Actual source code: matproduct.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  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_Basic(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_Basic(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_Basic;
 93:   return(0);
 94: }

 96: static PetscErrorCode MatProductNumeric_RARt_Basic(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_Basic(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_Basic;
139:   return(0);
140: }

142: static PetscErrorCode MatProductNumeric_ABC_Basic(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_Basic(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_Basic;
185:   return(0);
186: }

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

194:   switch (product->type) {
195:   case MATPRODUCT_PtAP:
196:     MatProductSymbolic_PtAP_Basic(mat);
197:     break;
198:   case MATPRODUCT_RARt:
199:     MatProductSymbolic_RARt_Basic(mat);
200:     break;
201:   case MATPRODUCT_ABC:
202:     MatProductSymbolic_ABC_Basic(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 (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
398:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
399:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
400:   else bname = Bnames[0];

402:   /* Check matrices sizes */
403:   Am = A->rmap->N;
404:   An = A->cmap->N;
405:   Bm = B->rmap->N;
406:   Bn = B->cmap->N;
407:   Cm = C ? C->rmap->N : 0;
408:   Cn = C ? C->cmap->N : 0;
409:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
410:   if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
411:   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);
412:   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);

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

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

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

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

474:       PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
475:       if (isdense) {

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

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

499:    Logically Collective on Mat

501:    Input Parameter:
502: .  mat - the matrix

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

507:    Level: intermediate

509: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
510: @*/
511: PetscErrorCode MatProductSetFromOptions(Mat mat)
512: {

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

528: /*@C
529:    MatProductView - View a MatProduct

531:    Logically Collective on Mat

533:    Input Parameter:
534: .  mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()

536:    Level: intermediate

538: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
539: @*/
540: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
541: {

546:   if (!mat->product) return(0);
547:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);}
550:   if (mat->product->view) {
551:     (*mat->product->view)(mat,viewer);
552:   }
553:   return(0);
554: }

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

566:   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);
567:   (*mat->ops->matmultnumeric)(A,B,mat);
568:   return(0);
569: }

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

578:   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);
579:   (*mat->ops->transposematmultnumeric)(A,B,mat);
580:   return(0);
581: }

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

590:   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);
591:   (*mat->ops->mattransposemultnumeric)(A,B,mat);
592:   return(0);
593: }

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

602:   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);
603:   (*mat->ops->ptapnumeric)(A,B,mat);
604:   return(0);
605: }

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

614:   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);
615:   (*mat->ops->rartnumeric)(A,B,mat);
616:   return(0);
617: }

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

626:   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);
627:   (*mat->ops->matmatmultnumeric)(A,B,C,mat);
628:   return(0);
629: }

631: /* ----------------------------------------------- */

633: /*@
634:    MatProductNumeric - Implement a matrix product with numerical values.

636:    Collective on Mat

638:    Input/Output Parameter:
639: .  mat - the matrix holding the product

641:    Level: intermediate

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

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

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

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

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

705: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
706: {
708:   Mat_Product    *product = mat->product;
709:   Mat            A=product->A,B=product->B;

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

718: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
719: {
721:   Mat_Product    *product = mat->product;
722:   Mat            A=product->A,B=product->B;

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

731: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
732: {
734:   Mat_Product    *product = mat->product;
735:   Mat            A=product->A,B=product->B,C=product->C;

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

744: /* ----------------------------------------------- */

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

749:    Collective on Mat

751:    Input/Output Parameter:
752: .  mat - the matrix to hold a product

754:    Level: intermediate

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

758: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
759: @*/
760: PetscErrorCode MatProductSymbolic(Mat mat)
761: {
763:   PetscLogEvent  eventtype=-1;

767:   MatCheckProduct(mat,1);
768:   if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot run symbolic phase. Product data not empty");

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

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

804: /*@
805:    MatProductSetFill - Set an expected fill of the matrix product.

807:    Collective on Mat

809:    Input Parameters:
810: +  mat - the matrix product
811: -  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.

813:    Level: intermediate

815: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
816: @*/
817: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
818: {
821:   MatCheckProduct(mat,1);
822:   if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
823:   else mat->product->fill = fill;
824:   return(0);
825: }

827: /*@
828:    MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.

830:    Collective on Mat

832:    Input Parameters:
833: +  mat - the matrix product
834: -  alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.

836:    Level: intermediate

838: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
839: @*/
840: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
841: {

846:   MatCheckProduct(mat,1);
847:   PetscFree(mat->product->alg);
848:   PetscStrallocpy(alg,&mat->product->alg);
849:   return(0);
850: }

852: /*@
853:    MatProductSetType - Sets a particular matrix product type

855:    Collective on Mat

857:    Input Parameters:
858: +  mat - the matrix
859: -  productype   - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.

861:    Level: intermediate

863: .seealso: MatProductCreate(), MatProductType
864: @*/
865: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
866: {

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

886: /*@
887:    MatProductClear - Clears matrix product internal structure.

889:    Collective on Mat

891:    Input Parameters:
892: .  mat - the product matrix

894:    Level: intermediate

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

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

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

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

941:   MatProductSetAlgorithm(D,MATPRODUCTALGORITHM_DEFAULT);
942:   MatProductSetFill(D,PETSC_DEFAULT);

944:   PetscObjectReference((PetscObject)A);
945:   PetscObjectReference((PetscObject)B);
946:   PetscObjectReference((PetscObject)C);
947:   return(0);
948: }

950: /*@
951:    MatProductCreateWithMat - Setup a given matrix as a matrix product.

953:    Collective on Mat

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

961:    Output Parameters:
962: .  D - the product matrix

964:    Notes:
965:      Any product data attached to D will be cleared

967:    Level: intermediate

969: .seealso: MatProductCreate(), MatProductClear()
970: @*/
971: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
972: {

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

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

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

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

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

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

1011:    Collective on Mat

1013:    Input Parameters:
1014: +  A - the first matrix
1015: .  B - the second matrix
1016: -  C - the third matrix (optional)

1018:    Output Parameters:
1019: .  D - the product matrix

1021:    Level: intermediate

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

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

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

1042:   if (C) {
1045:     MatCheckPreallocated(C,3);
1046:     if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1047:     if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1048:   }

1051:   MatCreate(PetscObjectComm((PetscObject)A),D);
1052:   MatProductCreate_Private(A,B,C,*D);
1053:   return(0);
1054: }