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;
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;
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: MatCreateMPISBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&sA);
35: if (bs == 1){
36: if (prob == 1){ /* tridiagonal matrix */
37: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
38: for (i=1; i<n-1; i++) {
39: col[0] = i-1; col[1] = i; col[2] = i+1;
40: /* PetscTrValid(0,0,0,0); */
41: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
42: }
43: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
44: value[0]= 0.1; value[1]=-1; value[2]=2;
45: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
47: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
48: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
49: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
50: }
51: else if (prob ==2){ /* matrix for the five point stencil */
52: n1 = (int) sqrt(n);
53: if (n1*n1 != n){
54: SETERRQ(PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
55: }
56:
57: for (i=0; i<n1; i++) {
58: for (j=0; j<n1; j++) {
59: I = j + n1*i;
60: if (i>0) {J = I - n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
61: if (i<n1-1) {J = I + n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
62: if (j>0) {J = I - 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
63: if (j<n1-1) {J = I + 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
64: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
65: }
66: }
67: }
68: } /* end of if (bs == 1) */
69: else { /* bs > 1 */
70: for (block=0; block<n/bs; block++){
71: /* diagonal blocks */
72: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
73: for (i=1+block*bs; i<bs-1+block*bs; i++) {
74: col[0] = i-1; col[1] = i; col[2] = i+1;
75: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
76: }
77: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
78: value[0]=-1.0; value[1]=4.0;
79: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
81: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
82: value[0]=4.0; value[1] = -1.0;
83: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
84: }
85: /* off-diagonal blocks */
86: value[0]=-1.0;
87: for (i=0; i<(n/bs-1)*bs; i++){
88: col[0]=i+bs;
89: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
90: col[0]=i; row=i+bs;
91: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
92: }
93: }
94: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
97: /* Test MatView() */
98: /*
99: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
100: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
101: */
102: /* Assemble MPIBAIJ matrix A */
103: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);
105: if (bs == 1){
106: if (prob == 1){ /* tridiagonal matrix */
107: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
108: for (i=1; i<n-1; i++) {
109: col[0] = i-1; col[1] = i; col[2] = i+1;
110: /* PetscTrValid(0,0,0,0); */
111: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
112: }
113: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
114: value[0]= 0.1; value[1]=-1; value[2]=2;
115: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
117: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
118: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
119: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
120: }
121: else if (prob ==2){ /* matrix for the five point stencil */
122: n1 = (int) sqrt(n);
123: for (i=0; i<n1; i++) {
124: for (j=0; j<n1; j++) {
125: I = j + n1*i;
126: if (i>0) {J = I - n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
127: if (i<n1-1) {J = I + n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
128: if (j>0) {J = I - 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
129: if (j<n1-1) {J = I + 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
130: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
131: }
132: }
133: }
134: } /* end of if (bs == 1) */
135: else { /* bs > 1 */
136: for (block=0; block<n/bs; block++){
137: /* diagonal blocks */
138: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
139: for (i=1+block*bs; i<bs-1+block*bs; i++) {
140: col[0] = i-1; col[1] = i; col[2] = i+1;
141: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
142: }
143: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
144: value[0]=-1.0; value[1]=4.0;
145: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
147: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
148: value[0]=4.0; value[1] = -1.0;
149: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
150: }
151: /* off-diagonal blocks */
152: value[0]=-1.0;
153: for (i=0; i<(n/bs-1)*bs; i++){
154: col[0]=i+bs;
155: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
156: col[0]=i; row=i+bs;
157: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
158: }
159: }
160: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
163: /* Test MatGetSize(), MatGetLocalSize() */
164: MatGetSize(sA, &i,&j); MatGetSize(A, &i2,&j2);
165: i -= i2; j -= j2;
166: if (i || j) {
167: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
168: PetscSynchronizedFlush(PETSC_COMM_WORLD);
169: }
170:
171: MatGetLocalSize(sA, &i,&j); MatGetLocalSize(A, &i2,&j2);
172: i2 -= i; j2 -= j;
173: if (i2 || j2) {
174: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
175: PetscSynchronizedFlush(PETSC_COMM_WORLD);
176: }
178: /* vectors */
179: /*--------------------*/
180: /* i is obtained from MatGetLocalSize() */
181: VecCreate(PETSC_COMM_WORLD,&x);
182: VecSetSizes(x,i,PETSC_DECIDE);
183: VecSetFromOptions(x);
184: VecDuplicate(x,&y);
185: VecDuplicate(x,&u);
186: VecDuplicate(x,&s1);
187: VecDuplicate(x,&s2);
189: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
190: VecSetRandom(rctx,x);
191: VecSet(&one,u);
193: /* Test MatNorm() */
194: MatNorm(A,NORM_FROBENIUS,&r1);
195: MatNorm(sA,NORM_FROBENIUS,&r2);
196: r1 -= r2;
197: if (r1<-tol || r1>tol){
198: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatNorm(), A_fnorm - sA_fnorm = %16.14e\n",rank,r1);
199: PetscSynchronizedFlush(PETSC_COMM_WORLD);
200: }
202: /* Test MatGetOwnershipRange() */
203: MatGetOwnershipRange(sA,&rstart,&rend);
204: MatGetOwnershipRange(A,&i2,&j2);
205: i2 -= rstart; j2 -= rend;
206: if (i2 || j2) {
207: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
208: PetscSynchronizedFlush(PETSC_COMM_WORLD);
209: }
211: /* Test MatGetRow(): can only obtain rows associated with the given processor */
212: for (i=rstart; i<rstart+1; i++) {
213: MatGetRow(sA,i,&ncols,&cols,&vr);
214: /*
215: PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"[%d] get row %d: ",rank,i);
216: for (j=0; j<ncols; j++) {
217: PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"%d %g ",cols[j],vr[j]);
218: }
219: PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"\n");
220: PetscSynchronizedFlush(PETSC_COMM_WORLD);
221: */
222: MatRestoreRow(sA,i,&ncols,&cols,&vr);
223: }
225: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
226: MatDiagonalScale(A,x,x);
227: MatDiagonalScale(sA,x,x);
228: MatGetDiagonal(A,s1);
229: MatGetDiagonal(sA,s2);
230: VecNorm(s1,NORM_1,&r1);
231: VecNorm(s2,NORM_1,&r2);
232: r1 -= r2;
233: if (r1<-tol || r1>tol) {
234: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,r1);
235: PetscSynchronizedFlush(PETSC_COMM_WORLD);
236: }
237:
238: MatScale(&alpha,A);
239: MatScale(&alpha,sA);
241: /* Test MatGetRowMax() */
242: MatGetRowMax(A,s1);
243: MatGetRowMax(sA,s2);
245: VecNorm(s1,NORM_1,&r1);
246: VecNorm(s2,NORM_1,&r2);
247: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD);
248: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], r1: %g, r2: %g\n",rank,r1,r2);
249: PetscSynchronizedFlush(PETSC_COMM_WORLD);
250: */
251: r1 -= r2;
252: if (r1<-tol || r1>tol) {
253: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMax() \n");
254: }
256: /* Test MatMult(), MatMultAdd() */
257: for (i=0; i<10; i++) {
258: VecSetRandom(rctx,x);
259: MatMult(A,x,s1);
260: MatMult(sA,x,s2);
261: VecNorm(s1,NORM_1,&r1);
262: VecNorm(s2,NORM_1,&r2);
263: r1 -= r2;
264: if (r1<-tol || r1>tol) {
265: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,r1);
266: PetscSynchronizedFlush(PETSC_COMM_WORLD);
267: }
268: }
270: for (i=0; i<10; i++) {
271: VecSetRandom(rctx,x);
272: VecSetRandom(rctx,y);
273: MatMultAdd(A,x,y,s1);
274: MatMultAdd(sA,x,y,s2);
275: VecNorm(s1,NORM_1,&r1);
276: VecNorm(s2,NORM_1,&r2);
277: r1 -= r2;
278: if (r1<-tol || r1>tol) {
279: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,r1);
280: PetscSynchronizedFlush(PETSC_COMM_WORLD);
281: }
282: }
284: /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD); */
285: /* MatView(sA, PETSC_VIEWER_DRAW_WORLD); */
286:
287: VecDestroy(u);
288: VecDestroy(x);
289: VecDestroy(y);
290: VecDestroy(s1);
291: VecDestroy(s2);
292: MatDestroy(sA);
293: MatDestroy(A);
294: PetscRandomDestroy(rctx);
295:
296: PetscFinalize();
297: return 0;
298: }