Actual source code: matproduct.c
1: /*
2: Routines for matrix products. Calling procedure:
4: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
5: MatProductSetType(D, MATPRODUCT_AB/AtB/ABt/PtAP/RARt/ABC)
6: MatProductSetAlgorithm(D, alg)
7: MatProductSetFill(D,fill)
8: MatProductSetFromOptions(D)
9: -> MatProductSetFromOptions_Private(D)
10: # Check matrix global sizes
11: if the matrices have the same setfromoptions routine, use it
12: if not, try:
13: -> Query MatProductSetFromOptions_Atype_Btype_Ctype_C(D) from A, B and C (in order)
14: if found -> run the specific setup that must set the symbolic operation (these callbacks should never fail)
15: if callback not found or no symbolic operation set
16: -> Query MatProductSetFromOptions_anytype_C(D) from A, B and C (in order) (e.g, matrices may have inner matrices like MATTRANSPOSEVIRTUAL)
17: if dispatch found but combination still not present do
18: -> check if B is dense and product type AtB or AB -> if true, basic looping of dense columns
19: -> check if triple product (PtAP, RARt or ABC) -> if true, set the Basic routines
21: # The setfromoptions calls MatProductSetFromOptions_Atype_Btype_Ctype should
22: # Check matrix local sizes for mpi matrices
23: # Set default algorithm
24: # Get runtime option
25: # Set D->ops->productsymbolic = MatProductSymbolic_productype_Atype_Btype_Ctype if found
27: MatProductSymbolic(D)
28: # Call MatProductSymbolic_productype_Atype_Btype_Ctype()
29: the callback must set the numeric phase D->ops->productnumeric = MatProductNumeric_productype_Atype_Btype_Ctype
31: MatProductNumeric(D)
32: # Call the numeric phase
34: # The symbolic phases are allowed to set extra data structures and attach those to the product
35: # this additional data can be reused between multiple numeric phases with the same matrices
36: # if not needed, call
37: MatProductClear(D)
38: */
40: #include <petsc/private/matimpl.h>
42: const char *const MatProductTypes[] = {"UNSPECIFIED", "AB", "AtB", "ABt", "PtAP", "RARt", "ABC"};
44: /* these are basic implementations relying on the old function pointers
45: * they are dangerous and should be removed in the future */
46: static PetscErrorCode MatProductNumeric_PtAP_Unsafe(Mat C)
47: {
48: Mat_Product *product = C->product;
49: Mat P = product->B, AP = product->Dwork;
51: PetscFunctionBegin;
52: /* AP = A*P */
53: PetscCall(MatProductNumeric(AP));
54: /* C = P^T*AP */
55: product->type = MATPRODUCT_AtB;
56: PetscCall((*C->ops->transposematmultnumeric)(P, AP, C));
57: product->type = MATPRODUCT_PtAP;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode MatProductSymbolic_PtAP_Unsafe(Mat C)
62: {
63: Mat_Product *product = C->product;
64: Mat A = product->A, P = product->B, AP;
65: PetscReal fill = product->fill;
67: PetscFunctionBegin;
68: PetscCall(PetscInfo(C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
69: /* AP = A*P */
70: PetscCall(MatProductCreate(A, P, NULL, &AP));
71: PetscCall(MatProductSetType(AP, MATPRODUCT_AB));
72: PetscCall(MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT));
73: PetscCall(MatProductSetFill(AP, fill));
74: PetscCall(MatProductSetFromOptions(AP));
75: PetscCall(MatProductSymbolic(AP));
77: /* C = P^T*AP */
78: PetscCall(MatProductSetType(C, MATPRODUCT_AtB));
79: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
80: product->A = P;
81: product->B = AP;
82: PetscCall(MatProductSetFromOptions(C));
83: PetscCall(MatProductSymbolic(C));
85: /* resume user's original input matrix setting for A and B */
86: product->type = MATPRODUCT_PtAP;
87: product->A = A;
88: product->B = P;
89: product->Dwork = AP;
91: C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
96: {
97: Mat_Product *product = C->product;
98: Mat R = product->B, RA = product->Dwork;
100: PetscFunctionBegin;
101: /* RA = R*A */
102: PetscCall(MatProductNumeric(RA));
103: /* C = RA*R^T */
104: product->type = MATPRODUCT_ABt;
105: PetscCall((*C->ops->mattransposemultnumeric)(RA, R, C));
106: product->type = MATPRODUCT_RARt;
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
111: {
112: Mat_Product *product = C->product;
113: Mat A = product->A, R = product->B, RA;
114: PetscReal fill = product->fill;
116: PetscFunctionBegin;
117: PetscCall(PetscInfo(C, "for A %s, R %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
118: /* RA = R*A */
119: PetscCall(MatProductCreate(R, A, NULL, &RA));
120: PetscCall(MatProductSetType(RA, MATPRODUCT_AB));
121: PetscCall(MatProductSetAlgorithm(RA, MATPRODUCTALGORITHMDEFAULT));
122: PetscCall(MatProductSetFill(RA, fill));
123: PetscCall(MatProductSetFromOptions(RA));
124: PetscCall(MatProductSymbolic(RA));
126: /* C = RA*R^T */
127: PetscCall(MatProductSetType(C, MATPRODUCT_ABt));
128: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
129: product->A = RA;
130: PetscCall(MatProductSetFromOptions(C));
131: PetscCall(MatProductSymbolic(C));
133: /* resume user's original input matrix setting for A */
134: product->type = MATPRODUCT_RARt;
135: product->A = A;
136: product->Dwork = RA; /* save here so it will be destroyed with product C */
137: C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
142: {
143: Mat_Product *product = mat->product;
144: Mat A = product->A, BC = product->Dwork;
146: PetscFunctionBegin;
147: /* Numeric BC = B*C */
148: PetscCall(MatProductNumeric(BC));
149: /* Numeric mat = A*BC */
150: product->type = MATPRODUCT_AB;
151: PetscCall((*mat->ops->matmultnumeric)(A, BC, mat));
152: product->type = MATPRODUCT_ABC;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
157: {
158: Mat_Product *product = mat->product;
159: Mat B = product->B, C = product->C, BC;
160: PetscReal fill = product->fill;
162: PetscFunctionBegin;
163: PetscCall(PetscInfo(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));
164: /* Symbolic BC = B*C */
165: PetscCall(MatProductCreate(B, C, NULL, &BC));
166: PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
167: PetscCall(MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT));
168: PetscCall(MatProductSetFill(BC, fill));
169: PetscCall(MatProductSetFromOptions(BC));
170: PetscCall(MatProductSymbolic(BC));
172: /* Symbolic mat = A*BC */
173: PetscCall(MatProductSetType(mat, MATPRODUCT_AB));
174: PetscCall(MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT));
175: product->B = BC;
176: product->Dwork = BC;
177: PetscCall(MatProductSetFromOptions(mat));
178: PetscCall(MatProductSymbolic(mat));
180: /* resume user's original input matrix setting for B */
181: product->type = MATPRODUCT_ABC;
182: product->B = B;
183: mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
188: {
189: Mat_Product *product = mat->product;
191: PetscFunctionBegin;
192: switch (product->type) {
193: case MATPRODUCT_PtAP:
194: PetscCall(MatProductSymbolic_PtAP_Unsafe(mat));
195: break;
196: case MATPRODUCT_RARt:
197: PetscCall(MatProductSymbolic_RARt_Unsafe(mat));
198: break;
199: case MATPRODUCT_ABC:
200: PetscCall(MatProductSymbolic_ABC_Unsafe(mat));
201: break;
202: default:
203: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: MatProductReplaceMats - Replace the input matrices for the matrix-matrix product operation inside the computed matrix
211: Collective
213: Input Parameters:
214: + A - the matrix or `NULL` if not being replaced
215: . B - the matrix or `NULL` if not being replaced
216: . C - the matrix or `NULL` if not being replaced
217: - D - the matrix whose values are computed via a matrix-matrix product operation
219: Level: intermediate
221: Note:
222: To reuse the symbolic phase, the input matrices must have exactly the same data structure as the replaced one.
223: If the type of any of the input matrices is different than what was previously used, or their symmetry flag changed but
224: the symbolic phase took advantage of their symmetry, the product is cleared and `MatProductSetFromOptions()`
225: and `MatProductSymbolic()` are invoked again.
227: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic()`, `MatProductClear()`
228: @*/
229: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
230: {
231: Mat_Product *product;
232: PetscBool flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;
234: PetscFunctionBegin;
236: MatCheckProduct(D, 4);
237: product = D->product;
238: if (A) {
240: PetscCall(PetscObjectReference((PetscObject)A));
241: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA));
242: PetscCall(MatIsSymmetricKnown(A, &isset, &issym));
243: 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 */
244: flgA = PETSC_FALSE;
245: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
246: }
247: PetscCall(MatDestroy(&product->A));
248: product->A = A;
249: }
250: if (B) {
252: PetscCall(PetscObjectReference((PetscObject)B));
253: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB));
254: PetscCall(MatIsSymmetricKnown(B, &isset, &issym));
255: if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
256: flgB = PETSC_FALSE;
257: product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
258: }
259: PetscCall(MatDestroy(&product->B));
260: product->B = B;
261: }
262: if (C) {
264: PetscCall(PetscObjectReference((PetscObject)C));
265: PetscCall(PetscObjectTypeCompare((PetscObject)product->C, ((PetscObject)C)->type_name, &flgC));
266: PetscCall(MatIsSymmetricKnown(C, &isset, &issym));
267: if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
268: flgC = PETSC_FALSE;
269: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
270: }
271: PetscCall(MatDestroy(&product->C));
272: product->C = C;
273: }
274: /* Any of the replaced mats is of a different type, reset */
275: if (!flgA || !flgB || !flgC) {
276: if (D->product->destroy) PetscCall((*D->product->destroy)(&D->product->data));
277: D->product->destroy = NULL;
278: D->product->data = NULL;
279: if (D->ops->productnumeric || D->ops->productsymbolic) {
280: PetscCall(MatProductSetFromOptions(D));
281: PetscCall(MatProductSymbolic(D));
282: }
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
288: {
289: Mat_Product *product = C->product;
290: Mat A = product->A, B = product->B;
291: PetscInt k, K = B->cmap->N;
292: PetscBool t = PETSC_TRUE, iscuda = PETSC_FALSE;
293: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
294: char *Btype = NULL, *Ctype = NULL;
296: PetscFunctionBegin;
297: switch (product->type) {
298: case MATPRODUCT_AB:
299: t = PETSC_FALSE;
300: case MATPRODUCT_AtB:
301: break;
302: default:
303: 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);
304: }
305: if (PetscDefined(HAVE_CUDA)) {
306: VecType vtype;
308: PetscCall(MatGetVecType(A, &vtype));
309: PetscCall(PetscStrcmp(vtype, VECCUDA, &iscuda));
310: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECSEQCUDA, &iscuda));
311: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECMPICUDA, &iscuda));
312: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
313: PetscCall(PetscStrallocpy(((PetscObject)B)->type_name, &Btype));
314: PetscCall(PetscStrallocpy(((PetscObject)C)->type_name, &Ctype));
315: PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
316: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
317: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
318: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
319: }
320: PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
321: } else { /* Make sure we have up-to-date data on the CPU */
322: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
323: Bcpu = B->boundtocpu;
324: Ccpu = C->boundtocpu;
325: #endif
326: PetscCall(MatBindToCPU(B, PETSC_TRUE));
327: PetscCall(MatBindToCPU(C, PETSC_TRUE));
328: }
329: }
330: for (k = 0; k < K; k++) {
331: Vec x, y;
333: PetscCall(MatDenseGetColumnVecRead(B, k, &x));
334: PetscCall(MatDenseGetColumnVecWrite(C, k, &y));
335: if (t) {
336: PetscCall(MatMultTranspose(A, x, y));
337: } else {
338: PetscCall(MatMult(A, x, y));
339: }
340: PetscCall(MatDenseRestoreColumnVecRead(B, k, &x));
341: PetscCall(MatDenseRestoreColumnVecWrite(C, k, &y));
342: }
343: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
344: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
345: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
346: if (PetscDefined(HAVE_CUDA)) {
347: if (iscuda) {
348: PetscCall(MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B));
349: PetscCall(MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C));
350: } else {
351: PetscCall(MatBindToCPU(B, Bcpu));
352: PetscCall(MatBindToCPU(C, Ccpu));
353: }
354: }
355: PetscCall(PetscFree(Btype));
356: PetscCall(PetscFree(Ctype));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
361: {
362: Mat_Product *product = C->product;
363: Mat A = product->A, B = product->B;
364: PetscBool isdense;
366: PetscFunctionBegin;
367: switch (product->type) {
368: case MATPRODUCT_AB:
369: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
370: break;
371: case MATPRODUCT_AtB:
372: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
373: break;
374: default:
375: 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);
376: }
377: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
378: if (!isdense) {
379: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
380: /* If matrix type of C was not set or not dense, we need to reset the pointer */
381: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
382: }
383: C->ops->productnumeric = MatProductNumeric_X_Dense;
384: PetscCall(MatSetUp(C));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /* a single driver to query the dispatching */
389: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
390: {
391: Mat_Product *product = mat->product;
392: PetscInt Am, An, Bm, Bn, Cm, Cn;
393: Mat A = product->A, B = product->B, C = product->C;
394: const char *const Bnames[] = {"B", "R", "P"};
395: const char *bname;
396: PetscErrorCode (*fA)(Mat);
397: PetscErrorCode (*fB)(Mat);
398: PetscErrorCode (*fC)(Mat);
399: PetscErrorCode (*f)(Mat) = NULL;
401: PetscFunctionBegin;
402: mat->ops->productsymbolic = NULL;
403: mat->ops->productnumeric = NULL;
404: if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(PETSC_SUCCESS);
405: PetscCheck(A, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing A mat");
406: PetscCheck(B, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing B mat");
407: PetscCheck(product->type != MATPRODUCT_ABC || C, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing C mat");
408: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
409: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
410: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
411: else bname = Bnames[0];
413: /* Check matrices sizes */
414: Am = A->rmap->N;
415: An = A->cmap->N;
416: Bm = B->rmap->N;
417: Bn = B->cmap->N;
418: Cm = C ? C->rmap->N : 0;
419: Cn = C ? C->cmap->N : 0;
420: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
421: PetscInt t = Bn;
422: Bn = Bm;
423: Bm = t;
424: }
425: if (product->type == MATPRODUCT_AtB) An = Am;
427: PetscCheck(An == Bm, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_SIZ, "Matrix dimensions of A and %s are incompatible for MatProductType %s: A %" PetscInt_FMT "x%" PetscInt_FMT ", %s %" PetscInt_FMT "x%" PetscInt_FMT, bname,
428: MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
429: PetscCheck(!Cm || Cm == Bn, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_SIZ, "Matrix dimensions of B and C are incompatible for MatProductType %s: B %" PetscInt_FMT "x%" PetscInt_FMT ", C %" PetscInt_FMT "x%" PetscInt_FMT,
430: MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);
432: fA = A->ops->productsetfromoptions;
433: fB = B->ops->productsetfromoptions;
434: fC = C ? C->ops->productsetfromoptions : fA;
435: if (C) {
436: PetscCall(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));
437: } else {
438: PetscCall(PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name));
439: }
440: if (fA == fB && fA == fC && fA) {
441: PetscCall(PetscInfo(mat, " matching op\n"));
442: PetscCall((*fA)(mat));
443: }
444: /* We may have found f but it did not succeed */
445: if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
446: char mtypes[256];
447: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes)));
448: PetscCall(PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes)));
449: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
450: PetscCall(PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes)));
451: if (C) {
452: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
453: PetscCall(PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes)));
454: }
455: PetscCall(PetscStrlcat(mtypes, "_C", sizeof(mtypes)));
456: #if defined(__clang__)
457: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
458: #elif defined(__GNUC__) || defined(__GNUG__)
459: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
460: #endif
461: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
462: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
463: if (!f) {
464: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
465: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
466: }
467: if (!f && C) {
468: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
469: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
470: }
471: if (f) PetscCall((*f)(mat));
473: /* We may have found f but it did not succeed */
474: /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
475: if (!mat->ops->productsymbolic) {
476: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes)));
477: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
478: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
479: if (!f) {
480: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
481: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
482: }
483: if (!f && C) {
484: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
485: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
486: }
487: }
488: if (f) PetscCall((*f)(mat));
489: }
490: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
491: /* We may have found f but it did not succeed */
492: if (!mat->ops->productsymbolic) {
493: /* we can still compute the product if B is of type dense */
494: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
495: PetscBool isdense;
497: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
498: if (isdense) {
499: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
500: PetscCall(PetscInfo(mat, " using basic looping over columns of a dense matrix\n"));
501: }
502: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
503: /*
504: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
505: the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
506: before computing the symbolic phase
507: */
508: PetscCall(PetscInfo(mat, " symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n"));
509: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
510: }
511: }
512: if (!mat->ops->productsymbolic) PetscCall(PetscInfo(mat, " symbolic product is not supported\n"));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product operation where the type,
518: the algorithm etc are determined from the options database.
520: Logically Collective
522: Input Parameter:
523: . mat - the matrix whose values are computed via a matrix-matrix product operation
525: Options Database Keys:
526: + -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
527: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm` for possible values
528: - -mat_product_algorithm_backend_cpu - Use the CPU to perform the computation even if the matrix is a GPU matrix
530: Level: intermediate
532: Note:
533: The `-mat_product_clear` option reduces memory usage but means that the matrix cannot be re-used for a matrix-matrix product operation
535: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`,
536: `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductAlgorithm`
537: @*/
538: PetscErrorCode MatProductSetFromOptions(Mat mat)
539: {
540: PetscFunctionBegin;
542: MatCheckProduct(mat, 1);
543: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot call MatProductSetFromOptions() with already present data");
544: mat->product->setfromoptionscalled = PETSC_TRUE;
545: PetscObjectOptionsBegin((PetscObject)mat);
546: PetscCall(PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL));
547: PetscCall(PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()"));
548: PetscOptionsEnd();
549: PetscCall(MatProductSetFromOptions_Private(mat));
550: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing product after setup phase");
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: /*@
555: MatProductView - View the private matrix-matrix algorithm object within a matrix
557: Logically Collective
559: Input Parameters:
560: + mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
561: - viewer - where the information on the matrix-matrix algorithm of `mat` should be reviewed
563: Level: intermediate
565: Developer Note:
566: Shouldn't this information be printed from an appropriate `MatView()` with perhaps certain formats set?
568: .seealso: [](ch_matrices), `MatProductType`, `Mat`, `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
569: @*/
570: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
571: {
572: PetscFunctionBegin;
574: if (!mat->product) PetscFunctionReturn(PETSC_SUCCESS);
575: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer));
577: PetscCheckSameComm(mat, 1, viewer, 2);
578: if (mat->product->view) PetscCall((*mat->product->view)(mat, viewer));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /* these are basic implementations relying on the old function pointers
583: * they are dangerous and should be removed in the future */
584: PetscErrorCode MatProductNumeric_AB(Mat mat)
585: {
586: Mat_Product *product = mat->product;
587: Mat A = product->A, B = product->B;
589: PetscFunctionBegin;
590: PetscCall((*mat->ops->matmultnumeric)(A, B, mat));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: PetscErrorCode MatProductNumeric_AtB(Mat mat)
595: {
596: Mat_Product *product = mat->product;
597: Mat A = product->A, B = product->B;
599: PetscFunctionBegin;
600: PetscCall((*mat->ops->transposematmultnumeric)(A, B, mat));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: PetscErrorCode MatProductNumeric_ABt(Mat mat)
605: {
606: Mat_Product *product = mat->product;
607: Mat A = product->A, B = product->B;
609: PetscFunctionBegin;
610: PetscCall((*mat->ops->mattransposemultnumeric)(A, B, mat));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
615: {
616: Mat_Product *product = mat->product;
617: Mat A = product->A, B = product->B;
619: PetscFunctionBegin;
620: PetscCall((*mat->ops->ptapnumeric)(A, B, mat));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: PetscErrorCode MatProductNumeric_RARt(Mat mat)
625: {
626: Mat_Product *product = mat->product;
627: Mat A = product->A, B = product->B;
629: PetscFunctionBegin;
630: PetscCall((*mat->ops->rartnumeric)(A, B, mat));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: PetscErrorCode MatProductNumeric_ABC(Mat mat)
635: {
636: Mat_Product *product = mat->product;
637: Mat A = product->A, B = product->B, C = product->C;
639: PetscFunctionBegin;
640: PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: /*@
645: MatProductNumeric - Compute a matrix-matrix product operation with the numerical values
647: Collective
649: Input/Output Parameter:
650: . mat - the matrix whose values are computed via a matrix-matrix product operation
652: Level: intermediate
654: Note:
655: `MatProductSymbolic()` must have been called on `mat` before calling this function
657: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
658: @*/
659: PetscErrorCode MatProductNumeric(Mat mat)
660: {
661: PetscLogEvent eventtype = -1;
663: PetscFunctionBegin;
665: MatCheckProduct(mat, 1);
666: switch (mat->product->type) {
667: case MATPRODUCT_AB:
668: eventtype = MAT_MatMultNumeric;
669: break;
670: case MATPRODUCT_AtB:
671: eventtype = MAT_TransposeMatMultNumeric;
672: break;
673: case MATPRODUCT_ABt:
674: eventtype = MAT_MatTransposeMultNumeric;
675: break;
676: case MATPRODUCT_PtAP:
677: eventtype = MAT_PtAPNumeric;
678: break;
679: case MATPRODUCT_RARt:
680: eventtype = MAT_RARtNumeric;
681: break;
682: case MATPRODUCT_ABC:
683: eventtype = MAT_MatMatMultNumeric;
684: break;
685: default:
686: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
687: }
689: if (mat->ops->productnumeric) {
690: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
691: PetscUseTypeMethod(mat, productnumeric);
692: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
693: } else if (mat->product) {
694: char errstr[256];
696: if (mat->product->type == MATPRODUCT_ABC) {
697: PetscCall(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));
698: } else {
699: PetscCall(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));
700: }
701: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified numeric phase for product %s", errstr);
702: }
703: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after numeric phase for product");
705: if (mat->product->clear) PetscCall(MatProductClear(mat));
706: PetscCall(PetscObjectStateIncrease((PetscObject)mat));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /* these are basic implementations relying on the old function pointers
711: * they are dangerous and should be removed in the future */
712: PetscErrorCode MatProductSymbolic_AB(Mat mat)
713: {
714: Mat_Product *product = mat->product;
715: Mat A = product->A, B = product->B;
717: PetscFunctionBegin;
718: PetscCall((*mat->ops->matmultsymbolic)(A, B, product->fill, mat));
719: mat->ops->productnumeric = MatProductNumeric_AB;
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
724: {
725: Mat_Product *product = mat->product;
726: Mat A = product->A, B = product->B;
728: PetscFunctionBegin;
729: PetscCall((*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat));
730: mat->ops->productnumeric = MatProductNumeric_AtB;
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
735: {
736: Mat_Product *product = mat->product;
737: Mat A = product->A, B = product->B;
739: PetscFunctionBegin;
740: PetscCall((*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat));
741: mat->ops->productnumeric = MatProductNumeric_ABt;
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
746: {
747: Mat_Product *product = mat->product;
748: Mat A = product->A, B = product->B, C = product->C;
750: PetscFunctionBegin;
751: PetscCall((*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat));
752: mat->ops->productnumeric = MatProductNumeric_ABC;
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: /*@
757: MatProductSymbolic - Perform the symbolic portion of a matrix-matrix product operation, this creates a data structure for use with the numerical
758: product to be done with `MatProductNumeric()`
760: Collective
762: Input/Output Parameter:
763: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
765: Level: intermediate
767: Note:
768: `MatProductSetFromOptions()` must have been called on `mat` before calling this function
770: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
771: @*/
772: PetscErrorCode MatProductSymbolic(Mat mat)
773: {
774: PetscLogEvent eventtype = -1;
775: PetscBool missing = PETSC_FALSE;
776: Mat_Product *product = mat->product;
777: Mat A = product->A;
778: Mat B = product->B;
779: Mat C = product->C;
781: PetscFunctionBegin;
783: MatCheckProduct(mat, 1);
784: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot run symbolic phase. Product data not empty");
785: switch (mat->product->type) {
786: case MATPRODUCT_AB:
787: eventtype = MAT_MatMultSymbolic;
788: break;
789: case MATPRODUCT_AtB:
790: eventtype = MAT_TransposeMatMultSymbolic;
791: break;
792: case MATPRODUCT_ABt:
793: eventtype = MAT_MatTransposeMultSymbolic;
794: break;
795: case MATPRODUCT_PtAP:
796: eventtype = MAT_PtAPSymbolic;
797: break;
798: case MATPRODUCT_RARt:
799: eventtype = MAT_RARtSymbolic;
800: break;
801: case MATPRODUCT_ABC:
802: eventtype = MAT_MatMatMultSymbolic;
803: break;
804: default:
805: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
806: }
807: mat->ops->productnumeric = NULL;
808: if (mat->ops->productsymbolic) {
809: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
810: PetscUseTypeMethod(mat, productsymbolic);
811: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
812: } else missing = PETSC_TRUE;
813: if (missing || !mat->product || !mat->ops->productnumeric) {
814: char errstr[256];
816: if (mat->product->type == MATPRODUCT_ABC) {
817: PetscCall(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));
818: } else {
819: PetscCall(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));
820: }
821: PetscCheck(mat->product->setfromoptionscalled, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() first", errstr);
822: PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unspecified symbolic phase for product %s. The product is not supported", errstr);
823: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after symbolic phase for product %s", errstr);
824: }
825: #if defined(PETSC_HAVE_DEVICE)
826: PetscBool bindingpropagates;
827: bindingpropagates = (PetscBool)((A->boundtocpu && A->bindingpropagates) || (B->boundtocpu && B->bindingpropagates));
828: if (C) bindingpropagates = (PetscBool)(bindingpropagates || (C->boundtocpu && C->bindingpropagates));
829: if (bindingpropagates) {
830: PetscCall(MatBindToCPU(mat, PETSC_TRUE));
831: PetscCall(MatSetBindingPropagates(mat, PETSC_TRUE));
832: }
833: #endif
834: /* set block sizes */
835: switch (product->type) {
836: case MATPRODUCT_PtAP:
837: if (B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->cmap->bs, B->cmap->bs));
838: break;
839: case MATPRODUCT_RARt:
840: if (B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->rmap->bs, B->rmap->bs));
841: break;
842: case MATPRODUCT_ABC:
843: PetscCall(MatSetBlockSizesFromMats(mat, A, C));
844: break;
845: case MATPRODUCT_AB:
846: PetscCall(MatSetBlockSizesFromMats(mat, A, B));
847: break;
848: case MATPRODUCT_AtB:
849: if (A->cmap->bs > 1 || B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, B->cmap->bs));
850: break;
851: case MATPRODUCT_ABt:
852: if (A->rmap->bs > 1 || B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, B->rmap->bs));
853: break;
854: default:
855: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
856: }
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: MatProductSetFill - Set an expected fill of the matrix whose values are computed via a matrix-matrix product operation
863: Collective
865: Input Parameters:
866: + mat - the matrix whose values are to be computed via a matrix-matrix product operation
867: - fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use `PETSC_DETERMINE` or `PETSC_CURRENT` if you do not have a good estimate.
868: If the product is a dense matrix, this value is not used.
870: Level: intermediate
872: Notes:
873: Use `fill` of `PETSC_DETERMINE` to use the default value.
875: The deprecated `PETSC_DEFAULT` is also supported to mean use the current value.
877: .seealso: [](ch_matrices), `MatProduct`, `PETSC_DETERMINE`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
878: @*/
879: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
880: {
881: PetscFunctionBegin;
883: MatCheckProduct(mat, 1);
884: if (fill == (PetscReal)PETSC_DETERMINE) mat->product->fill = mat->product->default_fill;
885: else if (fill != (PetscReal)PETSC_CURRENT) mat->product->fill = fill;
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: /*@
890: MatProductSetAlgorithm - Requests a particular algorithm for a matrix-matrix product operation that will perform to compute the given matrix
892: Collective
894: Input Parameters:
895: + mat - the matrix whose values are computed via a matrix-matrix product operation
896: - alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
898: Options Database Key:
899: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm`
901: Level: intermediate
903: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductClear()`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`, `MatProductGetAlgorithm()`
904: @*/
905: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
906: {
907: PetscFunctionBegin;
909: MatCheckProduct(mat, 1);
910: PetscCall(PetscFree(mat->product->alg));
911: PetscCall(PetscStrallocpy(alg, &mat->product->alg));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*@
916: MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation
918: Not Collective
920: Input Parameter:
921: . mat - the matrix whose values are computed via a matrix-matrix product operation
923: Output Parameter:
924: . alg - the selected algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
926: Level: intermediate
928: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`
929: @*/
930: PetscErrorCode MatProductGetAlgorithm(Mat mat, MatProductAlgorithm *alg)
931: {
932: PetscFunctionBegin;
934: PetscAssertPointer(alg, 2);
935: if (mat->product) *alg = mat->product->alg;
936: else *alg = NULL;
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: /*@
941: MatProductSetType - Sets a particular matrix-matrix product operation to be used to compute the values of the given matrix
943: Collective
945: Input Parameters:
946: + mat - the matrix whose values are computed via a matrix-matrix product operation
947: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`,
948: see `MatProductType`
950: Level: intermediate
952: Note:
953: The small t represents the transpose operation.
955: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductType`,
956: `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
957: @*/
958: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
959: {
960: PetscFunctionBegin;
962: MatCheckProduct(mat, 1);
964: if (productype != mat->product->type) {
965: if (mat->product->destroy) PetscCall((*mat->product->destroy)(&mat->product->data));
966: mat->product->destroy = NULL;
967: mat->product->data = NULL;
968: mat->ops->productsymbolic = NULL;
969: mat->ops->productnumeric = NULL;
970: }
971: mat->product->type = productype;
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: MatProductClear - Clears from the matrix any internal data structures related to the computation of the values of the matrix from matrix-matrix product operations
978: Collective
980: Input Parameter:
981: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
983: Options Database Key:
984: . -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
986: Level: intermediate
988: Notes:
989: This function should be called to remove any intermediate data used to compute the matrix to free up memory.
991: After having called this function, matrix-matrix product operations can no longer be used on `mat`
993: Developer Note:
994: This frees the `Mat_Product` context that was attached to the matrix during `MatProductCreate()` or `MatProductCreateWithMat()`
996: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`
997: @*/
998: PetscErrorCode MatProductClear(Mat mat)
999: {
1000: Mat_Product *product = mat->product;
1002: PetscFunctionBegin;
1004: if (product) {
1005: PetscCall(MatDestroy(&product->A));
1006: PetscCall(MatDestroy(&product->B));
1007: PetscCall(MatDestroy(&product->C));
1008: PetscCall(PetscFree(product->alg));
1009: PetscCall(MatDestroy(&product->Dwork));
1010: if (product->destroy) PetscCall((*product->destroy)(&product->data));
1011: }
1012: PetscCall(PetscFree(mat->product));
1013: mat->ops->productsymbolic = NULL;
1014: mat->ops->productnumeric = NULL;
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /* Create a supporting struct and attach it to the matrix product */
1019: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
1020: {
1021: Mat_Product *product = NULL;
1023: PetscFunctionBegin;
1025: PetscCheck(!D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product already present");
1026: PetscCall(PetscNew(&product));
1027: product->A = A;
1028: product->B = B;
1029: product->C = C;
1030: product->type = MATPRODUCT_UNSPECIFIED;
1031: product->Dwork = NULL;
1032: product->api_user = PETSC_FALSE;
1033: product->clear = PETSC_FALSE;
1034: product->setfromoptionscalled = PETSC_FALSE;
1035: PetscObjectParameterSetDefault(product, fill, 2);
1036: D->product = product;
1038: PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
1039: PetscCall(MatProductSetFill(D, PETSC_DEFAULT));
1041: PetscCall(PetscObjectReference((PetscObject)A));
1042: PetscCall(PetscObjectReference((PetscObject)B));
1043: PetscCall(PetscObjectReference((PetscObject)C));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*@
1048: MatProductCreateWithMat - Set a given matrix to have its values computed via matrix-matrix operations on other matrices.
1050: Collective
1052: Input Parameters:
1053: + A - the first matrix
1054: . B - the second matrix
1055: . C - the third matrix (optional, use `NULL` if not needed)
1056: - D - the matrix whose values are to be computed via a matrix-matrix product operation
1058: Level: intermediate
1060: Notes:
1061: Use `MatProductCreate()` if the matrix you wish computed `D` does not exist
1063: See `MatProductCreate()` for details on the usage of the matrix-matrix product operations
1065: Any product data currently attached to `D` will be freed
1067: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductType`, `MatProductSetType()`, `MatProductAlgorithm`,
1068: `MatProductSetAlgorithm`, `MatProductCreate()`, `MatProductClear()`
1069: @*/
1070: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
1071: {
1072: PetscFunctionBegin;
1075: MatCheckPreallocated(A, 1);
1076: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1077: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1081: MatCheckPreallocated(B, 2);
1082: PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1083: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1085: if (C) {
1088: MatCheckPreallocated(C, 3);
1089: PetscCheck(C->assembled, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1090: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1091: }
1095: MatCheckPreallocated(D, 4);
1096: PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1097: PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1099: /* Create a supporting struct and attach it to D */
1100: PetscCall(MatProductClear(D));
1101: PetscCall(MatProductCreate_Private(A, B, C, D));
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: /*@
1106: MatProductCreate - create a matrix to hold the result of a matrix-matrix (or matrix-matrix-matrix) product operation
1108: Collective
1110: Input Parameters:
1111: + A - the first matrix
1112: . B - the second matrix
1113: - C - the third matrix (or `NULL`)
1115: Output Parameter:
1116: . D - the matrix whose values are to be computed via a matrix-matrix product operation
1118: Level: intermediate
1120: Example:
1121: .vb
1122: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1123: MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1124: MatProductSetAlgorithm(D, alg)
1125: MatProductSetFill(D,fill)
1126: MatProductSetFromOptions(D)
1127: MatProductSymbolic(D)
1128: MatProductNumeric(D)
1129: Change numerical values in some of the matrices
1130: MatProductNumeric(D)
1131: .ve
1133: Notes:
1134: Use `MatProductCreateWithMat()` if `D` the matrix you wish computed already exists.
1136: The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure of the input matrices.
1138: Developer Notes:
1139: It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1140: Is there error checking for it?
1142: On this call, auxiliary data needed to compute the product is stored in `D` in a `Mat_Product` context. A call to `MatProductClear()` frees this
1143: information.
1145: Each `MatProductAlgorithm` associated with a particular `MatType` stores additional data needed for the product computation
1146: (generally this data is computed in `MatProductSymbolic()`) inside the `Mat_Product` context in a `MatProductCtx_XXX` data structure
1147: and provides a `MatProductCtxDestroy_XXX()` routine to free that data. The `MatProductAlgorithm` and `MatType` specific destroy routine is called by
1148: `MatProductClear()`.
1150: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`,
1151: `MatProductSymbolic()`, `MatProductNumeric()`, `MatProductAlgorithm`, `MatProductType`
1152: @*/
1153: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1154: {
1155: PetscFunctionBegin;
1160: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix A");
1161: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix B");
1163: if (C) {
1166: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix C");
1167: }
1169: PetscAssertPointer(D, 4);
1170: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), D));
1171: /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1172: PetscCall(MatProductCreate_Private(A, B, C, *D));
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: /*
1177: These are safe basic implementations of ABC, RARt and PtAP
1178: that do not rely on mat->ops->matmatop function pointers.
1179: They only use the MatProduct API and are currently used by
1180: cuSPARSE and KOKKOS-KERNELS backends
1181: */
1182: typedef struct {
1183: Mat BC;
1184: Mat ABC;
1185: } MatProductCtx_MatMatMatPrivate;
1187: static PetscErrorCode MatProductCtxDestroy_MatMatMatPrivate(void **data)
1188: {
1189: MatProductCtx_MatMatMatPrivate *mmdata = *(MatProductCtx_MatMatMatPrivate **)data;
1191: PetscFunctionBegin;
1192: PetscCall(MatDestroy(&mmdata->BC));
1193: PetscCall(MatDestroy(&mmdata->ABC));
1194: PetscCall(PetscFree(*data));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1199: {
1200: Mat_Product *product = mat->product;
1201: MatProductCtx_MatMatMatPrivate *mmabc;
1203: PetscFunctionBegin;
1204: MatCheckProduct(mat, 1);
1205: PetscCheck(mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data empty");
1206: mmabc = (MatProductCtx_MatMatMatPrivate *)mat->product->data;
1207: PetscCheck(mmabc->BC->ops->productnumeric, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing numeric stage");
1208: /* use function pointer directly to prevent logging */
1209: PetscCall((*mmabc->BC->ops->productnumeric)(mmabc->BC));
1210: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1211: mat->product = mmabc->ABC->product;
1212: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1213: /* use function pointer directly to prevent logging */
1214: PetscUseTypeMethod(mat, productnumeric);
1215: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1216: mat->product = product;
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1221: {
1222: Mat_Product *product = mat->product;
1223: Mat A, B, C;
1224: MatProductType p1, p2;
1225: MatProductCtx_MatMatMatPrivate *mmabc;
1226: const char *prefix;
1228: PetscFunctionBegin;
1229: MatCheckProduct(mat, 1);
1230: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data not empty");
1231: PetscCall(MatGetOptionsPrefix(mat, &prefix));
1232: PetscCall(PetscNew(&mmabc));
1233: product->data = mmabc;
1234: product->destroy = MatProductCtxDestroy_MatMatMatPrivate;
1235: switch (product->type) {
1236: case MATPRODUCT_PtAP:
1237: p1 = MATPRODUCT_AB;
1238: p2 = MATPRODUCT_AtB;
1239: A = product->B;
1240: B = product->A;
1241: C = product->B;
1242: if (A->cmap->bs > 0 && C->cmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, C->cmap->bs));
1243: break;
1244: case MATPRODUCT_RARt:
1245: p1 = MATPRODUCT_ABt;
1246: p2 = MATPRODUCT_AB;
1247: A = product->B;
1248: B = product->A;
1249: C = product->B;
1250: if (A->rmap->bs > 0 && C->rmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, C->rmap->bs));
1251: break;
1252: case MATPRODUCT_ABC:
1253: p1 = MATPRODUCT_AB;
1254: p2 = MATPRODUCT_AB;
1255: A = product->A;
1256: B = product->B;
1257: C = product->C;
1258: PetscCall(MatSetBlockSizesFromMats(mat, A, C));
1259: break;
1260: default:
1261: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1262: }
1263: PetscCall(MatProductCreate(B, C, NULL, &mmabc->BC));
1264: PetscCall(MatSetOptionsPrefix(mmabc->BC, prefix));
1265: PetscCall(MatAppendOptionsPrefix(mmabc->BC, "P1_"));
1266: PetscCall(MatProductSetType(mmabc->BC, p1));
1267: PetscCall(MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT));
1268: PetscCall(MatProductSetFill(mmabc->BC, product->fill));
1269: mmabc->BC->product->api_user = product->api_user;
1270: PetscCall(MatProductSetFromOptions(mmabc->BC));
1271: PetscCheck(mmabc->BC->ops->productsymbolic, 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);
1272: /* use function pointer directly to prevent logging */
1273: PetscCall((*mmabc->BC->ops->productsymbolic)(mmabc->BC));
1275: PetscCall(MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC));
1276: PetscCall(MatSetOptionsPrefix(mmabc->ABC, prefix));
1277: PetscCall(MatAppendOptionsPrefix(mmabc->ABC, "P2_"));
1278: PetscCall(MatProductSetType(mmabc->ABC, p2));
1279: PetscCall(MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT));
1280: PetscCall(MatProductSetFill(mmabc->ABC, product->fill));
1281: mmabc->ABC->product->api_user = product->api_user;
1282: PetscCall(MatProductSetFromOptions(mmabc->ABC));
1283: PetscCheck(mmabc->ABC->ops->productsymbolic, 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);
1284: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1285: mat->product = mmabc->ABC->product;
1286: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1287: /* use function pointer directly to prevent logging */
1288: PetscUseTypeMethod(mat, productsymbolic);
1289: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1290: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1291: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1292: mat->product = product;
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: /*@
1297: MatProductGetType - Returns the type of matrix-matrix product associated with computing values for the given matrix
1299: Not Collective
1301: Input Parameter:
1302: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1304: Output Parameter:
1305: . mtype - the `MatProductType`
1307: Level: intermediate
1309: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1310: @*/
1311: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1312: {
1313: PetscFunctionBegin;
1315: PetscAssertPointer(mtype, 2);
1316: *mtype = MATPRODUCT_UNSPECIFIED;
1317: if (mat->product) *mtype = mat->product->type;
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: /*@
1322: MatProductGetMats - Returns the matrices associated with the matrix-matrix product associated with computing values for the given matrix
1324: Not Collective
1326: Input Parameter:
1327: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1329: Output Parameters:
1330: + A - the first matrix
1331: . B - the second matrix
1332: - C - the third matrix (may be `NULL` for some `MatProductType`)
1334: Level: intermediate
1336: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1337: @*/
1338: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1339: {
1340: PetscFunctionBegin;
1342: if (A) *A = mat->product ? mat->product->A : NULL;
1343: if (B) *B = mat->product ? mat->product->B : NULL;
1344: if (C) *C = mat->product ? mat->product->C : NULL;
1345: PetscFunctionReturn(PETSC_SUCCESS);
1346: }