Actual source code: ex321.c
1: static char help[] = "Tests that the assembled views of a MATIS follow a change made through MatISGetLocalMat().\n\n";
3: #include <petscmat.h>
5: // Prints the diagonal of A, one line per process.
6: static PetscErrorCode DiagonalView(Mat A, const char *label)
7: {
8: Vec d;
9: PetscInt n;
10: const PetscScalar *vals;
12: PetscFunctionBeginUser;
13: PetscCall(MatCreateVecs(A, NULL, &d));
14: PetscCall(MatGetDiagonal(A, d));
15: PetscCall(VecGetLocalSize(d, &n));
16: PetscCall(VecGetArrayRead(d, &vals));
17: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%s", label));
18: for (PetscInt i = 0; i < n; i++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, " %g", (double)PetscRealPart(vals[i])));
19: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n"));
20: PetscCall(VecRestoreArrayRead(d, &vals));
21: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
22: PetscCall(VecDestroy(&d));
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode CheckViews(Mat A, PetscBool subfirst)
27: {
28: Mat B, dA, dB, *sub;
29: MatState before, after;
30: IS rows;
31: PetscInt start, end;
32: PetscBool equal;
34: PetscFunctionBeginUser;
35: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
36: PetscCall(MatGetDiagonalBlock(B, &dB));
37: PetscCall(MatGetOwnershipRange(A, &start, &end));
38: PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &rows));
39: if (!subfirst) PetscCall(MatGetDiagonalBlock(A, &dA));
40: PetscCall(MatCreateSubMatrices(A, 1, &rows, &rows, MAT_INITIAL_MATRIX, &sub));
41: if (subfirst) PetscCall(MatGetDiagonalBlock(A, &dA));
42: PetscCall(MatMultEqual(dA, dB, 3, &equal));
43: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stale diagonal block");
44: PetscCall(MatMultEqual(sub[0], dB, 3, &equal));
45: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stale assembled submatrix");
46: PetscCall(MatGetState(dA, &before));
47: PetscCall(MatGetDiagonalBlock(A, &dA));
48: PetscCall(MatGetState(dA, &after));
49: PetscCall(MatStateCompare(before, after, &equal));
50: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "An unchanged diagonal block was rebuilt");
51: PetscCall(MatDestroySubMatrices(1, &sub));
52: PetscCall(ISDestroy(&rows));
53: PetscCall(MatDestroy(&B));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode TestState(Mat A)
58: {
59: Mat lA, B, C, D;
60: MatState before, after, lbefore, lafter;
61: PetscInt row = 0, n;
62: PetscMPIInt rank;
63: PetscBool same;
65: PetscFunctionBeginUser;
66: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
67: PetscCall(CheckViews(A, PETSC_FALSE));
68: PetscCall(MatGetState(A, &before));
69: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
70: PetscCall(MatAXPY(A, 0.5, B, SAME_NONZERO_PATTERN));
71: PetscCall(MatGetState(A, &after));
72: PetscCheck(after.state > before.state && after.nonzerostate == before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAXPY() did not propagate a numerical change");
73: PetscCall(CheckViews(A, PETSC_FALSE));
74: before = after;
76: PetscCall(MatCopy(B, A, SAME_NONZERO_PATTERN));
77: PetscCall(MatGetState(A, &after));
78: PetscCheck(after.state > before.state && after.nonzerostate == before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy() did not propagate a numerical change");
79: PetscCall(MatMultEqual(A, B, 3, &same));
80: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatCopy() with the same nonzero pattern produced incorrect values");
81: PetscCall(CheckViews(A, PETSC_TRUE));
82: PetscCall(MatDestroy(&B));
83: before = after;
85: if (!rank) {
86: PetscCall(MatISGetLocalMat(A, &lA));
87: PetscCall(MatScale(lA, 2.0));
88: PetscCall(MatISRestoreLocalMat(A, &lA));
89: }
90: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
91: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
92: PetscCall(MatGetState(A, &after));
93: PetscCheck(after.state > before.state && after.nonzerostate == before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A numerical change was not propagated correctly");
94: PetscCall(CheckViews(A, PETSC_TRUE));
96: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
97: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
98: PetscCall(MatSetOption(C, MAT_KEEP_NONZERO_PATTERN, PETSC_FALSE));
99: PetscCall(MatZeroRows(C, rank ? 0 : 1, &row, 1.0, NULL, NULL));
100: before = after;
101: if (!rank) {
102: PetscCall(MatISGetLocalMat(A, &lA));
103: PetscCall(MatSetOption(lA, MAT_KEEP_NONZERO_PATTERN, PETSC_FALSE));
104: PetscCall(MatZeroRows(lA, 1, &row, 1.0, NULL, NULL));
105: PetscCall(MatISRestoreLocalMat(A, &lA));
106: }
107: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
108: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
109: PetscCall(MatGetState(A, &after));
110: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A structural change on one rank was not propagated");
111: PetscCall(MatMultEqual(A, C, 3, &same));
112: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Zeroing a local row without keeping its nonzero pattern produced incorrect values");
113: PetscCall(CheckViews(A, PETSC_FALSE));
114: PetscCall(MatDestroy(&C));
116: // Restore the dropped entries through MatCopy(), with both views already cached.
117: before = after;
118: PetscCall(MatCopy(B, A, DIFFERENT_NONZERO_PATTERN));
119: PetscCall(MatGetState(A, &after));
120: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy() did not propagate a structural change");
121: PetscCall(MatMultEqual(A, B, 3, &same));
122: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatCopy() with different nonzero patterns produced incorrect values");
123: PetscCall(CheckViews(A, PETSC_TRUE));
125: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &C));
126: PetscCall(MatSetOption(C, MAT_KEEP_NONZERO_PATTERN, PETSC_FALSE));
127: PetscCall(MatZeroRows(C, rank ? 0 : 1, &row, 1.0, NULL, NULL));
128: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_FALSE));
129: before = after;
130: PetscCall(MatZeroRows(A, rank ? 0 : 1, &row, 1.0, NULL, NULL));
131: PetscCall(MatGetState(A, &after));
132: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatZeroRows() did not propagate a structural change");
133: PetscCall(MatMultEqual(A, C, 3, &same));
134: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatZeroRows() without keeping the nonzero pattern produced incorrect values");
135: PetscCall(CheckViews(A, PETSC_FALSE));
137: // MatAXPY() must add the missing entries and invalidate both cached views.
138: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &D));
139: PetscCall(MatAXPY(C, 0.5, D, DIFFERENT_NONZERO_PATTERN));
140: before = after;
141: PetscCall(MatAXPY(A, 0.5, B, DIFFERENT_NONZERO_PATTERN));
142: PetscCall(MatGetState(A, &after));
143: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAXPY() did not propagate a structural change");
144: PetscCall(MatMultEqual(A, C, 3, &same));
145: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatAXPY() with different nonzero patterns produced incorrect values");
146: PetscCall(CheckViews(A, PETSC_TRUE));
147: PetscCall(MatDestroy(&D));
148: PetscCall(MatDestroy(&C));
149: PetscCall(MatDestroy(&B));
151: // Replace the local matrix on one rank, retaining its type and nonzero pattern.
152: for (PetscInt i = 0; i < 2; i++) {
153: PetscCall(MatGetState(A, &before));
154: if (!rank) {
155: PetscCall(MatISGetLocalMat(A, &lA));
156: PetscCall(MatDuplicate(lA, MAT_COPY_VALUES, &B));
157: PetscCall(MatISRestoreLocalMat(A, &lA));
158: PetscCall(MatScale(B, 2.0));
159: PetscCall(MatISSetLocalMat(A, B));
160: PetscCall(MatDestroy(&B));
161: }
162: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
163: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
164: PetscCall(MatGetState(A, &after));
165: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local matrix replacement was not propagated");
166: PetscCall(CheckViews(A, (PetscBool)i));
167: }
169: // Replacing an empty local matrix must be detected even when its nonzero state is unchanged.
170: for (PetscInt i = 0; i < 2; i++) {
171: PetscCall(MatGetState(A, &before));
172: if (!rank) {
173: PetscCall(MatISGetLocalMat(A, &lA));
174: PetscCall(MatGetState(lA, &lbefore));
175: PetscCall(MatGetSize(lA, &n, NULL));
176: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 0, NULL, &B));
177: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
178: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
179: PetscCall(MatGetState(B, &lafter));
180: if (i) PetscCheck(lbefore.id != lafter.id && lbefore.nonzerostate == lafter.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected different local matrices with equal nonzero states");
181: PetscCall(MatISRestoreLocalMat(A, &lA));
182: PetscCall(MatISSetLocalMat(A, B));
183: PetscCall(MatDestroy(&B));
184: }
185: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
187: PetscCall(MatGetState(A, &after));
188: PetscCheck(after.state > before.state && after.nonzerostate > before.nonzerostate, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local matrix replacement was not propagated");
189: PetscCall(CheckViews(A, PETSC_TRUE));
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode TestReuse(Mat A)
195: {
196: Mat C;
197: PetscBool equal;
199: PetscFunctionBeginUser;
200: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &C));
201: PetscCall(MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &C));
202: PetscCall(MatMultEqual(A, C, 3, &equal));
203: PetscCheck(equal, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Reuse with the same source failed");
204: PetscCall(MatDestroy(&C));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: int main(int argc, char **args)
209: {
210: Mat A, Aij, dA, lA;
211: ISLocalToGlobalMapping l2g;
212: PetscInt *gidx;
213: PetscInt row = 1, nl = 4, N;
214: PetscMPIInt rank, size;
215: PetscBool teststate = PETSC_FALSE;
216: const PetscScalar elem[] = {2.0, -1.0, -1.0, 2.0};
218: PetscFunctionBeginUser;
219: PetscCall(PetscInitialize(&argc, &args, NULL, help));
220: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
221: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
222: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_state", &teststate, NULL));
224: // Subdomain r holds the nl nodes r*(nl-1) ... r*(nl-1)+nl-1, so consecutive subdomains share a node.
225: N = size * (nl - 1) + 1;
226: PetscCall(PetscMalloc1(nl, &gidx));
227: for (PetscInt i = 0; i < nl; i++) gidx[i] = rank * (nl - 1) + i;
228: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, nl, gidx, PETSC_OWN_POINTER, &l2g));
229: PetscCall(MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE, N, N, l2g, l2g, &A));
230: PetscCall(MatSetFromOptions(A));
231: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
232: PetscCall(MatISSetPreallocation(A, 3, NULL, 0, NULL));
233: for (PetscInt e = 0; e < nl - 1; e++) {
234: const PetscInt erows[] = {e, e + 1};
236: PetscCall(MatSetValuesLocal(A, 2, erows, 2, erows, elem, ADD_VALUES));
237: }
238: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
239: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
241: // Assembled views built before the change, the way a solver caches them.
242: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Aij));
243: PetscCall(MatGetDiagonalBlock(A, &dA));
245: // Impose a unit diagonal on one row of every local matrix, keeping the nonzero pattern.
246: PetscCall(MatISGetLocalMat(A, &lA));
247: PetscCall(MatSetOption(lA, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
248: PetscCall(MatZeroRows(lA, 1, &row, 1.0, NULL, NULL));
249: PetscCall(MatISRestoreLocalMat(A, &lA));
250: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
251: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
253: PetscCall(MatISGetLocalMat(A, &lA));
254: PetscCall(DiagonalView(lA, "local (Neumann):"));
255: PetscCall(MatISRestoreLocalMat(A, &lA));
256: PetscCall(MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &Aij));
257: PetscCall(DiagonalView(Aij, "assembled (MAT_REUSE_MATRIX):"));
258: PetscCall(MatDestroy(&Aij));
259: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Aij));
260: PetscCall(DiagonalView(Aij, "assembled (MAT_INITIAL_MATRIX):"));
261: PetscCall(MatGetDiagonalBlock(A, &dA));
262: PetscCall(DiagonalView(dA, "diagonal block:"));
264: if (teststate) {
265: PetscCall(TestReuse(A));
266: PetscCall(TestState(A));
267: }
269: PetscCall(MatDestroy(&Aij));
270: PetscCall(MatDestroy(&A));
271: PetscCall(PetscFinalize());
272: return 0;
273: }
275: /*TEST
277: test:
278: nsize: 2
279: diff_args: -j
281: test:
282: suffix: state
283: nsize: 2
284: args: -test_state -mat_is_keepassembled {{0 1}}
285: output_file: output/ex321_1.out
286: diff_args: -j
288: TEST*/