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 MATTRANSPOSEVIRTUAL)
 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: {
 49:   Mat_Product *product = C->product;
 50:   Mat          P = product->B, AP = product->Dwork;

 52:   /* AP = A*P */
 53:   MatProductNumeric(AP);
 54:   /* C = P^T*AP */
 55:   (*C->ops->transposematmultnumeric)(P, AP, C);
 56:   return 0;
 57: }

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

 65:   PetscInfo((PetscObject)C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name);
 66:   /* AP = A*P */
 67:   MatProductCreate(A, P, NULL, &AP);
 68:   MatProductSetType(AP, MATPRODUCT_AB);
 69:   MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT);
 70:   MatProductSetFill(AP, fill);
 71:   MatProductSetFromOptions(AP);
 72:   MatProductSymbolic(AP);

 74:   /* C = P^T*AP */
 75:   MatProductSetType(C, MATPRODUCT_AtB);
 76:   MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT);
 77:   product->A = P;
 78:   product->B = AP;
 79:   MatProductSetFromOptions(C);
 80:   MatProductSymbolic(C);

 82:   /* resume user's original input matrix setting for A and B */
 83:   product->A     = A;
 84:   product->B     = P;
 85:   product->Dwork = AP;

 87:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
 88:   return 0;
 89: }

 91: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
 92: {
 93:   Mat_Product *product = C->product;
 94:   Mat          R = product->B, RA = product->Dwork;

 96:   /* RA = R*A */
 97:   MatProductNumeric(RA);
 98:   /* C = RA*R^T */
 99:   (*C->ops->mattransposemultnumeric)(RA, R, C);
100:   return 0;
101: }

103: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
104: {
105:   Mat_Product *product = C->product;
106:   Mat          A = product->A, R = product->B, RA;
107:   PetscReal    fill = product->fill;

109:   PetscInfo((PetscObject)C, "for A %s, R %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name);
110:   /* RA = R*A */
111:   MatProductCreate(R, A, NULL, &RA);
112:   MatProductSetType(RA, MATPRODUCT_AB);
113:   MatProductSetAlgorithm(RA, MATPRODUCTALGORITHMDEFAULT);
114:   MatProductSetFill(RA, fill);
115:   MatProductSetFromOptions(RA);
116:   MatProductSymbolic(RA);

118:   /* C = RA*R^T */
119:   MatProductSetType(C, MATPRODUCT_ABt);
120:   MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT);
121:   product->A = RA;
122:   MatProductSetFromOptions(C);
123:   MatProductSymbolic(C);

125:   /* resume user's original input matrix setting for A */
126:   product->A             = A;
127:   product->Dwork         = RA; /* save here so it will be destroyed with product C */
128:   C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
129:   return 0;
130: }

132: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
133: {
134:   Mat_Product *product = mat->product;
135:   Mat          A = product->A, BC = product->Dwork;

137:   /* Numeric BC = B*C */
138:   MatProductNumeric(BC);
139:   /* Numeric mat = A*BC */
140:   (*mat->ops->matmultnumeric)(A, BC, mat);
141:   return 0;
142: }

144: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
145: {
146:   Mat_Product *product = mat->product;
147:   Mat          B = product->B, C = product->C, BC;
148:   PetscReal    fill = product->fill;

150:   PetscInfo((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);
151:   /* Symbolic BC = B*C */
152:   MatProductCreate(B, C, NULL, &BC);
153:   MatProductSetType(BC, MATPRODUCT_AB);
154:   MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT);
155:   MatProductSetFill(BC, fill);
156:   MatProductSetFromOptions(BC);
157:   MatProductSymbolic(BC);

159:   /* Symbolic mat = A*BC */
160:   MatProductSetType(mat, MATPRODUCT_AB);
161:   MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT);
162:   product->B     = BC;
163:   product->Dwork = BC;
164:   MatProductSetFromOptions(mat);
165:   MatProductSymbolic(mat);

167:   /* resume user's original input matrix setting for B */
168:   product->B               = B;
169:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
170:   return 0;
171: }

173: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
174: {
175:   Mat_Product *product = mat->product;

177:   switch (product->type) {
178:   case MATPRODUCT_PtAP:
179:     MatProductSymbolic_PtAP_Unsafe(mat);
180:     break;
181:   case MATPRODUCT_RARt:
182:     MatProductSymbolic_RARt_Unsafe(mat);
183:     break;
184:   case MATPRODUCT_ABC:
185:     MatProductSymbolic_ABC_Unsafe(mat);
186:     break;
187:   default:
188:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
189:   }
190:   return 0;
191: }

