Actual source code: ex75.c
petsc-3.7.3 2016-08-01
2: static char help[] = "Tests the vatious routines in MatMPISBAIJ format.\n";
4: #include <petscmat.h>
8: int main(int argc,char **args)
9: {
10: Vec x,y,u,s1,s2;
11: Mat A,sA,sB;
12: PetscRandom rctx;
13: PetscReal r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
14: PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
15: 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=2;
17: PetscMPIInt size,rank;
18: PetscBool flg;
19: MatType type;
21: PetscInitialize(&argc,&args,(char*)0,help);
22: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
23: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: MPI_Comm_size(PETSC_COMM_WORLD,&size);
28: n = mbs*bs;
30: /* Assemble MPISBAIJ matrix sA */
31: MatCreate(PETSC_COMM_WORLD,&sA);
32: MatSetSizes(sA,PETSC_DECIDE,PETSC_DECIDE,n,n);
33: MatSetType(sA,MATSBAIJ);
34: MatSetFromOptions(sA);
35: MatGetType(sA,&type);
36: /* printf(" mattype: %s\n",type); */
37: MatMPISBAIJSetPreallocation(sA,bs,d_nz,NULL,o_nz,NULL);
38: MatSeqSBAIJSetPreallocation(sA,bs,d_nz,NULL);
39: MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
41: if (bs == 1) {
42: if (prob == 1) { /* tridiagonal matrix */
43: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
44: for (i=1; i<n-1; i++) {
45: col[0] = i-1; col[1] = i; col[2] = i+1;
46: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
47: }
48: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
49: value[0]= 0.1; value[1]=-1; value[2]=2;
50: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
52: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
53: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
54: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
55: } else if (prob ==2) { /* matrix for the five point stencil */
56: n1 = (int) PetscSqrtReal((PetscReal)n);
57: if (n1*n1 != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
59: for (i=0; i<n1; i++) {
60: for (j=0; j<n1; j++) {
61: Ii = j + n1*i;
62: if (i>0) {J = Ii - n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
63: if (i<n1-1) {J = Ii + n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
64: if (j>0) {J = Ii - 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
65: if (j<n1-1) {J = Ii + 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
66: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
67: }
68: }
69: }
70: /* end of if (bs == 1) */
71: } else { /* bs > 1 */
72: for (block=0; block<n/bs; block++) {
73: /* diagonal blocks */
74: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
75: for (i=1+block*bs; i<bs-1+block*bs; i++) {
76: col[0] = i-1; col[1] = i; col[2] = i+1;
77: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
78: }
79: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
80: value[0]=-1.0; value[1]=4.0;
81: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
83: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
84: value[0]=4.0; value[1] = -1.0;
85: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
86: }
87: /* off-diagonal blocks */
88: value[0]=-1.0;
89: for (i=0; i<(n/bs-1)*bs; i++) {
90: col[0]=i+bs;
91: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
92: col[0]=i; row=i+bs;
93: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
94: }
95: }
96: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
99: /* Test MatView() */
100: /*
101: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
102: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
103: */
104: /* Assemble MPIBAIJ matrix A */
105: MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,NULL,o_nz,NULL,&A);
106: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
108: if (bs == 1) {
109: if (prob == 1) { /* tridiagonal matrix */
110: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
111: for (i=1; i<n-1; i++) {
112: col[0] = i-1; col[1] = i; col[2] = i+1;
113: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
114: }
115: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
116: value[0]= 0.1; value[1]=-1; value[2]=2;
117: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
119: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
120: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
121: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
122: } else if (prob ==2) { /* matrix for the five point stencil */
123: n1 = (int) PetscSqrtReal((PetscReal)n);
124: for (i=0; i<n1; i++) {
125: for (j=0; j<n1; j++) {
126: Ii = j + n1*i;
127: if (i>0) {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
128: if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
129: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
130: if (j<n1-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
131: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
132: }
133: }
134: }
135: /* end of if (bs == 1) */
136: } else { /* bs > 1 */
137: for (block=0; block<n/bs; block++) {
138: /* diagonal blocks */
139: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
140: for (i=1+block*bs; i<bs-1+block*bs; i++) {
141: col[0] = i-1; col[1] = i; col[2] = i+1;
142: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
143: }
144: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
145: value[0]=-1.0; value[1]=4.0;
146: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
148: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
149: value[0]=4.0; value[1] = -1.0;
150: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
151: }
152: /* off-diagonal blocks */
153: value[0]=-1.0;
154: for (i=0; i<(n/bs-1)*bs; i++) {
155: col[0]=i+bs;
156: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
157: col[0]=i; row=i+bs;
158: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
159: }
160: }
161: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
164: /* Test MatGetSize(), MatGetLocalSize() */
165: MatGetSize(sA, &i,&j); MatGetSize(A, &i2,&j2);
166: i -= i2; j -= j2;
167: if (i || j) {
168: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
169: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
170: }
172: MatGetLocalSize(sA, &i,&j); MatGetLocalSize(A, &i2,&j2);
173: i2 -= i; j2 -= j;
174: if (i2 || j2) {
175: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
176: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
177: }
179: /* vectors */
180: /*--------------------*/
181: /* i is obtained from MatGetLocalSize() */
182: VecCreate(PETSC_COMM_WORLD,&x);
183: VecSetSizes(x,i,PETSC_DECIDE);
184: VecSetFromOptions(x);
185: VecDuplicate(x,&y);
186: VecDuplicate(x,&u);
187: VecDuplicate(x,&s1);
188: VecDuplicate(x,&s2);
190: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
191: PetscRandomSetFromOptions(rctx);
192: VecSetRandom(x,rctx);
193: VecSet(u,one);
195: /* Test MatNorm() */
196: MatNorm(A,NORM_FROBENIUS,&r1);
197: MatNorm(sA,NORM_FROBENIUS,&r2);
198: rnorm = PetscAbsReal(r1-r2)/r2;
199: if (rnorm > tol && !rank) {
200: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
201: }
202: MatNorm(A,NORM_INFINITY,&r1);
203: MatNorm(sA,NORM_INFINITY,&r2);
204: rnorm = PetscAbsReal(r1-r2)/r2;
205: if (rnorm > tol && !rank) {
206: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
207: }
208: MatNorm(A,NORM_1,&r1);
209: MatNorm(sA,NORM_1,&r2);
210: rnorm = PetscAbsReal(r1-r2)/r2;
211: if (rnorm > tol && !rank) {
212: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
213: }
215: /* Test MatGetOwnershipRange() */
216: MatGetOwnershipRange(sA,&rstart,&rend);
217: MatGetOwnershipRange(A,&i2,&j2);
218: i2 -= rstart; j2 -= rend;
219: if (i2 || j2) {
220: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
221: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
222: }
224: /* Test MatDiagonalScale() */
225: MatDiagonalScale(A,x,x);
226: MatDiagonalScale(sA,x,x);
227: MatMultEqual(A,sA,10,&flg);
228: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
230: /* Test MatGetDiagonal(), MatScale() */
231: MatGetDiagonal(A,s1);
232: MatGetDiagonal(sA,s2);
233: VecNorm(s1,NORM_1,&r1);
234: VecNorm(s2,NORM_1,&r2);
235: r1 -= r2;
236: if (r1<-tol || r1>tol) {
237: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
238: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
239: }
241: MatScale(A,alpha);
242: MatScale(sA,alpha);
244: /* Test MatGetRowMaxAbs() */
245: MatGetRowMaxAbs(A,s1,NULL);
246: MatGetRowMaxAbs(sA,s2,NULL);
248: VecNorm(s1,NORM_1,&r1);
249: VecNorm(s2,NORM_1,&r2);
250: r1 -= r2;
251: if (r1<-tol || r1>tol) {
252: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
253: }
255: /* Test MatMult(), MatMultAdd() */
256: MatMultEqual(A,sA,10,&flg);
257: if (!flg) {
258: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
259: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
260: }
262: MatMultAddEqual(A,sA,10,&flg);
263: if (!flg) {
264: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
265: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
266: }
268: /* Test MatMultTranspose(), MatMultTransposeAdd() */
269: for (i=0; i<10; i++) {
270: VecSetRandom(x,rctx);
271: MatMultTranspose(A,x,s1);
272: MatMultTranspose(sA,x,s2);
273: VecNorm(s1,NORM_1,&r1);
274: VecNorm(s2,NORM_1,&r2);
275: r1 -= r2;
276: if (r1<-tol || r1>tol) {
277: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
278: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
279: }
280: }
281: for (i=0; i<10; i++) {
282: VecSetRandom(x,rctx);
283: VecSetRandom(y,rctx);
284: MatMultTransposeAdd(A,x,y,s1);
285: MatMultTransposeAdd(sA,x,y,s2);
286: VecNorm(s1,NORM_1,&r1);
287: VecNorm(s2,NORM_1,&r2);
288: r1 -= r2;
289: if (r1<-tol || r1>tol) {
290: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
291: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
292: }
293: }
294: /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD); */
295: /* MatView(sA, PETSC_VIEWER_DRAW_WORLD); */
297: /* Test MatDuplicate() */
298: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
299: MatEqual(sA,sB,&flg);
300: if (!flg) {
301: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
302: }
303: MatMultEqual(sA,sB,5,&flg);
304: if (!flg) {
305: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
306: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
307: }
308: MatMultAddEqual(sA,sB,5,&flg);
309: if (!flg) {
310: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
311: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
312: }
313: MatDestroy(&sB);
315: VecDestroy(&u);
316: VecDestroy(&x);
317: VecDestroy(&y);
318: VecDestroy(&s1);
319: VecDestroy(&s2);
320: MatDestroy(&sA);
321: MatDestroy(&A);
322: PetscRandomDestroy(&rctx);
324: PetscFinalize();
325: return 0;
326: }