Actual source code: matproduct.c

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

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

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

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

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

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

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

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

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

 51:   PetscFunctionBegin;
 52:   /* AP = A*P */
 53:   PetscCall(MatProductNumeric(AP));
 54:   /* C = P^T*AP */
 55:   PetscCall((*C->ops->transposematmultnumeric)(P, AP, C));
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 66:   PetscCall(PetscInfo((PetscObject)C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
 67:   /* AP = A*P */
 68:   PetscCall(MatProductCreate(A, P, NULL, &AP));
 69:   PetscCall(MatProductSetType(AP, MATPRODUCT_AB));
 70:   PetscCall(MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT));
 71:   PetscCall(MatProductSetFill(AP, fill));
 72:   PetscCall(MatProductSetFromOptions(AP));
 73:   PetscCall(MatProductSymbolic(AP));

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

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

 88:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

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

 97:   PetscFunctionBegin;
 98:   /* RA = R*A */
 99:   PetscCall(MatProductNumeric(RA));
100:   /* C = RA*R^T */
101:   PetscCall((*C->ops->mattransposemultnumeric)(RA, R, C));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

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

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

121:   /* C = RA*R^T */
122:   PetscCall(MatProductSetType(C, MATPRODUCT_ABt));
123:   PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
124:   product->A = RA;
125:   PetscCall(MatProductSetFromOptions(C));
126:   PetscCall(MatProductSymbolic(C));

128:   /* resume user's original input matrix setting for A */
129:   product->A             = A;
130:   product->Dwork         = RA; /* save here so it will be destroyed with product C */
131:   C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
136: {
137:   Mat_Product *product = mat->product;
138:   Mat          A = product->A, BC = product->Dwork;

140:   PetscFunctionBegin;
141:   /* Numeric BC = B*C */
142:   PetscCall(MatProductNumeric(BC));
143:   /* Numeric mat = A*BC */
144:   PetscCall((*mat->ops->matmultnumeric)(A, BC, mat));
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
149: {
150:   Mat_Product *product = mat->product;
151:   Mat          B = product->B, C = product->C, BC;
152:   PetscReal    fill = product->fill;

154:   PetscFunctionBegin;
155:   PetscCall(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));
156:   /* Symbolic BC = B*C */
157:   PetscCall(MatProductCreate(B, C, NULL, &BC));
158:   PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
159:   PetscCall(MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT));
160:   PetscCall(MatProductSetFill(BC, fill));
161:   PetscCall(MatProductSetFromOptions(BC));
162:   PetscCall(MatProductSymbolic(BC));

164:   /* Symbolic mat = A*BC */
165:   PetscCall(MatProductSetType(mat, MATPRODUCT_AB));
166:   PetscCall(MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT));
167:   product->B     = BC;
168:   product->Dwork = BC;
169:   PetscCall(MatProductSetFromOptions(mat));
170:   PetscCall(MatProductSymbolic(mat));