193: /* ----------------------------------------------- */
194: /*@C
195:    MatProductReplaceMats - Replace input matrices for a matrix product.

197:    Collective on Mat

199:    Input Parameters:
200: +  A - the matrix or NULL if not being replaced
201: .  B - the matrix or NULL if not being replaced
202: .  C - the matrix or NULL if not being replaced
203: -  D - the matrix product

205:    Level: intermediate

207:    Note:
208:      To reuse the symbolic phase, the input matrices must have exactly the same data structure as the replaced one.
209:      If the type of any of the input matrices is different than what was previously used, or their symmetry flag changed but
210:      the symbolic phase took advantage of their symmetry, the product is cleared and `MatProductSetFromOptions()` and `MatProductSymbolic()` are invoked again.

212: .seealso: `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic().` `MatProductClear()`
213: @*/
214: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
215: {
216:   Mat_Product *product;
217:   PetscBool    flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;

220:   MatCheckProduct(D, 4);
221:   product = D->product;
222:   if (A) {
224:     PetscObjectReference((PetscObject)A);
225:     PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA);
226:     MatIsSymmetricKnown(A, &isset, &issym);
227:     if (product->symbolic_used_the_fact_A_is_symmetric && isset && !issym) { /* symbolic was built around a symmetric A, but the new A is not anymore */
228:       flgA                                           = PETSC_FALSE;
229:       product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
230:     }
231:     MatDestroy(&product->A);
232:     product->A = A;
233:   }
234:   if (B) {
236:     PetscObjectReference((PetscObject)B);
237:     PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB);
238:     MatIsSymmetricKnown(B, &isset, &issym);
239:     if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
240:       flgB                                           = PETSC_FALSE;
241:       product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
242:     }
243:     MatDestroy(&product->B);
244:     product->B = B;
245:   }
246:   if (C) {
248:     PetscObjectReference((PetscObject)C);
249:     PetscObjectTypeCompare((PetscObject)product->C, ((PetscObject)C)->type_name, &flgC);
250:     MatIsSymmetricKnown(C, &isset, &issym);
251:     if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
252:       flgC                                           = PETSC_FALSE;
253:       product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
254:     }
255:     MatDestroy(&product->C);
256:     product->C = C;
257:   }
258:   /* Any of the replaced mats is of a different type, reset */
259:   if (!flgA || !flgB || !flgC) {
260:     if (D->product->destroy) (*D->product->destroy)(D->product->data);
261:     D->product->destroy = NULL;
262:     D->product->data    = NULL;
263:     if (D->ops->productnumeric || D->ops->productsymbolic) {
264:       MatProductSetFromOptions(D);
265:       MatProductSymbolic(D);
266:     }
267:   }
268:   return 0;
269: }

271: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
272: {
273:   Mat_Product *product = C->product;
274:   Mat          A = product->A, B = product->B;
275:   PetscInt     k, K              = B->cmap->N;
276:   PetscBool    t = PETSC_TRUE, iscuda = PETSC_FALSE;
277:   PetscBool    Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
278:   char        *Btype = NULL, *Ctype = NULL;

280:   switch (product->type) {
281:   case MATPRODUCT_AB:
282:     t = PETSC_FALSE;
283:   case MATPRODUCT_AtB:
284:     break;
285:   default:
286:     SETERRQ(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);
287:   }
288:   if (PetscDefined(HAVE_CUDA)) {
289:     VecType vtype;

291:     MatGetVecType(A, &vtype);
292:     PetscStrcmp(vtype, VECCUDA, &iscuda);
293:     if (!iscuda) PetscStrcmp(vtype, VECSEQCUDA, &iscuda);
294:     if (!iscuda) PetscStrcmp(vtype, VECMPICUDA, &iscuda);
295:     if (iscuda) { /* Make sure we have up-to-date data on the GPU */
296:       PetscStrallocpy(((PetscObject)B)->type_name, &Btype);
297:       PetscStrallocpy(((PetscObject)C)->type_name, &Ctype);
298:       MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B);
299:       if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
300:         MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
301:         MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
302:       }
303:       MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C);
304:     } else { /* Make sure we have up-to-date data on the CPU */
305: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
306:       Bcpu = B->boundtocpu;
307:       Ccpu = C->boundtocpu;
308: #endif
309:       MatBindToCPU(B, PETSC_TRUE);
310:       MatBindToCPU(C, PETSC_TRUE);
311:     }
312:   }
313:   for (k = 0; k < K; k++) {
314:     Vec x, y;

316:     MatDenseGetColumnVecRead(B, k, &x);
317:     MatDenseGetColumnVecWrite(C, k, &y);
318:     if (t) {
319:       MatMultTranspose(A, x, y);
320:     } else {
321:       MatMult(A, x, y);
322:     }
323:     MatDenseRestoreColumnVecRead(B, k, &x);
324:     MatDenseRestoreColumnVecWrite(C, k, &y);
325:   }
326:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
327:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
328:   if (PetscDefined(HAVE_CUDA)) {
329:     if (iscuda) {
330:       MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B);
331:       MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C);
332:     } else {
333:       MatBindToCPU(B, Bcpu);
334:       MatBindToCPU(C, Ccpu);
335:     }
336:   }
337:   PetscFree(Btype);
338:   PetscFree(Ctype);
339:   return 0;
340: }

