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