Actual source code: matrart.c


  2: /*
  3:   Defines projective product routines where A is a SeqAIJ matrix
  4:           C = R * A * R^T
  5: */

  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>

 11: PetscErrorCode MatDestroy_SeqAIJ_RARt(void *data)
 12: {
 13:   Mat_RARt       *rart = (Mat_RARt*)data;

 15:   MatTransposeColoringDestroy(&rart->matcoloring);
 16:   MatDestroy(&rart->Rt);
 17:   MatDestroy(&rart->RARt);
 18:   MatDestroy(&rart->ARt);
 19:   PetscFree(rart->work);
 20:   if (rart->destroy) {
 21:     (*rart->destroy)(rart->data);
 22:   }
 23:   PetscFree(rart);
 24:   return 0;
 25: }

 27: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat C)
 28: {
 29:   Mat                  P;
 30:   PetscInt             *rti,*rtj;
 31:   Mat_RARt             *rart;
 32:   MatColoring          coloring;
 33:   MatTransposeColoring matcoloring;
 34:   ISColoring           iscoloring;
 35:   Mat                  Rt_dense,RARt_dense;

 37:   MatCheckProduct(C,4);
 39:   /* create symbolic P=Rt */
 40:   MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
 41:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);

 43:   /* get symbolic C=Pt*A*P */
 44:   MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
 45:   MatSetBlockSizes(C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
 46:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;

 48:   /* create a supporting struct */
 49:   PetscNew(&rart);
 50:   C->product->data    = rart;
 51:   C->product->destroy = MatDestroy_SeqAIJ_RARt;

 53:   /* ------ Use coloring ---------- */
 54:   /* inode causes memory problem */
 55:   MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);

 57:   /* Create MatTransposeColoring from symbolic C=R*A*R^T */
 58:   MatColoringCreate(C,&coloring);
 59:   MatColoringSetDistance(coloring,2);
 60:   MatColoringSetType(coloring,MATCOLORINGSL);
 61:   MatColoringSetFromOptions(coloring);
 62:   MatColoringApply(coloring,&iscoloring);
 63:   MatColoringDestroy(&coloring);
 64:   MatTransposeColoringCreate(C,iscoloring,&matcoloring);

 66:   rart->matcoloring = matcoloring;
 67:   ISColoringDestroy(&iscoloring);

 69:   /* Create Rt_dense */
 70:   MatCreate(PETSC_COMM_SELF,&Rt_dense);
 71:   MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
 72:   MatSetType(Rt_dense,MATSEQDENSE);
 73:   MatSeqDenseSetPreallocation(Rt_dense,NULL);

 75:   Rt_dense->assembled = PETSC_TRUE;
 76:   rart->Rt            = Rt_dense;

 78:   /* Create RARt_dense = R*A*Rt_dense */
 79:   MatCreate(PETSC_COMM_SELF,&RARt_dense);
 80:   MatSetSizes(RARt_dense,C->rmap->n,matcoloring->ncolors,C->rmap->n,matcoloring->ncolors);
 81:   MatSetType(RARt_dense,MATSEQDENSE);
 82:   MatSeqDenseSetPreallocation(RARt_dense,NULL);

 84:   rart->RARt = RARt_dense;

 86:   /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
 87:   PetscMalloc1(A->rmap->n*4,&rart->work);

 89:   /* clean up */
 90:   MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
 91:   MatDestroy(&P);

 93: #if defined(PETSC_USE_INFO)
 94:   {
 95:     Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
 96:     PetscReal density = (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
 97:     PetscInfo(C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
 98:     PetscInfo(C,"RARt_den %" PetscInt_FMT " %" PetscInt_FMT "; Rt %" PetscInt_FMT " %" PetscInt_FMT " (RARt->nz %" PetscInt_FMT ")/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,R->cmap->n,R->rmap->n,c->nz,density);
 99:   }
100: #endif
101:   return 0;
102: }

104: /*
105:  RAB = R * A * B, R and A in seqaij format, B in dense format;
106: */
107: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
108: {
109:   Mat_SeqAIJ        *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
110:   PetscScalar       r1,r2,r3,r4;
111:   const PetscScalar *b,*b1,*b2,*b3,*b4;
112:   MatScalar         *aa,*ra;
113:   PetscInt          cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
114:   PetscInt          am2=2*am,am3=3*am,bm4=4*bm;
115:   PetscScalar       *d,*c,*c2,*c3,*c4;
116:   PetscInt          *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
117:   PetscInt         rm2=2*rm,rm3=3*rm,colrm;

119:   if (!dm || !dn) return 0;

125:   { /*
126:      This approach is not as good as original ones (will be removed later), but it reveals that
127:      AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/tutorials/ex56.c
128:      */
129:     PetscBool via_matmatmult=PETSC_FALSE;
130:     PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
131:     if (via_matmatmult) {
132:       Mat AB_den = NULL;
133:       MatCreate(PetscObjectComm((PetscObject)A),&AB_den);
134:       MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,AB_den);
135:       MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
136:       MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
137:       MatDestroy(&AB_den);
138:       return 0;
139:     }
140:   }

142:   MatDenseGetArrayRead(B,&b);
143:   MatDenseGetArray(RAB,&d);
144:   b1   = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
145:   c    = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
146:   for (col=0; col<cn-4; col += 4) {  /* over columns of C */
147:     for (i=0; i<am; i++) {        /* over rows of A in those columns */
148:       r1 = r2 = r3 = r4 = 0.0;
149:       n  = ai[i+1] - ai[i];
150:       aj = a->j + ai[i];
151:       aa = a->a + ai[i];
152:       for (j=0; j<n; j++) {
153:         r1 += (*aa)*b1[*aj];
154:         r2 += (*aa)*b2[*aj];
155:         r3 += (*aa)*b3[*aj];
156:         r4 += (*aa++)*b4[*aj++];
157:       }
158:       c[i]       = r1;
159:       c[am  + i] = r2;
160:       c[am2 + i] = r3;
161:       c[am3 + i] = r4;
162:     }
163:     b1 += bm4;
164:     b2 += bm4;
165:     b3 += bm4;
166:     b4 += bm4;

168:     /* RAB[:,col] = R*C[:,col] */
169:     colrm = col*rm;
170:     for (i=0; i<rm; i++) {        /* over rows of R in those columns */
171:       r1 = r2 = r3 = r4 = 0.0;
172:       n  = r->i[i+1] - r->i[i];
173:       rj = r->j + r->i[i];
174:       ra = r->a + r->i[i];
175:       for (j=0; j<n; j++) {
176:         r1 += (*ra)*c[*rj];
177:         r2 += (*ra)*c2[*rj];
178:         r3 += (*ra)*c3[*rj];
179:         r4 += (*ra++)*c4[*rj++];
180:       }
181:       d[colrm + i]       = r1;
182:       d[colrm + rm + i]  = r2;
183:       d[colrm + rm2 + i] = r3;
184:       d[colrm + rm3 + i] = r4;
185:     }
186:   }
187:   for (; col<cn; col++) {     /* over extra columns of C */
188:     for (i=0; i<am; i++) {  /* over rows of A in those columns */
189:       r1 = 0.0;
190:       n  = a->i[i+1] - a->i[i];
191:       aj = a->j + a->i[i];
192:       aa = a->a + a->i[i];
193:       for (j=0; j<n; j++) {
194:         r1 += (*aa++)*b1[*aj++];
195:       }
196:       c[i] = r1;
197:     }
198:     b1 += bm;

200:     for (i=0; i<rm; i++) {  /* over rows of R in those columns */
201:       r1 = 0.0;
202:       n  = r->i[i+1] - r->i[i];
203:       rj = r->j + r->i[i];
204:       ra = r->a + r->i[i];
205:       for (j=0; j<n; j++) {
206:         r1 += (*ra++)*c[*rj++];
207:       }
208:       d[col*rm + i] = r1;
209:     }
210:   }
211:   PetscLogFlops(cn*2.0*(a->nz + r->nz));

213:   MatDenseRestoreArrayRead(B,&b);
214:   MatDenseRestoreArray(RAB,&d);
215:   MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
216:   MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
217:   return 0;
218: }

220: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
221: {
222:   Mat_RARt             *rart;
223:   MatTransposeColoring matcoloring;
224:   Mat                  Rt,RARt;

226:   MatCheckProduct(C,3);
228:   rart = (Mat_RARt*)C->product->data;

230:   /* Get dense Rt by Apply MatTransposeColoring to R */
231:   matcoloring = rart->matcoloring;
232:   Rt          = rart->Rt;
233:   MatTransColoringApplySpToDen(matcoloring,R,Rt);

235:   /* Get dense RARt = R*A*Rt -- dominates! */
236:   RARt = rart->RARt;
237:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);

239:   /* Recover C from C_dense */
240:   MatTransColoringApplyDenToSp(matcoloring,RARt,C);
241:   return 0;
242: }

