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: MPI_Comm comm;
12: PetscInt *rows;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &args, NULL, help));
16: comm = PETSC_COMM_WORLD;
17: PetscCallMPI(MPI_Comm_size(comm, &NP));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
19: PetscCheck(M >= 6, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Matrix has to have more than 6 columns");
20: /* Hypre matrix */
21: PetscCall(MatCreate(comm, &B));
22: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, M));
23: PetscCall(MatSetType(B, MATHYPRE));
24: PetscCall(MatHYPRESetPreallocation(B, 9, NULL, 9, NULL));
26: /* PETSc AIJ matrix */
27: PetscCall(MatCreate(comm, &A));
28: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, M));
29: PetscCall(MatSetType(A, MATAIJ));
30: PetscCall(MatSeqAIJSetPreallocation(A, 9, NULL));
31: PetscCall(MatMPIAIJSetPreallocation(A, 9, NULL, 9, NULL));
33: /*Set Values */
34: for (i = 0; i < M; i++) {
35: PetscInt cols[] = {0, 1, 2, 3, 4, 5};
36: PetscScalar vals[6] = {0};
37: PetscScalar value[] = {100};
38: for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;
40: PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, ADD_VALUES));
41: PetscCall(MatSetValues(B, 1, &i, 1, &i, value, ADD_VALUES));
42: PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, ADD_VALUES));
43: PetscCall(MatSetValues(A, 1, &i, 1, &i, value, ADD_VALUES));
44: }
46: /* MAT_FLUSH_ASSEMBLY currently not supported */
47: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
49: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
50: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
52: /* Compare A and B */
53: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
54: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
55: PetscCall(MatNorm(C, NORM_INFINITY, &err));
56: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues %g", err);
57: PetscCall(MatDestroy(&C));
59: /* MatZeroRows */
60: PetscCall(PetscMalloc1(M, &rows));
61: for (i = 0; i < M; i++) rows[i] = i;
62: PetscCall(MatZeroRows(B, M, rows, 10.0, NULL, NULL));
63: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
64: PetscCall(MatZeroRows(A, M, rows, 10.0, NULL, NULL));
65: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
66: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
67: PetscCall(MatNorm(C, NORM_INFINITY, &err));
68: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatZeroRows %g", err);
69: PetscCall(MatDestroy(&C));
70: PetscCall(PetscFree(rows));
72: /* Test MatZeroEntries */
73: PetscCall(MatZeroEntries(B));
74: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
75: PetscCall(MatNorm(C, NORM_INFINITY, &err));
76: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error MatZeroEntries %g", err);
77: PetscCall(MatDestroy(&C));
79: /* Insert Values */
80: for (i = 0; i < M; i++) {
81: PetscInt cols[] = {0, 1, 2, 3, 4, 5};
82: PetscScalar vals[6] = {0};
83: PetscScalar value[] = {100};
85: for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;
87: PetscCall(MatSetValues(B, 1, &i, 6, cols, vals, INSERT_VALUES));
88: PetscCall(MatSetValues(B, 1, &i, 1, &i, value, INSERT_VALUES));
89: PetscCall(MatSetValues(A, 1, &i, 6, cols, vals, INSERT_VALUES));
90: PetscCall(MatSetValues(A, 1, &i, 1, &i, value, INSERT_VALUES));
91: }
93: /* MAT_FLUSH_ASSEMBLY currently not supported */
94: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
95: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
96: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
97: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
99: /* Rows are not sorted with HYPRE so we need an intermediate sort
100: They use a temporary buffer, so we can sort inplace the const memory */
101: {
102: const PetscInt *idxA, *idxB;
103: const PetscScalar *vA, *vB;
104: PetscInt rstart, rend, nzA, nzB;
105: PetscInt cols[] = {0, 1, 2, 3, 4, -5};
106: PetscInt *rows;
107: PetscScalar *valuesA, *valuesB;
108: PetscBool flg;
110: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
111: for (i = rstart; i < rend; i++) {
112: PetscCall(MatGetRow(A, i, &nzA, &idxA, &vA));
113: PetscCall(MatGetRow(B, i, &nzB, &idxB, &vB));
114: PetscCheck(nzA == nzB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT, nzA - nzB);
115: PetscCall(PetscSortIntWithScalarArray(nzB, (PetscInt *)idxB, (PetscScalar *)vB));
116: PetscCall(PetscArraycmp(idxA, idxB, nzA, &flg));
117: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (indices)", i);
118: PetscCall(PetscArraycmp(vA, vB, nzA, &flg));
119: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetRow %" PetscInt_FMT " (values)", i);
120: PetscCall(MatRestoreRow(A, i, &nzA, &idxA, &vA));
121: PetscCall(MatRestoreRow(B, i, &nzB, &idxB, &vB));
122: }
124: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
125: PetscCall(PetscCalloc3((rend - rstart) * 6, &valuesA, (rend - rstart) * 6, &valuesB, rend - rstart, &rows));
126: for (i = rstart; i < rend; i++) rows[i - rstart] = i;
128: PetscCall(MatGetValues(A, rend - rstart, rows, 6, cols, valuesA));
129: PetscCall(MatGetValues(B, rend - rstart, rows, 6, cols, valuesB));
131: for (i = 0; i < (rend - rstart); i++) {
132: PetscCall(PetscArraycmp(valuesA + 6 * i, valuesB + 6 * i, 6, &flg));
133: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error MatGetValues %" PetscInt_FMT, i + rstart);
134: }
135: PetscCall(PetscFree3(valuesA, valuesB, rows));
136: }
138: /* Compare A and B */
139: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C));
140: PetscCall(MatAXPY(C, -1., A, SAME_NONZERO_PATTERN));
141: PetscCall(MatNorm(C, NORM_INFINITY, &err));
142: PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Error MatSetValues with INSERT_VALUES %g", err);
144: PetscCall(MatDestroy(&A));
145: PetscCall(MatDestroy(&B));
146: PetscCall(MatDestroy(&C));
148: PetscCall(PetscFinalize());
149: return 0;
150: }
152: /*TEST
154: build:
155: requires: hypre
157: test:
158: suffix: 1
159: requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
161: test:
162: suffix: 2
163: requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
164: output_file: output/ex225_1.out
165: nsize: 2
167: TEST*/