xref: /petsc/src/mat/tests/ex16.c (revision a69119a591a03a9d906b29c0a4e9802e4d7c9795)
1 static char help[] = "Tests MatDenseGetArray() and MatView()/MatLoad() with binary viewers.\n\n";
2 
3 #include <petscmat.h>
4 #include <petscviewer.h>
5 
6 static PetscErrorCode CheckValues(Mat A, PetscBool one) {
7   const PetscScalar *array;
8   PetscInt           M, N, rstart, rend, lda, i, j;
9 
10   PetscFunctionBegin;
11   PetscCall(MatDenseGetArrayRead(A, &array));
12   PetscCall(MatDenseGetLDA(A, &lda));
13   PetscCall(MatGetSize(A, &M, &N));
14   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
15   for (i = rstart; i < rend; i++) {
16     for (j = 0; j < N; j++) {
17       PetscInt  ii = i - rstart, jj = j;
18       PetscReal v = (PetscReal)(one ? 1 : (1 + i + j * M));
19       PetscReal w = PetscRealPart(array[ii + jj * lda]);
20       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);
21     }
22   }
23   PetscCall(MatDenseRestoreArrayRead(A, &array));
24   PetscFunctionReturn(0);
25 }
26 
27 #define CheckValuesIJ(A)  CheckValues(A, PETSC_FALSE)
28 #define CheckValuesOne(A) CheckValues(A, PETSC_TRUE)
29 
30 int main(int argc, char **args) {
31   Mat          A;
32   PetscInt     i, j, M = 4, N = 3, rstart, rend;
33   PetscScalar *array;
34   char         mattype[256];
35   PetscViewer  view;
36 
37   PetscFunctionBeginUser;
38   PetscCall(PetscInitialize(&argc, &args, NULL, help));
39   PetscCall(PetscStrcpy(mattype, MATMPIDENSE));
40   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_type", mattype, sizeof(mattype), NULL));
41   /*
42       Create a parallel dense matrix shared by all processors
43   */
44   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &A));
45   PetscCall(MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A));
46   /*
47      Set values into the matrix
48   */
49   for (i = 0; i < M; i++) {
50     for (j = 0; j < N; j++) {
51       PetscScalar v = (PetscReal)(1 + i + j * M);
52       PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
53     }
54   }
55   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
56   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
57   PetscCall(MatScale(A, 2.0));
58   PetscCall(MatScale(A, 1.0 / 2.0));
59 
60   /*
61       Store the binary matrix to a file
62   */
63   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
64   for (i = 0; i < 2; i++) {
65     PetscCall(MatView(A, view));
66     PetscCall(PetscViewerPushFormat(view, PETSC_VIEWER_NATIVE));
67     PetscCall(MatView(A, view));
68     PetscCall(PetscViewerPopFormat(view));
69   }
70   PetscCall(PetscViewerDestroy(&view));
71   PetscCall(MatDestroy(&A));
72 
73   /*
74       Now reload the matrix and check its values
75   */
76   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
77   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
78   PetscCall(MatSetType(A, mattype));
79   for (i = 0; i < 4; i++) {
80     if (i > 0) PetscCall(MatZeroEntries(A));
81     PetscCall(MatLoad(A, view));
82     PetscCall(CheckValuesIJ(A));
83   }
84   PetscCall(PetscViewerDestroy(&view));
85 
86   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
87   PetscCall(PetscMalloc1((rend - rstart) * N, &array));
88   for (i = 0; i < (rend - rstart) * N; i++) array[i] = (PetscReal)1;
89   PetscCall(MatDensePlaceArray(A, array));
90   PetscCall(MatScale(A, 2.0));
91   PetscCall(MatScale(A, 1.0 / 2.0));
92   PetscCall(CheckValuesOne(A));
93   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
94   PetscCall(MatView(A, view));
95   PetscCall(MatDenseResetArray(A));
96   PetscCall(PetscFree(array));
97   PetscCall(CheckValuesIJ(A));
98   PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
99   PetscCall(MatView(A, view));
100   PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
101   PetscCall(PetscViewerDestroy(&view));
102   PetscCall(MatDestroy(&A));
103 
104   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
105   PetscCall(MatSetType(A, mattype));
106   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
107   PetscCall(MatLoad(A, view));
108   PetscCall(CheckValuesOne(A));
109   PetscCall(MatZeroEntries(A));
110   PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
111   PetscCall(MatLoad(A, view));
112   PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
113   PetscCall(CheckValuesIJ(A));
114   PetscCall(PetscViewerDestroy(&view));
115   PetscCall(MatDestroy(&A));
116 
117   {
118     PetscInt m = PETSC_DECIDE, n = PETSC_DECIDE;
119     PetscCall(PetscSplitOwnership(PETSC_COMM_WORLD, &m, &M));
120     PetscCall(PetscSplitOwnership(PETSC_COMM_WORLD, &n, &N));
121     /* TODO: MatCreateDense requires data!=NULL at all processes! */
122     PetscCall(PetscMalloc1(m * N + 1, &array));
123 
124     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
125     PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, array, &A));
126     PetscCall(MatLoad(A, view));
127     PetscCall(CheckValuesOne(A));
128     PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
129     PetscCall(MatLoad(A, view));
130     PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
131     PetscCall(CheckValuesIJ(A));
132     PetscCall(MatDestroy(&A));
133     PetscCall(PetscViewerDestroy(&view));
134 
135     PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, array, &A));
136     PetscCall(CheckValuesIJ(A));
137     PetscCall(MatDestroy(&A));
138 
139     PetscCall(PetscFree(array));
140   }
141 
142   PetscCall(PetscFinalize());
143   return 0;
144 }
145 
146 /*TEST
147 
148    testset:
149       args: -viewer_binary_mpiio 0
150       output_file: output/ex16.out
151       test:
152         suffix: stdio_1
153         nsize: 1
154         args: -mat_type seqdense
155       test:
156         suffix: stdio_2
157         nsize: 2
158       test:
159         suffix: stdio_3
160         nsize: 3
161       test:
162         suffix: stdio_4
163         nsize: 4
164       test:
165         suffix: stdio_5
166         nsize: 5
167       test:
168         requires: cuda
169         args: -mat_type seqdensecuda
170         suffix: stdio_cuda_1
171         nsize: 1
172       test:
173         requires: cuda
174         args: -mat_type mpidensecuda
175         suffix: stdio_cuda_2
176         nsize: 2
177       test:
178         requires: cuda
179         args: -mat_type mpidensecuda
180         suffix: stdio_cuda_3
181         nsize: 3
182       test:
183         requires: cuda
184         args: -mat_type mpidensecuda
185         suffix: stdio_cuda_4
186         nsize: 4
187       test:
188         requires: cuda
189         args: -mat_type mpidensecuda
190         suffix: stdio_cuda_5
191         nsize: 5
192 
193    testset:
194       requires: mpiio
195       args: -viewer_binary_mpiio 1
196       output_file: output/ex16.out
197       test:
198         suffix: mpiio_1
199         nsize: 1
200       test:
201         suffix: mpiio_2
202         nsize: 2
203       test:
204         suffix: mpiio_3
205         nsize: 3
206       test:
207         suffix: mpiio_4
208         nsize: 4
209       test:
210         suffix: mpiio_5
211         nsize: 5
212       test:
213         requires: cuda
214         args: -mat_type mpidensecuda
215         suffix: mpiio_cuda_1
216         nsize: 1
217       test:
218         requires: cuda
219         args: -mat_type mpidensecuda
220         suffix: mpiio_cuda_2
221         nsize: 2
222       test:
223         requires: cuda
224         args: -mat_type mpidensecuda
225         suffix: mpiio_cuda_3
226         nsize: 3
227       test:
228         requires: cuda
229         args: -mat_type mpidensecuda
230         suffix: mpiio_cuda_4
231         nsize: 4
232       test:
233         requires: cuda
234         args: -mat_type mpidensecuda
235         suffix: mpiio_cuda_5
236         nsize: 5
237 
238 TEST*/
239