342: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
343: {
344:   Mat_Product *product = C->product;
345:   Mat          A = product->A, B = product->B;
346:   PetscBool    isdense;

348:   switch (product->type) {
349:   case MATPRODUCT_AB:
350:     MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N);
351:     break;
352:   case MATPRODUCT_AtB:
353:     MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N);
354:     break;
355:   default:
356:     SETERRQ(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);
357:   }
358:   PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, "");
359:   if (!isdense) {
360:     MatSetType(C, ((PetscObject)B)->type_name);
361:     /* If matrix type of C was not set or not dense, we need to reset the pointer */
362:     C->ops->productsymbolic = MatProductSymbolic_X_Dense;
363:   }
364:   C->ops->productnumeric = MatProductNumeric_X_Dense;
365:   MatSetUp(C);
366:   return 0;
367: }

369: /* a single driver to query the dispatching */
370: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
371: {
372:   Mat_Product      *product = mat->product;
373:   PetscInt          Am, An, Bm, Bn, Cm, Cn;
374:   Mat               A = product->A, B = product->B, C = product->C;
375:   const char *const Bnames[] = {"B", "R", "P"};
376:   const char       *bname;
377:   PetscErrorCode (*fA)(Mat);
378:   PetscErrorCode (*fB)(Mat);
379:   PetscErrorCode (*fC)(Mat);
380:   PetscErrorCode (*f)(Mat) = NULL;

382:   mat->ops->productsymbolic = NULL;
383:   mat->ops->productnumeric  = NULL;
384:   if (product->type == MATPRODUCT_UNSPECIFIED) return 0;
388:   if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
389:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
390:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
391:   else bname = Bnames[0];

393:   /* Check matrices sizes */
394:   Am = A->rmap->N;
395:   An = A->cmap->N;
396:   Bm = B->rmap->N;
397:   Bn = B->cmap->N;
398:   Cm = C ? C->rmap->N : 0;
399:   Cn = C ? C->cmap->N : 0;
400:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
401:     PetscInt t = Bn;
402:     Bn         = Bm;
403:     Bm         = t;
404:   }
405:   if (product->type == MATPRODUCT_AtB) {
406:     PetscInt t = An;
407:     An         = Am;
408:     Am         = t;
409:   }
411:              MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
413:              MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);

415:   fA = A->ops->productsetfromoptions;
416:   fB = B->ops->productsetfromoptions;
417:   fC = C ? C->ops->productsetfromoptions : fA;
418:   if (C) {
419:     PetscInfo(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);
420:   } else {
421:     PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name);
422:   }
423:   if (fA == fB && fA == fC && fA) {
424:     PetscInfo(mat, "  matching op\n");
425:     (*fA)(mat);
426:   }
427:   /* We may have found f but it did not succeed */
428:   if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
429:     char mtypes[256];
430:     PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes));
431:     PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes));
432:     PetscStrlcat(mtypes, "_", sizeof(mtypes));
433:     PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes));
434:     if (C) {
435:       PetscStrlcat(mtypes, "_", sizeof(mtypes));
436:       PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes));
437:     }
438:     PetscStrlcat(mtypes, "_C", sizeof(mtypes));

