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:   product->type = MATPRODUCT_AtB;
     56:   PetscCall((*C->ops->transposematmultnumeric)(P, AP, C));
     57:   product->type = MATPRODUCT_PtAP;
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

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

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

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

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

     91:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }

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

    100:   PetscFunctionBegin;
    101:   /* RA = R*A */
    102:   PetscCall(MatProductNumeric(RA));
    103:   /* C = RA*R^T */
    104:   product->type = MATPRODUCT_ABt;
    105:   PetscCall((*C->ops->mattransposemultnumeric)(RA, R, C));
    106:   product->type = MATPRODUCT_RARt;
    107:   PetscFunctionReturn(PETSC_SUCCESS);
    108: }

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

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

    126:   /* C = RA*R^T */
    127:   PetscCall(MatProductSetType(C, MATPRODUCT_ABt));
    128:   PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
    129:   product->A = RA;
    130:   PetscCall(MatProductSetFromOptions(C));
    131:   PetscCall(MatProductSymbolic(C));

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

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

    146:   PetscFunctionBegin;
    147:   /* Numeric BC = B*C */
    148:   PetscCall(MatProductNumeric(BC));
    149:   /* Numeric mat = A*BC */
    150:   product->type = MATPRODUCT_AB;
    151:   PetscCall((*mat->ops->matmultnumeric)(A, BC, mat));
    152:   product->type = MATPRODUCT_ABC;
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }

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

    162:   PetscFunctionBegin;
    163:   PetscCall(PetscInfo(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));
    164:   /* Symbolic BC = B*C */
    165:   PetscCall(MatProductCreate(B, C, NULL, &BC));
    166:   PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
    167:   PetscCall(MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT));
    168:   PetscCall(MatProductSetFill(BC, fill));
    169:   PetscCall(MatProductSetFromOptions(BC));
    170:   PetscCall(MatProductSymbolic(BC));

    172:   /* Symbolic mat = A*BC */
    173:   PetscCall(MatProductSetType(mat, MATPRODUCT_AB));
    174:   PetscCall(MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT));
    175:   product->B     = BC;
    176:   product->Dwork = BC;
    177:   PetscCall(MatProductSetFromOptions(mat));
    178:   PetscCall(MatProductSymbolic(mat));

    180:   /* resume user's original input matrix setting for B */
    181:   product->type            = MATPRODUCT_ABC;
    182:   product->B               = B;
    183:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

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

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

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

    211:   Collective

    213:   Input Parameters:
    214: + A - the matrix or `NULL` if not being replaced
    215: . B - the matrix or `NULL` if not being replaced
    216: . C - the matrix or `NULL` if not being replaced
    217: - D - the matrix whose values are computed via a matrix-matrix product operation

    219:   Level: intermediate

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

    227: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic()`, `MatProductClear()`
    228: @*/
    229: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
    230: {
    231:   Mat_Product *product;
    232:   PetscBool    flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;

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

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

    296:   PetscFunctionBegin;
    297:   switch (product->type) {
    298:   case MATPRODUCT_AB:
    299:     t = PETSC_FALSE;
    300:   case MATPRODUCT_AtB:
    301:     break;
    302:   default:
    303:     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);
    304:   }
    305:   if (PetscDefined(HAVE_CUDA)) {
    306:     VecType vtype;

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

    333:     PetscCall(MatDenseGetColumnVecRead(B, k, &x));
    334:     PetscCall(MatDenseGetColumnVecWrite(C, k, &y));
    335:     if (t) {
    336:       PetscCall(MatMultTranspose(A, x, y));
    337:     } else {
    338:       PetscCall(MatMult(A, x, y));
    339:     }
    340:     PetscCall(MatDenseRestoreColumnVecRead(B, k, &x));
    341:     PetscCall(MatDenseRestoreColumnVecWrite(C, k, &y));
    342:   }
    343:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
    344:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    345:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    346:   if (PetscDefined(HAVE_CUDA)) {
    347:     if (iscuda) {
    348:       PetscCall(MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B));
    349:       PetscCall(MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C));
    350:     } else {
    351:       PetscCall(MatBindToCPU(B, Bcpu));
    352:       PetscCall(MatBindToCPU(C, Ccpu));
    353:     }
    354:   }
    355:   PetscCall(PetscFree(Btype));
    356:   PetscCall(PetscFree(Ctype));
    357:   PetscFunctionReturn(PETSC_SUCCESS);
    358: }

    360: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
    361: {
    362:   Mat_Product *product = C->product;
    363:   Mat          A = product->A, B = product->B;
    364:   PetscBool    isdense;

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

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

    401:   PetscFunctionBegin;
    402:   mat->ops->productsymbolic = NULL;
    403:   mat->ops->productnumeric  = NULL;
    404:   if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(PETSC_SUCCESS);
    405:   PetscCheck(A, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing A mat");
    406:   PetscCheck(B, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing B mat");
    407:   PetscCheck(product->type != MATPRODUCT_ABC || C, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing C mat");
    408:   if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
    409:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
    410:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
    411:   else bname = Bnames[0];

    413:   /* Check matrices sizes */
    414:   Am = A->rmap->N;
    415:   An = A->cmap->N;
    416:   Bm = B->rmap->N;
    417:   Bn = B->cmap->N;
    418:   Cm = C ? C->rmap->N : 0;
    419:   Cn = C ? C->cmap->N : 0;
    420:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
    421:     PetscInt t = Bn;
    422:     Bn         = Bm;
    423:     Bm         = t;
    424:   }
    425:   if (product->type == MATPRODUCT_AtB) An = Am;

    427:   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,
    428:              MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
    429:   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,
    430:              MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);

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

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

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

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

    520:   Logically Collective

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

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

    530:   Level: intermediate

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

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

    554: /*@
    555:   MatProductView - View the private matrix-matrix algorithm object within a matrix

    557:   Logically Collective

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

    563:   Level: intermediate

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

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

    582: /* these are basic implementations relying on the old function pointers
    583:  * they are dangerous and should be removed in the future */
    584: PetscErrorCode MatProductNumeric_AB(Mat mat)
    585: {
    586:   Mat_Product *product = mat->product;
    587:   Mat          A = product->A, B = product->B;

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

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

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

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

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

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

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

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

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

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

    639:   PetscFunctionBegin;
    640:   PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }

    644: /*@
    645:   MatProductNumeric - Compute a matrix-matrix product operation with the numerical values

    647:   Collective

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

    652:   Level: intermediate

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

    657: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
    658: @*/
    659: PetscErrorCode MatProductNumeric(Mat mat)
    660: {
    661:   PetscLogEvent eventtype = -1;

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

    689:   if (mat->ops->productnumeric) {
    690:     PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
    691:     PetscUseTypeMethod(mat, productnumeric);
    692:     PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
    693:   } else if (mat->product) {
    694:     char errstr[256];

    696:     if (mat->product->type == MATPRODUCT_ABC) {
    697:       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));
    698:     } else {
    699:       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));
    700:     }
    701:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified numeric phase for product %s", errstr);
    702:   }
    703:   PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after numeric phase for product");

    705:   if (mat->product->clear) PetscCall(MatProductClear(mat));
    706:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
    707:   PetscFunctionReturn(PETSC_SUCCESS);
    708: }

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

    717:   PetscFunctionBegin;
    718:   PetscCall((*mat->ops->matmultsymbolic)(A, B, product->fill, mat));
    719:   mat->ops->productnumeric = MatProductNumeric_AB;
    720:   PetscFunctionReturn(PETSC_SUCCESS);
    721: }

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

    728:   PetscFunctionBegin;
    729:   PetscCall((*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat));
    730:   mat->ops->productnumeric = MatProductNumeric_AtB;
    731:   PetscFunctionReturn(PETSC_SUCCESS);
    732: }

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

    739:   PetscFunctionBegin;
    740:   PetscCall((*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat));
    741:   mat->ops->productnumeric = MatProductNumeric_ABt;
    742:   PetscFunctionReturn(PETSC_SUCCESS);
    743: }

    745: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
    746: {
    747:   Mat_Product *product = mat->product;
    748:   Mat          A = product->A, B = product->B, C = product->C;

    750:   PetscFunctionBegin;
    751:   PetscCall((*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat));
    752:   mat->ops->productnumeric = MatProductNumeric_ABC;
    753:   PetscFunctionReturn(PETSC_SUCCESS);
    754: }

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

    760:   Collective

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

    765:   Level: intermediate

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

    770: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
    771: @*/
    772: PetscErrorCode MatProductSymbolic(Mat mat)
    773: {
    774:   PetscLogEvent eventtype = -1;
    775:   PetscBool     missing   = PETSC_FALSE;
    776:   Mat_Product  *product   = mat->product;
    777:   Mat           A         = product->A;
    778:   Mat           B         = product->B;
    779:   Mat           C         = product->C;

    781:   PetscFunctionBegin;
    783:   MatCheckProduct(mat, 1);
    784:   PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot run symbolic phase. Product data not empty");
    785:   switch (mat->product->type) {
    786:   case MATPRODUCT_AB:
    787:     eventtype = MAT_MatMultSymbolic;
    788:     break;
    789:   case MATPRODUCT_AtB:
    790:     eventtype = MAT_TransposeMatMultSymbolic;
    791:     break;
    792:   case MATPRODUCT_ABt:
    793:     eventtype = MAT_MatTransposeMultSymbolic;
    794:     break;
    795:   case MATPRODUCT_PtAP:
    796:     eventtype = MAT_PtAPSymbolic;
    797:     break;
    798:   case MATPRODUCT_RARt:
    799:     eventtype = MAT_RARtSymbolic;
    800:     break;
    801:   case MATPRODUCT_ABC:
    802:     eventtype = MAT_MatMatMultSymbolic;
    803:     break;
    804:   default:
    805:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
    806:   }
    807:   mat->ops->productnumeric = NULL;
    808:   if (mat->ops->productsymbolic) {
    809:     PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
    810:     PetscUseTypeMethod(mat, productsymbolic);
    811:     PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
    812:   } else missing = PETSC_TRUE;
    813:   if (missing || !mat->product || !mat->ops->productnumeric) {
    814:     char errstr[256];

    816:     if (mat->product->type == MATPRODUCT_ABC) {
    817:       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));
    818:     } else {
    819:       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));
    820:     }
    821:     PetscCheck(mat->product->setfromoptionscalled, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() first", errstr);
    822:     PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unspecified symbolic phase for product %s. The product is not supported", errstr);
    823:     PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after symbolic phase for product %s", errstr);
    824:   }
    825: #if defined(PETSC_HAVE_DEVICE)
    826:   PetscBool bindingpropagates;
    827:   bindingpropagates = (PetscBool)((A->boundtocpu && A->bindingpropagates) || (B->boundtocpu && B->bindingpropagates));
    828:   if (C) bindingpropagates = (PetscBool)(bindingpropagates || (C->boundtocpu && C->bindingpropagates));
    829:   if (bindingpropagates) {
    830:     PetscCall(MatBindToCPU(mat, PETSC_TRUE));
    831:     PetscCall(MatSetBindingPropagates(mat, PETSC_TRUE));
    832:   }
    833: #endif
    834:   /* set block sizes */
    835:   switch (product->type) {
    836:   case MATPRODUCT_PtAP:
    837:     if (B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->cmap->bs, B->cmap->bs));
    838:     break;
    839:   case MATPRODUCT_RARt:
    840:     if (B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->rmap->bs, B->rmap->bs));
    841:     break;
    842:   case MATPRODUCT_ABC:
    843:     PetscCall(MatSetBlockSizesFromMats(mat, A, C));
    844:     break;
    845:   case MATPRODUCT_AB:
    846:     PetscCall(MatSetBlockSizesFromMats(mat, A, B));
    847:     break;
    848:   case MATPRODUCT_AtB:
    849:     if (A->cmap->bs > 1 || B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, B->cmap->bs));
    850:     break;
    851:   case MATPRODUCT_ABt:
    852:     if (A->rmap->bs > 1 || B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, B->rmap->bs));
    853:     break;
    854:   default:
    855:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
    856:   }
    857:   PetscFunctionReturn(PETSC_SUCCESS);
    858: }

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

    863:   Collective

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

    870:   Level: intermediate

    872:   Notes:
    873:   Use `fill` of `PETSC_DETERMINE` to use the default value.

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

    877: .seealso: [](ch_matrices), `MatProduct`, `PETSC_DETERMINE`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
    878: @*/
    879: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
    880: {
    881:   PetscFunctionBegin;
    883:   MatCheckProduct(mat, 1);
    884:   if (fill == (PetscReal)PETSC_DETERMINE) mat->product->fill = mat->product->default_fill;
    885:   else if (fill != (PetscReal)PETSC_CURRENT) mat->product->fill = fill;
    886:   PetscFunctionReturn(PETSC_SUCCESS);
    887: }

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

    892:   Collective

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

    898:   Options Database Key:
    899: . -mat_product_algorithm algorithm - Sets the algorithm, see `MatProductAlgorithm`

    901:   Level: intermediate

    903: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductClear()`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`, `MatProductGetAlgorithm()`
    904: @*/
    905: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
    906: {
    907:   PetscFunctionBegin;
    909:   MatCheckProduct(mat, 1);
    910:   PetscCall(PetscFree(mat->product->alg));
    911:   PetscCall(PetscStrallocpy(alg, &mat->product->alg));
    912:   PetscFunctionReturn(PETSC_SUCCESS);
    913: }

    915: /*@
    916:   MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation

    918:   Not Collective

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

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

    926:   Level: intermediate

    928: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`
    929: @*/
    930: PetscErrorCode MatProductGetAlgorithm(Mat mat, MatProductAlgorithm *alg)
    931: {
    932:   PetscFunctionBegin;
    934:   PetscAssertPointer(alg, 2);
    935:   if (mat->product) *alg = mat->product->alg;
    936:   else *alg = NULL;
    937:   PetscFunctionReturn(PETSC_SUCCESS);
    938: }

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

    943:   Collective

    945:   Input Parameters:
    946: + mat        - the matrix whose values are computed via a matrix-matrix product operation
    947: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`,
    948:                see `MatProductType`

    950:   Level: intermediate

    952:   Note:
    953:   The small t represents the transpose operation.

    955: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductType`,
    956:           `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
    957: @*/
    958: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
    959: {
    960:   PetscFunctionBegin;
    962:   MatCheckProduct(mat, 1);
    964:   if (productype != mat->product->type) {
    965:     if (mat->product->destroy) PetscCall((*mat->product->destroy)(&mat->product->data));
    966:     mat->product->destroy     = NULL;
    967:     mat->product->data        = NULL;
    968:     mat->ops->productsymbolic = NULL;
    969:     mat->ops->productnumeric  = NULL;
    970:   }
    971:   mat->product->type = productype;
    972:   PetscFunctionReturn(PETSC_SUCCESS);
    973: }

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

    978:   Collective

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

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

    986:   Level: intermediate

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

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

    993:   Developer Note:
    994:   This frees the `Mat_Product` context that was attached to the matrix during `MatProductCreate()` or `MatProductCreateWithMat()`

    996: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`
    997: @*/
    998: PetscErrorCode MatProductClear(Mat mat)
    999: {
    1000:   Mat_Product *product = mat->product;

    1002:   PetscFunctionBegin;
    1004:   if (product) {
    1005:     PetscCall(MatDestroy(&product->A));
    1006:     PetscCall(MatDestroy(&product->B));
    1007:     PetscCall(MatDestroy(&product->C));
    1008:     PetscCall(PetscFree(product->alg));
    1009:     PetscCall(MatDestroy(&product->Dwork));
    1010:     if (product->destroy) PetscCall((*product->destroy)(&product->data));
    1011:   }
    1012:   PetscCall(PetscFree(mat->product));
    1013:   mat->ops->productsymbolic = NULL;
    1014:   mat->ops->productnumeric  = NULL;
    1015:   PetscFunctionReturn(PETSC_SUCCESS);
    1016: }

    1018: /* Create a supporting struct and attach it to the matrix product */
    1019: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
    1020: {
    1021:   Mat_Product *product = NULL;

    1023:   PetscFunctionBegin;
    1025:   PetscCheck(!D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product already present");
    1026:   PetscCall(PetscNew(&product));
    1027:   product->A                    = A;
    1028:   product->B                    = B;
    1029:   product->C                    = C;
    1030:   product->type                 = MATPRODUCT_UNSPECIFIED;
    1031:   product->Dwork                = NULL;
    1032:   product->api_user             = PETSC_FALSE;
    1033:   product->clear                = PETSC_FALSE;
    1034:   product->setfromoptionscalled = PETSC_FALSE;
    1035:   PetscObjectParameterSetDefault(product, fill, 2);
    1036:   D->product = product;

    1038:   PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
    1039:   PetscCall(MatProductSetFill(D, PETSC_DEFAULT));

    1041:   PetscCall(PetscObjectReference((PetscObject)A));
    1042:   PetscCall(PetscObjectReference((PetscObject)B));
    1043:   PetscCall(PetscObjectReference((PetscObject)C));
    1044:   PetscFunctionReturn(PETSC_SUCCESS);
    1045: }

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

    1050:   Collective

    1052:   Input Parameters:
    1053: + A - the first matrix
    1054: . B - the second matrix
    1055: . C - the third matrix (optional, use `NULL` if not needed)
    1056: - D - the matrix whose values are to be computed via a matrix-matrix product operation

    1058:   Level: intermediate

    1060:   Notes:
    1061:   Use `MatProductCreate()` if the matrix you wish computed `D` does not exist

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

    1065:   Any product data currently attached to `D` will be freed

    1067: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductType`, `MatProductSetType()`, `MatProductAlgorithm`,
    1068:           `MatProductSetAlgorithm`, `MatProductCreate()`, `MatProductClear()`
    1069: @*/
    1070: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
    1071: {
    1072:   PetscFunctionBegin;
    1075:   MatCheckPreallocated(A, 1);
    1076:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1077:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    1081:   MatCheckPreallocated(B, 2);
    1082:   PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1083:   PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    1085:   if (C) {
    1088:     MatCheckPreallocated(C, 3);
    1089:     PetscCheck(C->assembled, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1090:     PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
    1091:   }

    1095:   MatCheckPreallocated(D, 4);
    1096:   PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    1097:   PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

    1099:   /* Create a supporting struct and attach it to D */
    1100:   PetscCall(MatProductClear(D));
    1101:   PetscCall(MatProductCreate_Private(A, B, C, D));
    1102:   PetscFunctionReturn(PETSC_SUCCESS);
    1103: }

    1105: /*@
    1106:   MatProductCreate - create a matrix to hold the result of a matrix-matrix (or matrix-matrix-matrix) product operation

    1108:   Collective

    1110:   Input Parameters:
    1111: + A - the first matrix
    1112: . B - the second matrix
    1113: - C - the third matrix (or `NULL`)

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

    1118:   Level: intermediate

    1120:   Example:
    1121: .vb
    1122:     MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
    1123:     MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
    1124:     MatProductSetAlgorithm(D, alg)
    1125:     MatProductSetFill(D,fill)
    1126:     MatProductSetFromOptions(D)
    1127:     MatProductSymbolic(D)
    1128:     MatProductNumeric(D)
    1129:     Change numerical values in some of the matrices
    1130:     MatProductNumeric(D)
    1131: .ve

    1133:   Notes:
    1134:   Use `MatProductCreateWithMat()` if `D` the matrix you wish computed already exists.

    1136:   The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure of the input matrices.

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

    1142:   On this call, auxiliary data needed to compute the product is stored in `D` in a `Mat_Product` context. A call to `MatProductClear()` frees this
    1143:   information.

    1145:   Each `MatProductAlgorithm` associated with a particular `MatType` stores additional data needed for the product computation
    1146:   (generally this data is computed in `MatProductSymbolic()`) inside the `Mat_Product` context in a `MatProductCtx_XXX` data structure
    1147:   and provides a `MatProductCtxDestroy_XXX()` routine to free that data. The `MatProductAlgorithm` and `MatType` specific destroy routine is called by
    1148:   `MatProductClear()`.

    1150: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`,
    1151:           `MatProductSymbolic()`, `MatProductNumeric()`, `MatProductAlgorithm`, `MatProductType`
    1152: @*/
    1153: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
    1154: {
    1155:   PetscFunctionBegin;
    1160:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix A");
    1161:   PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix B");

    1163:   if (C) {
    1166:     PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix C");
    1167:   }

    1169:   PetscAssertPointer(D, 4);
    1170:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), D));
    1171:   /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
    1172:   PetscCall(MatProductCreate_Private(A, B, C, *D));
    1173:   PetscFunctionReturn(PETSC_SUCCESS);
    1174: }

    1176: /*
    1177:    These are safe basic implementations of ABC, RARt and PtAP
    1178:    that do not rely on mat->ops->matmatop function pointers.
    1179:    They only use the MatProduct API and are currently used by
    1180:    cuSPARSE and KOKKOS-KERNELS backends
    1181: */
    1182: typedef struct {
    1183:   Mat BC;
    1184:   Mat ABC;
    1185: } MatProductCtx_MatMatMatPrivate;

    1187: static PetscErrorCode MatProductCtxDestroy_MatMatMatPrivate(PetscCtxRt data)
    1188: {
    1189:   MatProductCtx_MatMatMatPrivate *mmdata = *(MatProductCtx_MatMatMatPrivate **)data;

    1191:   PetscFunctionBegin;
    1192:   PetscCall(MatDestroy(&mmdata->BC));
    1193:   PetscCall(MatDestroy(&mmdata->ABC));
    1194:   PetscCall(PetscFree(mmdata));
    1195:   PetscFunctionReturn(PETSC_SUCCESS);
    1196: }

    1198: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
    1199: {
    1200:   Mat_Product                    *product = mat->product;
    1201:   MatProductCtx_MatMatMatPrivate *mmabc;

    1203:   PetscFunctionBegin;
    1204:   MatCheckProduct(mat, 1);
    1205:   PetscCheck(mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data empty");
    1206:   mmabc = (MatProductCtx_MatMatMatPrivate *)mat->product->data;
    1207:   PetscCheck(mmabc->BC->ops->productnumeric, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing numeric stage");
    1208:   /* use function pointer directly to prevent logging */
    1209:   PetscCall((*mmabc->BC->ops->productnumeric)(mmabc->BC));
    1210:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
    1211:   mat->product             = mmabc->ABC->product;
    1212:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
    1213:   /* use function pointer directly to prevent logging */
    1214:   PetscUseTypeMethod(mat, productnumeric);
    1215:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
    1216:   mat->product             = product;
    1217:   PetscFunctionReturn(PETSC_SUCCESS);
    1218: }

    1220: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
    1221: {
    1222:   Mat_Product                    *product = mat->product;
    1223:   Mat                             A, B, C;
    1224:   MatProductType                  p1, p2;
    1225:   MatProductCtx_MatMatMatPrivate *mmabc;
    1226:   const char                     *prefix;

    1228:   PetscFunctionBegin;
    1229:   MatCheckProduct(mat, 1);
    1230:   PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data not empty");
    1231:   PetscCall(MatGetOptionsPrefix(mat, &prefix));
    1232:   PetscCall(PetscNew(&mmabc));
    1233:   product->data    = mmabc;
    1234:   product->destroy = MatProductCtxDestroy_MatMatMatPrivate;
    1235:   switch (product->type) {
    1236:   case MATPRODUCT_PtAP:
    1237:     p1 = MATPRODUCT_AB;
    1238:     p2 = MATPRODUCT_AtB;
    1239:     A  = product->B;
    1240:     B  = product->A;
    1241:     C  = product->B;
    1242:     if (A->cmap->bs > 0 && C->cmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, C->cmap->bs));
    1243:     break;
    1244:   case MATPRODUCT_RARt:
    1245:     p1 = MATPRODUCT_ABt;
    1246:     p2 = MATPRODUCT_AB;
    1247:     A  = product->B;
    1248:     B  = product->A;
    1249:     C  = product->B;
    1250:     if (A->rmap->bs > 0 && C->rmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, C->rmap->bs));
    1251:     break;
    1252:   case MATPRODUCT_ABC:
    1253:     p1 = MATPRODUCT_AB;
    1254:     p2 = MATPRODUCT_AB;
    1255:     A  = product->A;
    1256:     B  = product->B;
    1257:     C  = product->C;
    1258:     PetscCall(MatSetBlockSizesFromMats(mat, A, C));
    1259:     break;
    1260:   default:
    1261:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
    1262:   }
    1263:   PetscCall(MatProductCreate(B, C, NULL, &mmabc->BC));
    1264:   PetscCall(MatSetOptionsPrefix(mmabc->BC, prefix));
    1265:   PetscCall(MatAppendOptionsPrefix(mmabc->BC, "P1_"));
    1266:   PetscCall(MatProductSetType(mmabc->BC, p1));
    1267:   PetscCall(MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT));
    1268:   PetscCall(MatProductSetFill(mmabc->BC, product->fill));
    1269:   mmabc->BC->product->api_user = product->api_user;
    1270:   PetscCall(MatProductSetFromOptions(mmabc->BC));
    1271:   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);
    1272:   /* use function pointer directly to prevent logging */
    1273:   PetscCall((*mmabc->BC->ops->productsymbolic)(mmabc->BC));

    1275:   PetscCall(MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC));
    1276:   PetscCall(MatSetOptionsPrefix(mmabc->ABC, prefix));
    1277:   PetscCall(MatAppendOptionsPrefix(mmabc->ABC, "P2_"));
    1278:   PetscCall(MatProductSetType(mmabc->ABC, p2));
    1279:   PetscCall(MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT));
    1280:   PetscCall(MatProductSetFill(mmabc->ABC, product->fill));
    1281:   mmabc->ABC->product->api_user = product->api_user;
    1282:   PetscCall(MatProductSetFromOptions(mmabc->ABC));
    1283:   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);
    1284:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
    1285:   mat->product              = mmabc->ABC->product;
    1286:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
    1287:   /* use function pointer directly to prevent logging */
    1288:   PetscUseTypeMethod(mat, productsymbolic);
    1289:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
    1290:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
    1291:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
    1292:   mat->product                    = product;
    1293:   PetscFunctionReturn(PETSC_SUCCESS);
    1294: }

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

    1299:   Not Collective

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

    1304:   Output Parameter:
    1305: . mtype - the `MatProductType`

    1307:   Level: intermediate

    1309: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
    1310: @*/
    1311: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
    1312: {
    1313:   PetscFunctionBegin;
    1315:   PetscAssertPointer(mtype, 2);
    1316:   *mtype = MATPRODUCT_UNSPECIFIED;
    1317:   if (mat->product) *mtype = mat->product->type;
    1318:   PetscFunctionReturn(PETSC_SUCCESS);
    1319: }

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

    1324:   Not Collective

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

    1329:   Output Parameters:
    1330: + A - the first matrix
    1331: . B - the second matrix
    1332: - C - the third matrix (may be `NULL` for some `MatProductType`)

    1334:   Level: intermediate

    1336: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
    1337: @*/
    1338: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
    1339: {
    1340:   PetscFunctionBegin;
    1342:   if (A) *A = mat->product ? mat->product->A : NULL;
    1343:   if (B) *B = mat->product ? mat->product->B : NULL;
    1344:   if (C) *C = mat->product ? mat->product->C : NULL;
    1345:   PetscFunctionReturn(PETSC_SUCCESS);
    1346: }