Actual source code: ex225.c
1: static char help[] = "Test Hypre matrix APIs\n";
3: #include <petscmathypre.h>
5: int main(int argc, char **args)
6: {
7: Mat A, B, C;
8: PetscReal err;
9: PetscInt i, j, M = 20;
10: PetscMPIInt NP;
11: PetscInt *rows;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &args, NULL, help));
15: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &NP));
16: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
17: PetscCheck(M >= 6, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Matrix has to have more than 6 columns");
18: /* Hypre matrix */
19: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
20: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, M));
21: PetscCall(MatSetType(B, MATHYPRE));
22: PetscCall(MatHYPRESetPreallocation(B, 9, NULL, 9, NULL));
24: /* PETSc AIJ matrix */
25: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
26: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, M));
27: PetscCall(MatSetType(A, MATAIJ));
28: PetscCall(MatSeqAIJSetPreallocation(A, 9, NULL));
29: PetscCall(MatMPIAIJSetPreallocation(A, 9, NULL, 9, NULL));
31: /*Set Values */
32: for (i = 0; i < M; i++) {
33: PetscInt cols[] = {0, 1, 2, 3, 4, 5};
34: PetscScalar vals[6] = {0};
35: PetscScalar value[] = {100};
36: for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;
38: PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, ADD_VALUES));
39: PetscCall(MatSetValues(B, 1, &i, 1, &i, value, ADD_VALUES));
40: PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, ADD_VALUES));
41: PetscCall(MatSetValues(A, 1, &i, 1, &i, value, ADD_VALUES));
42: }
44: /* MAT_FLUSH_ASSEMBLY currently not supported */
45: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
46: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
47: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
50: /* Compare A and B */
51: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
52: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
53: PetscCall(MatNorm(C, NORM_INFINITY, &err));
54: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues %g", err);
55: PetscCall(MatDestroy(&C));
57: /* MatZeroRows */
58: PetscCall(PetscMalloc1(M, &rows));
59: for (i = 0; i < M; i++) rows[i] = i;
60: PetscCall(MatZeroRows(B, M, rows, 10.0, NULL, NULL));
61: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
62: PetscCall(MatZeroRows(A, M, rows, 10.0, NULL, NULL));
63: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
64: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
65: PetscCall(MatNorm(C, NORM_INFINITY, &err));
66: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatZeroRows %g", err);
67: PetscCall(MatDestroy(&C));
68: PetscCall(PetscFree(rows));
70: /* Test MatZeroEntries */
71: PetscCall(MatZeroEntries(B));
72: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
73: PetscCall(MatNorm(C, NORM_INFINITY, &err));
74: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatZeroEntries %g", err);
75: PetscCall(MatDestroy(&C));
77: /* Insert Values */
78: for (i = 0; i < M; i++) {
79: PetscInt cols[] = {0, 1, 2, 3, 4, 5};
80: PetscScalar vals[6] = {0};
81: PetscScalar value[] = {100};
83: for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;
85: PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, INSERT_VALUES));
86: PetscCall(MatSetValues(B, 1, &i, 1, &i, value, INSERT_VALUES));
87: PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, INSERT_VALUES));
88: PetscCall(MatSetValues(A, 1, &i, 1, &i, value, INSERT_VALUES));
89: }
91: /* MAT_FLUSH_ASSEMBLY currently not supported */
92: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
93: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
94: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
95: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
97: /* Rows are not sorted with HYPRE so we need an intermediate sort
98: They use a temporary buffer, so we can sort inplace the const memory */
99: {
100: const PetscInt *idxA, *idxB;
101: const PetscScalar *vA, *vB;
102: PetscInt rstart, rend, nzA, nzB;
103: PetscInt cols[] = {0, 1, 2, 3, 4, -5};
104: PetscInt *rows;
105: PetscScalar *valuesA, *valuesB;
106: PetscBool flg;
108: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
109: for (i = rstart; i < rend; i++) {
110: PetscCall(MatGetRow(A, i, &nzA, &idxA, &vA));
111: PetscCall(MatGetRow(B, i, &nzB, &idxB, &vB));
112: PetscCheck(nzA == nzB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT, nzA - nzB);
113: PetscCall(PetscSortIntWithScalarArray(nzB, (PetscInt *)idxB, (PetscScalar *)vB));
114: PetscCall(PetscArraycmp(idxA, idxB, nzA, &flg));
115: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (indices)", i);
116: PetscCall(PetscArraycmp(vA, vB, nzA, &flg));
117: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (values)", i);
118: PetscCall(MatRestoreRow(A, i, &nzA, &idxA, &vA));
119: PetscCall(MatRestoreRow(B, i, &nzB, &idxB, &vB));
120: }
122: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
123: PetscCall(PetscCalloc3((rend - rstart) * 6, &valuesA, (rend - rstart) * 6, &valuesB, rend - rstart, &rows));
124: for (i = rstart; i < rend; i++) rows[i - rstart] = i;
125: PetscCall(MatGetValues(A, rend - rstart, rows, 6, cols, valuesA));
126: PetscCall(MatGetValues(B, rend - rstart, rows, 6, cols, valuesB));
127: PetscCall(PetscArraycmp(valuesA, valuesB, 6 * (rend - rstart), &flg));
128: if (!flg) {
129: PetscCall(PetscScalarView(6 * (rend - rstart), valuesA, PETSC_VIEWER_STDOUT_WORLD));
130: PetscCall(PetscScalarView(6 * (rend - rstart), valuesB, PETSC_VIEWER_STDOUT_WORLD));
131: }
132: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error unexpected values from MatGetValues");
133: PetscCall(PetscFree3(valuesA, valuesB, rows));
134: }
136: /* Compare A and B */
137: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
138: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
139: PetscCall(MatNorm(C, NORM_INFINITY, &err));
140: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues with INSERT_VALUES %g", err);
142: PetscCall(MatDestroy(&A));
143: PetscCall(MatDestroy(&B));
144: PetscCall(MatDestroy(&C));
146: PetscCall(PetscFinalize());
147: return 0;
148: }
150: /*TEST
152: build:
153: requires: hypre
155: test:
156: suffix: 1
157: requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
158: output_file: output/empty.out
160: test:
161: suffix: 2
162: requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
163: output_file: output/empty.out
164: nsize: 2
166: TEST*/