172:   /* resume user's original input matrix setting for B */
173:   product->B               = B;
174:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
179: {
180:   Mat_Product *product = mat->product;

182:   PetscFunctionBegin;
183:   switch (product->type) {
184:   case MATPRODUCT_PtAP:
185:     PetscCall(MatProductSymbolic_PtAP_Unsafe(mat));
186:     break;
187:   case MATPRODUCT_RARt:
188:     PetscCall(MatProductSymbolic_RARt_Unsafe(mat));
189:     break;
190:   case MATPRODUCT_ABC:
191:     PetscCall(MatProductSymbolic_ABC_Unsafe(mat));
192:     break;
193:   default:
194:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
195:   }
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: /*@
200:   MatProductReplaceMats - Replace the input matrices for the matrix-matrix product operation inside the computed matrix

202:   Collective

204:   Input Parameters:
205: + A - the matrix or `NULL` if not being replaced
206: . B - the matrix or `NULL` if not being replaced
207: . C - the matrix or `NULL` if not being replaced
208: - D - the matrix whose values are computed via a matrix-matrix product operation

210:   Level: intermediate

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

218: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic().` `MatProductClear()`
219: @*/
220: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
221: {
222:   Mat_Product *product;
223:   PetscBool    flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;

225:   PetscFunctionBegin;
227:   MatCheckProduct(D, 4);
228:   product = D->product;
229:   if (A) {
231:     PetscCall(PetscObjectReference((PetscObject)A));
232:     PetscCall(PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA));
233:     PetscCall(MatIsSymmetricKnown(A, &isset, &issym));
234:     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 */
235:       flgA                                           = PETSC_FALSE;
236:       product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
237:     }
238:     PetscCall(MatDestroy(&product->A));
239:     product->A = A;
240:   }
241:   if (B) {
243:     PetscCall(PetscObjectReference((PetscObject)B));
244:     PetscCall(PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB));
245:     PetscCall(MatIsSymmetricKnown(B, &isset, &issym));
246:     if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
247:       flgB                                           = PETSC_FALSE;
248:       product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
249:     }
250:     PetscCall(MatDestroy(&product->B));
251:     product->B = B;
252:   }
253:   if (C) {
255:     PetscCall(PetscObjectReference((PetscObject)C));
256:     PetscCall(PetscObjectTypeCompare((PetscObject)product->C, ((PetscObject)C)->type_name, &flgC));
257:     PetscCall(MatIsSymmetricKnown(C, &isset, &issym));
258:     if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
259:       flgC                                           = PETSC_FALSE;
260:       product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
261:     }
262:     PetscCall(MatDestroy(&product->C));
263:     product->C = C;
264:   }
265:   /* Any of the replaced mats is of a different type, reset */
266:   if (!flgA || !flgB || !flgC) {
267:     if (D->product->destroy) PetscCall((*D->product->destroy)(D->product->data));
268:     D->product->destroy = NULL;
269:     D->product->data    = NULL;
270:     if (D->ops->productnumeric || D->ops->productsymbolic) {
271:       PetscCall(MatProductSetFromOptions(D));
272:       PetscCall(MatProductSymbolic(D));
273:     }
274:   }
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
279: {
280:   Mat_Product *product = C->product;
281:   Mat          A = product->A, B = product->B;
282:   PetscInt     k, K              = B->cmap->N;
283:   PetscBool    t = PETSC_TRUE, iscuda = PETSC_FALSE;
284:   PetscBool    Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
285:   char        *Btype = NULL, *Ctype = NULL;

287:   PetscFunctionBegin;
288:   switch (product->type) {
289:   case MATPRODUCT_AB:
290:     t = PETSC_FALSE;
291:   case MATPRODUCT_AtB:
292:     break;
293:   default:
294:     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);
295:   }
296:   if (PetscDefined(HAVE_CUDA)) {
297:     VecType vtype;

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

324:     PetscCall(MatDenseGetColumnVecRead(B, k, &x));
325:     PetscCall(MatDenseGetColumnVecWrite(C, k, &y));
326:     if (t) {
327:       PetscCall(MatMultTranspose(A, x, y));
328:     } else {
329:       PetscCall(MatMult(A, x, y));
330:     }
331:     PetscCall(MatDenseRestoreColumnVecRead(B, k, &x));
332:     PetscCall(MatDenseRestoreColumnVecWrite(C, k, &y));
333:   }
334:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
335:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
336:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
337:   if (PetscDefined(HAVE_CUDA)) {
338:     if (iscuda) {
339:       PetscCall(MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B));
340:       PetscCall(MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C));
341:     } else {
342:       PetscCall(MatBindToCPU(B, Bcpu));
343:       PetscCall(MatBindToCPU(C, Ccpu));
344:     }
345:   }
346:   PetscCall(PetscFree(Btype));
347:   PetscCall(PetscFree(Ctype));
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
352: {
353:   Mat_Product *product = C->product;
354:   Mat          A = product->A, B = product->B;
355:   PetscBool    isdense;

357:   PetscFunctionBegin;
358:   switch (product->type) {
359:   case MATPRODUCT_AB:
360:     PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
361:     break;
362:   case MATPRODUCT_AtB:
363:     PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
364:     break;
365:   default:
366:     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);
367:   }
368:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
369:   if (!isdense) {
370:     PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
371:     /* If matrix type of C was not set or not dense, we need to reset the pointer */
372:     C->ops->productsymbolic = MatProductSymbolic_X_Dense;
373:   }
374:   C->ops->productnumeric = MatProductNumeric_X_Dense;
375:   PetscCall(MatSetUp(C));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: /* a single driver to query the dispatching */
380: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
381: {
382:   Mat_Product      *product = mat->product;
383:   PetscInt          Am, An, Bm, Bn, Cm, Cn;
384:   Mat               A = product->A, B = product->B, C = product->C;
385:   const char *const Bnames[] = {"B", "R", "P"};
386:   const char       *bname;
387:   PetscErrorCode (*fA)(Mat);
388:   PetscErrorCode (*fB)(Mat);
389:   PetscErrorCode (*fC)(Mat);
390:   PetscErrorCode (*f)(Mat) = NULL;

392:   PetscFunctionBegin;
393:   mat->ops->productsymbolic = NULL;
394:   mat->ops->productnumeric  = NULL;
395:   if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(PETSC_SUCCESS);
396:   PetscCheck(A, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing A mat");
397:   PetscCheck(B, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing B mat");
398:   PetscCheck(product->type != MATPRODUCT_ABC || C, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing C mat");
399:   if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
400:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
401:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
402:   else bname = Bnames[0];

404:   /* Check matrices sizes */
405:   Am = A->rmap->N;
406:   An = A->cmap->N;
407:   Bm = B->rmap->N;
408:   Bn = B->cmap->N;
409:   Cm = C ? C->rmap->N : 0;
410:   Cn = C ? C->cmap->N : 0;
411:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
412:     PetscInt t = Bn;
413:     Bn         = Bm;
414:     Bm         = t;
415:   }
416:   if (product->type == MATPRODUCT_AtB) {
417:     PetscInt t = An;
418:     An         = Am;
419:     Am         = t;
420:   }
421:   PetscCheck(An == Bm, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_SIZ, "Matrix dimensions of A and %s are incompatible for MatProductType %s: A %" PetscInt_FMT "x%" PetscInt_FMT ", %s %" PetscInt_FMT "x%" PetscInt_FMT, bname,
422:              MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
423:   PetscCheck(!Cm || Cm == Bn, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_SIZ, "Matrix dimensions of B and C are incompatible for MatProductType %s: B %" PetscInt_FMT "x%" PetscInt_FMT ", C %" PetscInt_FMT "x%" PetscInt_FMT,
424:              MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);

426:   fA = A->ops->productsetfromoptions;
427:   fB = B->ops->productsetfromoptions;
428:   fC = C ? C->ops->productsetfromoptions : fA;
429:   if (C) {
430:     PetscCall(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));
431:   } else {
432:     PetscCall(PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name));
433:   }
434:   if (fA == fB && fA == fC && fA) {
435:     PetscCall(PetscInfo(mat, "  matching op\n"));
436:     PetscCall((*fA)(mat));
437:   }
438:   /* We may have found f but it did not succeed */
439:   if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
440:     char mtypes[256];
441:     PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes)));
442:     PetscCall(PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes)));
443:     PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
444:     PetscCall(PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes)));
445:     if (C) {
446:       PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
447:       PetscCall(PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes)));
448:     }
449:     PetscCall(PetscStrlcat(mtypes, "_C", sizeof(mtypes)));
450: #if defined(__clang__)
451:     PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
452: #elif defined(__GNUC__) || defined(__GNUG__)
453:     PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
454: #endif
455:     PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
456:     PetscCall(PetscInfo(mat, "  querying %s from A? %p\n", mtypes, f));
457:     if (!f) {
458:       PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
459:       PetscCall(PetscInfo(mat, "  querying %s from %s? %p\n", mtypes, bname, f));
460:     }
461:     if (!f && C) {
462:       PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
463:       PetscCall(PetscInfo(mat, "  querying %s from C? %p\n", mtypes, f));
464:     }
465:     if (f) PetscCall((*f)(mat));

