1 2 /* 3 Defines projective product routines where A is a SeqAIJ matrix 4 C = R * A * R^T 5 */ 6 7 #include <../src/mat/impls/aij/seq/aij.h> 8 #include <../src/mat/utils/freespace.h> 9 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 10 #include <petsctime.h> 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatDestroy_SeqAIJ_RARt" 14 PetscErrorCode MatDestroy_SeqAIJ_RARt(Mat A) 15 { 16 PetscErrorCode ierr; 17 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 18 Mat_RARt *rart = a->rart; 19 20 PetscFunctionBegin; 21 ierr = MatTransposeColoringDestroy(&rart->matcoloring);CHKERRQ(ierr); 22 ierr = MatDestroy(&rart->Rt);CHKERRQ(ierr); 23 ierr = MatDestroy(&rart->RARt);CHKERRQ(ierr); 24 ierr = MatDestroy(&rart->ARt);CHKERRQ(ierr); 25 ierr = PetscFree(rart->work);CHKERRQ(ierr); 26 27 A->ops->destroy = rart->destroy; 28 if (A->ops->destroy) { 29 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 30 } 31 ierr = PetscFree(rart);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart" 37 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat *C) 38 { 39 PetscErrorCode ierr; 40 Mat P; 41 PetscInt *rti,*rtj; 42 Mat_RARt *rart; 43 MatColoring coloring; 44 MatTransposeColoring matcoloring; 45 ISColoring iscoloring; 46 Mat Rt_dense,RARt_dense; 47 Mat_SeqAIJ *c; 48 49 PetscFunctionBegin; 50 /* create symbolic P=Rt */ 51 ierr = MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr); 52 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);CHKERRQ(ierr); 53 54 /* get symbolic C=Pt*A*P */ 55 ierr = MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);CHKERRQ(ierr); 56 (*C)->rmap->bs = R->rmap->bs; 57 (*C)->cmap->bs = R->rmap->bs; 58 (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart; 59 60 /* create a supporting struct */ 61 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 62 c = (Mat_SeqAIJ*)(*C)->data; 63 c->rart = rart; 64 65 /* ------ Use coloring ---------- */ 66 /* inode causes memory problem, don't know why */ 67 if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'"); 68 69 /* Create MatTransposeColoring from symbolic C=R*A*R^T */ 70 ierr = MatColoringCreate(*C,&coloring);CHKERRQ(ierr); 71 ierr = MatColoringSetDistance(coloring,2);CHKERRQ(ierr); 72 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 73 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 74 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 75 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 76 ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr); 77 78 rart->matcoloring = matcoloring; 79 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 80 81 /* Create Rt_dense */ 82 ierr = MatCreate(PETSC_COMM_SELF,&Rt_dense);CHKERRQ(ierr); 83 ierr = MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr); 84 ierr = MatSetType(Rt_dense,MATSEQDENSE);CHKERRQ(ierr); 85 ierr = MatSeqDenseSetPreallocation(Rt_dense,NULL);CHKERRQ(ierr); 86 87 Rt_dense->assembled = PETSC_TRUE; 88 rart->Rt = Rt_dense; 89 90 /* Create RARt_dense = R*A*Rt_dense */ 91 ierr = MatCreate(PETSC_COMM_SELF,&RARt_dense);CHKERRQ(ierr); 92 ierr = MatSetSizes(RARt_dense,(*C)->rmap->n,matcoloring->ncolors,(*C)->rmap->n,matcoloring->ncolors);CHKERRQ(ierr); 93 ierr = MatSetType(RARt_dense,MATSEQDENSE);CHKERRQ(ierr); 94 ierr = MatSeqDenseSetPreallocation(RARt_dense,NULL);CHKERRQ(ierr); 95 96 rart->RARt = RARt_dense; 97 98 /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */ 99 ierr = PetscMalloc(A->rmap->n*4*sizeof(PetscScalar),&rart->work);CHKERRQ(ierr); 100 101 rart->destroy = (*C)->ops->destroy; 102 (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt; 103 104 /* clean up */ 105 ierr = MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr); 106 ierr = MatDestroy(&P);CHKERRQ(ierr); 107 108 #if defined(PETSC_USE_INFO) 109 { 110 PetscReal density= (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n); 111 ierr = PetscInfo(*C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");CHKERRQ(ierr); 112 ierr = PetscInfo6(*C,"RARt_den %D %D; Rt %D %D (RARt->nz %D)/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,R->cmap->n,R->rmap->n,c->nz,density);CHKERRQ(ierr); 113 } 114 #endif 115 PetscFunctionReturn(0); 116 } 117 118 /* 119 RAB = R * A * B, R and A in seqaij format, B in dense format; 120 */ 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense" 123 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work) 124 { 125 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data; 126 PetscErrorCode ierr; 127 PetscScalar *b,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 128 MatScalar *aa,*ra; 129 PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n; 130 PetscInt am2=2*am,am3=3*am,bm4=4*bm; 131 PetscScalar *d,*c,*c2,*c3,*c4; 132 PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n; 133 PetscInt rm2=2*rm,rm3=3*rm,colrm; 134 135 PetscFunctionBegin; 136 if (!dm || !dn) PetscFunctionReturn(0); 137 if (bm != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,bm); 138 if (am != R->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in R %D not equal rows in A %D\n",R->cmap->n,am); 139 if (R->rmap->n != RAB->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in RAB %D not equal rows in R %D\n",RAB->rmap->n,R->rmap->n); 140 if (B->cmap->n != RAB->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in RAB %D not equal columns in B %D\n",RAB->cmap->n,B->cmap->n); 141 142 { /* 143 This approach is not as good as original ones (will be removed later), but it reveals that 144 AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c 145 */ 146 PetscBool via_matmatmult=PETSC_FALSE; 147 ierr = PetscOptionsGetBool(NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);CHKERRQ(ierr); 148 if (via_matmatmult) { 149 Mat AB_den; 150 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);CHKERRQ(ierr); 151 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);CHKERRQ(ierr); 152 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);CHKERRQ(ierr); 153 ierr = MatDestroy(&AB_den);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 } 157 158 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 159 ierr = MatDenseGetArray(RAB,&d);CHKERRQ(ierr); 160 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 161 c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am; 162 for (col=0; col<cn-4; col += 4) { /* over columns of C */ 163 for (i=0; i<am; i++) { /* over rows of A in those columns */ 164 r1 = r2 = r3 = r4 = 0.0; 165 n = ai[i+1] - ai[i]; 166 aj = a->j + ai[i]; 167 aa = a->a + ai[i]; 168 for (j=0; j<n; j++) { 169 r1 += (*aa)*b1[*aj]; 170 r2 += (*aa)*b2[*aj]; 171 r3 += (*aa)*b3[*aj]; 172 r4 += (*aa++)*b4[*aj++]; 173 } 174 c[i] = r1; 175 c[am + i] = r2; 176 c[am2 + i] = r3; 177 c[am3 + i] = r4; 178 } 179 b1 += bm4; 180 b2 += bm4; 181 b3 += bm4; 182 b4 += bm4; 183 184 /* RAB[:,col] = R*C[:,col] */ 185 colrm = col*rm; 186 for (i=0; i<rm; i++) { /* over rows of R in those columns */ 187 r1 = r2 = r3 = r4 = 0.0; 188 n = r->i[i+1] - r->i[i]; 189 rj = r->j + r->i[i]; 190 ra = r->a + r->i[i]; 191 for (j=0; j<n; j++) { 192 r1 += (*ra)*c[*rj]; 193 r2 += (*ra)*c2[*rj]; 194 r3 += (*ra)*c3[*rj]; 195 r4 += (*ra++)*c4[*rj++]; 196 } 197 d[colrm + i] = r1; 198 d[colrm + rm + i] = r2; 199 d[colrm + rm2 + i] = r3; 200 d[colrm + rm3 + i] = r4; 201 } 202 } 203 for (; col<cn; col++) { /* over extra columns of C */ 204 for (i=0; i<am; i++) { /* over rows of A in those columns */ 205 r1 = 0.0; 206 n = a->i[i+1] - a->i[i]; 207 aj = a->j + a->i[i]; 208 aa = a->a + a->i[i]; 209 for (j=0; j<n; j++) { 210 r1 += (*aa++)*b1[*aj++]; 211 } 212 c[i] = r1; 213 } 214 b1 += bm; 215 216 for (i=0; i<rm; i++) { /* over rows of R in those columns */ 217 r1 = 0.0; 218 n = r->i[i+1] - r->i[i]; 219 rj = r->j + r->i[i]; 220 ra = r->a + r->i[i]; 221 for (j=0; j<n; j++) { 222 r1 += (*ra++)*c[*rj++]; 223 } 224 d[col*rm + i] = r1; 225 } 226 } 227 ierr = PetscLogFlops(cn*2.0*(a->nz + r->nz));CHKERRQ(ierr); 228 229 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 230 ierr = MatDenseRestoreArray(RAB,&d);CHKERRQ(ierr); 231 ierr = MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 232 ierr = MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart" 238 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C) 239 { 240 PetscErrorCode ierr; 241 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 242 Mat_RARt *rart=c->rart; 243 MatTransposeColoring matcoloring; 244 Mat Rt,RARt; 245 246 PetscFunctionBegin; 247 /* Get dense Rt by Apply MatTransposeColoring to R */ 248 matcoloring = rart->matcoloring; 249 Rt = rart->Rt; 250 ierr = MatTransColoringApplySpToDen(matcoloring,R,Rt);CHKERRQ(ierr); 251 252 /* Get dense RARt = R*A*Rt -- dominates! */ 253 RARt = rart->RARt; 254 ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);CHKERRQ(ierr); 255 256 /* Recover C from C_dense */ 257 ierr = MatTransColoringApplyDenToSp(matcoloring,RARt,C);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult" 263 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C) 264 { 265 PetscErrorCode ierr; 266 Mat ARt,RARt; 267 Mat_SeqAIJ *c; 268 Mat_RARt *rart; 269 270 PetscFunctionBegin; 271 /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */ 272 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);CHKERRQ(ierr); 273 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);CHKERRQ(ierr); 274 *C = RARt; 275 RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult; 276 277 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 278 c = (Mat_SeqAIJ*)(*C)->data; 279 c->rart = rart; 280 rart->ARt = ARt; 281 rart->destroy = RARt->ops->destroy; 282 RARt->ops->destroy = MatDestroy_SeqAIJ_RARt; 283 #if defined(PETSC_USE_INFO) 284 ierr = PetscInfo(*C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");CHKERRQ(ierr); 285 #endif 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult" 291 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C) 292 { 293 PetscErrorCode ierr; 294 Mat_SeqAIJ *c=(Mat_SeqAIJ*)C->data; 295 Mat_RARt *rart=c->rart; 296 Mat ARt=rart->ARt; 297 298 PetscFunctionBegin; 299 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt);CHKERRQ(ierr); /* dominate! */ 300 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ" 306 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat *C) 307 { 308 PetscErrorCode ierr; 309 Mat Rt; 310 Mat_SeqAIJ *c; 311 Mat_RARt *rart; 312 313 PetscFunctionBegin; 314 ierr = MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); 315 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);CHKERRQ(ierr); 316 317 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 318 rart->Rt = Rt; 319 c = (Mat_SeqAIJ*)(*C)->data; 320 c->rart = rart; 321 rart->destroy = (*C)->ops->destroy; 322 (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt; 323 (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ; 324 #if defined(PETSC_USE_INFO) 325 ierr = PetscInfo(*C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");CHKERRQ(ierr); 326 #endif 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ" 332 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C) 333 { 334 PetscErrorCode ierr; 335 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 336 Mat_RARt *rart = c->rart; 337 Mat Rt = rart->Rt; 338 339 PetscFunctionBegin; 340 ierr = MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);CHKERRQ(ierr); 341 ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344 345 #undef __FUNCT__ 346 #define __FUNCT__ "MatRARt_SeqAIJ_SeqAIJ" 347 PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C) 348 { 349 PetscErrorCode ierr; 350 const char *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"}; 351 PetscInt alg=0; /* set default algorithm */ 352 353 PetscFunctionBegin; 354 if (scall == MAT_INITIAL_MATRIX) { 355 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 356 ierr = PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);CHKERRQ(ierr); 357 ierr = PetscOptionsEnd();CHKERRQ(ierr); 358 359 ierr = PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr); 360 switch (alg) { 361 case 1: 362 /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */ 363 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);CHKERRQ(ierr); 364 break; 365 case 2: 366 /* via coloring_rart: apply coloring C = R*A*R^T */ 367 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);CHKERRQ(ierr); 368 break; 369 default: 370 /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */ 371 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);CHKERRQ(ierr); 372 break; 373 } 374 ierr = PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr); 375 } 376 377 ierr = PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr); 378 ierr = (*(*C)->ops->rartnumeric)(A,R,*C);CHKERRQ(ierr); 379 ierr = PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr); 380 PetscFunctionReturn(0); 381 } 382