1 static char help[] = "Test Hypre matrix APIs\n";
2
3 #include <petscmathypre.h>
4
main(int argc,char ** args)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;
13
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));
25
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));
32
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;
39
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 }
45
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));
51
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));
58
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));
71
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));
78
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};
84
85 for (j = 0; j < 6; j++) vals[j] = ((PetscReal)j) / NP;
86
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 }
92
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));
98
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;
109
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 }
123
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;
127
128 PetscCall(MatGetValues(A, rend - rstart, rows, 6, cols, valuesA));
129 PetscCall(MatGetValues(B, rend - rstart, rows, 6, cols, valuesB));
130
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 }
137
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);
143
144 PetscCall(MatDestroy(&A));
145 PetscCall(MatDestroy(&B));
146 PetscCall(MatDestroy(&C));
147
148 PetscCall(PetscFinalize());
149 return 0;
150 }
151
152 /*TEST
153
154 build:
155 requires: hypre
156
157 test:
158 suffix: 1
159 requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
160 output_file: output/empty.out
161
162 test:
163 suffix: 2
164 requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
165 output_file: output/empty.out
166 nsize: 2
167
168 TEST*/
169