Actual source code: ex75.c
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: PetscTruth 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);
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: }
56: else if (prob ==2){ /* matrix for the five point stencil */
57: n1 = (int) sqrt((double)n);
58: if (n1*n1 != n){
59: SETERRQ(PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
60: }
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: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);
110: if (bs == 1){
111: if (prob == 1){ /* tridiagonal matrix */
112: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
113: for (i=1; i<n-1; i++) {
114: col[0] = i-1; col[1] = i; col[2] = i+1;
115: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
116: }
117: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
118: value[0]= 0.1; value[1]=-1; value[2]=2;
119: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
121: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
122: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
123: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
124: }
125: else if (prob ==2){ /* matrix for the five point stencil */
126: n1 = (int) sqrt((double)n);
127: for (i=0; i<n1; i++) {
128: for (j=0; j<n1; j++) {
129: Ii = j + n1*i;
130: if (i>0) {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
131: if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
132: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
133: if (j<n1-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
134: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
135: }
136: }
137: }
138: } /* end of if (bs == 1) */
139: else { /* bs > 1 */
140: for (block=0; block<n/bs; block++){
141: /* diagonal blocks */
142: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
143: for (i=1+block*bs; i<bs-1+block*bs; i++) {
144: col[0] = i-1; col[1] = i; col[2] = i+1;
145: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
146: }
147: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
148: value[0]=-1.0; value[1]=4.0;
149: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
151: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
152: value[0]=4.0; value[1] = -1.0;
153: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
154: }
155: /* off-diagonal blocks */
156: value[0]=-1.0;
157: for (i=0; i<(n/bs-1)*bs; i++){
158: col[0]=i+bs;
159: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
160: col[0]=i; row=i+bs;
161: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
162: }
163: }
164: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
165: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
167: /* Test MatGetSize(), MatGetLocalSize() */
168: MatGetSize(sA, &i,&j); MatGetSize(A, &i2,&j2);
169: i -= i2; j -= j2;
170: if (i || j) {
171: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
172: PetscSynchronizedFlush(PETSC_COMM_WORLD);
173: }
174:
175: MatGetLocalSize(sA, &i,&j); MatGetLocalSize(A, &i2,&j2);
176: i2 -= i; j2 -= j;
177: if (i2 || j2) {
178: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
179: PetscSynchronizedFlush(PETSC_COMM_WORLD);
180: }
182: /* vectors */
183: /*--------------------*/
184: /* i is obtained from MatGetLocalSize() */
185: VecCreate(PETSC_COMM_WORLD,&x);
186: VecSetSizes(x,i,PETSC_DECIDE);
187: VecSetFromOptions(x);
188: VecDuplicate(x,&y);
189: VecDuplicate(x,&u);
190: VecDuplicate(x,&s1);
191: VecDuplicate(x,&s2);
193: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
194: PetscRandomSetFromOptions(rctx);
195: VecSetRandom(x,rctx);
196: VecSet(u,one);
198: /* Test MatNorm() */
199: MatNorm(A,NORM_FROBENIUS,&r1);
200: MatNorm(sA,NORM_FROBENIUS,&r2);
201: rnorm = PetscAbsScalar(r1-r2)/r2;
202: if (rnorm > tol && !rank){
203: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
204: }
205: MatNorm(A,NORM_INFINITY,&r1);
206: MatNorm(sA,NORM_INFINITY,&r2);
207: rnorm = PetscAbsScalar(r1-r2)/r2;
208: if (rnorm > tol && !rank){
209: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
210: }
211: MatNorm(A,NORM_1,&r1);
212: MatNorm(sA,NORM_1,&r2);
213: rnorm = PetscAbsScalar(r1-r2)/r2;
214: if (rnorm > tol && !rank){
215: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
216: }
218: /* Test MatGetOwnershipRange() */
219: MatGetOwnershipRange(sA,&rstart,&rend);
220: MatGetOwnershipRange(A,&i2,&j2);
221: i2 -= rstart; j2 -= rend;
222: if (i2 || j2) {
223: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
224: PetscSynchronizedFlush(PETSC_COMM_WORLD);
225: }
227: /* Test MatDiagonalScale() */
228: MatDiagonalScale(A,x,x);
229: MatDiagonalScale(sA,x,x);
230: MatMultEqual(A,sA,10,&flg);
231: if (!flg) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
232:
233: /* Test MatGetDiagonal(), MatScale() */
234: MatGetDiagonal(A,s1);
235: MatGetDiagonal(sA,s2);
236: VecNorm(s1,NORM_1,&r1);
237: VecNorm(s2,NORM_1,&r2);
238: r1 -= r2;
239: if (r1<-tol || r1>tol) {
240: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%G \n",rank,r1);
241: PetscSynchronizedFlush(PETSC_COMM_WORLD);
242: }
243:
244: MatScale(A,alpha);
245: MatScale(sA,alpha);
247: /* Test MatGetRowMaxAbs() */
248: MatGetRowMaxAbs(A,s1,PETSC_NULL);
249: MatGetRowMaxAbs(sA,s2,PETSC_NULL);
251: VecNorm(s1,NORM_1,&r1);
252: VecNorm(s2,NORM_1,&r2);
253: r1 -= r2;
254: if (r1<-tol || r1>tol) {
255: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
256: }
258: /* Test MatMult(), MatMultAdd() */
259: MatMultEqual(A,sA,10,&flg);
260: if (!flg){
261: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
262: PetscSynchronizedFlush(PETSC_COMM_WORLD);
263: }
265: MatMultAddEqual(A,sA,10,&flg);
266: if (!flg){
267: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
268: PetscSynchronizedFlush(PETSC_COMM_WORLD);
269: }
271: /* Test MatMultTranspose(), MatMultTransposeAdd() */
272: for (i=0; i<10; i++) {
273: VecSetRandom(x,rctx);
274: MatMultTranspose(A,x,s1);
275: MatMultTranspose(sA,x,s2);
276: VecNorm(s1,NORM_1,&r1);
277: VecNorm(s2,NORM_1,&r2);
278: r1 -= r2;
279: if (r1<-tol || r1>tol) {
280: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%G\n",rank,r1);
281: PetscSynchronizedFlush(PETSC_COMM_WORLD);
282: }
283: }
284: for (i=0; i<10; i++) {
285: VecSetRandom(x,rctx);
286: VecSetRandom(y,rctx);
287: MatMultTransposeAdd(A,x,y,s1);
288: MatMultTransposeAdd(sA,x,y,s2);
289: VecNorm(s1,NORM_1,&r1);
290: VecNorm(s2,NORM_1,&r2);
291: r1 -= r2;
292: if (r1<-tol || r1>tol) {
293: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%G \n",rank,r1);
294: PetscSynchronizedFlush(PETSC_COMM_WORLD);
295: }
296: }
297: /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD); */
298: /* MatView(sA, PETSC_VIEWER_DRAW_WORLD); */
300: /* Test MatDuplicate() */
301: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
302: MatEqual(sA,sB,&flg);
303: if (!flg){
304: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
305: }
306: MatMultEqual(sA,sB,5,&flg);
307: if (!flg){
308: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
309: PetscSynchronizedFlush(PETSC_COMM_WORLD);
310: }
311: MatMultAddEqual(sA,sB,5,&flg);
312: if (!flg){
313: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
314: PetscSynchronizedFlush(PETSC_COMM_WORLD);
315: }
316: MatDestroy(sB);
317:
318: VecDestroy(u);
319: VecDestroy(x);
320: VecDestroy(y);
321: VecDestroy(s1);
322: VecDestroy(s2);
323: MatDestroy(sA);
324: MatDestroy(A);
325: PetscRandomDestroy(rctx);
326:
327: PetscFinalize();
328: return 0;
329: }