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: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqAIJ(Mat);
43: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_MPIAIJ(Mat);
45: const char *const MatProductTypes[] = {"UNSPECIFIED", "AB", "AtB", "ABt", "PtAP", "RARt", "ABC"};
47: /* these are basic implementations relying on the old function pointers
48: * they are dangerous and should be removed in the future */
49: static PetscErrorCode MatProductNumeric_PtAP_Unsafe(Mat C)
50: {
51: Mat_Product *product = C->product;
52: Mat P = product->B, AP = product->Dwork;
54: PetscFunctionBegin;
55: /* AP = A*P */
56: PetscCall(MatProductNumeric(AP));
57: /* C = P^T*AP */
58: product->type = MATPRODUCT_AtB;
59: PetscCall((*C->ops->transposematmultnumeric)(P, AP, C));
60: product->type = MATPRODUCT_PtAP;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: static PetscErrorCode MatProductSymbolic_PtAP_Unsafe(Mat C)
65: {
66: Mat_Product *product = C->product;
67: Mat A = product->A, P = product->B, AP;
68: PetscReal fill = product->fill;
70: PetscFunctionBegin;
71: PetscCall(PetscInfo(C, "for A %s, P %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
72: /* AP = A*P */
73: PetscCall(MatProductCreate(A, P, NULL, &AP));
74: PetscCall(MatProductSetType(AP, MATPRODUCT_AB));
75: PetscCall(MatProductSetAlgorithm(AP, MATPRODUCTALGORITHMDEFAULT));
76: PetscCall(MatProductSetFill(AP, fill));
77: PetscCall(MatProductSetFromOptions(AP));
78: PetscCall(MatProductSymbolic(AP));
80: /* C = P^T*AP */
81: PetscCall(MatProductSetType(C, MATPRODUCT_AtB));
82: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
83: product->A = P;
84: product->B = AP;
85: PetscCall(MatProductSetFromOptions(C));
86: PetscCall(MatProductSymbolic(C));
88: /* resume user's original input matrix setting for A and B */
89: product->type = MATPRODUCT_PtAP;
90: product->A = A;
91: product->B = P;
92: product->Dwork = AP;
94: C->ops->productnumeric = MatProductNumeric_PtAP_Unsafe;
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode MatProductNumeric_RARt_Unsafe(Mat C)
99: {
100: Mat_Product *product = C->product;
101: Mat R = product->B, RA = product->Dwork;
103: PetscFunctionBegin;
104: /* RA = R*A */
105: PetscCall(MatProductNumeric(RA));
106: /* C = RA*R^T */
107: product->type = MATPRODUCT_ABt;
108: PetscCall((*C->ops->mattransposemultnumeric)(RA, R, C));
109: product->type = MATPRODUCT_RARt;
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode MatProductSymbolic_RARt_Unsafe(Mat C)
114: {
115: Mat_Product *product = C->product;
116: Mat A = product->A, R = product->B, RA;
117: PetscReal fill = product->fill;
119: PetscFunctionBegin;
120: PetscCall(PetscInfo(C, "for A %s, R %s is used\n", ((PetscObject)product->A)->type_name, ((PetscObject)product->B)->type_name));
121: /* RA = R*A */
122: PetscCall(MatProductCreate(R, A, NULL, &RA));
123: PetscCall(MatProductSetType(RA, MATPRODUCT_AB));
124: PetscCall(MatProductSetAlgorithm(RA, MATPRODUCTALGORITHMDEFAULT));
125: PetscCall(MatProductSetFill(RA, fill));
126: PetscCall(MatProductSetFromOptions(RA));
127: PetscCall(MatProductSymbolic(RA));
129: /* C = RA*R^T */
130: PetscCall(MatProductSetType(C, MATPRODUCT_ABt));
131: PetscCall(MatProductSetAlgorithm(C, MATPRODUCTALGORITHMDEFAULT));
132: product->A = RA;
133: PetscCall(MatProductSetFromOptions(C));
134: PetscCall(MatProductSymbolic(C));
136: /* resume user's original input matrix setting for A */
137: product->type = MATPRODUCT_RARt;
138: product->A = A;
139: product->Dwork = RA; /* save here so it will be destroyed with product C */
140: C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
145: {
146: Mat_Product *product = mat->product;
147: Mat A = product->A, BC = product->Dwork;
149: PetscFunctionBegin;
150: /* Numeric BC = B*C */
151: PetscCall(MatProductNumeric(BC));
152: /* Numeric mat = A*BC */
153: product->type = MATPRODUCT_AB;
154: PetscCall((*mat->ops->matmultnumeric)(A, BC, mat));
155: product->type = MATPRODUCT_ABC;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
160: {
161: Mat_Product *product = mat->product;
162: Mat B = product->B, C = product->C, BC;
163: PetscReal fill = product->fill;
165: PetscFunctionBegin;
166: 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));
167: /* Symbolic BC = B*C */
168: PetscCall(MatProductCreate(B, C, NULL, &BC));
169: PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
170: PetscCall(MatProductSetAlgorithm(BC, MATPRODUCTALGORITHMDEFAULT));
171: PetscCall(MatProductSetFill(BC, fill));
172: PetscCall(MatProductSetFromOptions(BC));
173: PetscCall(MatProductSymbolic(BC));
175: /* Symbolic mat = A*BC */
176: PetscCall(MatProductSetType(mat, MATPRODUCT_AB));
177: PetscCall(MatProductSetAlgorithm(mat, MATPRODUCTALGORITHMDEFAULT));
178: product->B = BC;
179: product->Dwork = BC;
180: PetscCall(MatProductSetFromOptions(mat));
181: PetscCall(MatProductSymbolic(mat));
183: /* resume user's original input matrix setting for B */
184: product->type = MATPRODUCT_ABC;
185: product->B = B;
186: mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
191: {
192: Mat_Product *product = mat->product;
194: PetscFunctionBegin;
195: switch (product->type) {
196: case MATPRODUCT_PtAP:
197: PetscCall(MatProductSymbolic_PtAP_Unsafe(mat));
198: break;
199: case MATPRODUCT_RARt:
200: PetscCall(MatProductSymbolic_RARt_Unsafe(mat));
201: break;
202: case MATPRODUCT_ABC:
203: PetscCall(MatProductSymbolic_ABC_Unsafe(mat));
204: break;
205: default:
206: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[product->type]);
207: }
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@
212: MatProductReplaceMats - Replace the input matrices for the matrix-matrix product operation inside the computed matrix
214: Collective
216: Input Parameters:
217: + A - the matrix or `NULL` if not being replaced
218: . B - the matrix or `NULL` if not being replaced
219: . C - the matrix or `NULL` if not being replaced
220: - D - the matrix whose values are computed via a matrix-matrix product operation
222: Level: intermediate
224: Note:
225: To reuse the symbolic phase, the input matrices must have exactly the same data structure as the replaced one.
226: If the type of any of the input matrices is different than what was previously used, or their symmetry flag changed but
227: the symbolic phase took advantage of their symmetry, the product is cleared and `MatProductSetFromOptions()`
228: and `MatProductSymbolic()` are invoked again.
230: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductSetFromOptions()`, `MatProductSymbolic()`, `MatProductClear()`
231: @*/
232: PetscErrorCode MatProductReplaceMats(Mat A, Mat B, Mat C, Mat D)
233: {
234: Mat_Product *product;
235: PetscBool flgA = PETSC_TRUE, flgB = PETSC_TRUE, flgC = PETSC_TRUE, isset, issym;
237: PetscFunctionBegin;
239: MatCheckProduct(D, 4);
240: product = D->product;
241: if (A) {
243: PetscCall(PetscObjectReference((PetscObject)A));
244: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, ((PetscObject)A)->type_name, &flgA));
245: PetscCall(MatIsSymmetricKnown(A, &isset, &issym));
246: 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 */
247: flgA = PETSC_FALSE;
248: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
249: }
250: PetscCall(MatDestroy(&product->A));
251: product->A = A;
252: }
253: if (B) {
255: PetscCall(PetscObjectReference((PetscObject)B));
256: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, ((PetscObject)B)->type_name, &flgB));
257: PetscCall(MatIsSymmetricKnown(B, &isset, &issym));
258: if (product->symbolic_used_the_fact_B_is_symmetric && isset && !issym) {
259: flgB = PETSC_FALSE;
260: product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
261: }
262: PetscCall(MatDestroy(&product->B));
263: product->B = B;
264: }
265: if (C) {
267: PetscCall(PetscObjectReference((PetscObject)C));
268: PetscCall(PetscObjectTypeCompare((PetscObject)product->C, ((PetscObject)C)->type_name, &flgC));
269: PetscCall(MatIsSymmetricKnown(C, &isset, &issym));
270: if (product->symbolic_used_the_fact_C_is_symmetric && isset && !issym) {
271: flgC = PETSC_FALSE;
272: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
273: }
274: PetscCall(MatDestroy(&product->C));
275: product->C = C;
276: }
277: /* Any of the replaced mats is of a different type, reset */
278: if (!flgA || !flgB || !flgC) {
279: if (D->product->destroy) PetscCall((*D->product->destroy)(&D->product->data));
280: D->product->destroy = NULL;
281: D->product->data = NULL;
282: if (D->ops->productnumeric || D->ops->productsymbolic) {
283: PetscCall(MatProductSetFromOptions(D));
284: PetscCall(MatProductSymbolic(D));
285: }
286: }
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
291: {
292: Mat_Product *product = C->product;
293: Mat A = product->A, B = product->B;
294: PetscInt k, K = B->cmap->N;
295: PetscBool t = PETSC_TRUE, iscuda = PETSC_FALSE;
296: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
297: char *Btype = NULL, *Ctype = NULL;
299: PetscFunctionBegin;
300: switch (product->type) {
301: case MATPRODUCT_AB:
302: t = PETSC_FALSE;
303: case MATPRODUCT_AtB:
304: break;
305: default:
306: 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);
307: }
308: if (PetscDefined(HAVE_CUDA)) {
309: VecType vtype;
311: PetscCall(MatGetVecType(A, &vtype));
312: PetscCall(PetscStrcmp(vtype, VECCUDA, &iscuda));
313: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECSEQCUDA, &iscuda));
314: if (!iscuda) PetscCall(PetscStrcmp(vtype, VECMPICUDA, &iscuda));
315: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
316: PetscCall(PetscStrallocpy(((PetscObject)B)->type_name, &Btype));
317: PetscCall(PetscStrallocpy(((PetscObject)C)->type_name, &Ctype));
318: PetscCall(MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B));
319: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
320: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
321: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
322: }
323: PetscCall(MatConvert(C, MATDENSECUDA, MAT_INPLACE_MATRIX, &C));
324: } else { /* Make sure we have up-to-date data on the CPU */
325: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_VIENNACL)
326: Bcpu = B->boundtocpu;
327: Ccpu = C->boundtocpu;
328: #endif
329: PetscCall(MatBindToCPU(B, PETSC_TRUE));
330: PetscCall(MatBindToCPU(C, PETSC_TRUE));
331: }
332: }
333: for (k = 0; k < K; k++) {
334: Vec x, y;
336: PetscCall(MatDenseGetColumnVecRead(B, k, &x));
337: PetscCall(MatDenseGetColumnVecWrite(C, k, &y));
338: if (t) {
339: PetscCall(MatMultTranspose(A, x, y));
340: } else {
341: PetscCall(MatMult(A, x, y));
342: }
343: PetscCall(MatDenseRestoreColumnVecRead(B, k, &x));
344: PetscCall(MatDenseRestoreColumnVecWrite(C, k, &y));
345: }
346: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
347: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
348: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
349: if (PetscDefined(HAVE_CUDA)) {
350: if (iscuda) {
351: PetscCall(MatConvert(B, Btype, MAT_INPLACE_MATRIX, &B));
352: PetscCall(MatConvert(C, Ctype, MAT_INPLACE_MATRIX, &C));
353: } else {
354: PetscCall(MatBindToCPU(B, Bcpu));
355: PetscCall(MatBindToCPU(C, Ccpu));
356: }
357: }
358: PetscCall(PetscFree(Btype));
359: PetscCall(PetscFree(Ctype));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
364: {
365: Mat_Product *product = C->product;
366: Mat A = product->A, B = product->B;
367: PetscBool isdense;
369: PetscFunctionBegin;
370: switch (product->type) {
371: case MATPRODUCT_AB:
372: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
373: break;
374: case MATPRODUCT_AtB:
375: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
376: break;
377: default:
378: 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);
379: }
380: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
381: if (!isdense) {
382: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
383: PetscCall(MatSetVecType(C, B->defaultvectype));
384: /* If matrix type of C was not set or not dense, we need to reset the pointer */
385: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
386: }
387: C->ops->productnumeric = MatProductNumeric_X_Dense;
388: PetscCall(MatSetUp(C));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /* a single driver to query the dispatching */
393: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
394: {
395: Mat_Product *product = mat->product;
396: PetscInt Am, An, Bm, Bn, Cm, Cn;
397: Mat A = product->A, B = product->B, C = product->C;
398: const char *const Bnames[] = {"B", "R", "P"};
399: const char *bname;
400: PetscErrorCode (*fA)(Mat);
401: PetscErrorCode (*fB)(Mat);
402: PetscErrorCode (*fC)(Mat);
403: PetscErrorCode (*f)(Mat) = NULL;
405: PetscFunctionBegin;
406: mat->ops->productsymbolic = NULL;
407: mat->ops->productnumeric = NULL;
408: if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(PETSC_SUCCESS);
409: PetscCheck(A, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing A mat");
410: PetscCheck(B, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing B mat");
411: PetscCheck(product->type != MATPRODUCT_ABC || C, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing C mat");
412: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
413: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
414: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
415: else bname = Bnames[0];
417: /* Check matrices sizes */
418: Am = A->rmap->N;
419: An = A->cmap->N;
420: Bm = B->rmap->N;
421: Bn = B->cmap->N;
422: Cm = C ? C->rmap->N : 0;
423: Cn = C ? C->cmap->N : 0;
424: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) {
425: PetscInt t = Bn;
426: Bn = Bm;
427: Bm = t;
428: }
429: if (product->type == MATPRODUCT_AtB) An = Am;
431: 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,
432: MatProductTypes[product->type], A->rmap->N, A->cmap->N, bname, B->rmap->N, B->cmap->N);
433: 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,
434: MatProductTypes[product->type], B->rmap->N, B->cmap->N, Cm, Cn);
436: /* graph products use the host CSR structure, including for AIJ device types */
437: if (mat->structure_only && (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB)) {
438: PetscBool aseq, bseq, ampi, bmpi;
440: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &aseq));
441: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &bseq));
442: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &i));
443: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATMPIAIJ, &bmpi));
444: if (aseq && bseq) {
445: PetscCall(MatProductSetFromOptions_SeqAIJ(mat));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: } else if (ampi && bmpi) {
448: PetscCall(MatProductSetFromOptions_MPIAIJ(mat));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
451: }
453: fA = A->ops->productsetfromoptions;
454: fB = B->ops->productsetfromoptions;
455: fC = C ? C->ops->productsetfromoptions : fA;
456: if (C) {
457: 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));
458: } else {
459: PetscCall(PetscInfo(mat, "MatProductType %s for A %s, %s %s\n", MatProductTypes[product->type], ((PetscObject)A)->type_name, bname, ((PetscObject)B)->type_name));
460: }
461: if (fA == fB && fA == fC && fA) {
462: PetscCall(PetscInfo(mat, " matching op\n"));
463: PetscCall((*fA)(mat));
464: }
465: /* We may have found f but it did not succeed */
466: if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
467: char mtypes[256];
468: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_", sizeof(mtypes)));
469: PetscCall(PetscStrlcat(mtypes, ((PetscObject)A)->type_name, sizeof(mtypes)));
470: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
471: PetscCall(PetscStrlcat(mtypes, ((PetscObject)B)->type_name, sizeof(mtypes)));
472: if (C) {
473: PetscCall(PetscStrlcat(mtypes, "_", sizeof(mtypes)));
474: PetscCall(PetscStrlcat(mtypes, ((PetscObject)C)->type_name, sizeof(mtypes)));
475: }
476: PetscCall(PetscStrlcat(mtypes, "_C", sizeof(mtypes)));
477: #if defined(__clang__)
478: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
479: #elif defined(__GNUC__) || defined(__GNUG__)
480: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
481: #endif
482: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
483: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
484: if (!f) {
485: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
486: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
487: }
488: if (!f && C) {
489: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
490: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
491: }
492: if (f) PetscCall((*f)(mat));
494: /* We may have found f but it did not succeed */
495: /* some matrices (i.e. MATTRANSPOSEVIRTUAL, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
496: if (!mat->ops->productsymbolic) {
497: PetscCall(PetscStrncpy(mtypes, "MatProductSetFromOptions_anytype_C", sizeof(mtypes)));
498: PetscCall(PetscObjectQueryFunction((PetscObject)A, mtypes, &f));
499: PetscCall(PetscInfo(mat, " querying %s from A? %p\n", mtypes, f));
500: if (!f) {
501: PetscCall(PetscObjectQueryFunction((PetscObject)B, mtypes, &f));
502: PetscCall(PetscInfo(mat, " querying %s from %s? %p\n", mtypes, bname, f));
503: }
504: if (!f && C) {
505: PetscCall(PetscObjectQueryFunction((PetscObject)C, mtypes, &f));
506: PetscCall(PetscInfo(mat, " querying %s from C? %p\n", mtypes, f));
507: }
508: }
509: if (f) PetscCall((*f)(mat));
510: }
511: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
512: /* We may have found f but it did not succeed */
513: if (!mat->ops->productsymbolic) {
514: /* we can still compute the product if B is of type dense */
515: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
516: PetscBool isdense;
518: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
519: if (isdense) {
520: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
521: PetscCall(PetscInfo(mat, " using basic looping over columns of a dense matrix\n"));
522: }
523: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
524: /*
525: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
526: the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
527: before computing the symbolic phase
528: */
529: PetscCall(PetscInfo(mat, " symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n"));
530: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
531: }
532: }
533: if (!mat->ops->productsymbolic) PetscCall(PetscInfo(mat, " symbolic product is not supported\n"));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: MatProductSetFromOptions - Sets the options for the computation of a matrix-matrix product operation where the type,
539: the algorithm etc are determined from the options database.
541: Logically Collective
543: Input Parameter:
544: . mat - the matrix whose values are computed via a matrix-matrix product operation
546: Options Database Keys:
547: + -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
548: . -mat_product_algorithm algorithm - Sets the algorithm, see `MatProductAlgorithm` for possible values
549: - -mat_product_algorithm_backend_cpu - Use the CPU to perform the computation even if the matrix is a GPU matrix
551: Level: intermediate
553: Note:
554: The `-mat_product_clear` option reduces memory usage but means that the matrix cannot be re-used for a matrix-matrix product operation
556: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatSetFromOptions()`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductNumeric()`,
557: `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductAlgorithm`
558: @*/
559: PetscErrorCode MatProductSetFromOptions(Mat mat)
560: {
561: PetscFunctionBegin;
563: MatCheckProduct(mat, 1);
564: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot call MatProductSetFromOptions() with already present data");
565: mat->product->setfromoptionscalled = PETSC_TRUE;
566: PetscObjectOptionsBegin((PetscObject)mat);
567: PetscCall(PetscOptionsBool("-mat_product_clear", "Clear intermediate data structures after MatProductNumeric() has been called", "MatProductClear", mat->product->clear, &mat->product->clear, NULL));
568: PetscCall(PetscOptionsDeprecated("-mat_freeintermediatedatastructures", "-mat_product_clear", "3.13", "Or call MatProductClear() after MatProductNumeric()"));
569: PetscOptionsEnd();
570: PetscCall(MatProductSetFromOptions_Private(mat));
571: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing product after setup phase");
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*@
576: MatProductView - View the private matrix-matrix algorithm object within a matrix
578: Logically Collective
580: Input Parameters:
581: + mat - the matrix obtained with `MatProductCreate()` or `MatProductCreateWithMat()`
582: - viewer - where the information on the matrix-matrix algorithm of `mat` should be reviewed
584: Level: intermediate
586: Developer Note:
587: Shouldn't this information be printed from an appropriate `MatView()` with perhaps certain formats set?
589: .seealso: [](ch_matrices), `MatProductType`, `Mat`, `MatProductSetFromOptions()`, `MatView()`, `MatProductCreate()`, `MatProductCreateWithMat()`
590: @*/
591: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
592: {
593: PetscFunctionBegin;
595: if (!mat->product) PetscFunctionReturn(PETSC_SUCCESS);
596: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat), &viewer));
598: PetscCheckSameComm(mat, 1, viewer, 2);
599: if (mat->product->view) PetscCall((*mat->product->view)(mat, viewer));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /* these are basic implementations relying on the old function pointers
604: * they are dangerous and should be removed in the future */
605: PetscErrorCode MatProductNumeric_AB(Mat mat)
606: {
607: Mat_Product *product = mat->product;
608: Mat A = product->A, B = product->B;
610: PetscFunctionBegin;
611: PetscCall((*mat->ops->matmultnumeric)(A, B, mat));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: PetscErrorCode MatProductNumeric_AtB(Mat mat)
616: {
617: Mat_Product *product = mat->product;
618: Mat A = product->A, B = product->B;
620: PetscFunctionBegin;
621: PetscCall((*mat->ops->transposematmultnumeric)(A, B, mat));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: PetscErrorCode MatProductNumeric_ABt(Mat mat)
626: {
627: Mat_Product *product = mat->product;
628: Mat A = product->A, B = product->B;
630: PetscFunctionBegin;
631: PetscCall((*mat->ops->mattransposemultnumeric)(A, B, mat));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
636: {
637: Mat_Product *product = mat->product;
638: Mat A = product->A, B = product->B;
640: PetscFunctionBegin;
641: PetscCall((*mat->ops->ptapnumeric)(A, B, mat));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: PetscErrorCode MatProductNumeric_RARt(Mat mat)
646: {
647: Mat_Product *product = mat->product;
648: Mat A = product->A, B = product->B;
650: PetscFunctionBegin;
651: PetscCall((*mat->ops->rartnumeric)(A, B, mat));
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: PetscErrorCode MatProductNumeric_ABC(Mat mat)
656: {
657: Mat_Product *product = mat->product;
658: Mat A = product->A, B = product->B, C = product->C;
660: PetscFunctionBegin;
661: PetscCall((*mat->ops->matmatmultnumeric)(A, B, C, mat));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: /*@
666: MatProductNumeric - Compute a matrix-matrix product operation with the numerical values
668: Collective
670: Input/Output Parameter:
671: . mat - the matrix whose values are computed via a matrix-matrix product operation
673: Level: intermediate
675: Note:
676: `MatProductSymbolic()` must have been called on `mat` before calling this function
678: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`, `MatProductSetType()`, `MatProductCreate()`, `MatSetType()`, `MatProductSymbolic()`
679: @*/
680: PetscErrorCode MatProductNumeric(Mat mat)
681: {
682: PetscLogEvent eventtype = -1;
684: PetscFunctionBegin;
686: MatCheckProduct(mat, 1);
687: PetscCheck(!mat->structure_only, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot compute numerical values of a MAT_STRUCTURE_ONLY product");
688: switch (mat->product->type) {
689: case MATPRODUCT_AB:
690: eventtype = MAT_MatMultNumeric;
691: break;
692: case MATPRODUCT_AtB:
693: eventtype = MAT_TransposeMatMultNumeric;
694: break;
695: case MATPRODUCT_ABt:
696: eventtype = MAT_MatTransposeMultNumeric;
697: break;
698: case MATPRODUCT_PtAP:
699: eventtype = MAT_PtAPNumeric;
700: break;
701: case MATPRODUCT_RARt:
702: eventtype = MAT_RARtNumeric;
703: break;
704: case MATPRODUCT_ABC:
705: eventtype = MAT_MatMatMultNumeric;
706: break;
707: default:
708: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
709: }
711: if (mat->ops->productnumeric) {
712: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
713: PetscUseTypeMethod(mat, productnumeric);
714: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
715: } else if (mat->product) {
716: char errstr[256];
718: if (mat->product->type == MATPRODUCT_ABC) {
719: 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));
720: } else {
721: 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));
722: }
723: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified numeric phase for product %s", errstr);
724: }
725: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after numeric phase for product");
727: if (mat->product->clear) PetscCall(MatProductClear(mat));
728: PetscCall(PetscObjectStateIncrease((PetscObject)mat));
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: /* these are basic implementations relying on the old function pointers
733: * they are dangerous and should be removed in the future */
734: PetscErrorCode MatProductSymbolic_AB(Mat mat)
735: {
736: Mat_Product *product = mat->product;
737: Mat A = product->A, B = product->B;
739: PetscFunctionBegin;
740: PetscCall((*mat->ops->matmultsymbolic)(A, B, product->fill, mat));
741: mat->ops->productnumeric = MatProductNumeric_AB;
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
746: {
747: Mat_Product *product = mat->product;
748: Mat A = product->A, B = product->B;
750: PetscFunctionBegin;
751: PetscCall((*mat->ops->transposematmultsymbolic)(A, B, product->fill, mat));
752: mat->ops->productnumeric = MatProductNumeric_AtB;
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
757: {
758: Mat_Product *product = mat->product;
759: Mat A = product->A, B = product->B;
761: PetscFunctionBegin;
762: PetscCall((*mat->ops->mattransposemultsymbolic)(A, B, product->fill, mat));
763: mat->ops->productnumeric = MatProductNumeric_ABt;
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
768: {
769: Mat_Product *product = mat->product;
770: Mat A = product->A, B = product->B, C = product->C;
772: PetscFunctionBegin;
773: PetscCall((*mat->ops->matmatmultsymbolic)(A, B, C, product->fill, mat));
774: mat->ops->productnumeric = MatProductNumeric_ABC;
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: /*@
779: MatProductSymbolic - Perform the symbolic portion of a matrix-matrix product operation, this creates a data structure for use with the numerical
780: product to be done with `MatProductNumeric()`
782: Collective
784: Input/Output Parameter:
785: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
787: Level: intermediate
789: Notes:
790: `MatProductSetFromOptions()` must have been called on `mat` before calling this function.
791: For `MATPRODUCT_AB` and `MATPRODUCT_AtB` with AIJ inputs, setting `MAT_STRUCTURE_ONLY` on `mat` before
792: `MatProductSetFromOptions()` requests an assembled graph without numerical values. These products use the
793: native CPU symbolic algorithms (`sorted` for sequential AB, `outerproduct` for sequential AtB, and `scalable`
794: for MPI products), overriding the selected algorithm, and return `MATSEQAIJ` or `MATMPIAIJ`, including for device AIJ inputs.
795: Call `MatProductClear()` when the product context is no longer needed. `MatProductNumeric()` cannot be used on a structure-only result.
797: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductCreateWithMat()`, `MatProductSetFromOptions()`, `MatProductNumeric()`, `MatProductSetType()`, `MatProductSetAlgorithm()`
798: @*/
799: PetscErrorCode MatProductSymbolic(Mat mat)
800: {
801: PetscLogEvent eventtype = -1;
802: PetscBool missing = PETSC_FALSE;
803: Mat_Product *product = mat->product;
804: Mat A = product->A;
805: Mat B = product->B;
806: Mat C = product->C;
808: PetscFunctionBegin;
810: MatCheckProduct(mat, 1);
811: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Cannot run symbolic phase. Product data not empty");
812: switch (mat->product->type) {
813: case MATPRODUCT_AB:
814: eventtype = MAT_MatMultSymbolic;
815: break;
816: case MATPRODUCT_AtB:
817: eventtype = MAT_TransposeMatMultSymbolic;
818: break;
819: case MATPRODUCT_ABt:
820: eventtype = MAT_MatTransposeMultSymbolic;
821: break;
822: case MATPRODUCT_PtAP:
823: eventtype = MAT_PtAPSymbolic;
824: break;
825: case MATPRODUCT_RARt:
826: eventtype = MAT_RARtSymbolic;
827: break;
828: case MATPRODUCT_ABC:
829: eventtype = MAT_MatMatMultSymbolic;
830: break;
831: default:
832: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[mat->product->type]);
833: }
834: mat->ops->productnumeric = NULL;
835: if (mat->ops->productsymbolic) {
836: PetscCall(PetscLogEventBegin(eventtype, mat, 0, 0, 0));
837: PetscUseTypeMethod(mat, productsymbolic);
838: PetscCall(PetscLogEventEnd(eventtype, mat, 0, 0, 0));
839: } else missing = PETSC_TRUE;
840: if (missing || !mat->product || !mat->ops->productnumeric) {
841: char errstr[256];
843: if (mat->product->type == MATPRODUCT_ABC) {
844: 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));
845: } else {
846: 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));
847: }
848: PetscCheck(mat->product->setfromoptionscalled, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() first", errstr);
849: PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unspecified symbolic phase for product %s. The product is not supported", errstr);
850: PetscCheck(mat->product, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing struct after symbolic phase for product %s", errstr);
851: }
852: #if PetscDefined(HAVE_DEVICE)
853: PetscBool bindingpropagates;
854: bindingpropagates = (PetscBool)((A->boundtocpu && A->bindingpropagates) || (B->boundtocpu && B->bindingpropagates));
855: if (C) bindingpropagates = (PetscBool)(bindingpropagates || (C->boundtocpu && C->bindingpropagates));
856: if (bindingpropagates) {
857: PetscCall(MatBindToCPU(mat, PETSC_TRUE));
858: PetscCall(MatSetBindingPropagates(mat, PETSC_TRUE));
859: }
860: #endif
861: /* set block sizes */
862: switch (product->type) {
863: case MATPRODUCT_PtAP:
864: if (B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->cmap->bs, B->cmap->bs));
865: break;
866: case MATPRODUCT_RARt:
867: if (B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, B->rmap->bs, B->rmap->bs));
868: break;
869: case MATPRODUCT_ABC:
870: PetscCall(MatSetBlockSizesFromMats(mat, A, C));
871: break;
872: case MATPRODUCT_AB:
873: PetscCall(MatSetBlockSizesFromMats(mat, A, B));
874: break;
875: case MATPRODUCT_AtB:
876: if (A->cmap->bs > 1 || B->cmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, B->cmap->bs));
877: break;
878: case MATPRODUCT_ABt:
879: if (A->rmap->bs > 1 || B->rmap->bs > 1) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, B->rmap->bs));
880: break;
881: default:
882: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
883: }
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: MatProductSetFill - Set an expected fill of the matrix whose values are computed via a matrix-matrix product operation
890: Collective
892: Input Parameters:
893: + mat - the matrix whose values are to be computed via a matrix-matrix product operation
894: - 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.
895: If the product is a dense matrix, this value is not used.
897: Level: intermediate
899: Notes:
900: Use `fill` of `PETSC_DETERMINE` to use the default value.
902: The deprecated `PETSC_DEFAULT` is also supported to mean use the current value.
904: .seealso: [](ch_matrices), `MatProduct`, `PETSC_DETERMINE`, `Mat`, `MatProductSetFromOptions()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
905: @*/
906: PetscErrorCode MatProductSetFill(Mat mat, PetscReal fill)
907: {
908: PetscFunctionBegin;
910: MatCheckProduct(mat, 1);
911: if (fill == (PetscReal)PETSC_DETERMINE) mat->product->fill = mat->product->default_fill;
912: else if (fill != (PetscReal)PETSC_CURRENT) mat->product->fill = fill;
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: /*@
917: MatProductSetAlgorithm - Requests a particular algorithm for a matrix-matrix product operation that will perform to compute the given matrix
919: Collective
921: Input Parameters:
922: + mat - the matrix whose values are computed via a matrix-matrix product operation
923: - alg - particular implementation algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
925: Options Database Key:
926: . -mat_product_algorithm algorithm - Sets the algorithm, see `MatProductAlgorithm`
928: Level: intermediate
930: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductClear()`, `MatProductSetType()`, `MatProductSetFill()`, `MatProductCreate()`, `MatProductAlgorithm`, `MatProductType`, `MatProductGetAlgorithm()`
931: @*/
932: PetscErrorCode MatProductSetAlgorithm(Mat mat, MatProductAlgorithm alg)
933: {
934: PetscFunctionBegin;
936: MatCheckProduct(mat, 1);
937: PetscCall(PetscFree(mat->product->alg));
938: PetscCall(PetscStrallocpy(alg, &mat->product->alg));
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: /*@
943: MatProductGetAlgorithm - Returns the selected algorithm for a matrix-matrix product operation
945: Not Collective
947: Input Parameter:
948: . mat - the matrix whose values are computed via a matrix-matrix product operation
950: Output Parameter:
951: . alg - the selected algorithm of the matrix product, e.g., `MATPRODUCTALGORITHMDEFAULT`.
953: Level: intermediate
955: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductSetAlgorithm()`
956: @*/
957: PetscErrorCode MatProductGetAlgorithm(Mat mat, MatProductAlgorithm *alg)
958: {
959: PetscFunctionBegin;
961: PetscAssertPointer(alg, 2);
962: if (mat->product) *alg = mat->product->alg;
963: else *alg = NULL;
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: /*@
968: MatProductSetType - Sets a particular matrix-matrix product operation to be used to compute the values of the given matrix
970: Collective
972: Input Parameters:
973: + mat - the matrix whose values are computed via a matrix-matrix product operation
974: - productype - matrix product type, e.g., `MATPRODUCT_AB`,`MATPRODUCT_AtB`,`MATPRODUCT_ABt`,`MATPRODUCT_PtAP`,`MATPRODUCT_RARt`,`MATPRODUCT_ABC`,
975: see `MatProductType`
977: Level: intermediate
979: Note:
980: The small t represents the transpose operation.
982: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`, `MatProductType`,
983: `MATPRODUCT_AB`, `MATPRODUCT_AtB`, `MATPRODUCT_ABt`, `MATPRODUCT_PtAP`, `MATPRODUCT_RARt`, `MATPRODUCT_ABC`
984: @*/
985: PetscErrorCode MatProductSetType(Mat mat, MatProductType productype)
986: {
987: PetscFunctionBegin;
989: MatCheckProduct(mat, 1);
991: if (productype != mat->product->type) {
992: if (mat->product->destroy) PetscCall((*mat->product->destroy)(&mat->product->data));
993: mat->product->destroy = NULL;
994: mat->product->data = NULL;
995: mat->ops->productsymbolic = NULL;
996: mat->ops->productnumeric = NULL;
997: }
998: mat->product->type = productype;
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: /*@
1003: MatProductClear - Clears from the matrix any internal data structures related to the computation of the values of the matrix from matrix-matrix product operations
1005: Collective
1007: Input Parameter:
1008: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1010: Options Database Key:
1011: . -mat_product_clear - Clear intermediate data structures after `MatProductNumeric()` has been called
1013: Level: intermediate
1015: Notes:
1016: This function should be called to remove any intermediate data used to compute the matrix to free up memory.
1018: After having called this function, matrix-matrix product operations can no longer be used on `mat`
1020: Developer Note:
1021: This frees the `Mat_Product` context that was attached to the matrix during `MatProductCreate()` or `MatProductCreateWithMat()`
1023: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreate()`
1024: @*/
1025: PetscErrorCode MatProductClear(Mat mat)
1026: {
1027: Mat_Product *product = mat->product;
1029: PetscFunctionBegin;
1031: if (product) {
1032: PetscCall(MatDestroy(&product->A));
1033: PetscCall(MatDestroy(&product->B));
1034: PetscCall(MatDestroy(&product->C));
1035: PetscCall(PetscFree(product->alg));
1036: PetscCall(MatDestroy(&product->Dwork));
1037: if (product->destroy) PetscCall((*product->destroy)(&product->data));
1038: }
1039: PetscCall(PetscFree(mat->product));
1040: mat->ops->productsymbolic = NULL;
1041: mat->ops->productnumeric = NULL;
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /* Create a supporting struct and attach it to the matrix product */
1046: PetscErrorCode MatProductCreate_Private(Mat A, Mat B, Mat C, Mat D)
1047: {
1048: Mat_Product *product = NULL;
1050: PetscFunctionBegin;
1052: PetscCheck(!D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product already present");
1053: PetscCall(PetscNew(&product));
1054: product->A = A;
1055: product->B = B;
1056: product->C = C;
1057: product->type = MATPRODUCT_UNSPECIFIED;
1058: product->Dwork = NULL;
1059: product->api_user = PETSC_FALSE;
1060: product->clear = PETSC_FALSE;
1061: product->setfromoptionscalled = PETSC_FALSE;
1062: PetscObjectParameterSetDefault(product, fill, 2);
1063: D->product = product;
1065: PetscCall(MatProductSetAlgorithm(D, MATPRODUCTALGORITHMDEFAULT));
1066: PetscCall(MatProductSetFill(D, PETSC_DEFAULT));
1068: PetscCall(PetscObjectReference((PetscObject)A));
1069: PetscCall(PetscObjectReference((PetscObject)B));
1070: PetscCall(PetscObjectReference((PetscObject)C));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@
1075: MatProductCreateWithMat - Set a given matrix to have its values computed via matrix-matrix operations on other matrices.
1077: Collective
1079: Input Parameters:
1080: + A - the first matrix
1081: . B - the second matrix
1082: . C - the third matrix (optional, use `NULL` if not needed)
1083: - D - the matrix whose values are to be computed via a matrix-matrix product operation
1085: Level: intermediate
1087: Notes:
1088: Use `MatProductCreate()` if the matrix you wish computed `D` does not exist
1090: See `MatProductCreate()` for details on the usage of the matrix-matrix product operations
1092: Any product data currently attached to `D` will be freed
1094: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductType`, `MatProductSetType()`, `MatProductAlgorithm`,
1095: `MatProductSetAlgorithm`, `MatProductCreate()`, `MatProductClear()`
1096: @*/
1097: PetscErrorCode MatProductCreateWithMat(Mat A, Mat B, Mat C, Mat D)
1098: {
1099: PetscFunctionBegin;
1102: MatCheckPreallocated(A, 1);
1103: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1104: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1108: MatCheckPreallocated(B, 2);
1109: PetscCheck(B->assembled, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1110: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1112: if (C) {
1115: MatCheckPreallocated(C, 3);
1116: PetscCheck(C->assembled, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1117: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1118: }
1122: MatCheckPreallocated(D, 4);
1123: PetscCheck(D->assembled, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
1124: PetscCheck(!D->factortype, PetscObjectComm((PetscObject)D), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1126: /* Create a supporting struct and attach it to D */
1127: PetscCall(MatProductClear(D));
1128: PetscCall(MatProductCreate_Private(A, B, C, D));
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }
1132: /*@
1133: MatProductCreate - create a matrix to hold the result of a matrix-matrix (or matrix-matrix-matrix) product operation
1135: Collective
1137: Input Parameters:
1138: + A - the first matrix
1139: . B - the second matrix
1140: - C - the third matrix (or `NULL`)
1142: Output Parameter:
1143: . D - the matrix whose values are to be computed via a matrix-matrix product operation
1145: Level: intermediate
1147: Example:
1148: .vb
1149: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
1150: MatProductSetType(D, MATPRODUCT_AB or MATPRODUCT_AtB or MATPRODUCT_ABt or MATPRODUCT_PtAP or MATPRODUCT_RARt or MATPRODUCT_ABC)
1151: MatProductSetAlgorithm(D, alg)
1152: MatProductSetFill(D,fill)
1153: MatProductSetFromOptions(D)
1154: MatProductSymbolic(D)
1155: MatProductNumeric(D)
1156: Change numerical values in some of the matrices
1157: MatProductNumeric(D)
1158: .ve
1160: Notes:
1161: Use `MatProductCreateWithMat()` if `D` the matrix you wish computed already exists.
1163: The information computed during the symbolic stage can be reused for new numerical computations with the same non-zero structure of the input matrices.
1165: Developer Notes:
1166: It is undocumented what happens if the nonzero structure of the input matrices changes. Is the symbolic stage automatically redone? Does it crash?
1167: Is there error checking for it?
1169: 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
1170: information.
1172: Each `MatProductAlgorithm` associated with a particular `MatType` stores additional data needed for the product computation
1173: (generally this data is computed in `MatProductSymbolic()`) inside the `Mat_Product` context in a `MatProductCtx_XXX` data structure
1174: and provides a `MatProductCtxDestroy_XXX()` routine to free that data. The `MatProductAlgorithm` and `MatType` specific destroy routine is called by
1175: `MatProductClear()`.
1177: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductClear()`,
1178: `MatProductSymbolic()`, `MatProductNumeric()`, `MatProductAlgorithm`, `MatProductType`
1179: @*/
1180: PetscErrorCode MatProductCreate(Mat A, Mat B, Mat C, Mat *D)
1181: {
1182: PetscFunctionBegin;
1187: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix A");
1188: PetscCheck(!B->factortype, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix B");
1190: if (C) {
1193: PetscCheck(!C->factortype, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix C");
1194: }
1196: PetscAssertPointer(D, 4);
1197: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), D));
1198: /* Delay setting type of D to the MatProduct symbolic phase, as we allow sparse A and dense B */
1199: PetscCall(MatProductCreate_Private(A, B, C, *D));
1200: PetscFunctionReturn(PETSC_SUCCESS);
1201: }
1203: /*
1204: These are safe basic implementations of ABC, RARt and PtAP
1205: that do not rely on mat->ops->matmatop function pointers.
1206: They only use the MatProduct API and are currently used by
1207: cuSPARSE and KOKKOS-KERNELS backends
1208: */
1209: typedef struct {
1210: Mat BC;
1211: Mat ABC;
1212: } MatProductCtx_MatMatMatPrivate;
1214: static PetscErrorCode MatProductCtxDestroy_MatMatMatPrivate(PetscCtxRt data)
1215: {
1216: MatProductCtx_MatMatMatPrivate *mmdata = *(MatProductCtx_MatMatMatPrivate **)data;
1218: PetscFunctionBegin;
1219: PetscCall(MatDestroy(&mmdata->BC));
1220: PetscCall(MatDestroy(&mmdata->ABC));
1221: PetscCall(PetscFree(mmdata));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1226: {
1227: Mat_Product *product = mat->product;
1228: MatProductCtx_MatMatMatPrivate *mmabc;
1230: PetscFunctionBegin;
1231: MatCheckProduct(mat, 1);
1232: PetscCheck(mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data empty");
1233: mmabc = (MatProductCtx_MatMatMatPrivate *)mat->product->data;
1234: PetscCheck(mmabc->BC->ops->productnumeric, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Missing numeric stage");
1235: /* use function pointer directly to prevent logging */
1236: PetscCall((*mmabc->BC->ops->productnumeric)(mmabc->BC));
1237: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1238: mat->product = mmabc->ABC->product;
1239: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1240: /* use function pointer directly to prevent logging */
1241: PetscUseTypeMethod(mat, productnumeric);
1242: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1243: mat->product = product;
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1248: {
1249: Mat_Product *product = mat->product;
1250: Mat A, B, C;
1251: MatProductType p1, p2;
1252: MatProductCtx_MatMatMatPrivate *mmabc;
1253: const char *prefix;
1255: PetscFunctionBegin;
1256: MatCheckProduct(mat, 1);
1257: PetscCheck(!mat->product->data, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Product data not empty");
1258: PetscCall(MatGetOptionsPrefix(mat, &prefix));
1259: PetscCall(PetscNew(&mmabc));
1260: product->data = mmabc;
1261: product->destroy = MatProductCtxDestroy_MatMatMatPrivate;
1262: switch (product->type) {
1263: case MATPRODUCT_PtAP:
1264: p1 = MATPRODUCT_AB;
1265: p2 = MATPRODUCT_AtB;
1266: A = product->B;
1267: B = product->A;
1268: C = product->B;
1269: if (A->cmap->bs > 0 && C->cmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->cmap->bs, C->cmap->bs));
1270: break;
1271: case MATPRODUCT_RARt:
1272: p1 = MATPRODUCT_ABt;
1273: p2 = MATPRODUCT_AB;
1274: A = product->B;
1275: B = product->A;
1276: C = product->B;
1277: if (A->rmap->bs > 0 && C->rmap->bs > 0) PetscCall(MatSetBlockSizes(mat, A->rmap->bs, C->rmap->bs));
1278: break;
1279: case MATPRODUCT_ABC:
1280: p1 = MATPRODUCT_AB;
1281: p2 = MATPRODUCT_AB;
1282: A = product->A;
1283: B = product->B;
1284: C = product->C;
1285: PetscCall(MatSetBlockSizesFromMats(mat, A, C));
1286: break;
1287: default:
1288: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "Not for ProductType %s", MatProductTypes[product->type]);
1289: }
1290: PetscCall(MatProductCreate(B, C, NULL, &mmabc->BC));
1291: PetscCall(MatSetOptionsPrefix(mmabc->BC, prefix));
1292: PetscCall(MatAppendOptionsPrefix(mmabc->BC, "P1_"));
1293: PetscCall(MatProductSetType(mmabc->BC, p1));
1294: PetscCall(MatProductSetAlgorithm(mmabc->BC, MATPRODUCTALGORITHMDEFAULT));
1295: PetscCall(MatProductSetFill(mmabc->BC, product->fill));
1296: mmabc->BC->product->api_user = product->api_user;
1297: PetscCall(MatProductSetFromOptions(mmabc->BC));
1298: 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);
1299: /* use function pointer directly to prevent logging */
1300: PetscCall((*mmabc->BC->ops->productsymbolic)(mmabc->BC));
1302: PetscCall(MatProductCreate(A, mmabc->BC, NULL, &mmabc->ABC));
1303: PetscCall(MatSetOptionsPrefix(mmabc->ABC, prefix));
1304: PetscCall(MatAppendOptionsPrefix(mmabc->ABC, "P2_"));
1305: PetscCall(MatProductSetType(mmabc->ABC, p2));
1306: PetscCall(MatProductSetAlgorithm(mmabc->ABC, MATPRODUCTALGORITHMDEFAULT));
1307: PetscCall(MatProductSetFill(mmabc->ABC, product->fill));
1308: mmabc->ABC->product->api_user = product->api_user;
1309: PetscCall(MatProductSetFromOptions(mmabc->ABC));
1310: 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);
1311: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1312: mat->product = mmabc->ABC->product;
1313: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1314: /* use function pointer directly to prevent logging */
1315: PetscUseTypeMethod(mat, productsymbolic);
1316: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1317: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1318: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1319: mat->product = product;
1320: PetscFunctionReturn(PETSC_SUCCESS);
1321: }
1323: /*@
1324: MatProductGetType - Returns the type of matrix-matrix product associated with computing values for the given matrix
1326: Not Collective
1328: Input Parameter:
1329: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1331: Output Parameter:
1332: . mtype - the `MatProductType`
1334: Level: intermediate
1336: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductCreate()`, `MatProductType`, `MatProductAlgorithm`
1337: @*/
1338: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1339: {
1340: PetscFunctionBegin;
1342: PetscAssertPointer(mtype, 2);
1343: *mtype = MATPRODUCT_UNSPECIFIED;
1344: if (mat->product) *mtype = mat->product->type;
1345: PetscFunctionReturn(PETSC_SUCCESS);
1346: }
1348: /*@
1349: MatProductGetMats - Returns the matrices associated with the matrix-matrix product associated with computing values for the given matrix
1351: Not Collective
1353: Input Parameter:
1354: . mat - the matrix whose values are to be computed via a matrix-matrix product operation
1356: Output Parameters:
1357: + A - the first matrix
1358: . B - the second matrix
1359: - C - the third matrix (may be `NULL` for some `MatProductType`)
1361: Level: intermediate
1363: .seealso: [](ch_matrices), `MatProduct`, `Mat`, `MatProductCreateWithMat()`, `MatProductSetType()`, `MatProductSetAlgorithm()`, `MatProductCreate()`
1364: @*/
1365: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1366: {
1367: PetscFunctionBegin;
1369: if (A) *A = mat->product ? mat->product->A : NULL;
1370: if (B) *B = mat->product ? mat->product->B : NULL;
1371: if (C) *C = mat->product ? mat->product->C : NULL;
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }