Actual source code: matproduct.c

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

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

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

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

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

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

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

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

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

     54:   /* AP = A*P */
     55:   MatProductNumeric(AP);
     56:   /* C = P^T*AP */
     57:   if (!C->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
     58:   (*C->ops->transposematmultnumeric)(P,AP,C);
     59:   return(0);
     60: }

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

     70:   PetscInfo2((PetscObject)C,"for A %s, P %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
     71:   /* AP = A*P */
     72:   MatProductCreate(A,P,NULL,&AP);
     73:   MatProductSetType(AP,MATPRODUCT_AB);
     74:   MatProductSetAlgorithm(AP,MATPRODUCTALGORITHM_DEFAULT);
     75:   MatProductSetFill(AP,fill);
     76:   MatProductSetFromOptions(AP);
     77:   MatProductSymbolic(AP);

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

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

     92:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
     93:   return(0);
     94: }

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

    103:   /* RA = R*A */
    104:   MatProductNumeric(RA);
    105:   /* C = RA*R^T */
    106:   if (!C->ops->mattransposemultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
    107:   (*C->ops->mattransposemultnumeric)(RA,R,C);
    108:   return(0);
    109: }

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

    119:   PetscInfo2((PetscObject)C,"for A %s, R %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
    120:   /* RA = R*A */
    121:   MatProductCreate(R,A,NULL,&RA);
    122:   MatProductSetType(RA,MATPRODUCT_AB);
    123:   MatProductSetAlgorithm(RA,MATPRODUCTALGORITHM_DEFAULT);
    124:   MatProductSetFill(RA,fill);
    125:   MatProductSetFromOptions(RA);
    126:   MatProductSymbolic(RA);

    128:   /* C = RA*R^T */
    129:   MatProductSetType(C,MATPRODUCT_ABt);
    130:   MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);
    131:   product->A = RA;
    132:   MatProductSetFromOptions(C);
    133:   MatProductSymbolic(C);

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

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

    149:   /* Numeric BC = B*C */
    150:   MatProductNumeric(BC);
    151:   /* Numeric mat = A*BC */
    152:   if (!mat->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
    153:   (*mat->ops->matmultnumeric)(A,BC,mat);
    154:   return(0);
    155: }

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

    165:   PetscInfo3((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);
    166:   /* Symbolic BC = B*C */
    167:   MatProductCreate(B,C,NULL,&BC);
    168:   MatProductSetType(BC,MATPRODUCT_AB);
    169:   MatProductSetAlgorithm(BC,MATPRODUCTALGORITHM_DEFAULT);
    170:   MatProductSetFill(BC,fill);
    171:   MatProductSetFromOptions(BC);
    172:   MatProductSymbolic(BC);

    174:   /* Symbolic mat = A*BC */
    175:   MatProductSetType(mat,MATPRODUCT_AB);
    176:   MatProductSetAlgorithm(mat,MATPRODUCTALGORITHM_DEFAULT);
    177:   product->B     = BC;
    178:   product->Dwork = BC;
    179:   MatProductSetFromOptions(mat);
    180:   MatProductSymbolic(mat);

    182:   /* resume user's original input matrix setting for B */
    183:   product->B = B;
    184:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
    185:   return(0);
    186: }

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

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

    209: /* ----------------------------------------------- */
    210: /*@C
    211:    MatProductReplaceMats - Replace input matrices for a matrix product.

    213:    Collective on Mat

    215:    Input Parameters:
    216: +  A - the matrix or NULL if not being replaced
    217: .  B - the matrix or NULL if not being replaced
    218: .  C - the matrix or NULL if not being replaced
    219: -  D - the matrix product

    221:    Level: intermediate

    223:    Notes:
    224:      To reuse the symbolic phase, input matrices must have exactly the same data structure as the replaced one.
    225:      If the type of any of the input matrices is different than what previously used, the product is cleared and MatProductSetFromOptions()/MatProductSymbolic() are invoked again.

    227: .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
    228: @*/
    229: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
    230: {
    232:   Mat_Product    *product;
    233:   PetscBool      flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;

    237:   MatCheckProduct(D,4);
    238:   product = D->product;
    239:   if (A) {
    241:     PetscObjectReference((PetscObject)A);
    242:     PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);
    243:     MatDestroy(&product->A);
    244:     product->A = A;
    245:   }
    246:   if (B) {
    248:     PetscObjectReference((PetscObject)B);
    249:     PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);
    250:     MatDestroy(&product->B);
    251:     product->B = B;
    252:   }
    253:   if (C) {
    255:     PetscObjectReference((PetscObject)C);
    256:     PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);
    257:     MatDestroy(&product->C);
    258:     product->C = C;
    259:   }
    260:   /* Any of the replaced mats is of a different type, reset */
    261:   if (!flgA || !flgB || !flgC) {
    262:     if (D->product->destroy) {
    263:       (*D->product->destroy)(D->product->data);
    264:     }
    265:     D->product->destroy = NULL;
    266:     D->product->data = NULL;
    267:     if (D->ops->productnumeric || D->ops->productsymbolic) {
    268:       MatProductSetFromOptions(D);
    269:       MatProductSymbolic(D);
    270:     }
    271:   }
    272:   return(0);
    273: }

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

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

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

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

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

    359:   switch (product->type) {
    360:   case MATPRODUCT_AB:
    361:     MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
    362:     break;
    363:   case MATPRODUCT_AtB:
    364:     MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
    365:     break;
    366:   default: SETERRQ3(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:   PetscObjectBaseTypeCompareAny((PetscObject)C,&isdense,MATSEQDENSE,MATMPIDENSE,"");
    369:   if (!isdense) {
    370:     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:   MatSetUp(C);
    376:   return(0);
    377: }

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

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

    405:   /* Check matrices sizes */
    406:   Am = A->rmap->N;
    407:   An = A->cmap->N;
    408:   Bm = B->rmap->N;
    409:   Bn = B->cmap->N;
    410:   Cm = C ? C->rmap->N : 0;
    411:   Cn = C ? C->cmap->N : 0;
    412:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
    413:   if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
    414:   if (An != Bm) SETERRQ7(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of A and %s are incompatible for MatProductType %s: A %Dx%D, %s %Dx%D",bname,MatProductTypes[product->type],A->rmap->N,A->cmap->N,bname,B->rmap->N,B->cmap->N);
    415:   if (Cm && Cm != Bn) SETERRQ5(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of B and C are incompatible for MatProductType %s: B %Dx%D, C %Dx%D",MatProductTypes[product->type],B->rmap->N,B->cmap->N,Cm,Cn);

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

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

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

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

    477:       PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
    478:       if (isdense) {

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

    499: /*@C
    500:    MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.

    502:    Logically Collective on Mat

    504:    Input Parameter:
    505: .  mat - the matrix

    507:    Options Database Keys:
    508: .    -mat_product_clear - Clear intermediate data structures after MatProductNumeric() has been called

    510:    Level: intermediate

    512: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
    513: @*/
    514: PetscErrorCode MatProductSetFromOptions(Mat mat)
    515: {

    520:   MatCheckProduct(mat,1);
    521:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot call MatProductSetFromOptions with already present data");
    522:   PetscObjectOptionsBegin((PetscObject)mat);
    523:   PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);
    524:   PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");
    525:   PetscOptionsEnd();
    526:   MatProductSetFromOptions_Private(mat);
    527:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after setup phase");
    528:   return(0);
    529: }

    531: /*@C
    532:    MatProductView - View a MatProduct

    534:    Logically Collective on Mat

    536:    Input Parameter:
    537: .  mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()

    539:    Level: intermediate

    541: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
    542: @*/
    543: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
    544: {

    549:   if (!mat->product) return(0);
    550:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);}
    553:   if (mat->product->view) {
    554:     (*mat->product->view)(mat,viewer);
    555:   }
    556:   return(0);
    557: }

    559: /* ----------------------------------------------- */
    560: /* these are basic implementations relying on the old function pointers
    561:  * they are dangerous and should be removed in the future */
    562: PetscErrorCode MatProductNumeric_AB(Mat mat)
    563: {
    565:   Mat_Product    *product = mat->product;
    566:   Mat            A=product->A,B=product->B;

    569:   if (!mat->ops->matmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    570:   (*mat->ops->matmultnumeric)(A,B,mat);
    571:   return(0);
    572: }

    574: PetscErrorCode MatProductNumeric_AtB(Mat mat)
    575: {
    577:   Mat_Product    *product = mat->product;
    578:   Mat            A=product->A,B=product->B;

    581:   if (!mat->ops->transposematmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    582:   (*mat->ops->transposematmultnumeric)(A,B,mat);
    583:   return(0);
    584: }

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

    593:   if (!mat->ops->mattransposemultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    594:   (*mat->ops->mattransposemultnumeric)(A,B,mat);
    595:   return(0);
    596: }

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

    605:   if (!mat->ops->ptapnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    606:   (*mat->ops->ptapnumeric)(A,B,mat);
    607:   return(0);
    608: }

    610: PetscErrorCode MatProductNumeric_RARt(Mat mat)
    611: {
    613:   Mat_Product    *product = mat->product;
    614:   Mat            A=product->A,B=product->B;

    617:   if (!mat->ops->rartnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    618:   (*mat->ops->rartnumeric)(A,B,mat);
    619:   return(0);
    620: }

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

    629:   if (!mat->ops->matmatmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
    630:   (*mat->ops->matmatmultnumeric)(A,B,C,mat);
    631:   return(0);
    632: }

    634: /* ----------------------------------------------- */

    636: /*@
    637:    MatProductNumeric - Implement a matrix product with numerical values.

    639:    Collective on Mat

    641:    Input/Output Parameter:
    642: .  mat - the matrix holding the product

    644:    Level: intermediate

    646:    Notes: MatProductSymbolic() must have been called on mat before calling this function

    648: .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
    649: @*/
    650: PetscErrorCode MatProductNumeric(Mat mat)
    651: {
    653:   PetscLogEvent  eventtype=-1;

    657:   MatCheckProduct(mat,1);
    658:   /* log event */
    659:   switch (mat->product->type) {
    660:   case MATPRODUCT_AB:
    661:     eventtype = MAT_MatMultNumeric;
    662:     break;
    663:   case MATPRODUCT_AtB:
    664:     eventtype = MAT_TransposeMatMultNumeric;
    665:     break;
    666:   case MATPRODUCT_ABt:
    667:     eventtype = MAT_MatTransposeMultNumeric;
    668:     break;
    669:   case MATPRODUCT_PtAP:
    670:     eventtype = MAT_PtAPNumeric;
    671:     break;
    672:   case MATPRODUCT_RARt:
    673:     eventtype = MAT_RARtNumeric;
    674:     break;
    675:   case MATPRODUCT_ABC:
    676:     eventtype = MAT_MatMatMultNumeric;
    677:     break;
    678:   default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
    679:   }
    680:   if (mat->ops->productnumeric) {
    681:     PetscLogEventBegin(eventtype,mat,0,0,0);
    682:     (*mat->ops->productnumeric)(mat);
    683:     PetscLogEventEnd(eventtype,mat,0,0,0);
    684:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSymbolic() first");
    685:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after numeric phase");
    686:   if (mat->product->clear) {
    687:     MatProductClear(mat);
    688:   }
    689:   PetscObjectStateIncrease((PetscObject)mat);
    690:   return(0);
    691: }

    693: /* ----------------------------------------------- */
    694: /* these are basic implementations relying on the old function pointers
    695:  * they are dangerous and should be removed in the future */
    696: PetscErrorCode MatProductSymbolic_AB(Mat mat)
    697: {
    699:   Mat_Product    *product = mat->product;
    700:   Mat            A=product->A,B=product->B;

    703:   if (!mat->ops->matmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
    704:   (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);
    705:   mat->ops->productnumeric = MatProductNumeric_AB;
    706:   return(0);
    707: }

    709: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
    710: {
    712:   Mat_Product    *product = mat->product;
    713:   Mat            A=product->A,B=product->B;

    716:   if (!mat->ops->transposematmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
    717:   (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
    718:   mat->ops->productnumeric = MatProductNumeric_AtB;
    719:   return(0);
    720: }

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

    729:   if (!mat->ops->mattransposemultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
    730:   (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
    731:   mat->ops->productnumeric = MatProductNumeric_ABt;
    732:   return(0);
    733: }

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

    742:   if (!mat->ops->matmatmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
    743:   (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
    744:   mat->ops->productnumeric = MatProductNumeric_ABC;
    745:   return(0);
    746: }

    748: /* ----------------------------------------------- */

    750: /*@
    751:    MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.

    753:    Collective on Mat

    755:    Input/Output Parameter:
    756: .  mat - the matrix to hold a product

    758:    Level: intermediate

    760:    Notes: MatProductSetFromOptions() must have been called on mat before calling this function

    762: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
    763: @*/
    764: PetscErrorCode MatProductSymbolic(Mat mat)
    765: {
    767:   PetscLogEvent  eventtype=-1;

    771:   MatCheckProduct(mat,1);
    772:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot run symbolic phase. Product data not empty");
    773:   /* log event */
    774:   switch (mat->product->type) {
    775:   case MATPRODUCT_AB:
    776:     eventtype = MAT_MatMultSymbolic;
    777:     break;
    778:   case MATPRODUCT_AtB:
    779:     eventtype = MAT_TransposeMatMultSymbolic;
    780:     break;
    781:   case MATPRODUCT_ABt:
    782:     eventtype = MAT_MatTransposeMultSymbolic;
    783:     break;
    784:   case MATPRODUCT_PtAP:
    785:     eventtype = MAT_PtAPSymbolic;
    786:     break;
    787:   case MATPRODUCT_RARt:
    788:     eventtype = MAT_RARtSymbolic;
    789:     break;
    790:   case MATPRODUCT_ABC:
    791:     eventtype = MAT_MatMatMultSymbolic;
    792:     break;
    793:   default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
    794:   }

    796:   mat->ops->productnumeric = NULL;
    797:   if (mat->ops->productsymbolic) {
    798:     PetscLogEventBegin(eventtype,mat,0,0,0);
    799:     (*mat->ops->productsymbolic)(mat);
    800:     PetscLogEventEnd(eventtype,mat,0,0,0);
    801:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSetFromOptions() first");
    802:   if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after symbolic phase");
    803:   if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Symbolic phase did not specify the numeric phase");
    804:   return(0);
    805: }

    807: /*@
    808:    MatProductSetFill - Set an expected fill of the matrix product.

    810:    Collective on Mat

    812:    Input Parameters:
    813: +  mat - the matrix product
    814: -  fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use PETSC_DEFAULT if you do not have a good estimate. If the product is a dense matrix, this is irrelevent.

    816:    Level: intermediate

    818: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
    819: @*/
    820: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
    821: {
    824:   MatCheckProduct(mat,1);
    825:   if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
    826:   else mat->product->fill = fill;
    827:   return(0);
    828: }

    830: /*@
    831:    MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.

    833:    Collective on Mat

    835:    Input Parameters:
    836: +  mat - the matrix product
    837: -  alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.

    839:    Level: intermediate

    841: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
    842: @*/
    843: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
    844: {

    849:   MatCheckProduct(mat,1);
    850:   PetscFree(mat->product->alg);
    851:   PetscStrallocpy(alg,&mat->product->alg);
    852:   return(0);
    853: }

    855: /*@
    856:    MatProductSetType - Sets a particular matrix product type

    858:    Collective on Mat

    860:    Input Parameters:
    861: +  mat - the matrix
    862: -  productype   - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.

    864:    Level: intermediate

    866: .seealso: MatProductCreate(), MatProductType
    867: @*/
    868: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
    869: {

    874:   MatCheckProduct(mat,1);
    876:   if (productype != mat->product->type) {
    877:     if (mat->product->destroy) {
    878:       (*mat->product->destroy)(mat->product->data);
    879:     }
    880:     mat->product->destroy = NULL;
    881:     mat->product->data = NULL;
    882:     mat->ops->productsymbolic = NULL;
    883:     mat->ops->productnumeric  = NULL;
    884:   }
    885:   mat->product->type = productype;
    886:   return(0);
    887: }

    889: /*@
    890:    MatProductClear - Clears matrix product internal structure.

    892:    Collective on Mat

    894:    Input Parameters:
    895: .  mat - the product matrix

    897:    Level: intermediate

    899:    Notes: this function should be called to remove any intermediate data used by the product
    900:           After having called this function, MatProduct operations can no longer be used on mat
    901: @*/
    902: PetscErrorCode MatProductClear(Mat mat)
    903: {
    905:   Mat_Product    *product = mat->product;

    909:   if (product) {
    910:     MatDestroy(&product->A);
    911:     MatDestroy(&product->B);
    912:     MatDestroy(&product->C);
    913:     PetscFree(product->alg);
    914:     MatDestroy(&product->Dwork);
    915:     if (product->destroy) {
    916:       (*product->destroy)(product->data);
    917:     }
    918:   }
    919:   PetscFree(mat->product);
    920:   mat->ops->productsymbolic = NULL;
    921:   mat->ops->productnumeric = NULL;
    922:   return(0);
    923: }

    925: /* Create a supporting struct and attach it to the matrix product */
    926: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
    927: {
    929:   Mat_Product    *product=NULL;

    933:   if (D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product already present");
    934:   PetscNewLog(D,&product);
    935:   product->A        = A;
    936:   product->B        = B;
    937:   product->C        = C;
    938:   product->type     = MATPRODUCT_UNSPECIFIED;
    939:   product->Dwork    = NULL;
    940:   product->api_user = PETSC_FALSE;
    941:   product->clear    = PETSC_FALSE;
    942:   D->product        = product;

    944:   MatProductSetAlgorithm(D,MATPRODUCTALGORITHM_DEFAULT);
    945:   MatProductSetFill(D,PETSC_DEFAULT);

    947:   PetscObjectReference((PetscObject)A);
    948:   PetscObjectReference((PetscObject)B);
    949:   PetscObjectReference((PetscObject)C);
    950:   return(0);
    951: }

    953: /*@
    954:    MatProductCreateWithMat - Setup a given matrix as a matrix product.

    956:    Collective on Mat

    958:    Input Parameters:
    959: +  A - the first matrix
    960: .  B - the second matrix
    961: .  C - the third matrix (optional)
    962: -  D - the matrix which will be used as a product

    964:    Output Parameters:
    965: .  D - the product matrix

    967:    Notes:
    968:      Any product data attached to D will be cleared

    970:    Level: intermediate

    972: .seealso: MatProductCreate(), MatProductClear()
    973: @*/
    974: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
    975: {

    981:   MatCheckPreallocated(A,1);
    982:   if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    983:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    987:   MatCheckPreallocated(B,2);
    988:   if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    989:   if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    991:   if (C) {
    994:     MatCheckPreallocated(C,3);
    995:     if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    996:     if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    997:   }

    1001:   MatCheckPreallocated(D,4);
    1002:   if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1003:   if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    1005:   /* Create a supporting struct and attach it to D */
    1006:   MatProductClear(D);
    1007:   MatProductCreate_Private(A,B,C,D);
    1008:   return(0);
    1009: }

    1011: /*@
    1012:    MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.

    1014:    Collective on Mat

    1016:    Input Parameters:
    1017: +  A - the first matrix
    1018: .  B - the second matrix
    1019: -  C - the third matrix (optional)

    1021:    Output Parameters:
    1022: .  D - the product matrix

    1024:    Level: intermediate

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

    1037:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix A");
    1038:   if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix B");

    1040:   if (C) {
    1043:     if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix C");
    1044:   }

    1047:   MatCreate(PetscObjectComm((PetscObject)A),D);
    1048:   MatProductCreate_Private(A,B,C,*D);
    1049:   return(0);
    1050: }

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

    1063: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
    1064: {
    1065:   PetscErrorCode   ierr;
    1066:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

    1069:   MatDestroy(&mmdata->BC);
    1070:   MatDestroy(&mmdata->ABC);
    1071:   PetscFree(data);
    1072:   return(0);
    1073: }

    1075: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
    1076: {
    1077:   PetscErrorCode   ierr;
    1078:   Mat_Product      *product = mat->product;
    1079:   MatMatMatPrivate *mmabc;

    1082:   MatCheckProduct(mat,1);
    1083:   if (!mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data empty");
    1084:   mmabc = (MatMatMatPrivate *)mat->product->data;
    1085:   if (!mmabc->BC->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
    1086:   /* use function pointer directly to prevent logging */
    1087:   (*mmabc->BC->ops->productnumeric)(mmabc->BC);
    1088:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
    1089:   mat->product = mmabc->ABC->product;
    1090:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
    1091:   if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
    1092:   /* use function pointer directly to prevent logging */
    1093:   (*mat->ops->productnumeric)(mat);
    1094:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
    1095:   mat->product = product;
    1096:   return(0);
    1097: }

    1099: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
    1100: {
    1101:   PetscErrorCode   ierr;
    1102:   Mat_Product      *product = mat->product;
    1103:   Mat              A, B ,C;
    1104:   MatProductType   p1,p2;
    1105:   MatMatMatPrivate *mmabc;
    1106:   const char       *prefix;

    1109:   MatCheckProduct(mat,1);
    1110:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data not empty");
    1111:   MatGetOptionsPrefix(mat,&prefix);
    1112:   PetscNew(&mmabc);
    1113:   product->data    = mmabc;
    1114:   product->destroy = MatDestroy_MatMatMatPrivate;
    1115:   switch (product->type) {
    1116:   case MATPRODUCT_PtAP:
    1117:     p1 = MATPRODUCT_AB;
    1118:     p2 = MATPRODUCT_AtB;
    1119:     A = product->B;
    1120:     B = product->A;
    1121:     C = product->B;
    1122:     break;
    1123:   case MATPRODUCT_RARt:
    1124:     p1 = MATPRODUCT_ABt;
    1125:     p2 = MATPRODUCT_AB;
    1126:     A = product->B;
    1127:     B = product->A;
    1128:     C = product->B;
    1129:     break;
    1130:   case MATPRODUCT_ABC:
    1131:     p1 = MATPRODUCT_AB;
    1132:     p2 = MATPRODUCT_AB;
    1133:     A = product->A;
    1134:     B = product->B;
    1135:     C = product->C;
    1136:     break;
    1137:   default:
    1138:     SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Not for ProductType %s",MatProductTypes[product->type]);
    1139:   }
    1140:   MatProductCreate(B,C,NULL,&mmabc->BC);
    1141:   MatSetOptionsPrefix(mmabc->BC,prefix);
    1142:   MatAppendOptionsPrefix(mmabc->BC,"P1_");
    1143:   MatProductSetType(mmabc->BC,p1);
    1144:   MatProductSetAlgorithm(mmabc->BC,MATPRODUCTALGORITHM_DEFAULT);
    1145:   MatProductSetFill(mmabc->BC,product->fill);
    1146:   mmabc->BC->product->api_user = product->api_user;
    1147:   MatProductSetFromOptions(mmabc->BC);
    1148:   if (!mmabc->BC->ops->productsymbolic) SETERRQ3(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);
    1149:   /* use function pointer directly to prevent logging */
    1150:   (*mmabc->BC->ops->productsymbolic)(mmabc->BC);

    1152:   MatProductCreate(A,mmabc->BC,NULL,&mmabc->ABC);
    1153:   MatSetOptionsPrefix(mmabc->ABC,prefix);
    1154:   MatAppendOptionsPrefix(mmabc->ABC,"P2_");
    1155:   MatProductSetType(mmabc->ABC,p2);
    1156:   MatProductSetAlgorithm(mmabc->ABC,MATPRODUCTALGORITHM_DEFAULT);
    1157:   MatProductSetFill(mmabc->ABC,product->fill);
    1158:   mmabc->ABC->product->api_user = product->api_user;
    1159:   MatProductSetFromOptions(mmabc->ABC);
    1160:   if (!mmabc->ABC->ops->productsymbolic) SETERRQ3(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);
    1161:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
    1162:   mat->product = mmabc->ABC->product;
    1163:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
    1164:   /* use function pointer directly to prevent logging */
    1165:   (*mat->ops->productsymbolic)(mat);
    1166:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
    1167:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
    1168:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
    1169:   mat->product                    = product;
    1170:   return(0);
    1171: }