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>
MakeValue(PetscInt i,PetscInt j,PetscInt M)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
CheckValuesAIJ(Mat A)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 PetscErrorCode ierr;
20
21 PetscFunctionBegin;
22 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
23 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
24 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&seqsbaij);CHKERRQ(ierr);
25 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&mpisbaij);CHKERRQ(ierr);
26 sbaij = (seqsbaij||mpisbaij) ? PETSC_TRUE : PETSC_FALSE;
27 for (i=rstart; i<rend; i++) {
28 for (j=(sbaij?i:0); j<N; j++) {
29 ierr = MatGetValue(A,i,j,&val);CHKERRQ(ierr);
30 v = MakeValue(i,j,M); w = PetscRealPart(val);
31 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);
32 }
33 }
34 PetscFunctionReturn(0);
35 }
36
main(int argc,char ** args)37 int main(int argc,char **args)
38 {
39 Mat A;
40 PetscInt M = 24,N = 24,bs = 3;
41 PetscInt rstart,rend,i,j;
42 PetscErrorCode ierr;
43 PetscViewer view;
44
45 ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
46 /*
47 Create a parallel SBAIJ matrix shared by all processors
48 */
49 ierr = MatCreateSBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,M,N,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A);CHKERRQ(ierr);
50
51 /*
52 Set values into the matrix
53 */
54 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
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,MATSBAIJ);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 SEQSBAIJ 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,MATSEQSBAIJ);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 MPISBAIJ 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,MATMPISBAIJ);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 /*TEST
126
127 testset:
128 args: -viewer_binary_mpiio 0
129 output_file: output/ex50.out
130 test:
131 suffix: stdio_1
132 nsize: 1
133 test:
134 suffix: stdio_2
135 nsize: 2
136 test:
137 suffix: stdio_3
138 nsize: 3
139 test:
140 suffix: stdio_4
141 nsize: 4
142 test:
143 suffix: stdio_5
144 nsize: 4
145
146 testset:
147 requires: mpiio
148 args: -viewer_binary_mpiio 1
149 output_file: output/ex50.out
150 test:
151 suffix: mpiio_1
152 nsize: 1
153 test:
154 suffix: mpiio_2
155 nsize: 2
156 test:
157 suffix: mpiio_3
158 nsize: 3
159 test:
160 suffix: mpiio_4
161 nsize: 4
162 test:
163 suffix: mpiio_5
164 nsize: 5
165
166 TEST*/
167