467:     /* We may have found f but it did not succeed */
468:     /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
469:     if (!mat->ops->productsymbolic) {
470:       PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes)));
471:       PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
472:       PetscCall(PetscInfo(mat, "  querying %s from A? %p\n", mtypes, f));
473:       if (!f) {
474:         PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
475:         PetscCall(PetscInfo(mat, "  querying %s from %s? %p\n", mtypes, bname, f));
476:       }
477:       if (!f && C) {
478:         PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
479:         PetscCall(PetscInfo(mat, "  querying %s from C? %p\n", mtypes, f));
480:       }
481:     }
482:     if (f) PetscCall((*f)(mat));
483:   }
484:   PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
485:   /* We may have found f but it did not succeed */
486:   if (!mat->ops->productsymbolic) {
487:     /* we can still compute the product if B is of type dense */
488:     if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
489:       PetscBool isdense;

491:       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
492:       if (isdense) {
493:         mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
494:         PetscCall(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 combination 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:       PetscCall(PetscInfo(mat, "  symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n"));
503:       mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
504:     }
505:   }
506:   if (!mat->ops->productsymbolic) PetscCall(PetscInfo(mat, "  symbolic product is not supported\n"));
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

510: /*@
511:   MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product operation where the type,
512:   the algorithm etc are determined from the options database.

514:   Logically Collective

516:   Input Parameter:
517: . mat - the matrix whose values are computed via a matrix-matrix product operation

519:   Options Database Keys:
520: + -mat_product_clear                 - Clear intermediate data structures after `MatProductNumeric()` has been called
521: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm` for possible values
522: - -mat_product_algorithm_backend_cpu - Use the CPU to perform the computation even if the matrix is a GPU matrix

524:   Level: intermediate

526:   Note:
527:   The `-mat_product_clear` option reduces memory usage but means that the matrix cannot be re-used for a matrix-matrix product operation

529: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`,
530:           `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductAlgorithm`
531: @*/
532: PetscErrorCode MatProductSetFromOptions(Mat mat)
533: {
534:   PetscFunctionBegin;
536:   MatCheckProduct(mat, 1);
537:   PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot call MatProductSetFromOptions() with already present data");
538:   mat->product->setfromoptionscalled = PETSC_TRUE;
539:   PetscObjectOptionsBegin((PetscObject)mat);
540:   PetscCall(PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL));
541:   PetscCall(PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()"));
542:   PetscOptionsEnd();
543:   PetscCall(MatProductSetFromOptions_Private(mat));
544:   PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing product after setup phase");
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: /*@
549:   MatProductView - View the private matrix-matrix algorithm object within a matrix

551:   Logically Collective

553:   Input Parameters:
554: + mat    - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
555: - viewer - where the information on the matrix-matrix algorithm of `mat` should be reviewed

557:   Level: intermediate

559:   Developer Note:
560:   Shouldn't this information be printed from an approriate `MatView()` with perhaps certain formats set?

562: .seealso: [](ch_matrices), `MatProductType`, `Mat`, `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
563: @*/
564: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
565: {
566:   PetscFunctionBegin;
568:   if (!mat->product) PetscFunctionReturn(PETSC_SUCCESS);
569:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer));
571:   PetscCheckSameComm(mat, 1, viewer, 2);
572:   if (mat->product->view) PetscCall((*mat->product->view)(mat, viewer));
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

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

583:   PetscFunctionBegin;
584:   PetscCall((*mat->ops->matmultnumeric)(A, B, mat));
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

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

593:   PetscFunctionBegin;
594:   PetscCall((*mat->ops->transposematmultnumeric)(A, B, mat));
595:   PetscFunctionReturn(PETSC_SUCCESS);
596: }

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

