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