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