xref: /petsc/src/mat/impls/aij/seq/matrart.c (revision 55bea0eb881ce1ddd4964f044db496e4529865de)
1 
2 /*
3   Defines projective product routines where A is a SeqAIJ matrix
4           C = R * A * R^T
5 */
6 
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*/
10 #include <petsctime.h>
11 
12 /* #define RART_PROFILE */
13 
14 #undef __FUNCT__
15 #define __FUNCT__ "MatDestroy_SeqAIJ_RARt"
16 PetscErrorCode MatDestroy_SeqAIJ_RARt(Mat A)
17 {
18   PetscErrorCode ierr;
19   Mat_SeqAIJ     *a    = (Mat_SeqAIJ*)A->data;
20   Mat_RARt       *rart = a->rart;
21 
22   PetscFunctionBegin;
23   ierr = MatTransposeColoringDestroy(&rart->matcoloring);CHKERRQ(ierr);
24   ierr = MatDestroy(&rart->Rt);CHKERRQ(ierr);
25   ierr = MatDestroy(&rart->RARt);CHKERRQ(ierr);
26   ierr = MatDestroy(&rart->ARt);CHKERRQ(ierr);
27   ierr = PetscFree(rart->work);CHKERRQ(ierr);
28 
29   A->ops->destroy = rart->destroy;
30   if (A->ops->destroy) {
31     ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
32   }
33   ierr = PetscFree(rart);CHKERRQ(ierr);
34   PetscFunctionReturn(0);
35 }
36 
37 #undef __FUNCT__
38 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart"
39 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat *C)
40 {
41   PetscErrorCode       ierr;
42   Mat                  P;
43   PetscInt             *rti,*rtj;
44   Mat_RARt             *rart;
45   MatTransposeColoring matcoloring;
46   ISColoring           iscoloring;
47   Mat                  Rt_dense,RARt_dense;
48 #if defined(RART_PROFILE)
49   PetscLogDouble       GColor=0.0,MCCreate=0.0,MDenCreate=0.0,t0,tf,etime=0.0;
50 #endif
51   Mat_SeqAIJ           *c;
52 
53   PetscFunctionBegin;
54 #if defined(RART_PROFILE)
55   ierr = PetscTime(&t0);CHKERRQ(ierr);
56 #endif
57   /* create symbolic P=Rt */
58   ierr = MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr);
59   ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);CHKERRQ(ierr);
60 
61   /* get symbolic C=Pt*A*P */
62   ierr           = MatPtAPSymbolic_SeqAIJ_SeqAIJ(A,P,fill,C);CHKERRQ(ierr);
63   (*C)->rmap->bs = R->rmap->bs;
64   (*C)->cmap->bs = R->rmap->bs;
65   (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;
66 
67   /* create a supporting struct */
68   ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr);
69   c       = (Mat_SeqAIJ*)(*C)->data;
70   c->rart = rart;
71 #if defined(RART_PROFILE)
72   ierr   = PetscTime(&tf);CHKERRQ(ierr);
73   etime += tf - t0;
74 #endif
75 
76   /* ------ Use coloring ---------- */
77   /* Create MatTransposeColoring from symbolic C=R*A*R^T */
78 #if defined(RART_PROFILE)
79   ierr    = PetscTime(&t0);CHKERRQ(ierr);
80 #endif
81   ierr    = MatGetColoring(*C,MATCOLORINGLF,&iscoloring);CHKERRQ(ierr);
82 #if defined(RART_PROFILE)
83   ierr    = PetscTime(&tf);CHKERRQ(ierr);
84   GColor += tf - t0;
85   ierr = PetscTime(&t0);CHKERRQ(ierr);
86 #endif
87   ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr);
88 
89   rart->matcoloring = matcoloring;
90 
91   ierr      = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
92 #if defined(RART_PROFILE)
93   ierr      = PetscTime(&tf);CHKERRQ(ierr);
94   MCCreate += tf - t0;
95   ierr = PetscTime(&t0);CHKERRQ(ierr);
96 #endif
97 
98   /* Create Rt_dense */
99   ierr = MatCreate(PETSC_COMM_SELF,&Rt_dense);CHKERRQ(ierr);
100   ierr = MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr);
101   ierr = MatSetType(Rt_dense,MATSEQDENSE);CHKERRQ(ierr);
102   ierr = MatSeqDenseSetPreallocation(Rt_dense,NULL);CHKERRQ(ierr);
103 
104   Rt_dense->assembled = PETSC_TRUE;
105   rart->Rt            = Rt_dense;
106 
107   /* Create RARt_dense = R*A*Rt_dense */
108   ierr = MatCreate(PETSC_COMM_SELF,&RARt_dense);CHKERRQ(ierr);
109   ierr = MatSetSizes(RARt_dense,(*C)->rmap->n,matcoloring->ncolors,(*C)->rmap->n,matcoloring->ncolors);CHKERRQ(ierr);
110   ierr = MatSetType(RARt_dense,MATSEQDENSE);CHKERRQ(ierr);
111   ierr = MatSeqDenseSetPreallocation(RARt_dense,NULL);CHKERRQ(ierr);
112 
113   rart->RARt = RARt_dense;
114 
115   /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
116   ierr = PetscMalloc(A->rmap->n*4*sizeof(PetscScalar),&rart->work);CHKERRQ(ierr);
117 
118 #if defined(RART_PROFILE)
119   ierr        = PetscTime(&tf);CHKERRQ(ierr);
120   MDenCreate += tf - t0;
121 #endif
122 
123   rart->destroy      = (*C)->ops->destroy;
124   (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt;
125 
126   /* clean up */
127   ierr = MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr);
128   ierr = MatDestroy(&P);CHKERRQ(ierr);
129 
130 #if defined(PETSC_USE_INFO)
131   {
132     PetscReal density= (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
133     ierr = PetscInfo6(*C,"RARt_den %D %D; Rt_den %D %D, (RARt->nz %D)/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,Rt_dense->rmap->n,Rt_dense->cmap->n,c->nz,density);CHKERRQ(ierr);
134 #if defined(RART_PROFILE)
135     ierr = PetscInfo5(*C,"Sym = GetColor %g + MColorCreate %g + MDenCreate %g + other %g = %g\n",GColor,MCCreate,MDenCreate,etime,GColor+MCCreate+MDenCreate+etime);CHKERRQ(ierr);
136 #endif
137   }
138 #endif
139   PetscFunctionReturn(0);
140 }
141 
142 /*
143  RAB = R * A * B, R and A in seqaij format, B in dense format;
144 */
145 #undef __FUNCT__
146 #define __FUNCT__ "MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense"
147 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
148 {
149   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
150   PetscErrorCode ierr;
151   PetscScalar    *b,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
152   MatScalar      *aa,*ra;
153   PetscInt       cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
154   PetscInt       am2=2*am,am3=3*am,bm4=4*bm;
155   PetscScalar    *d,*c,*c2,*c3,*c4;
156   PetscInt       *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
157   PetscInt       rm2=2*rm,rm3=3*rm,colrm;
158 #if defined(RART_PROFILE)
159   PetscLogDouble t0,tf,Mult_sp_den1=0.0,Mult_sp_den2=0.0;
160 #endif
161 
162   PetscFunctionBegin;
163   if (!dm || !dn) PetscFunctionReturn(0);
164   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);
165   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);
166   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);
167   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);
168 
169   { /*
170      This approach is not as good as original ones (will be removed later), but it reveals that
171      AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c
172      */
173     PetscBool via_matmatmult=PETSC_FALSE;
174     ierr = PetscOptionsGetBool(NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);CHKERRQ(ierr);
175     if (via_matmatmult) {
176       Mat AB_den;
177       ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);CHKERRQ(ierr);
178 
179 #if defined(RART_PROFILE)
180       ierr = PetscTime(&t0);CHKERRQ(ierr);
181 #endif
182       ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);CHKERRQ(ierr);
183 #if defined(RART_PROFILE)
184       ierr  = PetscTime(&tf);CHKERRQ(ierr);
185       Mult_sp_den1 += tf - t0;
186       ierr = PetscTime(&t0);CHKERRQ(ierr);
187 #endif
188       ierr = MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);CHKERRQ(ierr);
189 #if defined(RART_PROFILE)
190       ierr  = PetscTime(&tf);CHKERRQ(ierr);
191       Mult_sp_den2 += tf - t0;
192       printf(" MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense...time = %g + %g = %g; ncolor: %d\n",Mult_sp_den1, Mult_sp_den2, Mult_sp_den1+ Mult_sp_den2,B->cmap->n);
193 #endif
194       ierr = MatDestroy(&AB_den);CHKERRQ(ierr);
195       PetscFunctionReturn(0);
196     }
197   }
198 
199   ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr);
200   ierr = MatDenseGetArray(RAB,&d);CHKERRQ(ierr);
201   b1   = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
202   c    = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
203   for (col=0; col<cn-4; col += 4) {  /* over columns of C */
204     for (i=0; i<am; i++) {        /* over rows of A in those columns */
205       r1 = r2 = r3 = r4 = 0.0;
206       n  = ai[i+1] - ai[i];
207       aj = a->j + ai[i];
208       aa = a->a + ai[i];
209       for (j=0; j<n; j++) {
210         r1 += (*aa)*b1[*aj];
211         r2 += (*aa)*b2[*aj];
212         r3 += (*aa)*b3[*aj];
213         r4 += (*aa++)*b4[*aj++];
214       }
215       c[i]       = r1;
216       c[am  + i] = r2;
217       c[am2 + i] = r3;
218       c[am3 + i] = r4;
219     }
220     b1 += bm4;
221     b2 += bm4;
222     b3 += bm4;
223     b4 += bm4;
224 
225     /* RAB[:,col] = R*C[:,col] */
226     colrm = col*rm;
227     for (i=0; i<rm; i++) {        /* over rows of R in those columns */
228       r1 = r2 = r3 = r4 = 0.0;
229       n  = r->i[i+1] - r->i[i];
230       rj = r->j + r->i[i];
231       ra = r->a + r->i[i];
232       for (j=0; j<n; j++) {
233         r1 += (*ra)*c[*rj];
234         r2 += (*ra)*c2[*rj];
235         r3 += (*ra)*c3[*rj];
236         r4 += (*ra++)*c4[*rj++];
237       }
238       d[colrm + i]       = r1;
239       d[colrm + rm + i]  = r2;
240       d[colrm + rm2 + i] = r3;
241       d[colrm + rm3 + i] = r4;
242     }
243   }
244   for (; col<cn; col++) {     /* over extra columns of C */
245     for (i=0; i<am; i++) {  /* over rows of A in those columns */
246       r1 = 0.0;
247       n  = a->i[i+1] - a->i[i];
248       aj = a->j + a->i[i];
249       aa = a->a + a->i[i];
250       for (j=0; j<n; j++) {
251         r1 += (*aa++)*b1[*aj++];
252       }
253       c[i] = r1;
254     }
255     b1 += bm;
256 
257     for (i=0; i<rm; i++) {  /* over rows of R in those columns */
258       r1 = 0.0;
259       n  = r->i[i+1] - r->i[i];
260       rj = r->j + r->i[i];
261       ra = r->a + r->i[i];
262       for (j=0; j<n; j++) {
263         r1 += (*ra++)*c[*rj++];
264       }
265       d[col*rm + i] = r1;
266     }
267   }
268   ierr = PetscLogFlops(cn*2.0*(a->nz + r->nz));CHKERRQ(ierr);
269 
270   ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr);
271   ierr = MatDenseRestoreArray(RAB,&d);CHKERRQ(ierr);
272   ierr = MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
273   ierr = MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
274   PetscFunctionReturn(0);
275 }
276 
277 #undef __FUNCT__
278 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart"
279 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
280 {
281   PetscErrorCode       ierr;
282   Mat_SeqAIJ           *c = (Mat_SeqAIJ*)C->data;
283   Mat_RARt             *rart=c->rart;
284   MatTransposeColoring matcoloring;
285   Mat                  Rt,RARt;
286 #if defined(RART_PROFILE)
287   PetscLogDouble       Mult_sp_den=0.0,app1=0.0,app2=0.0,t0,tf;
288 #endif
289 
290   PetscFunctionBegin;
291   /* Get dense Rt by Apply MatTransposeColoring to R */
292   matcoloring = rart->matcoloring;
293   Rt          = rart->Rt;
294 
295 #if defined(RART_PROFILE)
296   ierr  = PetscTime(&t0);CHKERRQ(ierr);
297 #endif
298   ierr  = MatTransColoringApplySpToDen(matcoloring,R,Rt);CHKERRQ(ierr);
299 #if defined(RART_PROFILE)
300   ierr  = PetscTime(&tf);CHKERRQ(ierr);
301   app1 += tf - t0;
302 #endif
303 
304   /* Get dense RARt = R*A*Rt -- dominates! */
305 #if defined(RART_PROFILE)
306   ierr = PetscTime(&t0);CHKERRQ(ierr);
307 #endif
308   RARt = rart->RARt;
309   ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);CHKERRQ(ierr);
310 #if defined(RART_PROFILE)
311   ierr = PetscTime(&tf);CHKERRQ(ierr);
312   Mult_sp_den += tf - t0;
313 #endif
314 
315   /* Recover C from C_dense */
316 #if defined(RART_PROFILE)
317   ierr = PetscTime(&t0);CHKERRQ(ierr);
318 #endif
319   ierr = MatTransColoringApplyDenToSp(matcoloring,RARt,C);CHKERRQ(ierr);
320 #if defined(RART_PROFILE)
321   ierr = PetscTime(&tf);CHKERRQ(ierr);
322   app2 += tf - t0;
323 #endif
324 
325 #if defined(RART_PROFILE)
326   printf("MatRARtNumeric_SeqAIJ_SeqAIJ: ColorApp %g + %g + Mult_sp_den %g  = %g\n",app1,app2,Mult_sp_den,app1+app2+Mult_sp_den);CHKERRQ(ierr);
327 #endif
328   PetscFunctionReturn(0);
329 }
330 
331 #undef __FUNCT__
332 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult"
333 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C)
334 {
335   PetscErrorCode  ierr;
336   Mat             ARt,RARt;
337   Mat_SeqAIJ     *c;
338   Mat_RARt       *rart;
339 
340   PetscFunctionBegin;
341   /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */
342   ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);CHKERRQ(ierr);
343   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);CHKERRQ(ierr);
344   *C                     = RARt;
345   RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;
346 
347   ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr);
348   c         = (Mat_SeqAIJ*)(*C)->data;
349   c->rart   = rart;
350   rart->ARt = ARt;
351   rart->destroy      = RARt->ops->destroy;
352   RARt->ops->destroy = MatDestroy_SeqAIJ_RARt;
353   PetscFunctionReturn(0);
354 }
355 
356 #undef __FUNCT__
357 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult"
358 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
359 {
360   PetscErrorCode  ierr;
361   Mat_SeqAIJ      *c=(Mat_SeqAIJ*)C->data;
362   Mat_RARt        *rart=c->rart;
363   Mat             ARt=rart->ARt;
364 #if defined(RART_PROFILE)
365   PetscLogDouble  t0,t1,t2;
366 #endif
367 
368   PetscFunctionBegin;
369 #if defined(RART_PROFILE)
370     ierr = PetscTime(&t0);CHKERRQ(ierr);
371 #endif
372     ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt);CHKERRQ(ierr); /* dominate! */
373 #if defined(RART_PROFILE)
374     ierr = PetscTime(&t1);CHKERRQ(ierr);
375 #endif
376     ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);CHKERRQ(ierr);
377 #if defined(RART_PROFILE)
378     ierr = PetscTime(&t2);CHKERRQ(ierr);
379     printf(" matrart_color_art_num = %g + %g = %g\n",t1-t0,t2-t1,t2-t0);
380 #endif
381 
382 #if defined(PETSC_USE_INFO)
383     ierr = PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");CHKERRQ(ierr);
384 #endif
385   PetscFunctionReturn(0);
386 }
387 
388 #undef __FUNCT__
389 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ"
390 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat *C)
391 {
392   PetscErrorCode  ierr;
393   Mat             Rt;
394   Mat_SeqAIJ      *c;
395   Mat_RARt        *rart;
396 
397   PetscFunctionBegin;
398   ierr = MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr);
399   ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);CHKERRQ(ierr);
400 
401   ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr);
402   rart->Rt = Rt;
403   c        = (Mat_SeqAIJ*)(*C)->data;
404   c->rart  = rart;
405   rart->destroy          = (*C)->ops->destroy;
406   (*C)->ops->destroy     = MatDestroy_SeqAIJ_RARt;
407   (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
408   PetscFunctionReturn(0);
409 }
410 
411 #undef __FUNCT__
412 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ"
413 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
414 {
415   PetscErrorCode  ierr;
416   Mat_SeqAIJ      *c = (Mat_SeqAIJ*)C->data;
417   Mat_RARt        *rart = c->rart;
418   Mat             Rt = rart->Rt;
419 
420   PetscFunctionBegin;
421   ierr = MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);CHKERRQ(ierr);
422   ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);CHKERRQ(ierr);
423 #if defined(PETSC_USE_INFO)
424   ierr = PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");CHKERRQ(ierr);
425 #endif
426   PetscFunctionReturn(0);
427 }
428 
429 #undef __FUNCT__
430 #define __FUNCT__ "MatRARt_SeqAIJ_SeqAIJ"
431 PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
432 {
433   PetscErrorCode ierr;
434   const char     *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
435   PetscInt       alg=0; /* set default algorithm */
436 
437   PetscFunctionBegin;
438   if (scall == MAT_INITIAL_MATRIX) {
439     /* ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); -- prefix ? */
440     ierr = PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);CHKERRQ(ierr);
441     /* ierr = PetscOptionsEnd();CHKERRQ(ierr); */
442     ierr = PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr);
443     switch (alg) {
444     case 1:
445       /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
446       ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);CHKERRQ(ierr);
447       break;
448     case 2:
449       /* via coloring_rart: apply coloring C = R*A*R^T                          */
450       ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);CHKERRQ(ierr);
451       break;
452     default:
453       /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
454       ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);CHKERRQ(ierr);
455       break;
456     }
457     ierr = PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr);
458   }
459 
460   ierr = PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr);
461   ierr = (*(*C)->ops->rartnumeric)(A,R,*C);CHKERRQ(ierr);
462   ierr = PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr);
463   PetscFunctionReturn(0);
464 }
465