Actual source code: matrart.c
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: {
13: Mat_RARt *rart = (Mat_RARt*)data;
15: MatTransposeColoringDestroy(&rart->matcoloring);
16: MatDestroy(&rart->Rt);
17: MatDestroy(&rart->RARt);
18: MatDestroy(&rart->ARt);
19: PetscFree(rart->work);
20: if (rart->destroy) {
21: (*rart->destroy)(rart->data);
22: }
23: PetscFree(rart);
24: return 0;
25: }
27: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat C)
28: {
29: Mat P;
30: PetscInt *rti,*rtj;
31: Mat_RARt *rart;
32: MatColoring coloring;
33: MatTransposeColoring matcoloring;
34: ISColoring iscoloring;
35: Mat Rt_dense,RARt_dense;
37: MatCheckProduct(C,4);
39: /* create symbolic P=Rt */
40: MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
41: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);
43: /* get symbolic C=Pt*A*P */
44: MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
45: MatSetBlockSizes(C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
46: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;
48: /* create a supporting struct */
49: PetscNew(&rart);
50: C->product->data = rart;
51: C->product->destroy = MatDestroy_SeqAIJ_RARt;
53: /* ------ Use coloring ---------- */
54: /* inode causes memory problem */
55: MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);
57: /* Create MatTransposeColoring from symbolic C=R*A*R^T */
58: MatColoringCreate(C,&coloring);
59: MatColoringSetDistance(coloring,2);
60: MatColoringSetType(coloring,MATCOLORINGSL);
61: MatColoringSetFromOptions(coloring);
62: MatColoringApply(coloring,&iscoloring);
63: MatColoringDestroy(&coloring);
64: MatTransposeColoringCreate(C,iscoloring,&matcoloring);
66: rart->matcoloring = matcoloring;
67: ISColoringDestroy(&iscoloring);
69: /* Create Rt_dense */
70: MatCreate(PETSC_COMM_SELF,&Rt_dense);
71: MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
72: MatSetType(Rt_dense,MATSEQDENSE);
73: MatSeqDenseSetPreallocation(Rt_dense,NULL);
75: Rt_dense->assembled = PETSC_TRUE;
76: rart->Rt = Rt_dense;
78: /* Create RARt_dense = R*A*Rt_dense */
79: MatCreate(PETSC_COMM_SELF,&RARt_dense);
80: MatSetSizes(RARt_dense,C->rmap->n,matcoloring->ncolors,C->rmap->n,matcoloring->ncolors);
81: MatSetType(RARt_dense,MATSEQDENSE);
82: MatSeqDenseSetPreallocation(RARt_dense,NULL);
84: rart->RARt = RARt_dense;
86: /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
87: PetscMalloc1(A->rmap->n*4,&rart->work);
89: /* clean up */
90: MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
91: MatDestroy(&P);
93: #if defined(PETSC_USE_INFO)
94: {
95: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
96: PetscReal density = (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
97: PetscInfo(C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
98: PetscInfo(C,"RARt_den %" PetscInt_FMT " %" PetscInt_FMT "; Rt %" PetscInt_FMT " %" PetscInt_FMT " (RARt->nz %" PetscInt_FMT ")/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,R->cmap->n,R->rmap->n,c->nz,density);
99: }
100: #endif
101: return 0;
102: }
104: /*
105: RAB = R * A * B, R and A in seqaij format, B in dense format;
106: */
107: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
108: {
109: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
110: PetscScalar r1,r2,r3,r4;
111: const PetscScalar *b,*b1,*b2,*b3,*b4;
112: MatScalar *aa,*ra;
113: PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
114: PetscInt am2=2*am,am3=3*am,bm4=4*bm;
115: PetscScalar *d,*c,*c2,*c3,*c4;
116: PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
117: PetscInt rm2=2*rm,rm3=3*rm,colrm;
119: if (!dm || !dn) return 0;
125: { /*
126: This approach is not as good as original ones (will be removed later), but it reveals that
127: AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/tutorials/ex56.c
128: */
129: PetscBool via_matmatmult=PETSC_FALSE;
130: PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
131: if (via_matmatmult) {
132: Mat AB_den = NULL;
133: MatCreate(PetscObjectComm((PetscObject)A),&AB_den);
134: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,AB_den);
135: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
136: MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
137: MatDestroy(&AB_den);
138: return 0;
139: }
140: }
142: MatDenseGetArrayRead(B,&b);
143: MatDenseGetArray(RAB,&d);
144: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
145: c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
146: for (col=0; col<cn-4; col += 4) { /* over columns of C */
147: for (i=0; i<am; i++) { /* over rows of A in those columns */
148: r1 = r2 = r3 = r4 = 0.0;
149: n = ai[i+1] - ai[i];
150: aj = a->j + ai[i];
151: aa = a->a + ai[i];
152: for (j=0; j<n; j++) {
153: r1 += (*aa)*b1[*aj];
154: r2 += (*aa)*b2[*aj];
155: r3 += (*aa)*b3[*aj];
156: r4 += (*aa++)*b4[*aj++];
157: }
158: c[i] = r1;
159: c[am + i] = r2;
160: c[am2 + i] = r3;
161: c[am3 + i] = r4;
162: }
163: b1 += bm4;
164: b2 += bm4;
165: b3 += bm4;
166: b4 += bm4;
168: /* RAB[:,col] = R*C[:,col] */
169: colrm = col*rm;
170: for (i=0; i<rm; i++) { /* over rows of R in those columns */
171: r1 = r2 = r3 = r4 = 0.0;
172: n = r->i[i+1] - r->i[i];
173: rj = r->j + r->i[i];
174: ra = r->a + r->i[i];
175: for (j=0; j<n; j++) {
176: r1 += (*ra)*c[*rj];
177: r2 += (*ra)*c2[*rj];
178: r3 += (*ra)*c3[*rj];
179: r4 += (*ra++)*c4[*rj++];
180: }
181: d[colrm + i] = r1;
182: d[colrm + rm + i] = r2;
183: d[colrm + rm2 + i] = r3;
184: d[colrm + rm3 + i] = r4;
185: }
186: }
187: for (; col<cn; col++) { /* over extra columns of C */
188: for (i=0; i<am; i++) { /* over rows of A in those columns */
189: r1 = 0.0;
190: n = a->i[i+1] - a->i[i];
191: aj = a->j + a->i[i];
192: aa = a->a + a->i[i];
193: for (j=0; j<n; j++) {
194: r1 += (*aa++)*b1[*aj++];
195: }
196: c[i] = r1;
197: }
198: b1 += bm;
200: for (i=0; i<rm; i++) { /* over rows of R in those columns */
201: r1 = 0.0;
202: n = r->i[i+1] - r->i[i];
203: rj = r->j + r->i[i];
204: ra = r->a + r->i[i];
205: for (j=0; j<n; j++) {
206: r1 += (*ra++)*c[*rj++];
207: }
208: d[col*rm + i] = r1;
209: }
210: }
211: PetscLogFlops(cn*2.0*(a->nz + r->nz));
213: MatDenseRestoreArrayRead(B,&b);
214: MatDenseRestoreArray(RAB,&d);
215: MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
216: MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
217: return 0;
218: }
220: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
221: {
222: Mat_RARt *rart;
223: MatTransposeColoring matcoloring;
224: Mat Rt,RARt;
226: MatCheckProduct(C,3);
228: rart = (Mat_RARt*)C->product->data;
230: /* Get dense Rt by Apply MatTransposeColoring to R */
231: matcoloring = rart->matcoloring;
232: Rt = rart->Rt;
233: MatTransColoringApplySpToDen(matcoloring,R,Rt);
235: /* Get dense RARt = R*A*Rt -- dominates! */
236: RARt = rart->RARt;
237: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);
239: /* Recover C from C_dense */
240: MatTransColoringApplyDenToSp(matcoloring,RARt,C);
241: return 0;
242: }
244: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat C)
245: {
246: Mat ARt;
247: Mat_RARt *rart;
248: char *alg;
250: MatCheckProduct(C,4);
252: /* create symbolic ARt = A*R^T */
253: MatProductCreate(A,R,NULL,&ARt);
254: MatProductSetType(ARt,MATPRODUCT_ABt);
255: MatProductSetAlgorithm(ARt,"sorted");
256: MatProductSetFill(ARt,fill);
257: MatProductSetFromOptions(ARt);
258: MatProductSymbolic(ARt);
260: /* compute symbolic C = R*ARt */
261: /* set algorithm for C = R*ARt */
262: PetscStrallocpy(C->product->alg,&alg);
263: MatProductSetAlgorithm(C,"sorted");
264: MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,C);
265: /* resume original algorithm for C */
266: MatProductSetAlgorithm(C,alg);
267: PetscFree(alg);
269: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
271: PetscNew(&rart);
272: rart->ARt = ARt;
273: C->product->data = rart;
274: C->product->destroy = MatDestroy_SeqAIJ_RARt;
275: PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
276: return 0;
277: }
279: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
280: {
281: Mat_RARt *rart;
283: MatCheckProduct(C,3);
285: rart = (Mat_RARt*)C->product->data;
286: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,rart->ARt); /* dominate! */
287: MatMatMultNumeric_SeqAIJ_SeqAIJ(R,rart->ARt,C);
288: return 0;
289: }
291: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat C)
292: {
293: Mat Rt;
294: Mat_RARt *rart;
296: MatCheckProduct(C,4);
298: MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
299: MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);
301: PetscNew(&rart);
302: rart->data = C->product->data;
303: rart->destroy = C->product->destroy;
304: rart->Rt = Rt;
305: C->product->data = rart;
306: C->product->destroy = MatDestroy_SeqAIJ_RARt;
307: C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
308: PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
309: return 0;
310: }
312: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
313: {
314: Mat_RARt *rart;
316: MatCheckProduct(C,3);
318: rart = (Mat_RARt*)C->product->data;
319: MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&rart->Rt);
320: /* MatMatMatMultSymbolic used a different data */
321: C->product->data = rart->data;
322: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,rart->Rt,C);
323: C->product->data = rart;
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 */
333: if (scall == MAT_INITIAL_MATRIX) {
334: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatRARt","Mat");
335: PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
336: PetscOptionsEnd();
338: PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
339: MatCreate(PETSC_COMM_SELF,C);
340: switch (alg) {
341: case 1:
342: /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
343: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,*C);
344: break;
345: case 2:
346: /* via coloring_rart: apply coloring C = R*A*R^T */
347: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,*C);
348: break;
349: default:
350: /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
351: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,*C);
352: break;
353: }
354: PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
355: }
357: PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
358: ((*C)->ops->rartnumeric)(A,R,*C);
359: PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
360: return 0;
361: }
363: /* ------------------------------------------------------------- */
364: PetscErrorCode MatProductSymbolic_RARt_SeqAIJ_SeqAIJ(Mat C)
365: {
366: Mat_Product *product = C->product;
367: Mat A=product->A,R=product->B;
368: MatProductAlgorithm alg=product->alg;
369: PetscReal fill=product->fill;
370: PetscBool flg;
372: PetscStrcmp(alg,"r*a*rt",&flg);
373: if (flg) {
374: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
375: goto next;
376: }
378: PetscStrcmp(alg,"r*art",&flg);
379: if (flg) {
380: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
381: goto next;
382: }
384: PetscStrcmp(alg,"coloring_rart",&flg);
385: if (flg) {
386: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
387: goto next;
388: }
390: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductAlgorithm is not supported");
392: next:
393: C->ops->productnumeric = MatProductNumeric_RARt;
394: return 0;
395: }