1 static char help[] = "Tests MatView()/MatLoad() with binary viewers for SBAIJ matrices.\n\n";
2
3 #include <petscmat.h>
4 #include <petscviewer.h>
5
6 #include <petsc/private/hashtable.h>
MakeValue(PetscInt i,PetscInt j,PetscInt M)7 static PetscReal MakeValue(PetscInt i, PetscInt j, PetscInt M)
8 {
9 PetscHash_t h = PetscHashCombine(PetscHashInt(i), PetscHashInt(j));
10 return (PetscReal)((h % 5 == 0) ? (1 + i + j * M) : 0);
11 }
12
CheckValuesAIJ(Mat A)13 static PetscErrorCode CheckValuesAIJ(Mat A)
14 {
15 PetscInt M, N, rstart, rend, i, j;
16 PetscReal v, w;
17 PetscScalar val;
18 PetscBool seqsbaij, mpisbaij, sbaij;
19
20 PetscFunctionBegin;
21 PetscCall(MatGetSize(A, &M, &N));
22 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
23 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &seqsbaij));
24 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &mpisbaij));
25 sbaij = (seqsbaij || mpisbaij) ? PETSC_TRUE : PETSC_FALSE;
26 for (i = rstart; i < rend; i++) {
27 for (j = (sbaij ? i : 0); j < N; j++) {
28 PetscCall(MatGetValue(A, i, j, &val));
29 v = MakeValue(i, j, M);
30 w = PetscRealPart(val);
31 PetscCheck(PetscAbsReal(v - w) <= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Matrix entry (%" PetscInt_FMT ",%" PetscInt_FMT ") should be %g, got %g", i, j, (double)v, (double)w);
32 }
33 }
34 PetscFunctionReturn(PETSC_SUCCESS);
35 }
36
main(int argc,char ** args)37 int main(int argc, char **args)
38 {
39 Mat A;
40 PetscInt M = 24, N = 24, bs = 3;
41 PetscInt rstart, rend, i, j;
42 PetscViewer view;
43
44 PetscFunctionBeginUser;
45 PetscCall(PetscInitialize(&argc, &args, NULL, help));
46 /*
47 Create a parallel SBAIJ matrix shared by all processors
48 */
49 PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A));
50
51 /*
52 Set values into the matrix
53 */
54 PetscCall(MatGetSize(A, &M, &N));
55 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
56 for (i = rstart; i < rend; i++) {
57 for (j = 0; j < N; j++) {
58 PetscReal v = MakeValue(i, j, M);
59 if (PetscAbsReal(v) > 0) PetscCall(MatSetValue(A, i, j, v, INSERT_VALUES));
60 }
61 }
62 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
63 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
64 PetscCall(MatViewFromOptions(A, NULL, "-mat_base_view"));
65
66 /*
67 Store the binary matrix to a file
68 */
69 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
70 for (i = 0; i < 3; i++) PetscCall(MatView(A, view));
71 PetscCall(PetscViewerDestroy(&view));
72 PetscCall(MatDestroy(&A));
73
74 /*
75 Now reload the matrix and check its values
76 */
77 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
78 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
79 PetscCall(MatSetType(A, MATSBAIJ));
80 for (i = 0; i < 3; i++) {
81 if (i > 0) PetscCall(MatZeroEntries(A));
82 PetscCall(MatLoad(A, view));
83 PetscCall(CheckValuesAIJ(A));
84 }
85 PetscCall(PetscViewerDestroy(&view));
86 PetscCall(MatViewFromOptions(A, NULL, "-mat_load_view"));
87 PetscCall(MatDestroy(&A));
88
89 /*
90 Reload in SEQSBAIJ matrix and check its values
91 */
92 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_READ, &view));
93 PetscCall(MatCreate(PETSC_COMM_SELF, &A));
94 PetscCall(MatSetType(A, MATSEQSBAIJ));
95 for (i = 0; i < 3; i++) {
96 if (i > 0) PetscCall(MatZeroEntries(A));
97 PetscCall(MatLoad(A, view));
98 PetscCall(CheckValuesAIJ(A));
99 }
100 PetscCall(PetscViewerDestroy(&view));
101 PetscCall(MatDestroy(&A));
102
103 /*
104 Reload in MPISBAIJ matrix and check its values
105 */
106 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
107 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
108 PetscCall(MatSetType(A, MATMPISBAIJ));
109 for (i = 0; i < 3; i++) {
110 if (i > 0) PetscCall(MatZeroEntries(A));
111 PetscCall(MatLoad(A, view));
112 PetscCall(CheckValuesAIJ(A));
113 }
114 PetscCall(PetscViewerDestroy(&view));
115 PetscCall(MatDestroy(&A));
116
117 PetscCall(PetscFinalize());
118 return 0;
119 }
120
121 /*TEST
122
123 testset:
124 args: -viewer_binary_mpiio 0
125 output_file: output/empty.out
126 test:
127 suffix: stdio_1
128 nsize: 1
129 test:
130 suffix: stdio_2
131 nsize: 2
132 test:
133 suffix: stdio_3
134 nsize: 3
135 test:
136 suffix: stdio_4
137 nsize: 4
138 test:
139 suffix: stdio_5
140 nsize: 4
141
142 testset:
143 requires: mpiio
144 args: -viewer_binary_mpiio 1
145 output_file: output/empty.out
146 test:
147 suffix: mpiio_1
148 nsize: 1
149 test:
150 suffix: mpiio_2
151 nsize: 2
152 test:
153 suffix: mpiio_3
154 nsize: 3
155 test:
156 suffix: mpiio_4
157 nsize: 4
158 test:
159 suffix: mpiio_5
160 nsize: 5
161
162 TEST*/
163