603:   PetscFunctionBegin;
604:   PetscCall((*mat->ops->mattransposemultnumeric)(A, B, mat));
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
609: {
610:   Mat_Product *product = mat->product;
611:   Mat          A = product->A, B = product->B;

613:   PetscFunctionBegin;
614:   PetscCall((*mat->ops->ptapnumeric)(A, B, mat));
615:   PetscFunctionReturn(PETSC_SUCCESS);
616: }

618: PetscErrorCode MatProductNumeric_RARt(Mat mat)
619: {
620:   Mat_Product *product = mat->product;
621:   Mat          A = product->A, B = product->B;

623:   PetscFunctionBegin;
624:   PetscCall((*mat->ops->rartnumeric)(A, B, mat));
625:   PetscFunctionReturn(PETSC_SUCCESS);
626: }

628: PetscErrorCode MatProductNumeric_ABC(Mat mat)
629: {
630:   Mat_Product *product = mat->product;
631:   Mat          A = product->A, B = product->B, C = product->C;

633:   PetscFunctionBegin;
634:   PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }

638: /*@
639:   MatProductNumeric - Compute a matrix-matrix product operation with the numerical values

641:   Collective

643:   Input/Output Parameter:
644: . mat - the matrix whose values are computed via a matrix-matrix product operation

646:   Level: intermediate

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

651: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
652: @*/
653: PetscErrorCode MatProductNumeric(Mat mat)
654: {
655:   PetscLogEvent eventtype = -1;

657:   PetscFunctionBegin;
659:   MatCheckProduct(mat, 1);
660:   switch (mat->product->type) {
661:   case MATPRODUCT_AB:
662:     eventtype = MAT_MatMultNumeric;
663:     break;
664:   case MATPRODUCT_AtB:
665:     eventtype = MAT_TransposeMatMultNumeric;
666:     break;
667:   case MATPRODUCT_ABt:
668:     eventtype = MAT_MatTransposeMultNumeric;
669:     break;
670:   case MATPRODUCT_PtAP:
671:     eventtype = MAT_PtAPNumeric;
672:     break;
673:   case MATPRODUCT_RARt:
674:     eventtype = MAT_RARtNumeric;
675:     break;
676:   case MATPRODUCT_ABC:
677:     eventtype = MAT_MatMatMultNumeric;
678:     break;
679:   default:
680:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
681:   }

683:   if (mat->ops->productnumeric) {
684:     PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
685:     PetscUseTypeMethod(mat, productnumeric);
686:     PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
687:   } else if (mat->product) {
688:     char errstr[256];

690:     if (mat->product->type == MATPRODUCT_ABC) {
691:       PetscCall(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));
692:     } else {
693:       PetscCall(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));
694:     }
695:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified numeric phase for product %s", errstr);
696:   }
697:   PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after numeric phase for product");

699:   if (mat->product->clear) PetscCall(MatProductClear(mat));
700:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
701:   PetscFunctionReturn(PETSC_SUCCESS);
702: }

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

