Actual source code: matproduct.c


  2: /*
  3:     Routines for matrix products. Calling procedure:

  5:     MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
  6:     MatProductSetType(D, MATPRODUCT_AB/AtB/ABt/PtAP/RARt/ABC)
  7:     MatProductSetAlgorithm(D, alg)
  8:     MatProductSetFill(D,fill)
  9:     MatProductSetFromOptions(D)
 10:       -> MatProductSetFromOptions_Private(D)
 11:            # Check matrix global sizes
 12:            if the matrices have the same setfromoptions routine, use it
 13:            if not, try:
 14:              -> Query MatProductSetFromOptions_Atype_Btype_Ctype_C(D) from A, B and C (in order)
 15:              if found -> run the specific setup that must set the symbolic operation (these callbacks should never fail)
 16:            if callback not found or no symbolic operation set
 17:              -> Query MatProductSetFromOptions_anytype_C(D) from A, B and C (in order) (e.g, matrices may have inner matrices like MATTRANSPOSEMAT)
 18:            if dispatch found but combination still not present do
 19:              -> check if B is dense and product type AtB or AB -> if true, basic looping of dense columns
 20:              -> check if triple product (PtAP, RARt or ABC) -> if true, set the Basic routines

 22:     #  The setfromoptions calls MatProductSetFromOptions_Atype_Btype_Ctype should
 23:     #    Check matrix local sizes for mpi matrices
 24:     #    Set default algorithm
 25:     #    Get runtime option
 26:     #    Set D->ops->productsymbolic = MatProductSymbolic_productype_Atype_Btype_Ctype if found

 28:     MatProductSymbolic(D)
 29:       # Call MatProductSymbolic_productype_Atype_Btype_Ctype()
 30:         the callback must set the numeric phase D->ops->productnumeric = MatProductNumeric_productype_Atype_Btype_Ctype

 32:     MatProductNumeric(D)
 33:       # Call the numeric phase

 35:     # The symbolic phases are allowed to set extra data structures and attach those to the product
 36:     # this additional data can be reused between multiple numeric phases with the same matrices
 37:     # if not needed, call
 38:     MatProductClear(D)
 39: */

 41: #include <petsc/private/matimpl.h>

 43: const char *const MatProductTypes[] = {"UNSPECIFIED","AB","AtB","ABt","PtAP","RARt","ABC"};

 45: /* these are basic implementations relying on the old function pointers
 46:  * they are dangerous and should be removed in the future */
 47: static PetscErrorCode MatProductNumeric_PtAP_Unsafe(Mat C)
 48: {
 50:   Mat_Product    *product = C->product;
 51:   Mat            P = product->B,AP = product->Dwork;

 54:   /* AP = A*P */
 55:   MatProductNumeric(AP);
 56:   /* C = P^T*AP */
 57:   if (!C->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
 58:   (*C->ops->transposematmultnumeric)(P,AP,C);
 59:   return(0);
 60: }

 62: static PetscErrorCode MatProductSymbolic_PtAP_Unsafe(Mat C)
 63: {
 65:   Mat_Product    *product = C->product;
 66:   Mat            A=product->A,P=product->B,AP;
 67:   PetscReal      fill=product->fill;

 70:   PetscInfo2((PetscObject)C,"for A %s, P %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
 71:   /* AP = A*P */
 72:   MatProductCreate(A,P,NULL,&AP);
 73:   MatProductSetType(AP,MATPRODUCT_AB);
 74:   MatProductSetAlgorithm(AP,MATPRODUCTALGORITHM_DEFAULT);
 75:   MatProductSetFill(AP,fill);
 76:   MatProductSetFromOptions(AP);
 77:   MatProductSymbolic(AP);

 79:   /* C = P^T*AP */
 80:   MatProductSetType(C,MATPRODUCT_AtB);
 81:   MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
 82:   product->A = P;
 83:   product->B = AP;
 84:   MatProductSetFromOptions(C);
 85:   MatProductSymbolic(C);

 87:   /* resume user's original input matrix setting for A and B */
 88:   product->A     = A;
 89:   product->B     = P;
 90:   product->Dwork = AP;

 92:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
 93:   return(0);
 94: }

 96: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
 97: {
 99:   Mat_Product    *product = C->product;
100:   Mat            R=product->B,RA=product->Dwork;

103:   /* RA = R*A */
104:   MatProductNumeric(RA);
105:   /* C = RA*R^T */
106:   if (!C->ops->mattransposemultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
107:   (*C->ops->mattransposemultnumeric)(RA,R,C);
108:   return(0);
109: }

111: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
112: {
114:   Mat_Product    *product = C->product;
115:   Mat            A=product->A,R=product->B,RA;
116:   PetscReal      fill=product->fill;

119:   PetscInfo2((PetscObject)C,"for A %s, R %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
120:   /* RA = R*A */
121:   MatProductCreate(R,A,NULL,&RA);
122:   MatProductSetType(RA,MATPRODUCT_AB);
123:   MatProductSetAlgorithm(RA,MATPRODUCTALGORITHM_DEFAULT);
124:   MatProductSetFill(RA,fill);
125:   MatProductSetFromOptions(RA);
126:   MatProductSymbolic(RA);

128:   /* C = RA*R^T */
129:   MatProductSetType(C,MATPRODUCT_ABt);
130:   MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
131:   product->A = RA;
132:   MatProductSetFromOptions(C);
133:   MatProductSymbolic(C);

135:   /* resume user's original input matrix setting for A */
136:   product->A     = A;
137:   product->Dwork = RA; /* save here so it will be destroyed with product C */
138:   C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
139:   return(0);
140: }

142: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
143: {
145:   Mat_Product    *product = mat->product;
146:   Mat            A=product->A,BC=product->Dwork;

149:   /* Numeric BC = B*C */
150:   MatProductNumeric(BC);
151:   /* Numeric mat = A*BC */
152:   if (!mat->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
153:   (*mat->ops->matmultnumeric)(A,BC,mat);
154:   return(0);
155: }

157: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
158: {
160:   Mat_Product    *product = mat->product;
161:   Mat            B=product->B,C=product->C,BC;
162:   PetscReal      fill=product->fill;

165:   PetscInfo3((PetscObject)mat,"for A %s, B %s, C %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name,((PetscObject)product->C)->type_name);
166:   /* Symbolic BC = B*C */
167:   MatProductCreate(B,C,NULL,&BC);
168:   MatProductSetType(BC,MATPRODUCT_AB);
169:   MatProductSetAlgorithm(BC,MATPRODUCTALGORITHM_DEFAULT);
170:   MatProductSetFill(BC,fill);
171:   MatProductSetFromOptions(BC);
172:   MatProductSymbolic(BC);

174:   /* Symbolic mat = A*BC */
175:   MatProductSetType(mat,MATPRODUCT_AB);
176:   MatProductSetAlgorithm(mat,MATPRODUCTALGORITHM_DEFAULT);
177:   product->B     = BC;
178:   product->Dwork = BC;
179:   MatProductSetFromOptions(mat);
180:   MatProductSymbolic(mat);

182:   /* resume user's original input matrix setting for B */
183:   product->B = B;
184:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
185:   return(0);
186: }

188: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
189: {
191:   Mat_Product    *product = mat->product;

194:   switch (product->type) {
195:   case MATPRODUCT_PtAP:
196:     MatProductSymbolic_PtAP_Unsafe(mat);
197:     break;
198:   case MATPRODUCT_RARt:
199:     MatProductSymbolic_RARt_Unsafe(mat);
200:     break;
201:   case MATPRODUCT_ABC:
202:     MatProductSymbolic_ABC_Unsafe(mat);
203:     break;
204:   default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]);
205:   }
206:   return(0);
207: }

209: /* ----------------------------------------------- */
210: /*@C
211:    MatProductReplaceMats - Replace input matrices for a matrix product.

213:    Collective on Mat

215:    Input Parameters:
216: +  A - the matrix or NULL if not being replaced
217: .  B - the matrix or NULL if not being replaced
218: .  C - the matrix or NULL if not being replaced
219: -  D - the matrix product

221:    Level: intermediate

223:    Notes:
224:      To reuse the symbolic phase, input matrices must have exactly the same data structure as the replaced one.
225:      If the type of any of the input matrices is different than what was previously used, or their symmetry changed but
226:      the symbolic phase took advantage of their symmetry, the product is cleared and MatProductSetFromOptions()/MatProductSymbolic() are invoked again.

228: .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
229: @*/
230: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
231: {
233:   Mat_Product    *product;
234:   PetscBool      flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;

238:   MatCheckProduct(D,4);
239:   product = D->product;
240:   if (A) {
242:     PetscObjectReference((PetscObject)A);
243:     PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);
244:     if (product->symbolic_used_the_fact_A_is_symmetric && !A->symmetric) { /* symbolic was built around a symmetric A, but the new A is not anymore */
245:       flgA = PETSC_FALSE;
246:       product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
247:     }
248:     MatDestroy(&product->A);
249:     product->A = A;
250:   }
251:   if (B) {
253:     PetscObjectReference((PetscObject)B);
254:     PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);
255:     if (product->symbolic_used_the_fact_B_is_symmetric && !B->symmetric) {
256:       flgB = PETSC_FALSE;
257:       product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
258:     }
259:     MatDestroy(&product->B);
260:     product->B = B;
261:   }
262:   if (C) {
264:     PetscObjectReference((PetscObject)C);
265:     PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);
266:     if (product->symbolic_used_the_fact_C_is_symmetric && !C->symmetric) {
267:       flgC = PETSC_FALSE;
268:       product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
269:     }
270:     MatDestroy(&product->C);
271:     product->C = C;
272:   }
273:   /* Any of the replaced mats is of a different type, reset */
274:   if (!flgA || !flgB || !flgC) {
275:     if (D->product->destroy) {
276:       (*D->product->destroy)(D->product->data);
277:     }
278:     D->product->destroy = NULL;
279:     D->product->data = NULL;
280:     if (D->ops->productnumeric || D->ops->productsymbolic) {
281:       MatProductSetFromOptions(D);
282:       MatProductSymbolic(D);
283:     }
284:   }
285:   return(0);
286: }

288: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
289: {
291:   Mat_Product    *product = C->product;
292:   Mat            A = product->A, B = product->B;
293:   PetscInt       k, K = B->cmap->N;
294:   PetscBool      t = PETSC_TRUE,iscuda = PETSC_FALSE;
295:   PetscBool      Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
296:   char           *Btype = NULL,*Ctype = NULL;

299:   switch (product->type) {
300:   case MATPRODUCT_AB:
301:     t = PETSC_FALSE;
302:   case MATPRODUCT_AtB:
303:     break;
304:   default: SETERRQ3(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProductNumeric type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
305:   }
306:   if (PetscDefined(HAVE_CUDA)) {
307:     VecType vtype;

309:     MatGetVecType(A,&vtype);
310:     PetscStrcmp(vtype,VECCUDA,&iscuda);
311:     if (!iscuda) {
312:       PetscStrcmp(vtype,VECSEQCUDA,&iscuda);
313:     }
314:     if (!iscuda) {
315:       PetscStrcmp(vtype,VECMPICUDA,&iscuda);
316:     }
317:     if (iscuda) { /* Make sure we have up-to-date data on the GPU */
318:       PetscStrallocpy(((PetscObject)B)->type_name,&Btype);
319:       PetscStrallocpy(((PetscObject)C)->type_name,&Ctype);
320:       MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B);
321:       if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
322:         MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
323:         MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
324:       }
325:       MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C);
326:     } else { /* Make sure we have up-to-date data on the CPU */
327: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
328:       Bcpu = B->boundtocpu;
329:       Ccpu = C->boundtocpu;
330: #endif
331:       MatBindToCPU(B,PETSC_TRUE);
332:       MatBindToCPU(C,PETSC_TRUE);
333:     }
334:   }
335:   for (k=0;k<K;k++) {
336:     Vec x,y;

338:     MatDenseGetColumnVecRead(B,k,&x);
339:     MatDenseGetColumnVecWrite(C,k,&y);
340:     if (t) {
341:       MatMultTranspose(A,x,y);
342:     } else {
343:       MatMult(A,x,y);
344:     }
345:     MatDenseRestoreColumnVecRead(B,k,&x);
346:     MatDenseRestoreColumnVecWrite(C,k,&y);
347:   }
348:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
349:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
350:   if (PetscDefined(HAVE_CUDA)) {
351:     if (iscuda) {
352:       MatConvert(B,Btype,MAT_INPLACE_MATRIX,&B);
353:       MatConvert(C,Ctype,MAT_INPLACE_MATRIX,&C);
354:     } else {
355:       MatBindToCPU(B,Bcpu);
356:       MatBindToCPU(C,Ccpu);
357:     }
358:   }
359:   PetscFree(Btype);
360:   PetscFree(Ctype);
361:   return(0);
362: }

