Actual source code: matrart.c
petsc-3.6.1 2015-08-06
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> /*I "petscmat.h" I*/
13: PetscErrorCode MatDestroy_SeqAIJ_RARt(Mat A)
14: {
16: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
17: Mat_RARt *rart = a->rart;
20: MatTransposeColoringDestroy(&rart->matcoloring);
21: MatDestroy(&rart->Rt);
22: MatDestroy(&rart->RARt);
23: MatDestroy(&rart->ARt);
24: PetscFree(rart->work);
26: A->ops->destroy = rart->destroy;
27: if (A->ops->destroy) {
28: (*A->ops->destroy)(A);
29: }
30: PetscFree(rart);
31: return(0);
32: }
36: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat *C)
37: {
38: PetscErrorCode ierr;
39: Mat P;
40: PetscInt *rti,*rtj;
41: Mat_RARt *rart;
42: MatColoring coloring;
43: MatTransposeColoring matcoloring;
44: ISColoring iscoloring;
45: Mat Rt_dense,RARt_dense;
46: Mat_SeqAIJ *c;
49: /* create symbolic P=Rt */
50: MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
51: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);
53: /* get symbolic C=Pt*A*P */
54: MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
55: MatSetBlockSizes(*C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
56: (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;
58: /* create a supporting struct */
59: PetscNew(&rart);
60: c = (Mat_SeqAIJ*)(*C)->data;
61: c->rart = rart;
63: /* ------ Use coloring ---------- */
64: /* inode causes memory problem, don't know why */
65: if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'");
67: /* Create MatTransposeColoring from symbolic C=R*A*R^T */
68: MatColoringCreate(*C,&coloring);
69: MatColoringSetDistance(coloring,2);
70: MatColoringSetType(coloring,MATCOLORINGSL);
71: MatColoringSetFromOptions(coloring);
72: MatColoringApply(coloring,&iscoloring);
73: MatColoringDestroy(&coloring);
74: MatTransposeColoringCreate(*C,iscoloring,&matcoloring);
76: rart->matcoloring = matcoloring;
77: ISColoringDestroy(&iscoloring);
79: /* Create Rt_dense */
80: MatCreate(PETSC_COMM_SELF,&Rt_dense);
81: MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
82: MatSetType(Rt_dense,MATSEQDENSE);
83: MatSeqDenseSetPreallocation(Rt_dense,NULL);
85: Rt_dense->assembled = PETSC_TRUE;
86: rart->Rt = Rt_dense;
88: /* Create RARt_dense = R*A*Rt_dense */
89: MatCreate(PETSC_COMM_SELF,&RARt_dense);
90: MatSetSizes(RARt_dense,(*C)->rmap->n,matcoloring->ncolors,(*C)->rmap->n,matcoloring->ncolors);
91: MatSetType(RARt_dense,MATSEQDENSE);
92: MatSeqDenseSetPreallocation(RARt_dense,NULL);
94: rart->RARt = RARt_dense;
96: /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
97: PetscMalloc1(A->rmap->n*4,&rart->work);
99: rart->destroy = (*C)->ops->destroy;
100: (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt;
102: /* clean up */
103: MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
104: MatDestroy(&P);
106: #if defined(PETSC_USE_INFO)
107: {
108: PetscReal density= (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
109: PetscInfo(*C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
110: 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);
111: }
112: #endif
113: return(0);
114: }
116: /*
117: RAB = R * A * B, R and A in seqaij format, B in dense format;
118: */
121: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
122: {
123: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
125: PetscScalar *b,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
126: MatScalar *aa,*ra;
127: PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
128: PetscInt am2=2*am,am3=3*am,bm4=4*bm;
129: PetscScalar *d,*c,*c2,*c3,*c4;
130: PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
131: PetscInt rm2=2*rm,rm3=3*rm,colrm;
134: if (!dm || !dn) return(0);
135: 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);
136: 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);
137: 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);
138: 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);
140: { /*
141: This approach is not as good as original ones (will be removed later), but it reveals that
142: AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c
143: */
144: PetscBool via_matmatmult=PETSC_FALSE;
145: PetscOptionsGetBool(NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
146: if (via_matmatmult) {
147: Mat AB_den;
148: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);
149: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
150: MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
151: MatDestroy(&AB_den);
152: return(0);
153: }
154: }
156: MatDenseGetArray(B,&b);
157: MatDenseGetArray(RAB,&d);
158: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
159: c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
160: for (col=0; col<cn-4; col += 4) { /* over columns of C */
161: for (i=0; i<am; i++) { /* over rows of A in those columns */
162: r1 = r2 = r3 = r4 = 0.0;
163: n = ai[i+1] - ai[i];
164: aj = a->j + ai[i];
165: aa = a->a + ai[i];
166: for (j=0; j<n; j++) {
167: r1 += (*aa)*b1[*aj];
168: r2 += (*aa)*b2[*aj];
169: r3 += (*aa)*b3[*aj];
170: r4 += (*aa++)*b4[*aj++];
171: }
172: c[i] = r1;
173: c[am + i] = r2;
174: c[am2 + i] = r3;
175: c[am3 + i] = r4;
176: }
177: b1 += bm4;
178: b2 += bm4;
179: b3 += bm4;
180: b4 += bm4;
182: /* RAB[:,col] = R*C[:,col] */
183: colrm = col*rm;
184: for (i=0; i<rm; i++) { /* over rows of R in those columns */
185: r1 = r2 = r3 = r4 = 0.0;
186: n = r->i[i+1] - r->i[i];
187: rj = r->j + r->i[i];
188: ra = r->a + r->i[i];
189: for (j=0; j<n; j++) {
190: r1 += (*ra)*c[*rj];
191: r2 += (*ra)*c2[*rj];
192: r3 += (*ra)*c3[*rj];
193: r4 += (*ra++)*c4[*rj++];
194: }
195: d[colrm + i] = r1;
196: d[colrm + rm + i] = r2;
197: d[colrm + rm2 + i] = r3;
198: d[colrm + rm3 + i] = r4;
199: }
200: }
201: for (; col<cn; col++) { /* over extra columns of C */
202: for (i=0; i<am; i++) { /* over rows of A in those columns */
203: r1 = 0.0;
204: n = a->i[i+1] - a->i[i];
205: aj = a->j + a->i[i];
206: aa = a->a + a->i[i];
207: for (j=0; j<n; j++) {
208: r1 += (*aa++)*b1[*aj++];
209: }
210: c[i] = r1;
211: }
212: b1 += bm;
214: for (i=0; i<rm; i++) { /* over rows of R in those columns */
215: r1 = 0.0;
216: n = r->i[i+1] - r->i[i];
217: rj = r->j + r->i[i];
218: ra = r->a + r->i[i];
219: for (j=0; j<n; j++) {
220: r1 += (*ra++)*c[*rj++];
221: }
222: d[col*rm + i] = r1;
223: }
224: }
225: PetscLogFlops(cn*2.0*(a->nz + r->nz));
227: MatDenseRestoreArray(B,&b);
228: MatDenseRestoreArray(RAB,&d);
229: MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
230: MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
231: return(0);
232: }
236: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
237: {
238: PetscErrorCode ierr;
239: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
240: Mat_RARt *rart=c->rart;
241: MatTransposeColoring matcoloring;
242: Mat Rt,RARt;
245: /* Get dense Rt by Apply MatTransposeColoring to R */
246: matcoloring = rart->matcoloring;
247: Rt = rart->Rt;
248: MatTransColoringApplySpToDen(matcoloring,R,Rt);
250: /* Get dense RARt = R*A*Rt -- dominates! */
251: RARt = rart->RARt;
252: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);
254: /* Recover C from C_dense */
255: MatTransColoringApplyDenToSp(matcoloring,RARt,C);
256: return(0);
257: }
261: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C)
262: {
263: PetscErrorCode ierr;
264: Mat ARt,RARt;
265: Mat_SeqAIJ *c;
266: Mat_RARt *rart;
269: /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */
270: MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);
271: MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);
272: *C = RARt;
273: RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
275: PetscNew(&rart);
276: c = (Mat_SeqAIJ*)(*C)->data;
277: c->rart = rart;
278: rart->ARt = ARt;
279: rart->destroy = RARt->ops->destroy;
280: RARt->ops->destroy = MatDestroy_SeqAIJ_RARt;
281: #if defined(PETSC_USE_INFO)
282: PetscInfo(*C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
283: #endif
284: return(0);
285: }
289: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
290: {
291: PetscErrorCode ierr;
292: Mat_SeqAIJ *c=(Mat_SeqAIJ*)C->data;
293: Mat_RARt *rart=c->rart;
294: Mat ARt=rart->ARt;
295:
297: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt); /* dominate! */
298: MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);
299: return(0);
300: }
304: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat *C)
305: {
306: PetscErrorCode ierr;
307: Mat Rt;
308: Mat_SeqAIJ *c;
309: Mat_RARt *rart;
312: MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
313: MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);
315: PetscNew(&rart);
316: rart->Rt = Rt;
317: c = (Mat_SeqAIJ*)(*C)->data;
318: c->rart = rart;
319: rart->destroy = (*C)->ops->destroy;
320: (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt;
321: (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
322: #if defined(PETSC_USE_INFO)
323: PetscInfo(*C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
324: #endif
325: return(0);
326: }
330: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
331: {
332: PetscErrorCode ierr;
333: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
334: Mat_RARt *rart = c->rart;
335: Mat Rt = rart->Rt;
338: MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);
339: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);
340: return(0);
341: }
345: PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
346: {
348: const char *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
349: PetscInt alg=0; /* set default algorithm */
350:
352: if (scall == MAT_INITIAL_MATRIX) {
353: PetscObjectOptionsBegin((PetscObject)A);
354: PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
355: PetscOptionsEnd();
357: PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
358: switch (alg) {
359: case 1:
360: /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
361: MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
362: break;
363: case 2:
364: /* via coloring_rart: apply coloring C = R*A*R^T */
365: MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
366: break;
367: default:
368: /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
369: MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
370: break;
371: }
372: PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
373: }
375: PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
376: (*(*C)->ops->rartnumeric)(A,R,*C);
377: PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
378: return(0);
379: }