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 MATTRANSPOSEMAT)
 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: {
 49:   Mat_Product    *product = C->product;
 50:   Mat            P = product->B,AP = product->Dwork;

 52:   /* AP = A*P */
 53:   MatProductNumeric(AP);
 54:   /* C = P^T*AP */
 56:   (*C->ops->transposematmultnumeric)(P,AP,C);
 57:   return 0;
 58: }

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

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

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

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

 88:   C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
 89:   return 0;
 90: }

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

 97:   /* RA = R*A */
 98:   MatProductNumeric(RA);
 99:   /* C = RA*R^T */
101:   (*C->ops->mattransposemultnumeric)(RA,R,C);
102:   return 0;
103: }

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

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

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

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

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

139:   /* Numeric BC = B*C */
140:   MatProductNumeric(BC);
141:   /* Numeric mat = A*BC */
143:   (*mat->ops->matmultnumeric)(A,BC,mat);
144:   return 0;
145: }

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

153:   PetscInfo((PetscObject)mat,"for A %s, B %s, C %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name,((PetscObject)product->C)->type_name);
154:   /* Symbolic BC = B*C */
155:   MatProductCreate(B,C,NULL,&BC);
156:   MatProductSetType(BC,MATPRODUCT_AB);
157:   MatProductSetAlgorithm(BC,MATPRODUCTALGORITHMDEFAULT);
158:   MatProductSetFill(BC,fill);
159:   MatProductSetFromOptions(BC);
160:   MatProductSymbolic(BC);

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

170:   /* resume user's original input matrix setting for B */
171:   product->B = B;
172:   mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
173:   return 0;
174: }

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

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

195: /* ----------------------------------------------- */
196: /*@C
197:    MatProductReplaceMats - Replace input matrices for a matrix product.

199:    Collective on Mat

201:    Input Parameters:
202: +  A - the matrix or NULL if not being replaced
203: .  B - the matrix or NULL if not being replaced
204: .  C - the matrix or NULL if not being replaced
205: -  D - the matrix product

207:    Level: intermediate

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

214: .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
215: @*/
216: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
217: {
218:   Mat_Product    *product;
219:   PetscBool      flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;

222:   MatCheckProduct(D,4);
223:   product = D->product;
224:   if (A) {
226:     PetscObjectReference((PetscObject)A);
227:     PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);
228:     if (product->symbolic_used_the_fact_A_is_symmetric && !A->symmetric) { /* symbolic was built around a symmetric A, but the new A is not anymore */
229:       flgA = PETSC_FALSE;
230:       product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
231:     }
232:     MatDestroy(&product->A);
233:     product->A = A;
234:   }
235:   if (B) {
237:     PetscObjectReference((PetscObject)B);
238:     PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);
239:     if (product->symbolic_used_the_fact_B_is_symmetric && !B->symmetric) {
240:       flgB = PETSC_FALSE;
241:       product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
242:     }
243:     MatDestroy(&product->B);
244:     product->B = B;
245:   }
246:   if (C) {
248:     PetscObjectReference((PetscObject)C);
249:     PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);
250:     if (product->symbolic_used_the_fact_C_is_symmetric && !C->symmetric) {
251:       flgC = PETSC_FALSE;
252:       product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
253:     }
254:     MatDestroy(&product->C);
255:     product->C = C;
256:   }
257:   /* Any of the replaced mats is of a different type, reset */
258:   if (!flgA || !flgB || !flgC) {
259:     if (D->product->destroy) {
260:       (*D->product->destroy)(D->product->data);
261:     }
262:     D->product->destroy = NULL;
263:     D->product->data = NULL;
264:     if (D->ops->productnumeric || D->ops->productsymbolic) {
265:       MatProductSetFromOptions(D);
266:       MatProductSymbolic(D);
267:     }
268:   }
269:   return 0;
270: }

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