244: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat C)
245: {
246:   Mat            ARt;
247:   Mat_RARt       *rart;
248:   char           *alg;

250:   MatCheckProduct(C,4);
252:   /* create symbolic ARt = A*R^T  */
253:   MatProductCreate(A,R,NULL,&ARt);
254:   MatProductSetType(ARt,MATPRODUCT_ABt);
255:   MatProductSetAlgorithm(ARt,"sorted");
256:   MatProductSetFill(ARt,fill);
257:   MatProductSetFromOptions(ARt);
258:   MatProductSymbolic(ARt);

260:   /* compute symbolic C = R*ARt */
261:   /* set algorithm for C = R*ARt */
262:   PetscStrallocpy(C->product->alg,&alg);
263:   MatProductSetAlgorithm(C,"sorted");
264:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,C);
265:   /* resume original algorithm for C */
266:   MatProductSetAlgorithm(C,alg);
267:   PetscFree(alg);

269:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;

271:   PetscNew(&rart);
272:   rart->ARt = ARt;
273:   C->product->data    = rart;
274:   C->product->destroy = MatDestroy_SeqAIJ_RARt;
275:   PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
276:   return 0;
277: }

279: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
280: {
281:   Mat_RARt       *rart;

283:   MatCheckProduct(C,3);
285:   rart = (Mat_RARt*)C->product->data;
286:   MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,rart->ARt); /* dominate! */
287:   MatMatMultNumeric_SeqAIJ_SeqAIJ(R,rart->ARt,C);
288:   return 0;
289: }

