Actual source code: normm.c


  2: #include <petsc/private/matimpl.h>

  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;

 11: PetscErrorCode MatScale_Normal(Mat inA, PetscScalar scale)
 12: {
 13:   Mat_Normal *a = (Mat_Normal *)inA->data;

 15:   PetscFunctionBegin;
 16:   a->scale *= scale;
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

 20: PetscErrorCode MatDiagonalScale_Normal(Mat inA, Vec left, Vec right)
 21: {
 22:   Mat_Normal *a = (Mat_Normal *)inA->data;

 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(PETSC_SUCCESS);
 42: }

 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;

 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(PETSC_SUCCESS);
 58: }

 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;

 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(PETSC_SUCCESS);
 84: }

 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;

 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(PETSC_SUCCESS);
101: }

103: PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
104: {
105:   Mat_Normal *a = (Mat_Normal *)A->data;
106:   Mat         C;

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(PETSC_SUCCESS);
115: }

117: PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str)
118: {
119:   Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data;

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(PETSC_SUCCESS);
130: }

132: PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y)
133: {
134:   Mat_Normal *Na = (Mat_Normal *)N->data;
135:   Vec         in;

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(PETSC_SUCCESS);
149: }

151: PetscErrorCode MatMultAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
152: {
153:   Mat_Normal *Na = (Mat_Normal *)N->data;
154:   Vec         in;

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(PETSC_SUCCESS);
173: }

175: PetscErrorCode MatMultTranspose_Normal(Mat N, Vec x, Vec y)
176: {
177:   Mat_Normal *Na = (Mat_Normal *)N->data;
178:   Vec         in;

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(PETSC_SUCCESS);
192: }

194: PetscErrorCode MatMultTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
195: {
196:   Mat_Normal *Na = (Mat_Normal *)N->data;
197:   Vec         in;

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(PETSC_SUCCESS);
216: }