281:   switch (product->type) {
282:   case MATPRODUCT_AB:
283:     t = PETSC_FALSE;
284:   case MATPRODUCT_AtB:
285:     break;
286:   default: 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);
287:   }
288:   if (PetscDefined(HAVE_CUDA)) {
289:     VecType vtype;

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

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

346: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
347: {
348:   Mat_Product    *product = C->product;
349:   Mat            A = product->A, B = product->B;
350:   PetscBool      isdense;

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

372: /* a single driver to query the dispatching */
373: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
374: {
375:   Mat_Product       *product = mat->product;
376:   PetscInt          Am,An,Bm,Bn,Cm,Cn;
377:   Mat               A = product->A,B = product->B,C = product->C;
378:   const char* const Bnames[] = { "B", "R", "P" };
379:   const char*       bname;
380:   PetscErrorCode    (*fA)(Mat);
381:   PetscErrorCode    (*fB)(Mat);
382:   PetscErrorCode    (*fC)(Mat);
383:   PetscErrorCode    (*f)(Mat)=NULL;

385:   mat->ops->productsymbolic = NULL;
386:   mat->ops->productnumeric = NULL;
387:   if (product->type == MATPRODUCT_UNSPECIFIED) return 0;
391:   if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
392:   if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
393:   else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
394:   else bname = Bnames[0];

396:   /* Check matrices sizes */
397:   Am = A->rmap->N;
398:   An = A->cmap->N;
399:   Bm = B->rmap->N;
400:   Bn = B->cmap->N;
401:   Cm = C ? C->rmap->N : 0;
402:   Cn = C ? C->cmap->N : 0;
403:   if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
404:   if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }

408:   fA = A->ops->productsetfromoptions;
409:   fB = B->ops->productsetfromoptions;
410:   fC = C ? C->ops->productsetfromoptions : fA;
411:   if (C) {
412:     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);
413:   } else {
414:     PetscInfo(mat,"MatProductType %s for A %s, %s %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name);
415:   }
416:   if (fA == fB && fA == fC && fA) {
417:     PetscInfo(mat,"  matching op\n");
418:     (*fA)(mat);
419:   }
420:   /* We may have found f but it did not succeed */
421:   if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
422:     char  mtypes[256];
423:     PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
424:     PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
425:     PetscStrlcat(mtypes,"_",sizeof(mtypes));
426:     PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
427:     if (C) {
428:       PetscStrlcat(mtypes,"_",sizeof(mtypes));
429:       PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));
430:     }
431:     PetscStrlcat(mtypes,"_C",sizeof(mtypes));

433:     PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
434:     PetscInfo(mat,"  querying %s from A? %p\n",mtypes,f);
435:     if (!f) {
436:       PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
437:       PetscInfo(mat,"  querying %s from %s? %p\n",mtypes,bname,f);
438:     }
439:     if (!f && C) {
440:       PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
441:       PetscInfo(mat,"  querying %s from C? %p\n",mtypes,f);
442:     }
443:     if (f) (*f)(mat);

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

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

469:       PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
470:       if (isdense) {

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

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

494:    Logically Collective on Mat

496:    Input Parameter:
497: .  mat - the matrix

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

502:    Level: intermediate

504: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
505: @*/
506: PetscErrorCode MatProductSetFromOptions(Mat mat)
507: {

511:   MatCheckProduct(mat,1);
513:   PetscObjectOptionsBegin((PetscObject)mat);
514:   PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);
515:   PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");
516:   PetscOptionsEnd();
517:   MatProductSetFromOptions_Private(mat);
519:   return 0;
520: }

522: /*@C
523:    MatProductView - View a MatProduct

525:    Logically Collective on Mat

527:    Input Parameter:
528: .  mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()

530:    Level: intermediate

532: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
533: @*/
534: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
535: {
537:   if (!mat->product) return 0;
538:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);
541:   if (mat->product->view) {
542:     (*mat->product->view)(mat,viewer);
543:   }
544:   return 0;
545: }

547: /* ----------------------------------------------- */
548: /* these are basic implementations relying on the old function pointers
549:  * they are dangerous and should be removed in the future */
550: PetscErrorCode MatProductNumeric_AB(Mat mat)
551: {
552:   Mat_Product    *product = mat->product;
553:   Mat            A=product->A,B=product->B;

556:   (*mat->ops->matmultnumeric)(A,B,mat);
557:   return 0;
558: }

