xref: /petsc/src/mat/impls/normal/normm.c (revision c3e4dd79e17d3f0a51a524340fae4661630fd09e)
1 
2 #include <petsc/private/matimpl.h>          /*I "petscmat.h" I*/
3 
4 typedef struct {
5   Mat         A;
6   Mat         D; /* local submatrix for diagonal part */
7   Vec         w,left,right,leftwork,rightwork;
8   PetscScalar scale;
9 } Mat_Normal;
10 
11 PetscErrorCode MatScale_Normal(Mat inA,PetscScalar scale)
12 {
13   Mat_Normal *a = (Mat_Normal*)inA->data;
14 
15   PetscFunctionBegin;
16   a->scale *= scale;
17   PetscFunctionReturn(0);
18 }
19 
20 PetscErrorCode MatDiagonalScale_Normal(Mat inA,Vec left,Vec right)
21 {
22   Mat_Normal     *a = (Mat_Normal*)inA->data;
23 
24   PetscFunctionBegin;
25   if (left) {
26     if (!a->left) {
27       PetscCall(VecDuplicate(left,&a->left));
28       PetscCall(VecCopy(left,a->left));
29     } else {
30       PetscCall(VecPointwiseMult(a->left,left,a->left));
31     }
32   }
33   if (right) {
34     if (!a->right) {
35       PetscCall(VecDuplicate(right,&a->right));
36       PetscCall(VecCopy(right,a->right));
37     } else {
38       PetscCall(VecPointwiseMult(a->right,right,a->right));
39     }
40   }
41   PetscFunctionReturn(0);
42 }
43 
44 PetscErrorCode MatIncreaseOverlap_Normal(Mat A,PetscInt is_max,IS is[],PetscInt ov)
45 {
46   Mat_Normal     *a = (Mat_Normal*)A->data;
47   Mat            pattern;
48 
49   PetscFunctionBegin;
50   PetscCheck(ov >= 0,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified");
51   PetscCall(MatProductCreate(a->A,a->A,NULL,&pattern));
52   PetscCall(MatProductSetType(pattern,MATPRODUCT_AtB));
53   PetscCall(MatProductSetFromOptions(pattern));
54   PetscCall(MatProductSymbolic(pattern));
55   PetscCall(MatIncreaseOverlap(pattern,is_max,is,ov));
56   PetscCall(MatDestroy(&pattern));
57   PetscFunctionReturn(0);
58 }
59 
60 PetscErrorCode MatCreateSubMatrices_Normal(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[])
61 {
62   Mat_Normal     *a = (Mat_Normal*)mat->data;
63   Mat            B = a->A, *suba;
64   IS             *row;
65   PetscInt       M;
66 
67   PetscFunctionBegin;
68   PetscCheck(!a->left && !a->right && irow == icol,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Not implemented");
69   if (scall != MAT_REUSE_MATRIX) {
70     PetscCall(PetscCalloc1(n,submat));
71   }
72   PetscCall(MatGetSize(B,&M,NULL));
73   PetscCall(PetscMalloc1(n,&row));
74   PetscCall(ISCreateStride(PETSC_COMM_SELF,M,0,1,&row[0]));
75   PetscCall(ISSetIdentity(row[0]));
76   for (M = 1; M < n; ++M) row[M] = row[0];
77   PetscCall(MatCreateSubMatrices(B,n,row,icol,MAT_INITIAL_MATRIX,&suba));
78   for (M = 0; M < n; ++M) {
79     PetscCall(MatCreateNormal(suba[M],*submat+M));
80     ((Mat_Normal*)(*submat)[M]->data)->scale = a->scale;
81   }
82   PetscCall(ISDestroy(&row[0]));
83   PetscCall(PetscFree(row));
84   PetscCall(MatDestroySubMatrices(n,&suba));
85   PetscFunctionReturn(0);
86 }
87 
88 PetscErrorCode MatPermute_Normal(Mat A,IS rowp,IS colp,Mat *B)
89 {
90   Mat_Normal     *a = (Mat_Normal*)A->data;
91   Mat            C,Aa = a->A;
92   IS             row;
93 
94   PetscFunctionBegin;
95   PetscCheck(rowp == colp,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Row permutation and column permutation must be the same");
96   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa),Aa->rmap->n,Aa->rmap->rstart,1,&row));
97   PetscCall(ISSetIdentity(row));
98   PetscCall(MatPermute(Aa,row,colp,&C));
99   PetscCall(ISDestroy(&row));
100   PetscCall(MatCreateNormal(C,B));
101   PetscCall(MatDestroy(&C));
102   PetscFunctionReturn(0);
103 }
104 
105 PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
106 {
107   Mat_Normal     *a = (Mat_Normal*)A->data;
108   Mat            C;
109 
110   PetscFunctionBegin;
111   PetscCheck(!a->left && !a->right,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
112   PetscCall(MatDuplicate(a->A,op,&C));
113   PetscCall(MatCreateNormal(C,B));
114   PetscCall(MatDestroy(&C));
115   if (op == MAT_COPY_VALUES) ((Mat_Normal*)(*B)->data)->scale = a->scale;
116   PetscFunctionReturn(0);
117 }
118 
119 PetscErrorCode MatCopy_Normal(Mat A,Mat B,MatStructure str)
120 {
121   Mat_Normal     *a = (Mat_Normal*)A->data,*b = (Mat_Normal*)B->data;
122 
123   PetscFunctionBegin;
124   PetscCheck(!a->left && !a->right,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
125   PetscCall(MatCopy(a->A,b->A,str));
126   b->scale = a->scale;
127   PetscCall(VecDestroy(&b->left));
128   PetscCall(VecDestroy(&b->right));
129   PetscCall(VecDestroy(&b->leftwork));
130   PetscCall(VecDestroy(&b->rightwork));
131   PetscFunctionReturn(0);
132 }
133 
134 PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
135 {
136   Mat_Normal     *Na = (Mat_Normal*)N->data;
137   Vec            in;
138 
139   PetscFunctionBegin;
140   in = x;
141   if (Na->right) {
142     if (!Na->rightwork) {
143       PetscCall(VecDuplicate(Na->right,&Na->rightwork));
144     }
145     PetscCall(VecPointwiseMult(Na->rightwork,Na->right,in));
146     in   = Na->rightwork;
147   }
148   PetscCall(MatMult(Na->A,in,Na->w));
149   PetscCall(MatMultTranspose(Na->A,Na->w,y));
150   if (Na->left) PetscCall(VecPointwiseMult(y,Na->left,y));
151   PetscCall(VecScale(y,Na->scale));
152   PetscFunctionReturn(0);
153 }
154 
155 PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
156 {
157   Mat_Normal     *Na = (Mat_Normal*)N->data;
158   Vec            in;
159 
160   PetscFunctionBegin;
161   in = v1;
162   if (Na->right) {
163     if (!Na->rightwork) {
164       PetscCall(VecDuplicate(Na->right,&Na->rightwork));
165     }
166     PetscCall(VecPointwiseMult(Na->rightwork,Na->right,in));
167     in   = Na->rightwork;
168   }
169   PetscCall(MatMult(Na->A,in,Na->w));
170   PetscCall(VecScale(Na->w,Na->scale));
171   if (Na->left) {
172     PetscCall(MatMultTranspose(Na->A,Na->w,v3));
173     PetscCall(VecPointwiseMult(v3,Na->left,v3));
174     PetscCall(VecAXPY(v3,1.0,v2));
175   } else {
176     PetscCall(MatMultTransposeAdd(Na->A,Na->w,v2,v3));
177   }
178   PetscFunctionReturn(0);
179 }
180 
181 PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y)
182 {
183   Mat_Normal     *Na = (Mat_Normal*)N->data;
184   Vec            in;
185 
186   PetscFunctionBegin;
187   in = x;
188   if (Na->left) {
189     if (!Na->leftwork) {
190       PetscCall(VecDuplicate(Na->left,&Na->leftwork));
191     }
192     PetscCall(VecPointwiseMult(Na->leftwork,Na->left,in));
193     in   = Na->leftwork;
194   }
195   PetscCall(MatMult(Na->A,in,Na->w));
196   PetscCall(MatMultTranspose(Na->A,Na->w,y));
197   if (Na->right) PetscCall(VecPointwiseMult(y,Na->right,y));
198   PetscCall(VecScale(y,Na->scale));
199   PetscFunctionReturn(0);
200 }
201 
202 PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
203 {
204   Mat_Normal     *Na = (Mat_Normal*)N->data;
205   Vec            in;
206 
207   PetscFunctionBegin;
208   in = v1;
209   if (Na->left) {
210     if (!Na->leftwork) {
211       PetscCall(VecDuplicate(Na->left,&Na->leftwork));
212     }
213     PetscCall(VecPointwiseMult(Na->leftwork,Na->left,in));
214     in   = Na->leftwork;
215   }
216   PetscCall(MatMult(Na->A,in,Na->w));
217   PetscCall(VecScale(Na->w,Na->scale));
218   if (Na->right) {
219     PetscCall(MatMultTranspose(Na->A,Na->w,v3));
220     PetscCall(VecPointwiseMult(v3,Na->right,v3));
221     PetscCall(VecAXPY(v3,1.0,v2));
222   } else {
223     PetscCall(MatMultTransposeAdd(Na->A,Na->w,v2,v3));
224   }
225   PetscFunctionReturn(0);
226 }
227 
228 PetscErrorCode MatDestroy_Normal(Mat N)
229 {
230   Mat_Normal     *Na = (Mat_Normal*)N->data;
231 
232   PetscFunctionBegin;
233   PetscCall(MatDestroy(&Na->A));
234   PetscCall(MatDestroy(&Na->D));
235   PetscCall(VecDestroy(&Na->w));
236   PetscCall(VecDestroy(&Na->left));
237   PetscCall(VecDestroy(&Na->right));
238   PetscCall(VecDestroy(&Na->leftwork));
239   PetscCall(VecDestroy(&Na->rightwork));
240   PetscCall(PetscFree(N->data));
241   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMat_C",NULL));
242   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_seqaij_C",NULL));
243   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_mpiaij_C",NULL));
244   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_seqdense_C",NULL));
245   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_mpidense_C",NULL));
246   PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_dense_C",NULL));
247   PetscFunctionReturn(0);
248 }
249 
250 /*
251       Slow, nonscalable version
252 */
253 PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
254 {
255   Mat_Normal        *Na = (Mat_Normal*)N->data;
256   Mat               A   = Na->A;
257   PetscInt          i,j,rstart,rend,nnz;
258   const PetscInt    *cols;
259   PetscScalar       *diag,*work,*values;
260   const PetscScalar *mvalues;
261 
262   PetscFunctionBegin;
263   PetscCall(PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work));
264   PetscCall(PetscArrayzero(work,A->cmap->N));
265   PetscCall(MatGetOwnershipRange(A,&rstart,&rend));
266   for (i=rstart; i<rend; i++) {
267     PetscCall(MatGetRow(A,i,&nnz,&cols,&mvalues));
268     for (j=0; j<nnz; j++) {
269       work[cols[j]] += mvalues[j]*mvalues[j];
270     }
271     PetscCall(MatRestoreRow(A,i,&nnz,&cols,&mvalues));
272   }
273   PetscCall(MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N)));
274   rstart = N->cmap->rstart;
275   rend   = N->cmap->rend;
276   PetscCall(VecGetArray(v,&values));
277   PetscCall(PetscArraycpy(values,diag+rstart,rend-rstart));
278   PetscCall(VecRestoreArray(v,&values));
279   PetscCall(PetscFree2(diag,work));
280   PetscCall(VecScale(v,Na->scale));
281   PetscFunctionReturn(0);
282 }
283 
284 PetscErrorCode MatGetDiagonalBlock_Normal(Mat N,Mat *D)
285 {
286   Mat_Normal *Na = (Mat_Normal*)N->data;
287   Mat        M,A = Na->A;
288 
289   PetscFunctionBegin;
290   PetscCall(MatGetDiagonalBlock(A,&M));
291   PetscCall(MatCreateNormal(M,&Na->D));
292   *D = Na->D;
293   PetscFunctionReturn(0);
294 }
295 
296 PetscErrorCode MatNormalGetMat_Normal(Mat A,Mat *M)
297 {
298   Mat_Normal *Aa = (Mat_Normal*)A->data;
299 
300   PetscFunctionBegin;
301   *M = Aa->A;
302   PetscFunctionReturn(0);
303 }
304 
305 /*@
306       MatNormalGetMat - Gets the Mat object stored inside a MATNORMAL
307 
308    Logically collective on Mat
309 
310    Input Parameter:
311 .   A  - the MATNORMAL matrix
312 
313    Output Parameter:
314 .   M - the matrix object stored inside A
315 
316    Level: intermediate
317 
318 .seealso: `MatCreateNormal()`
319 
320 @*/
321 PetscErrorCode MatNormalGetMat(Mat A,Mat *M)
322 {
323   PetscFunctionBegin;
324   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
325   PetscValidType(A,1);
326   PetscValidPointer(M,2);
327   PetscUseMethod(A,"MatNormalGetMat_C",(Mat,Mat*),(A,M));
328   PetscFunctionReturn(0);
329 }
330 
331 PetscErrorCode MatConvert_Normal_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
332 {
333   Mat_Normal     *Aa = (Mat_Normal*)A->data;
334   Mat            B;
335   PetscInt       m,n,M,N;
336 
337   PetscFunctionBegin;
338   PetscCall(MatGetSize(A,&M,&N));
339   PetscCall(MatGetLocalSize(A,&m,&n));
340   if (reuse == MAT_REUSE_MATRIX) {
341     B = *newmat;
342     PetscCall(MatProductReplaceMats(Aa->A,Aa->A,NULL,B));
343   } else {
344     PetscCall(MatProductCreate(Aa->A,Aa->A,NULL,&B));
345     PetscCall(MatProductSetType(B,MATPRODUCT_AtB));
346     PetscCall(MatProductSetFromOptions(B));
347     PetscCall(MatProductSymbolic(B));
348     PetscCall(MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE));
349   }
350   PetscCall(MatProductNumeric(B));
351   if (reuse == MAT_INPLACE_MATRIX) {
352     PetscCall(MatHeaderReplace(A,&B));
353   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
354   PetscCall(MatConvert(*newmat,MATAIJ,MAT_INPLACE_MATRIX,newmat));
355   PetscFunctionReturn(0);
356 }
357 
358 typedef struct {
359   Mat work[2];
360 } Normal_Dense;
361 
362 PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
363 {
364   Mat            A,B;
365   Normal_Dense   *contents;
366   Mat_Normal     *a;
367   PetscScalar    *array;
368 
369   PetscFunctionBegin;
370   MatCheckProduct(C,3);
371   A = C->product->A;
372   a = (Mat_Normal*)A->data;
373   B = C->product->B;
374   contents = (Normal_Dense*)C->product->data;
375   PetscCheck(contents,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
376   if (a->right) {
377     PetscCall(MatCopy(B,C,SAME_NONZERO_PATTERN));
378     PetscCall(MatDiagonalScale(C,a->right,NULL));
379   }
380   PetscCall(MatProductNumeric(contents->work[0]));
381   PetscCall(MatDenseGetArrayWrite(C,&array));
382   PetscCall(MatDensePlaceArray(contents->work[1],array));
383   PetscCall(MatProductNumeric(contents->work[1]));
384   PetscCall(MatDenseRestoreArrayWrite(C,&array));
385   PetscCall(MatDenseResetArray(contents->work[1]));
386   PetscCall(MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE));
387   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
388   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
389   PetscCall(MatScale(C,a->scale));
390   PetscFunctionReturn(0);
391 }
392 
393 PetscErrorCode MatNormal_DenseDestroy(void *ctx)
394 {
395   Normal_Dense   *contents = (Normal_Dense*)ctx;
396 
397   PetscFunctionBegin;
398   PetscCall(MatDestroy(contents->work));
399   PetscCall(MatDestroy(contents->work+1));
400   PetscCall(PetscFree(contents));
401   PetscFunctionReturn(0);
402 }
403 
404 PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
405 {
406   Mat            A,B;
407   Normal_Dense   *contents = NULL;
408   Mat_Normal     *a;
409   PetscScalar    *array;
410   PetscInt       n,N,m,M;
411 
412   PetscFunctionBegin;
413   MatCheckProduct(C,4);
414   PetscCheck(!C->product->data,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
415   A = C->product->A;
416   a = (Mat_Normal*)A->data;
417   PetscCheck(!a->left,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Not implemented");
418   B = C->product->B;
419   PetscCall(MatGetLocalSize(C,&m,&n));
420   PetscCall(MatGetSize(C,&M,&N));
421   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
422     PetscCall(MatGetLocalSize(B,NULL,&n));
423     PetscCall(MatGetSize(B,NULL,&N));
424     PetscCall(MatGetLocalSize(A,&m,NULL));
425     PetscCall(MatGetSize(A,&M,NULL));
426     PetscCall(MatSetSizes(C,m,n,M,N));
427   }
428   PetscCall(MatSetType(C,((PetscObject)B)->type_name));
429   PetscCall(MatSetUp(C));
430   PetscCall(PetscNew(&contents));
431   C->product->data = contents;
432   C->product->destroy = MatNormal_DenseDestroy;
433   if (a->right) {
434     PetscCall(MatProductCreate(a->A,C,NULL,contents->work));
435   } else {
436     PetscCall(MatProductCreate(a->A,B,NULL,contents->work));
437   }
438   PetscCall(MatProductSetType(contents->work[0],MATPRODUCT_AB));
439   PetscCall(MatProductSetFromOptions(contents->work[0]));
440   PetscCall(MatProductSymbolic(contents->work[0]));
441   PetscCall(MatProductCreate(a->A,contents->work[0],NULL,contents->work+1));
442   PetscCall(MatProductSetType(contents->work[1],MATPRODUCT_AtB));
443   PetscCall(MatProductSetFromOptions(contents->work[1]));
444   PetscCall(MatProductSymbolic(contents->work[1]));
445   PetscCall(MatDenseGetArrayWrite(C,&array));
446   PetscCall(MatSeqDenseSetPreallocation(contents->work[1],array));
447   PetscCall(MatMPIDenseSetPreallocation(contents->work[1],array));
448   PetscCall(MatDenseRestoreArrayWrite(C,&array));
449   C->ops->productnumeric = MatProductNumeric_Normal_Dense;
450   PetscFunctionReturn(0);
451 }
452 
453 PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
454 {
455   PetscFunctionBegin;
456   C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
457   PetscFunctionReturn(0);
458 }
459 
460 PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
461 {
462   Mat_Product    *product = C->product;
463 
464   PetscFunctionBegin;
465   if (product->type == MATPRODUCT_AB) {
466     PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C));
467   }
468   PetscFunctionReturn(0);
469 }
470 
471 /*@
472       MatCreateNormal - Creates a new matrix object that behaves like A'*A.
473 
474    Collective on Mat
475 
476    Input Parameter:
477 .   A  - the (possibly rectangular) matrix
478 
479    Output Parameter:
480 .   N - the matrix that represents A'*A
481 
482    Level: intermediate
483 
484    Notes:
485     The product A'*A is NOT actually formed! Rather the new matrix
486           object performs the matrix-vector product by first multiplying by
487           A and then A'
488 @*/
489 PetscErrorCode  MatCreateNormal(Mat A,Mat *N)
490 {
491   PetscInt       n,nn;
492   Mat_Normal     *Na;
493   VecType        vtype;
494 
495   PetscFunctionBegin;
496   PetscCall(MatGetSize(A,NULL,&nn));
497   PetscCall(MatGetLocalSize(A,NULL,&n));
498   PetscCall(MatCreate(PetscObjectComm((PetscObject)A),N));
499   PetscCall(MatSetSizes(*N,n,n,nn,nn));
500   PetscCall(PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL));
501   PetscCall(PetscLayoutReference(A->cmap,&(*N)->rmap));
502   PetscCall(PetscLayoutReference(A->cmap,&(*N)->cmap));
503 
504   PetscCall(PetscNewLog(*N,&Na));
505   (*N)->data = (void*) Na;
506   PetscCall(PetscObjectReference((PetscObject)A));
507   Na->A      = A;
508   Na->scale  = 1.0;
509 
510   PetscCall(MatCreateVecs(A,NULL,&Na->w));
511 
512   (*N)->ops->destroy           = MatDestroy_Normal;
513   (*N)->ops->mult              = MatMult_Normal;
514   (*N)->ops->multtranspose     = MatMultTranspose_Normal;
515   (*N)->ops->multtransposeadd  = MatMultTransposeAdd_Normal;
516   (*N)->ops->multadd           = MatMultAdd_Normal;
517   (*N)->ops->getdiagonal       = MatGetDiagonal_Normal;
518   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_Normal;
519   (*N)->ops->scale             = MatScale_Normal;
520   (*N)->ops->diagonalscale     = MatDiagonalScale_Normal;
521   (*N)->ops->increaseoverlap   = MatIncreaseOverlap_Normal;
522   (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
523   (*N)->ops->permute           = MatPermute_Normal;
524   (*N)->ops->duplicate         = MatDuplicate_Normal;
525   (*N)->ops->copy              = MatCopy_Normal;
526   (*N)->assembled              = PETSC_TRUE;
527   (*N)->preallocated           = PETSC_TRUE;
528 
529   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMat_C",MatNormalGetMat_Normal));
530   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_seqaij_C",MatConvert_Normal_AIJ));
531   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_mpiaij_C",MatConvert_Normal_AIJ));
532   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_seqdense_C",MatProductSetFromOptions_Normal_Dense));
533   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_mpidense_C",MatProductSetFromOptions_Normal_Dense));
534   PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_dense_C",MatProductSetFromOptions_Normal_Dense));
535   PetscCall(MatSetOption(*N,MAT_SYMMETRIC,PETSC_TRUE));
536   PetscCall(MatGetVecType(A,&vtype));
537   PetscCall(MatSetVecType(*N,vtype));
538 #if defined(PETSC_HAVE_DEVICE)
539   PetscCall(MatBindToCPU(*N,A->boundtocpu));
540 #endif
541   PetscFunctionReturn(0);
542 }
543