218: PetscErrorCode MatDestroy_Normal(Mat N)
219: {
220:   Mat_Normal *Na = (Mat_Normal *)N->data;

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: #if defined(PETSC_HAVE_HYPRE)
235:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_hypre_C", NULL));
236: #endif
237:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL));
238:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL));
239:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_dense_C", NULL));
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*
244:       Slow, nonscalable version
245: */
246: PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
247: {
248:   Mat_Normal        *Na = (Mat_Normal *)N->data;
249:   Mat                A  = Na->A;
250:   PetscInt           i, j, rstart, rend, nnz;
251:   const PetscInt    *cols;
252:   PetscScalar       *diag, *work, *values;
253:   const PetscScalar *mvalues;

255:   PetscFunctionBegin;
256:   PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
257:   PetscCall(PetscArrayzero(work, A->cmap->N));
258:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
259:   for (i = rstart; i < rend; i++) {
260:     PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
261:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
262:     PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
263:   }
264:   PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
265:   rstart = N->cmap->rstart;
266:   rend   = N->cmap->rend;
267:   PetscCall(VecGetArray(v, &values));
268:   PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
269:   PetscCall(VecRestoreArray(v, &values));
270:   PetscCall(PetscFree2(diag, work));
271:   PetscCall(VecScale(v, Na->scale));
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
276: {
277:   Mat_Normal *Na = (Mat_Normal *)N->data;
278:   Mat         M, A = Na->A;

280:   PetscFunctionBegin;
281:   PetscCall(MatGetDiagonalBlock(A, &M));
282:   PetscCall(MatCreateNormal(M, &Na->D));
283:   *D = Na->D;
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
288: {
289:   Mat_Normal *Aa = (Mat_Normal *)A->data;

291:   PetscFunctionBegin;
292:   *M = Aa->A;
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:       MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`

299:    Logically Collective

301:    Input Parameter:
302: .   A  - the `MATNORMAL` matrix

304:    Output Parameter:
305: .   M - the matrix object stored inside `A`

307:    Level: intermediate

309: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
310: @*/
311: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
312: {
313:   PetscFunctionBegin;
317:   PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
318:   PetscFunctionReturn(PETSC_SUCCESS);
319: }

321: PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
322: {
323:   Mat_Normal *Aa = (Mat_Normal *)A->data;
324:   Mat         B;
325:   PetscInt    m, n, M, N;

327:   PetscFunctionBegin;
328:   PetscCall(MatGetSize(A, &M, &N));
329:   PetscCall(MatGetLocalSize(A, &m, &n));
330:   if (reuse == MAT_REUSE_MATRIX) {
331:     B = *newmat;
332:     PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
333:   } else {
334:     PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
335:     PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
336:     PetscCall(MatProductSetFromOptions(B));
337:     PetscCall(MatProductSymbolic(B));
338:     PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
339:   }
340:   PetscCall(MatProductNumeric(B));
341:   if (reuse == MAT_INPLACE_MATRIX) {
342:     PetscCall(MatHeaderReplace(A, &B));
343:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
344:   PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
345:   PetscFunctionReturn(PETSC_SUCCESS);
346: }

348: #if defined(PETSC_HAVE_HYPRE)
349: PetscErrorCode MatConvert_Normal_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
350: {
351:   PetscFunctionBegin;
352:   if (reuse == MAT_INITIAL_MATRIX) {
353:     PetscCall(MatConvert(A, MATAIJ, reuse, B));
354:     PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
355:   } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }
358: #endif

360: typedef struct {
361:   Mat work[2];
362: } Normal_Dense;

364: PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
365: {
366:   Mat           A, B;
367:   Normal_Dense *contents;
368:   Mat_Normal   *a;
369:   PetscScalar  *array;

371:   PetscFunctionBegin;
372:   MatCheckProduct(C, 1);
373:   A        = C->product->A;
374:   a        = (Mat_Normal *)A->data;
375:   B        = C->product->B;
376:   contents = (Normal_Dense *)C->product->data;
377:   PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
378:   if (a->right) {
379:     PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
380:     PetscCall(MatDiagonalScale(C, a->right, NULL));
381:   }
382:   PetscCall(MatProductNumeric(contents->work[0]));
383:   PetscCall(MatDenseGetArrayWrite(C, &array));
384:   PetscCall(MatDensePlaceArray(contents->work[1], array));
385:   PetscCall(MatProductNumeric(contents->work[1]));
386:   PetscCall(MatDenseRestoreArrayWrite(C, &array));
387:   PetscCall(MatDenseResetArray(contents->work[1]));
388:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
389:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
390:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
391:   PetscCall(MatScale(C, a->scale));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: PetscErrorCode MatNormal_DenseDestroy(void *ctx)
396: {
397:   Normal_Dense *contents = (Normal_Dense *)ctx;

399:   PetscFunctionBegin;
400:   PetscCall(MatDestroy(contents->work));
401:   PetscCall(MatDestroy(contents->work + 1));
402:   PetscCall(PetscFree(contents));
403:   PetscFunctionReturn(PETSC_SUCCESS);
404: }

406: PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
407: {
408:   Mat           A, B;
409:   Normal_Dense *contents = NULL;
410:   Mat_Normal   *a;
411:   PetscScalar  *array;
412:   PetscInt      n, N, m, M;

414:   PetscFunctionBegin;
415:   MatCheckProduct(C, 1);
416:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
417:   A = C->product->A;
418:   a = (Mat_Normal *)A->data;
419:   PetscCheck(!a->left, PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "Not implemented");
420:   B = C->product->B;
421:   PetscCall(MatGetLocalSize(C, &m, &n));
422:   PetscCall(MatGetSize(C, &M, &N));
423:   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
424:     PetscCall(MatGetLocalSize(B, NULL, &n));
425:     PetscCall(MatGetSize(B, NULL, &N));
426:     PetscCall(MatGetLocalSize(A, &m, NULL));
427:     PetscCall(MatGetSize(A, &M, NULL));
428:     PetscCall(MatSetSizes(C, m, n, M, N));
429:   }
430:   PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
431:   PetscCall(MatSetUp(C));
432:   PetscCall(PetscNew(&contents));
433:   C->product->data    = contents;
434:   C->product->destroy = MatNormal_DenseDestroy;
435:   if (a->right) {
436:     PetscCall(MatProductCreate(a->A, C, NULL, contents->work));
437:   } else {
438:     PetscCall(MatProductCreate(a->A, B, NULL, contents->work));
439:   }
440:   PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB));
441:   PetscCall(MatProductSetFromOptions(contents->work[0]));
442:   PetscCall(MatProductSymbolic(contents->work[0]));
443:   PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1));
444:   PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB));
445:   PetscCall(MatProductSetFromOptions(contents->work[1]));
446:   PetscCall(MatProductSymbolic(contents->work[1]));
447:   PetscCall(MatDenseGetArrayWrite(C, &array));
448:   PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array));
449:   PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array));
450:   PetscCall(MatDenseRestoreArrayWrite(C, &array));
451:   C->ops->productnumeric = MatProductNumeric_Normal_Dense;
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

455: PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
456: {
457:   PetscFunctionBegin;
458:   C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
459:   PetscFunctionReturn(PETSC_SUCCESS);
460: }

462: PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
463: {
464:   Mat_Product *product = C->product;

466:   PetscFunctionBegin;
467:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C));
468:   PetscFunctionReturn(PETSC_SUCCESS);
469: }

471: /*MC
472:   MATNORMAL - a matrix that behaves like A'*A for `MatMult()` while only containing A

474:   Level: intermediate

476: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
477: M*/

479: /*@
480:       MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A.

482:    Collective

484:    Input Parameter:
485: .   A  - the (possibly rectangular) matrix

487:    Output Parameter:
488: .   N - the matrix that represents A'*A

490:    Level: intermediate

492:    Notes:
493:     The product A'*A is NOT actually formed! Rather the new matrix
494:           object performs the matrix-vector product, `MatMult()`, by first multiplying by
495:           A and then A'

497: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
498: @*/
499: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
500: {
501:   PetscInt    n, nn;
502:   Mat_Normal *Na;
503:   VecType     vtype;

505:   PetscFunctionBegin;
506:   PetscCall(MatGetSize(A, NULL, &nn));
507:   PetscCall(MatGetLocalSize(A, NULL, &n));
508:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
509:   PetscCall(MatSetSizes(*N, n, n, nn, nn));
510:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL));
511:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
512:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));

514:   PetscCall(PetscNew(&Na));
515:   (*N)->data = (void *)Na;
516:   PetscCall(PetscObjectReference((PetscObject)A));
517:   Na->A     = A;
518:   Na->scale = 1.0;

520:   PetscCall(MatCreateVecs(A, NULL, &Na->w));

522:   (*N)->ops->destroy           = MatDestroy_Normal;
523:   (*N)->ops->mult              = MatMult_Normal;
524:   (*N)->ops->multtranspose     = MatMultTranspose_Normal;
525:   (*N)->ops->multtransposeadd  = MatMultTransposeAdd_Normal;
526:   (*N)->ops->multadd           = MatMultAdd_Normal;
527:   (*N)->ops->getdiagonal       = MatGetDiagonal_Normal;
528:   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_Normal;
529:   (*N)->ops->scale             = MatScale_Normal;
530:   (*N)->ops->diagonalscale     = MatDiagonalScale_Normal;
531:   (*N)->ops->increaseoverlap   = MatIncreaseOverlap_Normal;
532:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
533:   (*N)->ops->permute           = MatPermute_Normal;
534:   (*N)->ops->duplicate         = MatDuplicate_Normal;
535:   (*N)->ops->copy              = MatCopy_Normal;
536:   (*N)->assembled              = PETSC_TRUE;
537:   (*N)->preallocated           = PETSC_TRUE;

539:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMat_C", MatNormalGetMat_Normal));
540:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ));
541:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ));
542: #if defined(PETSC_HAVE_HYPRE)
543:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_hypre_C", MatConvert_Normal_HYPRE));
544: #endif
545:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense));
546:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense));
547:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_dense_C", MatProductSetFromOptions_Normal_Dense));
548:   PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE));
549:   PetscCall(MatGetVecType(A, &vtype));
550:   PetscCall(MatSetVecType(*N, vtype));
551: #if defined(PETSC_HAVE_DEVICE)
552:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
553: #endif
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }