1
2 static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
3
4 #include <petscmat.h>
5
main(int argc,char ** args)6 int main(int argc,char **args)
7 {
8 Vec x,y,u,s1,s2;
9 Mat A,sA,sB;
10 PetscRandom rctx;
11 PetscReal r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
12 PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
13 PetscInt n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=1;
14 PetscErrorCode ierr;
15 PetscMPIInt size,rank;
16 PetscBool flg;
17
18 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
19 ierr = PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr);
20 ierr = PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);CHKERRQ(ierr);
21 ierr = PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL);CHKERRQ(ierr);
22
23 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
24 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
25
26 /* Create a BAIJ matrix A */
27 n = mbs*bs;
28 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
29 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
30 ierr = MatSetType(A,MATBAIJ);CHKERRQ(ierr);
31 ierr = MatSetFromOptions(A);CHKERRQ(ierr);
32 ierr = MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL);CHKERRQ(ierr);
33 ierr = MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL);CHKERRQ(ierr);
34 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
35
36 if (bs == 1) {
37 if (prob == 1) { /* tridiagonal matrix */
38 value[0] = -1.0; value[1] = 0.0; value[2] = -1.0;
39 for (i=1; i<n-1; i++) {
40 col[0] = i-1; col[1] = i; col[2] = i+1;
41 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
42 }
43 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
44 value[0]= 0.1; value[1]=-1.0; value[2]=0.0;
45 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
46
47 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
48 value[0] = 0.0; value[1] = -1.0; value[2]=0.1;
49 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
50 } else if (prob ==2) { /* matrix for the five point stencil */
51 n1 = (int) PetscSqrtReal((PetscReal)n);
52 for (i=0; i<n1; i++) {
53 for (j=0; j<n1; j++) {
54 Ii = j + n1*i;
55 if (i>0) {J = Ii - n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
56 if (i<n1-1) {J = Ii + n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
57 if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
58 if (j<n1-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
59 ierr = MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr);
60 }
61 }
62 }
63 /* end of if (bs == 1) */
64 } else { /* bs > 1 */
65 for (block=0; block<n/bs; block++) {
66 /* diagonal blocks */
67 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
68 for (i=1+block*bs; i<bs-1+block*bs; i++) {
69 col[0] = i-1; col[1] = i; col[2] = i+1;
70 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
71 }
72 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
73 value[0]=-1.0; value[1]=4.0;
74 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
75
76 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
77 value[0]=4.0; value[1] = -1.0;
78 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
79 }
80 /* off-diagonal blocks */
81 value[0]=-1.0;
82 for (i=0; i<(n/bs-1)*bs; i++) {
83 col[0]=i+bs;
84 ierr = MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
85 col[0]=i; row=i+bs;
86 ierr = MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
87 }
88 }
89 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
90 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
91 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
92
93 /* Get SBAIJ matrix sA from A */
94 ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
95
96 /* Test MatGetSize(), MatGetLocalSize() */
97 ierr = MatGetSize(sA, &i,&j);CHKERRQ(ierr);
98 ierr = MatGetSize(A, &i2,&j2);CHKERRQ(ierr);
99 i -= i2; j -= j2;
100 if (i || j) {
101 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);CHKERRQ(ierr);
102 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
103 }
104
105 ierr = MatGetLocalSize(sA, &i,&j);CHKERRQ(ierr);
106 ierr = MatGetLocalSize(A, &i2,&j2);CHKERRQ(ierr);
107 i2 -= i; j2 -= j;
108 if (i2 || j2) {
109 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);CHKERRQ(ierr);
110 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
111 }
112
113 /* vectors */
114 /*--------------------*/
115 /* i is obtained from MatGetLocalSize() */
116 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
117 ierr = VecSetSizes(x,i,PETSC_DECIDE);CHKERRQ(ierr);
118 ierr = VecSetFromOptions(x);CHKERRQ(ierr);
119 ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
120 ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
121 ierr = VecDuplicate(x,&s1);CHKERRQ(ierr);
122 ierr = VecDuplicate(x,&s2);CHKERRQ(ierr);
123
124 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
125 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
126 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
127 ierr = VecSet(u,one);CHKERRQ(ierr);
128
129 /* Test MatNorm() */
130 ierr = MatNorm(A,NORM_FROBENIUS,&r1);CHKERRQ(ierr);
131 ierr = MatNorm(sA,NORM_FROBENIUS,&r2);CHKERRQ(ierr);
132 rnorm = PetscAbsReal(r1-r2)/r2;
133 if (rnorm > tol && !rank) {
134 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
135 }
136 ierr = MatNorm(A,NORM_INFINITY,&r1);CHKERRQ(ierr);
137 ierr = MatNorm(sA,NORM_INFINITY,&r2);CHKERRQ(ierr);
138 rnorm = PetscAbsReal(r1-r2)/r2;
139 if (rnorm > tol && !rank) {
140 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
141 }
142 ierr = MatNorm(A,NORM_1,&r1);CHKERRQ(ierr);
143 ierr = MatNorm(sA,NORM_1,&r2);CHKERRQ(ierr);
144 rnorm = PetscAbsReal(r1-r2)/r2;
145 if (rnorm > tol && !rank) {
146 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
147 }
148
149 /* Test MatGetOwnershipRange() */
150 ierr = MatGetOwnershipRange(sA,&rstart,&rend);CHKERRQ(ierr);
151 ierr = MatGetOwnershipRange(A,&i2,&j2);CHKERRQ(ierr);
152 i2 -= rstart; j2 -= rend;
153 if (i2 || j2) {
154 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);CHKERRQ(ierr);
155 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
156 }
157
158 /* Test MatDiagonalScale() */
159 ierr = MatDiagonalScale(A,x,x);CHKERRQ(ierr);
160 ierr = MatDiagonalScale(sA,x,x);CHKERRQ(ierr);
161 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
162 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
163
164 /* Test MatGetDiagonal(), MatScale() */
165 ierr = MatGetDiagonal(A,s1);CHKERRQ(ierr);
166 ierr = MatGetDiagonal(sA,s2);CHKERRQ(ierr);
167 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
168 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
169 r1 -= r2;
170 if (r1<-tol || r1>tol) {
171 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);CHKERRQ(ierr);
172 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
173 }
174
175 ierr = MatScale(A,alpha);CHKERRQ(ierr);
176 ierr = MatScale(sA,alpha);CHKERRQ(ierr);
177
178 /* Test MatGetRowMaxAbs() */
179 ierr = MatGetRowMaxAbs(A,s1,NULL);CHKERRQ(ierr);
180 ierr = MatGetRowMaxAbs(sA,s2,NULL);CHKERRQ(ierr);
181
182 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
183 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
184 r1 -= r2;
185 if (r1<-tol || r1>tol) {
186 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");CHKERRQ(ierr);
187 }
188
189 /* Test MatMult(), MatMultAdd() */
190 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
191 if (!flg) {
192 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);CHKERRQ(ierr);
193 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
194 }
195
196 ierr = MatMultAddEqual(A,sA,10,&flg);CHKERRQ(ierr);
197 if (!flg) {
198 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);CHKERRQ(ierr);
199 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
200 }
201
202 /* Test MatMultTranspose(), MatMultTransposeAdd() */
203 for (i=0; i<10; i++) {
204 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
205 ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr);
206 ierr = MatMultTranspose(sA,x,s2);CHKERRQ(ierr);
207 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
208 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
209 r1 -= r2;
210 if (r1<-tol || r1>tol) {
211 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);CHKERRQ(ierr);
212 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
213 }
214 }
215 for (i=0; i<10; i++) {
216 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
217 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr);
218 ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr);
219 ierr = MatMultTransposeAdd(sA,x,y,s2);CHKERRQ(ierr);
220 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
221 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
222 r1 -= r2;
223 if (r1<-tol || r1>tol) {
224 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);CHKERRQ(ierr);
225 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
226 }
227 }
228
229 /* Test MatDuplicate() */
230 ierr = MatDuplicate(sA,MAT_COPY_VALUES,&sB);CHKERRQ(ierr);
231 ierr = MatEqual(sA,sB,&flg);CHKERRQ(ierr);
232 if (!flg) {
233 ierr = PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");CHKERRQ(ierr);
234 }
235 ierr = MatMultEqual(sA,sB,5,&flg);CHKERRQ(ierr);
236 if (!flg) {
237 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);CHKERRQ(ierr);
238 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
239 }
240 ierr = MatMultAddEqual(sA,sB,5,&flg);CHKERRQ(ierr);
241 if (!flg) {
242 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);CHKERRQ(ierr);
243 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
244 }
245 ierr = MatDestroy(&sB);CHKERRQ(ierr);
246 ierr = VecDestroy(&u);CHKERRQ(ierr);
247 ierr = VecDestroy(&x);CHKERRQ(ierr);
248 ierr = VecDestroy(&y);CHKERRQ(ierr);
249 ierr = VecDestroy(&s1);CHKERRQ(ierr);
250 ierr = VecDestroy(&s2);CHKERRQ(ierr);
251 ierr = MatDestroy(&sA);CHKERRQ(ierr);
252 ierr = MatDestroy(&A);CHKERRQ(ierr);
253 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
254 ierr = PetscFinalize();
255 return ierr;
256 }
257
258 /*TEST
259
260 test:
261 nsize: {{1 3}}
262 args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
263
264 TEST*/
265