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: PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
182: PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatAXPY() has been called on the input Mat"); // TODO FIXME
183: PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME
184: PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME
185: PetscCall(MatShellGetContext(N, &Na));
186: A = Na->A;
187: PetscCall(MatGetDiagonalBlock(A, &M));
188: PetscCall(MatCreateNormal(M, &Na->D));
189: *D = Na->D;
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
194: {
195: Mat_Normal *Aa;
197: PetscFunctionBegin;
198: PetscCall(MatShellGetContext(A, &Aa));
199: *M = Aa->A;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`
206: Logically Collective
208: Input Parameter:
209: . A - the `MATNORMAL` matrix
211: Output Parameter:
212: . M - the matrix object stored inside `A`
214: Level: intermediate
216: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
217: @*/
218: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
219: {
220: PetscFunctionBegin;
223: PetscAssertPointer(M, 2);
224: PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
229: {
230: Mat_Normal *Aa;
231: Mat B;
232: PetscInt m, n, M, N;
234: PetscFunctionBegin;
235: PetscCall(MatShellGetContext(A, &Aa));
236: PetscCall(MatGetSize(A, &M, &N));
237: PetscCall(MatGetLocalSize(A, &m, &n));
238: if (reuse == MAT_REUSE_MATRIX) {
239: B = *newmat;
240: PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
241: } else {
242: PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
243: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
244: PetscCall(MatProductSetFromOptions(B));
245: PetscCall(MatProductSymbolic(B));
246: PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
247: }
248: PetscCall(MatProductNumeric(B));
249: if (reuse == MAT_INPLACE_MATRIX) {
250: PetscCall(MatHeaderReplace(A, &B));
251: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
252: PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: #if defined(PETSC_HAVE_HYPRE)
257: static PetscErrorCode MatConvert_Normal_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
258: {
259: PetscFunctionBegin;
260: if (reuse == MAT_INITIAL_MATRIX) {
261: PetscCall(MatConvert(A, MATAIJ, reuse, B));
262: PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
263: } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
266: #endif
268: typedef struct {
269: Mat work[2];
270: } Normal_Dense;
272: static PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
273: {
274: Mat A, B;
275: Normal_Dense *contents;
276: Mat_Normal *a;
277: PetscScalar *array;
279: PetscFunctionBegin;
280: MatCheckProduct(C, 1);
281: A = C->product->A;
282: B = C->product->B;
283: PetscCall(MatShellGetContext(A, &a));
284: contents = (Normal_Dense *)C->product->data;
285: PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
286: if (((Mat_Shell *)A->data)->right) {
287: PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
288: PetscCall(MatDiagonalScale(C, ((Mat_Shell *)A->data)->right, NULL));
289: }
290: PetscCall(MatProductNumeric(contents->work[0]));
291: PetscCall(MatDenseGetArrayWrite(C, &array));
292: PetscCall(MatDensePlaceArray(contents->work[1], array));
293: PetscCall(MatProductNumeric(contents->work[1]));
294: PetscCall(MatDenseRestoreArrayWrite(C, &array));
295: PetscCall(MatDenseResetArray(contents->work[1]));
296: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
297: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
298: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
299: PetscCall(MatScale(C, ((Mat_Shell *)A->data)->vscale));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode MatNormal_DenseDestroy(void *ctx)
304: {
305: Normal_Dense *contents = (Normal_Dense *)ctx;
307: PetscFunctionBegin;
308: PetscCall(MatDestroy(contents->work));
309: PetscCall(MatDestroy(contents->work + 1));
310: PetscCall(PetscFree(contents));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: static PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
315: {
316: Mat A, B;
317: Normal_Dense *contents = NULL;
318: Mat_Normal *a;
319: PetscScalar *array;
320: PetscInt n, N, m, M;
322: PetscFunctionBegin;
323: MatCheckProduct(C, 1);
324: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
325: A = C->product->A;
326: B = C->product->B;
327: PetscCall(MatShellGetContext(A, &a));
328: 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
329: 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
330: 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
331: 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
332: PetscCall(MatGetLocalSize(C, &m, &n));
333: PetscCall(MatGetSize(C, &M, &N));
334: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
335: PetscCall(MatGetLocalSize(B, NULL, &n));
336: PetscCall(MatGetSize(B, NULL, &N));
337: PetscCall(MatGetLocalSize(A, &m, NULL));
338: PetscCall(MatGetSize(A, &M, NULL));
339: PetscCall(MatSetSizes(C, m, n, M, N));
340: }
341: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
342: PetscCall(MatSetUp(C));
343: PetscCall(PetscNew(&contents));
344: C->product->data = contents;
345: C->product->destroy = MatNormal_DenseDestroy;
346: if (((Mat_Shell *)A->data)->right) {
347: PetscCall(MatProductCreate(a->A, C, NULL, contents->work));
348: } else {
349: PetscCall(MatProductCreate(a->A, B, NULL, contents->work));
350: }
351: PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB));
352: PetscCall(MatProductSetFromOptions(contents->work[0]));
353: PetscCall(MatProductSymbolic(contents->work[0]));
354: PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1));
355: PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB));
356: PetscCall(MatProductSetFromOptions(contents->work[1]));
357: PetscCall(MatProductSymbolic(contents->work[1]));
358: PetscCall(MatDenseGetArrayWrite(C, &array));
359: PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array));
360: PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array));
361: PetscCall(MatDenseRestoreArrayWrite(C, &array));
362: C->ops->productnumeric = MatProductNumeric_Normal_Dense;
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
367: {
368: PetscFunctionBegin;
369: C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
374: {
375: Mat_Product *product = C->product;
377: PetscFunctionBegin;
378: if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*MC
383: MATNORMAL - a matrix that behaves like A'*A for `MatMult()` while only containing A
385: Level: intermediate
387: Developer Notes:
388: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
390: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
392: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
393: M*/
395: /*@
396: MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A.
398: Collective
400: Input Parameter:
401: . A - the (possibly rectangular) matrix
403: Output Parameter:
404: . N - the matrix that represents A'*A
406: Level: intermediate
408: Notes:
409: The product A'*A is NOT actually formed! Rather the new matrix
410: object performs the matrix-vector product, `MatMult()`, by first multiplying by
411: A and then A'
413: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
414: @*/
415: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
416: {
417: Mat_Normal *Na;
418: VecType vtype;
420: PetscFunctionBegin;
421: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
422: PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
423: PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
424: PetscCall(MatSetType(*N, MATSHELL));
425: PetscCall(PetscNew(&Na));
426: PetscCall(MatShellSetContext(*N, Na));
427: PetscCall(PetscObjectReference((PetscObject)A));
428: Na->A = A;
429: PetscCall(MatCreateVecs(A, NULL, &Na->w));
431: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
432: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Normal));
433: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Normal));
434: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Normal));
435: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Normal));
436: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Normal));
437: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Normal));
438: (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_Normal;
439: (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal;
440: (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
441: (*N)->ops->permute = MatPermute_Normal;
443: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalGetMat_Normal));
444: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ));
445: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ));
446: #if defined(PETSC_HAVE_HYPRE)
447: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_hypre_C", MatConvert_Normal_HYPRE));
448: #endif
449: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense));
450: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense));
451: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
452: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
453: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
454: PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE));
455: PetscCall(MatGetVecType(A, &vtype));
456: PetscCall(MatSetVecType(*N, vtype));
457: #if defined(PETSC_HAVE_DEVICE)
458: PetscCall(MatBindToCPU(*N, A->boundtocpu));
459: #endif
460: PetscCall(MatSetUp(*N));
461: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }