Actual source code: normm.c

  1: #include <../src/mat/impls/shell/shell.h>

  3: typedef struct {
  4:   Mat A;
  5:   Mat D; /* local submatrix for diagonal part */
  6:   Vec w;
  7: } Mat_Normal;

  9: static PetscErrorCode MatIncreaseOverlap_Normal(Mat A, PetscInt is_max, IS is[], PetscInt ov)
 10: {
 11:   Mat_Normal *a;
 12:   Mat         pattern;

 14:   PetscFunctionBegin;
 15:   PetscCheck(ov >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap specified");
 16:   PetscCall(MatShellGetContext(A, &a));
 17:   PetscCall(MatProductCreate(a->A, a->A, NULL, &pattern));
 18:   PetscCall(MatProductSetType(pattern, MATPRODUCT_AtB));
 19:   PetscCall(MatProductSetFromOptions(pattern));
 20:   PetscCall(MatProductSymbolic(pattern));
 21:   PetscCall(MatIncreaseOverlap(pattern, is_max, is, ov));
 22:   PetscCall(MatDestroy(&pattern));
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

 26: static PetscErrorCode MatCreateSubMatrices_Normal(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
 27: {
 28:   Mat_Normal *a;
 29:   Mat         B, *suba;
 30:   IS         *row;
 31:   PetscInt    M;

 33:   PetscFunctionBegin;
 34:   PetscCheck(!((Mat_Shell *)mat->data)->zrows && !((Mat_Shell *)mat->data)->zcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatZeroRows()/MatZeroRowsColumns() after the SubMatrices creation
 35:   PetscCheck(!((Mat_Shell *)mat->data)->axpy, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatAXPY() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatAXPY() after the SubMatrices creation
 36:   PetscCheck(!((Mat_Shell *)mat->data)->left && !((Mat_Shell *)mat->data)->right, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalScale() after the SubMatrices creation with a SubVector
 37:   PetscCheck(!((Mat_Shell *)mat->data)->dshift, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalSet() after the SubMatrices creation with a SubVector
 38:   PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
 39:   PetscCall(MatShellGetContext(mat, &a));
 40:   B = a->A;
 41:   if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
 42:   PetscCall(MatGetSize(B, &M, NULL));
 43:   PetscCall(PetscMalloc1(n, &row));
 44:   PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
 45:   PetscCall(ISSetIdentity(row[0]));
 46:   for (M = 1; M < n; ++M) row[M] = row[0];
 47:   PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
 48:   for (M = 0; M < n; ++M) {
 49:     PetscCall(MatCreateNormal(suba[M], *submat + M));
 50:     ((Mat_Shell *)(*submat)[M]->data)->vscale = ((Mat_Shell *)mat->data)->vscale;
 51:     ((Mat_Shell *)(*submat)[M]->data)->vshift = ((Mat_Shell *)mat->data)->vshift;
 52:   }
 53:   PetscCall(ISDestroy(&row[0]));
 54:   PetscCall(PetscFree(row));
 55:   PetscCall(MatDestroySubMatrices(n, &suba));
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }

 59: static PetscErrorCode MatPermute_Normal(Mat A, IS rowp, IS colp, Mat *B)
 60: {
 61:   Mat_Normal *a;
 62:   Mat         C, Aa;
 63:   IS          row;

 65:   PetscFunctionBegin;
 66:   PetscCheck(!((Mat_Shell *)A->data)->zrows && !((Mat_Shell *)A->data)->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatPermute() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatZeroRows()/MatZeroRowsColumns() after the permutation with a permuted zrows and zcols
 67:   PetscCheck(!((Mat_Shell *)A->data)->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatPermute() if MatAXPY() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatAXPY() after the permutation with a permuted axpy
 68:   PetscCheck(!((Mat_Shell *)A->data)->left && !((Mat_Shell *)A->data)->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatPermute() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalScale() after the permutation with a permuted left and right
 69:   PetscCheck(!((Mat_Shell *)A->data)->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatPermute() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalSet() after the permutation with a permuted dshift
 70:   PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
 71:   PetscCall(MatShellGetContext(A, &a));
 72:   Aa = a->A;
 73:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
 74:   PetscCall(ISSetIdentity(row));
 75:   PetscCall(MatPermute(Aa, row, colp, &C));
 76:   PetscCall(ISDestroy(&row));
 77:   PetscCall(MatCreateNormal(C, B));
 78:   PetscCall(MatDestroy(&C));
 79:   ((Mat_Shell *)(*B)->data)->vscale = ((Mat_Shell *)A->data)->vscale;
 80:   ((Mat_Shell *)(*B)->data)->vshift = ((Mat_Shell *)A->data)->vshift;
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: static PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
 85: {
 86:   Mat_Normal *a;
 87:   Mat         C;

 89:   PetscFunctionBegin;
 90:   PetscCall(MatShellGetContext(A, &a));
 91:   PetscCall(MatDuplicate(a->A, op, &C));
 92:   PetscCall(MatCreateNormal(C, B));
 93:   PetscCall(MatDestroy(&C));
 94:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
 95:   PetscFunctionReturn(PETSC_SUCCESS);
 96: }

 98: static PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str)
 99: {
100:   Mat_Normal *a, *b;

102:   PetscFunctionBegin;
103:   PetscCall(MatShellGetContext(A, &a));
104:   PetscCall(MatShellGetContext(B, &b));
105:   PetscCall(MatCopy(a->A, b->A, str));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: static PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y)
110: {
111:   Mat_Normal *Na;

113:   PetscFunctionBegin;
114:   PetscCall(MatShellGetContext(N, &Na));
115:   PetscCall(MatMult(Na->A, x, Na->w));
116:   PetscCall(MatMultTranspose(Na->A, Na->w, y));
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: static PetscErrorCode MatDestroy_Normal(Mat N)
121: {
122:   Mat_Normal *Na;

124:   PetscFunctionBegin;
125:   PetscCall(MatShellGetContext(N, &Na));
126:   PetscCall(MatDestroy(&Na->A));
127:   PetscCall(MatDestroy(&Na->D));
128:   PetscCall(VecDestroy(&Na->w));
129:   PetscCall(PetscFree(Na));
130:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL));
131:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_seqaij_C", NULL));
132:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_mpiaij_C", NULL));
133: #if defined(PETSC_HAVE_HYPRE)
134:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_hypre_C", NULL));
135: #endif
136:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL));
137:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL));
138:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: /*
143:       Slow, nonscalable version
144: */
145: static PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
146: {
147:   Mat_Normal        *Na;
148:   Mat                A;
149:   PetscInt           i, j, rstart, rend, nnz;
150:   const PetscInt    *cols;
151:   PetscScalar       *diag, *work, *values;
152:   const PetscScalar *mvalues;

154:   PetscFunctionBegin;
155:   PetscCall(MatShellGetContext(N, &Na));
156:   A = Na->A;
157:   PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
158:   PetscCall(PetscArrayzero(work, A->cmap->N));
159:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
160:   for (i = rstart; i < rend; i++) {
161:     PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
162:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
163:     PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
164:   }
165:   PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
166:   rstart = N->cmap->rstart;
167:   rend   = N->cmap->rend;
168:   PetscCall(VecGetArray(v, &values));
169:   PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
170:   PetscCall(VecRestoreArray(v, &values));
171:   PetscCall(PetscFree2(diag, work));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: static PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
176: {
177:   Mat_Normal *Na;
178:   Mat         M, A;

180:   PetscFunctionBegin;
181:   PetscCall(MatShellGetContext(N, &Na));
182:   A = Na->A;
183:   PetscCall(MatGetDiagonalBlock(A, &M));
184:   PetscCall(MatCreateNormal(M, &Na->D));
185:   *D = Na->D;
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: static PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
190: {
191:   Mat_Normal *Aa;

193:   PetscFunctionBegin;
194:   PetscCall(MatShellGetContext(A, &Aa));
195:   *M = Aa->A;
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: /*@
200:   MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`

202:   Logically Collective

204:   Input Parameter:
205: . A - the `MATNORMAL` matrix

207:   Output Parameter:
208: . M - the matrix object stored inside `A`

210:   Level: intermediate

212: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
213: @*/
214: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
215: {
216:   PetscFunctionBegin;
219:   PetscAssertPointer(M, 2);
220:   PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
221:   PetscFunctionReturn(PETSC_SUCCESS);
222: }

224: static PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
225: {
226:   Mat_Normal *Aa;
227:   Mat         B;
228:   PetscInt    m, n, M, N;

230:   PetscFunctionBegin;
231:   PetscCall(MatShellGetContext(A, &Aa));
232:   PetscCall(MatGetSize(A, &M, &N));
233:   PetscCall(MatGetLocalSize(A, &m, &n));
234:   if (reuse == MAT_REUSE_MATRIX) {
235:     B = *newmat;
236:     PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
237:   } else {
238:     PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
239:     PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
240:     PetscCall(MatProductSetFromOptions(B));
241:     PetscCall(MatProductSymbolic(B));
242:     PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
243:   }
244:   PetscCall(MatProductNumeric(B));
245:   if (reuse == MAT_INPLACE_MATRIX) {
246:     PetscCall(MatHeaderReplace(A, &B));
247:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
248:   PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: #if defined(PETSC_HAVE_HYPRE)
253: static PetscErrorCode MatConvert_Normal_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
254: {
255:   PetscFunctionBegin;
256:   if (reuse == MAT_INITIAL_MATRIX) {
257:     PetscCall(MatConvert(A, MATAIJ, reuse, B));
258:     PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
259:   } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }
262: #endif

264: typedef struct {
265:   Mat work[2];
266: } Normal_Dense;

268: static PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
269: {
270:   Mat           A, B;
271:   Normal_Dense *contents;
272:   Mat_Normal   *a;
273:   PetscScalar  *array;

275:   PetscFunctionBegin;
276:   MatCheckProduct(C, 1);
277:   A = C->product->A;
278:   B = C->product->B;
279:   PetscCall(MatShellGetContext(A, &a));
280:   contents = (Normal_Dense *)C->product->data;
281:   PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
282:   if (((Mat_Shell *)A->data)->right) {
283:     PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
284:     PetscCall(MatDiagonalScale(C, ((Mat_Shell *)A->data)->right, NULL));
285:   }
286:   PetscCall(MatProductNumeric(contents->work[0]));
287:   PetscCall(MatDenseGetArrayWrite(C, &array));
288:   PetscCall(MatDensePlaceArray(contents->work[1], array));
289:   PetscCall(MatProductNumeric(contents->work[1]));
290:   PetscCall(MatDenseRestoreArrayWrite(C, &array));
291:   PetscCall(MatDenseResetArray(contents->work[1]));
292:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
293:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
294:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
295:   PetscCall(MatScale(C, ((Mat_Shell *)A->data)->vscale));
296:   PetscFunctionReturn(PETSC_SUCCESS);
297: }

299: static PetscErrorCode MatNormal_DenseDestroy(void *ctx)
300: {
301:   Normal_Dense *contents = (Normal_Dense *)ctx;

303:   PetscFunctionBegin;
304:   PetscCall(MatDestroy(contents->work));
305:   PetscCall(MatDestroy(contents->work + 1));
306:   PetscCall(PetscFree(contents));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: static PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
311: {
312:   Mat           A, B;
313:   Normal_Dense *contents = NULL;
314:   Mat_Normal   *a;
315:   PetscScalar  *array;
316:   PetscInt      n, N, m, M;

318:   PetscFunctionBegin;
319:   MatCheckProduct(C, 1);
320:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
321:   A = C->product->A;
322:   B = C->product->B;
323:   PetscCall(MatShellGetContext(A, &a));
324:   PetscCheck(!((Mat_Shell *)A->data)->zrows && !((Mat_Shell *)A->data)->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatProductSymbolic() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
325:   PetscCheck(!((Mat_Shell *)A->data)->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatProductSymbolic() if MatAXPY() has been called on the input Mat");          // TODO FIXME
326:   PetscCheck(!((Mat_Shell *)A->data)->left, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatProductSymbolic() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME
327:   PetscCheck(!((Mat_Shell *)A->data)->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatProductSymbolic() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME
328:   PetscCall(MatGetLocalSize(C, &m, &n));
329:   PetscCall(MatGetSize(C, &M, &N));
330:   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
331:     PetscCall(MatGetLocalSize(B, NULL, &n));
332:     PetscCall(MatGetSize(B, NULL, &N));
333:     PetscCall(MatGetLocalSize(A, &m, NULL));
334:     PetscCall(MatGetSize(A, &M, NULL));
335:     PetscCall(MatSetSizes(C, m, n, M, N));
336:   }
337:   PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
338:   PetscCall(MatSetUp(C));
339:   PetscCall(PetscNew(&contents));
340:   C->product->data    = contents;
341:   C->product->destroy = MatNormal_DenseDestroy;
342:   if (((Mat_Shell *)A->data)->right) {
343:     PetscCall(MatProductCreate(a->A, C, NULL, contents->work));
344:   } else {
345:     PetscCall(MatProductCreate(a->A, B, NULL, contents->work));
346:   }
347:   PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB));
348:   PetscCall(MatProductSetFromOptions(contents->work[0]));
349:   PetscCall(MatProductSymbolic(contents->work[0]));
350:   PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1));
351:   PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB));
352:   PetscCall(MatProductSetFromOptions(contents->work[1]));
353:   PetscCall(MatProductSymbolic(contents->work[1]));
354:   PetscCall(MatDenseGetArrayWrite(C, &array));
355:   PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array));
356:   PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array));
357:   PetscCall(MatDenseRestoreArrayWrite(C, &array));
358:   C->ops->productnumeric = MatProductNumeric_Normal_Dense;
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: static PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
363: {
364:   PetscFunctionBegin;
365:   C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: static PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
370: {
371:   Mat_Product *product = C->product;

373:   PetscFunctionBegin;
374:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C));
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

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

381:   Level: intermediate

383:   Developer Notes:
384:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

386:   Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage

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

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

394:   Collective

396:   Input Parameter:
397: . A - the (possibly rectangular) matrix

399:   Output Parameter:
400: . N - the matrix that represents A'*A

402:   Level: intermediate

404:   Notes:
405:   The product A'*A is NOT actually formed! Rather the new matrix
406:   object performs the matrix-vector product, `MatMult()`, by first multiplying by
407:   A and then A'

409: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
410: @*/
411: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
412: {
413:   Mat_Normal *Na;
414:   VecType     vtype;

416:   PetscFunctionBegin;
417:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
418:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
419:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
420:   PetscCall(MatSetType(*N, MATSHELL));
421:   PetscCall(PetscNew(&Na));
422:   PetscCall(MatShellSetContext(*N, Na));
423:   PetscCall(PetscObjectReference((PetscObject)A));
424:   Na->A = A;
425:   PetscCall(MatCreateVecs(A, NULL, &Na->w));

427:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
428:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Normal));
429:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Normal));
430:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Normal));
431:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Normal));
432:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Normal));
433:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (void (*)(void))MatGetDiagonalBlock_Normal));
434:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Normal));
435:   (*N)->ops->increaseoverlap   = MatIncreaseOverlap_Normal;
436:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
437:   (*N)->ops->permute           = MatPermute_Normal;

439:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalGetMat_Normal));
440:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ));
441:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ));
442: #if defined(PETSC_HAVE_HYPRE)
443:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_hypre_C", MatConvert_Normal_HYPRE));
444: #endif
445:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense));
446:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense));
447:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
448:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
449:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
450:   PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE));
451:   PetscCall(MatGetVecType(A, &vtype));
452:   PetscCall(MatSetVecType(*N, vtype));
453: #if defined(PETSC_HAVE_DEVICE)
454:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
455: #endif
456:   PetscCall(MatSetUp(*N));
457:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL));
458:   PetscFunctionReturn(PETSC_SUCCESS);
459: }