Actual source code: matrart.c
petsc-3.13.6 2020-09-29
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 */
61: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);
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/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 = NULL;
143: MatCreate(PetscObjectComm((PetscObject)A),&AB_den);
144: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,AB_den);
145: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
146: MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
147: MatDestroy(&AB_den);
148: return(0);
149: }
150: }
152: MatDenseGetArrayRead(B,&b);
153: MatDenseGetArray(RAB,&d);
154: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
155: c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
156: for (col=0; col<cn-4; col += 4) { /* over columns of C */
157: for (i=0; i<am; i++) { /* over rows of A in those columns */
158: r1 = r2 = r3 = r4 = 0.0;
159: n = ai[i+1] - ai[i];
160: aj = a->j + ai[i];
161: aa = a->a + ai[i];
162: for (j=0; j<n; j++) {
163: r1 += (*aa)*b1[*aj];
164: r2 += (*aa)*b2[*aj];
165: r3 += (*aa)*b3[*aj];
166: r4 += (*aa++)*b4[*aj++];
167: }
168: c[i] = r1;
169: c[am + i] = r2;
170: c[am2 + i] = r3;
171: c[am3 + i] = r4;
172: }
173: b1 += bm4;
174: b2 += bm4;
175: b3 += bm4;
176: b4 += bm4;
178: /* RAB[:,col] = R*C[:,col] */
179: colrm = col*rm;
180: for (i=0; i<rm; i++) { /* over rows of R in those columns */
181: r1 = r2 = r3 = r4 = 0.0;
182: n = r->i[i+1] - r->i[i];
183: rj = r->j + r->i[i];
184: ra = r->a + r->i[i];
185: for (j=0; j<n; j++) {
186: r1 += (*ra)*c[*rj];
187: r2 += (*ra)*c2[*rj];
188: r3 += (*ra)*c3[*rj];
189: r4 += (*ra++)*c4[*rj++];
190: }
191: d[colrm + i] = r1;
192: d[colrm + rm + i] = r2;
193: d[colrm + rm2 + i] = r3;
194: d[colrm + rm3 + i] = r4;
195: }
196: }
197: for (; col<cn; col++) { /* over extra columns of C */
198: for (i=0; i<am; i++) { /* over rows of A in those columns */
199: r1 = 0.0;
200: n = a->i[i+1] - a->i[i];
201: aj = a->j + a->i[i];
202: aa = a->a + a->i[i];
203: for (j=0; j<n; j++) {
204: r1 += (*aa++)*b1[*aj++];
205: }
206: c[i] = r1;
207: }
208: b1 += bm;
210: for (i=0; i<rm; i++) { /* over rows of R in those columns */
211: r1 = 0.0;
212: n = r->i[i+1] - r->i[i];
213: rj = r->j + r->i[i];
214: ra = r->a + r->i[i];
215: for (j=0; j<n; j++) {
216: r1 += (*ra++)*c[*rj++];
217: }
218: d[col*rm + i] = r1;
219: }
220: }
221: PetscLogFlops(cn*2.0*(a->nz + r->nz));
223: MatDenseRestoreArrayRead(B,&b);
224: MatDenseRestoreArray(RAB,&d);
225: MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
226: MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
227: return(0);
228: }
230: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
231: {
232: PetscErrorCode ierr;
233: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
234: Mat_RARt *rart=c->rart;
235: MatTransposeColoring matcoloring;
236: Mat Rt,RARt;
239: /* Get dense Rt by Apply MatTransposeColoring to R */
240: matcoloring = rart->matcoloring;
241: Rt = rart->Rt;
242: MatTransColoringApplySpToDen(matcoloring,R,Rt);
244: /* Get dense RARt = R*A*Rt -- dominates! */
245: RARt = rart->RARt;
246: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);
248: /* Recover C from C_dense */
249: MatTransColoringApplyDenToSp(matcoloring,RARt,C);
250: return(0);
251: }
253: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat C)
254: {
256: Mat ARt;
257: Mat_SeqAIJ *c;
258: Mat_RARt *rart;
259: Mat_Product *product = C->product;
260: MatProductAlgorithm alg=product->alg;
263: /* create symbolic ARt = A*R^T */
264: MatProductCreate(A,R,NULL,&ARt);
265: MatProductSetType(ARt,MATPRODUCT_ABt);
266: MatProductSetAlgorithm(ARt,"sorted");
267: MatProductSetFill(ARt,fill);
268: MatProductSetFromOptions(ARt);
269: MatProductSymbolic(ARt);
271: /* compute symbolic C = R*ARt */
272: C->product->alg = "sorted"; /* set algorithm for C = R*ARt */
273: MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,C);
274: C->product->alg = alg; /* resume original algorithm for C */
276: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
278: PetscNew(&rart);
279: c = (Mat_SeqAIJ*)C->data;
280: c->rart = rart;
281: rart->ARt = ARt;
282: rart->destroy = C->ops->destroy;
283: C->ops->destroy = MatDestroy_SeqAIJ_RARt;
284: #if defined(PETSC_USE_INFO)
285: PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
286: #endif
287: return(0);
288: }
290: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
291: {
292: PetscErrorCode ierr;
293: Mat_SeqAIJ *c=(Mat_SeqAIJ*)C->data;
294: Mat_RARt *rart=c->rart;
295: Mat ARt=rart->ARt;
298: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt); /* dominate! */
299: MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);
300: return(0);
301: }
303: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat C)
304: {
306: Mat Rt;
307: Mat_SeqAIJ *c;
308: Mat_RARt *rart;
311: MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
312: MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);
314: PetscNew(&rart);
315: rart->Rt = Rt;
316: c = (Mat_SeqAIJ*)C->data;
317: c->rart = rart;
318: rart->destroy = C->ops->destroy;
319: C->ops->destroy = MatDestroy_SeqAIJ_RARt;
320: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
321: #if defined(PETSC_USE_INFO)
322: PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
323: #endif
324: return(0);
325: }
327: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
328: {
329: PetscErrorCode ierr;
330: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
331: Mat_RARt *rart = c->rart;
332: Mat Rt = rart->Rt;
335: MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);
336: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);
337: return(0);
338: }
340: PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
341: {
343: const char *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
344: PetscInt alg=0; /* set default algorithm */
347: if (scall == MAT_INITIAL_MATRIX) {
348: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatRARt","Mat");
349: PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
350: PetscOptionsEnd();
352: PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
353: MatCreate(PETSC_COMM_SELF,C);
354: switch (alg) {
355: case 1:
356: /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
357: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,*C);
358: break;
359: case 2:
360: /* via coloring_rart: apply coloring C = R*A*R^T */
361: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,*C);
362: break;
363: default:
364: /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
365: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,*C);
366: break;
367: }
368: PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
369: }
371: PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
372: ((*C)->ops->rartnumeric)(A,R,*C);
373: PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
374: return(0);
375: }
377: /* ------------------------------------------------------------- */
378: PetscErrorCode MatProductSymbolic_RARt_SeqAIJ_SeqAIJ(Mat C)
379: {
380: PetscErrorCode ierr;
381: Mat_Product *product = C->product;
382: Mat A=product->A,R=product->B;
383: MatProductAlgorithm alg=product->alg;
384: PetscReal fill=product->fill;
385: PetscBool flg;
388: PetscStrcmp(alg,"r*a*rt",&flg);
389: if (flg) {
390: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
391: goto next;
392: }
394: PetscStrcmp(alg,"r*art",&flg);
395: if (flg) {
396: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
397: goto next;
398: }
400: PetscStrcmp(alg,"coloring_rart",&flg);
401: if (flg) {
402: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
403: goto next;
404: }
406: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductAlgorithm is not supported");
408: next:
409: C->ops->productnumeric = MatProductNumeric_RARt;
410: return(0);
411: }