560: PetscErrorCode MatProductNumeric_AtB(Mat mat)
561: {
562:   Mat_Product    *product = mat->product;
563:   Mat            A=product->A,B=product->B;

566:   (*mat->ops->transposematmultnumeric)(A,B,mat);
567:   return 0;
568: }

570: PetscErrorCode MatProductNumeric_ABt(Mat mat)
571: {
572:   Mat_Product    *product = mat->product;
573:   Mat            A=product->A,B=product->B;

576:   (*mat->ops->mattransposemultnumeric)(A,B,mat);
577:   return 0;
578: }

580: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
581: {
582:   Mat_Product    *product = mat->product;
583:   Mat            A=product->A,B=product->B;

586:   (*mat->ops->ptapnumeric)(A,B,mat);
587:   return 0;
588: }

590: PetscErrorCode MatProductNumeric_RARt(Mat mat)
591: {
592:   Mat_Product    *product = mat->product;
593:   Mat            A=product->A,B=product->B;

596:   (*mat->ops->rartnumeric)(A,B,mat);
597:   return 0;
598: }

600: PetscErrorCode MatProductNumeric_ABC(Mat mat)
601: {
602:   Mat_Product    *product = mat->product;
603:   Mat            A=product->A,B=product->B,C=product->C;

606:   (*mat->ops->matmatmultnumeric)(A,B,C,mat);
607:   return 0;
608: }

610: /* ----------------------------------------------- */

612: /*@
613:    MatProductNumeric - Implement a matrix product with numerical values.

615:    Collective on Mat

617:    Input/Output Parameter:
618: .  mat - the matrix holding the product

620:    Level: intermediate

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

624: .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
625: @*/
626: PetscErrorCode MatProductNumeric(Mat mat)
627: {
628:   PetscLogEvent  eventtype = -1;
629:   PetscBool      missing = PETSC_FALSE;

632:   MatCheckProduct(mat,1);
633:   switch (mat->product->type) {
634:   case MATPRODUCT_AB:
635:     eventtype = MAT_MatMultNumeric;
636:     break;
637:   case MATPRODUCT_AtB:
638:     eventtype = MAT_TransposeMatMultNumeric;
639:     break;
640:   case MATPRODUCT_ABt:
641:     eventtype = MAT_MatTransposeMultNumeric;
642:     break;
643:   case MATPRODUCT_PtAP:
644:     eventtype = MAT_PtAPNumeric;
645:     break;
646:   case MATPRODUCT_RARt:
647:     eventtype = MAT_RARtNumeric;
648:     break;
649:   case MATPRODUCT_ABC:
650:     eventtype = MAT_MatMatMultNumeric;
651:     break;
652:   default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
653:   }

655:   if (mat->ops->productnumeric) {
656:     PetscLogEventBegin(eventtype,mat,0,0,0);
657:     (*mat->ops->productnumeric)(mat);
658:     PetscLogEventEnd(eventtype,mat,0,0,0);
659:   } else missing = PETSC_TRUE;

661:   if (missing || !mat->product) {
662:     char errstr[256];

664:     if (mat->product->type == MATPRODUCT_ABC) {
665:       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);
666:     } else {
667:       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);
668:     }
671:   }

673:   if (mat->product->clear) {
674:     MatProductClear(mat);
675:   }
676:   PetscObjectStateIncrease((PetscObject)mat);
677:   return 0;
678: }

680: /* ----------------------------------------------- */
681: /* these are basic implementations relying on the old function pointers
682:  * they are dangerous and should be removed in the future */
683: PetscErrorCode MatProductSymbolic_AB(Mat mat)
684: {
685:   Mat_Product    *product = mat->product;
686:   Mat            A=product->A,B=product->B;

689:   (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);
690:   mat->ops->productnumeric = MatProductNumeric_AB;
691:   return 0;
692: }

694: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
695: {
696:   Mat_Product    *product = mat->product;
697:   Mat            A=product->A,B=product->B;

700:   (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
701:   mat->ops->productnumeric = MatProductNumeric_AtB;
702:   return 0;
703: }

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

711:   (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
712:   mat->ops->productnumeric = MatProductNumeric_ABt;
713:   return 0;
714: }

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

722:   (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
723:   mat->ops->productnumeric = MatProductNumeric_ABC;
724:   return 0;
725: }

727: /* ----------------------------------------------- */

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

732:    Collective on Mat

734:    Input/Output Parameter:
735: .  mat - the matrix to hold a product

737:    Level: intermediate

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

741: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
742: @*/
743: PetscErrorCode MatProductSymbolic(Mat mat)
744: {
745:   PetscLogEvent  eventtype = -1;
746:   PetscBool      missing = PETSC_FALSE;

749:   MatCheckProduct(mat,1);
751:   switch (mat->product->type) {
752:   case MATPRODUCT_AB:
753:     eventtype = MAT_MatMultSymbolic;
754:     break;
755:   case MATPRODUCT_AtB:
756:     eventtype = MAT_TransposeMatMultSymbolic;
757:     break;
758:   case MATPRODUCT_ABt:
759:     eventtype = MAT_MatTransposeMultSymbolic;
760:     break;
761:   case MATPRODUCT_PtAP:
762:     eventtype = MAT_PtAPSymbolic;
763:     break;
764:   case MATPRODUCT_RARt:
765:     eventtype = MAT_RARtSymbolic;
766:     break;
767:   case MATPRODUCT_ABC:
768:     eventtype = MAT_MatMatMultSymbolic;
769:     break;
770:   default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
771:   }
772:   mat->ops->productnumeric = NULL;
773:   if (mat->ops->productsymbolic) {
774:     PetscLogEventBegin(eventtype,mat,0,0,0);
775:     (*mat->ops->productsymbolic)(mat);
776:     PetscLogEventEnd(eventtype,mat,0,0,0);
777:   } else missing = PETSC_TRUE;

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

782:     if (mat->product->type == MATPRODUCT_ABC) {
783:       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);
784:     } else {
785:       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);
786:     }
790:   }
791:   return 0;
792: }

794: /*@
795:    MatProductSetFill - Set an expected fill of the matrix product.

797:    Collective on Mat

799:    Input Parameters:
800: +  mat - the matrix product
801: -  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 irrelevant.

803:    Level: intermediate

805: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
806: @*/
807: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
808: {
810:   MatCheckProduct(mat,1);
811:   if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
812:   else mat->product->fill = fill;
813:   return 0;
814: }

816: /*@
817:    MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.

819:    Collective on Mat

821:    Input Parameters:
822: +  mat - the matrix product
823: -  alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHMDEFAULT.

825:    Options Database Key:
826: .  -mat_product_algorithm <algorithm> - Sets the algorithm; use -help for a list
827:     of available algorithms (for instance, scalable, outerproduct, etc.)

829:    Level: intermediate

831: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
832: @*/
833: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
834: {
836:   MatCheckProduct(mat,1);
837:   PetscFree(mat->product->alg);
838:   PetscStrallocpy(alg,&mat->product->alg);
839:   return 0;
840: }

842: /*@
843:    MatProductSetType - Sets a particular matrix product type

845:    Collective on Mat

847:    Input Parameters:
848: +  mat - the matrix
849: -  productype   - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.

851:    Level: intermediate

853: .seealso: MatProductCreate(), MatProductType
854: @*/
855: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
856: {
858:   MatCheckProduct(mat,1);
860:   if (productype != mat->product->type) {
861:     if (mat->product->destroy) {
862:       (*mat->product->destroy)(mat->product->data);
863:     }
864:     mat->product->destroy = NULL;
865:     mat->product->data = NULL;
866:     mat->ops->productsymbolic = NULL;
867:     mat->ops->productnumeric  = NULL;
868:   }
869:   mat->product->type = productype;
870:   return 0;
871: }

873: /*@
874:    MatProductClear - Clears matrix product internal structure.

876:    Collective on Mat

878:    Input Parameters:
879: .  mat - the product matrix

881:    Level: intermediate

883:    Notes: this function should be called to remove any intermediate data used by the product
884:           After having called this function, MatProduct operations can no longer be used on mat
885: @*/
886: PetscErrorCode MatProductClear(Mat mat)
887: {
888:   Mat_Product    *product = mat->product;

891:   if (product) {
892:     MatDestroy(&product->A);
893:     MatDestroy(&product->B);
894:     MatDestroy(&product->C);
895:     PetscFree(product->alg);
896:     MatDestroy(&product->Dwork);
897:     if (product->destroy) {
898:       (*product->destroy)(product->data);
899:     }
900:   }
901:   PetscFree(mat->product);
902:   mat->ops->productsymbolic = NULL;
903:   mat->ops->productnumeric = NULL;
904:   return 0;
905: }

907: /* Create a supporting struct and attach it to the matrix product */
908: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
909: {
910:   Mat_Product    *product=NULL;

914:   PetscNewLog(D,&product);
915:   product->A        = A;
916:   product->B        = B;
917:   product->C        = C;
918:   product->type     = MATPRODUCT_UNSPECIFIED;
919:   product->Dwork    = NULL;
920:   product->api_user = PETSC_FALSE;
921:   product->clear    = PETSC_FALSE;
922:   D->product        = product;

924:   MatProductSetAlgorithm(D,MATPRODUCTALGORITHMDEFAULT);
925:   MatProductSetFill(D,PETSC_DEFAULT);

927:   PetscObjectReference((PetscObject)A);
928:   PetscObjectReference((PetscObject)B);
929:   PetscObjectReference((PetscObject)C);
930:   return 0;
931: }

933: /*@
934:    MatProductCreateWithMat - Setup a given matrix as a matrix product.

936:    Collective on Mat

938:    Input Parameters:
939: +  A - the first matrix
940: .  B - the second matrix
941: .  C - the third matrix (optional)
942: -  D - the matrix which will be used as a product

944:    Output Parameters:
945: .  D - the product matrix

947:    Notes:
948:      Any product data attached to D will be cleared

950:    Level: intermediate

952: .seealso: MatProductCreate(), MatProductClear()
953: @*/
954: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
955: {
958:   MatCheckPreallocated(A,1);

964:   MatCheckPreallocated(B,2);

968:   if (C) {
971:     MatCheckPreallocated(C,3);
974:   }

978:   MatCheckPreallocated(D,4);

982:   /* Create a supporting struct and attach it to D */
983:   MatProductClear(D);
984:   MatProductCreate_Private(A,B,C,D);
985:   return 0;
986: }

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

991:    Collective on Mat

993:    Input Parameters:
994: +  A - the first matrix
995: .  B - the second matrix
996: -  C - the third matrix (optional)

998:    Output Parameters:
999: .  D - the product matrix

1001:    Level: intermediate

1003: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1004: @*/
1005: PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1006: {

1014:   if (C) {
1018:   }

1021:   MatCreate(PetscObjectComm((PetscObject)A),D);
1022:   MatProductCreate_Private(A,B,C,*D);
1023:   return 0;
1024: }

1026: /*
1027:    These are safe basic implementations of ABC, RARt and PtAP
1028:    that do not rely on mat->ops->matmatop function pointers.
1029:    They only use the MatProduct API and are currently used by
1030:    cuSPARSE and KOKKOS-KERNELS backends
1031: */
1032: typedef struct {
1033:   Mat BC;
1034:   Mat ABC;
1035: } MatMatMatPrivate;

1037: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1038: {
1039:   MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;

1041:   MatDestroy(&mmdata->BC);
1042:   MatDestroy(&mmdata->ABC);
1043:   PetscFree(data);
1044:   return 0;
1045: }

