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: {
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 was previously used, or their symmetry changed but
226: the symbolic phase took advantage of their symmetry, the product is cleared and MatProductSetFromOptions()/MatProductSymbolic() are invoked again.
228: .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
229: @*/
230: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
231: {
233: Mat_Product *product;
234: PetscBool flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;
238: MatCheckProduct(D,4);
239: product = D->product;
240: if (A) {
242: PetscObjectReference((PetscObject)A);
243: PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);
244: 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 */
245: flgA = PETSC_FALSE;
246: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
247: }
248: MatDestroy(&product->A);
249: product->A = A;
250: }
251: if (B) {
253: PetscObjectReference((PetscObject)B);
254: PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);
255: if (product->symbolic_used_the_fact_B_is_symmetric && !B->symmetric) {
256: flgB = PETSC_FALSE;
257: product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
258: }
259: MatDestroy(&product->B);
260: product->B = B;
261: }
262: if (C) {
264: PetscObjectReference((PetscObject)C);
265: PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);
266: if (product->symbolic_used_the_fact_C_is_symmetric && !C->symmetric) {
267: flgC = PETSC_FALSE;
268: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
269: }
270: MatDestroy(&product->C);
271: product->C = C;
272: }
273: /* Any of the replaced mats is of a different type, reset */
274: if (!flgA || !flgB || !flgC) {
275: if (D->product->destroy) {
276: (*D->product->destroy)(D->product->data);
277: }
278: D->product->destroy = NULL;
279: D->product->data = NULL;
280: if (D->ops->productnumeric || D->ops->productsymbolic) {
281: MatProductSetFromOptions(D);
282: MatProductSymbolic(D);
283: }
284: }
285: return(0);
286: }
288: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
289: {
291: Mat_Product *product = C->product;
292: Mat A = product->A, B = product->B;
293: PetscInt k, K = B->cmap->N;
294: PetscBool t = PETSC_TRUE,iscuda = PETSC_FALSE;
295: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
296: char *Btype = NULL,*Ctype = NULL;
299: switch (product->type) {
300: case MATPRODUCT_AB:
301: t = PETSC_FALSE;
302: case MATPRODUCT_AtB:
303: break;
304: 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);
305: }
306: if (PetscDefined(HAVE_CUDA)) {
307: VecType vtype;
309: MatGetVecType(A,&vtype);
310: PetscStrcmp(vtype,VECCUDA,&iscuda);
311: if (!iscuda) {
312: PetscStrcmp(vtype,VECSEQCUDA,&iscuda);
313: }
314: if (!iscuda) {
315: PetscStrcmp(vtype,VECMPICUDA,&iscuda);
316: }
317: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
318: PetscStrallocpy(((PetscObject)B)->type_name,&Btype);
319: PetscStrallocpy(((PetscObject)C)->type_name,&Ctype);
320: MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B);
321: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
322: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
324: }
325: MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C);
326: } else { /* Make sure we have up-to-date data on the CPU */
327: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
328: Bcpu = B->boundtocpu;
329: Ccpu = C->boundtocpu;
330: #endif
331: MatBindToCPU(B,PETSC_TRUE);
332: MatBindToCPU(C,PETSC_TRUE);
333: }
334: }
335: for (k=0;k<K;k++) {
336: Vec x,y;
338: MatDenseGetColumnVecRead(B,k,&x);
339: MatDenseGetColumnVecWrite(C,k,&y);
340: if (t) {
341: MatMultTranspose(A,x,y);
342: } else {
343: MatMult(A,x,y);
344: }
345: MatDenseRestoreColumnVecRead(B,k,&x);
346: MatDenseRestoreColumnVecWrite(C,k,&y);
347: }
348: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
349: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
350: if (PetscDefined(HAVE_CUDA)) {
351: if (iscuda) {
352: MatConvert(B,Btype,MAT_INPLACE_MATRIX,&B);
353: MatConvert(C,Ctype,MAT_INPLACE_MATRIX,&C);
354: } else {
355: MatBindToCPU(B,Bcpu);
356: MatBindToCPU(C,Ccpu);
357: }
358: }
359: PetscFree(Btype);
360: PetscFree(Ctype);
361: return(0);
362: }
364: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
365: {
367: Mat_Product *product = C->product;
368: Mat A = product->A, B = product->B;
369: PetscBool isdense;
372: switch (product->type) {
373: case MATPRODUCT_AB:
374: MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
375: break;
376: case MATPRODUCT_AtB:
377: MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
378: break;
379: 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);
380: }
381: PetscObjectBaseTypeCompareAny((PetscObject)C,&isdense,MATSEQDENSE,MATMPIDENSE,"");
382: if (!isdense) {
383: MatSetType(C,((PetscObject)B)->type_name);
384: /* If matrix type of C was not set or not dense, we need to reset the pointer */
385: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
386: }
387: C->ops->productnumeric = MatProductNumeric_X_Dense;
388: MatSetUp(C);
389: return(0);
390: }
392: /* a single driver to query the dispatching */
393: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
394: {
395: PetscErrorCode ierr;
396: Mat_Product *product = mat->product;
397: PetscInt Am,An,Bm,Bn,Cm,Cn;
398: Mat A = product->A,B = product->B,C = product->C;
399: const char* const Bnames[] = { "B", "R", "P" };
400: const char* bname;
401: PetscErrorCode (*fA)(Mat);
402: PetscErrorCode (*fB)(Mat);
403: PetscErrorCode (*fC)(Mat);
404: PetscErrorCode (*f)(Mat)=NULL;
407: mat->ops->productsymbolic = NULL;
408: mat->ops->productnumeric = NULL;
409: if (product->type == MATPRODUCT_UNSPECIFIED) return(0);
410: if (!A) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing A mat");
411: if (!B) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing B mat");
412: if (product->type == MATPRODUCT_ABC && !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing C mat");
413: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
414: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
415: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
416: else bname = Bnames[0];
418: /* Check matrices sizes */
419: Am = A->rmap->N;
420: An = A->cmap->N;
421: Bm = B->rmap->N;
422: Bn = B->cmap->N;
423: Cm = C ? C->rmap->N : 0;
424: Cn = C ? C->cmap->N : 0;
425: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
426: if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
427: 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);
428: 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);
430: fA = A->ops->productsetfromoptions;
431: fB = B->ops->productsetfromoptions;
432: fC = C ? C->ops->productsetfromoptions : fA;
433: if (C) {
434: 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);
435: } else {
436: PetscInfo4(mat,"MatProductType %s for A %s, %s %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name);
437: }
438: if (fA == fB && fA == fC && fA) {
439: PetscInfo(mat," matching op\n");
440: (*fA)(mat);
441: } else {
442: /* query MatProductSetFromOptions_Atype_Btype_Ctype */
443: char mtypes[256];
444: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
445: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
446: PetscStrlcat(mtypes,"_",sizeof(mtypes));
447: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
448: if (C) {
449: PetscStrlcat(mtypes,"_",sizeof(mtypes));
450: PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));
451: }
452: PetscStrlcat(mtypes,"_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: if (f) { (*f)(mat); }
466: /* We may have found f but it did not succeed */
467: /* some matrices (i.e. MATTRANSPOSE, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
468: if (!mat->ops->productsymbolic) {
469: PetscStrncpy(mtypes,"MatProductSetFromOptions_anytype_C",sizeof(mtypes));
470: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
471: PetscInfo2(mat," querying %s from A? %p\n",mtypes,f);
472: if (!f) {
473: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
474: PetscInfo3(mat," querying %s from %s? %p\n",mtypes,bname,f);
475: }
476: if (!f && C) {
477: PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
478: PetscInfo2(mat," querying %s from C? %p\n",mtypes,f);
479: }
480: }
481: if (f) { (*f)(mat); }
482: }
484: /* We may have found f but it did not succeed */
485: if (!mat->ops->productsymbolic) {
486: /* we can still compute the product if B is of type dense */
487: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
488: PetscBool isdense;
490: PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
491: if (isdense) {
493: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
494: PetscInfo(mat," using basic looping over columns of a dense matrix\n");
495: }
496: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
497: /*
498: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
499: the compination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
500: before computing the symbolic phase
501: */
502: PetscInfo(mat," symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n");
503: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
504: }
505: }
506: if (!mat->ops->productsymbolic) {
507: PetscInfo(mat," symbolic product is not supported\n");
508: }
509: return(0);
510: }
512: /*@C
513: MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.
515: Logically Collective on Mat
517: Input Parameter:
518: . mat - the matrix
520: Options Database Keys:
521: . -mat_product_clear - Clear intermediate data structures after MatProductNumeric() has been called
523: Level: intermediate
525: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
526: @*/
527: PetscErrorCode MatProductSetFromOptions(Mat mat)
528: {
533: MatCheckProduct(mat,1);
534: if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot call MatProductSetFromOptions with already present data");
535: PetscObjectOptionsBegin((PetscObject)mat);
536: PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);
537: PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");
538: PetscOptionsEnd();
539: MatProductSetFromOptions_Private(mat);
540: if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after setup phase");
541: return(0);
542: }
544: /*@C
545: MatProductView - View a MatProduct
547: Logically Collective on Mat
549: Input Parameter:
550: . mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()
552: Level: intermediate
554: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
555: @*/
556: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
557: {
562: if (!mat->product) return(0);
563: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);}
566: if (mat->product->view) {
567: (*mat->product->view)(mat,viewer);
568: }
569: return(0);
570: }
572: /* ----------------------------------------------- */
573: /* these are basic implementations relying on the old function pointers
574: * they are dangerous and should be removed in the future */
575: PetscErrorCode MatProductNumeric_AB(Mat mat)
576: {
578: Mat_Product *product = mat->product;
579: Mat A=product->A,B=product->B;
582: 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);
583: (*mat->ops->matmultnumeric)(A,B,mat);
584: return(0);
585: }
587: PetscErrorCode MatProductNumeric_AtB(Mat mat)
588: {
590: Mat_Product *product = mat->product;
591: Mat A=product->A,B=product->B;
594: 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);
595: (*mat->ops->transposematmultnumeric)(A,B,mat);
596: return(0);
597: }
599: PetscErrorCode MatProductNumeric_ABt(Mat mat)
600: {
602: Mat_Product *product = mat->product;
603: Mat A=product->A,B=product->B;
606: 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);
607: (*mat->ops->mattransposemultnumeric)(A,B,mat);
608: return(0);
609: }
611: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
612: {
614: Mat_Product *product = mat->product;
615: Mat A=product->A,B=product->B;
618: 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);
619: (*mat->ops->ptapnumeric)(A,B,mat);
620: return(0);
621: }
623: PetscErrorCode MatProductNumeric_RARt(Mat mat)
624: {
626: Mat_Product *product = mat->product;
627: Mat A=product->A,B=product->B;
630: 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);
631: (*mat->ops->rartnumeric)(A,B,mat);
632: return(0);
633: }
635: PetscErrorCode MatProductNumeric_ABC(Mat mat)
636: {
638: Mat_Product *product = mat->product;
639: Mat A=product->A,B=product->B,C=product->C;
642: 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);
643: (*mat->ops->matmatmultnumeric)(A,B,C,mat);
644: return(0);
645: }
647: /* ----------------------------------------------- */
649: /*@
650: MatProductNumeric - Implement a matrix product with numerical values.
652: Collective on Mat
654: Input/Output Parameter:
655: . mat - the matrix holding the product
657: Level: intermediate
659: Notes: MatProductSymbolic() must have been called on mat before calling this function
661: .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
662: @*/
663: PetscErrorCode MatProductNumeric(Mat mat)
664: {
666: PetscLogEvent eventtype=-1;
670: MatCheckProduct(mat,1);
671: /* log event */
672: switch (mat->product->type) {
673: case MATPRODUCT_AB:
674: eventtype = MAT_MatMultNumeric;
675: break;
676: case MATPRODUCT_AtB:
677: eventtype = MAT_TransposeMatMultNumeric;
678: break;
679: case MATPRODUCT_ABt:
680: eventtype = MAT_MatTransposeMultNumeric;
681: break;
682: case MATPRODUCT_PtAP:
683: eventtype = MAT_PtAPNumeric;
684: break;
685: case MATPRODUCT_RARt:
686: eventtype = MAT_RARtNumeric;
687: break;
688: case MATPRODUCT_ABC:
689: eventtype = MAT_MatMatMultNumeric;
690: break;
691: default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
692: }
693: if (mat->ops->productnumeric) {
694: PetscLogEventBegin(eventtype,mat,0,0,0);
695: (*mat->ops->productnumeric)(mat);
696: PetscLogEventEnd(eventtype,mat,0,0,0);
697: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSymbolic() first");
698: if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after numeric phase");
699: if (mat->product->clear) {
700: MatProductClear(mat);
701: }
702: PetscObjectStateIncrease((PetscObject)mat);
703: return(0);
704: }
706: /* ----------------------------------------------- */
707: /* these are basic implementations relying on the old function pointers
708: * they are dangerous and should be removed in the future */
709: PetscErrorCode MatProductSymbolic_AB(Mat mat)
710: {
712: Mat_Product *product = mat->product;
713: Mat A=product->A,B=product->B;
716: if (!mat->ops->matmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
717: (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);
718: mat->ops->productnumeric = MatProductNumeric_AB;
719: return(0);
720: }
722: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
723: {
725: Mat_Product *product = mat->product;
726: Mat A=product->A,B=product->B;
729: if (!mat->ops->transposematmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
730: (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
731: mat->ops->productnumeric = MatProductNumeric_AtB;
732: return(0);
733: }
735: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
736: {
738: Mat_Product *product = mat->product;
739: Mat A=product->A,B=product->B;
742: if (!mat->ops->mattransposemultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
743: (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
744: mat->ops->productnumeric = MatProductNumeric_ABt;
745: return(0);
746: }
748: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
749: {
751: Mat_Product *product = mat->product;
752: Mat A=product->A,B=product->B,C=product->C;
755: if (!mat->ops->matmatmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
756: (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
757: mat->ops->productnumeric = MatProductNumeric_ABC;
758: return(0);
759: }
761: /* ----------------------------------------------- */
763: /*@
764: MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.
766: Collective on Mat
768: Input/Output Parameter:
769: . mat - the matrix to hold a product
771: Level: intermediate
773: Notes: MatProductSetFromOptions() must have been called on mat before calling this function
775: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
776: @*/
777: PetscErrorCode MatProductSymbolic(Mat mat)
778: {
780: PetscLogEvent eventtype=-1;
784: MatCheckProduct(mat,1);
785: if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot run symbolic phase. Product data not empty");
786: /* log event */
787: switch (mat->product->type) {
788: case MATPRODUCT_AB:
789: eventtype = MAT_MatMultSymbolic;
790: break;
791: case MATPRODUCT_AtB:
792: eventtype = MAT_TransposeMatMultSymbolic;
793: break;
794: case MATPRODUCT_ABt:
795: eventtype = MAT_MatTransposeMultSymbolic;
796: break;
797: case MATPRODUCT_PtAP:
798: eventtype = MAT_PtAPSymbolic;
799: break;
800: case MATPRODUCT_RARt:
801: eventtype = MAT_RARtSymbolic;
802: break;
803: case MATPRODUCT_ABC:
804: eventtype = MAT_MatMatMultSymbolic;
805: break;
806: default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
807: }
809: mat->ops->productnumeric = NULL;
810: if (mat->ops->productsymbolic) {
811: PetscLogEventBegin(eventtype,mat,0,0,0);
812: (*mat->ops->productsymbolic)(mat);
813: PetscLogEventEnd(eventtype,mat,0,0,0);
814: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSetFromOptions() first");
815: if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after symbolic phase");
816: if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Symbolic phase did not specify the numeric phase");
817: return(0);
818: }
820: /*@
821: MatProductSetFill - Set an expected fill of the matrix product.
823: Collective on Mat
825: Input Parameters:
826: + mat - the matrix product
827: - 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.
829: Level: intermediate
831: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
832: @*/
833: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
834: {
837: MatCheckProduct(mat,1);
838: if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
839: else mat->product->fill = fill;
840: return(0);
841: }
843: /*@
844: MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.
846: Collective on Mat
848: Input Parameters:
849: + mat - the matrix product
850: - alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.
852: Level: intermediate
854: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
855: @*/
856: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
857: {
862: MatCheckProduct(mat,1);
863: PetscFree(mat->product->alg);
864: PetscStrallocpy(alg,&mat->product->alg);
865: return(0);
866: }
868: /*@
869: MatProductSetType - Sets a particular matrix product type
871: Collective on Mat
873: Input Parameters:
874: + mat - the matrix
875: - productype - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.
877: Level: intermediate
879: .seealso: MatProductCreate(), MatProductType
880: @*/
881: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
882: {
887: MatCheckProduct(mat,1);
889: if (productype != mat->product->type) {
890: if (mat->product->destroy) {
891: (*mat->product->destroy)(mat->product->data);
892: }
893: mat->product->destroy = NULL;
894: mat->product->data = NULL;
895: mat->ops->productsymbolic = NULL;
896: mat->ops->productnumeric = NULL;
897: }
898: mat->product->type = productype;
899: return(0);
900: }
902: /*@
903: MatProductClear - Clears matrix product internal structure.
905: Collective on Mat
907: Input Parameters:
908: . mat - the product matrix
910: Level: intermediate
912: Notes: this function should be called to remove any intermediate data used by the product
913: After having called this function, MatProduct operations can no longer be used on mat
914: @*/
915: PetscErrorCode MatProductClear(Mat mat)
916: {
918: Mat_Product *product = mat->product;
922: if (product) {
923: MatDestroy(&product->A);
924: MatDestroy(&product->B);
925: MatDestroy(&product->C);
926: PetscFree(product->alg);
927: MatDestroy(&product->Dwork);
928: if (product->destroy) {
929: (*product->destroy)(product->data);
930: }
931: }
932: PetscFree(mat->product);
933: mat->ops->productsymbolic = NULL;
934: mat->ops->productnumeric = NULL;
935: return(0);
936: }
938: /* Create a supporting struct and attach it to the matrix product */
939: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
940: {
942: Mat_Product *product=NULL;
946: if (D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product already present");
947: PetscNewLog(D,&product);
948: product->A = A;
949: product->B = B;
950: product->C = C;
951: product->type = MATPRODUCT_UNSPECIFIED;
952: product->Dwork = NULL;
953: product->api_user = PETSC_FALSE;
954: product->clear = PETSC_FALSE;
955: D->product = product;
957: MatProductSetAlgorithm(D,MATPRODUCTALGORITHM_DEFAULT);
958: MatProductSetFill(D,PETSC_DEFAULT);
960: PetscObjectReference((PetscObject)A);
961: PetscObjectReference((PetscObject)B);
962: PetscObjectReference((PetscObject)C);
963: return(0);
964: }
966: /*@
967: MatProductCreateWithMat - Setup a given matrix as a matrix product.
969: Collective on Mat
971: Input Parameters:
972: + A - the first matrix
973: . B - the second matrix
974: . C - the third matrix (optional)
975: - D - the matrix which will be used as a product
977: Output Parameters:
978: . D - the product matrix
980: Notes:
981: Any product data attached to D will be cleared
983: Level: intermediate
985: .seealso: MatProductCreate(), MatProductClear()
986: @*/
987: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
988: {
994: MatCheckPreallocated(A,1);
995: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
996: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1000: MatCheckPreallocated(B,2);
1001: if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1002: if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1004: if (C) {
1007: MatCheckPreallocated(C,3);
1008: if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1009: if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1010: }
1014: MatCheckPreallocated(D,4);
1015: if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1016: if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1018: /* Create a supporting struct and attach it to D */
1019: MatProductClear(D);
1020: MatProductCreate_Private(A,B,C,D);
1021: return(0);
1022: }
1024: /*@
1025: MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.
1027: Collective on Mat
1029: Input Parameters:
1030: + A - the first matrix
1031: . B - the second matrix
1032: - C - the third matrix (optional)
1034: Output Parameters:
1035: . D - the product matrix
1037: Level: intermediate
1039: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1040: @*/
1041: PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1042: {
1050: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix A");
1051: if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix B");
1053: if (C) {
1056: if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix C");
1057: }
1060: MatCreate(PetscObjectComm((PetscObject)A),D);
1061: MatProductCreate_Private(A,B,C,*D);
1062: return(0);
1063: }
1065: /*
1066: These are safe basic implementations of ABC, RARt and PtAP
1067: that do not rely on mat->ops->matmatop function pointers.
1068: They only use the MatProduct API and are currently used by
1069: cuSPARSE and KOKKOS-KERNELS backends
1070: */
1071: typedef struct {
1072: Mat BC;
1073: Mat ABC;
1074: } MatMatMatPrivate;
1076: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1077: {
1078: PetscErrorCode ierr;
1079: MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;
1082: MatDestroy(&mmdata->BC);
1083: MatDestroy(&mmdata->ABC);
1084: PetscFree(data);
1085: return(0);
1086: }
1088: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1089: {
1090: PetscErrorCode ierr;
1091: Mat_Product *product = mat->product;
1092: MatMatMatPrivate *mmabc;
1095: MatCheckProduct(mat,1);
1096: if (!mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data empty");
1097: mmabc = (MatMatMatPrivate *)mat->product->data;
1098: if (!mmabc->BC->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
1099: /* use function pointer directly to prevent logging */
1100: (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1101: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1102: mat->product = mmabc->ABC->product;
1103: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1104: if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
1105: /* use function pointer directly to prevent logging */
1106: (*mat->ops->productnumeric)(mat);
1107: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1108: mat->product = product;
1109: return(0);
1110: }
1112: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1113: {
1114: PetscErrorCode ierr;
1115: Mat_Product *product = mat->product;
1116: Mat A, B ,C;
1117: MatProductType p1,p2;
1118: MatMatMatPrivate *mmabc;
1119: const char *prefix;
1122: MatCheckProduct(mat,1);
1123: if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Product data not empty");
1124: MatGetOptionsPrefix(mat,&prefix);
1125: PetscNew(&mmabc);
1126: product->data = mmabc;
1127: product->destroy = MatDestroy_MatMatMatPrivate;
1128: switch (product->type) {
1129: case MATPRODUCT_PtAP:
1130: p1 = MATPRODUCT_AB;
1131: p2 = MATPRODUCT_AtB;
1132: A = product->B;
1133: B = product->A;
1134: C = product->B;
1135: break;
1136: case MATPRODUCT_RARt:
1137: p1 = MATPRODUCT_ABt;
1138: p2 = MATPRODUCT_AB;
1139: A = product->B;
1140: B = product->A;
1141: C = product->B;
1142: break;
1143: case MATPRODUCT_ABC:
1144: p1 = MATPRODUCT_AB;
1145: p2 = MATPRODUCT_AB;
1146: A = product->A;
1147: B = product->B;
1148: C = product->C;
1149: break;
1150: default:
1151: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Not for ProductType %s",MatProductTypes[product->type]);
1152: }
1153: MatProductCreate(B,C,NULL,&mmabc->BC);
1154: MatSetOptionsPrefix(mmabc->BC,prefix);
1155: MatAppendOptionsPrefix(mmabc->BC,"P1_");
1156: MatProductSetType(mmabc->BC,p1);
1157: MatProductSetAlgorithm(mmabc->BC,MATPRODUCTALGORITHM_DEFAULT);
1158: MatProductSetFill(mmabc->BC,product->fill);
1159: mmabc->BC->product->api_user = product->api_user;
1160: MatProductSetFromOptions(mmabc->BC);
1161: 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);
1162: /* use function pointer directly to prevent logging */
1163: (*mmabc->BC->ops->productsymbolic)(mmabc->BC);
1165: MatProductCreate(A,mmabc->BC,NULL,&mmabc->ABC);
1166: MatSetOptionsPrefix(mmabc->ABC,prefix);
1167: MatAppendOptionsPrefix(mmabc->ABC,"P2_");
1168: MatProductSetType(mmabc->ABC,p2);
1169: MatProductSetAlgorithm(mmabc->ABC,MATPRODUCTALGORITHM_DEFAULT);
1170: MatProductSetFill(mmabc->ABC,product->fill);
1171: mmabc->ABC->product->api_user = product->api_user;
1172: MatProductSetFromOptions(mmabc->ABC);
1173: 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);
1174: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1175: mat->product = mmabc->ABC->product;
1176: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1177: /* use function pointer directly to prevent logging */
1178: (*mat->ops->productsymbolic)(mat);
1179: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1180: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1181: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1182: mat->product = product;
1183: return(0);
1184: }