291: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat C)
292: {
293:   Mat            Rt;
294:   Mat_RARt       *rart;

296:   MatCheckProduct(C,4);
298:   MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
299:   MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);

301:   PetscNew(&rart);
302:   rart->data = C->product->data;
303:   rart->destroy = C->product->destroy;
304:   rart->Rt = Rt;
305:   C->product->data    = rart;
306:   C->product->destroy = MatDestroy_SeqAIJ_RARt;
307:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
308:   PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
309:   return 0;
310: }

312: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
313: {
314:   Mat_RARt       *rart;

316:   MatCheckProduct(C,3);
318:   rart = (Mat_RARt*)C->product->data;
319:   MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&rart->Rt);
320:   /* MatMatMatMultSymbolic used a different data */
321:   C->product->data = rart->data;
322:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,rart->Rt,C);
323:   C->product->data = rart;
324:   return 0;
325: }

327: PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
328: {
330:   const char     *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
331:   PetscInt       alg=0; /* set default algorithm */

333:   if (scall == MAT_INITIAL_MATRIX) {
334:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatRARt","Mat");
335:     PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
336:     PetscOptionsEnd();

338:     PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
339:     MatCreate(PETSC_COMM_SELF,C);
340:     switch (alg) {
341:     case 1:
342:       /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
343:       MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,*C);
344:       break;
345:     case 2:
346:       /* via coloring_rart: apply coloring C = R*A*R^T                          */
347:       MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,*C);
348:       break;
349:     default:
350:       /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
351:       MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,*C);
352:       break;
353:     }
354:     PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
355:   }

357:   PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
358:   ((*C)->ops->rartnumeric)(A,R,*C);
359:   PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
360:   return 0;
361: }

363: /* ------------------------------------------------------------- */
364: PetscErrorCode MatProductSymbolic_RARt_SeqAIJ_SeqAIJ(Mat C)
365: {
366:   Mat_Product         *product = C->product;
367:   Mat                 A=product->A,R=product->B;
368:   MatProductAlgorithm alg=product->alg;
369:   PetscReal           fill=product->fill;
370:   PetscBool           flg;

372:   PetscStrcmp(alg,"r*a*rt",&flg);
373:   if (flg) {
374:     MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
375:     goto next;
376:   }

378:   PetscStrcmp(alg,"r*art",&flg);
379:   if (flg) {
380:     MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
381:     goto next;
382:   }

384:   PetscStrcmp(alg,"coloring_rart",&flg);
385:   if (flg) {
386:     MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
387:     goto next;
388:   }

390:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductAlgorithm is not supported");

392: next:
393:   C->ops->productnumeric = MatProductNumeric_RARt;
394:   return 0;
395: }