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