440:     PetscObjectQueryFunction((PetscObject)A, mtypes, &f);
441:     PetscInfo(mat, "  querying %s from A? %p\n", mtypes, f);
442:     if (!f) {
443:       PetscObjectQueryFunction((PetscObject)B, mtypes, &f);
444:       PetscInfo(mat, "  querying %s from %s? %p\n", mtypes, bname, f);
445:     }
446:     if (!f && C) {
447:       PetscObjectQueryFunction((PetscObject)C, mtypes, &f);
448:       PetscInfo(mat, "  querying %s from C? %p\n", mtypes, f);
449:     }
450:     if (f) (*f)(mat);

452:     /* We may have found f but it did not succeed */
453:     /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
454:     if (!mat->ops->productsymbolic) {
455:       PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes));
456:       PetscObjectQueryFunction((PetscObject)A, mtypes, &f);
457:       PetscInfo(mat, "  querying %s from A? %p\n", mtypes, f);
458:       if (!f) {
459:         PetscObjectQueryFunction((PetscObject)B, mtypes, &f);
460:         PetscInfo(mat, "  querying %s from %s? %p\n", mtypes, bname, f);
461:       }
462:       if (!f && C) {
463:         PetscObjectQueryFunction((PetscObject)C, mtypes, &f);
464:         PetscInfo(mat, "  querying %s from C? %p\n", mtypes, f);
465:       }
466:     }
467:     if (f) (*f)(mat);
468:   }

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

476:       PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, "");
477:       if (isdense) {
478:         mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
479:         PetscInfo(mat, "  using basic looping over columns of a dense matrix\n");
480:       }
481:     } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
482:       /*
483:          TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
484:                the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
485:                before computing the symbolic phase
486:       */
487:       PetscInfo(mat, "  symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n");
488:       mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
489:     }
490:   }
491:   if (!mat->ops->productsymbolic) PetscInfo(mat, "  symbolic product is not supported\n");
492:   return 0;
493: }

495: /*@C
496:    MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product where the type, the algorithm etc are determined from the options database.

498:    Logically Collective on Mat

500:    Input Parameter:
501: .  mat - the matrix

503:    Options Database Keys:
504: .    -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called

506:    Level: intermediate

508: .seealso: `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
509: @*/
510: PetscErrorCode MatProductSetFromOptions(Mat mat)
511: {
513:   MatCheckProduct(mat, 1);
515:   PetscObjectOptionsBegin((PetscObject)mat);
516:   PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL);
517:   PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()");
518:   PetscOptionsEnd();
519:   MatProductSetFromOptions_Private(mat);
521:   return 0;
522: }

524: /*@C
525:    MatProductView - View the private `Mat_Product` algorithm object within a matrix

527:    Logically Collective on mat

529:    Input Parameter:
530: .  mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`

532:    Level: intermediate

534: .seealso: `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
535: @*/
536: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
537: {
539:   if (!mat->product) return 0;
540:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer);
543:   if (mat->product->view) (*mat->product->view)(mat, viewer);
544:   return 0;
545: }

547: /* ----------------------------------------------- */
548: /* these are basic implementations relying on the old function pointers
549:  * they are dangerous and should be removed in the future */
550: PetscErrorCode MatProductNumeric_AB(Mat mat)
551: {
552:   Mat_Product *product = mat->product;
553:   Mat          A = product->A, B = product->B;

555:   (*mat->ops->matmultnumeric)(A, B, mat);
556:   return 0;
557: }

559: PetscErrorCode MatProductNumeric_AtB(Mat mat)
560: {
561:   Mat_Product *product = mat->product;
562:   Mat          A = product->A, B = product->B;

564:   (*mat->ops->transposematmultnumeric)(A, B, mat);
565:   return 0;
566: }

568: PetscErrorCode MatProductNumeric_ABt(Mat mat)
569: {
570:   Mat_Product *product = mat->product;
571:   Mat          A = product->A, B = product->B;

573:   (*mat->ops->mattransposemultnumeric)(A, B, mat);
574:   return 0;
575: }

577: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
578: {
579:   Mat_Product *product = mat->product;
580:   Mat          A = product->A, B = product->B;

582:   (*mat->ops->ptapnumeric)(A, B, mat);
583:   return 0;
584: }

586: PetscErrorCode MatProductNumeric_RARt(Mat mat)
587: {
588:   Mat_Product *product = mat->product;
589:   Mat          A = product->A, B = product->B;

591:   (*mat->ops->rartnumeric)(A, B, mat);
592:   return 0;
593: }

595: PetscErrorCode MatProductNumeric_ABC(Mat mat)
596: {
597:   Mat_Product *product = mat->product;
598:   Mat          A = product->A, B = product->B, C = product->C;

600:   (*mat->ops->matmatmultnumeric)(A, B, C, mat);
601:   return 0;
602: }

