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) An = Am;

    418:   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,
    419:              MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
    420:   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,
    421:              MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);

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

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

    488:       PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
    489:       if (isdense) {
    490:         mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
    491:         PetscCall(PetscInfo(mat, "  using basic looping over columns of a dense matrix\n"));
    492:       }
    493:     } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
    494:       /*
    495:          TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
    496:                the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
    497:                before computing the symbolic phase
    498:       */
    499:       PetscCall(PetscInfo(mat, "  symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n"));
    500:       mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
    501:     }
    502:   }
    503:   if (!mat->ops->productsymbolic) PetscCall(PetscInfo(mat, "  symbolic product is not supported\n"));
    504:   PetscFunctionReturn(PETSC_SUCCESS);
    505: }

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

    511:   Logically Collective

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

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

    521:   Level: intermediate

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

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

    545: /*@
    546:   MatProductView - View the private matrix-matrix algorithm object within a matrix

    548:   Logically Collective

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

    554:   Level: intermediate

    556:   Developer Note:
    557:   Shouldn't this information be printed from an appropriate `MatView()` with perhaps certain formats set?

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

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

    580:   PetscFunctionBegin;
    581:   PetscCall((*mat->ops->matmultnumeric)(A, B, mat));
    582:   PetscFunctionReturn(PETSC_SUCCESS);
    583: }

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

    590:   PetscFunctionBegin;
    591:   PetscCall((*mat->ops->transposematmultnumeric)(A, B, mat));
    592:   PetscFunctionReturn(PETSC_SUCCESS);
    593: }

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

    600:   PetscFunctionBegin;
    601:   PetscCall((*mat->ops->mattransposemultnumeric)(A, B, mat));
    602:   PetscFunctionReturn(PETSC_SUCCESS);
    603: }

    605: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
    606: {
    607:   Mat_Product *product = mat->product;
    608:   Mat          A = product->A, B = product->B;

    610:   PetscFunctionBegin;
    611:   PetscCall((*mat->ops->ptapnumeric)(A, B, mat));
    612:   PetscFunctionReturn(PETSC_SUCCESS);
    613: }

    615: PetscErrorCode MatProductNumeric_RARt(Mat mat)
    616: {
    617:   Mat_Product *product = mat->product;
    618:   Mat          A = product->A, B = product->B;

    620:   PetscFunctionBegin;
    621:   PetscCall((*mat->ops->rartnumeric)(A, B, mat));
    622:   PetscFunctionReturn(PETSC_SUCCESS);
    623: }

    625: PetscErrorCode MatProductNumeric_ABC(Mat mat)
    626: {
    627:   Mat_Product *product = mat->product;
    628:   Mat          A = product->A, B = product->B, C = product->C;

    630:   PetscFunctionBegin;
    631:   PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
    632:   PetscFunctionReturn(PETSC_SUCCESS);
    633: }

    635: /*@
    636:   MatProductNumeric - Compute a matrix-matrix product operation with the numerical values

    638:   Collective

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

    643:   Level: intermediate

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

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

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

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

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

    696:   if (mat->product->clear) PetscCall(MatProductClear(mat));
    697:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
    698:   PetscFunctionReturn(PETSC_SUCCESS);
    699: }

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

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

    714: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
    715: {
    716:   Mat_Product *product = mat->product;
    717:   Mat          A = product->A, B = product->B;

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

    725: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
    726: {
    727:   Mat_Product *product = mat->product;
    728:   Mat          A = product->A, B = product->B;

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

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

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

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

    751:   Collective

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

    756:   Level: intermediate

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

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

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

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

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

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

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

    832:   Collective

    834:   Input Parameters:
    835: + mat  - the matrix whose values are to be computed via a matrix-matrix product operation
    836: - fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use `PETSC_DETERMINE` or `PETSC_CURRENT` if you do not have a good estimate.
    837:          If the product is a dense matrix, this value is not used.

    839:   Level: intermediate

    841:   Notes:
    842:   Use `fill` of `PETSC_DETERMINE` to use the default value.

    844:   The deprecated `PETSC_DEFAULT` is also supported to mean use the current value.

    846: .seealso: [](ch_matrices), `MatProduct`, `PETSC_DETERMINE`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
    847: @*/
    848: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
    849: {
    850:   PetscFunctionBegin;
    852:   MatCheckProduct(mat, 1);
    853:   if (fill == (PetscReal)PETSC_DETERMINE) mat->product->fill = mat->product->default_fill;
    854:   else if (fill != (PetscReal)PETSC_CURRENT) mat->product->fill = fill;
    855:   PetscFunctionReturn(PETSC_SUCCESS);
    856: }

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

    861:   Collective

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

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

    870:   Level: intermediate

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

    884: /*@
    885:   MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation

    887:   Not Collective

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

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

    895:   Level: intermediate

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

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

    912:   Collective

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

    919:   Level: intermediate

    921:   Note:
    922:   The small t represents the transpose operation.

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

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

    947:   Collective

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

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

    955:   Level: intermediate

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

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

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

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

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

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

    1004:   PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
    1005:   PetscCall(MatProductSetFill(D, PETSC_DEFAULT));

    1007:   PetscCall(PetscObjectReference((PetscObject)A));
    1008:   PetscCall(PetscObjectReference((PetscObject)B));
    1009:   PetscCall(PetscObjectReference((PetscObject)C));
    1010:   PetscFunctionReturn(PETSC_SUCCESS);
    1011: }

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

    1016:   Collective

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

    1024:   Level: intermediate

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

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

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

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

    1047:   MatCheckPreallocated(B, 2);
    1048:   PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1049:   PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

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

    1061:   MatCheckPreallocated(D, 4);
    1062:   PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1063:   PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

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

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

    1074:   Collective

    1076:   Input Parameters:
    1077: + A - the first matrix
    1078: . B - the second matrix
    1079: - C - the third matrix (or `NULL`)

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

    1084:   Level: intermediate

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

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

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

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

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

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

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

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

    1144: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
    1145: {
    1146:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

    1148:   PetscFunctionBegin;
    1149:   PetscCall(MatDestroy(&mmdata->BC));
    1150:   PetscCall(MatDestroy(&mmdata->ABC));
    1151:   PetscCall(PetscFree(data));
    1152:   PetscFunctionReturn(PETSC_SUCCESS);
    1153: }

    1155: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
    1156: {
    1157:   Mat_Product      *product = mat->product;
    1158:   MatMatMatPrivate *mmabc;

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

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

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

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

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

    1253:   Not Collective

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

    1258:   Output Parameter:
    1259: . mtype - the `MatProductType`

    1261:   Level: intermediate

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

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

    1278:   Not Collective

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

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

    1288:   Level: intermediate

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