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