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