364: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
365: {
367:   Mat_Product    *product = C->product;
368:   Mat            A = product->A, B = product->B;
369:   PetscBool      isdense;

372:   switch (product->type) {
373:   case MATPRODUCT_AB:
374:     MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
375:     break;
376:   case MATPRODUCT_AtB:
377:     MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
378:     break;
379:   default: SETERRQ3(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
380:   }
381:   PetscObjectBaseTypeCompareAny((PetscObject)C,&isdense,MATSEQDENSE,MATMPIDENSE,"");
382:   if (!isdense) {
383:     MatSetType(C,((PetscObject)B)->type_name);
384:     /* If matrix type of C was not set or not dense, we need to reset the pointer */
385:     C->ops->productsymbolic = MatProductSymbolic_X_Dense;
386:   }
387:   C->ops->productnumeric = MatProductNumeric_X_Dense;
388:   MatSetUp(C);
389:   return(0);
390: }

392: /* a single driver to query the dispatching */
393: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
394: {
395:   PetscErrorCode    ierr;
396:   Mat_Product       *product = mat->product;
397:   PetscInt          Am,An,Bm,Bn,Cm,Cn;
398:   Mat               A = product->A,B = product->B,C = product->C;
399:   const char* const Bnames[] = { "B", "R", "P" };
400:   const char*       bname;
401:   PetscErrorCode    (*fA)(Mat);
402:   PetscErrorCode    (*fB)(Mat);
403:   PetscErrorCode    (*fC)(Mat);
404:   PetscErrorCode    (*f)(Mat)=NULL;

407:   mat->ops->productsymbolic = NULL;
408:   mat->ops->productnumeric = NULL;
409:   if (product->type == MATPRODUCT_UNSPECIFIED) return(0);
410:   if (!A) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing A mat");
411:   if (!B) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing B mat");
412:   if (product->type == MATPRODUCT_ABC && !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing C mat");
413:   if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
414:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
415:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
416:   else bname = Bnames[0];

418:   /* Check matrices sizes */
419:   Am = A->rmap->N;
420:   An = A->cmap->N;
421:   Bm = B->rmap->N;
422:   Bn = B->cmap->N;
423:   Cm = C ? C->rmap->N : 0;
424:   Cn = C ? C->cmap->N : 0;
425:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
426:   if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
427:   if (An != Bm) SETERRQ7(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of A and %s are incompatible for MatProductType %s: A %Dx%D, %s %Dx%D",bname,MatProductTypes[product->type],A->rmap->N,A->cmap->N,bname,B->rmap->N,B->cmap->N);
428:   if (Cm && Cm != Bn) SETERRQ5(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of B and C are incompatible for MatProductType %s: B %Dx%D, C %Dx%D",MatProductTypes[product->type],B->rmap->N,B->cmap->N,Cm,Cn);

430:   fA = A->ops->productsetfromoptions;
431:   fB = B->ops->productsetfromoptions;
432:   fC = C ? C->ops->productsetfromoptions : fA;
433:   if (C) {
434:     PetscInfo5(mat,"MatProductType %s for A %s, %s %s, C %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name,((PetscObject)C)->type_name);
435:   } else {
436:     PetscInfo4(mat,"MatProductType %s for A %s, %s %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name);
437:   }
438:   if (fA == fB && fA == fC && fA) {
439:     PetscInfo(mat,"  matching op\n");
440:     (*fA)(mat);
441:   } else {
442:     /* query MatProductSetFromOptions_Atype_Btype_Ctype */
443:     char  mtypes[256];
444:     PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
445:     PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
446:     PetscStrlcat(mtypes,"_",sizeof(mtypes));
447:     PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
448:     if (C) {
449:       PetscStrlcat(mtypes,"_",sizeof(mtypes));
450:       PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));
451:     }
452:     PetscStrlcat(mtypes,"_C",sizeof(mtypes));

454:     PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
455:     PetscInfo2(mat,"  querying %s from A? %p\n",mtypes,f);
456:     if (!f) {
457:       PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
458:       PetscInfo3(mat,"  querying %s from %s? %p\n",mtypes,bname,f);
459:     }
460:     if (!f && C) {
461:       PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
462:       PetscInfo2(mat,"  querying %s from C? %p\n",mtypes,f);
463:     }
464:     if (f) { (*f)(mat); }

466:     /* We may have found f but it did not succeed */
467:     /* some matrices (i.e. MATTRANSPOSE, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
468:     if (!mat->ops->productsymbolic) {
469:       PetscStrncpy(mtypes,"MatProductSetFromOptions_anytype_C",sizeof(mtypes));
470:       PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
471:       PetscInfo2(mat,"  querying %s from A? %p\n",mtypes,f);
472:       if (!f) {
473:         PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
474:         PetscInfo3(mat,"  querying %s from %s? %p\n",mtypes,bname,f);
475:       }
476:       if (!f && C) {
477:         PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
478:         PetscInfo2(mat,"  querying %s from C? %p\n",mtypes,f);
479:       }
480:     }
481:     if (f) { (*f)(mat); }
482:   }

484:   /* We may have found f but it did not succeed */
485:   if (!mat->ops->productsymbolic) {
486:     /* we can still compute the product if B is of type dense */
487:     if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
488:       PetscBool isdense;

490:       PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
491:       if (isdense) {

493:         mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
494:         PetscInfo(mat,"  using basic looping over columns of a dense matrix\n");
495:       }
496:     } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
497:       /*
498:          TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
499:                the compination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
500:                before computing the symbolic phase
501:       */
502:       PetscInfo(mat,"  symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n");
503:       mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
504:     }
505:   }
506:   if (!mat->ops->productsymbolic) {
507:     PetscInfo(mat,"  symbolic product is not supported\n");
508:   }
509:   return(0);
510: }

512: /*@C
513:    MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.

515:    Logically Collective on Mat

517:    Input Parameter:
518: .  mat - the matrix

520:    Options Database Keys:
521: .    -mat_product_clear - Clear intermediate data structures after MatProductNumeric() has been called

523:    Level: intermediate

525: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
526: @*/
527: PetscErrorCode MatProductSetFromOptions(Mat mat)
528: {

533:   MatCheckProduct(mat,1);
534:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot call MatProductSetFromOptions with already present data");
535:   PetscObjectOptionsBegin((PetscObject)mat);
536:   PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);
537:   PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");
538:   PetscOptionsEnd();
539:   MatProductSetFromOptions_Private(mat);
540:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after setup phase");
541:   return(0);
542: }

544: /*@C
545:    MatProductView - View a MatProduct

547:    Logically Collective on Mat

549:    Input Parameter:
550: .  mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()

552:    Level: intermediate

554: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
555: @*/
556: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
557: {

562:   if (!mat->product) return(0);
563:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);}
566:   if (mat->product->view) {
567:     (*mat->product->view)(mat,viewer);
568:   }
569:   return(0);
570: }

572: /* ----------------------------------------------- */
573: /* these are basic implementations relying on the old function pointers
574:  * they are dangerous and should be removed in the future */
575: PetscErrorCode MatProductNumeric_AB(Mat mat)
576: {
578:   Mat_Product    *product = mat->product;
579:   Mat            A=product->A,B=product->B;

582:   if (!mat->ops->matmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
583:   (*mat->ops->matmultnumeric)(A,B,mat);
584:   return(0);
585: }

587: PetscErrorCode MatProductNumeric_AtB(Mat mat)
588: {
590:   Mat_Product    *product = mat->product;
591:   Mat            A=product->A,B=product->B;

594:   if (!mat->ops->transposematmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
595:   (*mat->ops->transposematmultnumeric)(A,B,mat);
596:   return(0);
597: }

599: PetscErrorCode MatProductNumeric_ABt(Mat mat)
600: {
602:   Mat_Product    *product = mat->product;
603:   Mat            A=product->A,B=product->B;

606:   if (!mat->ops->mattransposemultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
607:   (*mat->ops->mattransposemultnumeric)(A,B,mat);
608:   return(0);
609: }

611: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
612: {
614:   Mat_Product    *product = mat->product;
615:   Mat            A=product->A,B=product->B;

618:   if (!mat->ops->ptapnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
619:   (*mat->ops->ptapnumeric)(A,B,mat);
620:   return(0);
621: }

623: PetscErrorCode MatProductNumeric_RARt(Mat mat)
624: {
626:   Mat_Product    *product = mat->product;
627:   Mat            A=product->A,B=product->B;

630:   if (!mat->ops->rartnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
631:   (*mat->ops->rartnumeric)(A,B,mat);
632:   return(0);
633: }

635: PetscErrorCode MatProductNumeric_ABC(Mat mat)
636: {
638:   Mat_Product    *product = mat->product;
639:   Mat            A=product->A,B=product->B,C=product->C;

642:   if (!mat->ops->matmatmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
643:   (*mat->ops->matmatmultnumeric)(A,B,C,mat);
644:   return(0);
645: }

647: /* ----------------------------------------------- */

649: /*@
650:    MatProductNumeric - Implement a matrix product with numerical values.

652:    Collective on Mat

654:    Input/Output Parameter:
655: .  mat - the matrix holding the product

657:    Level: intermediate

659:    Notes: MatProductSymbolic() must have been called on mat before calling this function

661: .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
662: @*/
663: PetscErrorCode MatProductNumeric(Mat mat)
664: {
666:   PetscLogEvent  eventtype=-1;

670:   MatCheckProduct(mat,1);
671:   /* log event */
672:   switch (mat->product->type) {
673:   case MATPRODUCT_AB:
674:     eventtype = MAT_MatMultNumeric;
675:     break;
676:   case MATPRODUCT_AtB:
677:     eventtype = MAT_TransposeMatMultNumeric;
678:     break;
679:   case MATPRODUCT_ABt:
680:     eventtype = MAT_MatTransposeMultNumeric;
681:     break;
682:   case MATPRODUCT_PtAP:
683:     eventtype = MAT_PtAPNumeric;
684:     break;
685:   case MATPRODUCT_RARt:
686:     eventtype = MAT_RARtNumeric;
687:     break;
688:   case MATPRODUCT_ABC:
689:     eventtype = MAT_MatMatMultNumeric;
690:     break;
691:   default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
692:   }
693:   if (mat->ops->productnumeric) {
694:     PetscLogEventBegin(eventtype,mat,0,0,0);
695:     (*mat->ops->productnumeric)(mat);
696:     PetscLogEventEnd(eventtype,mat,0,0,0);
697:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSymbolic() first");
698:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after numeric phase");
699:   if (mat->product->clear) {
700:     MatProductClear(mat);
701:   }
702:   PetscObjectStateIncrease((PetscObject)mat);
703:   return(0);
704: }

706: /* ----------------------------------------------- */
707: /* these are basic implementations relying on the old function pointers
708:  * they are dangerous and should be removed in the future */
709: PetscErrorCode MatProductSymbolic_AB(Mat mat)
710: {
712:   Mat_Product    *product = mat->product;
713:   Mat            A=product->A,B=product->B;

716:   if (!mat->ops->matmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
717:   (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);
718:   mat->ops->productnumeric = MatProductNumeric_AB;
719:   return(0);
720: }

722: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
723: {
725:   Mat_Product    *product = mat->product;
726:   Mat            A=product->A,B=product->B;

729:   if (!mat->ops->transposematmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
730:   (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
731:   mat->ops->productnumeric = MatProductNumeric_AtB;
732:   return(0);
733: }

735: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
736: {
738:   Mat_Product    *product = mat->product;
739:   Mat            A=product->A,B=product->B;

742:   if (!mat->ops->mattransposemultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
743:   (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
744:   mat->ops->productnumeric = MatProductNumeric_ABt;
745:   return(0);
746: }

748: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
749: {
751:   Mat_Product    *product = mat->product;
752:   Mat            A=product->A,B=product->B,C=product->C;

755:   if (!mat->ops->matmatmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
756:   (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
757:   mat->ops->productnumeric = MatProductNumeric_ABC;
758:   return(0);
759: }

761: /* ----------------------------------------------- */

763: /*@
764:    MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.

766:    Collective on Mat

768:    Input/Output Parameter:
769: .  mat - the matrix to hold a product

771:    Level: intermediate

773:    Notes: MatProductSetFromOptions() must have been called on mat before calling this function

775: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
776: @*/
777: PetscErrorCode MatProductSymbolic(Mat mat)
778: {
780:   PetscLogEvent  eventtype=-1;

784:   MatCheckProduct(mat,1);
785:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot run symbolic phase. Product data not empty");
786:   /* log event */
787:   switch (mat->product->type) {
788:   case MATPRODUCT_AB:
789:     eventtype = MAT_MatMultSymbolic;
790:     break;
791:   case MATPRODUCT_AtB:
792:     eventtype = MAT_TransposeMatMultSymbolic;
793:     break;
794:   case MATPRODUCT_ABt:
795:     eventtype = MAT_MatTransposeMultSymbolic;
796:     break;
797:   case MATPRODUCT_PtAP:
798:     eventtype = MAT_PtAPSymbolic;
799:     break;
800:   case MATPRODUCT_RARt:
801:     eventtype = MAT_RARtSymbolic;
802:     break;
803:   case MATPRODUCT_ABC:
804:     eventtype = MAT_MatMatMultSymbolic;
805:     break;
806:   default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
807:   }

809:   mat->ops->productnumeric = NULL;
810:   if (mat->ops->productsymbolic) {
811:     PetscLogEventBegin(eventtype,mat,0,0,0);
812:     (*mat->ops->productsymbolic)(mat);
813:     PetscLogEventEnd(eventtype,mat,0,0,0);
814:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSetFromOptions() first");
815:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after symbolic phase");
816:   if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Symbolic phase did not specify the numeric phase");
817:   return(0);
818: }

820: /*@
821:    MatProductSetFill - Set an expected fill of the matrix product.

823:    Collective on Mat

825:    Input Parameters:
826: +  mat - the matrix product
827: -  fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use PETSC_DEFAULT if you do not have a good estimate. If the product is a dense matrix, this is irrelevant.

829:    Level: intermediate

831: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
832: @*/
833: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
834: {
837:   MatCheckProduct(mat,1);
838:   if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
839:   else mat->product->fill = fill;
840:   return(0);
841: }

843: /*@
844:    MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.

846:    Collective on Mat

848:    Input Parameters:
849: +  mat - the matrix product
850: -  alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.

852:    Level: intermediate

854: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
855: @*/
856: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
857: {

862:   MatCheckProduct(mat,1);
863:   PetscFree(mat->product->alg);
864:   PetscStrallocpy(alg,&mat->product->alg);
865:   return(0);
866: }

868: /*@
869:    MatProductSetType - Sets a particular matrix product type

871:    Collective on Mat

873:    Input Parameters:
874: +  mat - the matrix
875: -  productype   - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.

877:    Level: intermediate

879: .seealso: MatProductCreate(), MatProductType
880: @*/
881: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
882: {

887:   MatCheckProduct(mat,1);
889:   if (productype != mat->product->type) {
890:     if (mat->product->destroy) {
891:       (*mat->product->destroy)(mat->product->data);
892:     }
893:     mat->product->destroy = NULL;
894:     mat->product->data = NULL;
895:     mat->ops->productsymbolic = NULL;
896:     mat->ops->productnumeric  = NULL;
897:   }
898:   mat->product->type = productype;
899:   return(0);
900: }

902: /*@
903:    MatProductClear - Clears matrix product internal structure.

905:    Collective on Mat

907:    Input Parameters:
908: .  mat - the product matrix

910:    Level: intermediate

912:    Notes: this function should be called to remove any intermediate data used by the product
913:           After having called this function, MatProduct operations can no longer be used on mat
914: @*/
915: PetscErrorCode MatProductClear(Mat mat)
916: {
918:   Mat_Product    *product = mat->product;

922:   if (product) {
923:     MatDestroy(&product->A);
924:     MatDestroy(&product->B);
925:     MatDestroy(&product->C);
926:     PetscFree(product->alg);
927:     MatDestroy(&product->Dwork);
928:     if (product->destroy) {
929:       (*product->destroy)(product->data);
930:     }
931:   }
932:   PetscFree(mat->product);
933:   mat->ops->productsymbolic = NULL;
934:   mat->ops->productnumeric = NULL;
935:   return(0);
936: }

938: /* Create a supporting struct and attach it to the matrix product */
939: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
940: {
942:   Mat_Product    *product=NULL;

946:   if (D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product already present");
947:   PetscNewLog(D,&product);
948:   product->A        = A;
949:   product->B        = B;
950:   product->C        = C;
951:   product->type     = MATPRODUCT_UNSPECIFIED;
952:   product->Dwork    = NULL;
953:   product->api_user = PETSC_FALSE;
954:   product->clear    = PETSC_FALSE;
955:   D->product        = product;

957:   MatProductSetAlgorithm(D,MATPRODUCTALGORITHM_DEFAULT);
958:   MatProductSetFill(D,PETSC_DEFAULT);

960:   PetscObjectReference((PetscObject)A);
961:   PetscObjectReference((PetscObject)B);
962:   PetscObjectReference((PetscObject)C);
963:   return(0);
964: }

966: /*@
967:    MatProductCreateWithMat - Setup a given matrix as a matrix product.

969:    Collective on Mat

971:    Input Parameters:
972: +  A - the first matrix
973: .  B - the second matrix
974: .  C - the third matrix (optional)
975: -  D - the matrix which will be used as a product

977:    Output Parameters:
978: .  D - the product matrix

980:    Notes:
981:      Any product data attached to D will be cleared

983:    Level: intermediate

985: .seealso: MatProductCreate(), MatProductClear()
986: @*/
987: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
988: {

994:   MatCheckPreallocated(A,1);
995:   if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
996:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

1000:   MatCheckPreallocated(B,2);
1001:   if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1002:   if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

1004:   if (C) {
1007:     MatCheckPreallocated(C,3);
1008:     if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1009:     if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1010:   }

1014:   MatCheckPreallocated(D,4);
1015:   if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1016:   if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

1018:   /* Create a supporting struct and attach it to D */
1019:   MatProductClear(D);
1020:   MatProductCreate_Private(A,B,C,D);
1021:   return(0);
1022: }

1024: /*@
1025:    MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.

1027:    Collective on Mat

1029:    Input Parameters:
1030: +  A - the first matrix
1031: .  B - the second matrix
1032: -  C - the third matrix (optional)

1034:    Output Parameters:
1035: .  D - the product matrix

1037:    Level: intermediate

1039: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1040: @*/
1041: PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1042: {

1050:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix A");
1051:   if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix B");

1053:   if (C) {
1056:     if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix C");
1057:   }

1060:   MatCreate(PetscObjectComm((PetscObject)A),D);
1061:   MatProductCreate_Private(A,B,C,*D);
1062:   return(0);
1063: }

1065: /*
1066:    These are safe basic implementations of ABC, RARt and PtAP
1067:    that do not rely on mat->ops->matmatop function pointers.
1068:    They only use the MatProduct API and are currently used by
1069:    cuSPARSE and KOKKOS-KERNELS backends
1070: */
1071: typedef struct {
1072:   Mat BC;
1073:   Mat ABC;
1074: } MatMatMatPrivate;

1076: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1077: {
1078:   PetscErrorCode   ierr;
1079:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

1082:   MatDestroy(&mmdata->BC);
1083:   MatDestroy(&mmdata->ABC);
1084:   PetscFree(data);
1085:   return(0);
1086: }

1088: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1089: {
1090:   PetscErrorCode   ierr;
1091:   Mat_Product      *product = mat->product;
1092:   MatMatMatPrivate *mmabc;

1095:   MatCheckProduct(mat,1);
1096:   if (!mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data empty");
1097:   mmabc = (MatMatMatPrivate *)mat->product->data;
1098:   if (!mmabc->BC->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
1099:   /* use function pointer directly to prevent logging */
1100:   (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1101:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1102:   mat->product = mmabc->ABC->product;
1103:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1104:   if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
1105:   /* use function pointer directly to prevent logging */
1106:   (*mat->ops->productnumeric)(mat);
1107:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1108:   mat->product = product;
1109:   return(0);
1110: }

1112: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1113: {
1114:   PetscErrorCode   ierr;
1115:   Mat_Product      *product = mat->product;
1116:   Mat              A, B ,C;
1117:   MatProductType   p1,p2;
1118:   MatMatMatPrivate *mmabc;
1119:   const char       *prefix;

1122:   MatCheckProduct(mat,1);
1123:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data not empty");
1124:   MatGetOptionsPrefix(mat,&prefix);
1125:   PetscNew(&mmabc);
1126:   product->data    = mmabc;
1127:   product->destroy = MatDestroy_MatMatMatPrivate;
1128:   switch (product->type) {
1129:   case MATPRODUCT_PtAP:
1130:     p1 = MATPRODUCT_AB;
1131:     p2 = MATPRODUCT_AtB;
1132:     A = product->B;
1133:     B = product->A;
1134:     C = product->B;
1135:     break;
1136:   case MATPRODUCT_RARt:
1137:     p1 = MATPRODUCT_ABt;
1138:     p2 = MATPRODUCT_AB;
1139:     A = product->B;
1140:     B = product->A;
1141:     C = product->B;
1142:     break;
1143:   case MATPRODUCT_ABC:
1144:     p1 = MATPRODUCT_AB;
1145:     p2 = MATPRODUCT_AB;
1146:     A = product->A;
1147:     B = product->B;
1148:     C = product->C;
1149:     break;
1150:   default:
1151:     SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Not for ProductType %s",MatProductTypes[product->type]);
1152:   }
1153:   MatProductCreate(B,C,NULL,&mmabc->BC);
1154:   MatSetOptionsPrefix(mmabc->BC,prefix);
1155:   MatAppendOptionsPrefix(mmabc->BC,"P1_");
1156:   MatProductSetType(mmabc->BC,p1);
1157:   MatProductSetAlgorithm(mmabc->BC,MATPRODUCTALGORITHM_DEFAULT);
1158:   MatProductSetFill(mmabc->BC,product->fill);
1159:   mmabc->BC->product->api_user = product->api_user;
1160:   MatProductSetFromOptions(mmabc->BC);
1161:   if (!mmabc->BC->ops->productsymbolic) SETERRQ3(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Symbolic ProductType %s not supported with %s and %s",MatProductTypes[p1],((PetscObject)B)->type_name,((PetscObject)C)->type_name);
1162:   /* use function pointer directly to prevent logging */
1163:   (*mmabc->BC->ops->productsymbolic)(mmabc->BC);

1165:   MatProductCreate(A,mmabc->BC,NULL,&mmabc->ABC);
1166:   MatSetOptionsPrefix(mmabc->ABC,prefix);
1167:   MatAppendOptionsPrefix(mmabc->ABC,"P2_");
1168:   MatProductSetType(mmabc->ABC,p2);
1169:   MatProductSetAlgorithm(mmabc->ABC,MATPRODUCTALGORITHM_DEFAULT);
1170:   MatProductSetFill(mmabc->ABC,product->fill);
1171:   mmabc->ABC->product->api_user = product->api_user;
1172:   MatProductSetFromOptions(mmabc->ABC);
1173:   if (!mmabc->ABC->ops->productsymbolic) SETERRQ3(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Symbolic ProductType %s not supported with %s and %s",MatProductTypes[p2],((PetscObject)A)->type_name,((PetscObject)mmabc->BC)->type_name);
1174:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1175:   mat->product = mmabc->ABC->product;
1176:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1177:   /* use function pointer directly to prevent logging */
1178:   (*mat->ops->productsymbolic)(mat);
1179:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1180:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
1181:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
1182:   mat->product                    = product;
1183:   return(0);
1184: }