1047: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1048: {
1049:   Mat_Product      *product = mat->product;
1050:   MatMatMatPrivate *mmabc;

1052:   MatCheckProduct(mat,1);
1054:   mmabc = (MatMatMatPrivate *)mat->product->data;
1056:   /* use function pointer directly to prevent logging */
1057:   (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1058:   /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1059:   mat->product = mmabc->ABC->product;
1060:   mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1062:   /* use function pointer directly to prevent logging */
1063:   (*mat->ops->productnumeric)(mat);
1064:   mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1065:   mat->product = product;
1066:   return 0;
1067: }

1069: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1070: {
1071:   Mat_Product      *product = mat->product;
1072:   Mat              A, B ,C;
1073:   MatProductType   p1,p2;
1074:   MatMatMatPrivate *mmabc;
1075:   const char       *prefix;

1077:   MatCheckProduct(mat,1);
1079:   MatGetOptionsPrefix(mat,&prefix);
1080:   PetscNew(&mmabc);
1081:   product->data    = mmabc;
1082:   product->destroy = MatDestroy_MatMatMatPrivate;
1083:   switch (product->type) {
1084:   case MATPRODUCT_PtAP:
1085:     p1 = MATPRODUCT_AB;
1086:     p2 = MATPRODUCT_AtB;
1087:     A = product->B;
1088:     B = product->A;
1089:     C = product->B;
1090:     break;
1091:   case MATPRODUCT_RARt:
1092:     p1 = MATPRODUCT_ABt;
1093:     p2 = MATPRODUCT_AB;
1094:     A = product->B;
1095:     B = product->A;
1096:     C = product->B;
1097:     break;
1098:   case MATPRODUCT_ABC:
1099:     p1 = MATPRODUCT_AB;
1100:     p2 = MATPRODUCT_AB;
1101:     A = product->A;
1102:     B = product->B;
1103:     C = product->C;
1104:     break;
1105:   default:
1106:     SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Not for ProductType %s",MatProductTypes[product->type]);
1107:   }
1108:   MatProductCreate(B,C,NULL,&mmabc->BC);
1109:   MatSetOptionsPrefix(mmabc->BC,prefix);
1110:   MatAppendOptionsPrefix(mmabc->BC,"P1_");
1111:   MatProductSetType(mmabc->BC,p1);
1112:   MatProductSetAlgorithm(mmabc->BC,MATPRODUCTALGORITHMDEFAULT);
1113:   MatProductSetFill(mmabc->BC,product->fill);
1114:   mmabc->BC->product->api_user = product->api_user;
1115:   MatProductSetFromOptions(mmabc->BC);
1117:   /* use function pointer directly to prevent logging */
1118:   (*mmabc->BC->ops->productsymbolic)(mmabc->BC);

1120:   MatProductCreate(A,mmabc->BC,NULL,&mmabc->ABC);
1121:   MatSetOptionsPrefix(mmabc->ABC,prefix);
1122:   MatAppendOptionsPrefix(mmabc->ABC,"P2_");
1123:   MatProductSetType(mmabc->ABC,p2);
1124:   MatProductSetAlgorithm(mmabc->ABC,MATPRODUCTALGORITHMDEFAULT);
1125:   MatProductSetFill(mmabc->ABC,product->fill);
1126:   mmabc->ABC->product->api_user = product->api_user;
1127:   MatProductSetFromOptions(mmabc->ABC);
1129:   /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1130:   mat->product = mmabc->ABC->product;
1131:   mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1132:   /* use function pointer directly to prevent logging */
1133:   (*mat->ops->productsymbolic)(mat);
1134:   mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1135:   mat->ops->productsymbolic       = MatProductSymbolic_ABC_Basic;
1136:   mat->ops->productnumeric        = MatProductNumeric_ABC_Basic;
1137:   mat->product                    = product;
1138:   return 0;
1139: }

1141: /*@
1142:    MatProductGetType - Returns the type of the MatProduct.

1144:    Not collective

1146:    Input Parameter:
1147: .  mat - the matrix

1149:    Output Parameter:
1150: .  mtype - the MatProduct type

1152:    Level: intermediate

1154: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductCreate()
1155: @*/
1156: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1157: {
1160:   *mtype = MATPRODUCT_UNSPECIFIED;
1161:   if (mat->product) *mtype = mat->product->type;
1162:   return 0;
1163: }

1165: /*@
1166:    MatProductGetMats - Returns the matrices associated with the MatProduct.

1168:    Not collective

1170:    Input Parameter:
1171: .  mat - the product matrix

1173:    Output Parameters:
1174: +  A - the first matrix
1175: .  B - the second matrix
1176: -  C - the third matrix (optional)

1178:    Level: intermediate

1180: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1181: @*/
1182: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1183: {
1185:   if (A) *A = mat->product ? mat->product->A : NULL;
1186:   if (B) *B = mat->product ? mat->product->B : NULL;
1187:   if (C) *C = mat->product ? mat->product->C : NULL;
1188:   return 0;
1189: }