711:   PetscFunctionBegin;
712:   PetscCall((*mat->ops->matmultsymbolic)(A, B, product->fill, mat));
713:   mat->ops->productnumeric = MatProductNumeric_AB;
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
718: {
719:   Mat_Product *product = mat->product;
720:   Mat          A = product->A, B = product->B;

722:   PetscFunctionBegin;
723:   PetscCall((*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat));
724:   mat->ops->productnumeric = MatProductNumeric_AtB;
725:   PetscFunctionReturn(PETSC_SUCCESS);
726: }

728: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
729: {
730:   Mat_Product *product = mat->product;
731:   Mat          A = product->A, B = product->B;

733:   PetscFunctionBegin;
734:   PetscCall((*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat));
735:   mat->ops->productnumeric = MatProductNumeric_ABt;
736:   PetscFunctionReturn(PETSC_SUCCESS);
737: }

739: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
740: {
741:   Mat_Product *product = mat->product;
742:   Mat          A = product->A, B = product->B, C = product->C;

744:   PetscFunctionBegin;
745:   PetscCall((*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat));
746:   mat->ops->productnumeric = MatProductNumeric_ABC;
747:   PetscFunctionReturn(PETSC_SUCCESS);
748: }

750: /*@
751:   MatProductSymbolic - Perform the symbolic portion of a matrix-matrix product operation, this creates a data structure for use with the numerical
752:   product to be done with `MatProductNumeric()`

754:   Collective

756:   Input/Output Parameter:
757: . mat - the matrix whose values are to be computed via a matrix-matrix product operation

759:   Level: intermediate

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

764: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
765: @*/
766: PetscErrorCode MatProductSymbolic(Mat mat)
767: {
768:   PetscLogEvent eventtype = -1;
769:   PetscBool     missing   = PETSC_FALSE;

771:   PetscFunctionBegin;
773:   MatCheckProduct(mat, 1);
774:   PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot run symbolic phase. Product data not empty");
775:   switch (mat->product->type) {
776:   case MATPRODUCT_AB:
777:     eventtype = MAT_MatMultSymbolic;
778:     break;
779:   case MATPRODUCT_AtB:
780:     eventtype = MAT_TransposeMatMultSymbolic;
781:     break;
782:   case MATPRODUCT_ABt:
783:     eventtype = MAT_MatTransposeMultSymbolic;
784:     break;
785:   case MATPRODUCT_PtAP:
786:     eventtype = MAT_PtAPSymbolic;
787:     break;
788:   case MATPRODUCT_RARt:
789:     eventtype = MAT_RARtSymbolic;
790:     break;
791:   case MATPRODUCT_ABC:
792:     eventtype = MAT_MatMatMultSymbolic;
793:     break;
794:   default:
795:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
796:   }
797:   mat->ops->productnumeric = NULL;
798:   if (mat->ops->productsymbolic) {
799:     PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
800:     PetscUseTypeMethod(mat, productsymbolic);
801:     PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
802:   } else missing = PETSC_TRUE;

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

807:     if (mat->product->type == MATPRODUCT_ABC) {
808:       PetscCall(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));
809:     } else {
810:       PetscCall(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));
811:     }
812:     PetscCheck(mat->product->setfromoptionscalled, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() first", errstr);
813:     PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unspecified symbolic phase for product %s. The product is not supported", errstr);
814:     PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after symbolic phase for product %s", errstr);
815:   }

817: #if defined(PETSC_HAVE_DEVICE)
818:   Mat       A = mat->product->A;
819:   Mat       B = mat->product->B;
820:   Mat       C = mat->product->C;
821:   PetscBool bindingpropagates;
822:   bindingpropagates = (PetscBool)((A->boundtocpu && A->bindingpropagates) || (B->boundtocpu && B->bindingpropagates));
823:   if (C) bindingpropagates = (PetscBool)(bindingpropagates || (C->boundtocpu && C->bindingpropagates));
824:   if (bindingpropagates) {
825:     PetscCall(MatBindToCPU(mat, PETSC_TRUE));
826:     PetscCall(MatSetBindingPropagates(mat, PETSC_TRUE));
827:   }
828: #endif
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: /*@
833:   MatProductSetFill - Set an expected fill of the matrix whose values are computed via a matrix-matrix product operation

835:   Collective

837:   Input Parameters:
838: + mat  - the matrix whose values are to be computed via a matrix-matrix product operation
839: - 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.
840:           If the product is a dense matrix, this value is not used.

842:   Level: intermediate

844: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
845: @*/
846: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
847: {
848:   PetscFunctionBegin;
850:   MatCheckProduct(mat, 1);
851:   if (fill == (PetscReal)PETSC_DEFAULT || fill == (PetscReal)PETSC_DECIDE) mat->product->fill = 2.0;
852:   else mat->product->fill = fill;
853:   PetscFunctionReturn(PETSC_SUCCESS);
854: }

