Actual source code: ex315.c
1: static char help[] = "Tests that MATSELL honors MAT_NEW_NONZERO_LOCATIONS and MAT_IGNORE_ZERO_ENTRIES.\n\n";
3: /*
4: MatSetValues() must silently ignore an entry that would create a new nonzero location once
5: MAT_NEW_NONZERO_LOCATIONS is false, and must ignore a zero value once MAT_IGNORE_ZERO_ENTRIES is
6: true. On one process this exercises MatSetValues_SeqSELL(); on two it exercises
7: MatSetValues_MPISELL(), probing the diagonal block, an off-diagonal column that garray already
8: knows, and an off-diagonal column it does not. A MATAIJ matrix receives the same calls and is
9: used as the reference for the resulting values.
11: MAT_IGNORE_ZERO_ENTRIES leaves nothing visible behind, so it is probed indirectly: the ignored
12: zero is followed by a nonzero write to the same location with MAT_NEW_NONZERO_LOCATIONS false. If
13: the zero had wrongly created the location, that second write would land and the values would differ.
15: MAT_IGNORE_ZERO_ENTRIES does not apply to the diagonal, which a zero must still create. Two
16: diagonal locations are left out of the initial structure so the same indirect probe can be run
17: there with the opposite expectation, once with ADD_VALUES and once with INSERT_VALUES, since the
18: guards for the two insert modes are separate.
20: The last probe repeats the off-diagonal one after MatDisAssemble_MPISELL() has replaced the
21: off-diagonal block, which happens once a new off-diagonal column arrives while new locations are
22: still allowed. The replacement block must still honor the option.
24: CheckDiagonalInOffDiagonalBlock() covers the one case where the two rules meet. The diagonal the
25: exception protects is the one MatInvertDiagonalForSOR_SeqSELL() reads, which lives in the
26: diagonal block, so only that block exempts it. A matrix whose row and column layouts differ can
27: place a global (i,i) in the off-diagonal block instead, where it must be dropped like any other
28: zero, and dropped the same way whether or not garray already knows the column.
30: CheckOffDiagonalInDiagonalBlock() covers the reverse, which the same layout produces on the other
31: rank: a global (i,j) with i != j that lands at a block-local (r,r) of the diagonal block. Whether
32: the exemption applies is a question about the global indices, so this zero must be dropped too.
33: */
35: #include <petscmat.h>
37: static PetscErrorCode CheckEntry(Mat A, Mat B, PetscInt row, PetscInt col, const char *where)
38: {
39: PetscScalar va, vb;
41: PetscFunctionBeginUser;
42: PetscCall(MatGetValues(A, 1, &row, 1, &col, &va));
43: PetscCall(MatGetValues(B, 1, &row, 1, &col, &vb));
44: PetscCheck(va == vb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MATSELL has %g at (%" PetscInt_FMT ",%" PetscInt_FMT ") but MATAIJ has %g: %s", (double)PetscRealPart(va), row, col, (double)PetscRealPart(vb), where);
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: /* set the same entry in both matrices, then reassemble both */
49: static PetscErrorCode SetAndAssemble(Mat A, Mat B, PetscInt row, PetscInt col, PetscScalar value, InsertMode addv)
50: {
51: PetscFunctionBeginUser;
52: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, addv));
53: PetscCall(MatSetValues(B, 1, &row, 1, &col, &value, addv));
54: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
55: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
56: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
57: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: /* only the owning rank sets the entry, but the assembly is collective */
62: static PetscErrorCode SetOnOwnerAndAssemble(Mat A, Mat B, PetscBool owner, PetscInt row, PetscInt col, PetscScalar value)
63: {
64: PetscFunctionBeginUser;
65: if (owner) {
66: PetscCall(MatSetValues(A, 1, &row, 1, &col, &value, INSERT_VALUES));
67: PetscCall(MatSetValues(B, 1, &row, 1, &col, &value, INSERT_VALUES));
68: }
69: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
71: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /* An 8 by 8 MATSELL and MATAIJ pair over two ranks whose row and column layouts differ: rank 0 owns
77: rows 0-3 but only columns 0-1, and rank 1 owns rows 4-7 and columns 2-7. A global (2,2) therefore
78: belongs to rank 0's off-diagonal block, while a global (4,2) reaches rank 1 as block-local (0,0)
79: of its diagonal block. Neither location is filled here, so both are free for the probes below. */
80: static PetscErrorCode CreateMismatchedLayoutPair(Mat *A, Mat *B)
81: {
82: PetscInt i, rstart, rend, ncol, col;
83: PetscMPIInt rank;
84: PetscScalar value = 1.0;
86: PetscFunctionBeginUser;
87: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
88: ncol = rank ? 6 : 2;
90: PetscCall(MatCreate(PETSC_COMM_WORLD, A));
91: PetscCall(MatSetSizes(*A, 4, ncol, 8, 8));
92: PetscCall(MatSetType(*A, MATSELL));
93: PetscCall(MatSetFromOptions(*A));
94: PetscCall(MatSeqSELLSetPreallocation(*A, 8, NULL));
95: PetscCall(MatMPISELLSetPreallocation(*A, 8, NULL, 8, NULL));
97: PetscCall(MatCreate(PETSC_COMM_WORLD, B));
98: PetscCall(MatSetSizes(*B, 4, ncol, 8, 8));
99: PetscCall(MatSetType(*B, MATAIJ));
100: PetscCall(MatSeqAIJSetPreallocation(*B, 8, NULL));
101: PetscCall(MatMPIAIJSetPreallocation(*B, 8, NULL, 8, NULL));
103: PetscCall(MatGetOwnershipRange(*A, &rstart, &rend));
104: for (i = rstart; i < rend; i++) {
105: col = (i + 4) % 8;
106: PetscCall(MatSetValues(*A, 1, &i, 1, &col, &value, INSERT_VALUES));
107: PetscCall(MatSetValues(*B, 1, &i, 1, &col, &value, INSERT_VALUES));
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /* A zero at a global (i,i) that the column layout puts in the off-diagonal block must be dropped,
113: not exempted. ingarray selects which of the two suppression paths in MatSetValues_MPISELL() the
114: probe reaches; both must reach the same answer, and it must be the MATAIJ answer. */
115: static PetscErrorCode CheckDiagonalInOffDiagonalBlock(PetscBool ingarray)
116: {
117: Mat A, B;
118: PetscInt i, probe = 2;
119: PetscMPIInt rank;
120: PetscBool owner;
121: PetscScalar value = 1.0;
123: PetscFunctionBeginUser;
124: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
125: owner = (PetscBool)(rank == 0);
126: PetscCall(CreateMismatchedLayoutPair(&A, &B));
127: if (ingarray && owner) {
128: i = 0;
129: PetscCall(MatSetValues(A, 1, &i, 1, &probe, &value, INSERT_VALUES));
130: PetscCall(MatSetValues(B, 1, &i, 1, &probe, &value, INSERT_VALUES));
131: }
132: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
133: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
134: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
135: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
137: /* new locations stay allowed, or nonew would suppress the insert before the zero test is reached */
138: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
139: PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
140: PetscCall(SetOnOwnerAndAssemble(A, B, owner, probe, probe, 0.0));
141: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
142: PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
143: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
144: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
145: PetscCall(SetOnOwnerAndAssemble(A, B, owner, probe, probe, 5.0));
146: if (owner) PetscCall(CheckEntry(A, B, probe, probe, ingarray ? "global diagonal in the off-diagonal block, column in garray" : "global diagonal in the off-diagonal block, column not in garray"));
148: PetscCall(MatDestroy(&A));
149: PetscCall(MatDestroy(&B));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /* The mirror case. A global (4,2) is not on the diagonal, but rank 1 owns rows from 4 and columns
154: from 2, so it arrives at block-local (0,0) of the diagonal block. The exemption belongs to the
155: global diagonal, the one MatInvertDiagonalForSOR_SeqAIJ() and MatInvertDiagonalForSOR_SeqSELL()
156: read, so a zero here must be dropped like any other. */
157: static PetscErrorCode CheckOffDiagonalInDiagonalBlock(void)
158: {
159: Mat A, B;
160: PetscInt prow = 4, pcol = 2;
161: PetscMPIInt rank;
162: PetscBool owner;
164: PetscFunctionBeginUser;
165: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
166: owner = (PetscBool)(rank == 1);
167: PetscCall(CreateMismatchedLayoutPair(&A, &B));
168: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
169: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
170: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
171: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
173: /* new locations stay allowed, or nonew would suppress the insert before the zero test is reached */
174: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
175: PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
176: PetscCall(SetOnOwnerAndAssemble(A, B, owner, prow, pcol, 0.0));
177: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
178: PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
179: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
180: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
181: PetscCall(SetOnOwnerAndAssemble(A, B, owner, prow, pcol, 5.0));
182: if (owner) PetscCall(CheckEntry(A, B, prow, pcol, "global off-diagonal at a block-local diagonal position"));
184: PetscCall(MatDestroy(&A));
185: PetscCall(MatDestroy(&B));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: int main(int argc, char **args)
190: {
191: Mat A, B;
192: PetscInt i, rstart, rend, col, known, unknown, rebuilt, dadd, dins;
193: PetscMPIInt rank, size;
194: PetscScalar value = 1.0;
196: PetscFunctionBeginUser;
197: PetscCall(PetscInitialize(&argc, &args, NULL, help));
198: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
199: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
200: PetscCheck(size == 1 || size == 2, PETSC_COMM_WORLD, PETSC_ERR_USER, "This test requires 1 or 2 processes");
202: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
203: PetscCall(MatSetSizes(A, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE));
204: PetscCall(MatSetType(A, MATSELL));
205: PetscCall(MatSetFromOptions(A));
206: PetscCall(MatSeqSELLSetPreallocation(A, 4, NULL));
207: PetscCall(MatMPISELLSetPreallocation(A, 4, NULL, 4, NULL));
209: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
210: PetscCall(MatSetSizes(B, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE));
211: PetscCall(MatSetType(B, MATAIJ));
212: PetscCall(MatSeqAIJSetPreallocation(B, 4, NULL));
213: PetscCall(MatMPIAIJSetPreallocation(B, 4, NULL, 4, NULL));
215: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
216: /* known enters garray during the first assembly; unknown never does */
217: known = rank ? 0 : 4;
218: unknown = rank ? 1 : 5;
219: /* rebuilt is the off-diagonal column whose arrival forces MatDisAssemble_MPISELL() */
220: rebuilt = rank ? 2 : 6;
221: /* the diagonal is filled except at these two rows, which the zero-value probes below must create */
222: dadd = rstart + 2;
223: dins = rstart + 3;
225: for (i = rstart; i < rend; i++) {
226: if (i == dadd || i == dins) continue;
227: PetscCall(MatSetValues(A, 1, &i, 1, &i, &value, INSERT_VALUES));
228: PetscCall(MatSetValues(B, 1, &i, 1, &i, &value, INSERT_VALUES));
229: }
230: if (size > 1) {
231: PetscCall(MatSetValues(A, 1, &rstart, 1, &known, &value, INSERT_VALUES));
232: PetscCall(MatSetValues(B, 1, &rstart, 1, &known, &value, INSERT_VALUES));
233: }
234: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
236: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
237: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
239: /* from here on, a new nonzero location must be ignored rather than created */
240: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
241: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
243: i = rstart + 1;
244: PetscCall(SetAndAssemble(A, B, i, rstart + 2, 7.0, INSERT_VALUES)); /* diagonal block */
245: PetscCall(CheckEntry(A, B, i, rstart + 2, "new location in the diagonal block"));
246: if (size > 1) {
247: PetscCall(SetAndAssemble(A, B, i, known, 7.0, INSERT_VALUES)); /* off-diagonal, column in garray */
248: PetscCall(CheckEntry(A, B, i, known, "new location in an off-diagonal column garray knows"));
249: PetscCall(SetAndAssemble(A, B, i, unknown, 7.0, INSERT_VALUES)); /* off-diagonal, column not in garray */
250: PetscCall(CheckEntry(A, B, i, unknown, "new location in an off-diagonal column garray does not know"));
251: }
253: /* an existing location must still be writable */
254: PetscCall(SetAndAssemble(A, B, i, i, 9.0, INSERT_VALUES));
255: PetscCall(CheckEntry(A, B, i, i, "overwrite of an existing diagonal entry"));
257: /* a zero value must not create a location either */
258: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
259: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
260: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
261: PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
262: col = rstart + 3;
263: PetscCall(SetAndAssemble(A, B, i, col, 0.0, INSERT_VALUES));
264: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
265: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
266: PetscCall(SetAndAssemble(A, B, i, col, 5.0, INSERT_VALUES));
267: PetscCall(CheckEntry(A, B, i, col, "location that an ignored zero must not have created"));
269: /* but a zero on the diagonal must create its location, for both insert modes */
270: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
271: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
272: PetscCall(SetAndAssemble(A, B, dadd, dadd, 0.0, ADD_VALUES));
273: PetscCall(SetAndAssemble(A, B, dins, dins, 0.0, INSERT_VALUES));
274: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
275: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
276: PetscCall(SetAndAssemble(A, B, dadd, dadd, 5.0, INSERT_VALUES));
277: PetscCall(CheckEntry(A, B, dadd, dadd, "diagonal location that a zero added under MAT_IGNORE_ZERO_ENTRIES must have created"));
278: PetscCall(SetAndAssemble(A, B, dins, dins, 5.0, INSERT_VALUES));
279: PetscCall(CheckEntry(A, B, dins, dins, "diagonal location that a zero inserted under MAT_IGNORE_ZERO_ENTRIES must have created"));
281: /* the same off-diagonal probe, but against the block MatDisAssemble_MPISELL() builds. New
282: locations must be allowed for the zero to reach the rebuilt block at all, since a column that
283: garray does not know is what triggers the disassembly. */
284: if (size > 1) {
285: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
286: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
287: PetscCall(SetAndAssemble(A, B, i, rebuilt, 0.0, INSERT_VALUES));
288: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
289: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
290: PetscCall(SetAndAssemble(A, B, i, rebuilt, 5.0, INSERT_VALUES));
291: PetscCall(CheckEntry(A, B, i, rebuilt, "off-diagonal location that an ignored zero must not have created in the rebuilt block"));
292: }
294: PetscCall(MatDestroy(&A));
295: PetscCall(MatDestroy(&B));
297: if (size > 1) {
298: PetscCall(CheckDiagonalInOffDiagonalBlock(PETSC_TRUE));
299: PetscCall(CheckDiagonalInOffDiagonalBlock(PETSC_FALSE));
300: PetscCall(CheckOffDiagonalInDiagonalBlock());
301: }
302: PetscCall(PetscFinalize());
303: return 0;
304: }
306: /*TEST
308: testset:
309: output_file: output/empty.out
311: test:
312: suffix: 1
313: nsize: {{1 2}}
314: args: -mat_type sell
316: test:
317: suffix: cuda
318: nsize: {{1 2}}
319: requires: cuda !complex
320: args: -mat_type sellcuda
322: test:
323: suffix: hip
324: nsize: {{1 2}}
325: requires: hip !complex
326: args: -mat_type sellhip
328: TEST*/