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: `Mat`, `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: `Mat`, `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

529:    Input Parameters:
530: +  mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
531: -  viewer - where `mat` should be reviewed

533:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

605: /* ----------------------------------------------- */

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

610:    Collective

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

615:    Level: intermediate

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

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

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

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

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

661:     if (mat->product->type == MATPRODUCT_ABC) {
662:       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);
663:     } else {
664:       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);
665:     }
668:   }

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

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

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

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

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

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

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

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

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

718: /* ----------------------------------------------- */

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

724:    Collective

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

729:    Level: intermediate

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

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

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

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

776:     if (mat->product->type == MATPRODUCT_ABC) {
777:       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);
778:     } else {
779:       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);
780:     }
783:   }
784:   return 0;
785: }

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

790:    Collective on Mat

792:    Input Parameters:
793: +  mat - the matrix product result matrix
794: -  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.

796:    Level: intermediate

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

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

812:    Collective

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

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

822:    Level: intermediate

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

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

838:    Collective

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

844:    Level: intermediate

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

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

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

871:    Collective

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

876:    Level: intermediate

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

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

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

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

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

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

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

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

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

933:    Collective on Mat

935:    Input Parameters:
936: +  A - the first matrix
937: .  B - the second matrix
938: .  C - the third matrix (optional, use `NULL` if not needed)
939: -  D - the matrix which will be used to hold the product

941:    Level: intermediate

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: .seealso: `Mat`, `MatProductCreate()`, `MatProductClear()`
951: @*/
952: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
953: {
956:   MatCheckPreallocated(A, 1);

962:   MatCheckPreallocated(B, 2);

966:   if (C) {
969:     MatCheckPreallocated(C, 3);
972:   }

976:   MatCheckPreallocated(D, 4);

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

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

989:    Collective on A

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

996:    Output Parameters:
997: .  D - the product matrix

999:    Level: intermediate

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

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

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

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

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

1034:   if (C) {
1038:   }

1041:   MatCreate(PetscObjectComm((PetscObject)A), D);
1042:   /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1043:   MatProductCreate_Private(A, B, C, *D);
1044:   return 0;
1045: }

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

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

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

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

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

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

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

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

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

1164:    Not collective

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

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

1172:    Level: intermediate

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

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

1188:    Not collective

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

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

1198:    Level: intermediate

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