1
2 /*
3 Routines for matrix products. Calling procedure:
4
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 MATTRANPOSEMAT)
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
21
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
27
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
31
32 MatProductNumeric(D)
33 # Call the numeric phase
34
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 */
40
41 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
42
43 const char *const MatProductTypes[] = {"UNSPECIFIED","AB","AtB","ABt","PtAP","RARt","ABC"};
44
45 /* these are basic implementations relying on the old function pointers
46 * they are dangerous and should be removed in the future */
MatProductNumeric_PtAP_Basic(Mat C)47 static PetscErrorCode MatProductNumeric_PtAP_Basic(Mat C)
48 {
49 PetscErrorCode ierr;
50 Mat_Product *product = C->product;
51 Mat P = product->B,AP = product->Dwork;
52
53 PetscFunctionBegin;
54 /* AP = A*P */
55 ierr = MatProductNumeric(AP);CHKERRQ(ierr);
56 /* C = P^T*AP */
57 if (!C->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
58 ierr = (*C->ops->transposematmultnumeric)(P,AP,C);CHKERRQ(ierr);
59 PetscFunctionReturn(0);
60 }
61
MatProductSymbolic_PtAP_Basic(Mat C)62 static PetscErrorCode MatProductSymbolic_PtAP_Basic(Mat C)
63 {
64 PetscErrorCode ierr;
65 Mat_Product *product = C->product;
66 Mat A=product->A,P=product->B,AP;
67 PetscReal fill=product->fill;
68
69 PetscFunctionBegin;
70 ierr = PetscInfo2((PetscObject)C,"for A %s, P %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);CHKERRQ(ierr);
71 /* AP = A*P */
72 ierr = MatProductCreate(A,P,NULL,&AP);CHKERRQ(ierr);
73 ierr = MatProductSetType(AP,MATPRODUCT_AB);CHKERRQ(ierr);
74 ierr = MatProductSetAlgorithm(AP,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
75 ierr = MatProductSetFill(AP,fill);CHKERRQ(ierr);
76 ierr = MatProductSetFromOptions(AP);CHKERRQ(ierr);
77 ierr = MatProductSymbolic(AP);CHKERRQ(ierr);
78
79 /* C = P^T*AP */
80 ierr = MatProductSetType(C,MATPRODUCT_AtB);CHKERRQ(ierr);
81 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
82 product->A = P;
83 product->B = AP;
84 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
85 ierr = MatProductSymbolic(C);CHKERRQ(ierr);
86
87 /* resume user's original input matrix setting for A and B */
88 product->A = A;
89 product->B = P;
90 product->Dwork = AP;
91
92 C->ops->productnumeric = MatProductNumeric_PtAP_Basic;
93 PetscFunctionReturn(0);
94 }
95
MatProductNumeric_RARt_Basic(Mat C)96 static PetscErrorCode MatProductNumeric_RARt_Basic(Mat C)
97 {
98 PetscErrorCode ierr;
99 Mat_Product *product = C->product;
100 Mat R=product->B,RA=product->Dwork;
101
102 PetscFunctionBegin;
103 /* RA = R*A */
104 ierr = MatProductNumeric(RA);CHKERRQ(ierr);
105 /* C = RA*R^T */
106 if (!C->ops->mattransposemultnumeric) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric stage");
107 ierr = (*C->ops->mattransposemultnumeric)(RA,R,C);CHKERRQ(ierr);
108 PetscFunctionReturn(0);
109 }
110
MatProductSymbolic_RARt_Basic(Mat C)111 static PetscErrorCode MatProductSymbolic_RARt_Basic(Mat C)
112 {
113 PetscErrorCode ierr;
114 Mat_Product *product = C->product;
115 Mat A=product->A,R=product->B,RA;
116 PetscReal fill=product->fill;
117
118 PetscFunctionBegin;
119 ierr = PetscInfo2((PetscObject)C,"for A %s, R %s is used\n",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name);CHKERRQ(ierr);
120 /* RA = R*A */
121 ierr = MatProductCreate(R,A,NULL,&RA);CHKERRQ(ierr);
122 ierr = MatProductSetType(RA,MATPRODUCT_AB);CHKERRQ(ierr);
123 ierr = MatProductSetAlgorithm(RA,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
124 ierr = MatProductSetFill(RA,fill);CHKERRQ(ierr);
125 ierr = MatProductSetFromOptions(RA);CHKERRQ(ierr);
126 ierr = MatProductSymbolic(RA);CHKERRQ(ierr);
127
128 /* C = RA*R^T */
129 ierr = MatProductSetType(C,MATPRODUCT_ABt);CHKERRQ(ierr);
130 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
131 product->A = RA;
132 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr);
133 ierr = MatProductSymbolic(C);CHKERRQ(ierr);
134
135 /* resume user's original input matrix setting for A */
136 product->A = A;
137 product->Dwork = RA; /* save here so it will be destroyed with product C */
138 C->ops->productnumeric = MatProductNumeric_RARt_Basic;
139 PetscFunctionReturn(0);
140 }
141
MatProductNumeric_ABC_Basic(Mat mat)142 static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat)
143 {
144 PetscErrorCode ierr;
145 Mat_Product *product = mat->product;
146 Mat A=product->A,BC=product->Dwork;
147
148 PetscFunctionBegin;
149 /* Numeric BC = B*C */
150 ierr = MatProductNumeric(BC);CHKERRQ(ierr);
151 /* Numeric mat = A*BC */
152 if (!mat->ops->transposematmultnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric stage");
153 ierr = (*mat->ops->matmultnumeric)(A,BC,mat);CHKERRQ(ierr);
154 PetscFunctionReturn(0);
155 }
156
MatProductSymbolic_ABC_Basic(Mat mat)157 static PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat)
158 {
159 PetscErrorCode ierr;
160 Mat_Product *product = mat->product;
161 Mat B=product->B,C=product->C,BC;
162 PetscReal fill=product->fill;
163
164 PetscFunctionBegin;
165 ierr = PetscInfo3((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);CHKERRQ(ierr);
166 /* Symbolic BC = B*C */
167 ierr = MatProductCreate(B,C,NULL,&BC);CHKERRQ(ierr);
168 ierr = MatProductSetType(BC,MATPRODUCT_AB);CHKERRQ(ierr);
169 ierr = MatProductSetAlgorithm(BC,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
170 ierr = MatProductSetFill(BC,fill);CHKERRQ(ierr);
171 ierr = MatProductSetFromOptions(BC);CHKERRQ(ierr);
172 ierr = MatProductSymbolic(BC);CHKERRQ(ierr);
173
174 /* Symbolic mat = A*BC */
175 ierr = MatProductSetType(mat,MATPRODUCT_AB);CHKERRQ(ierr);
176 ierr = MatProductSetAlgorithm(mat,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
177 product->B = BC;
178 product->Dwork = BC;
179 ierr = MatProductSetFromOptions(mat);CHKERRQ(ierr);
180 ierr = MatProductSymbolic(mat);CHKERRQ(ierr);
181
182 /* resume user's original input matrix setting for B */
183 product->B = B;
184 mat->ops->productnumeric = MatProductNumeric_ABC_Basic;
185 PetscFunctionReturn(0);
186 }
187
MatProductSymbolic_Basic(Mat mat)188 static PetscErrorCode MatProductSymbolic_Basic(Mat mat)
189 {
190 PetscErrorCode ierr;
191 Mat_Product *product = mat->product;
192
193 PetscFunctionBegin;
194 switch (product->type) {
195 case MATPRODUCT_PtAP:
196 ierr = MatProductSymbolic_PtAP_Basic(mat);CHKERRQ(ierr);
197 break;
198 case MATPRODUCT_RARt:
199 ierr = MatProductSymbolic_RARt_Basic(mat);CHKERRQ(ierr);
200 break;
201 case MATPRODUCT_ABC:
202 ierr = MatProductSymbolic_ABC_Basic(mat);CHKERRQ(ierr);
203 break;
204 default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[product->type]);
205 }
206 PetscFunctionReturn(0);
207 }
208
209 /* ----------------------------------------------- */
210 /*@C
211 MatProductReplaceMats - Replace input matrices for a matrix product.
212
213 Collective on Mat
214
215 Input Parameters:
216 + A - the matrix or NULL if not being replaced
217 . B - the matrix or NULL if not being replaced
218 . C - the matrix or NULL if not being replaced
219 - D - the matrix product
220
221 Level: intermediate
222
223 Notes:
224 To reuse the symbolic phase, input matrices must have exactly the same data structure as the replaced one.
225 If the type of any of the input matrices is different than what previously used, the product is cleared and MatProductSetFromOptions()/MatProductSymbolic() are invoked again.
226
227 .seealso: MatProductCreate(), MatProductSetFromOptions(), MatProductSymbolic(). MatProductClear()
228 @*/
MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)229 PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D)
230 {
231 PetscErrorCode ierr;
232 Mat_Product *product;
233 PetscBool flgA = PETSC_TRUE,flgB = PETSC_TRUE,flgC = PETSC_TRUE;
234
235 PetscFunctionBegin;
236 PetscValidHeaderSpecific(D,MAT_CLASSID,4);
237 MatCheckProduct(D,4);
238 product = D->product;
239 if (A) {
240 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
241 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
242 ierr = PetscObjectTypeCompare((PetscObject)product->A,((PetscObject)A)->type_name,&flgA);CHKERRQ(ierr);
243 ierr = MatDestroy(&product->A);CHKERRQ(ierr);
244 product->A = A;
245 }
246 if (B) {
247 PetscValidHeaderSpecific(B,MAT_CLASSID,2);
248 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);
249 ierr = PetscObjectTypeCompare((PetscObject)product->B,((PetscObject)B)->type_name,&flgB);CHKERRQ(ierr);
250 ierr = MatDestroy(&product->B);CHKERRQ(ierr);
251 product->B = B;
252 }
253 if (C) {
254 PetscValidHeaderSpecific(C,MAT_CLASSID,3);
255 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
256 ierr = PetscObjectTypeCompare((PetscObject)product->C,((PetscObject)C)->type_name,&flgC);CHKERRQ(ierr);
257 ierr = MatDestroy(&product->C);CHKERRQ(ierr);
258 product->C = C;
259 }
260 /* Any of the replaced mats is of a different type, reset */
261 if (!flgA || !flgB || !flgC) {
262 if (D->product->destroy) {
263 ierr = (*D->product->destroy)(D->product->data);CHKERRQ(ierr);
264 }
265 D->product->destroy = NULL;
266 D->product->data = NULL;
267 if (D->ops->productnumeric || D->ops->productsymbolic) {
268 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr);
269 ierr = MatProductSymbolic(D);CHKERRQ(ierr);
270 }
271 }
272 PetscFunctionReturn(0);
273 }
274
MatProductNumeric_X_Dense(Mat C)275 static PetscErrorCode MatProductNumeric_X_Dense(Mat C)
276 {
277 PetscErrorCode ierr;
278 Mat_Product *product = C->product;
279 Mat A = product->A, B = product->B;
280 PetscInt k, K = B->cmap->N;
281 PetscBool t = PETSC_TRUE,iscuda = PETSC_FALSE;
282 PetscBool Bcpu = PETSC_TRUE, Ccpu = PETSC_TRUE;
283 char *Btype = NULL,*Ctype = NULL;
284
285 PetscFunctionBegin;
286 switch (product->type) {
287 case MATPRODUCT_AB:
288 t = PETSC_FALSE;
289 case MATPRODUCT_AtB:
290 break;
291 default: SETERRQ3(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);
292 }
293 if (PetscDefined(HAVE_CUDA)) {
294 VecType vtype;
295
296 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr);
297 ierr = PetscStrcmp(vtype,VECCUDA,&iscuda);CHKERRQ(ierr);
298 if (!iscuda) {
299 ierr = PetscStrcmp(vtype,VECSEQCUDA,&iscuda);CHKERRQ(ierr);
300 }
301 if (!iscuda) {
302 ierr = PetscStrcmp(vtype,VECMPICUDA,&iscuda);CHKERRQ(ierr);
303 }
304 if (iscuda) { /* Make sure we have up-to-date data on the GPU */
305 ierr = PetscStrallocpy(((PetscObject)B)->type_name,&Btype);CHKERRQ(ierr);
306 ierr = PetscStrallocpy(((PetscObject)C)->type_name,&Ctype);CHKERRQ(ierr);
307 ierr = MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr);
308 if (!C->assembled) { /* need to flag the matrix as assembled, otherwise MatConvert will complain */
309 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
310 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
311 }
312 ierr = MatConvert(C,MATDENSECUDA,MAT_INPLACE_MATRIX,&C);CHKERRQ(ierr);
313 } else { /* Make sure we have up-to-date data on the CPU */
314 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
315 Bcpu = B->boundtocpu;
316 Ccpu = C->boundtocpu;
317 #endif
318 ierr = MatBindToCPU(B,PETSC_TRUE);CHKERRQ(ierr);
319 ierr = MatBindToCPU(C,PETSC_TRUE);CHKERRQ(ierr);
320 }
321 }
322 for (k=0;k<K;k++) {
323 Vec x,y;
324
325 ierr = MatDenseGetColumnVecRead(B,k,&x);CHKERRQ(ierr);
326 ierr = MatDenseGetColumnVecWrite(C,k,&y);CHKERRQ(ierr);
327 if (t) {
328 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
329 } else {
330 ierr = MatMult(A,x,y);CHKERRQ(ierr);
331 }
332 ierr = MatDenseRestoreColumnVecRead(B,k,&x);CHKERRQ(ierr);
333 ierr = MatDenseRestoreColumnVecWrite(C,k,&y);CHKERRQ(ierr);
334 }
335 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
336 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
337 if (PetscDefined(HAVE_CUDA)) {
338 if (iscuda) {
339 ierr = MatConvert(B,Btype,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr);
340 ierr = MatConvert(C,Ctype,MAT_INPLACE_MATRIX,&C);CHKERRQ(ierr);
341 } else {
342 ierr = MatBindToCPU(B,Bcpu);CHKERRQ(ierr);
343 ierr = MatBindToCPU(C,Ccpu);CHKERRQ(ierr);
344 }
345 }
346 ierr = PetscFree(Btype);CHKERRQ(ierr);
347 ierr = PetscFree(Ctype);CHKERRQ(ierr);
348 PetscFunctionReturn(0);
349 }
350
MatProductSymbolic_X_Dense(Mat C)351 static PetscErrorCode MatProductSymbolic_X_Dense(Mat C)
352 {
353 PetscErrorCode ierr;
354 Mat_Product *product = C->product;
355 Mat A = product->A, B = product->B;
356 PetscBool isdense;
357
358 PetscFunctionBegin;
359 switch (product->type) {
360 case MATPRODUCT_AB:
361 ierr = MatSetSizes(C,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);CHKERRQ(ierr);
362 break;
363 case MATPRODUCT_AtB:
364 ierr = MatSetSizes(C,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);CHKERRQ(ierr);
365 break;
366 default: SETERRQ3(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
367 }
368 ierr = PetscObjectBaseTypeCompareAny((PetscObject)C,&isdense,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr);
369 if (!isdense) {
370 ierr = MatSetType(C,((PetscObject)B)->type_name);CHKERRQ(ierr);
371 /* If matrix type of C was not set or not dense, we need to reset the pointer */
372 C->ops->productsymbolic = MatProductSymbolic_X_Dense;
373 }
374 C->ops->productnumeric = MatProductNumeric_X_Dense;
375 ierr = MatSetUp(C);CHKERRQ(ierr);
376 PetscFunctionReturn(0);
377 }
378
379 /* a single driver to query the dispatching */
MatProductSetFromOptions_Private(Mat mat)380 static PetscErrorCode MatProductSetFromOptions_Private(Mat mat)
381 {
382 PetscErrorCode ierr;
383 Mat_Product *product = mat->product;
384 PetscInt Am,An,Bm,Bn,Cm,Cn;
385 Mat A = product->A,B = product->B,C = product->C;
386 const char* const Bnames[] = { "B", "R", "P" };
387 const char* bname;
388 PetscErrorCode (*fA)(Mat);
389 PetscErrorCode (*fB)(Mat);
390 PetscErrorCode (*fC)(Mat);
391 PetscErrorCode (*f)(Mat)=NULL;
392
393 PetscFunctionBegin;
394 mat->ops->productsymbolic = NULL;
395 mat->ops->productnumeric = NULL;
396 if (product->type == MATPRODUCT_UNSPECIFIED) PetscFunctionReturn(0);
397 if (product->type != MATPRODUCT_ABC) C = NULL; /* do not use C if not needed */
398 if (product->type == MATPRODUCT_RARt) bname = Bnames[1];
399 else if (product->type == MATPRODUCT_PtAP) bname = Bnames[2];
400 else bname = Bnames[0];
401
402 /* Check matrices sizes */
403 Am = A->rmap->N;
404 An = A->cmap->N;
405 Bm = B->rmap->N;
406 Bn = B->cmap->N;
407 Cm = C ? C->rmap->N : 0;
408 Cn = C ? C->cmap->N : 0;
409 if (product->type == MATPRODUCT_RARt || product->type == MATPRODUCT_ABt) { PetscInt t = Bn; Bn = Bm; Bm = t; }
410 if (product->type == MATPRODUCT_AtB) { PetscInt t = An; An = Am; Am = t; }
411 if (An != Bm) SETERRQ7(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of A and %s are incompatible for MatProductType %s: A %Dx%D, %s %Dx%D",bname,MatProductTypes[product->type],A->rmap->N,A->cmap->N,bname,B->rmap->N,B->cmap->N);
412 if (Cm && Cm != Bn) SETERRQ5(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions of B and C are incompatible for MatProductType %s: B %Dx%D, C %Dx%D",MatProductTypes[product->type],B->rmap->N,B->cmap->N,Cm,Cn);
413
414 fA = A->ops->productsetfromoptions;
415 fB = B->ops->productsetfromoptions;
416 fC = C ? C->ops->productsetfromoptions : fA;
417 if (C) {
418 ierr = PetscInfo5(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);CHKERRQ(ierr);
419 } else {
420 ierr = PetscInfo4(mat,"MatProductType %s for A %s, %s %s\n",MatProductTypes[product->type],((PetscObject)A)->type_name,bname,((PetscObject)B)->type_name);CHKERRQ(ierr);
421 }
422 if (fA == fB && fA == fC && fA) {
423 ierr = PetscInfo(mat," matching op\n");CHKERRQ(ierr);
424 ierr = (*fA)(mat);CHKERRQ(ierr);
425 } else {
426 /* query MatProductSetFromOptions_Atype_Btype_Ctype */
427 char mtypes[256];
428 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr);
429 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr);
430 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr);
431 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr);
432 if (C) {
433 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr);
434 ierr = PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));CHKERRQ(ierr);
435 }
436 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr);
437
438 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr);
439 ierr = PetscInfo2(mat," querying %s from A? %p\n",mtypes,f);CHKERRQ(ierr);
440 if (!f) {
441 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr);
442 ierr = PetscInfo3(mat," querying %s from %s? %p\n",mtypes,bname,f);CHKERRQ(ierr);
443 }
444 if (!f && C) {
445 ierr = PetscObjectQueryFunction((PetscObject)C,mtypes,&f);CHKERRQ(ierr);
446 ierr = PetscInfo2(mat," querying %s from C? %p\n",mtypes,f);CHKERRQ(ierr);
447 }
448 if (f) { ierr = (*f)(mat);CHKERRQ(ierr); }
449
450 /* We may have found f but it did not succeed */
451 /* some matrices (i.e. MATTRANSPOSE, MATSHELL constructed from MatConvert), knows what to do with their inner matrices */
452 if (!mat->ops->productsymbolic) {
453 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_anytype_C",sizeof(mtypes));CHKERRQ(ierr);
454 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr);
455 ierr = PetscInfo2(mat," querying %s from A? %p\n",mtypes,f);CHKERRQ(ierr);
456 if (!f) {
457 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr);
458 ierr = PetscInfo3(mat," querying %s from %s? %p\n",mtypes,bname,f);CHKERRQ(ierr);
459 }
460 if (!f && C) {
461 ierr = PetscObjectQueryFunction((PetscObject)C,mtypes,&f);CHKERRQ(ierr);
462 ierr = PetscInfo2(mat," querying %s from C? %p\n",mtypes,f);CHKERRQ(ierr);
463 }
464 }
465 if (f) { ierr = (*f)(mat);CHKERRQ(ierr); }
466 }
467
468 /* We may have found f but it did not succeed */
469 if (!mat->ops->productsymbolic) {
470 /* we can still compute the product if B is of type dense */
471 if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) {
472 PetscBool isdense;
473
474 ierr = PetscObjectBaseTypeCompareAny((PetscObject)B,&isdense,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr);
475 if (isdense) {
476
477 mat->ops->productsymbolic = MatProductSymbolic_X_Dense;
478 ierr = PetscInfo(mat," using basic looping over columns of a dense matrix\n");CHKERRQ(ierr);
479 }
480 } else if (product->type != MATPRODUCT_ABt) { /* use MatProductSymbolic/Numeric_Basic() for triple products only */
481 /*
482 TODO: this should be changed to a proper setfromoptions, not setting the symbolic pointer here, because we do not know if
483 the compination will succeed. In order to be sure, we need MatProductGetProductType to return the type of the result
484 before computing the symbolic phase
485 */
486 ierr = PetscInfo(mat," symbolic product not supported, using MatProductSymbolic_Basic() implementation\n");CHKERRQ(ierr);
487 mat->ops->productsymbolic = MatProductSymbolic_Basic;
488 }
489 }
490 if (!mat->ops->productsymbolic) {
491 ierr = PetscInfo(mat," symbolic product is not supported\n");CHKERRQ(ierr);
492 }
493 PetscFunctionReturn(0);
494 }
495
496 /*@C
497 MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database.
498
499 Logically Collective on Mat
500
501 Input Parameter:
502 . mat - the matrix
503
504 Options Database Keys:
505 . -mat_product_clear - Clear intermediate data structures after MatProductNumeric() has been called
506
507 Level: intermediate
508
509 .seealso: MatSetFromOptions(), MatProductCreate(), MatProductCreateWithMat()
510 @*/
MatProductSetFromOptions(Mat mat)511 PetscErrorCode MatProductSetFromOptions(Mat mat)
512 {
513 PetscErrorCode ierr;
514
515 PetscFunctionBegin;
516 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
517 MatCheckProduct(mat,1);
518 if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot call MatProductSetFromOptions with already present data");
519 ierr = PetscObjectOptionsBegin((PetscObject)mat);CHKERRQ(ierr);
520 ierr = PetscOptionsBool("-mat_product_clear","Clear intermediate data structures after MatProductNumeric() has been called","MatProductClear",mat->product->clear,&mat->product->clear,NULL);CHKERRQ(ierr);
521 ierr = PetscOptionsDeprecated("-mat_freeintermediatedatastructures","-mat_product_clear","3.13","Or call MatProductClear() after MatProductNumeric()");CHKERRQ(ierr);
522 ierr = PetscOptionsEnd();CHKERRQ(ierr);
523 ierr = MatProductSetFromOptions_Private(mat);CHKERRQ(ierr);
524 if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after setup phase");
525 PetscFunctionReturn(0);
526 }
527
528 /*@C
529 MatProductView - View a MatProduct
530
531 Logically Collective on Mat
532
533 Input Parameter:
534 . mat - the matrix obtained with MatProductCreate() or MatProductCreateWithMat()
535
536 Level: intermediate
537
538 .seealso: MatView(), MatProductCreate(), MatProductCreateWithMat()
539 @*/
MatProductView(Mat mat,PetscViewer viewer)540 PetscErrorCode MatProductView(Mat mat, PetscViewer viewer)
541 {
542 PetscErrorCode ierr;
543
544 PetscFunctionBegin;
545 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
546 if (!mat->product) PetscFunctionReturn(0);
547 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mat),&viewer);CHKERRQ(ierr);}
548 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
549 PetscCheckSameComm(mat,1,viewer,2);
550 if (mat->product->view) {
551 ierr = (*mat->product->view)(mat,viewer);CHKERRQ(ierr);
552 }
553 PetscFunctionReturn(0);
554 }
555
556 /* ----------------------------------------------- */
557 /* these are basic implementations relying on the old function pointers
558 * they are dangerous and should be removed in the future */
MatProductNumeric_AB(Mat mat)559 PetscErrorCode MatProductNumeric_AB(Mat mat)
560 {
561 PetscErrorCode ierr;
562 Mat_Product *product = mat->product;
563 Mat A=product->A,B=product->B;
564
565 PetscFunctionBegin;
566 if (!mat->ops->matmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
567 ierr = (*mat->ops->matmultnumeric)(A,B,mat);CHKERRQ(ierr);
568 PetscFunctionReturn(0);
569 }
570
MatProductNumeric_AtB(Mat mat)571 PetscErrorCode MatProductNumeric_AtB(Mat mat)
572 {
573 PetscErrorCode ierr;
574 Mat_Product *product = mat->product;
575 Mat A=product->A,B=product->B;
576
577 PetscFunctionBegin;
578 if (!mat->ops->transposematmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
579 ierr = (*mat->ops->transposematmultnumeric)(A,B,mat);CHKERRQ(ierr);
580 PetscFunctionReturn(0);
581 }
582
MatProductNumeric_ABt(Mat mat)583 PetscErrorCode MatProductNumeric_ABt(Mat mat)
584 {
585 PetscErrorCode ierr;
586 Mat_Product *product = mat->product;
587 Mat A=product->A,B=product->B;
588
589 PetscFunctionBegin;
590 if (!mat->ops->mattransposemultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
591 ierr = (*mat->ops->mattransposemultnumeric)(A,B,mat);CHKERRQ(ierr);
592 PetscFunctionReturn(0);
593 }
594
MatProductNumeric_PtAP(Mat mat)595 PetscErrorCode MatProductNumeric_PtAP(Mat mat)
596 {
597 PetscErrorCode ierr;
598 Mat_Product *product = mat->product;
599 Mat A=product->A,B=product->B;
600
601 PetscFunctionBegin;
602 if (!mat->ops->ptapnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
603 ierr = (*mat->ops->ptapnumeric)(A,B,mat);CHKERRQ(ierr);
604 PetscFunctionReturn(0);
605 }
606
MatProductNumeric_RARt(Mat mat)607 PetscErrorCode MatProductNumeric_RARt(Mat mat)
608 {
609 PetscErrorCode ierr;
610 Mat_Product *product = mat->product;
611 Mat A=product->A,B=product->B;
612
613 PetscFunctionBegin;
614 if (!mat->ops->rartnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
615 ierr = (*mat->ops->rartnumeric)(A,B,mat);CHKERRQ(ierr);
616 PetscFunctionReturn(0);
617 }
618
MatProductNumeric_ABC(Mat mat)619 PetscErrorCode MatProductNumeric_ABC(Mat mat)
620 {
621 PetscErrorCode ierr;
622 Mat_Product *product = mat->product;
623 Mat A=product->A,B=product->B,C=product->C;
624
625 PetscFunctionBegin;
626 if (!mat->ops->matmatmultnumeric) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing numeric implementation of product %s for mat %s",MatProductTypes[product->type],((PetscObject)mat)->type_name);
627 ierr = (*mat->ops->matmatmultnumeric)(A,B,C,mat);CHKERRQ(ierr);
628 PetscFunctionReturn(0);
629 }
630
631 /* ----------------------------------------------- */
632
633 /*@
634 MatProductNumeric - Implement a matrix product with numerical values.
635
636 Collective on Mat
637
638 Input/Output Parameter:
639 . mat - the matrix holding the product
640
641 Level: intermediate
642
643 Notes: MatProductSymbolic() must have been called on mat before calling this function
644
645 .seealso: MatProductCreate(), MatSetType(), MatProductSymbolic()
646 @*/
MatProductNumeric(Mat mat)647 PetscErrorCode MatProductNumeric(Mat mat)
648 {
649 PetscErrorCode ierr;
650 PetscLogEvent eventtype=-1;
651
652 PetscFunctionBegin;
653 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
654 MatCheckProduct(mat,1);
655 /* log event */
656 switch (mat->product->type) {
657 case MATPRODUCT_AB:
658 eventtype = MAT_MatMultNumeric;
659 break;
660 case MATPRODUCT_AtB:
661 eventtype = MAT_TransposeMatMultNumeric;
662 break;
663 case MATPRODUCT_ABt:
664 eventtype = MAT_MatTransposeMultNumeric;
665 break;
666 case MATPRODUCT_PtAP:
667 eventtype = MAT_PtAPNumeric;
668 break;
669 case MATPRODUCT_RARt:
670 eventtype = MAT_RARtNumeric;
671 break;
672 case MATPRODUCT_ABC:
673 eventtype = MAT_MatMatMultNumeric;
674 break;
675 default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
676 }
677 if (mat->ops->productnumeric) {
678 ierr = PetscLogEventBegin(eventtype,mat,0,0,0);CHKERRQ(ierr);
679 ierr = (*mat->ops->productnumeric)(mat);CHKERRQ(ierr);
680 ierr = PetscLogEventEnd(eventtype,mat,0,0,0);CHKERRQ(ierr);
681 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSymbolic() first");
682 if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after numeric phase");
683 if (mat->product->clear) {
684 ierr = MatProductClear(mat);CHKERRQ(ierr);
685 }
686 PetscFunctionReturn(0);
687 }
688
689 /* ----------------------------------------------- */
690 /* these are basic implementations relying on the old function pointers
691 * they are dangerous and should be removed in the future */
MatProductSymbolic_AB(Mat mat)692 PetscErrorCode MatProductSymbolic_AB(Mat mat)
693 {
694 PetscErrorCode ierr;
695 Mat_Product *product = mat->product;
696 Mat A=product->A,B=product->B;
697
698 PetscFunctionBegin;
699 if (!mat->ops->matmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
700 ierr = (*mat->ops->matmultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr);
701 mat->ops->productnumeric = MatProductNumeric_AB;
702 PetscFunctionReturn(0);
703 }
704
MatProductSymbolic_AtB(Mat mat)705 PetscErrorCode MatProductSymbolic_AtB(Mat mat)
706 {
707 PetscErrorCode ierr;
708 Mat_Product *product = mat->product;
709 Mat A=product->A,B=product->B;
710
711 PetscFunctionBegin;
712 if (!mat->ops->transposematmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
713 ierr = (*mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr);
714 mat->ops->productnumeric = MatProductNumeric_AtB;
715 PetscFunctionReturn(0);
716 }
717
MatProductSymbolic_ABt(Mat mat)718 PetscErrorCode MatProductSymbolic_ABt(Mat mat)
719 {
720 PetscErrorCode ierr;
721 Mat_Product *product = mat->product;
722 Mat A=product->A,B=product->B;
723
724 PetscFunctionBegin;
725 if (!mat->ops->mattransposemultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
726 ierr = (*mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr);
727 mat->ops->productnumeric = MatProductNumeric_ABt;
728 PetscFunctionReturn(0);
729 }
730
MatProductSymbolic_ABC(Mat mat)731 PetscErrorCode MatProductSymbolic_ABC(Mat mat)
732 {
733 PetscErrorCode ierr;
734 Mat_Product *product = mat->product;
735 Mat A=product->A,B=product->B,C=product->C;
736
737 PetscFunctionBegin;
738 if (!mat->ops->matmatmultsymbolic) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing symbolic implementation of product %s",MatProductTypes[product->type]);
739 ierr = (*mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);CHKERRQ(ierr);
740 mat->ops->productnumeric = MatProductNumeric_ABC;
741 PetscFunctionReturn(0);
742 }
743
744 /* ----------------------------------------------- */
745
746 /*@
747 MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce.
748
749 Collective on Mat
750
751 Input/Output Parameter:
752 . mat - the matrix to hold a product
753
754 Level: intermediate
755
756 Notes: MatProductSetFromOptions() must have been called on mat before calling this function
757
758 .seealso: MatProductCreate(), MatProductCreateWithMat(), MatProductSetFromOptions(), MatProductNumeric(), MatProductSetType(), MatProductSetAlgorithm()
759 @*/
MatProductSymbolic(Mat mat)760 PetscErrorCode MatProductSymbolic(Mat mat)
761 {
762 PetscErrorCode ierr;
763 PetscLogEvent eventtype=-1;
764
765 PetscFunctionBegin;
766 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
767 MatCheckProduct(mat,1);
768 if (mat->product->data) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Cannot run symbolic phase. Product data not empty");
769
770 /* log event */
771 switch (mat->product->type) {
772 case MATPRODUCT_AB:
773 eventtype = MAT_MatMultSymbolic;
774 break;
775 case MATPRODUCT_AtB:
776 eventtype = MAT_TransposeMatMultSymbolic;
777 break;
778 case MATPRODUCT_ABt:
779 eventtype = MAT_MatTransposeMultSymbolic;
780 break;
781 case MATPRODUCT_PtAP:
782 eventtype = MAT_PtAPSymbolic;
783 break;
784 case MATPRODUCT_RARt:
785 eventtype = MAT_RARtSymbolic;
786 break;
787 case MATPRODUCT_ABC:
788 eventtype = MAT_MatMatMultSymbolic;
789 break;
790 default: SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[mat->product->type]);
791 }
792
793 mat->ops->productnumeric = NULL;
794 if (mat->ops->productsymbolic) {
795 ierr = PetscLogEventBegin(eventtype,mat,0,0,0);CHKERRQ(ierr);
796 ierr = (*mat->ops->productsymbolic)(mat);CHKERRQ(ierr);
797 ierr = PetscLogEventEnd(eventtype,mat,0,0,0);CHKERRQ(ierr);
798 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ORDER,"Call MatProductSetFromOptions() first");
799 if (!mat->product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Missing product after symbolic phase");
800 if (!mat->ops->productnumeric) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_PLIB,"Symbolic phase did not specify the numeric phase");
801 PetscFunctionReturn(0);
802 }
803
804 /*@
805 MatProductSetFill - Set an expected fill of the matrix product.
806
807 Collective on Mat
808
809 Input Parameters:
810 + mat - the matrix product
811 - 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.
812
813 Level: intermediate
814
815 .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate()
816 @*/
MatProductSetFill(Mat mat,PetscReal fill)817 PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill)
818 {
819 PetscFunctionBegin;
820 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
821 MatCheckProduct(mat,1);
822 if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) mat->product->fill = 2.0;
823 else mat->product->fill = fill;
824 PetscFunctionReturn(0);
825 }
826
827 /*@
828 MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation.
829
830 Collective on Mat
831
832 Input Parameters:
833 + mat - the matrix product
834 - alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT.
835
836 Level: intermediate
837
838 .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate(), MatProductAlgorithm
839 @*/
MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)840 PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg)
841 {
842 PetscErrorCode ierr;
843
844 PetscFunctionBegin;
845 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
846 MatCheckProduct(mat,1);
847 ierr = PetscFree(mat->product->alg);CHKERRQ(ierr);
848 ierr = PetscStrallocpy(alg,&mat->product->alg);CHKERRQ(ierr);
849 PetscFunctionReturn(0);
850 }
851
852 /*@
853 MatProductSetType - Sets a particular matrix product type
854
855 Collective on Mat
856
857 Input Parameters:
858 + mat - the matrix
859 - productype - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC.
860
861 Level: intermediate
862
863 .seealso: MatProductCreate(), MatProductType
864 @*/
MatProductSetType(Mat mat,MatProductType productype)865 PetscErrorCode MatProductSetType(Mat mat,MatProductType productype)
866 {
867 PetscErrorCode ierr;
868
869 PetscFunctionBegin;
870 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
871 MatCheckProduct(mat,1);
872 PetscValidLogicalCollectiveEnum(mat,productype,2);
873 if (productype != mat->product->type) {
874 if (mat->product->destroy) {
875 ierr = (*mat->product->destroy)(mat->product->data);CHKERRQ(ierr);
876 }
877 mat->product->destroy = NULL;
878 mat->product->data = NULL;
879 mat->ops->productsymbolic = NULL;
880 mat->ops->productnumeric = NULL;
881 }
882 mat->product->type = productype;
883 PetscFunctionReturn(0);
884 }
885
886 /*@
887 MatProductClear - Clears matrix product internal structure.
888
889 Collective on Mat
890
891 Input Parameters:
892 . mat - the product matrix
893
894 Level: intermediate
895
896 Notes: this function should be called to remove any intermediate data used by the product
897 After having called this function, MatProduct operations can no longer be used on mat
898 @*/
MatProductClear(Mat mat)899 PetscErrorCode MatProductClear(Mat mat)
900 {
901 PetscErrorCode ierr;
902 Mat_Product *product = mat->product;
903
904 PetscFunctionBegin;
905 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
906 if (product) {
907 ierr = MatDestroy(&product->A);CHKERRQ(ierr);
908 ierr = MatDestroy(&product->B);CHKERRQ(ierr);
909 ierr = MatDestroy(&product->C);CHKERRQ(ierr);
910 ierr = PetscFree(product->alg);CHKERRQ(ierr);
911 ierr = MatDestroy(&product->Dwork);CHKERRQ(ierr);
912 if (product->destroy) {
913 ierr = (*product->destroy)(product->data);CHKERRQ(ierr);
914 }
915 }
916 ierr = PetscFree(mat->product);CHKERRQ(ierr);
917 mat->ops->productsymbolic = NULL;
918 mat->ops->productnumeric = NULL;
919 PetscFunctionReturn(0);
920 }
921
922 /* Create a supporting struct and attach it to the matrix product */
MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)923 PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D)
924 {
925 PetscErrorCode ierr;
926 Mat_Product *product=NULL;
927
928 PetscFunctionBegin;
929 PetscValidHeaderSpecific(D,MAT_CLASSID,4);
930 if (D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product already present");
931 ierr = PetscNewLog(D,&product);CHKERRQ(ierr);
932 product->A = A;
933 product->B = B;
934 product->C = C;
935 product->type = MATPRODUCT_UNSPECIFIED;
936 product->Dwork = NULL;
937 product->api_user = PETSC_FALSE;
938 product->clear = PETSC_FALSE;
939 D->product = product;
940
941 ierr = MatProductSetAlgorithm(D,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr);
942 ierr = MatProductSetFill(D,PETSC_DEFAULT);CHKERRQ(ierr);
943
944 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
945 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);
946 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
947 PetscFunctionReturn(0);
948 }
949
950 /*@
951 MatProductCreateWithMat - Setup a given matrix as a matrix product.
952
953 Collective on Mat
954
955 Input Parameters:
956 + A - the first matrix
957 . B - the second matrix
958 . C - the third matrix (optional)
959 - D - the matrix which will be used as a product
960
961 Output Parameters:
962 . D - the product matrix
963
964 Notes:
965 Any product data attached to D will be cleared
966
967 Level: intermediate
968
969 .seealso: MatProductCreate(), MatProductClear()
970 @*/
MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)971 PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D)
972 {
973 PetscErrorCode ierr;
974
975 PetscFunctionBegin;
976 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
977 PetscValidType(A,1);
978 MatCheckPreallocated(A,1);
979 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
980 if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
981
982 PetscValidHeaderSpecific(B,MAT_CLASSID,2);
983 PetscValidType(B,2);
984 MatCheckPreallocated(B,2);
985 if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
986 if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
987
988 if (C) {
989 PetscValidHeaderSpecific(C,MAT_CLASSID,3);
990 PetscValidType(C,3);
991 MatCheckPreallocated(C,3);
992 if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
993 if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
994 }
995
996 PetscValidHeaderSpecific(D,MAT_CLASSID,4);
997 PetscValidType(D,4);
998 MatCheckPreallocated(D,4);
999 if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1000 if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1001
1002 /* Create a supporting struct and attach it to D */
1003 ierr = MatProductClear(D);CHKERRQ(ierr);
1004 ierr = MatProductCreate_Private(A,B,C,D);CHKERRQ(ierr);
1005 PetscFunctionReturn(0);
1006 }
1007
1008 /*@
1009 MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations.
1010
1011 Collective on Mat
1012
1013 Input Parameters:
1014 + A - the first matrix
1015 . B - the second matrix
1016 - C - the third matrix (optional)
1017
1018 Output Parameters:
1019 . D - the product matrix
1020
1021 Level: intermediate
1022
1023 .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm()
1024 @*/
MatProductCreate(Mat A,Mat B,Mat C,Mat * D)1025 PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D)
1026 {
1027 PetscErrorCode ierr;
1028
1029 PetscFunctionBegin;
1030 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1031 PetscValidType(A,1);
1032 MatCheckPreallocated(A,1);
1033 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1034 if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1035
1036 PetscValidHeaderSpecific(B,MAT_CLASSID,2);
1037 PetscValidType(B,2);
1038 MatCheckPreallocated(B,2);
1039 if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1040 if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1041
1042 if (C) {
1043 PetscValidHeaderSpecific(C,MAT_CLASSID,3);
1044 PetscValidType(C,3);
1045 MatCheckPreallocated(C,3);
1046 if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1047 if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1048 }
1049
1050 PetscValidPointer(D,4);
1051 ierr = MatCreate(PetscObjectComm((PetscObject)A),D);CHKERRQ(ierr);
1052 ierr = MatProductCreate_Private(A,B,C,*D);CHKERRQ(ierr);
1053 PetscFunctionReturn(0);
1054 }
1055