Actual source code: ex270.c
1: #include <petscsys.h>
2: #include <petscmat.h>
3: #include <petscoptions.h>
5: static const char help[] = "Test MatGetValue/Row for hypre matrix on device\n";
7: static PetscErrorCode CheckBandedMatrix(MPI_Comm comm, PetscInt N, PetscInt bw, const char fill[])
8: {
9: Mat A;
10: PetscInt rstart, rend;
11: PetscBool iscoo, isconvert, ok = PETSC_TRUE;
12: PetscScalar expected;
14: PetscFunctionBeginUser;
15: PetscCall(PetscStrcmp(fill, "coo", &iscoo));
16: PetscCall(PetscStrcmp(fill, "convert", &isconvert));
18: PetscCall(MatCreate(comm, &A));
19: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N));
20: PetscCall(MatSetType(A, isconvert ? MATAIJ : MATHYPRE));
21: PetscCall(MatSeqAIJSetPreallocation(A, 2 * bw + 1, NULL));
22: PetscCall(MatMPIAIJSetPreallocation(A, 2 * bw + 1, NULL, 2 * bw + 1, NULL));
23: PetscCall(MatSetUp(A));
24: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
26: if (iscoo) {
27: PetscCount ncoo = 0;
28: PetscInt *coo_i, *coo_j;
29: PetscInt nmax = (rend - rstart) * (2 * bw + 1); /* an upper bound; rows at the ends hold fewer */
30: PetscScalar *coo_v;
32: PetscCall(PetscMalloc3(nmax, &coo_i, nmax, &coo_j, nmax, &coo_v));
33: for (PetscInt i = rstart; i < rend; i++) {
34: for (PetscInt j = PetscMax(0, i - bw); j <= PetscMin(N - 1, i + bw); j++) {
35: coo_i[ncoo] = i;
36: coo_j[ncoo] = j;
37: coo_v[ncoo] = (PetscScalar)(100 * (i + 1) + j + 1);
38: ncoo++;
39: }
40: }
41: PetscCall(MatSetPreallocationCOO(A, ncoo, coo_i, coo_j));
42: PetscCall(MatSetValuesCOO(A, coo_v, INSERT_VALUES));
43: PetscCall(PetscFree3(coo_i, coo_j, coo_v));
44: } else {
45: for (PetscInt i = rstart; i < rend; i++) {
46: for (PetscInt j = PetscMax(0, i - bw); j <= PetscMin(N - 1, i + bw); j++) {
47: expected = (PetscScalar)(100 * (i + 1) + j + 1);
48: PetscCall(MatSetValues(A, 1, &i, 1, &j, &expected, INSERT_VALUES));
49: }
50: }
51: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
52: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
53: }
55: if (isconvert) {
56: Mat H;
58: PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &H));
59: PetscCall(MatDestroy(&A));
60: A = H;
61: /* the conversion leaves the matrix on the host; the reordering only matters once it moves */
62: PetscCall(MatBindToCPU(A, PETSC_FALSE));
63: }
65: for (PetscInt i = rstart; i < rend && ok; i++) {
66: const PetscInt *cols;
67: const PetscScalar *rowvals;
68: PetscInt ncols;
70: /* MatGetValues() over the whole row range, so that structural zeros are covered too */
71: for (PetscInt j = 0; j < N; j++) {
72: PetscScalar got;
74: expected = (j >= i - bw && j <= i + bw) ? (PetscScalar)(100 * (i + 1) + j + 1) : 0.0;
75: PetscCall(MatGetValues(A, 1, &i, 1, &j, &got));
76: if (PetscAbsScalar(got - expected) > PETSC_SMALL) ok = PETSC_FALSE;
77: }
79: PetscCall(MatGetRow(A, i, &ncols, &cols, &rowvals));
80: for (PetscInt k = 0; k < ncols; k++) {
81: expected = (PetscScalar)(100 * (i + 1) + cols[k] + 1);
82: if (PetscAbsScalar(rowvals[k] - expected) > PETSC_SMALL) ok = PETSC_FALSE;
83: }
84: PetscCall(MatRestoreRow(A, i, &ncols, &cols, &rowvals));
85: }
87: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ok, 1, MPI_C_BOOL, MPI_LAND, comm));
88: PetscCall(PetscPrintf(comm, "Banded matrix check (%s fill): %s\n", fill, ok ? "OK" : "FAILED"));
90: PetscCall(MatDestroy(&A));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: int main(int argc, char **argv)
95: {
96: PetscInt n = 10; // elements
97: PetscReal L = 10.0; // domain length
99: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
101: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
102: PetscCall(PetscOptionsGetReal(NULL, NULL, "-L", &L, NULL));
104: const PetscInt N = n + 1; // nodes
105: const PetscReal h = L / (PetscReal)n;
107: Mat M;
108: PetscCall(MatCreate(PETSC_COMM_WORLD, &M));
109: PetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, N, N));
110: PetscCall(MatSetType(M, MATAIJ));
111: PetscCall(MatSetFromOptions(M));
112: // Tridiagonal pattern: up to 3 nnz/row (ends have 2). Use same for on/off diag for simplicity.
113: PetscCall(MatSetUp(M));
114: PetscCall(MatMPIAIJSetPreallocation(M, 3, NULL, 3, NULL));
115: PetscCall(MatSeqAIJSetPreallocation(M, 3, NULL));
117: // Ownership range for rows (nodes)
118: PetscInt rstart, rend;
119: PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
121: // Element matrix (scaled)
122: const PetscScalar s = (PetscScalar)(h / 6.0);
123: const PetscScalar ae[4] = {2 * s, 1 * s, 1 * s, 2 * s};
125: // Assemble by looping over elements that start at locally-owned row i
126: PetscInt idx[2];
127: for (PetscInt i = PetscMax(rstart, 0); i < PetscMin(rend, n); ++i) {
128: idx[0] = i;
129: idx[1] = i + 1;
130: PetscCall(MatSetValues(M, 2, idx, 2, idx, ae, ADD_VALUES));
131: }
133: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
134: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
136: // --- Verification of MatGetRow ------------------------------------------
137: {
138: const PetscReal tol = 100 * PETSC_MACHINE_EPSILON; // tight but safe
139: PetscBool ok = PETSC_TRUE;
141: for (PetscInt i = rstart; i < rend; ++i) {
142: const PetscInt *cols;
143: const PetscScalar *getrowvals;
144: PetscScalar *getvals;
145: PetscInt ncols = 0, expected_ncols = 0;
146: PetscInt expected_cols[3];
147: PetscScalar expected_vals[3];
149: // Build expected data
150: if (i > 0) {
151: expected_cols[expected_ncols] = i - 1;
152: expected_vals[expected_ncols] = s;
153: ++expected_ncols;
154: }
155: expected_cols[expected_ncols] = i;
156: expected_vals[expected_ncols] = (i == 0 || i == N - 1) ? (2 * s) : (4 * s); // diag: h/3 at ends, 2h/3 interior
157: ++expected_ncols;
158: if (i + 1 < N) {
159: expected_cols[expected_ncols] = i + 1;
160: expected_vals[expected_ncols] = s;
161: ++expected_ncols;
162: }
164: PetscCall(PetscMalloc1(expected_ncols, &getvals));
165: PetscCall(MatGetRow(M, i, &ncols, &cols, &getrowvals));
166: PetscCall(MatGetValues(M, 1, &i, expected_ncols, (PetscInt *)expected_cols, getvals));
168: // Compare counts
169: if (ncols != expected_ncols) {
170: ok = PETSC_FALSE;
171: goto rowdone;
172: }
174: // Compare values (match by column)
175: for (PetscInt k = 0; k < ncols; ++k) {
176: PetscInt expected_k = -1;
177: /* Matrix is small. Just do a linear search */
178: for (PetscInt l = 0; l < expected_ncols; ++l) {
179: if (expected_cols[l] == cols[k]) {
180: expected_k = l;
181: break;
182: }
183: }
184: if (expected_k < 0) {
185: ok = PETSC_FALSE;
186: goto rowdone;
187: }
188: if ((PetscAbsScalar(getrowvals[k] - expected_vals[expected_k]) > tol) || (PetscAbsScalar(getvals[expected_k] - expected_vals[expected_k]) > tol)) {
189: ok = PETSC_FALSE;
190: goto rowdone;
191: }
192: }
194: rowdone:
195: PetscCall(MatRestoreRow(M, i, &ncols, &cols, &getrowvals));
196: PetscCall(PetscFree(getvals));
197: if (!ok) break;
198: }
200: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ok, 1, MPI_C_BOOL, MPI_LAND, PETSC_COMM_WORLD));
201: if (ok) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mass matrix check: OK\n"));
202: else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mass matrix check: FAILED\n"));
203: }
204: // --------------------------------------------------------------------------
206: PetscCall(MatDestroy(&M));
208: {
209: PetscInt bw = 3;
210: PetscBool flg;
211: char fill[16];
213: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bw", &bw, NULL));
214: PetscCall(PetscOptionsGetString(NULL, NULL, "-fill", fill, sizeof(fill), &flg));
215: if (!flg) PetscCall(PetscStrncpy(fill, "setvalues", sizeof(fill)));
216: PetscCall(CheckBandedMatrix(PETSC_COMM_WORLD, N, bw, fill));
217: }
219: PetscCall(PetscFinalize());
220: return 0;
221: }
223: /*TEST
225: test:
226: requires: hypre
227: suffix: 1
228: args: -mat_type hypre
230: test:
231: requires: hypre
232: suffix: 2
233: args: -mat_type hypre
234: nsize: 4
236: test:
237: requires: hypre
238: suffix: coo
239: output_file: output/ex270_coo.out
240: args: -mat_type hypre -fill coo
241: nsize: {{1 4}}
243: test:
244: requires: hypre
245: suffix: convert
246: output_file: output/ex270_convert.out
247: args: -mat_type hypre -fill convert
248: nsize: {{1 4}}
250: TEST*/