Actual source code: ex75.c
1: /*$Id: ex75.c,v 1.28 2001/09/11 16:32:50 bsmith Exp $*/
3: /* Program usage: mpirun -np <procs> ex75 [-help] [all PETSc options] */
5: static char help[] = "Tests the vatious routines in MatMPISBAIJ format.\n";
7: #include petscmat.h
11: int main(int argc,char **args)
12: {
13: int bs=1, mbs=16, d_nz=3, o_nz=3, prob=2;
14: Vec x,y,u,s1,s2;
15: Mat A,sA,sB;
16: PetscRandom rctx;
17: double r1,r2,tol=1.e-10;
18: int i,j,i2,j2,I,J,ierr;
19: PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1,*vr;
20: int n,rank,size,col[3],n1,block,row;
21: int ncols,*cols,rstart,rend;
22: PetscTruth flg;
24: PetscInitialize(&argc,&args,(char *)0,help);
25: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
27:
28: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30:
31: n = mbs*bs;
32:
33: /* Assemble MPISBAIJ matrix sA */
34: MatCreateMPISBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&sA);
36: if (bs == 1){
37: if (prob == 1){ /* tridiagonal matrix */
38: value[0] = -1.0; value[1] = 2.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: /* PetscTrValid(0,0,0,0); */
42: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
43: }
44: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
45: value[0]= 0.1; value[1]=-1; value[2]=2;
46: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
48: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
49: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
50: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
51: }
52: else if (prob ==2){ /* matrix for the five point stencil */
53: n1 = (int) sqrt((double)n);
54: if (n1*n1 != n){
55: SETERRQ(PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
56: }
57:
58: for (i=0; i<n1; i++) {
59: for (j=0; j<n1; j++) {
60: I = j + n1*i;
61: if (i>0) {J = I - n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
62: if (i<n1-1) {J = I + n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
63: if (j>0) {J = I - 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
64: if (j<n1-1) {J = I + 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
65: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
66: }
67: }
68: }
69: } /* end of if (bs == 1) */
70: else { /* bs > 1 */
71: for (block=0; block<n/bs; block++){
72: /* diagonal blocks */
73: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
74: for (i=1+block*bs; i<bs-1+block*bs; i++) {
75: col[0] = i-1; col[1] = i; col[2] = i+1;
76: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
77: }
78: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
79: value[0]=-1.0; value[1]=4.0;
80: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
82: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
83: value[0]=4.0; value[1] = -1.0;
84: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
85: }
86: /* off-diagonal blocks */
87: value[0]=-1.0;
88: for (i=0; i<(n/bs-1)*bs; i++){
89: col[0]=i+bs;
90: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
91: col[0]=i; row=i+bs;
92: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
93: }
94: }
95: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
96: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
98: /* Test MatView() */
99: /*
100: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
101: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
102: */
103: /* Assemble MPIBAIJ matrix A */
104: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);
106: if (bs == 1){
107: if (prob == 1){ /* tridiagonal matrix */
108: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
109: for (i=1; i<n-1; i++) {
110: col[0] = i-1; col[1] = i; col[2] = i+1;
111: /* PetscTrValid(0,0,0,0); */
112: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
113: }
114: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
115: value[0]= 0.1; value[1]=-1; value[2]=2;
116: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
118: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
119: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
120: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
121: }
122: else if (prob ==2){ /* matrix for the five point stencil */
123: n1 = (int) sqrt((double)n);
124: for (i=0; i<n1; i++) {
125: for (j=0; j<n1; j++) {
126: I = j + n1*i;
127: if (i>0) {J = I - n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
128: if (i<n1-1) {J = I + n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
129: if (j>0) {J = I - 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
130: if (j<n1-1) {J = I + 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
131: MatSetValues(A,1,&I,1,&I,&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);
170: }
171:
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);
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,RANDOM_DEFAULT,&rctx);
191: VecSetRandom(rctx,x);
192: VecSet(&one,u);
194: /* Test MatNorm() */
195: MatNorm(A,NORM_FROBENIUS,&r1);
196: MatNorm(sA,NORM_FROBENIUS,&r2);
197: r1 -= r2;
198: if (r1<-tol || r1>tol){
199: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatNorm(), A_fnorm - sA_fnorm = %16.14e\n",rank,r1);
200: PetscSynchronizedFlush(PETSC_COMM_WORLD);
201: }
203: /* Test MatGetOwnershipRange() */
204: MatGetOwnershipRange(sA,&rstart,&rend);
205: MatGetOwnershipRange(A,&i2,&j2);
206: i2 -= rstart; j2 -= rend;
207: if (i2 || j2) {
208: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
209: PetscSynchronizedFlush(PETSC_COMM_WORLD);
210: }
212: /* Test MatGetRow(): can only obtain rows associated with the given processor */
213: for (i=rstart; i<rstart+1; i++) {
214: MatGetRow(sA,i,&ncols,&cols,&vr);
215: MatRestoreRow(sA,i,&ncols,&cols,&vr);
216: }
218: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
219: MatDiagonalScale(A,x,x);
220: MatDiagonalScale(sA,x,x);
221: MatGetDiagonal(A,s1);
222: MatGetDiagonal(sA,s2);
223: VecNorm(s1,NORM_1,&r1);
224: VecNorm(s2,NORM_1,&r2);
225: r1 -= r2;
226: if (r1<-tol || r1>tol) {
227: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,r1);
228: PetscSynchronizedFlush(PETSC_COMM_WORLD);
229: }
230:
231: MatScale(&alpha,A);
232: MatScale(&alpha,sA);
234: /* Test MatGetRowMax() */
235: MatGetRowMax(A,s1);
236: MatGetRowMax(sA,s2);
238: VecNorm(s1,NORM_1,&r1);
239: VecNorm(s2,NORM_1,&r2);
240: r1 -= r2;
241: if (r1<-tol || r1>tol) {
242: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMax() \n");
243: }
245: /* Test MatMult(), MatMultAdd() */
246: for (i=0; i<10; i++) {
247: VecSetRandom(rctx,x);
248: MatMult(A,x,s1);
249: MatMult(sA,x,s2);
250: VecNorm(s1,NORM_1,&r1);
251: VecNorm(s2,NORM_1,&r2);
252: r1 -= r2;
253: if (r1<-tol || r1>tol) {
254: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,r1);
255: PetscSynchronizedFlush(PETSC_COMM_WORLD);
256: }
257: }
259: for (i=0; i<10; i++) {
260: VecSetRandom(rctx,x);
261: VecSetRandom(rctx,y);
262: MatMultAdd(A,x,y,s1);
263: MatMultAdd(sA,x,y,s2);
264: VecNorm(s1,NORM_1,&r1);
265: VecNorm(s2,NORM_1,&r2);
266: r1 -= r2;
267: if (r1<-tol || r1>tol) {
268: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,r1);
269: PetscSynchronizedFlush(PETSC_COMM_WORLD);
270: }
271: }
273: /* Test MatMultTranspose(), MatMultTransposeAdd() */
274: for (i=0; i<10; i++) {
275: VecSetRandom(rctx,x);
276: MatMultTranspose(A,x,s1);
277: MatMultTranspose(sA,x,s2);
278: VecNorm(s1,NORM_1,&r1);
279: VecNorm(s2,NORM_1,&r2);
280: r1 -= r2;
281: if (r1<-tol || r1>tol) {
282: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,r1);
283: PetscSynchronizedFlush(PETSC_COMM_WORLD);
284: }
285: }
286: for (i=0; i<10; i++) {
287: VecSetRandom(rctx,x);
288: VecSetRandom(rctx,y);
289: MatMultTransposeAdd(A,x,y,s1);
290: MatMultTransposeAdd(sA,x,y,s2);
291: VecNorm(s1,NORM_1,&r1);
292: VecNorm(s2,NORM_1,&r2);
293: r1 -= r2;
294: if (r1<-tol || r1>tol) {
295: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,r1);
296: PetscSynchronizedFlush(PETSC_COMM_WORLD);
297: }
298: }
299: /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD); */
300: /* MatView(sA, PETSC_VIEWER_DRAW_WORLD); */
302: /* Test MatDuplicate() */
303: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
304: MatEqual(sA,sB,&flg);
305: if (!flg){
306: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
307: }
308: for (i=0; i<5; i++) {
309: VecSetRandom(rctx,x);
310: MatMult(A,x,s1);
311: MatMult(sB,x,s2);
312: VecNorm(s1,NORM_1,&r1);
313: VecNorm(s2,NORM_1,&r2);
314: r1 -= r2;
315: if (r1<-tol || r1>tol) {
316: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult(), err=%g\n",rank,r1);
317: PetscSynchronizedFlush(PETSC_COMM_WORLD);
318: }
319: }
321: for (i=0; i<5; i++) {
322: VecSetRandom(rctx,x);
323: VecSetRandom(rctx,y);
324: MatMultAdd(A,x,y,s1);
325: MatMultAdd(sB,x,y,s2);
326: VecNorm(s1,NORM_1,&r1);
327: VecNorm(s2,NORM_1,&r2);
328: r1 -= r2;
329: if (r1<-tol || r1>tol) {
330: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(), err=%g \n",rank,r1);
331: PetscSynchronizedFlush(PETSC_COMM_WORLD);
332: }
333: }
334: MatDestroy(sB);
335:
336: VecDestroy(u);
337: VecDestroy(x);
338: VecDestroy(y);
339: VecDestroy(s1);
340: VecDestroy(s2);
341: MatDestroy(sA);
342: MatDestroy(A);
343: PetscRandomDestroy(rctx);
344:
345: PetscFinalize();
346: return 0;
347: }