856: /*@
857:   MatProductSetAlgorithm - Requests a particular algorithm for a matrix-matrix product operation that will perform to compute the given matrix

859:   Collective

861:   Input Parameters:
862: + mat - the matrix whose values are computed via a matrix-matrix product operation
863: - alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.

865:   Options Database Key:
866: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm`

868:   Level: intermediate

870: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductClear()`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`, `MatProductGetAlgorithm()`
871: @*/
872: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
873: {
874:   PetscFunctionBegin;
876:   MatCheckProduct(mat, 1);
877:   PetscCall(PetscFree(mat->product->alg));
878:   PetscCall(PetscStrallocpy(alg, &mat->product->alg));
879:   PetscFunctionReturn(PETSC_SUCCESS);
880: }

882: /*@
883:   MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation

885:   Not Collective

887:   Input Parameter:
888: . mat - the matrix whose values are computed via a matrix-matrix product operation

890:   Output Parameter:
891: . alg - the selected algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.

893:   Level: intermediate

895: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`
896: @*/
897: PetscErrorCode MatProductGetAlgorithm(Mat mat, MatProductAlgorithm *alg)
898: {
899:   PetscFunctionBegin;
901:   PetscAssertPointer(alg, 2);
902:   if (mat->product) *alg = mat->product->alg;
903:   else *alg = NULL;
904:   PetscFunctionReturn(PETSC_SUCCESS);
905: }

907: /*@
908:   MatProductSetType - Sets a particular matrix-matrix product operation to be used to compute the values of the given matrix

910:   Collective

912:   Input Parameters:
913: + mat        - the matrix whose values are computed via a matrix-matrix product operation
914: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`,
915:                   see `MatProductType`

917:   Level: intermediate

919:   Note:
920:   The small t represents the transpose operation.

922: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductType`,
923:           `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
924: @*/
925: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
926: {
927:   PetscFunctionBegin;
929:   MatCheckProduct(mat, 1);
931:   if (productype != mat->product->type) {
932:     if (mat->product->destroy) PetscCall((*mat->product->destroy)(mat->product->data));
933:     mat->product->destroy     = NULL;
934:     mat->product->data        = NULL;
935:     mat->ops->productsymbolic = NULL;
936:     mat->ops->productnumeric  = NULL;
937:   }
938:   mat->product->type = productype;
939:   PetscFunctionReturn(PETSC_SUCCESS);
940: }

942: /*@
943:   MatProductClear - Clears from the matrix any internal data structures related to the computation of the values of the matrix from matrix-matrix product operations

945:   Collective

947:   Input Parameter:
948: . mat - the matrix whose values are to be computed via a matrix-matrix product operation

950:   Options Database Key:
951: . -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called

953:   Level: intermediate

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

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

960: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`
961: @*/
962: PetscErrorCode MatProductClear(Mat mat)
963: {
964:   Mat_Product *product = mat->product;

966:   PetscFunctionBegin;
968:   if (product) {
969:     PetscCall(MatDestroy(&product->A));
970:     PetscCall(MatDestroy(&product->B));
971:     PetscCall(MatDestroy(&product->C));
972:     PetscCall(PetscFree(product->alg));
973:     PetscCall(MatDestroy(&product->Dwork));
974:     if (product->destroy) PetscCall((*product->destroy)(product->data));
975:   }
976:   PetscCall(PetscFree(mat->product));
977:   mat->ops->productsymbolic = NULL;
978:   mat->ops->productnumeric  = NULL;
979:   PetscFunctionReturn(PETSC_SUCCESS);
980: }

