Actual source code: matrart.c
petsc-3.12.5 2020-03-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, 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: }