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 MATTRANSPOSEVIRTUAL)
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 */
55: (*C->ops->transposematmultnumeric)(P, AP, C);
56: return 0;
57: }
59: static PetscErrorCode MatProductSymbolic_PtAP_Unsafe(Mat C)
60: {
61: Mat_Product *product = C->product;
62: Mat A = product->A, P = product->B, AP;
63: PetscReal fill = product->fill;
65: PetscInfo((PetscObject)C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name);
66: /* AP = A*P */
67: MatProductCreate(A, P, NULL, &AP);
68: MatProductSetType(AP, MATPRODUCT_AB);
69: MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT);
70: MatProductSetFill(AP, fill);
71: MatProductSetFromOptions(AP);
72: MatProductSymbolic(AP);
74: /* C = P^T*AP */
75: MatProductSetType(C, MATPRODUCT_AtB);
76: MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT);
77: product->A = P;
78: product->B = AP;
79: MatProductSetFromOptions(C);
80: MatProductSymbolic(C);
82: /* resume user's original input matrix setting for A and B */
83: product->A = A;
84: product->B = P;
85: product->Dwork = AP;
87: C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
88: return 0;
89: }
91: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
92: {
93: Mat_Product *product = C->product;
94: Mat R = product->B, RA = product->Dwork;
96: /* RA = R*A */
97: MatProductNumeric(RA);
98: /* C = RA*R^T */
99: (*C->ops->mattransposemultnumeric)(RA, R, C);
100: return 0;
101: }
103: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
104: {
105: Mat_Product *product = C->product;
106: Mat A = product->A, R = product->B, RA;
107: PetscReal fill = product->fill;
109: PetscInfo((PetscObject)C, "for A %s, R %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name);
110: /* RA = R*A */
111: MatProductCreate(R, A, NULL, &RA);
112: MatProductSetType(RA, MATPRODUCT_AB);
113: MatProductSetAlgorithm(RA, MATPRODUCTALGORITHMDEFAULT);
114: MatProductSetFill(RA, fill);
115: MatProductSetFromOptions(RA);
116: MatProductSymbolic(RA);
118: /* C = RA*R^T */
119: MatProductSetType(C, MATPRODUCT_ABt);
120: MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT);
121: product->A = RA;
122: MatProductSetFromOptions(C);
123: MatProductSymbolic(C);
125: /* resume user's original input matrix setting for A */
126: product->A = A;
127: product->Dwork = RA; /* save here so it will be destroyed with product C */
128: C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
129: return 0;
130: }
132: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
133: {
134: Mat_Product *product = mat->product;
135: Mat A = product->A, BC = product->Dwork;
137: /* Numeric BC = B*C */
138: MatProductNumeric(BC);
139: /* Numeric mat = A*BC */
140: (*mat->ops->matmultnumeric)(A, BC, mat);
141: return 0;
142: }
144: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
145: {
146: Mat_Product *product = mat->product;
147: Mat B = product->B, C = product->C, BC;
148: PetscReal fill = product->fill;
150: 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);
151: /* Symbolic BC = B*C */
152: MatProductCreate(B, C, NULL, &BC);
153: MatProductSetType(BC, MATPRODUCT_AB);
154: MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT);
155: MatProductSetFill(BC, fill);
156: MatProductSetFromOptions(BC);
157: MatProductSymbolic(BC);
159: /* Symbolic mat = A*BC */
160: MatProductSetType(mat, MATPRODUCT_AB);
161: MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT);
162: product->B = BC;
163: product->Dwork = BC;
164: MatProductSetFromOptions(mat);
165: MatProductSymbolic(mat);
167: /* resume user's original input matrix setting for B */
168: product->B = B;
169: mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
170: return 0;
171: }
173: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
174: {
175: Mat_Product *product = mat->product;
177: switch (product->type) {
178: case MATPRODUCT_PtAP:
179: MatProductSymbolic_PtAP_Unsafe(mat);
180: break;
181: case MATPRODUCT_RARt:
182: MatProductSymbolic_RARt_Unsafe(mat);
183: break;
184: case MATPRODUCT_ABC:
185: MatProductSymbolic_ABC_Unsafe(mat);
186: break;
187: default:
188: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
189: }
190: return 0;
191: }
193: /* ----------------------------------------------- */
194: /*@C
195: MatProductReplaceMats - Replace input matrices for a matrix product.
197: Collective on Mat
199: Input Parameters:
200: + A - the matrix or NULL if not being replaced
201: . B - the matrix or NULL if not being replaced
202: . C - the matrix or NULL if not being replaced
203: - D - the matrix product
205: Level: intermediate
207: Note:
208: To reuse the symbolic phase, the input matrices must have exactly the same data structure as the replaced one.
209: If the type of any of the input matrices is different than what was previously used, or their symmetry flag changed but
210: the symbolic phase took advantage of their symmetry, the product is cleared and `MatProductSetFromOptions()` and `MatProductSymbolic()` are invoked again.
212: .seealso: `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic().` `MatProductClear()`
213: @*/
214: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
215: {
216: Mat_Product *product;
217: PetscBool flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;
220: MatCheckProduct(D, 4);
221: product = D->product;
222: if (A) {
224: PetscObjectReference((PetscObject)A);
225: PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA);
226: MatIsSymmetricKnown(A, &isset, &issym);
227: if (product->symbolic_used_the_fact_A_is_symmetric && isset && !issym) { /* symbolic was built around a symmetric A, but the new A is not anymore */
228: flgA = PETSC_FALSE;
229: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
230: }
231: MatDestroy(&product->A);
232: product->A = A;
233: }
234: if (B) {
236: PetscObjectReference((PetscObject)B);
237: PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB);
238: MatIsSymmetricKnown(B, &isset, &issym);
239: if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
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: MatIsSymmetricKnown(C, &isset, &issym);
251: if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
252: flgC = PETSC_FALSE;
253: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
254: }
255: MatDestroy(&product->C);
256: product->C = C;
257: }
258: /* Any of the replaced mats is of a different type, reset */
259: if (!flgA || !flgB || !flgC) {
260: if (D->product->destroy) (*D->product->destroy)(D->product->data);
261: D->product->destroy = NULL;
262: D->product->data = NULL;
263: if (D->ops->productnumeric || D->ops->productsymbolic) {
264: MatProductSetFromOptions(D);
265: MatProductSymbolic(D);
266: }
267: }
268: return 0;
269: }
271: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
272: {
273: Mat_Product *product = C->product;
274: Mat A = product->A, B = product->B;
275: PetscInt k, K = B->cmap->N;
276: PetscBool t = PETSC_TRUE, iscuda = PETSC_FALSE;
277: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
278: char *Btype = NULL, *Ctype = NULL;
280: switch (product->type) {
281: case MATPRODUCT_AB:
282: t = PETSC_FALSE;
283: case MATPRODUCT_AtB:
284: break;
285: default:
286: 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) PetscStrcmp(vtype, VECSEQCUDA, &iscuda);
294: if (!iscuda) PetscStrcmp(vtype, VECMPICUDA, &iscuda);
295: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
296: PetscStrallocpy(((PetscObject)B)->type_name, &Btype);
297: PetscStrallocpy(((PetscObject)C)->type_name, &Ctype);
298: MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B);
299: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
300: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
301: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
302: }
303: MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C);
304: } else { /* Make sure we have up-to-date data on the CPU */
305: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
306: Bcpu = B->boundtocpu;
307: Ccpu = C->boundtocpu;
308: #endif
309: MatBindToCPU(B, PETSC_TRUE);
310: MatBindToCPU(C, PETSC_TRUE);
311: }
312: }
313: for (k = 0; k < K; k++) {
314: Vec x, y;
316: MatDenseGetColumnVecRead(B, k, &x);
317: MatDenseGetColumnVecWrite(C, k, &y);
318: if (t) {
319: MatMultTranspose(A, x, y);
320: } else {
321: MatMult(A, x, y);
322: }
323: MatDenseRestoreColumnVecRead(B, k, &x);
324: MatDenseRestoreColumnVecWrite(C, k, &y);
325: }
326: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
327: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
328: if (PetscDefined(HAVE_CUDA)) {
329: if (iscuda) {
330: MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B);
331: MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C);
332: } else {
333: MatBindToCPU(B, Bcpu);
334: MatBindToCPU(C, Ccpu);
335: }
336: }
337: PetscFree(Btype);
338: PetscFree(Ctype);
339: return 0;
340: }
342: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
343: {
344: Mat_Product *product = C->product;
345: Mat A = product->A, B = product->B;
346: PetscBool isdense;
348: switch (product->type) {
349: case MATPRODUCT_AB:
350: MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N);
351: break;
352: case MATPRODUCT_AtB:
353: MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N);
354: break;
355: default:
356: 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);
357: }
358: PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, "");
359: if (!isdense) {
360: MatSetType(C, ((PetscObject)B)->type_name);
361: /* If matrix type of C was not set or not dense, we need to reset the pointer */
362: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
363: }
364: C->ops->productnumeric = MatProductNumeric_X_Dense;
365: MatSetUp(C);
366: return 0;
367: }
369: /* a single driver to query the dispatching */
370: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
371: {
372: Mat_Product *product = mat->product;
373: PetscInt Am, An, Bm, Bn, Cm, Cn;
374: Mat A = product->A, B = product->B, C = product->C;
375: const char *const Bnames[] = {"B", "R", "P"};
376: const char *bname;
377: PetscErrorCode (*fA)(Mat);
378: PetscErrorCode (*fB)(Mat);
379: PetscErrorCode (*fC)(Mat);
380: PetscErrorCode (*f)(Mat) = NULL;
382: mat->ops->productsymbolic = NULL;
383: mat->ops->productnumeric = NULL;
384: if (product->type == MATPRODUCT_UNSPECIFIED) return 0;
388: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
389: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
390: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
391: else bname = Bnames[0];
393: /* Check matrices sizes */
394: Am = A->rmap->N;
395: An = A->cmap->N;
396: Bm = B->rmap->N;
397: Bn = B->cmap->N;
398: Cm = C ? C->rmap->N : 0;
399: Cn = C ? C->cmap->N : 0;
400: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
401: PetscInt t = Bn;
402: Bn = Bm;
403: Bm = t;
404: }
405: if (product->type == MATPRODUCT_AtB) {
406: PetscInt t = An;
407: An = Am;
408: Am = t;
409: }
411: MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
413: MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);
415: fA = A->ops->productsetfromoptions;
416: fB = B->ops->productsetfromoptions;
417: fC = C ? C->ops->productsetfromoptions : fA;
418: if (C) {
419: 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);
420: } else {
421: PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name);
422: }
423: if (fA == fB && fA == fC && fA) {
424: PetscInfo(mat, " matching op\n");
425: (*fA)(mat);
426: }
427: /* We may have found f but it did not succeed */
428: if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
429: char mtypes[256];
430: PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes));
431: PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes));
432: PetscStrlcat(mtypes, "_", sizeof(mtypes));
433: PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes));
434: if (C) {
435: PetscStrlcat(mtypes, "_", sizeof(mtypes));
436: PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes));
437: }
438: PetscStrlcat(mtypes, "_C", sizeof(mtypes));
440: PetscObjectQueryFunction((PetscObject)A, mtypes, &f);
441: PetscInfo(mat, " querying %s from A? %p\n", mtypes, f);
442: if (!f) {
443: PetscObjectQueryFunction((PetscObject)B, mtypes, &f);
444: PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f);
445: }
446: if (!f && C) {
447: PetscObjectQueryFunction((PetscObject)C, mtypes, &f);
448: PetscInfo(mat, " querying %s from C? %p\n", mtypes, f);
449: }
450: if (f) (*f)(mat);
452: /* We may have found f but it did not succeed */
453: /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
454: if (!mat->ops->productsymbolic) {
455: PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes));
456: PetscObjectQueryFunction((PetscObject)A, mtypes, &f);
457: PetscInfo(mat, " querying %s from A? %p\n", mtypes, f);
458: if (!f) {
459: PetscObjectQueryFunction((PetscObject)B, mtypes, &f);
460: PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f);
461: }
462: if (!f && C) {
463: PetscObjectQueryFunction((PetscObject)C, mtypes, &f);
464: PetscInfo(mat, " querying %s from C? %p\n", mtypes, f);
465: }
466: }
467: if (f) (*f)(mat);
468: }
470: /* We may have found f but it did not succeed */
471: if (!mat->ops->productsymbolic) {
472: /* we can still compute the product if B is of type dense */
473: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
474: PetscBool isdense;
476: PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, "");
477: if (isdense) {
478: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
479: PetscInfo(mat, " using basic looping over columns of a dense matrix\n");
480: }
481: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
482: /*
483: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
484: the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
485: before computing the symbolic phase
486: */
487: PetscInfo(mat, " symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n");
488: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
489: }
490: }
491: if (!mat->ops->productsymbolic) PetscInfo(mat, " symbolic product is not supported\n");
492: return 0;
493: }
495: /*@C
496: MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product where the type, the algorithm etc are determined from the options database.
498: Logically Collective on Mat
500: Input Parameter:
501: . mat - the matrix
503: Options Database Keys:
504: . -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
506: Level: intermediate
508: .seealso: `Mat`, `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
509: @*/
510: PetscErrorCode MatProductSetFromOptions(Mat mat)
511: {
513: MatCheckProduct(mat, 1);
515: PetscObjectOptionsBegin((PetscObject)mat);
516: PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL);
517: PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()");
518: PetscOptionsEnd();
519: MatProductSetFromOptions_Private(mat);
521: return 0;
522: }
524: /*@C
525: MatProductView - View the private `Mat_Product` algorithm object within a matrix
527: Logically Collective
529: Input Parameters:
530: + mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
531: - viewer - where `mat` should be reviewed
533: Level: intermediate
535: .seealso: `Mat`, `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
536: @*/
537: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
538: {
540: if (!mat->product) return 0;
541: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer);
544: if (mat->product->view) (*mat->product->view)(mat, viewer);
545: return 0;
546: }
548: /* ----------------------------------------------- */
549: /* these are basic implementations relying on the old function pointers
550: * they are dangerous and should be removed in the future */
551: PetscErrorCode MatProductNumeric_AB(Mat mat)
552: {
553: Mat_Product *product = mat->product;
554: 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;
565: (*mat->ops->transposematmultnumeric)(A, B, mat);
566: return 0;
567: }
569: PetscErrorCode MatProductNumeric_ABt(Mat mat)
570: {
571: Mat_Product *product = mat->product;
572: Mat A = product->A, B = product->B;
574: (*mat->ops->mattransposemultnumeric)(A, B, mat);
575: return 0;
576: }
578: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
579: {
580: Mat_Product *product = mat->product;
581: Mat A = product->A, B = product->B;
583: (*mat->ops->ptapnumeric)(A, B, mat);
584: return 0;
585: }
587: PetscErrorCode MatProductNumeric_RARt(Mat mat)
588: {
589: Mat_Product *product = mat->product;
590: Mat A = product->A, B = product->B;
592: (*mat->ops->rartnumeric)(A, B, mat);
593: return 0;
594: }
596: PetscErrorCode MatProductNumeric_ABC(Mat mat)
597: {
598: Mat_Product *product = mat->product;
599: Mat A = product->A, B = product->B, C = product->C;
601: (*mat->ops->matmatmultnumeric)(A, B, C, mat);
602: return 0;
603: }
605: /* ----------------------------------------------- */
607: /*@
608: MatProductNumeric - Compute a matrix product with numerical values.
610: Collective
612: Input/Output Parameter:
613: . mat - the matrix holding the product
615: Level: intermediate
617: Note:
618: `MatProductSymbolic()` must have been called on mat before calling this function
620: .seealso: `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
621: @*/
622: PetscErrorCode MatProductNumeric(Mat mat)
623: {
624: PetscLogEvent eventtype = -1;
625: PetscBool missing = PETSC_FALSE;
628: MatCheckProduct(mat, 1);
629: switch (mat->product->type) {
630: case MATPRODUCT_AB:
631: eventtype = MAT_MatMultNumeric;
632: break;
633: case MATPRODUCT_AtB:
634: eventtype = MAT_TransposeMatMultNumeric;
635: break;
636: case MATPRODUCT_ABt:
637: eventtype = MAT_MatTransposeMultNumeric;
638: break;
639: case MATPRODUCT_PtAP:
640: eventtype = MAT_PtAPNumeric;
641: break;
642: case MATPRODUCT_RARt:
643: eventtype = MAT_RARtNumeric;
644: break;
645: case MATPRODUCT_ABC:
646: eventtype = MAT_MatMatMultNumeric;
647: break;
648: default:
649: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
650: }
652: if (mat->ops->productnumeric) {
653: PetscLogEventBegin(eventtype, mat, 0, 0, 0);
654: PetscUseTypeMethod(mat, productnumeric);
655: PetscLogEventEnd(eventtype, mat, 0, 0, 0);
656: } else missing = PETSC_TRUE;
658: if (missing || !mat->product) {
659: char errstr[256];
661: if (mat->product->type == MATPRODUCT_ABC) {
662: 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);
663: } else {
664: 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);
665: }
668: }
670: if (mat->product->clear) MatProductClear(mat);
671: PetscObjectStateIncrease((PetscObject)mat);
672: return 0;
673: }
675: /* ----------------------------------------------- */
676: /* these are basic implementations relying on the old function pointers
677: * they are dangerous and should be removed in the future */
678: PetscErrorCode MatProductSymbolic_AB(Mat mat)
679: {
680: Mat_Product *product = mat->product;
681: Mat A = product->A, B = product->B;
683: (*mat->ops->matmultsymbolic)(A, B, product->fill, mat);
684: mat->ops->productnumeric = MatProductNumeric_AB;
685: return 0;
686: }
688: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
689: {
690: Mat_Product *product = mat->product;
691: Mat A = product->A, B = product->B;
693: (*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat);
694: mat->ops->productnumeric = MatProductNumeric_AtB;
695: return 0;
696: }
698: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
699: {
700: Mat_Product *product = mat->product;
701: Mat A = product->A, B = product->B;
703: (*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat);
704: mat->ops->productnumeric = MatProductNumeric_ABt;
705: return 0;
706: }
708: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
709: {
710: Mat_Product *product = mat->product;
711: Mat A = product->A, B = product->B, C = product->C;
713: (*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat);
714: mat->ops->productnumeric = MatProductNumeric_ABC;
715: return 0;
716: }
718: /* ----------------------------------------------- */
720: /*@
721: MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical product done with
722: `MatProductNumeric()`
724: Collective
726: Input/Output Parameter:
727: . mat - the matrix to hold a product
729: Level: intermediate
731: Note:
732: `MatProductSetFromOptions()` must have been called on mat before calling this function
734: .seealso: `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
735: @*/
736: PetscErrorCode MatProductSymbolic(Mat mat)
737: {
738: PetscLogEvent eventtype = -1;
739: PetscBool missing = PETSC_FALSE;
742: MatCheckProduct(mat, 1);
744: switch (mat->product->type) {
745: case MATPRODUCT_AB:
746: eventtype = MAT_MatMultSymbolic;
747: break;
748: case MATPRODUCT_AtB:
749: eventtype = MAT_TransposeMatMultSymbolic;
750: break;
751: case MATPRODUCT_ABt:
752: eventtype = MAT_MatTransposeMultSymbolic;
753: break;
754: case MATPRODUCT_PtAP:
755: eventtype = MAT_PtAPSymbolic;
756: break;
757: case MATPRODUCT_RARt:
758: eventtype = MAT_RARtSymbolic;
759: break;
760: case MATPRODUCT_ABC:
761: eventtype = MAT_MatMatMultSymbolic;
762: break;
763: default:
764: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
765: }
766: mat->ops->productnumeric = NULL;
767: if (mat->ops->productsymbolic) {
768: PetscLogEventBegin(eventtype, mat, 0, 0, 0);
769: PetscUseTypeMethod(mat, productsymbolic);
770: PetscLogEventEnd(eventtype, mat, 0, 0, 0);
771: } else missing = PETSC_TRUE;
773: if (missing || !mat->product || !mat->ops->productnumeric) {
774: char errstr[256];
776: if (mat->product->type == MATPRODUCT_ABC) {
777: 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);
778: } else {
779: 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);
780: }
783: }
784: return 0;
785: }
787: /*@
788: MatProductSetFill - Set an expected fill of the matrix product.
790: Collective on Mat
792: Input Parameters:
793: + mat - the matrix product result matrix
794: - 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 value is not used.
796: Level: intermediate
798: .seealso: `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
799: @*/
800: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
801: {
803: MatCheckProduct(mat, 1);
804: if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
805: else mat->product->fill = fill;
806: return 0;
807: }
809: /*@
810: MatProductSetAlgorithm - Requests a particular algorithm for a matrix product computation that will perform to compute the given matrix
812: Collective
814: Input Parameters:
815: + mat - the matrix product
816: - alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
818: Options Database Key:
819: . -mat_product_algorithm <algorithm> - Sets the algorithm; use -help for a list
820: of available algorithms (for instance, scalable, outerproduct, etc.)
822: Level: intermediate
824: .seealso: `Mat`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`
825: @*/
826: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
827: {
829: MatCheckProduct(mat, 1);
830: PetscFree(mat->product->alg);
831: PetscStrallocpy(alg, &mat->product->alg);
832: return 0;
833: }
835: /*@
836: MatProductSetType - Sets a particular matrix product type to be used to compute the given matrix
838: Collective
840: Input Parameters:
841: + mat - the matrix
842: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`.
844: Level: intermediate
846: Note:
847: The small t represents the transpose operation.
849: .seealso: `Mat`, `MatProductCreate()`, `MatProductType`, `MatProductType`,
850: `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
851: @*/
852: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
853: {
855: MatCheckProduct(mat, 1);
857: if (productype != mat->product->type) {
858: if (mat->product->destroy) (*mat->product->destroy)(mat->product->data);
859: mat->product->destroy = NULL;
860: mat->product->data = NULL;
861: mat->ops->productsymbolic = NULL;
862: mat->ops->productnumeric = NULL;
863: }
864: mat->product->type = productype;
865: return 0;
866: }
868: /*@
869: MatProductClear - Clears matrix product internal datastructures.
871: Collective
873: Input Parameters:
874: . mat - the product matrix
876: Level: intermediate
878: Notes:
879: This function should be called to remove any intermediate data used to compute the matrix to free up memory.
881: After having called this function, matrix-matrix operations can no longer be used on mat
883: .seealso: `Mat`, `MatProductCreate()`
884: @*/
885: PetscErrorCode MatProductClear(Mat mat)
886: {
887: Mat_Product *product = mat->product;
890: if (product) {
891: MatDestroy(&product->A);
892: MatDestroy(&product->B);
893: MatDestroy(&product->C);
894: PetscFree(product->alg);
895: MatDestroy(&product->Dwork);
896: if (product->destroy) (*product->destroy)(product->data);
897: }
898: PetscFree(mat->product);
899: mat->ops->productsymbolic = NULL;
900: mat->ops->productnumeric = NULL;
901: return 0;
902: }
904: /* Create a supporting struct and attach it to the matrix product */
905: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
906: {
907: Mat_Product *product = NULL;
911: PetscNew(&product);
912: product->A = A;
913: product->B = B;
914: product->C = C;
915: product->type = MATPRODUCT_UNSPECIFIED;
916: product->Dwork = NULL;
917: product->api_user = PETSC_FALSE;
918: product->clear = PETSC_FALSE;
919: D->product = product;
921: MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT);
922: MatProductSetFill(D, PETSC_DEFAULT);
924: PetscObjectReference((PetscObject)A);
925: PetscObjectReference((PetscObject)B);
926: PetscObjectReference((PetscObject)C);
927: return 0;
928: }
930: /*@
931: MatProductCreateWithMat - Setup a given matrix as a matrix product of other matrices
933: Collective on Mat
935: Input Parameters:
936: + A - the first matrix
937: . B - the second matrix
938: . C - the third matrix (optional, use `NULL` if not needed)
939: - D - the matrix which will be used to hold the product
941: Level: intermediate
943: Notes:
944: Use `MatProductCreate()` if the matrix you wish computed (the D matrix) does not already exist
946: See `MatProductCreate()` for details on the usage of the MatProduct routines
948: Any product data currently attached to `D` will be cleared
950: .seealso: `Mat`, `MatProductCreate()`, `MatProductClear()`
951: @*/
952: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
953: {
956: MatCheckPreallocated(A, 1);
962: MatCheckPreallocated(B, 2);
966: if (C) {
969: MatCheckPreallocated(C, 3);
972: }
976: MatCheckPreallocated(D, 4);
980: /* Create a supporting struct and attach it to D */
981: MatProductClear(D);
982: MatProductCreate_Private(A, B, C, D);
983: return 0;
984: }
986: /*@
987: MatProductCreate - create a matrix to hold the result of a matrix-matrix product operation
989: Collective on A
991: Input Parameters:
992: + A - the first matrix
993: . B - the second matrix
994: - C - the third matrix (optional)
996: Output Parameters:
997: . D - the product matrix
999: Level: intermediate
1001: Example of Usage:
1002: .vb
1003: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1004: MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1005: MatProductSetAlgorithm(D, alg)
1006: MatProductSetFill(D,fill)
1007: MatProductSetFromOptions(D)
1008: MatProductSymbolic(D)
1009: MatProductNumeric(D)
1010: Change numerical values in some of the matrices
1011: MatProductNumeric(D)
1012: .ve
1014: Notes:
1015: Use `MatProductCreateWithMat()` if the matrix you wish computed, the D matrix, already exists.
1017: The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure
1019: Developer Note:
1020: It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1021: Is there error checking for it?
1023: .seealso: `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`
1024: @*/
1025: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1026: {
1034: if (C) {
1038: }
1041: MatCreate(PetscObjectComm((PetscObject)A), D);
1042: /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1043: MatProductCreate_Private(A, B, C, *D);
1044: return 0;
1045: }
1047: /*
1048: These are safe basic implementations of ABC, RARt and PtAP
1049: that do not rely on mat->ops->matmatop function pointers.
1050: They only use the MatProduct API and are currently used by
1051: cuSPARSE and KOKKOS-KERNELS backends
1052: */
1053: typedef struct {
1054: Mat BC;
1055: Mat ABC;
1056: } MatMatMatPrivate;
1058: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1059: {
1060: MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;
1062: MatDestroy(&mmdata->BC);
1063: MatDestroy(&mmdata->ABC);
1064: PetscFree(data);
1065: return 0;
1066: }
1068: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1069: {
1070: Mat_Product *product = mat->product;
1071: MatMatMatPrivate *mmabc;
1073: MatCheckProduct(mat, 1);
1075: mmabc = (MatMatMatPrivate *)mat->product->data;
1077: /* use function pointer directly to prevent logging */
1078: (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1079: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1080: mat->product = mmabc->ABC->product;
1081: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1082: /* use function pointer directly to prevent logging */
1083: PetscUseTypeMethod(mat, productnumeric);
1084: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1085: mat->product = product;
1086: return 0;
1087: }
1089: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1090: {
1091: Mat_Product *product = mat->product;
1092: Mat A, B, C;
1093: MatProductType p1, p2;
1094: MatMatMatPrivate *mmabc;
1095: const char *prefix;
1097: MatCheckProduct(mat, 1);
1099: MatGetOptionsPrefix(mat, &prefix);
1100: PetscNew(&mmabc);
1101: product->data = mmabc;
1102: product->destroy = MatDestroy_MatMatMatPrivate;
1103: switch (product->type) {
1104: case MATPRODUCT_PtAP:
1105: p1 = MATPRODUCT_AB;
1106: p2 = MATPRODUCT_AtB;
1107: A = product->B;
1108: B = product->A;
1109: C = product->B;
1110: break;
1111: case MATPRODUCT_RARt:
1112: p1 = MATPRODUCT_ABt;
1113: p2 = MATPRODUCT_AB;
1114: A = product->B;
1115: B = product->A;
1116: C = product->B;
1117: break;
1118: case MATPRODUCT_ABC:
1119: p1 = MATPRODUCT_AB;
1120: p2 = MATPRODUCT_AB;
1121: A = product->A;
1122: B = product->B;
1123: C = product->C;
1124: break;
1125: default:
1126: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1127: }
1128: MatProductCreate(B, C, NULL, &mmabc->BC);
1129: MatSetOptionsPrefix(mmabc->BC, prefix);
1130: MatAppendOptionsPrefix(mmabc->BC, "P1_");
1131: MatProductSetType(mmabc->BC, p1);
1132: MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT);
1133: MatProductSetFill(mmabc->BC, product->fill);
1134: mmabc->BC->product->api_user = product->api_user;
1135: MatProductSetFromOptions(mmabc->BC);
1137: /* use function pointer directly to prevent logging */
1138: (*mmabc->BC->ops->productsymbolic)(mmabc->BC);
1140: MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC);
1141: MatSetOptionsPrefix(mmabc->ABC, prefix);
1142: MatAppendOptionsPrefix(mmabc->ABC, "P2_");
1143: MatProductSetType(mmabc->ABC, p2);
1144: MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT);
1145: MatProductSetFill(mmabc->ABC, product->fill);
1146: mmabc->ABC->product->api_user = product->api_user;
1147: MatProductSetFromOptions(mmabc->ABC);
1149: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1150: mat->product = mmabc->ABC->product;
1151: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1152: /* use function pointer directly to prevent logging */
1153: PetscUseTypeMethod(mat, productsymbolic);
1154: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1155: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1156: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1157: mat->product = product;
1158: return 0;
1159: }
1161: /*@
1162: MatProductGetType - Returns the type of matrix-matrix product associated with the given matrix.
1164: Not collective
1166: Input Parameter:
1167: . mat - the matrix
1169: Output Parameter:
1170: . mtype - the `MatProductType`
1172: Level: intermediate
1174: .seealso: `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1175: @*/
1176: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1177: {
1180: *mtype = MATPRODUCT_UNSPECIFIED;
1181: if (mat->product) *mtype = mat->product->type;
1182: return 0;
1183: }
1185: /*@
1186: MatProductGetMats - Returns the matrices associated with the matrix-matrix product this matrix can receive
1188: Not collective
1190: Input Parameter:
1191: . mat - the product matrix
1193: Output Parameters:
1194: + A - the first matrix
1195: . B - the second matrix
1196: - C - the third matrix (optional)
1198: Level: intermediate
1200: .seealso: `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1201: @*/
1202: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1203: {
1205: if (A) *A = mat->product ? mat->product->A : NULL;
1206: if (B) *B = mat->product ? mat->product->B : NULL;
1207: if (C) *C = mat->product ? mat->product->C : NULL;
1208: return 0;
1209: }