604: /* ----------------------------------------------- */

606: /*@
607:    MatProductNumeric - Compute a matrix product with numerical values.

609:    Collective on mat

611:    Input/Output Parameter:
612: .  mat - the matrix holding the product

614:    Level: intermediate

616:    Note:
617:    `MatProductSymbolic()` must have been called on mat before calling this function

619: .seealso: `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
620: @*/
621: PetscErrorCode MatProductNumeric(Mat mat)
622: {
623:   PetscLogEvent eventtype = -1;
624:   PetscBool     missing   = PETSC_FALSE;

627:   MatCheckProduct(mat, 1);
628:   switch (mat->product->type) {
629:   case MATPRODUCT_AB:
630:     eventtype = MAT_MatMultNumeric;
631:     break;
632:   case MATPRODUCT_AtB:
633:     eventtype = MAT_TransposeMatMultNumeric;
634:     break;
635:   case MATPRODUCT_ABt:
636:     eventtype = MAT_MatTransposeMultNumeric;
637:     break;
638:   case MATPRODUCT_PtAP:
639:     eventtype = MAT_PtAPNumeric;
640:     break;
641:   case MATPRODUCT_RARt:
642:     eventtype = MAT_RARtNumeric;
643:     break;
644:   case MATPRODUCT_ABC:
645:     eventtype = MAT_MatMatMultNumeric;
646:     break;
647:   default:
648:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
649:   }

651:   if (mat->ops->productnumeric) {
652:     PetscLogEventBegin(eventtype, mat, 0, 0, 0);
653:     PetscUseTypeMethod(mat, productnumeric);
654:     PetscLogEventEnd(eventtype, mat, 0, 0, 0);
655:   } else missing = PETSC_TRUE;

657:   if (missing || !mat->product) {
658:     char errstr[256];

660:     if (mat->product->type == MATPRODUCT_ABC) {
661:       PetscSNPrintf(errstr, 256, "%s with A %s, B %s, C %s", MatProductTypes[mat->product->type], ((PetscObject)mat->product->A)->type_name, ((PetscObject)mat->product->B)->type_name, ((PetscObject)mat->product->C)->type_name);
662:     } else {
663:       PetscSNPrintf(errstr, 256, "%s with A %s, B %s", MatProductTypes[mat->product->type], ((PetscObject)mat->product->A)->type_name, ((PetscObject)mat->product->B)->type_name);
664:     }
667:   }

669:   if (mat->product->clear) MatProductClear(mat);
670:   PetscObjectStateIncrease((PetscObject)mat);
671:   return 0;
672: }

674: /* ----------------------------------------------- */
675: /* these are basic implementations relying on the old function pointers
676:  * they are dangerous and should be removed in the future */
677: PetscErrorCode MatProductSymbolic_AB(Mat mat)
678: {
679:   Mat_Product *product = mat->product;
680:   Mat          A = product->A, B = product->B;

682:   (*mat->ops->matmultsymbolic)(A, B, product->fill, mat);
683:   mat->ops->productnumeric = MatProductNumeric_AB;
684:   return 0;
685: }

687: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
688: {
689:   Mat_Product *product = mat->product;
690:   Mat          A = product->A, B = product->B;

692:   (*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat);
693:   mat->ops->productnumeric = MatProductNumeric_AtB;
694:   return 0;
695: }

697: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
698: {
699:   Mat_Product *product = mat->product;
700:   Mat          A = product->A, B = product->B;

702:   (*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat);
703:   mat->ops->productnumeric = MatProductNumeric_ABt;
704:   return 0;
705: }

707: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
708: {
709:   Mat_Product *product = mat->product;
710:   Mat          A = product->A, B = product->B, C = product->C;

712:   (*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat);
713:   mat->ops->productnumeric = MatProductNumeric_ABC;
714:   return 0;
715: }

717: /* ----------------------------------------------- */

