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: PetscCall((*C->ops->transposematmultnumeric)(P, AP, C));
56: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
66: PetscCall(PetscInfo(C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
67: /* AP = A*P */
68: PetscCall(MatProductCreate(A, P, NULL, &AP));
69: PetscCall(MatProductSetType(AP, MATPRODUCT_AB));
70: PetscCall(MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT));
71: PetscCall(MatProductSetFill(AP, fill));
72: PetscCall(MatProductSetFromOptions(AP));
73: PetscCall(MatProductSymbolic(AP));
75: /* C = P^T*AP */
76: PetscCall(MatProductSetType(C, MATPRODUCT_AtB));
77: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
78: product->A = P;
79: product->B = AP;
80: PetscCall(MatProductSetFromOptions(C));
81: PetscCall(MatProductSymbolic(C));
83: /* resume user's original input matrix setting for A and B */
84: product->A = A;
85: product->B = P;
86: product->Dwork = AP;
88: C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
93: {
94: Mat_Product *product = C->product;
95: Mat R = product->B, RA = product->Dwork;
97: PetscFunctionBegin;
98: /* RA = R*A */
99: PetscCall(MatProductNumeric(RA));
100: /* C = RA*R^T */
101: PetscCall((*C->ops->mattransposemultnumeric)(RA, R, C));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
106: {
107: Mat_Product *product = C->product;
108: Mat A = product->A, R = product->B, RA;
109: PetscReal fill = product->fill;
111: PetscFunctionBegin;
112: PetscCall(PetscInfo(C, "for A %s, R %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
113: /* RA = R*A */
114: PetscCall(MatProductCreate(R, A, NULL, &RA));
115: PetscCall(MatProductSetType(RA, MATPRODUCT_AB));
116: PetscCall(MatProductSetAlgorithm(RA, MATPRODUCTALGORITHMDEFAULT));
117: PetscCall(MatProductSetFill(RA, fill));
118: PetscCall(MatProductSetFromOptions(RA));
119: PetscCall(MatProductSymbolic(RA));
121: /* C = RA*R^T */
122: PetscCall(MatProductSetType(C, MATPRODUCT_ABt));
123: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
124: product->A = RA;
125: PetscCall(MatProductSetFromOptions(C));
126: PetscCall(MatProductSymbolic(C));
128: /* resume user's original input matrix setting for A */
129: product->A = A;
130: product->Dwork = RA; /* save here so it will be destroyed with product C */
131: C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
136: {
137: Mat_Product *product = mat->product;
138: Mat A = product->A, BC = product->Dwork;
140: PetscFunctionBegin;
141: /* Numeric BC = B*C */
142: PetscCall(MatProductNumeric(BC));
143: /* Numeric mat = A*BC */
144: PetscCall((*mat->ops->matmultnumeric)(A, BC, mat));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
149: {
150: Mat_Product *product = mat->product;
151: Mat B = product->B, C = product->C, BC;
152: PetscReal fill = product->fill;
154: PetscFunctionBegin;
155: 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));
156: /* Symbolic BC = B*C */
157: PetscCall(MatProductCreate(B, C, NULL, &BC));
158: PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
159: PetscCall(MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT));
160: PetscCall(MatProductSetFill(BC, fill));
161: PetscCall(MatProductSetFromOptions(BC));
162: PetscCall(MatProductSymbolic(BC));
164: /* Symbolic mat = A*BC */
165: PetscCall(MatProductSetType(mat, MATPRODUCT_AB));
166: PetscCall(MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT));
167: product->B = BC;
168: product->Dwork = BC;
169: PetscCall(MatProductSetFromOptions(mat));
170: PetscCall(MatProductSymbolic(mat));
172: /* resume user's original input matrix setting for B */
173: product->B = B;
174: mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
179: {
180: Mat_Product *product = mat->product;
182: PetscFunctionBegin;
183: switch (product->type) {
184: case MATPRODUCT_PtAP:
185: PetscCall(MatProductSymbolic_PtAP_Unsafe(mat));
186: break;
187: case MATPRODUCT_RARt:
188: PetscCall(MatProductSymbolic_RARt_Unsafe(mat));
189: break;
190: case MATPRODUCT_ABC:
191: PetscCall(MatProductSymbolic_ABC_Unsafe(mat));
192: break;
193: default:
194: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: MatProductReplaceMats - Replace the input matrices for the matrix-matrix product operation inside the computed matrix
202: Collective
204: Input Parameters:
205: + A - the matrix or `NULL` if not being replaced
206: . B - the matrix or `NULL` if not being replaced
207: . C - the matrix or `NULL` if not being replaced
208: - D - the matrix whose values are computed via a matrix-matrix product operation
210: Level: intermediate
212: Note:
213: To reuse the symbolic phase, the input matrices must have exactly the same data structure as the replaced one.
214: If the type of any of the input matrices is different than what was previously used, or their symmetry flag changed but
215: the symbolic phase took advantage of their symmetry, the product is cleared and `MatProductSetFromOptions()`
216: and `MatProductSymbolic()` are invoked again.
218: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic()`, `MatProductClear()`
219: @*/
220: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
221: {
222: Mat_Product *product;
223: PetscBool flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;
225: PetscFunctionBegin;
227: MatCheckProduct(D, 4);
228: product = D->product;
229: if (A) {
231: PetscCall(PetscObjectReference((PetscObject)A));
232: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA));
233: PetscCall(MatIsSymmetricKnown(A, &isset, &issym));
234: 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 */
235: flgA = PETSC_FALSE;
236: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
237: }
238: PetscCall(MatDestroy(&product->A));
239: product->A = A;
240: }
241: if (B) {
243: PetscCall(PetscObjectReference((PetscObject)B));
244: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB));
245: PetscCall(MatIsSymmetricKnown(B, &isset, &issym));
246: if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
247: flgB = PETSC_FALSE;
248: product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
249: }
250: PetscCall(MatDestroy(&product->B));
251: product->B = B;
252: }
253: if (C) {
255: PetscCall(PetscObjectReference((PetscObject)C));
256: PetscCall(PetscObjectTypeCompare((PetscObject)product->C, ((PetscObject)C)->type_name, &flgC));
257: PetscCall(MatIsSymmetricKnown(C, &isset, &issym));
258: if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
259: flgC = PETSC_FALSE;
260: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
261: }
262: PetscCall(MatDestroy(&product->C));
263: product->C = C;
264: }
265: /* Any of the replaced mats is of a different type, reset */
266: if (!flgA || !flgB || !flgC) {
267: if (D->product->destroy) PetscCall((*D->product->destroy)(D->product->data));
268: D->product->destroy = NULL;
269: D->product->data = NULL;
270: if (D->ops->productnumeric || D->ops->productsymbolic) {
271: PetscCall(MatProductSetFromOptions(D));
272: PetscCall(MatProductSymbolic(D));
273: }
274: }
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
279: {
280: Mat_Product *product = C->product;
281: Mat A = product->A, B = product->B;
282: PetscInt k, K = B->cmap->N;
283: PetscBool t = PETSC_TRUE, iscuda = PETSC_FALSE;
284: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
285: char *Btype = NULL, *Ctype = NULL;
287: PetscFunctionBegin;
288: switch (product->type) {
289: case MATPRODUCT_AB:
290: t = PETSC_FALSE;
291: case MATPRODUCT_AtB:
292: break;
293: default:
294: 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);
295: }
296: if (PetscDefined(HAVE_CUDA)) {
297: VecType vtype;
299: PetscCall(MatGetVecType(A, &vtype));
300: PetscCall(PetscStrcmp(vtype, VECCUDA, &iscuda));
301: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECSEQCUDA, &iscuda));
302: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECMPICUDA, &iscuda));
303: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
304: PetscCall(PetscStrallocpy(((PetscObject)B)->type_name, &Btype));
305: PetscCall(PetscStrallocpy(((PetscObject)C)->type_name, &Ctype));
306: PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
307: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
308: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
309: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
310: }
311: PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
312: } else { /* Make sure we have up-to-date data on the CPU */
313: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
314: Bcpu = B->boundtocpu;
315: Ccpu = C->boundtocpu;
316: #endif
317: PetscCall(MatBindToCPU(B, PETSC_TRUE));
318: PetscCall(MatBindToCPU(C, PETSC_TRUE));
319: }
320: }
321: for (k = 0; k < K; k++) {
322: Vec x, y;
324: PetscCall(MatDenseGetColumnVecRead(B, k, &x));
325: PetscCall(MatDenseGetColumnVecWrite(C, k, &y));
326: if (t) {
327: PetscCall(MatMultTranspose(A, x, y));
328: } else {
329: PetscCall(MatMult(A, x, y));
330: }
331: PetscCall(MatDenseRestoreColumnVecRead(B, k, &x));
332: PetscCall(MatDenseRestoreColumnVecWrite(C, k, &y));
333: }
334: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
335: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
336: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
337: if (PetscDefined(HAVE_CUDA)) {
338: if (iscuda) {
339: PetscCall(MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B));
340: PetscCall(MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C));
341: } else {
342: PetscCall(MatBindToCPU(B, Bcpu));
343: PetscCall(MatBindToCPU(C, Ccpu));
344: }
345: }
346: PetscCall(PetscFree(Btype));
347: PetscCall(PetscFree(Ctype));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
352: {
353: Mat_Product *product = C->product;
354: Mat A = product->A, B = product->B;
355: PetscBool isdense;
357: PetscFunctionBegin;
358: switch (product->type) {
359: case MATPRODUCT_AB:
360: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
361: break;
362: case MATPRODUCT_AtB:
363: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
364: break;
365: default:
366: 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);
367: }
368: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
369: if (!isdense) {
370: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
371: /* If matrix type of C was not set or not dense, we need to reset the pointer */
372: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
373: }
374: C->ops->productnumeric = MatProductNumeric_X_Dense;
375: PetscCall(MatSetUp(C));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /* a single driver to query the dispatching */
380: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
381: {
382: Mat_Product *product = mat->product;
383: PetscInt Am, An, Bm, Bn, Cm, Cn;
384: Mat A = product->A, B = product->B, C = product->C;
385: const char *const Bnames[] = {"B", "R", "P"};
386: const char *bname;
387: PetscErrorCode (*fA)(Mat);
388: PetscErrorCode (*fB)(Mat);
389: PetscErrorCode (*fC)(Mat);
390: PetscErrorCode (*f)(Mat) = NULL;
392: PetscFunctionBegin;
393: mat->ops->productsymbolic = NULL;
394: mat->ops->productnumeric = NULL;
395: if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(PETSC_SUCCESS);
396: PetscCheck(A, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing A mat");
397: PetscCheck(B, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing B mat");
398: PetscCheck(product->type != MATPRODUCT_ABC || C, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing C mat");
399: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
400: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
401: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
402: else bname = Bnames[0];
404: /* Check matrices sizes */
405: Am = A->rmap->N;
406: An = A->cmap->N;
407: Bm = B->rmap->N;
408: Bn = B->cmap->N;
409: Cm = C ? C->rmap->N : 0;
410: Cn = C ? C->cmap->N : 0;
411: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
412: PetscInt t = Bn;
413: Bn = Bm;
414: Bm = t;
415: }
416: if (product->type == MATPRODUCT_AtB) An = Am;
418: 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,
419: MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
420: 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,
421: MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);
423: fA = A->ops->productsetfromoptions;
424: fB = B->ops->productsetfromoptions;
425: fC = C ? C->ops->productsetfromoptions : fA;
426: if (C) {
427: 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));
428: } else {
429: PetscCall(PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name));
430: }
431: if (fA == fB && fA == fC && fA) {
432: PetscCall(PetscInfo(mat, " matching op\n"));
433: PetscCall((*fA)(mat));
434: }
435: /* We may have found f but it did not succeed */
436: if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
437: char mtypes[256];
438: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes)));
439: PetscCall(PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes)));
440: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
441: PetscCall(PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes)));
442: if (C) {
443: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
444: PetscCall(PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes)));
445: }
446: PetscCall(PetscStrlcat(mtypes, "_C", sizeof(mtypes)));
447: #if defined(__clang__)
448: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
449: #elif defined(__GNUC__) || defined(__GNUG__)
450: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
451: #endif
452: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
453: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
454: if (!f) {
455: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
456: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
457: }
458: if (!f && C) {
459: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
460: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
461: }
462: if (f) PetscCall((*f)(mat));
464: /* We may have found f but it did not succeed */
465: /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
466: if (!mat->ops->productsymbolic) {
467: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes)));
468: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
469: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
470: if (!f) {
471: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
472: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
473: }
474: if (!f && C) {
475: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
476: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
477: }
478: }
479: if (f) PetscCall((*f)(mat));
480: }
481: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
482: /* We may have found f but it did not succeed */
483: if (!mat->ops->productsymbolic) {
484: /* we can still compute the product if B is of type dense */
485: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
486: PetscBool isdense;
488: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
489: if (isdense) {
490: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
491: PetscCall(PetscInfo(mat, " using basic looping over columns of a dense matrix\n"));
492: }
493: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
494: /*
495: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
496: the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
497: before computing the symbolic phase
498: */
499: PetscCall(PetscInfo(mat, " symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n"));
500: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
501: }
502: }
503: if (!mat->ops->productsymbolic) PetscCall(PetscInfo(mat, " symbolic product is not supported\n"));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@
508: MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product operation where the type,
509: the algorithm etc are determined from the options database.
511: Logically Collective
513: Input Parameter:
514: . mat - the matrix whose values are computed via a matrix-matrix product operation
516: Options Database Keys:
517: + -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
518: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm` for possible values
519: - -mat_product_algorithm_backend_cpu - Use the CPU to perform the computation even if the matrix is a GPU matrix
521: Level: intermediate
523: Note:
524: The `-mat_product_clear` option reduces memory usage but means that the matrix cannot be re-used for a matrix-matrix product operation
526: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`,
527: `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductAlgorithm`
528: @*/
529: PetscErrorCode MatProductSetFromOptions(Mat mat)
530: {
531: PetscFunctionBegin;
533: MatCheckProduct(mat, 1);
534: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot call MatProductSetFromOptions() with already present data");
535: mat->product->setfromoptionscalled = PETSC_TRUE;
536: PetscObjectOptionsBegin((PetscObject)mat);
537: PetscCall(PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL));
538: PetscCall(PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()"));
539: PetscOptionsEnd();
540: PetscCall(MatProductSetFromOptions_Private(mat));
541: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing product after setup phase");
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: /*@
546: MatProductView - View the private matrix-matrix algorithm object within a matrix
548: Logically Collective
550: Input Parameters:
551: + mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
552: - viewer - where the information on the matrix-matrix algorithm of `mat` should be reviewed
554: Level: intermediate
556: Developer Note:
557: Shouldn't this information be printed from an appropriate `MatView()` with perhaps certain formats set?
559: .seealso: [](ch_matrices), `MatProductType`, `Mat`, `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
560: @*/
561: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
562: {
563: PetscFunctionBegin;
565: if (!mat->product) PetscFunctionReturn(PETSC_SUCCESS);
566: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer));
568: PetscCheckSameComm(mat, 1, viewer, 2);
569: if (mat->product->view) PetscCall((*mat->product->view)(mat, viewer));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /* these are basic implementations relying on the old function pointers
574: * they are dangerous and should be removed in the future */
575: PetscErrorCode MatProductNumeric_AB(Mat mat)
576: {
577: Mat_Product *product = mat->product;
578: Mat A = product->A, B = product->B;
580: PetscFunctionBegin;
581: PetscCall((*mat->ops->matmultnumeric)(A, B, mat));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: PetscErrorCode MatProductNumeric_AtB(Mat mat)
586: {
587: Mat_Product *product = mat->product;
588: Mat A = product->A, B = product->B;
590: PetscFunctionBegin;
591: PetscCall((*mat->ops->transposematmultnumeric)(A, B, mat));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: PetscErrorCode MatProductNumeric_ABt(Mat mat)
596: {
597: Mat_Product *product = mat->product;
598: Mat A = product->A, B = product->B;
600: PetscFunctionBegin;
601: PetscCall((*mat->ops->mattransposemultnumeric)(A, B, mat));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
606: {
607: Mat_Product *product = mat->product;
608: Mat A = product->A, B = product->B;
610: PetscFunctionBegin;
611: PetscCall((*mat->ops->ptapnumeric)(A, B, mat));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: PetscErrorCode MatProductNumeric_RARt(Mat mat)
616: {
617: Mat_Product *product = mat->product;
618: Mat A = product->A, B = product->B;
620: PetscFunctionBegin;
621: PetscCall((*mat->ops->rartnumeric)(A, B, mat));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: PetscErrorCode MatProductNumeric_ABC(Mat mat)
626: {
627: Mat_Product *product = mat->product;
628: Mat A = product->A, B = product->B, C = product->C;
630: PetscFunctionBegin;
631: PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*@
636: MatProductNumeric - Compute a matrix-matrix product operation with the numerical values
638: Collective
640: Input/Output Parameter:
641: . mat - the matrix whose values are computed via a matrix-matrix product operation
643: Level: intermediate
645: Note:
646: `MatProductSymbolic()` must have been called on `mat` before calling this function
648: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
649: @*/
650: PetscErrorCode MatProductNumeric(Mat mat)
651: {
652: PetscLogEvent eventtype = -1;
654: PetscFunctionBegin;
656: MatCheckProduct(mat, 1);
657: switch (mat->product->type) {
658: case MATPRODUCT_AB:
659: eventtype = MAT_MatMultNumeric;
660: break;
661: case MATPRODUCT_AtB:
662: eventtype = MAT_TransposeMatMultNumeric;
663: break;
664: case MATPRODUCT_ABt:
665: eventtype = MAT_MatTransposeMultNumeric;
666: break;
667: case MATPRODUCT_PtAP:
668: eventtype = MAT_PtAPNumeric;
669: break;
670: case MATPRODUCT_RARt:
671: eventtype = MAT_RARtNumeric;
672: break;
673: case MATPRODUCT_ABC:
674: eventtype = MAT_MatMatMultNumeric;
675: break;
676: default:
677: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
678: }
680: if (mat->ops->productnumeric) {
681: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
682: PetscUseTypeMethod(mat, productnumeric);
683: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
684: } else if (mat->product) {
685: char errstr[256];
687: if (mat->product->type == MATPRODUCT_ABC) {
688: 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));
689: } else {
690: 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));
691: }
692: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified numeric phase for product %s", errstr);
693: }
694: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after numeric phase for product");
696: if (mat->product->clear) PetscCall(MatProductClear(mat));
697: PetscCall(PetscObjectStateIncrease((PetscObject)mat));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: /* these are basic implementations relying on the old function pointers
702: * they are dangerous and should be removed in the future */
703: PetscErrorCode MatProductSymbolic_AB(Mat mat)
704: {
705: Mat_Product *product = mat->product;
706: Mat A = product->A, B = product->B;
708: PetscFunctionBegin;
709: PetscCall((*mat->ops->matmultsymbolic)(A, B, product->fill, mat));
710: mat->ops->productnumeric = MatProductNumeric_AB;
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
715: {
716: Mat_Product *product = mat->product;
717: Mat A = product->A, B = product->B;
719: PetscFunctionBegin;
720: PetscCall((*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat));
721: mat->ops->productnumeric = MatProductNumeric_AtB;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
726: {
727: Mat_Product *product = mat->product;
728: Mat A = product->A, B = product->B;
730: PetscFunctionBegin;
731: PetscCall((*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat));
732: mat->ops->productnumeric = MatProductNumeric_ABt;
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
737: {
738: Mat_Product *product = mat->product;
739: Mat A = product->A, B = product->B, C = product->C;
741: PetscFunctionBegin;
742: PetscCall((*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat));
743: mat->ops->productnumeric = MatProductNumeric_ABC;
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /*@
748: MatProductSymbolic - Perform the symbolic portion of a matrix-matrix product operation, this creates a data structure for use with the numerical
749: product to be done with `MatProductNumeric()`
751: Collective
753: Input/Output Parameter:
754: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
756: Level: intermediate
758: Note:
759: `MatProductSetFromOptions()` must have been called on `mat` before calling this function
761: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
762: @*/
763: PetscErrorCode MatProductSymbolic(Mat mat)
764: {
765: PetscLogEvent eventtype = -1;
766: PetscBool missing = PETSC_FALSE;
768: PetscFunctionBegin;
770: MatCheckProduct(mat, 1);
771: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot run symbolic phase. Product data not empty");
772: switch (mat->product->type) {
773: case MATPRODUCT_AB:
774: eventtype = MAT_MatMultSymbolic;
775: break;
776: case MATPRODUCT_AtB:
777: eventtype = MAT_TransposeMatMultSymbolic;
778: break;
779: case MATPRODUCT_ABt:
780: eventtype = MAT_MatTransposeMultSymbolic;
781: break;
782: case MATPRODUCT_PtAP:
783: eventtype = MAT_PtAPSymbolic;
784: break;
785: case MATPRODUCT_RARt:
786: eventtype = MAT_RARtSymbolic;
787: break;
788: case MATPRODUCT_ABC:
789: eventtype = MAT_MatMatMultSymbolic;
790: break;
791: default:
792: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
793: }
794: mat->ops->productnumeric = NULL;
795: if (mat->ops->productsymbolic) {
796: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
797: PetscUseTypeMethod(mat, productsymbolic);
798: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
799: } else missing = PETSC_TRUE;
801: if (missing || !mat->product || !mat->ops->productnumeric) {
802: char errstr[256];
804: if (mat->product->type == MATPRODUCT_ABC) {
805: 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));
806: } else {
807: 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));
808: }
809: PetscCheck(mat->product->setfromoptionscalled, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() first", errstr);
810: PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unspecified symbolic phase for product %s. The product is not supported", errstr);
811: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after symbolic phase for product %s", errstr);
812: }
814: #if defined(PETSC_HAVE_DEVICE)
815: Mat A = mat->product->A;
816: Mat B = mat->product->B;
817: Mat C = mat->product->C;
818: PetscBool bindingpropagates;
819: bindingpropagates = (PetscBool)((A->boundtocpu && A->bindingpropagates) || (B->boundtocpu && B->bindingpropagates));
820: if (C) bindingpropagates = (PetscBool)(bindingpropagates || (C->boundtocpu && C->bindingpropagates));
821: if (bindingpropagates) {
822: PetscCall(MatBindToCPU(mat, PETSC_TRUE));
823: PetscCall(MatSetBindingPropagates(mat, PETSC_TRUE));
824: }
825: #endif
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: /*@
830: MatProductSetFill - Set an expected fill of the matrix whose values are computed via a matrix-matrix product operation
832: Collective
834: Input Parameters:
835: + mat - the matrix whose values are to be computed via a matrix-matrix product operation
836: - 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.
837: If the product is a dense matrix, this value is not used.
839: Level: intermediate
841: Notes:
842: Use `fill` of `PETSC_DETERMINE` to use the default value.
844: The deprecated `PETSC_DEFAULT` is also supported to mean use the current value.
846: .seealso: [](ch_matrices), `MatProduct`, `PETSC_DETERMINE`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
847: @*/
848: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
849: {
850: PetscFunctionBegin;
852: MatCheckProduct(mat, 1);
853: if (fill == (PetscReal)PETSC_DETERMINE) mat->product->fill = mat->product->default_fill;
854: else if (fill != (PetscReal)PETSC_CURRENT) mat->product->fill = fill;
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*@
859: MatProductSetAlgorithm - Requests a particular algorithm for a matrix-matrix product operation that will perform to compute the given matrix
861: Collective
863: Input Parameters:
864: + mat - the matrix whose values are computed via a matrix-matrix product operation
865: - alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
867: Options Database Key:
868: . -mat_product_algorithm <algorithm> - Sets the algorithm, see `MatProductAlgorithm`
870: Level: intermediate
872: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductClear()`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`, `MatProductGetAlgorithm()`
873: @*/
874: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
875: {
876: PetscFunctionBegin;
878: MatCheckProduct(mat, 1);
879: PetscCall(PetscFree(mat->product->alg));
880: PetscCall(PetscStrallocpy(alg, &mat->product->alg));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*@
885: MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation
887: Not Collective
889: Input Parameter:
890: . mat - the matrix whose values are computed via a matrix-matrix product operation
892: Output Parameter:
893: . alg - the selected algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
895: Level: intermediate
897: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`
898: @*/
899: PetscErrorCode MatProductGetAlgorithm(Mat mat, MatProductAlgorithm *alg)
900: {
901: PetscFunctionBegin;
903: PetscAssertPointer(alg, 2);
904: if (mat->product) *alg = mat->product->alg;
905: else *alg = NULL;
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }
909: /*@
910: MatProductSetType - Sets a particular matrix-matrix product operation to be used to compute the values of the given matrix
912: Collective
914: Input Parameters:
915: + mat - the matrix whose values are computed via a matrix-matrix product operation
916: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`,
917: see `MatProductType`
919: Level: intermediate
921: Note:
922: The small t represents the transpose operation.
924: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductType`,
925: `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
926: @*/
927: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
928: {
929: PetscFunctionBegin;
931: MatCheckProduct(mat, 1);
933: if (productype != mat->product->type) {
934: if (mat->product->destroy) PetscCall((*mat->product->destroy)(mat->product->data));
935: mat->product->destroy = NULL;
936: mat->product->data = NULL;
937: mat->ops->productsymbolic = NULL;
938: mat->ops->productnumeric = NULL;
939: }
940: mat->product->type = productype;
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@
945: MatProductClear - Clears from the matrix any internal data structures related to the computation of the values of the matrix from matrix-matrix product operations
947: Collective
949: Input Parameter:
950: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
952: Options Database Key:
953: . -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
955: Level: intermediate
957: Notes:
958: This function should be called to remove any intermediate data used to compute the matrix to free up memory.
960: After having called this function, matrix-matrix product operations can no longer be used on `mat`
962: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`
963: @*/
964: PetscErrorCode MatProductClear(Mat mat)
965: {
966: Mat_Product *product = mat->product;
968: PetscFunctionBegin;
970: if (product) {
971: PetscCall(MatDestroy(&product->A));
972: PetscCall(MatDestroy(&product->B));
973: PetscCall(MatDestroy(&product->C));
974: PetscCall(PetscFree(product->alg));
975: PetscCall(MatDestroy(&product->Dwork));
976: if (product->destroy) PetscCall((*product->destroy)(product->data));
977: }
978: PetscCall(PetscFree(mat->product));
979: mat->ops->productsymbolic = NULL;
980: mat->ops->productnumeric = NULL;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /* Create a supporting struct and attach it to the matrix product */
985: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
986: {
987: Mat_Product *product = NULL;
989: PetscFunctionBegin;
991: PetscCheck(!D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product already present");
992: PetscCall(PetscNew(&product));
993: product->A = A;
994: product->B = B;
995: product->C = C;
996: product->type = MATPRODUCT_UNSPECIFIED;
997: product->Dwork = NULL;
998: product->api_user = PETSC_FALSE;
999: product->clear = PETSC_FALSE;
1000: product->setfromoptionscalled = PETSC_FALSE;
1001: PetscObjectParameterSetDefault(product, fill, 2);
1002: D->product = product;
1004: PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
1005: PetscCall(MatProductSetFill(D, PETSC_DEFAULT));
1007: PetscCall(PetscObjectReference((PetscObject)A));
1008: PetscCall(PetscObjectReference((PetscObject)B));
1009: PetscCall(PetscObjectReference((PetscObject)C));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*@
1014: MatProductCreateWithMat - Set a given matrix to have its values computed via matrix-matrix operations on other matrices.
1016: Collective
1018: Input Parameters:
1019: + A - the first matrix
1020: . B - the second matrix
1021: . C - the third matrix (optional, use `NULL` if not needed)
1022: - D - the matrix whose values are to be computed via a matrix-matrix product operation
1024: Level: intermediate
1026: Notes:
1027: Use `MatProductCreate()` if the matrix you wish computed (the `D` matrix) does not already exist
1029: See `MatProductCreate()` for details on the usage of the matrix-matrix product operations
1031: Any product data currently attached to `D` will be cleared
1033: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductType`, `MatProductSetType()`, `MatProductAlgorithm`,
1034: `MatProductSetAlgorithm`, `MatProductCreate()`, `MatProductClear()`
1035: @*/
1036: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
1037: {
1038: PetscFunctionBegin;
1041: MatCheckPreallocated(A, 1);
1042: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1043: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1047: MatCheckPreallocated(B, 2);
1048: PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1049: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1051: if (C) {
1054: MatCheckPreallocated(C, 3);
1055: PetscCheck(C->assembled, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1056: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1057: }
1061: MatCheckPreallocated(D, 4);
1062: PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1063: PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1065: /* Create a supporting struct and attach it to D */
1066: PetscCall(MatProductClear(D));
1067: PetscCall(MatProductCreate_Private(A, B, C, D));
1068: PetscFunctionReturn(PETSC_SUCCESS);
1069: }
1071: /*@
1072: MatProductCreate - create a matrix to hold the result of a matrix-matrix product operation
1074: Collective
1076: Input Parameters:
1077: + A - the first matrix
1078: . B - the second matrix
1079: - C - the third matrix (or `NULL`)
1081: Output Parameter:
1082: . D - the matrix whose values are to be computed via a matrix-matrix product operation
1084: Level: intermediate
1086: Example:
1087: .vb
1088: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1089: MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1090: MatProductSetAlgorithm(D, alg)
1091: MatProductSetFill(D,fill)
1092: MatProductSetFromOptions(D)
1093: MatProductSymbolic(D)
1094: MatProductNumeric(D)
1095: Change numerical values in some of the matrices
1096: MatProductNumeric(D)
1097: .ve
1099: Notes:
1100: Use `MatProductCreateWithMat()` if the matrix you wish computed, the `D` matrix, already exists.
1102: The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure
1104: Developer Notes:
1105: It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1106: Is there error checking for it?
1108: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`
1109: @*/
1110: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1111: {
1112: PetscFunctionBegin;
1117: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix A");
1118: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix B");
1120: if (C) {
1123: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix C");
1124: }
1126: PetscAssertPointer(D, 4);
1127: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), D));
1128: /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1129: PetscCall(MatProductCreate_Private(A, B, C, *D));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /*
1134: These are safe basic implementations of ABC, RARt and PtAP
1135: that do not rely on mat->ops->matmatop function pointers.
1136: They only use the MatProduct API and are currently used by
1137: cuSPARSE and KOKKOS-KERNELS backends
1138: */
1139: typedef struct {
1140: Mat BC;
1141: Mat ABC;
1142: } MatMatMatPrivate;
1144: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1145: {
1146: MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;
1148: PetscFunctionBegin;
1149: PetscCall(MatDestroy(&mmdata->BC));
1150: PetscCall(MatDestroy(&mmdata->ABC));
1151: PetscCall(PetscFree(data));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1156: {
1157: Mat_Product *product = mat->product;
1158: MatMatMatPrivate *mmabc;
1160: PetscFunctionBegin;
1161: MatCheckProduct(mat, 1);
1162: PetscCheck(mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data empty");
1163: mmabc = (MatMatMatPrivate *)mat->product->data;
1164: PetscCheck(mmabc->BC->ops->productnumeric, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing numeric stage");
1165: /* use function pointer directly to prevent logging */
1166: PetscCall((*mmabc->BC->ops->productnumeric)(mmabc->BC));
1167: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1168: mat->product = mmabc->ABC->product;
1169: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1170: /* use function pointer directly to prevent logging */
1171: PetscUseTypeMethod(mat, productnumeric);
1172: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1173: mat->product = product;
1174: PetscFunctionReturn(PETSC_SUCCESS);
1175: }
1177: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1178: {
1179: Mat_Product *product = mat->product;
1180: Mat A, B, C;
1181: MatProductType p1, p2;
1182: MatMatMatPrivate *mmabc;
1183: const char *prefix;
1185: PetscFunctionBegin;
1186: MatCheckProduct(mat, 1);
1187: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data not empty");
1188: PetscCall(MatGetOptionsPrefix(mat, &prefix));
1189: PetscCall(PetscNew(&mmabc));
1190: product->data = mmabc;
1191: product->destroy = MatDestroy_MatMatMatPrivate;
1192: switch (product->type) {
1193: case MATPRODUCT_PtAP:
1194: p1 = MATPRODUCT_AB;
1195: p2 = MATPRODUCT_AtB;
1196: A = product->B;
1197: B = product->A;
1198: C = product->B;
1199: break;
1200: case MATPRODUCT_RARt:
1201: p1 = MATPRODUCT_ABt;
1202: p2 = MATPRODUCT_AB;
1203: A = product->B;
1204: B = product->A;
1205: C = product->B;
1206: break;
1207: case MATPRODUCT_ABC:
1208: p1 = MATPRODUCT_AB;
1209: p2 = MATPRODUCT_AB;
1210: A = product->A;
1211: B = product->B;
1212: C = product->C;
1213: break;
1214: default:
1215: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1216: }
1217: PetscCall(MatProductCreate(B, C, NULL, &mmabc->BC));
1218: PetscCall(MatSetOptionsPrefix(mmabc->BC, prefix));
1219: PetscCall(MatAppendOptionsPrefix(mmabc->BC, "P1_"));
1220: PetscCall(MatProductSetType(mmabc->BC, p1));
1221: PetscCall(MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT));
1222: PetscCall(MatProductSetFill(mmabc->BC, product->fill));
1223: mmabc->BC->product->api_user = product->api_user;
1224: PetscCall(MatProductSetFromOptions(mmabc->BC));
1225: 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);
1226: /* use function pointer directly to prevent logging */
1227: PetscCall((*mmabc->BC->ops->productsymbolic)(mmabc->BC));
1229: PetscCall(MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC));
1230: PetscCall(MatSetOptionsPrefix(mmabc->ABC, prefix));
1231: PetscCall(MatAppendOptionsPrefix(mmabc->ABC, "P2_"));
1232: PetscCall(MatProductSetType(mmabc->ABC, p2));
1233: PetscCall(MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT));
1234: PetscCall(MatProductSetFill(mmabc->ABC, product->fill));
1235: mmabc->ABC->product->api_user = product->api_user;
1236: PetscCall(MatProductSetFromOptions(mmabc->ABC));
1237: 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);
1238: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1239: mat->product = mmabc->ABC->product;
1240: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1241: /* use function pointer directly to prevent logging */
1242: PetscUseTypeMethod(mat, productsymbolic);
1243: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1244: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1245: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1246: mat->product = product;
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /*@
1251: MatProductGetType - Returns the type of matrix-matrix product associated with computing values for the given matrix
1253: Not Collective
1255: Input Parameter:
1256: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1258: Output Parameter:
1259: . mtype - the `MatProductType`
1261: Level: intermediate
1263: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1264: @*/
1265: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1266: {
1267: PetscFunctionBegin;
1269: PetscAssertPointer(mtype, 2);
1270: *mtype = MATPRODUCT_UNSPECIFIED;
1271: if (mat->product) *mtype = mat->product->type;
1272: PetscFunctionReturn(PETSC_SUCCESS);
1273: }
1275: /*@
1276: MatProductGetMats - Returns the matrices associated with the matrix-matrix product associated with computing values for the given matrix
1278: Not Collective
1280: Input Parameter:
1281: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1283: Output Parameters:
1284: + A - the first matrix
1285: . B - the second matrix
1286: - C - the third matrix (may be `NULL` for some `MatProductType`)
1288: Level: intermediate
1290: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1291: @*/
1292: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1293: {
1294: PetscFunctionBegin;
1296: if (A) *A = mat->product ? mat->product->A : NULL;
1297: if (B) *B = mat->product ? mat->product->B : NULL;
1298: if (C) *C = mat->product ? mat->product->C : NULL;
1299: PetscFunctionReturn(PETSC_SUCCESS);
1300: }