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