719: /*@
720:    MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical product done with
721:   `MatProductNumeric()`

723:    Collective on mat

725:    Input/Output Parameter:
726: .  mat - the matrix to hold a product

728:    Level: intermediate

730:    Note:
731:    `MatProductSetFromOptions()` must have been called on mat before calling this function

733: .seealso: `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
734: @*/
735: PetscErrorCode MatProductSymbolic(Mat mat)
736: {
737:   PetscLogEvent eventtype = -1;
738:   PetscBool     missing   = PETSC_FALSE;

741:   MatCheckProduct(mat, 1);
743:   switch (mat->product->type) {
744:   case MATPRODUCT_AB:
745:     eventtype = MAT_MatMultSymbolic;
746:     break;
747:   case MATPRODUCT_AtB:
748:     eventtype = MAT_TransposeMatMultSymbolic;
749:     break;
750:   case MATPRODUCT_ABt:
751:     eventtype = MAT_MatTransposeMultSymbolic;
752:     break;
753:   case MATPRODUCT_PtAP:
754:     eventtype = MAT_PtAPSymbolic;
755:     break;
756:   case MATPRODUCT_RARt:
757:     eventtype = MAT_RARtSymbolic;
758:     break;
759:   case MATPRODUCT_ABC:
760:     eventtype = MAT_MatMatMultSymbolic;
761:     break;
762:   default:
763:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
764:   }
765:   mat->ops->productnumeric = NULL;
766:   if (mat->ops->productsymbolic) {
767:     PetscLogEventBegin(eventtype, mat, 0, 0, 0);
768:     PetscUseTypeMethod(mat, productsymbolic);
769:     PetscLogEventEnd(eventtype, mat, 0, 0, 0);
770:   } else missing = PETSC_TRUE;

772:   if (missing || !mat->product || !mat->ops->productnumeric) {
773:     char errstr[256];

775:     if (mat->product->type == MATPRODUCT_ABC) {
776:       PetscSNPrintf(errstr, 256, "%s with A %s, B %s, C %s", MatProductTypes[mat->product->type], ((PetscObject)mat->product->A)->type_name, ((PetscObject)mat->product->B)->type_name, ((PetscObject)mat->product->C)->type_name);
777:     } else {
778:       PetscSNPrintf(errstr, 256, "%s with A %s, B %s", MatProductTypes[mat->product->type], ((PetscObject)mat->product->A)->type_name, ((PetscObject)mat->product->B)->type_name);
779:     }
782:   }
783:   return 0;
784: }

786: /*@
787:    MatProductSetFill - Set an expected fill of the matrix product.

789:    Collective on Mat

791:    Input Parameters:
792: +  mat - the matrix product result matrix
793: -  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 value is not used.

795:    Level: intermediate

797: .seealso: `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
798: @*/
799: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
800: {
802:   MatCheckProduct(mat, 1);
803:   if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
804:   else mat->product->fill = fill;
805:   return 0;
806: }

808: /*@
809:    MatProductSetAlgorithm - Requests a particular algorithm for a matrix product computation that will perform to compute the given matrix

811:    Collective on mat

813:    Input Parameters:
814: +  mat - the matrix product
815: -  alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.

817:    Options Database Key:
818: .  -mat_product_algorithm <algorithm> - Sets the algorithm; use -help for a list
819:     of available algorithms (for instance, scalable, outerproduct, etc.)

821:    Level: intermediate

823: .seealso: `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`
824: @*/
825: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
826: {
828:   MatCheckProduct(mat, 1);
829:   PetscFree(mat->product->alg);
830:   PetscStrallocpy(alg, &mat->product->alg);
831:   return 0;
832: }

834: /*@
835:    MatProductSetType - Sets a particular matrix product type to be used to compute the given matrix

837:    Collective on mat

839:    Input Parameters:
840: +  mat - the matrix
841: -  productype   - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`.

843:    Level: intermediate

845:    Note:
846:    The small t represents the transpose operation.

848: .seealso: `MatProductCreate()`, `MatProductType`, `MatProductType`,
849:           `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
850: @*/
851: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
852: {
854:   MatCheckProduct(mat, 1);
856:   if (productype != mat->product->type) {
857:     if (mat->product->destroy) (*mat->product->destroy)(mat->product->data);
858:     mat->product->destroy     = NULL;
859:     mat->product->data        = NULL;
860:     mat->ops->productsymbolic = NULL;
861:     mat->ops->productnumeric  = NULL;
862:   }
863:   mat->product->type = productype;
864:   return 0;
865: }

