Actual source code: matproduct.c
2: /*
3: Routines for matrix products. Calling procedure:
5: MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D)
6: MatProductSetType(D, MATPRODUCT_AB/AtB/ABt/PtAP/RARt/ABC)
7: MatProductSetAlgorithm(D, alg)
8: MatProductSetFill(D,fill)
9: MatProductSetFromOptions(D)
10: -> MatProductSetFromOptions_Private(D)
11: # Check matrix global sizes
12: if the matrices have the same setfromoptions routine, use it
13: if not, try:
14: -> Query MatProductSetFromOptions_Atype_Btype_Ctype_C(D) from A, B and C (in order)
15: if found -> run the specific setup that must set the symbolic operation (these callbacks should never fail)
16: if callback not found or no symbolic operation set
17: -> Query MatProductSetFromOptions_anytype_C(D) from A, B and C (in order) (e.g, matrices may have inner matrices like MATTRANSPOSEMAT)
18: if dispatch found but combination still not present do
19: -> check if B is dense and product type AtB or AB -> if true, basic looping of dense columns
20: -> check if triple product (PtAP, RARt or ABC) -> if true, set the Basic routines
22: # The setfromoptions calls MatProductSetFromOptions_Atype_Btype_Ctype should
23: # Check matrix local sizes for mpi matrices
24: # Set default algorithm
25: # Get runtime option
26: # Set D->ops->productsymbolic = MatProductSymbolic_productype_Atype_Btype_Ctype if found
28: MatProductSymbolic(D)
29: # Call MatProductSymbolic_productype_Atype_Btype_Ctype()
30: the callback must set the numeric phase D->ops->productnumeric = MatProductNumeric_productype_Atype_Btype_Ctype
32: MatProductNumeric(D)
33: # Call the numeric phase
35: # The symbolic phases are allowed to set extra data structures and attach those to the product
36: # this additional data can be reused between multiple numeric phases with the same matrices
37: # if not needed, call
38: MatProductClear(D)
39: */
41: #include <petsc/private/matimpl.h>
43: const char *const MatProductTypes[] = {"UNSPECIFIED","AB","AtB","ABt","PtAP","RARt","ABC"};
45: /* these are basic implementations relying on the old function pointers
46: * they are dangerous and should be removed in the future */
47: static PetscErrorCode MatProductNumeric_PtAP_Unsafe(Mat C)
48: {
49: Mat_Product *product = C->product;
50: Mat P = product->B,AP = product->Dwork;
52: /* AP = A*P */
53: MatProductNumeric(AP);
54: /* C = P^T*AP */
56: (*C->ops->transposematmultnumeric)(P,AP,C);
57: return 0;
58: }
60: static PetscErrorCode MatProductSymbolic_PtAP_Unsafe(Mat C)
61: {
62: Mat_Product *product = C->product;
63: Mat A=product->A,P=product->B,AP;
64: PetscReal fill=product->fill;
66: PetscInfo((PetscObject)C,"for A %s, P %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
67: /* AP = A*P */
68: MatProductCreate(A,P,NULL,&AP);
69: MatProductSetType(AP,MATPRODUCT_AB);
70: MatProductSetAlgorithm(AP,MATPRODUCTALGORITHMDEFAULT);
71: MatProductSetFill(AP,fill);
72: MatProductSetFromOptions(AP);
73: MatProductSymbolic(AP);
75: /* C = P^T*AP */
76: MatProductSetType(C,MATPRODUCT_AtB);
77: MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT);
78: product->A = P;
79: product->B = AP;
80: MatProductSetFromOptions(C);
81: 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: return 0;
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: /* RA = R*A */
98: MatProductNumeric(RA);
99: /* C = RA*R^T */
101: (*C->ops->mattransposemultnumeric)(RA,R,C);
102: return 0;
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: PetscInfo((PetscObject)C,"for A %s, R %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
112: /* RA = R*A */
113: MatProductCreate(R,A,NULL,&RA);
114: MatProductSetType(RA,MATPRODUCT_AB);
115: MatProductSetAlgorithm(RA,MATPRODUCTALGORITHMDEFAULT);
116: MatProductSetFill(RA,fill);
117: MatProductSetFromOptions(RA);
118: MatProductSymbolic(RA);
120: /* C = RA*R^T */
121: MatProductSetType(C,MATPRODUCT_ABt);
122: MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT);
123: product->A = RA;
124: MatProductSetFromOptions(C);
125: MatProductSymbolic(C);
127: /* resume user's original input matrix setting for A */
128: product->A = A;
129: product->Dwork = RA; /* save here so it will be destroyed with product C */
130: C->ops->productnumeric = MatProductNumeric_RARt_Unsafe;
131: return 0;
132: }
134: static PetscErrorCode MatProductNumeric_ABC_Unsafe(Mat mat)
135: {
136: Mat_Product *product = mat->product;
137: Mat A=product->A,BC=product->Dwork;
139: /* Numeric BC = B*C */
140: MatProductNumeric(BC);
141: /* Numeric mat = A*BC */
143: (*mat->ops->matmultnumeric)(A,BC,mat);
144: return 0;
145: }
147: static PetscErrorCode MatProductSymbolic_ABC_Unsafe(Mat mat)
148: {
149: Mat_Product *product = mat->product;
150: Mat B=product->B,C=product->C,BC;
151: PetscReal fill=product->fill;
153: PetscInfo((PetscObject)mat,"for A %s, B %s, C %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name,((PetscObject)product->C)->type_name);
154: /* Symbolic BC = B*C */
155: MatProductCreate(B,C,NULL,&BC);
156: MatProductSetType(BC,MATPRODUCT_AB);
157: MatProductSetAlgorithm(BC,MATPRODUCTALGORITHMDEFAULT);
158: MatProductSetFill(BC,fill);
159: MatProductSetFromOptions(BC);
160: MatProductSymbolic(BC);
162: /* Symbolic mat = A*BC */
163: MatProductSetType(mat,MATPRODUCT_AB);
164: MatProductSetAlgorithm(mat,MATPRODUCTALGORITHMDEFAULT);
165: product->B = BC;
166: product->Dwork = BC;
167: MatProductSetFromOptions(mat);
168: MatProductSymbolic(mat);
170: /* resume user's original input matrix setting for B */
171: product->B = B;
172: mat->ops->productnumeric = MatProductNumeric_ABC_Unsafe;
173: return 0;
174: }
176: static PetscErrorCode MatProductSymbolic_Unsafe(Mat mat)
177: {
178: Mat_Product *product = mat->product;
180: switch (product->type) {
181: case MATPRODUCT_PtAP:
182: MatProductSymbolic_PtAP_Unsafe(mat);
183: break;
184: case MATPRODUCT_RARt:
185: MatProductSymbolic_RARt_Unsafe(mat);
186: break;
187: case MATPRODUCT_ABC:
188: MatProductSymbolic_ABC_Unsafe(mat);
189: break;
190: default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]);
191: }
192: return 0;
193: }
195: /* ----------------------------------------------- */
196: /*@C
197: MatProductReplaceMats - Replace input matrices for a matrix product.
199: Collective on Mat
201: Input Parameters:
202: + A - the matrix or NULL if not being replaced
203: . B - the matrix or NULL if not being replaced
204: . C - the matrix or NULL if not being replaced
205: - D - the matrix product
207: Level: intermediate
209: Notes:
210: To reuse the symbolic phase, input matrices must have exactly the same data structure as the replaced one.
211: If the type of any of the input matrices is different than what was previously used, or their symmetry changed but
212: the symbolic phase took advantage of their symmetry, the product is cleared and MatProductSetFromOptions()/MatProductSymbolic() are invoked again.
214: .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
215: @*/
216: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
217: {
218: Mat_Product *product;
219: PetscBool flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;
222: MatCheckProduct(D,4);
223: product = D->product;
224: if (A) {
226: PetscObjectReference((PetscObject)A);
227: PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);
228: if (product->symbolic_used_the_fact_A_is_symmetric && !A->symmetric) { /* symbolic was built around a symmetric A, but the new A is not anymore */
229: flgA = PETSC_FALSE;
230: product->symbolic_used_the_fact_A_is_symmetric = PETSC_FALSE; /* reinit */
231: }
232: MatDestroy(&product->A);
233: product->A = A;
234: }
235: if (B) {
237: PetscObjectReference((PetscObject)B);
238: PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);
239: if (product->symbolic_used_the_fact_B_is_symmetric && !B->symmetric) {
240: flgB = PETSC_FALSE;
241: product->symbolic_used_the_fact_B_is_symmetric = PETSC_FALSE; /* reinit */
242: }
243: MatDestroy(&product->B);
244: product->B = B;
245: }
246: if (C) {
248: PetscObjectReference((PetscObject)C);
249: PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);
250: if (product->symbolic_used_the_fact_C_is_symmetric && !C->symmetric) {
251: flgC = PETSC_FALSE;
252: product->symbolic_used_the_fact_C_is_symmetric = PETSC_FALSE; /* reinit */
253: }
254: MatDestroy(&product->C);
255: product->C = C;
256: }
257: /* Any of the replaced mats is of a different type, reset */
258: if (!flgA || !flgB || !flgC) {
259: if (D->product->destroy) {
260: (*D->product->destroy)(D->product->data);
261: }
262: D->product->destroy = NULL;
263: D->product->data = NULL;
264: if (D->ops->productnumeric || D->ops->productsymbolic) {
265: MatProductSetFromOptions(D);
266: MatProductSymbolic(D);
267: }
268: }
269: return 0;
270: }
272: static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
273: {
274: Mat_Product *product = C->product;
275: Mat A = product->A, B = product->B;
276: PetscInt k, K = B->cmap->N;
277: PetscBool t = PETSC_TRUE,iscuda = PETSC_FALSE;
278: PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
279: char *Btype = NULL,*Ctype = NULL;
281: switch (product->type) {
282: case MATPRODUCT_AB:
283: t = PETSC_FALSE;
284: case MATPRODUCT_AtB:
285: break;
286: default: SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProductNumeric type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
287: }
288: if (PetscDefined(HAVE_CUDA)) {
289: VecType vtype;
291: MatGetVecType(A,&vtype);
292: PetscStrcmp(vtype,VECCUDA,&iscuda);
293: if (!iscuda) {
294: PetscStrcmp(vtype,VECSEQCUDA,&iscuda);
295: }
296: if (!iscuda) {
297: PetscStrcmp(vtype,VECMPICUDA,&iscuda);
298: }
299: if (iscuda) { /* Make sure we have up-to-date data on the GPU */
300: PetscStrallocpy(((PetscObject)B)->type_name,&Btype);
301: PetscStrallocpy(((PetscObject)C)->type_name,&Ctype);
302: MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B);
303: if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
304: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
306: }
307: MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C);
308: } else { /* Make sure we have up-to-date data on the CPU */
309: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
310: Bcpu = B->boundtocpu;
311: Ccpu = C->boundtocpu;
312: #endif
313: MatBindToCPU(B,PETSC_TRUE);
314: MatBindToCPU(C,PETSC_TRUE);
315: }
316: }
317: for (k=0;k<K;k++) {
318: Vec x,y;
320: MatDenseGetColumnVecRead(B,k,&x);
321: MatDenseGetColumnVecWrite(C,k,&y);
322: if (t) {
323: MatMultTranspose(A,x,y);
324: } else {
325: MatMult(A,x,y);
326: }
327: MatDenseRestoreColumnVecRead(B,k,&x);
328: MatDenseRestoreColumnVecWrite(C,k,&y);
329: }
330: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
331: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
332: if (PetscDefined(HAVE_CUDA)) {
333: if (iscuda) {
334: MatConvert(B,Btype,MAT_INPLACE_MATRIX,&B);
335: MatConvert(C,Ctype,MAT_INPLACE_MATRIX,&C);
336: } else {
337: MatBindToCPU(B,Bcpu);
338: MatBindToCPU(C,Ccpu);
339: }
340: }
341: PetscFree(Btype);
342: PetscFree(Ctype);
343: return 0;
344: }
346: static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
347: {
348: Mat_Product *product = C->product;
349: Mat A = product->A, B = product->B;
350: PetscBool isdense;
352: switch (product->type) {
353: case MATPRODUCT_AB:
354: MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
355: break;
356: case MATPRODUCT_AtB:
357: MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
358: break;
359: default: 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);
360: }
361: PetscObjectBaseTypeCompareAny((PetscObject)C,&isdense,MATSEQDENSE,MATMPIDENSE,"");
362: if (!isdense) {
363: MatSetType(C,((PetscObject)B)->type_name);
364: /* If matrix type of C was not set or not dense, we need to reset the pointer */
365: C->ops->productsymbolic = MatProductSymbolic_X_Dense;
366: }
367: C->ops->productnumeric = MatProductNumeric_X_Dense;
368: MatSetUp(C);
369: return 0;
370: }
372: /* a single driver to query the dispatching */
373: static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
374: {
375: Mat_Product *product = mat->product;
376: PetscInt Am,An,Bm,Bn,Cm,Cn;
377: Mat A = product->A,B = product->B,C = product->C;
378: const char* const Bnames[] = { "B", "R", "P" };
379: const char* bname;
380: PetscErrorCode (*fA)(Mat);
381: PetscErrorCode (*fB)(Mat);
382: PetscErrorCode (*fC)(Mat);
383: PetscErrorCode (*f)(Mat)=NULL;
385: mat->ops->productsymbolic = NULL;
386: mat->ops->productnumeric = NULL;
387: if (product->type == MATPRODUCT_UNSPECIFIED) return 0;
391: if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
392: if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
393: else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
394: else bname = Bnames[0];
396: /* Check matrices sizes */
397: Am = A->rmap->N;
398: An = A->cmap->N;
399: Bm = B->rmap->N;
400: Bn = B->cmap->N;
401: Cm = C ? C->rmap->N : 0;
402: Cn = C ? C->cmap->N : 0;
403: if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
404: if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
408: fA = A->ops->productsetfromoptions;
409: fB = B->ops->productsetfromoptions;
410: fC = C ? C->ops->productsetfromoptions : fA;
411: if (C) {
412: 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);
413: } else {
414: PetscInfo(mat,"MatProductType %s for A %s, %s %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name);
415: }
416: if (fA == fB && fA == fC && fA) {
417: PetscInfo(mat," matching op\n");
418: (*fA)(mat);
419: }
420: /* We may have found f but it did not succeed */
421: if (!mat->ops->productsymbolic) { /* query MatProductSetFromOptions_Atype_Btype_Ctype */
422: char mtypes[256];
423: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
424: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
425: PetscStrlcat(mtypes,"_",sizeof(mtypes));
426: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
427: if (C) {
428: PetscStrlcat(mtypes,"_",sizeof(mtypes));
429: PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));
430: }
431: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
433: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
434: PetscInfo(mat," querying %s from A? %p\n",mtypes,f);
435: if (!f) {
436: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
437: PetscInfo(mat," querying %s from %s? %p\n",mtypes,bname,f);
438: }
439: if (!f && C) {
440: PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
441: PetscInfo(mat," querying %s from C? %p\n",mtypes,f);
442: }
443: if (f) (*f)(mat);
445: /* We may have found f but it did not succeed */
446: /* some matrices (i.e. MATTRANSPOSE, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
447: if (!mat->ops->productsymbolic) {
448: PetscStrncpy(mtypes,"MatProductSetFromOptions_anytype_C",sizeof(mtypes));
449: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
450: PetscInfo(mat," querying %s from A? %p\n",mtypes,f);
451: if (!f) {
452: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
453: PetscInfo(mat," querying %s from %s? %p\n",mtypes,bname,f);
454: }
455: if (!f && C) {
456: PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
457: PetscInfo(mat," querying %s from C? %p\n",mtypes,f);
458: }
459: }
460: if (f) (*f)(mat);
461: }
463: /* We may have found f but it did not succeed */
464: if (!mat->ops->productsymbolic) {
465: /* we can still compute the product if B is of type dense */
466: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
467: PetscBool isdense;
469: PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");
470: if (isdense) {
472: mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
473: PetscInfo(mat," using basic looping over columns of a dense matrix\n");
474: }
475: } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Unsafe() for triple products only */
476: /*
477: TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
478: the combination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
479: before computing the symbolic phase
480: */
481: PetscInfo(mat," symbolic product not supported, using MatProductSymbolic_Unsafe() implementation\n");
482: mat->ops->productsymbolic = MatProductSymbolic_Unsafe;
483: }
484: }
485: if (!mat->ops->productsymbolic) {
486: PetscInfo(mat," symbolic product is not supported\n");
487: }
488: return 0;
489: }
491: /*@C
492: MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.
494: Logically Collective on Mat
496: Input Parameter:
497: . mat - the matrix
499: Options Database Keys:
500: . -mat_product_clear - Clear intermediate data structures after MatProductNumeric() has been called
502: Level: intermediate
504: .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
505: @*/
506: PetscErrorCode MatProductSetFromOptions(Mat mat)
507: {
511: MatCheckProduct(mat,1);
513: PetscObjectOptionsBegin((PetscObject)mat);
514: PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);
515: PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");
516: PetscOptionsEnd();
517: MatProductSetFromOptions_Private(mat);
519: return 0;
520: }
522: /*@C
523: MatProductView - View a MatProduct
525: Logically Collective on Mat
527: Input Parameter:
528: . mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()
530: Level: intermediate
532: .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
533: @*/
534: PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
535: {
537: if (!mat->product) return 0;
538: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);
541: if (mat->product->view) {
542: (*mat->product->view)(mat,viewer);
543: }
544: return 0;
545: }
547: /* ----------------------------------------------- */
548: /* these are basic implementations relying on the old function pointers
549: * they are dangerous and should be removed in the future */
550: PetscErrorCode MatProductNumeric_AB(Mat mat)
551: {
552: Mat_Product *product = mat->product;
553: Mat A=product->A,B=product->B;
556: (*mat->ops->matmultnumeric)(A,B,mat);
557: return 0;
558: }
560: PetscErrorCode MatProductNumeric_AtB(Mat mat)
561: {
562: Mat_Product *product = mat->product;
563: Mat A=product->A,B=product->B;
566: (*mat->ops->transposematmultnumeric)(A,B,mat);
567: return 0;
568: }
570: PetscErrorCode MatProductNumeric_ABt(Mat mat)
571: {
572: Mat_Product *product = mat->product;
573: Mat A=product->A,B=product->B;
576: (*mat->ops->mattransposemultnumeric)(A,B,mat);
577: return 0;
578: }
580: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
581: {
582: Mat_Product *product = mat->product;
583: Mat A=product->A,B=product->B;
586: (*mat->ops->ptapnumeric)(A,B,mat);
587: return 0;
588: }
590: PetscErrorCode MatProductNumeric_RARt(Mat mat)
591: {
592: Mat_Product *product = mat->product;
593: Mat A=product->A,B=product->B;
596: (*mat->ops->rartnumeric)(A,B,mat);
597: return 0;
598: }
600: PetscErrorCode MatProductNumeric_ABC(Mat mat)
601: {
602: Mat_Product *product = mat->product;
603: Mat A=product->A,B=product->B,C=product->C;
606: (*mat->ops->matmatmultnumeric)(A,B,C,mat);
607: return 0;
608: }
610: /* ----------------------------------------------- */
612: /*@
613: MatProductNumeric - Implement a matrix product with numerical values.
615: Collective on Mat
617: Input/Output Parameter:
618: . mat - the matrix holding the product
620: Level: intermediate
622: Notes: MatProductSymbolic() must have been called on mat before calling this function
624: .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
625: @*/
626: PetscErrorCode MatProductNumeric(Mat mat)
627: {
628: PetscLogEvent eventtype = -1;
629: PetscBool missing = PETSC_FALSE;
632: MatCheckProduct(mat,1);
633: switch (mat->product->type) {
634: case MATPRODUCT_AB:
635: eventtype = MAT_MatMultNumeric;
636: break;
637: case MATPRODUCT_AtB:
638: eventtype = MAT_TransposeMatMultNumeric;
639: break;
640: case MATPRODUCT_ABt:
641: eventtype = MAT_MatTransposeMultNumeric;
642: break;
643: case MATPRODUCT_PtAP:
644: eventtype = MAT_PtAPNumeric;
645: break;
646: case MATPRODUCT_RARt:
647: eventtype = MAT_RARtNumeric;
648: break;
649: case MATPRODUCT_ABC:
650: eventtype = MAT_MatMatMultNumeric;
651: break;
652: default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
653: }
655: if (mat->ops->productnumeric) {
656: PetscLogEventBegin(eventtype,mat,0,0,0);
657: (*mat->ops->productnumeric)(mat);
658: PetscLogEventEnd(eventtype,mat,0,0,0);
659: } else missing = PETSC_TRUE;
661: if (missing || !mat->product) {
662: char errstr[256];
664: if (mat->product->type == MATPRODUCT_ABC) {
665: 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);
666: } else {
667: 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);
668: }
671: }
673: if (mat->product->clear) {
674: MatProductClear(mat);
675: }
676: PetscObjectStateIncrease((PetscObject)mat);
677: return 0;
678: }
680: /* ----------------------------------------------- */
681: /* these are basic implementations relying on the old function pointers
682: * they are dangerous and should be removed in the future */
683: PetscErrorCode MatProductSymbolic_AB(Mat mat)
684: {
685: Mat_Product *product = mat->product;
686: Mat A=product->A,B=product->B;
689: (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);
690: mat->ops->productnumeric = MatProductNumeric_AB;
691: return 0;
692: }
694: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
695: {
696: Mat_Product *product = mat->product;
697: Mat A=product->A,B=product->B;
700: (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
701: mat->ops->productnumeric = MatProductNumeric_AtB;
702: return 0;
703: }
705: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
706: {
707: Mat_Product *product = mat->product;
708: Mat A=product->A,B=product->B;
711: (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
712: mat->ops->productnumeric = MatProductNumeric_ABt;
713: return 0;
714: }
716: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
717: {
718: Mat_Product *product = mat->product;
719: Mat A=product->A,B=product->B,C=product->C;
722: (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
723: mat->ops->productnumeric = MatProductNumeric_ABC;
724: return 0;
725: }
727: /* ----------------------------------------------- */
729: /*@
730: MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.
732: Collective on Mat
734: Input/Output Parameter:
735: . mat - the matrix to hold a product
737: Level: intermediate
739: Notes: MatProductSetFromOptions() must have been called on mat before calling this function
741: .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
742: @*/
743: PetscErrorCode MatProductSymbolic(Mat mat)
744: {
745: PetscLogEvent eventtype = -1;
746: PetscBool missing = PETSC_FALSE;
749: MatCheckProduct(mat,1);
751: switch (mat->product->type) {
752: case MATPRODUCT_AB:
753: eventtype = MAT_MatMultSymbolic;
754: break;
755: case MATPRODUCT_AtB:
756: eventtype = MAT_TransposeMatMultSymbolic;
757: break;
758: case MATPRODUCT_ABt:
759: eventtype = MAT_MatTransposeMultSymbolic;
760: break;
761: case MATPRODUCT_PtAP:
762: eventtype = MAT_PtAPSymbolic;
763: break;
764: case MATPRODUCT_RARt:
765: eventtype = MAT_RARtSymbolic;
766: break;
767: case MATPRODUCT_ABC:
768: eventtype = MAT_MatMatMultSymbolic;
769: break;
770: default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
771: }
772: mat->ops->productnumeric = NULL;
773: if (mat->ops->productsymbolic) {
774: PetscLogEventBegin(eventtype,mat,0,0,0);
775: (*mat->ops->productsymbolic)(mat);
776: PetscLogEventEnd(eventtype,mat,0,0,0);
777: } else missing = PETSC_TRUE;
779: if (missing || !mat->product || !mat->ops->productnumeric) {
780: char errstr[256];
782: if (mat->product->type == MATPRODUCT_ABC) {
783: 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);
784: } else {
785: 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);
786: }
790: }
791: return 0;
792: }
794: /*@
795: MatProductSetFill - Set an expected fill of the matrix product.
797: Collective on Mat
799: Input Parameters:
800: + mat - the matrix product
801: - fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use PETSC_DEFAULT if you do not have a good estimate. If the product is a dense matrix, this is irrelevant.
803: Level: intermediate
805: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
806: @*/
807: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
808: {
810: MatCheckProduct(mat,1);
811: if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
812: else mat->product->fill = fill;
813: return 0;
814: }
816: /*@
817: MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.
819: Collective on Mat
821: Input Parameters:
822: + mat - the matrix product
823: - alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHMDEFAULT.
825: Options Database Key:
826: . -mat_product_algorithm <algorithm> - Sets the algorithm; use -help for a list
827: of available algorithms (for instance, scalable, outerproduct, etc.)
829: Level: intermediate
831: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
832: @*/
833: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
834: {
836: MatCheckProduct(mat,1);
837: PetscFree(mat->product->alg);
838: PetscStrallocpy(alg,&mat->product->alg);
839: return 0;
840: }
842: /*@
843: MatProductSetType - Sets a particular matrix product type
845: Collective on Mat
847: Input Parameters:
848: + mat - the matrix
849: - productype - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.
851: Level: intermediate
853: .seealso: MatProductCreate(), MatProductType
854: @*/
855: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
856: {
858: MatCheckProduct(mat,1);
860: if (productype != mat->product->type) {
861: if (mat->product->destroy) {
862: (*mat->product->destroy)(mat->product->data);
863: }
864: mat->product->destroy = NULL;
865: mat->product->data = NULL;
866: mat->ops->productsymbolic = NULL;
867: mat->ops->productnumeric = NULL;
868: }
869: mat->product->type = productype;
870: return 0;
871: }
873: /*@
874: MatProductClear - Clears matrix product internal structure.
876: Collective on Mat
878: Input Parameters:
879: . mat - the product matrix
881: Level: intermediate
883: Notes: this function should be called to remove any intermediate data used by the product
884: After having called this function, MatProduct operations can no longer be used on mat
885: @*/
886: PetscErrorCode MatProductClear(Mat mat)
887: {
888: Mat_Product *product = mat->product;
891: if (product) {
892: MatDestroy(&product->A);
893: MatDestroy(&product->B);
894: MatDestroy(&product->C);
895: PetscFree(product->alg);
896: MatDestroy(&product->Dwork);
897: if (product->destroy) {
898: (*product->destroy)(product->data);
899: }
900: }
901: PetscFree(mat->product);
902: mat->ops->productsymbolic = NULL;
903: mat->ops->productnumeric = NULL;
904: return 0;
905: }
907: /* Create a supporting struct and attach it to the matrix product */
908: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
909: {
910: Mat_Product *product=NULL;
914: PetscNewLog(D,&product);
915: product->A = A;
916: product->B = B;
917: product->C = C;
918: product->type = MATPRODUCT_UNSPECIFIED;
919: product->Dwork = NULL;
920: product->api_user = PETSC_FALSE;
921: product->clear = PETSC_FALSE;
922: D->product = product;
924: MatProductSetAlgorithm(D,MATPRODUCTALGORITHMDEFAULT);
925: MatProductSetFill(D,PETSC_DEFAULT);
927: PetscObjectReference((PetscObject)A);
928: PetscObjectReference((PetscObject)B);
929: PetscObjectReference((PetscObject)C);
930: return 0;
931: }
933: /*@
934: MatProductCreateWithMat - Setup a given matrix as a matrix product.
936: Collective on Mat
938: Input Parameters:
939: + A - the first matrix
940: . B - the second matrix
941: . C - the third matrix (optional)
942: - D - the matrix which will be used as a product
944: Output Parameters:
945: . D - the product matrix
947: Notes:
948: Any product data attached to D will be cleared
950: Level: intermediate
952: .seealso: MatProductCreate(), MatProductClear()
953: @*/
954: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
955: {
958: MatCheckPreallocated(A,1);
964: MatCheckPreallocated(B,2);
968: if (C) {
971: MatCheckPreallocated(C,3);
974: }
978: MatCheckPreallocated(D,4);
982: /* Create a supporting struct and attach it to D */
983: MatProductClear(D);
984: MatProductCreate_Private(A,B,C,D);
985: return 0;
986: }
988: /*@
989: MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.
991: Collective on Mat
993: Input Parameters:
994: + A - the first matrix
995: . B - the second matrix
996: - C - the third matrix (optional)
998: Output Parameters:
999: . D - the product matrix
1001: Level: intermediate
1003: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1004: @*/
1005: PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1006: {
1014: if (C) {
1018: }
1021: MatCreate(PetscObjectComm((PetscObject)A),D);
1022: MatProductCreate_Private(A,B,C,*D);
1023: return 0;
1024: }
1026: /*
1027: These are safe basic implementations of ABC, RARt and PtAP
1028: that do not rely on mat->ops->matmatop function pointers.
1029: They only use the MatProduct API and are currently used by
1030: cuSPARSE and KOKKOS-KERNELS backends
1031: */
1032: typedef struct {
1033: Mat BC;
1034: Mat ABC;
1035: } MatMatMatPrivate;
1037: static PetscErrorCode MatDestroy_MatMatMatPrivate(void *data)
1038: {
1039: MatMatMatPrivate *mmdata = (MatMatMatPrivate *)data;
1041: MatDestroy(&mmdata->BC);
1042: MatDestroy(&mmdata->ABC);
1043: PetscFree(data);
1044: return 0;
1045: }
1047: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
1048: {
1049: Mat_Product *product = mat->product;
1050: MatMatMatPrivate *mmabc;
1052: MatCheckProduct(mat,1);
1054: mmabc = (MatMatMatPrivate *)mat->product->data;
1056: /* use function pointer directly to prevent logging */
1057: (*mmabc->BC->ops->productnumeric)(mmabc->BC);
1058: /* swap ABC product stuff with that of ABC for the numeric phase on mat */
1059: mat->product = mmabc->ABC->product;
1060: mat->ops->productnumeric = mmabc->ABC->ops->productnumeric;
1062: /* use function pointer directly to prevent logging */
1063: (*mat->ops->productnumeric)(mat);
1064: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1065: mat->product = product;
1066: return 0;
1067: }
1069: PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
1070: {
1071: Mat_Product *product = mat->product;
1072: Mat A, B ,C;
1073: MatProductType p1,p2;
1074: MatMatMatPrivate *mmabc;
1075: const char *prefix;
1077: MatCheckProduct(mat,1);
1079: MatGetOptionsPrefix(mat,&prefix);
1080: PetscNew(&mmabc);
1081: product->data = mmabc;
1082: product->destroy = MatDestroy_MatMatMatPrivate;
1083: switch (product->type) {
1084: case MATPRODUCT_PtAP:
1085: p1 = MATPRODUCT_AB;
1086: p2 = MATPRODUCT_AtB;
1087: A = product->B;
1088: B = product->A;
1089: C = product->B;
1090: break;
1091: case MATPRODUCT_RARt:
1092: p1 = MATPRODUCT_ABt;
1093: p2 = MATPRODUCT_AB;
1094: A = product->B;
1095: B = product->A;
1096: C = product->B;
1097: break;
1098: case MATPRODUCT_ABC:
1099: p1 = MATPRODUCT_AB;
1100: p2 = MATPRODUCT_AB;
1101: A = product->A;
1102: B = product->B;
1103: C = product->C;
1104: break;
1105: default:
1106: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Not for ProductType %s",MatProductTypes[product->type]);
1107: }
1108: MatProductCreate(B,C,NULL,&mmabc->BC);
1109: MatSetOptionsPrefix(mmabc->BC,prefix);
1110: MatAppendOptionsPrefix(mmabc->BC,"P1_");
1111: MatProductSetType(mmabc->BC,p1);
1112: MatProductSetAlgorithm(mmabc->BC,MATPRODUCTALGORITHMDEFAULT);
1113: MatProductSetFill(mmabc->BC,product->fill);
1114: mmabc->BC->product->api_user = product->api_user;
1115: MatProductSetFromOptions(mmabc->BC);
1117: /* use function pointer directly to prevent logging */
1118: (*mmabc->BC->ops->productsymbolic)(mmabc->BC);
1120: MatProductCreate(A,mmabc->BC,NULL,&mmabc->ABC);
1121: MatSetOptionsPrefix(mmabc->ABC,prefix);
1122: MatAppendOptionsPrefix(mmabc->ABC,"P2_");
1123: MatProductSetType(mmabc->ABC,p2);
1124: MatProductSetAlgorithm(mmabc->ABC,MATPRODUCTALGORITHMDEFAULT);
1125: MatProductSetFill(mmabc->ABC,product->fill);
1126: mmabc->ABC->product->api_user = product->api_user;
1127: MatProductSetFromOptions(mmabc->ABC);
1129: /* swap ABC product stuff with that of ABC for the symbolic phase on mat */
1130: mat->product = mmabc->ABC->product;
1131: mat->ops->productsymbolic = mmabc->ABC->ops->productsymbolic;
1132: /* use function pointer directly to prevent logging */
1133: (*mat->ops->productsymbolic)(mat);
1134: mmabc->ABC->ops->productnumeric = mat->ops->productnumeric;
1135: mat->ops->productsymbolic = MatProductSymbolic_ABC_Basic;
1136: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
1137: mat->product = product;
1138: return 0;
1139: }
1141: /*@
1142: MatProductGetType - Returns the type of the MatProduct.
1144: Not collective
1146: Input Parameter:
1147: . mat - the matrix
1149: Output Parameter:
1150: . mtype - the MatProduct type
1152: Level: intermediate
1154: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductCreate()
1155: @*/
1156: PetscErrorCode MatProductGetType(Mat mat, MatProductType *mtype)
1157: {
1160: *mtype = MATPRODUCT_UNSPECIFIED;
1161: if (mat->product) *mtype = mat->product->type;
1162: return 0;
1163: }
1165: /*@
1166: MatProductGetMats - Returns the matrices associated with the MatProduct.
1168: Not collective
1170: Input Parameter:
1171: . mat - the product matrix
1173: Output Parameters:
1174: + A - the first matrix
1175: . B - the second matrix
1176: - C - the third matrix (optional)
1178: Level: intermediate
1180: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1181: @*/
1182: PetscErrorCode MatProductGetMats(Mat mat, Mat *A, Mat *B, Mat *C)
1183: {
1185: if (A) *A = mat->product ? mat->product->A : NULL;
1186: if (B) *B = mat->product ? mat->product->B : NULL;
1187: if (C) *C = mat->product ? mat->product->C : NULL;
1188: return 0;
1189: }