Actual source code: matproduct.c
petsc-3.13.6 2020-09-29
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_producttype(D):
11: # Check matrix global sizes
12: -> MatProductSetFromOptions_Atype_Btype_Ctype(D);
13: ->MatProductSetFromOptions_Atype_Btype_Ctype_productype(D):
14: # Check matrix local sizes for mpi matrices
15: # Set default algorithm
16: # Get runtime option
17: # Set D->ops->productsymbolic = MatProductSymbolic_productype_Atype_Btype_Ctype;
19: PetscLogEventBegin()
20: MatProductSymbolic(D):
21: # Call MatxxxSymbolic_Atype_Btype_Ctype();
22: # Set D->ops->productnumeric = MatProductNumeric_productype_Atype_Btype_Ctype;
23: PetscLogEventEnd()
25: PetscLogEventBegin()
26: MatProductNumeric(D);
27: # Call (D->ops->matxxxnumeric)();
28: PetscLogEventEnd()
29: */
31: #include <petsc/private/matimpl.h>
33: const char *const MatProductTypes[] = {"AB","AtB","ABt","PtAP","RARt","ABC","MatProductType","MAT_Product_",0};
35: static PetscErrorCode MatProductNumeric_PtAP_Basic(Mat C)
36: {
38: Mat_Product *product = C->product;
39: Mat P = product->B,AP = product->Dwork;
42: /* AP = A*P */
43: MatProductNumeric(AP);
44: /* C = P^T*AP */
45: (C->ops->transposematmultnumeric)(P,AP,C);
46: return(0);
47: }
49: static PetscErrorCode MatProductSymbolic_PtAP_Basic(Mat C)
50: {
52: Mat_Product *product = C->product;
53: Mat A=product->A,P=product->B,AP;
54: PetscReal fill=product->fill;
57: /* AP = A*P */
58: MatProductCreate(A,P,NULL,&AP);
59: MatProductSetType(AP,MATPRODUCT_AB);
60: MatProductSetAlgorithm(AP,"default");
61: MatProductSetFill(AP,fill);
62: MatProductSetFromOptions(AP);
63: MatProductSymbolic(AP);
65: /* C = P^T*AP */
66: MatProductSetType(C,MATPRODUCT_AtB);
67: product->alg = "default";
68: product->A = P;
69: product->B = AP;
70: MatProductSetFromOptions(C);
71: MatProductSymbolic(C);
73: /* resume user's original input matrix setting for A and B */
74: product->A = A;
75: product->B = P;
76: product->Dwork = AP;
78: C->ops->productnumeric = MatProductNumeric_PtAP_Basic;
79: return(0);
80: }
82: static PetscErrorCode MatProductNumeric_RARt_Basic(Mat C)
83: {
85: Mat_Product *product = C->product;
86: Mat R=product->B,RA=product->Dwork;
89: /* RA = R*A */
90: MatProductNumeric(RA);
91: /* C = RA*R^T */
92: (C->ops->mattransposemultnumeric)(RA,R,C);
93: return(0);
94: }
96: static PetscErrorCode MatProductSymbolic_RARt_Basic(Mat C)
97: {
99: Mat_Product *product = C->product;
100: Mat A=product->A,R=product->B,RA;
101: PetscReal fill=product->fill;
104: /* RA = R*A */
105: MatProductCreate(R,A,NULL,&RA);
106: MatProductSetType(RA,MATPRODUCT_AB);
107: MatProductSetAlgorithm(RA,"default");
108: MatProductSetFill(RA,fill);
109: MatProductSetFromOptions(RA);
110: MatProductSymbolic(RA);
112: /* C = RA*R^T */
113: MatProductSetType(C,MATPRODUCT_ABt);
114: product->alg = "default";
115: product->A = RA;
116: MatProductSetFromOptions(C);
117: MatProductSymbolic(C);
119: /* resume user's original input matrix setting for A */
120: product->A = A;
121: product->Dwork = RA; /* save here so it will be destroyed with product C */
122: C->ops->productnumeric = MatProductNumeric_RARt_Basic;
123: return(0);
124: }
126: static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
127: {
129: Mat_Product *product = mat->product;
130: Mat A=product->A,BC=product->Dwork;
133: /* Numeric BC = B*C */
134: MatProductNumeric(BC);
135: /* Numeric mat = A*BC */
136: (mat->ops->matmultnumeric)(A,BC,mat);
137: return(0);
138: }
140: static PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
141: {
143: Mat_Product *product = mat->product;
144: Mat B=product->B,C=product->C,BC;
145: PetscReal fill=product->fill;
148: /* Symbolic BC = B*C */
149: MatProductCreate(B,C,NULL,&BC);
150: MatProductSetType(BC,MATPRODUCT_AB);
151: MatProductSetAlgorithm(BC,"default");
152: MatProductSetFill(BC,fill);
153: MatProductSetFromOptions(BC);
154: MatProductSymbolic(BC);
156: /* Symbolic mat = A*BC */
157: MatProductSetType(mat,MATPRODUCT_AB);
158: product->alg = "default";
159: product->B = BC;
160: product->Dwork = BC;
161: MatProductSetFromOptions(mat);
162: MatProductSymbolic(mat);
164: /* resume user's original input matrix setting for B */
165: product->B = B;
166: mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
167: return(0);
168: }
170: PetscErrorCode MatProductSymbolic_Basic(Mat mat)
171: {
173: Mat_Product *product = mat->product;
176: switch (product->type) {
177: case MATPRODUCT_PtAP:
178: PetscInfo2((PetscObject)mat, "MatProduct_Basic_PtAP() for A %s, P %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
179: MatProductSymbolic_PtAP_Basic(mat);
180: break;
181: case MATPRODUCT_RARt:
182: PetscInfo2((PetscObject)mat, "MatProduct_Basic_RARt() for A %s, R %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);
183: MatProductSymbolic_RARt_Basic(mat);
184: break;
185: case MATPRODUCT_ABC:
186: PetscInfo3((PetscObject)mat, "MatProduct_Basic_ABC() for A %s, B %s, C %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name,((PetscObject)product->C)->type_name);
187: MatProductSymbolic_ABC_Basic(mat);
188: break;
189: default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType is not supported");
190: }
191: return(0);
192: }
194: /* ----------------------------------------------- */
195: /*@C
196: MatProductReplaceMats - Replace input matrices for a matrix product.
198: Collective on Mat
200: Input Parameters:
201: + A - the matrix or NULL if not being replaced
202: . B - the matrix or NULL if not being replaced
203: . C - the matrix or NULL if not being replaced
204: - D - the matrix product
206: Level: intermediate
208: Notes:
209: Input matrix must have exactly same data structure as replaced one.
211: .seealso: MatProductCreate()
212: @*/
213: PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
214: {
216: Mat_Product *product=D->product;
219: if (!product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_NULL,"Mat D does not have struct 'product'. Call MatProductReplaceProduct(). \n");
220: if (A) {
221: if (!product->Areplaced) {
222: PetscObjectReference((PetscObject)A); /* take ownership of input */
223: MatDestroy(&product->A); /* release old reference */
224: product->A = A;
225: } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"Matrix A was changed by a PETSc internal routine, cannot be replaced");
226: }
227: if (B) {
228: if (!product->Breplaced) {
229: PetscObjectReference((PetscObject)B); /* take ownership of input */
230: MatDestroy(&product->B); /* release old reference */
231: product->B = B;
232: } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"Matrix B was changed by a PETSc internal routine, cannot be replaced");
233: }
234: if (C) {
235: PetscObjectReference((PetscObject)C); /* take ownership of input */
236: MatDestroy(&product->C); /* release old reference */
237: product->C = C;
238: }
239: return(0);
240: }
242: /* ----------------------------------------------- */
243: static PetscErrorCode MatProductSetFromOptions_AB(Mat mat)
244: {
246: Mat_Product *product = mat->product;
247: Mat A=product->A,B=product->B;
248: PetscBool sametype;
249: PetscErrorCode (*fA)(Mat);
250: PetscErrorCode (*fB)(Mat);
251: PetscErrorCode (*f)(Mat)=NULL;
252: PetscBool A_istrans,B_istrans;
255: /* Check matrix global sizes */
256: if (B->rmap->N!=A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N);
258: fA = A->ops->productsetfromoptions;
259: fB = B->ops->productsetfromoptions;
261: PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);
262: PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&A_istrans);
263: PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&B_istrans);
264: PetscInfo2(mat,"for A %s, B %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
266: if (fB == fA && sametype && (!A_istrans || !B_istrans)) {
267: PetscInfo(mat," sametype and matching op\n");
268: f = fB;
269: } else {
270: char mtypes[256];
271: PetscBool At_istrans=PETSC_TRUE,Bt_istrans=PETSC_TRUE;
272: Mat At = NULL,Bt = NULL;
274: if (A_istrans && !B_istrans) {
275: MatTransposeGetMat(A,&At);
276: PetscObjectTypeCompare((PetscObject)At,MATTRANSPOSEMAT,&At_istrans);
277: if (At_istrans) { /* mat = ATT * B */
278: Mat Att = NULL;
279: MatTransposeGetMat(At,&Att);
280: PetscObjectReference((PetscObject)Att);
281: MatDestroy(&product->A);
282: A = Att;
283: product->A = Att; /* use Att for matproduct */
284: product->Areplaced = PETSC_TRUE; /* Att = A, but has native matrix type */
285: } else { /* !At_istrans: mat = At^T*B */
286: PetscObjectReference((PetscObject)At);
287: MatDestroy(&product->A);
288: A = At;
289: product->A = At;
290: product->Areplaced = PETSC_TRUE;
291: product->type = MATPRODUCT_AtB;
292: }
293: } else if (!A_istrans && B_istrans) {
294: MatTransposeGetMat(B,&Bt);
295: PetscObjectTypeCompare((PetscObject)Bt,MATTRANSPOSEMAT,&Bt_istrans);
296: if (Bt_istrans) { /* mat = A * BTT */
297: Mat Btt = NULL;
298: MatTransposeGetMat(Bt,&Btt);
299: PetscObjectReference((PetscObject)Btt);
300: MatDestroy(&product->B);
301: B = Btt;
302: product->B = Btt; /* use Btt for matproduct */
303: product->Breplaced = PETSC_TRUE;
304: } else { /* !Bt_istrans */
305: /* mat = A*Bt^T */
306: PetscObjectReference((PetscObject)Bt);
307: MatDestroy(&product->B);
308: B = Bt;
309: product->B = Bt;
310: product->Breplaced = PETSC_TRUE;
311: product->type = MATPRODUCT_ABt;
312: }
313: } else if (A_istrans && B_istrans) { /* mat = At^T * Bt^T */
314: MatTransposeGetMat(A,&At);
315: PetscObjectTypeCompare((PetscObject)At,MATTRANSPOSEMAT,&At_istrans);
316: MatTransposeGetMat(B,&Bt);
317: PetscObjectTypeCompare((PetscObject)Bt,MATTRANSPOSEMAT,&Bt_istrans);
318: if (At_istrans && Bt_istrans) {
319: Mat Att= NULL,Btt = NULL;
320: MatTransposeGetMat(At,&Att);
321: MatTransposeGetMat(Bt,&Btt);
322: PetscObjectReference((PetscObject)Att);
323: PetscObjectReference((PetscObject)Btt);
324: MatDestroy(&product->A);
325: MatDestroy(&product->B);
326: A = Att;
327: product->A = Att; product->Areplaced = PETSC_TRUE;
328: B = Btt;
329: product->B = Btt; product->Breplaced = PETSC_TRUE;
330: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Not supported yet");
331: }
333: /* query MatProductSetFromOptions_Atype_Btype */
334: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
335: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
336: PetscStrlcat(mtypes,"_",sizeof(mtypes));
337: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
338: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
339: PetscInfo1(mat," querying %s\n",f);
340: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
341: if (!f) {
342: PetscInfo1(mat," querying %s\n",f);
343: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
344: }
345: }
347: if (!f) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
348: (*f)(mat);
349: return(0);
350: }
352: static PetscErrorCode MatProductSetFromOptions_AtB(Mat mat)
353: {
355: Mat_Product *product = mat->product;
356: Mat A=product->A,B=product->B;
357: PetscBool sametype;
358: PetscErrorCode (*fA)(Mat);
359: PetscErrorCode (*fB)(Mat);
360: PetscErrorCode (*f)(Mat)=NULL;
363: /* Check matrix global sizes */
364: if (B->rmap->N!=A->rmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->rmap->N);
366: fA = A->ops->productsetfromoptions;
367: fB = B->ops->productsetfromoptions;
369: PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);
370: PetscInfo2(mat,"for A %s, B %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
372: if (fB == fA && sametype) {
373: PetscInfo(mat," sametype and matching op\n");
374: f = fB;
375: } else {
376: char mtypes[256];
377: PetscBool istrans;
378: PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&istrans);
379: if (!istrans) {
380: /* query MatProductSetFromOptions_Atype_Btype */
381: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
382: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
383: PetscStrlcat(mtypes,"_",sizeof(mtypes));
384: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
385: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
386: PetscInfo1(mat," querying %s\n",f);
387: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
388: } else {
389: Mat T = NULL;
390: SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AtB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
392: MatTransposeGetMat(A,&T);
393: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
394: PetscStrlcat(mtypes,((PetscObject)T)->type_name,sizeof(mtypes));
395: PetscStrlcat(mtypes,"_",sizeof(mtypes));
396: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
397: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
399: product->type = MATPRODUCT_AtB;
400: PetscInfo1(mat," querying %s\n",f);
401: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
402: }
404: if (!f) {
405: PetscInfo1(mat," querying %s\n",f);
406: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
407: }
408: }
409: if (!f) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
411: (*f)(mat);
412: return(0);
413: }
415: static PetscErrorCode MatProductSetFromOptions_ABt(Mat mat)
416: {
418: Mat_Product *product = mat->product;
419: Mat A=product->A,B=product->B;
420: PetscBool sametype;
421: PetscErrorCode (*fA)(Mat);
422: PetscErrorCode (*fB)(Mat);
423: PetscErrorCode (*f)(Mat)=NULL;
426: /* Check matrix global sizes */
427: if (B->cmap->N!=A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, AN %D != BN %D",A->cmap->N,B->cmap->N);
429: fA = A->ops->productsetfromoptions;
430: fB = B->ops->productsetfromoptions;
432: PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);
433: PetscInfo2(mat,"for A %s, B %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
434: if (fB == fA && sametype) {
435: PetscInfo(mat," sametype and matching op\n");
436: f = fB;
437: } else {
438: char mtypes[256];
439: PetscBool istrans;
440: PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&istrans);
441: if (!istrans) {
442: /* query MatProductSetFromOptions_Atype_Btype */
443: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
444: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
445: PetscStrlcat(mtypes,"_",sizeof(mtypes));
446: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
447: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
448: PetscInfo1(mat," querying %s\n",f);
449: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
450: } else {
451: Mat T = NULL;
452: SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_ABt for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
454: MatTransposeGetMat(A,&T);
455: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
456: PetscStrlcat(mtypes,((PetscObject)T)->type_name,sizeof(mtypes));
457: PetscStrlcat(mtypes,"_",sizeof(mtypes));
458: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
459: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
461: product->type = MATPRODUCT_ABt;
462: PetscInfo1(mat," querying %s (ABt)\n",f);
463: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
464: }
466: if (!f) {
467: PetscInfo1(mat," querying %s\n",f);
468: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
469: }
470: }
471: if (!f) {
472: SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
473: }
475: (*f)(mat);
476: return(0);
477: }
479: static PetscErrorCode MatProductSetFromOptions_PtAP(Mat mat)
480: {
482: Mat_Product *product = mat->product;
483: Mat A=product->A,B=product->B;
484: PetscBool sametype;
485: PetscErrorCode (*fA)(Mat);
486: PetscErrorCode (*fB)(Mat);
487: PetscErrorCode (*f)(Mat)=NULL;
490: /* Check matrix global sizes */
491: if (A->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N);
492: if (B->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N);
494: fA = A->ops->productsetfromoptions;
495: fB = B->ops->productsetfromoptions;
497: PetscInfo2(mat,"for A %s, P %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
498: PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);
499: if (fB == fA && sametype) {
500: PetscInfo(mat," sametype and matching op\n");
501: f = fB;
502: } else {
503: /* query MatProductSetFromOptions_Atype_Btype */
504: char mtypes[256];
505: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
506: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
507: PetscStrlcat(mtypes,"_",sizeof(mtypes));
508: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
509: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
510: PetscInfo1(mat," querying %s\n",f);
511: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
513: if (!f) {
514: PetscInfo1(mat," querying %s\n",f);
515: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
516: }
517: }
519: if (f) {
520: (*f)(mat);
521: } else {
522: mat->ops->productsymbolic = MatProductSymbolic_Basic;
523: PetscInfo2(mat," for A %s, P %s uses MatProduct_Basic() implementation",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
524: }
525: return(0);
526: }
528: static PetscErrorCode MatProductSetFromOptions_RARt(Mat mat)
529: {
531: Mat_Product *product = mat->product;
532: Mat A=product->A,B=product->B;
533: PetscBool sametype;
534: PetscErrorCode (*fA)(Mat);
535: PetscErrorCode (*fB)(Mat);
536: PetscErrorCode (*f)(Mat)=NULL;
539: /* Check matrix global sizes */
540: if (A->rmap->N != B->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N);
541: if (B->cmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->cmap->N,A->cmap->N);
543: fA = A->ops->productsetfromoptions;
544: fB = B->ops->productsetfromoptions;
546: PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);
547: PetscInfo2(mat,"for A %s, B %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
548: if (fB == fA && sametype) {
549: PetscInfo(mat," sametype and matching op\n");
550: f = fB;
551: } else {
552: /* query MatProductSetFromOptions_Atype_Btype */
553: char mtypes[256];
554: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
555: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
556: PetscStrlcat(mtypes,"_",sizeof(mtypes));
557: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
558: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
559: PetscInfo1(mat," querying %s\n",f);
560: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
562: if (!f) {
563: PetscInfo1(mat," querying %s\n",f);
564: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
565: }
566: }
568: if (f) {
569: (*f)(mat);
570: } else {
571: PetscInfo2(mat," for A %s, P %s uses MatProduct_Basic() implementation",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
572: mat->ops->productsymbolic = MatProductSymbolic_Basic;
573: }
574: return(0);
575: }
577: static PetscErrorCode MatProductSetFromOptions_ABC(Mat mat)
578: {
580: Mat_Product *product = mat->product;
581: Mat A=product->A,B=product->B,C=product->C;
582: PetscErrorCode (*fA)(Mat);
583: PetscErrorCode (*fB)(Mat);
584: PetscErrorCode (*fC)(Mat);
585: PetscErrorCode (*f)(Mat)=NULL;
588: /* Check matrix global sizes */
589: if (B->rmap->N!= A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N);
590: if (C->rmap->N!= B->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",C->rmap->N,B->cmap->N);
592: fA = A->ops->productsetfromoptions;
593: fB = B->ops->productsetfromoptions;
594: fC = C->ops->productsetfromoptions;
595: PetscInfo3(mat,"for A %s, B %s and C %s\n",((PetscObject)A)->type_name,((PetscObject)B)->type_name,((PetscObject)C)->type_name);
596: if (fA == fB && fA == fC && fA) {
597: PetscInfo(mat," matching op\n");
598: f = fA;
599: } else {
600: /* query MatProductSetFromOptions_Atype_Btype_Ctype */
601: char mtypes[256];
602: PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));
603: PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));
604: PetscStrlcat(mtypes,"_",sizeof(mtypes));
605: PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));
606: PetscStrlcat(mtypes,"_",sizeof(mtypes));
607: PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));
608: PetscStrlcat(mtypes,"_C",sizeof(mtypes));
610: PetscInfo1(mat," querying %s\n",f);
611: PetscObjectQueryFunction((PetscObject)A,mtypes,&f);
612: if (!f) {
613: PetscInfo1(mat," querying %s\n",f);
614: PetscObjectQueryFunction((PetscObject)B,mtypes,&f);
615: }
616: if (!f) {
617: PetscInfo1(mat," querying %s\n",f);
618: PetscObjectQueryFunction((PetscObject)C,mtypes,&f);
619: }
620: }
622: if (f) {
623: (*f)(mat);
624: } else { /* use MatProductSymbolic/Numeric_Basic() */
625: PetscInfo3(mat," for A %s, B %s and C %s uses MatProduct_Basic() implementation",((PetscObject)A)->type_name,((PetscObject)B)->type_name,((PetscObject)C)->type_name);
626: mat->ops->productsymbolic = MatProductSymbolic_Basic;
627: }
628: return(0);
629: }
631: /*@C
632: MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.
634: Logically Collective on Mat
636: Input Parameter:
637: . mat - the matrix
639: Level: beginner
641: .seealso: MatSetFromOptions()
642: @*/
643: PetscErrorCode MatProductSetFromOptions(Mat mat)
644: {
650: if (mat->ops->productsetfromoptions) {
651: (*mat->ops->productsetfromoptions)(mat);
652: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Call MatProductSetType() first");
653: return(0);
654: }
656: /* ----------------------------------------------- */
657: PetscErrorCode MatProductNumeric_AB(Mat mat)
658: {
660: Mat_Product *product = mat->product;
661: Mat A=product->A,B=product->B;
664: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
665: (mat->ops->matmultnumeric)(A,B,mat);
666: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
667: return(0);
668: }
670: PetscErrorCode MatProductNumeric_AtB(Mat mat)
671: {
673: Mat_Product *product = mat->product;
674: Mat A=product->A,B=product->B;
677: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
678: (mat->ops->transposematmultnumeric)(A,B,mat);
679: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
680: return(0);
681: }
683: PetscErrorCode MatProductNumeric_ABt(Mat mat)
684: {
686: Mat_Product *product = mat->product;
687: Mat A=product->A,B=product->B;
690: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
691: (mat->ops->mattransposemultnumeric)(A,B,mat);
692: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
693: return(0);
694: }
696: PetscErrorCode MatProductNumeric_PtAP(Mat mat)
697: {
699: Mat_Product *product = mat->product;
700: Mat A=product->A,B=product->B;
703: PetscLogEventBegin(MAT_PtAPNumeric,mat,0,0,0);
704: (mat->ops->ptapnumeric)(A,B,mat);
705: PetscLogEventEnd(MAT_PtAPNumeric,mat,0,0,0);
706: return(0);
707: }
709: PetscErrorCode MatProductNumeric_RARt(Mat mat)
710: {
712: Mat_Product *product = mat->product;
713: Mat A=product->A,B=product->B;
716: PetscLogEventBegin(MAT_RARtNumeric,A,B,0,0);
717: (mat->ops->rartnumeric)(A,B,mat);
718: PetscLogEventEnd(MAT_RARtNumeric,A,B,0,0);
719: return(0);
720: }
722: PetscErrorCode MatProductNumeric_ABC(Mat mat)
723: {
725: Mat_Product *product = mat->product;
726: Mat A=product->A,B=product->B,C=product->C;
729: PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);
730: (mat->ops->matmatmultnumeric)(A,B,C,mat);
731: PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);
732: return(0);
733: }
735: /*@
736: MatProductNumeric - Implement a matrix product with numerical values.
738: Collective on Mat
740: Input Parameters:
741: . mat - the matrix to hold a product
743: Output Parameters:
744: . mat - the matrix product
746: Level: intermediate
748: .seealso: MatProductCreate(), MatSetType()
749: @*/
750: PetscErrorCode MatProductNumeric(Mat mat)
751: {
758: if (mat->ops->productnumeric) {
759: (*mat->ops->productnumeric)(mat);
760: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_USER,"Call MatProductSymbolic() first");
761: return(0);
762: }
764: /* ----------------------------------------------- */
765: PetscErrorCode MatProductSymbolic_AB(Mat mat)
766: {
768: Mat_Product *product = mat->product;
769: Mat A=product->A,B=product->B;
772: (mat->ops->matmultsymbolic)(A,B,product->fill,mat);
773: mat->ops->productnumeric = MatProductNumeric_AB;
774: return(0);
775: }
777: PetscErrorCode MatProductSymbolic_AtB(Mat mat)
778: {
780: Mat_Product *product = mat->product;
781: Mat A=product->A,B=product->B;
784: (mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);
785: mat->ops->productnumeric = MatProductNumeric_AtB;
786: return(0);
787: }
789: PetscErrorCode MatProductSymbolic_ABt(Mat mat)
790: {
792: Mat_Product *product = mat->product;
793: Mat A=product->A,B=product->B;
796: (mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);
797: mat->ops->productnumeric = MatProductNumeric_ABt;
798: return(0);
799: }
801: PetscErrorCode MatProductSymbolic_ABC(Mat mat)
802: {
804: Mat_Product *product = mat->product;
805: Mat A=product->A,B=product->B,C=product->C;
808: (mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);
809: mat->ops->productnumeric = MatProductNumeric_ABC;
810: return(0);
811: }
813: /*@
814: MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.
816: Collective on Mat
818: Input Parameters:
819: . mat - the matrix to hold a product
821: Output Parameters:
822: . mat - the matrix product data structure
824: Level: intermediate
826: .seealso: MatProductCreate(), MatSetType(), MatProductNumeric(), MatProductType, MatProductAlgorithm
827: @*/
828: PetscErrorCode MatProductSymbolic(Mat mat)
829: {
831: Mat_Product *product = mat->product;
832: MatProductType productype = product->type;
833: PetscLogEvent eventtype=-1;
838: /* log event */
839: switch (productype) {
840: case MATPRODUCT_AB:
841: eventtype = MAT_MatMultSymbolic;
842: break;
843: case MATPRODUCT_AtB:
844: eventtype = MAT_TransposeMatMultSymbolic;
845: break;
846: case MATPRODUCT_ABt:
847: eventtype = MAT_MatTransposeMultSymbolic;
848: break;
849: case MATPRODUCT_PtAP:
850: eventtype = MAT_PtAPSymbolic;
851: break;
852: case MATPRODUCT_RARt:
853: eventtype = MAT_RARtSymbolic;
854: break;
855: case MATPRODUCT_ABC:
856: eventtype = MAT_MatMatMultSymbolic;
857: break;
858: default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MATPRODUCT type is not supported");
859: }
861: if (mat->ops->productsymbolic) {
862: PetscLogEventBegin(eventtype,mat,0,0,0);
863: (*mat->ops->productsymbolic)(mat);
864: PetscLogEventEnd(eventtype,mat,0,0,0);
865: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_USER,"Call MatProductSetFromOptions() first");
866: return(0);
867: }
869: /*@
870: MatProductSetFill - Set an expected fill of the matrix product.
872: Collective on Mat
874: Input Parameters:
875: + mat - the matrix product
876: - 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 irrelevent.
878: Level: intermediate
880: .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
881: @*/
882: PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
883: {
884: Mat_Product *product = mat->product;
889: if (!product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first");
890: if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) {
891: product->fill = 2.0;
892: } else product->fill = fill;
893: return(0);
894: }
896: /*@
897: MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.
899: Collective on Mat
901: Input Parameters:
902: + mat - the matrix product
903: - alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.
905: Level: intermediate
907: .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate()
908: @*/
909: PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
910: {
911: Mat_Product *product = mat->product;
916: if (!product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first");
917: product->alg = alg;
918: return(0);
919: }
921: /*@
922: MatProductSetType - Sets a particular matrix product type, for example Mat*Mat.
924: Collective on Mat
926: Input Parameters:
927: + mat - the matrix
928: - productype - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.
930: Level: intermediate
932: .seealso: MatProductCreate(), MatProductType, MatProductAlgorithm
933: @*/
934: PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
935: {
937: Mat_Product *product = mat->product;
938: MPI_Comm comm;
943: PetscObjectGetComm((PetscObject)mat,&comm);
944: if (!product) SETERRQ(comm,PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first");
945: product->type = productype;
947: switch (productype) {
948: case MATPRODUCT_AB:
949: mat->ops->productsetfromoptions = MatProductSetFromOptions_AB;
950: break;
951: case MATPRODUCT_AtB:
952: mat->ops->productsetfromoptions = MatProductSetFromOptions_AtB;
953: break;
954: case MATPRODUCT_ABt:
955: mat->ops->productsetfromoptions = MatProductSetFromOptions_ABt;
956: break;
957: case MATPRODUCT_PtAP:
958: mat->ops->productsetfromoptions = MatProductSetFromOptions_PtAP;
959: break;
960: case MATPRODUCT_RARt:
961: mat->ops->productsetfromoptions = MatProductSetFromOptions_RARt;
962: break;
963: case MATPRODUCT_ABC:
964: mat->ops->productsetfromoptions = MatProductSetFromOptions_ABC;
965: break;
966: default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]);
967: }
968: return(0);
969: }
971: /*@
972: MatProductClear - Clears matrix product internal structure.
974: Collective on Mat
976: Input Parameters:
977: . mat - the product matrix
979: Level: intermediate
980: @*/
981: PetscErrorCode MatProductClear(Mat mat)
982: {
984: Mat_Product *product = mat->product;
987: if (product) {
988: /* release reference */
989: MatDestroy(&product->A);
990: MatDestroy(&product->B);
991: MatDestroy(&product->C);
992: MatDestroy(&product->Dwork);
993: PetscFree(mat->product);
994: }
995: return(0);
996: }
998: /* Create a supporting struct and attach it to the matrix product */
999: PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
1000: {
1002: Mat_Product *product=NULL;
1005: PetscNewLog(D,&product);
1006: product->A = A;
1007: product->B = B;
1008: product->C = C;
1009: product->Dwork = NULL;
1010: product->alg = MATPRODUCTALGORITHM_DEFAULT;
1011: product->fill = 2.0; /* PETSC_DEFAULT */
1012: product->Areplaced = PETSC_FALSE;
1013: product->Breplaced = PETSC_FALSE;
1014: product->api_user = PETSC_FALSE;
1015: D->product = product;
1017: /* take ownership */
1018: PetscObjectReference((PetscObject)A);
1019: PetscObjectReference((PetscObject)B);
1020: PetscObjectReference((PetscObject)C);
1021: return(0);
1022: }
1024: /*@
1025: MatProductCreateWithMat - Setup a given matrix as a matrix product.
1027: Collective on Mat
1029: Input Parameters:
1030: + A - the first matrix
1031: . B - the second matrix
1032: . C - the third matrix (optional)
1033: - D - the matrix which will be used as a product
1035: Output Parameters:
1036: . D - the product matrix
1038: Level: intermediate
1040: .seealso: MatProductCreate()
1041: @*/
1042: PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
1043: {
1049: MatCheckPreallocated(A,1);
1050: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1051: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1055: MatCheckPreallocated(B,2);
1056: if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1057: if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1059: if (C) {
1062: MatCheckPreallocated(C,3);
1063: if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1064: if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1065: }
1069: MatCheckPreallocated(D,4);
1070: if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1071: if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1073: /* Create a supporting struct and attach it to D */
1074: MatProductCreate_Private(A,B,C,D);
1075: return(0);
1076: }
1078: /*@
1079: MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.
1081: Collective on Mat
1083: Input Parameters:
1084: + A - the first matrix
1085: . B - the second matrix
1086: - C - the third matrix (optional)
1088: Output Parameters:
1089: . D - the product matrix
1091: Level: intermediate
1093: .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1094: @*/
1095: PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1096: {
1102: MatCheckPreallocated(A,1);
1103: if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1104: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1108: MatCheckPreallocated(B,2);
1109: if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1110: if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1112: if (C) {
1115: MatCheckPreallocated(C,3);
1116: if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1117: if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1118: }
1122: MatCreate(PetscObjectComm((PetscObject)A),D);
1123: MatProductCreate_Private(A,B,C,*D);
1124: return(0);
1125: }