Actual source code: ex75.c
2: /* Program usage: mpirun -np <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: PetscTruth flg;
21: 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: /* -mat_type <seqsbaij_derived type>, e.g., mpisbaijspooles, sbaijmumps */
37: MatSetFromOptions(sA);
38: MatGetType(sA,&type);
39: /* printf(" mattype: %s\n",type); */
40: MatMPISBAIJSetPreallocation(sA,bs,d_nz,PETSC_NULL,o_nz,PETSC_NULL);
42: if (bs == 1){
43: if (prob == 1){ /* tridiagonal matrix */
44: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
45: for (i=1; i<n-1; i++) {
46: col[0] = i-1; col[1] = i; col[2] = i+1;
47: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
48: }
49: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
50: value[0]= 0.1; value[1]=-1; value[2]=2;
51: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
53: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
54: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
55: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
56: }
57: else if (prob ==2){ /* matrix for the five point stencil */
58: n1 = (int) sqrt((double)n);
59: if (n1*n1 != n){
60: SETERRQ(PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
61: }
62:
63: for (i=0; i<n1; i++) {
64: for (j=0; j<n1; j++) {
65: Ii = j + n1*i;
66: if (i>0) {J = Ii - n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
67: if (i<n1-1) {J = Ii + n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
68: if (j>0) {J = Ii - 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
69: if (j<n1-1) {J = Ii + 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
70: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
71: }
72: }
73: }
74: } /* end of if (bs == 1) */
75: else { /* bs > 1 */
76: for (block=0; block<n/bs; block++){
77: /* diagonal blocks */
78: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
79: for (i=1+block*bs; i<bs-1+block*bs; i++) {
80: col[0] = i-1; col[1] = i; col[2] = i+1;
81: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
82: }
83: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
84: value[0]=-1.0; value[1]=4.0;
85: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
87: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
88: value[0]=4.0; value[1] = -1.0;
89: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
90: }
91: /* off-diagonal blocks */
92: value[0]=-1.0;
93: for (i=0; i<(n/bs-1)*bs; i++){
94: col[0]=i+bs;
95: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
96: col[0]=i; row=i+bs;
97: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
98: }
99: }
100: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
103: /* Test MatView() */
104: /*
105: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
106: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
107: */
108: /* Assemble MPIBAIJ matrix A */
109: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);
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) sqrt((double)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_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 MatGetRowMax() */
249: MatGetRowMax(A,s1);
250: MatGetRowMax(sA,s2);
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: MatGetRowMax() \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: }