867: /*@
868:    MatProductClear - Clears matrix product internal datastructures.

870:    Collective on mat

872:    Input Parameters:
873: .  mat - the product matrix

875:    Level: intermediate

877:    Notes:
878:    This function should be called to remove any intermediate data used to compute the matrix to free up memory.

880:    After having called this function, matrix-matrix operations can no longer be used on mat

882: .seealso: `MatProductCreate()`
883: @*/
884: PetscErrorCode MatProductClear(Mat mat)
885: {
886:   Mat_Product *product = mat->product;

889:   if (product) {
890:     MatDestroy(&product->A);
891:     MatDestroy(&product->B);
892:     MatDestroy(&product->C);
893:     PetscFree(product->alg);
894:     MatDestroy(&product->Dwork);
895:     if (product->destroy) (*product->destroy)(product->data);
896:   }
897:   PetscFree(mat->product);
898:   mat->ops->productsymbolic = NULL;
899:   mat->ops->productnumeric  = NULL;
900:   return 0;
901: }

903: /* Create a supporting struct and attach it to the matrix product */
904: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
905: {
906:   Mat_Product *product = NULL;

910:   PetscNew(&product);
911:   product->A        = A;
912:   product->B        = B;
913:   product->C        = C;
914:   product->type     = MATPRODUCT_UNSPECIFIED;
915:   product->Dwork    = NULL;
916:   product->api_user = PETSC_FALSE;
917:   product->clear    = PETSC_FALSE;
918:   D->product        = product;

920:   MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT);
921:   MatProductSetFill(D, PETSC_DEFAULT);

923:   PetscObjectReference((PetscObject)A);
924:   PetscObjectReference((PetscObject)B);
925:   PetscObjectReference((PetscObject)C);
926:   return 0;
927: }

929: /*@
930:    MatProductCreateWithMat - Setup a given matrix as a matrix product of other matrices

932:    Collective on Mat

934:    Input Parameters:
935: +  A - the first matrix
936: .  B - the second matrix
937: .  C - the third matrix (optional)
938: -  D - the matrix which will be used to hold the product

940:    Output Parameters:
941: .  D - the product matrix

943:    Notes:
944:    Use `MatProductCreate()` if the matrix you wish computed (the D matrix) does not already exist

946:    See `MatProductCreate()` for details on the usage of the MatProduct routines

948:    Any product data currently attached to D will be cleared

950:    Level: intermediate

952: .seealso: `MatProductCreate()`, `MatProductClear()`
953: @*/
954: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
955: {
958:   MatCheckPreallocated(A, 1);

964:   MatCheckPreallocated(B, 2);

968:   if (C) {
971:     MatCheckPreallocated(C, 3);
974:   }

978:   MatCheckPreallocated(D, 4);

982:   /* Create a supporting struct and attach it to D */
983:   MatProductClear(D);
984:   MatProductCreate_Private(A, B, C, D);
985:   return 0;
986: }

988: /*@
989:    MatProductCreate - create a matrix to hold the result of a matrix-matrix product operation

991:    Collective on A

993:    Input Parameters:
994: +  A - the first matrix
995: .  B - the second matrix
996: -  C - the third matrix (optional)

998:    Output Parameters:
999: .  D - the product matrix

1001:    Level: intermediate

1003:    Example of Usage:
1004: .vb
1005:     MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1006:     MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1007:     MatProductSetAlgorithm(D, alg)
1008:     MatProductSetFill(D,fill)
1009:     MatProductSetFromOptions(D)
1010:     MatProductSymbolic(D)
1011:     MatProductNumeric(D)
1012:     Change numerical values in some of the matrices
1013:     MatProductNumeric(D)
1014: .ve

1016:    Notes:
1017:    Use `MatProductCreateWithMat()` if the matrix you wish computed, the D matrix, already exists.

1019:    The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure

1021:    Developer Note:
1022:    It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1023:    Is there error checking for it?

1025: .seealso: `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`
1026: @*/
1027: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1028: {

1036:   if (C) {
1040:   }

1043:   MatCreate(PetscObjectComm((PetscObject)A), D);
1044:   MatProductCreate_Private(A, B, C, *D);
1045:   return 0;
1046: }

1048: /*
1049:    These are safe basic implementations of ABC, RARt and PtAP
1050:    that do not rely on mat->ops->matmatop function pointers.
1051:    They only use the MatProduct API and are currently used by
1052:    cuSPARSE and KOKKOS-KERNELS backends
1053: */
1054: typedef struct {
1055:   Mat BC;
1056:   Mat ABC;
1057: } MatMatMatPrivate;

1059: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1060: {
1061:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

1063:   MatDestroy(&mmdata->BC);
1064:   MatDestroy(&mmdata->ABC);
1065:   PetscFree(data);
1066:   return 0;
1067: }

