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