xref: /petsc/src/mat/tests/ex16.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94)
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   PetscErrorCode    ierr;
11 
12   PetscFunctionBegin;
13   ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr);
14   ierr = MatDenseGetLDA(A,&lda);CHKERRQ(ierr);
15   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
16   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
17   for (i=rstart; i<rend; i++) {
18     for (j=0; j<N; j++) {
19       PetscInt ii = i - rstart, jj = j;
20       PetscReal v = (PetscReal)(one ? 1 : (1 + i + j*M));
21       PetscReal w = PetscRealPart(array[ii + jj*lda]);
22       if (PetscAbsReal(v-w) > 0) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix entry (%D,%D) should be %g, got %g",i,j,(double)v,(double)w);
23     }
24   }
25   ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr);
26   PetscFunctionReturn(0);
27 }
28 
29 #define CheckValuesIJ(A)  CheckValues(A,PETSC_FALSE)
30 #define CheckValuesOne(A) CheckValues(A,PETSC_TRUE)
31 
32 int main(int argc,char **args)
33 {
34   Mat            A;
35   PetscInt       i,j,M = 4,N = 3,rstart,rend;
36   PetscErrorCode ierr;
37   PetscScalar    *array;
38   PetscViewer    view;
39 
40   ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
41   /*
42       Create a parallel dense matrix shared by all processors
43   */
44   ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,N,NULL,&A);CHKERRQ(ierr);
45 
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       ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);
53     }
54   }
55   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
56   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
57 
58   /*
59       Store the binary matrix to a file
60   */
61   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view);CHKERRQ(ierr);
62   for (i=0; i<2; i++) {
63     ierr = MatView(A,view);CHKERRQ(ierr);
64     ierr = PetscViewerPushFormat(view,PETSC_VIEWER_NATIVE);CHKERRQ(ierr);
65     ierr = MatView(A,view);CHKERRQ(ierr);
66     ierr = PetscViewerPopFormat(view);CHKERRQ(ierr);
67   }
68   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
69   ierr = MatDestroy(&A);CHKERRQ(ierr);
70 
71   /*
72       Now reload the matrix and check its values
73   */
74   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_READ,&view);CHKERRQ(ierr);
75   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
76   ierr = MatSetType(A,MATMPIDENSE);CHKERRQ(ierr);
77   for (i=0; i<4; i++) {
78     if (i > 0) {ierr = MatZeroEntries(A);CHKERRQ(ierr);}
79     ierr = MatLoad(A,view);CHKERRQ(ierr);
80     ierr = CheckValuesIJ(A);CHKERRQ(ierr);
81   }
82   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
83 
84   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
85   ierr = PetscMalloc1((rend-rstart)*N,&array);CHKERRQ(ierr);
86   for (i=0; i<(rend-rstart)*N; i++) array[i] = (PetscReal)1;
87   ierr = MatDensePlaceArray(A,array);CHKERRQ(ierr);
88   ierr = CheckValuesOne(A);CHKERRQ(ierr);
89   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_WRITE,&view);CHKERRQ(ierr);
90   ierr = MatView(A,view);CHKERRQ(ierr);
91   ierr = MatDenseResetArray(A);CHKERRQ(ierr);
92   ierr = PetscFree(array);CHKERRQ(ierr);
93   ierr = CheckValuesIJ(A);CHKERRQ(ierr);
94   ierr = PetscViewerBinarySetSkipHeader(view,PETSC_TRUE);CHKERRQ(ierr);
95   ierr = MatView(A,view);CHKERRQ(ierr);
96   ierr = PetscViewerBinarySetSkipHeader(view,PETSC_FALSE);CHKERRQ(ierr);
97   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
98   ierr = MatDestroy(&A);CHKERRQ(ierr);
99 
100   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
101   ierr = MatSetType(A,MATMPIDENSE);CHKERRQ(ierr);
102   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_READ,&view);CHKERRQ(ierr);
103   ierr = MatLoad(A,view);CHKERRQ(ierr);
104   ierr = CheckValuesOne(A);CHKERRQ(ierr);
105   ierr = MatZeroEntries(A);CHKERRQ(ierr);
106   ierr = PetscViewerBinarySetSkipHeader(view,PETSC_TRUE);CHKERRQ(ierr);
107   ierr = MatLoad(A,view);CHKERRQ(ierr);
108   ierr = PetscViewerBinarySetSkipHeader(view,PETSC_FALSE);CHKERRQ(ierr);
109   ierr = CheckValuesIJ(A);CHKERRQ(ierr);
110   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
111   ierr = MatDestroy(&A);CHKERRQ(ierr);
112 
113   {
114     PetscInt m = PETSC_DECIDE, n = PETSC_DECIDE;
115     ierr = PetscSplitOwnership(PETSC_COMM_WORLD,&m,&M);CHKERRQ(ierr);
116     ierr = PetscSplitOwnership(PETSC_COMM_WORLD,&n,&N);CHKERRQ(ierr);
117     /* TODO: MatCreateDense requires data!=NULL at all processes! */
118     ierr = PetscMalloc1(m*N+1,&array);CHKERRQ(ierr);
119 
120     ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"matrix.dat",FILE_MODE_READ,&view);CHKERRQ(ierr);
121     ierr = MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,array,&A);CHKERRQ(ierr);
122     ierr = MatLoad(A,view);CHKERRQ(ierr);
123     ierr = CheckValuesOne(A);CHKERRQ(ierr);
124     ierr = PetscViewerBinarySetSkipHeader(view,PETSC_TRUE);CHKERRQ(ierr);
125     ierr = MatLoad(A,view);CHKERRQ(ierr);
126     ierr = PetscViewerBinarySetSkipHeader(view,PETSC_FALSE);CHKERRQ(ierr);
127     ierr = CheckValuesIJ(A);CHKERRQ(ierr);
128     ierr = MatDestroy(&A);CHKERRQ(ierr);
129     ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
130 
131     ierr = MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,array,&A);CHKERRQ(ierr);
132     ierr = CheckValuesIJ(A);CHKERRQ(ierr);
133     ierr = MatDestroy(&A);CHKERRQ(ierr);
134 
135     ierr = PetscFree(array);CHKERRQ(ierr);
136   }
137 
138   ierr = PetscFinalize();
139   return ierr;
140 }
141 
142 
143 /*TEST
144 
145    testset:
146       args: -viewer_binary_mpiio 0
147       output_file: output/ex16.out
148       test:
149         suffix: stdio_1
150         nsize: 1
151       test:
152         suffix: stdio_2
153         nsize: 2
154       test:
155         suffix: stdio_3
156         nsize: 3
157       test:
158         suffix: stdio_4
159         nsize: 4
160       test:
161         suffix: stdio_5
162         nsize: 5
163 
164    testset:
165       requires: mpiio
166       args: -viewer_binary_mpiio 1
167       output_file: output/ex16.out
168       test:
169         suffix: mpiio_1
170         nsize: 1
171       test:
172         suffix: mpiio_2
173         nsize: 2
174       test:
175         suffix: mpiio_3
176         nsize: 3
177       test:
178         suffix: mpiio_4
179         nsize: 4
180       test:
181         suffix: mpiio_5
182         nsize: 5
183 
184 
185 TEST*/
186