Actual source code: matrart.c

petsc-master 2019-06-12
Report Typos and Errors

  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(Mat A)
 12: {
 14:   Mat_SeqAIJ     *a    = (Mat_SeqAIJ*)A->data;
 15:   Mat_RARt       *rart = a->rart;

 18:   MatTransposeColoringDestroy(&rart->matcoloring);
 19:   MatDestroy(&rart->Rt);
 20:   MatDestroy(&rart->RARt);
 21:   MatDestroy(&rart->ARt);
 22:   PetscFree(rart->work);

 24:   A->ops->destroy = rart->destroy;
 25:   if (A->ops->destroy) {
 26:     (*A->ops->destroy)(A);
 27:   }
 28:   PetscFree(rart);
 29:   return(0);
 30: }

 32: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat *C)
 33: {
 34:   PetscErrorCode       ierr;
 35:   Mat                  P;
 36:   PetscInt             *rti,*rtj;
 37:   Mat_RARt             *rart;
 38:   MatColoring          coloring;
 39:   MatTransposeColoring matcoloring;
 40:   ISColoring           iscoloring;
 41:   Mat                  Rt_dense,RARt_dense;
 42:   Mat_SeqAIJ           *c;

 45:   /* create symbolic P=Rt */
 46:   MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
 47:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);

 49:   /* get symbolic C=Pt*A*P */
 50:   MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
 51:   MatSetBlockSizes(*C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
 52:   (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;

 54:   /* create a supporting struct */
 55:   PetscNew(&rart);
 56:   c       = (Mat_SeqAIJ*)(*C)->data;
 57:   c->rart = rart;

 59:   /* ------ Use coloring ---------- */
 60:   /* inode causes memory problem, don't know why */
 61:   if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'");

 63:   /* Create MatTransposeColoring from symbolic C=R*A*R^T */
 64:   MatColoringCreate(*C,&coloring);
 65:   MatColoringSetDistance(coloring,2);
 66:   MatColoringSetType(coloring,MATCOLORINGSL);
 67:   MatColoringSetFromOptions(coloring);
 68:   MatColoringApply(coloring,&iscoloring);
 69:   MatColoringDestroy(&coloring);
 70:   MatTransposeColoringCreate(*C,iscoloring,&matcoloring);

 72:   rart->matcoloring = matcoloring;
 73:   ISColoringDestroy(&iscoloring);

 75:   /* Create Rt_dense */
 76:   MatCreate(PETSC_COMM_SELF,&Rt_dense);
 77:   MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
 78:   MatSetType(Rt_dense,MATSEQDENSE);
 79:   MatSeqDenseSetPreallocation(Rt_dense,NULL);

 81:   Rt_dense->assembled = PETSC_TRUE;
 82:   rart->Rt            = Rt_dense;

 84:   /* Create RARt_dense = R*A*Rt_dense */
 85:   MatCreate(PETSC_COMM_SELF,&RARt_dense);
 86:   MatSetSizes(RARt_dense,(*C)->rmap->n,matcoloring->ncolors,(*C)->rmap->n,matcoloring->ncolors);
 87:   MatSetType(RARt_dense,MATSEQDENSE);
 88:   MatSeqDenseSetPreallocation(RARt_dense,NULL);

 90:   rart->RARt = RARt_dense;

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

 95:   rart->destroy      = (*C)->ops->destroy;
 96:   (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt;

 98:   /* clean up */
 99:   MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
100:   MatDestroy(&P);

102: #if defined(PETSC_USE_INFO)
103:   {
104:     PetscReal density= (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
105:     PetscInfo(*C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
106:     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);
107:   }
108: #endif
109:   return(0);
110: }

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

129:   if (!dm || !dn) return(0);
130:   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);
131:   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);
132:   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);
133:   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);

135:   { /* 
136:      This approach is not as good as original ones (will be removed later), but it reveals that
137:      AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c
138:      */
139:     PetscBool via_matmatmult=PETSC_FALSE;
140:     PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
141:     if (via_matmatmult) {
142:       Mat AB_den;
143:       MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);
144:       MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
145:       MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
146:       MatDestroy(&AB_den);
147:       return(0);
148:     }
149:   }

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

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

209:     for (i=0; i<rm; i++) {  /* over rows of R in those columns */
210:       r1 = 0.0;
211:       n  = r->i[i+1] - r->i[i];
212:       rj = r->j + r->i[i];
213:       ra = r->a + r->i[i];
214:       for (j=0; j<n; j++) {
215:         r1 += (*ra++)*c[*rj++];
216:       }
217:       d[col*rm + i] = r1;
218:     }
219:   }
220:   PetscLogFlops(cn*2.0*(a->nz + r->nz));

222:   MatDenseRestoreArrayRead(B,&b);
223:   MatDenseRestoreArray(RAB,&d);
224:   MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
225:   MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
226:   return(0);
227: }

229: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
230: {
231:   PetscErrorCode       ierr;
232:   Mat_SeqAIJ           *c = (Mat_SeqAIJ*)C->data;
233:   Mat_RARt             *rart=c->rart;
234:   MatTransposeColoring matcoloring;
235:   Mat                  Rt,RARt;

238:   /* Get dense Rt by Apply MatTransposeColoring to R */
239:   matcoloring = rart->matcoloring;
240:   Rt          = rart->Rt;
241:   MatTransColoringApplySpToDen(matcoloring,R,Rt);

243:   /* Get dense RARt = R*A*Rt -- dominates! */
244:   RARt = rart->RARt;
245:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);

247:   /* Recover C from C_dense */
248:   MatTransColoringApplyDenToSp(matcoloring,RARt,C);
249:   return(0);
250: }

252: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C)
253: {
254:   PetscErrorCode  ierr;
255:   Mat             ARt,RARt;
256:   Mat_SeqAIJ     *c;
257:   Mat_RARt       *rart;

260:   /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */
261:   MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);
262:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);
263:   *C                     = RARt;
264:   RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;

266:   PetscNew(&rart);
267:   c         = (Mat_SeqAIJ*)(*C)->data;
268:   c->rart   = rart;
269:   rart->ARt = ARt;
270:   rart->destroy      = RARt->ops->destroy;
271:   RARt->ops->destroy = MatDestroy_SeqAIJ_RARt;
272: #if defined(PETSC_USE_INFO)
273:   PetscInfo(*C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
274: #endif
275:   return(0);
276: }

278: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
279: {
280:   PetscErrorCode  ierr;
281:   Mat_SeqAIJ      *c=(Mat_SeqAIJ*)C->data;
282:   Mat_RARt        *rart=c->rart;
283:   Mat             ARt=rart->ARt;
284: 
286:   MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt); /* dominate! */
287:   MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);
288:   return(0);
289: }

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

299:   MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
300:   MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);

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

315: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
316: {
317:   PetscErrorCode  ierr;
318:   Mat_SeqAIJ      *c = (Mat_SeqAIJ*)C->data;
319:   Mat_RARt        *rart = c->rart;
320:   Mat             Rt = rart->Rt;

323:   MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);
324:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);
325:   return(0);
326: }

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

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

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

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