982: /* Create a supporting struct and attach it to the matrix product */
983: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
984: {
985:   Mat_Product *product = NULL;

987:   PetscFunctionBegin;
989:   PetscCheck(!D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product already present");
990:   PetscCall(PetscNew(&product));
991:   product->A                    = A;
992:   product->B                    = B;
993:   product->C                    = C;
994:   product->type                 = MATPRODUCT_UNSPECIFIED;
995:   product->Dwork                = NULL;
996:   product->api_user             = PETSC_FALSE;
997:   product->clear                = PETSC_FALSE;
998:   product->setfromoptionscalled = PETSC_FALSE;
999:   D->product                    = product;

1001:   PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
1002:   PetscCall(MatProductSetFill(D, PETSC_DEFAULT));

1004:   PetscCall(PetscObjectReference((PetscObject)A));
1005:   PetscCall(PetscObjectReference((PetscObject)B));
1006:   PetscCall(PetscObjectReference((PetscObject)C));
1007:   PetscFunctionReturn(PETSC_SUCCESS);
1008: }

1010: /*@
1011:   MatProductCreateWithMat - Set a given matrix to have its values computed via matrix-matrix operations on other matrices.

1013:   Collective

1015:   Input Parameters:
1016: + A - the first matrix
1017: . B - the second matrix
1018: . C - the third matrix (optional, use `NULL` if not needed)
1019: - D - the matrix whose values are to be computed via a matrix-matrix product operation

1021:   Level: intermediate

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

1026:   See `MatProductCreate()` for details on the usage of the matrix-matrix product operations

1028:   Any product data currently attached to `D` will be cleared

1030: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductType`, `MatProductSetType()`, `MatProductAlgorithm`,
1031:           `MatProductSetAlgorithm`, `MatProductCreate()`, `MatProductClear()`
1032: @*/
1033: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
1034: {
1035:   PetscFunctionBegin;
1038:   MatCheckPreallocated(A, 1);
1039:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1040:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

1044:   MatCheckPreallocated(B, 2);
1045:   PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1046:   PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

1048:   if (C) {
1051:     MatCheckPreallocated(C, 3);
1052:     PetscCheck(C->assembled, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1053:     PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1054:   }

1058:   MatCheckPreallocated(D, 4);
1059:   PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1060:   PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

1062:   /* Create a supporting struct and attach it to D */
1063:   PetscCall(MatProductClear(D));
1064:   PetscCall(MatProductCreate_Private(A, B, C, D));
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }

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

1071:   Collective

1073:   Input Parameters:
1074: + A - the first matrix
1075: . B - the second matrix
1076: - C - the third matrix (or `NULL`)

1078:   Output Parameter:
1079: . D - the matrix whose values are to be computed via a matrix-matrix product operation

1081:   Level: intermediate

1083:   Example:
1084: .vb
1085:     MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1086:     MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1087:     MatProductSetAlgorithm(D, alg)
1088:     MatProductSetFill(D,fill)
1089:     MatProductSetFromOptions(D)
1090:     MatProductSymbolic(D)
1091:     MatProductNumeric(D)
1092:     Change numerical values in some of the matrices
1093:     MatProductNumeric(D)
1094: .ve

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

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

1101:   Developer Notes:
1102:   It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1103:   Is there error checking for it?

1105: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`
1106: @*/
1107: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1108: {
1109:   PetscFunctionBegin;
1114:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix A");
1115:   PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix B");

1117:   if (C) {
1120:     PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix C");
1121:   }

1123:   PetscAssertPointer(D, 4);
1124:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), D));
1125:   /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1126:   PetscCall(MatProductCreate_Private(A, B, C, *D));
1127:   PetscFunctionReturn(PETSC_SUCCESS);
1128: }

1130: /*
1131:    These are safe basic implementations of ABC, RARt and PtAP
1132:    that do not rely on mat->ops->matmatop function pointers.
1133:    They only use the MatProduct API and are currently used by
1134:    cuSPARSE and KOKKOS-KERNELS backends
1135: */
1136: typedef struct {
1137:   Mat BC;
1138:   Mat ABC;
1139: } MatMatMatPrivate;

1141: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1142: {
1143:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

1145:   PetscFunctionBegin;
1146:   PetscCall(MatDestroy(&mmdata->BC));
1147:   PetscCall(MatDestroy(&mmdata->ABC));
1148:   PetscCall(PetscFree(data));
1149:   PetscFunctionReturn(PETSC_SUCCESS);
1150: }

1152: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1153: {
1154:   Mat_Product      *product = mat->product;
1155:   MatMatMatPrivate *mmabc;

1157:   PetscFunctionBegin;
1158:   MatCheckProduct(mat, 1);
1159:   PetscCheck(mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data empty");
1160:   mmabc = (MatMatMatPrivate *)mat->product->data;
1161:   PetscCheck(mmabc->BC->ops->productnumeric, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing numeric stage");
1162:   /* use function pointer directly to prevent logging */
1163:   PetscCall((*mmabc->BC->ops->productnumeric)(mmabc->BC));
1164:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1165:   mat->product             = mmabc->ABC->product;
1166:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1167:   /* use function pointer directly to prevent logging */
1168:   PetscUseTypeMethod(mat, productnumeric);
1169:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1170:   mat->product             = product;
1171:   PetscFunctionReturn(PETSC_SUCCESS);
1172: }

1174: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1175: {
1176:   Mat_Product      *product = mat->product;
1177:   Mat               A, B, C;
1178:   MatProductType    p1, p2;
1179:   MatMatMatPrivate *mmabc;
1180:   const char       *prefix;

1182:   PetscFunctionBegin;
1183:   MatCheckProduct(mat, 1);
1184:   PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data not empty");
1185:   PetscCall(MatGetOptionsPrefix(mat, &prefix));
1186:   PetscCall(PetscNew(&mmabc));
1187:   product->data    = mmabc;
1188:   product->destroy = MatDestroy_MatMatMatPrivate;
1189:   switch (product->type) {
1190:   case MATPRODUCT_PtAP:
1191:     p1 = MATPRODUCT_AB;
1192:     p2 = MATPRODUCT_AtB;
1193:     A  = product->B;
1194:     B  = product->A;
1195:     C  = product->B;
1196:     break;
1197:   case MATPRODUCT_RARt:
1198:     p1 = MATPRODUCT_ABt;
1199:     p2 = MATPRODUCT_AB;
1200:     A  = product->B;
1201:     B  = product->A;
1202:     C  = product->B;
1203:     break;
1204:   case MATPRODUCT_ABC:
1205:     p1 = MATPRODUCT_AB;
1206:     p2 = MATPRODUCT_AB;
1207:     A  = product->A;
1208:     B  = product->B;
1209:     C  = product->C;
1210:     break;
1211:   default:
1212:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1213:   }
1214:   PetscCall(MatProductCreate(B, C, NULL, &mmabc->BC));
1215:   PetscCall(MatSetOptionsPrefix(mmabc->BC, prefix));
1216:   PetscCall(MatAppendOptionsPrefix(mmabc->BC, "P1_"));
1217:   PetscCall(MatProductSetType(mmabc->BC, p1));
1218:   PetscCall(MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT));
1219:   PetscCall(MatProductSetFill(mmabc->BC, product->fill));
1220:   mmabc->BC->product->api_user = product->api_user;
1221:   PetscCall(MatProductSetFromOptions(mmabc->BC));
1222:   PetscCheck(mmabc->BC->ops->productsymbolic, 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);
1223:   /* use function pointer directly to prevent logging */
1224:   PetscCall((*mmabc->BC->ops->productsymbolic)(mmabc->BC));

1226:   PetscCall(MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC));
1227:   PetscCall(MatSetOptionsPrefix(mmabc->ABC, prefix));
1228:   PetscCall(MatAppendOptionsPrefix(mmabc->ABC, "P2_"));
1229:   PetscCall(MatProductSetType(mmabc->ABC, p2));
1230:   PetscCall(MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT));
1231:   PetscCall(MatProductSetFill(mmabc->ABC, product->fill));
1232:   mmabc->ABC->product->api_user = product->api_user;
1233:   PetscCall(MatProductSetFromOptions(mmabc->ABC));
1234:   PetscCheck(mmabc->ABC->ops->productsymbolic, 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);
1235:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1236:   mat->product              = mmabc->ABC->product;
1237:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1238:   /* use function pointer directly to prevent logging */
1239:   PetscUseTypeMethod(mat, productsymbolic);
1240:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1241:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
1242:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
1243:   mat->product                    = product;
1244:   PetscFunctionReturn(PETSC_SUCCESS);
1245: }

1247: /*@
1248:   MatProductGetType - Returns the type of matrix-matrix product associated with computing values for the given matrix

1250:   Not Collective

1252:   Input Parameter:
1253: . mat - the matrix whose values are to be computed via a matrix-matrix product operation

1255:   Output Parameter:
1256: . mtype - the `MatProductType`

1258:   Level: intermediate

1260: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1261: @*/
1262: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1263: {
1264:   PetscFunctionBegin;
1266:   PetscAssertPointer(mtype, 2);
1267:   *mtype = MATPRODUCT_UNSPECIFIED;
1268:   if (mat->product) *mtype = mat->product->type;
1269:   PetscFunctionReturn(PETSC_SUCCESS);
1270: }

1272: /*@
1273:   MatProductGetMats - Returns the matrices associated with the matrix-matrix product associated with computing values for the given matrix

1275:   Not Collective

1277:   Input Parameter:
1278: . mat - the matrix whose values are to be computed via a matrix-matrix product operation

1280:   Output Parameters:
1281: + A - the first matrix
1282: . B - the second matrix
1283: - C - the third matrix (may be `NULL` for some `MatProductType`)

1285:   Level: intermediate

1287: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1288: @*/
1289: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1290: {
1291:   PetscFunctionBegin;
1293:   if (A) *A = mat->product ? mat->product->A : NULL;
1294:   if (B) *B = mat->product ? mat->product->B : NULL;
1295:   if (C) *C = mat->product ? mat->product->C : NULL;
1296:   PetscFunctionReturn(PETSC_SUCCESS);
1297: }