Actual source code: matrart.c
petsc-3.9.4 2018-09-11
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;
119: PetscScalar *b,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
120: MatScalar *aa,*ra;
121: PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
122: PetscInt am2=2*am,am3=3*am,bm4=4*bm;
123: PetscScalar *d,*c,*c2,*c3,*c4;
124: PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
125: PetscInt rm2=2*rm,rm3=3*rm,colrm;
128: if (!dm || !dn) return(0);
129: 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);
130: 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);
131: 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);
132: 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);
134: { /*
135: This approach is not as good as original ones (will be removed later), but it reveals that
136: AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c
137: */
138: PetscBool via_matmatmult=PETSC_FALSE;
139: PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
140: if (via_matmatmult) {
141: Mat AB_den;
142: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);
143: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
144: MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
145: MatDestroy(&AB_den);
146: return(0);
147: }
148: }
150: MatDenseGetArray(B,&b);
151: MatDenseGetArray(RAB,&d);
152: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
153: c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
154: for (col=0; col<cn-4; col += 4) { /* over columns of C */
155: for (i=0; i<am; i++) { /* over rows of A in those columns */
156: r1 = r2 = r3 = r4 = 0.0;
157: n = ai[i+1] - ai[i];
158: aj = a->j + ai[i];
159: aa = a->a + ai[i];
160: for (j=0; j<n; j++) {
161: r1 += (*aa)*b1[*aj];
162: r2 += (*aa)*b2[*aj];
163: r3 += (*aa)*b3[*aj];
164: r4 += (*aa++)*b4[*aj++];
165: }
166: c[i] = r1;
167: c[am + i] = r2;
168: c[am2 + i] = r3;
169: c[am3 + i] = r4;
170: }
171: b1 += bm4;
172: b2 += bm4;
173: b3 += bm4;
174: b4 += bm4;
176: /* RAB[:,col] = R*C[:,col] */
177: colrm = col*rm;
178: for (i=0; i<rm; i++) { /* over rows of R in those columns */
179: r1 = r2 = r3 = r4 = 0.0;
180: n = r->i[i+1] - r->i[i];
181: rj = r->j + r->i[i];
182: ra = r->a + r->i[i];
183: for (j=0; j<n; j++) {
184: r1 += (*ra)*c[*rj];
185: r2 += (*ra)*c2[*rj];
186: r3 += (*ra)*c3[*rj];
187: r4 += (*ra++)*c4[*rj++];
188: }
189: d[colrm + i] = r1;
190: d[colrm + rm + i] = r2;
191: d[colrm + rm2 + i] = r3;
192: d[colrm + rm3 + i] = r4;
193: }
194: }
195: for (; col<cn; col++) { /* over extra columns of C */
196: for (i=0; i<am; i++) { /* over rows of A in those columns */
197: r1 = 0.0;
198: n = a->i[i+1] - a->i[i];
199: aj = a->j + a->i[i];
200: aa = a->a + a->i[i];
201: for (j=0; j<n; j++) {
202: r1 += (*aa++)*b1[*aj++];
203: }
204: c[i] = r1;
205: }
206: b1 += bm;
208: for (i=0; i<rm; i++) { /* over rows of R in those columns */
209: r1 = 0.0;
210: n = r->i[i+1] - r->i[i];
211: rj = r->j + r->i[i];
212: ra = r->a + r->i[i];
213: for (j=0; j<n; j++) {
214: r1 += (*ra++)*c[*rj++];
215: }
216: d[col*rm + i] = r1;
217: }
218: }
219: PetscLogFlops(cn*2.0*(a->nz + r->nz));
221: MatDenseRestoreArray(B,&b);
222: MatDenseRestoreArray(RAB,&d);
223: MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
224: MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
225: return(0);
226: }
228: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
229: {
230: PetscErrorCode ierr;
231: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
232: Mat_RARt *rart=c->rart;
233: MatTransposeColoring matcoloring;
234: Mat Rt,RARt;
237: /* Get dense Rt by Apply MatTransposeColoring to R */
238: matcoloring = rart->matcoloring;
239: Rt = rart->Rt;
240: MatTransColoringApplySpToDen(matcoloring,R,Rt);
242: /* Get dense RARt = R*A*Rt -- dominates! */
243: RARt = rart->RARt;
244: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);
246: /* Recover C from C_dense */
247: MatTransColoringApplyDenToSp(matcoloring,RARt,C);
248: return(0);
249: }
251: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C)
252: {
253: PetscErrorCode ierr;
254: Mat ARt,RARt;
255: Mat_SeqAIJ *c;
256: Mat_RARt *rart;
259: /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */
260: MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);
261: MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);
262: *C = RARt;
263: RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
265: PetscNew(&rart);
266: c = (Mat_SeqAIJ*)(*C)->data;
267: c->rart = rart;
268: rart->ARt = ARt;
269: rart->destroy = RARt->ops->destroy;
270: RARt->ops->destroy = MatDestroy_SeqAIJ_RARt;
271: #if defined(PETSC_USE_INFO)
272: PetscInfo(*C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
273: #endif
274: return(0);
275: }
277: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
278: {
279: PetscErrorCode ierr;
280: Mat_SeqAIJ *c=(Mat_SeqAIJ*)C->data;
281: Mat_RARt *rart=c->rart;
282: Mat ARt=rart->ARt;
283:
285: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt); /* dominate! */
286: MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);
287: return(0);
288: }
290: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat *C)
291: {
292: PetscErrorCode ierr;
293: Mat Rt;
294: Mat_SeqAIJ *c;
295: Mat_RARt *rart;
298: MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
299: MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);
301: PetscNew(&rart);
302: rart->Rt = Rt;
303: c = (Mat_SeqAIJ*)(*C)->data;
304: c->rart = rart;
305: rart->destroy = (*C)->ops->destroy;
306: (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt;
307: (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
308: #if defined(PETSC_USE_INFO)
309: PetscInfo(*C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
310: #endif
311: return(0);
312: }
314: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
315: {
316: PetscErrorCode ierr;
317: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
318: Mat_RARt *rart = c->rart;
319: Mat Rt = rart->Rt;
322: MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);
323: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);
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 */
334: if (scall == MAT_INITIAL_MATRIX) {
335: PetscObjectOptionsBegin((PetscObject)A);
336: PetscOptionsObject->alreadyprinted = PETSC_FALSE; /* a hack to ensure the option shows in '-help' */
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: }