Actual source code: matproduct.c
petsc-3.14.6 2021-03-30
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: }