1069: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1070: {
1071:   Mat_Product      *product = mat->product;
1072:   MatMatMatPrivate *mmabc;

1074:   MatCheckProduct(mat, 1);
1076:   mmabc = (MatMatMatPrivate *)mat->product->data;
1078:   /* use function pointer directly to prevent logging */
1079:   (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1080:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1081:   mat->product             = mmabc->ABC->product;
1082:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1083:   /* use function pointer directly to prevent logging */
1084:   PetscUseTypeMethod(mat, productnumeric);
1085:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1086:   mat->product             = product;
1087:   return 0;
1088: }

1090: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1091: {
1092:   Mat_Product      *product = mat->product;
1093:   Mat               A, B, C;
1094:   MatProductType    p1, p2;
1095:   MatMatMatPrivate *mmabc;
1096:   const char       *prefix;

1098:   MatCheckProduct(mat, 1);
1100:   MatGetOptionsPrefix(mat, &prefix);
1101:   PetscNew(&mmabc);
1102:   product->data    = mmabc;
1103:   product->destroy = MatDestroy_MatMatMatPrivate;
1104:   switch (product->type) {
1105:   case MATPRODUCT_PtAP:
1106:     p1 = MATPRODUCT_AB;
1107:     p2 = MATPRODUCT_AtB;
1108:     A  = product->B;
1109:     B  = product->A;
1110:     C  = product->B;
1111:     break;
1112:   case MATPRODUCT_RARt:
1113:     p1 = MATPRODUCT_ABt;
1114:     p2 = MATPRODUCT_AB;
1115:     A  = product->B;
1116:     B  = product->A;
1117:     C  = product->B;
1118:     break;
1119:   case MATPRODUCT_ABC:
1120:     p1 = MATPRODUCT_AB;
1121:     p2 = MATPRODUCT_AB;
1122:     A  = product->A;
1123:     B  = product->B;
1124:     C  = product->C;
1125:     break;
1126:   default:
1127:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1128:   }
1129:   MatProductCreate(B, C, NULL, &mmabc->BC);
1130:   MatSetOptionsPrefix(mmabc->BC, prefix);
1131:   MatAppendOptionsPrefix(mmabc->BC, "P1_");
1132:   MatProductSetType(mmabc->BC, p1);
1133:   MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT);
1134:   MatProductSetFill(mmabc->BC, product->fill);
1135:   mmabc->BC->product->api_user = product->api_user;
1136:   MatProductSetFromOptions(mmabc->BC);
1138:   /* use function pointer directly to prevent logging */
1139:   (*mmabc->BC->ops->productsymbolic)(mmabc->BC);

1141:   MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC);
1142:   MatSetOptionsPrefix(mmabc->ABC, prefix);
1143:   MatAppendOptionsPrefix(mmabc->ABC, "P2_");
1144:   MatProductSetType(mmabc->ABC, p2);
1145:   MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT);
1146:   MatProductSetFill(mmabc->ABC, product->fill);
1147:   mmabc->ABC->product->api_user = product->api_user;
1148:   MatProductSetFromOptions(mmabc->ABC);
1150:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1151:   mat->product              = mmabc->ABC->product;
1152:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1153:   /* use function pointer directly to prevent logging */
1154:   PetscUseTypeMethod(mat, productsymbolic);
1155:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1156:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
1157:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
1158:   mat->product                    = product;
1159:   return 0;
1160: }

1162: /*@
1163:    MatProductGetType - Returns the type of matrix-matrix product associated with the given matrix.

1165:    Not collective

1167:    Input Parameter:
1168: .  mat - the matrix

1170:    Output Parameter:
1171: .  mtype - the `MatProductType`

1173:    Level: intermediate

1175: .seealso: `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1176: @*/
1177: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1178: {
1181:   *mtype = MATPRODUCT_UNSPECIFIED;
1182:   if (mat->product) *mtype = mat->product->type;
1183:   return 0;
1184: }

1186: /*@
1187:    MatProductGetMats - Returns the matrices associated with the matrix-matrix product this matrix can receive

1189:    Not collective

1191:    Input Parameter:
1192: .  mat - the product matrix

1194:    Output Parameters:
1195: +  A - the first matrix
1196: .  B - the second matrix
1197: -  C - the third matrix (optional)

1199:    Level: intermediate

1201: .seealso: `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1202: @*/
1203: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1204: {
1206:   if (A) *A = mat->product ? mat->product->A : NULL;
1207:   if (B) *B = mat->product ? mat->product->B : NULL;
1208:   if (C) *C = mat->product ? mat->product->C : NULL;
1209:   return 0;
1210: }