Actual source code: matrart.c
petsc-3.14.6 2021-03-30
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: {
14: Mat_RARt *rart = (Mat_RARt*)data;
17: MatTransposeColoringDestroy(&rart->matcoloring);
18: MatDestroy(&rart->Rt);
19: MatDestroy(&rart->RARt);
20: MatDestroy(&rart->ARt);
21: PetscFree(rart->work);
22: if (rart->destroy) {
23: (*rart->destroy)(rart->data);
24: }
25: PetscFree(rart);
26: return(0);
27: }
29: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat C)
30: {
31: PetscErrorCode ierr;
32: Mat P;
33: PetscInt *rti,*rtj;
34: Mat_RARt *rart;
35: MatColoring coloring;
36: MatTransposeColoring matcoloring;
37: ISColoring iscoloring;
38: Mat Rt_dense,RARt_dense;
41: MatCheckProduct(C,4);
42: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
43: /* create symbolic P=Rt */
44: MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
45: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);
47: /* get symbolic C=Pt*A*P */
48: MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
49: MatSetBlockSizes(C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
50: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;
52: /* create a supporting struct */
53: PetscNew(&rart);
54: C->product->data = rart;
55: C->product->destroy = MatDestroy_SeqAIJ_RARt;
57: /* ------ Use coloring ---------- */
58: /* inode causes memory problem */
59: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);
61: /* Create MatTransposeColoring from symbolic C=R*A*R^T */
62: MatColoringCreate(C,&coloring);
63: MatColoringSetDistance(coloring,2);
64: MatColoringSetType(coloring,MATCOLORINGSL);
65: MatColoringSetFromOptions(coloring);
66: MatColoringApply(coloring,&iscoloring);
67: MatColoringDestroy(&coloring);
68: MatTransposeColoringCreate(C,iscoloring,&matcoloring);
70: rart->matcoloring = matcoloring;
71: ISColoringDestroy(&iscoloring);
73: /* Create Rt_dense */
74: MatCreate(PETSC_COMM_SELF,&Rt_dense);
75: MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
76: MatSetType(Rt_dense,MATSEQDENSE);
77: MatSeqDenseSetPreallocation(Rt_dense,NULL);
79: Rt_dense->assembled = PETSC_TRUE;
80: rart->Rt = Rt_dense;
82: /* Create RARt_dense = R*A*Rt_dense */
83: MatCreate(PETSC_COMM_SELF,&RARt_dense);
84: MatSetSizes(RARt_dense,C->rmap->n,matcoloring->ncolors,C->rmap->n,matcoloring->ncolors);
85: MatSetType(RARt_dense,MATSEQDENSE);
86: MatSeqDenseSetPreallocation(RARt_dense,NULL);
88: rart->RARt = RARt_dense;
90: /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
91: PetscMalloc1(A->rmap->n*4,&rart->work);
93: /* clean up */
94: MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
95: MatDestroy(&P);
97: #if defined(PETSC_USE_INFO)
98: {
99: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
100: PetscReal density = (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
101: PetscInfo(C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
102: 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);
103: }
104: #endif
105: return(0);
106: }
108: /*
109: RAB = R * A * B, R and A in seqaij format, B in dense format;
110: */
111: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
112: {
113: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
114: PetscErrorCode ierr;
115: PetscScalar r1,r2,r3,r4;
116: const PetscScalar *b,*b1,*b2,*b3,*b4;
117: MatScalar *aa,*ra;
118: PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
119: PetscInt am2=2*am,am3=3*am,bm4=4*bm;
120: PetscScalar *d,*c,*c2,*c3,*c4;
121: PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
122: PetscInt rm2=2*rm,rm3=3*rm,colrm;
125: if (!dm || !dn) return(0);
126: 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);
127: 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);
128: 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);
129: 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);
131: { /*
132: This approach is not as good as original ones (will be removed later), but it reveals that
133: AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/tutorials/ex56.c
134: */
135: PetscBool via_matmatmult=PETSC_FALSE;
136: PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
137: if (via_matmatmult) {
138: Mat AB_den = NULL;
139: MatCreate(PetscObjectComm((PetscObject)A),&AB_den);
140: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,AB_den);
141: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
142: MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
143: MatDestroy(&AB_den);
144: return(0);
145: }
146: }
148: MatDenseGetArrayRead(B,&b);
149: MatDenseGetArray(RAB,&d);
150: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
151: c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
152: for (col=0; col<cn-4; col += 4) { /* over columns of C */
153: for (i=0; i<am; i++) { /* over rows of A in those columns */
154: r1 = r2 = r3 = r4 = 0.0;
155: n = ai[i+1] - ai[i];
156: aj = a->j + ai[i];
157: aa = a->a + ai[i];
158: for (j=0; j<n; j++) {
159: r1 += (*aa)*b1[*aj];
160: r2 += (*aa)*b2[*aj];
161: r3 += (*aa)*b3[*aj];
162: r4 += (*aa++)*b4[*aj++];
163: }
164: c[i] = r1;
165: c[am + i] = r2;
166: c[am2 + i] = r3;
167: c[am3 + i] = r4;
168: }
169: b1 += bm4;
170: b2 += bm4;
171: b3 += bm4;
172: b4 += bm4;
174: /* RAB[:,col] = R*C[:,col] */
175: colrm = col*rm;
176: for (i=0; i<rm; i++) { /* over rows of R in those columns */
177: r1 = r2 = r3 = r4 = 0.0;
178: n = r->i[i+1] - r->i[i];
179: rj = r->j + r->i[i];
180: ra = r->a + r->i[i];
181: for (j=0; j<n; j++) {
182: r1 += (*ra)*c[*rj];
183: r2 += (*ra)*c2[*rj];
184: r3 += (*ra)*c3[*rj];
185: r4 += (*ra++)*c4[*rj++];
186: }
187: d[colrm + i] = r1;
188: d[colrm + rm + i] = r2;
189: d[colrm + rm2 + i] = r3;
190: d[colrm + rm3 + i] = r4;
191: }
192: }
193: for (; col<cn; col++) { /* over extra columns of C */
194: for (i=0; i<am; i++) { /* over rows of A in those columns */
195: r1 = 0.0;
196: n = a->i[i+1] - a->i[i];
197: aj = a->j + a->i[i];
198: aa = a->a + a->i[i];
199: for (j=0; j<n; j++) {
200: r1 += (*aa++)*b1[*aj++];
201: }
202: c[i] = r1;
203: }
204: b1 += bm;
206: for (i=0; i<rm; i++) { /* over rows of R in those columns */
207: r1 = 0.0;
208: n = r->i[i+1] - r->i[i];
209: rj = r->j + r->i[i];
210: ra = r->a + r->i[i];
211: for (j=0; j<n; j++) {
212: r1 += (*ra++)*c[*rj++];
213: }
214: d[col*rm + i] = r1;
215: }
216: }
217: PetscLogFlops(cn*2.0*(a->nz + r->nz));
219: MatDenseRestoreArrayRead(B,&b);
220: MatDenseRestoreArray(RAB,&d);
221: MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
222: MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
223: return(0);
224: }
226: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
227: {
228: PetscErrorCode ierr;
229: Mat_RARt *rart;
230: MatTransposeColoring matcoloring;
231: Mat Rt,RARt;
234: MatCheckProduct(C,3);
235: if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
236: rart = (Mat_RARt*)C->product->data;
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: {
255: Mat ARt;
256: Mat_RARt *rart;
257: char *alg;
260: MatCheckProduct(C,4);
261: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
262: /* create symbolic ARt = A*R^T */
263: MatProductCreate(A,R,NULL,&ARt);
264: MatProductSetType(ARt,MATPRODUCT_ABt);
265: MatProductSetAlgorithm(ARt,"sorted");
266: MatProductSetFill(ARt,fill);
267: MatProductSetFromOptions(ARt);
268: MatProductSymbolic(ARt);
270: /* compute symbolic C = R*ARt */
271: /* set algorithm for C = R*ARt */
272: PetscStrallocpy(C->product->alg,&alg);
273: MatProductSetAlgorithm(C,"sorted");
274: MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,C);
275: /* resume original algorithm for C */
276: MatProductSetAlgorithm(C,alg);
277: PetscFree(alg);
279: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
281: PetscNew(&rart);
282: rart->ARt = ARt;
283: C->product->data = rart;
284: C->product->destroy = MatDestroy_SeqAIJ_RARt;
285: PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
286: return(0);
287: }
289: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
290: {
292: Mat_RARt *rart;
295: MatCheckProduct(C,3);
296: if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
297: rart = (Mat_RARt*)C->product->data;
298: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,rart->ARt); /* dominate! */
299: MatMatMultNumeric_SeqAIJ_SeqAIJ(R,rart->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_RARt *rart;
310: MatCheckProduct(C,4);
311: if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
312: MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
313: MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);
315: PetscNew(&rart);
316: rart->data = C->product->data;
317: rart->destroy = C->product->destroy;
318: rart->Rt = Rt;
319: C->product->data = rart;
320: C->product->destroy = MatDestroy_SeqAIJ_RARt;
321: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
322: PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
323: return(0);
324: }
326: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
327: {
329: Mat_RARt *rart;
332: MatCheckProduct(C,3);
333: if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
334: rart = (Mat_RARt*)C->product->data;
335: MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&rart->Rt);
336: /* MatMatMatMultSymbolic used a different data */
337: C->product->data = rart->data;
338: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,rart->Rt,C);
339: C->product->data = rart;
340: return(0);
341: }
343: PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
344: {
346: const char *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
347: PetscInt alg=0; /* set default algorithm */
350: if (scall == MAT_INITIAL_MATRIX) {
351: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatRARt","Mat");
352: PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
353: PetscOptionsEnd();
355: PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
356: MatCreate(PETSC_COMM_SELF,C);
357: switch (alg) {
358: case 1:
359: /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
360: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,*C);
361: break;
362: case 2:
363: /* via coloring_rart: apply coloring C = R*A*R^T */
364: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,*C);
365: break;
366: default:
367: /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
368: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,*C);
369: break;
370: }
371: PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
372: }
374: PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
375: ((*C)->ops->rartnumeric)(A,R,*C);
376: PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
377: return(0);
378: }
380: /* ------------------------------------------------------------- */
381: PetscErrorCode MatProductSymbolic_RARt_SeqAIJ_SeqAIJ(Mat C)
382: {
383: PetscErrorCode ierr;
384: Mat_Product *product = C->product;
385: Mat A=product->A,R=product->B;
386: MatProductAlgorithm alg=product->alg;
387: PetscReal fill=product->fill;
388: PetscBool flg;
391: PetscStrcmp(alg,"r*a*rt",&flg);
392: if (flg) {
393: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
394: goto next;
395: }
397: PetscStrcmp(alg,"r*art",&flg);
398: if (flg) {
399: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
400: goto next;
401: }
403: PetscStrcmp(alg,"coloring_rart",&flg);
404: if (flg) {
405: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
406: goto next;
407: }
409: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductAlgorithm is not supported");
411: next:
412: C->ops->productnumeric = MatProductNumeric_RARt;
413: return(0);
414: }