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