Actual source code: ex75.c
2: static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Vec x,y,u,s1,s2;
9: Mat A,sA,sB;
10: PetscRandom rctx;
11: PetscReal r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
12: PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
13: 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=1;
14: PetscMPIInt size,rank;
15: PetscBool flg;
17: PetscInitialize(&argc,&args,(char*)0,help);
18: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
19: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL);
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: /* Create a BAIJ matrix A */
26: n = mbs*bs;
27: MatCreate(PETSC_COMM_WORLD,&A);
28: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
29: MatSetType(A,MATBAIJ);
30: MatSetFromOptions(A);
31: MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL);
32: MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL);
33: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
35: if (bs == 1) {
36: if (prob == 1) { /* tridiagonal matrix */
37: value[0] = -1.0; value[1] = 0.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: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
41: }
42: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
43: value[0]= 0.1; value[1]=-1.0; value[2]=0.0;
44: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
46: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
47: value[0] = 0.0; value[1] = -1.0; value[2]=0.1;
48: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
49: } else if (prob ==2) { /* matrix for the five point stencil */
50: n1 = (int) PetscSqrtReal((PetscReal)n);
51: for (i=0; i<n1; i++) {
52: for (j=0; j<n1; j++) {
53: Ii = j + n1*i;
54: if (i>0) {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
55: if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
56: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
57: if (j<n1-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
58: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
59: }
60: }
61: }
62: /* end of if (bs == 1) */
63: } else { /* bs > 1 */
64: for (block=0; block<n/bs; block++) {
65: /* diagonal blocks */
66: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
67: for (i=1+block*bs; i<bs-1+block*bs; i++) {
68: col[0] = i-1; col[1] = i; col[2] = i+1;
69: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
70: }
71: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
72: value[0]=-1.0; value[1]=4.0;
73: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
75: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
76: value[0]=4.0; value[1] = -1.0;
77: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
78: }
79: /* off-diagonal blocks */
80: value[0]=-1.0;
81: for (i=0; i<(n/bs-1)*bs; i++) {
82: col[0]=i+bs;
83: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
84: col[0]=i; row=i+bs;
85: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
86: }
87: }
88: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
90: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
92: /* Get SBAIJ matrix sA from A */
93: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
95: /* Test MatGetSize(), MatGetLocalSize() */
96: MatGetSize(sA, &i,&j);
97: MatGetSize(A, &i2,&j2);
98: i -= i2; j -= j2;
99: if (i || j) {
100: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
101: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
102: }
104: MatGetLocalSize(sA, &i,&j);
105: MatGetLocalSize(A, &i2,&j2);
106: i2 -= i; j2 -= j;
107: if (i2 || j2) {
108: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
109: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
110: }
112: /* vectors */
113: /*--------------------*/
114: /* i is obtained from MatGetLocalSize() */
115: VecCreate(PETSC_COMM_WORLD,&x);
116: VecSetSizes(x,i,PETSC_DECIDE);
117: VecSetFromOptions(x);
118: VecDuplicate(x,&y);
119: VecDuplicate(x,&u);
120: VecDuplicate(x,&s1);
121: VecDuplicate(x,&s2);
123: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
124: PetscRandomSetFromOptions(rctx);
125: VecSetRandom(x,rctx);
126: VecSet(u,one);
128: /* Test MatNorm() */
129: MatNorm(A,NORM_FROBENIUS,&r1);
130: MatNorm(sA,NORM_FROBENIUS,&r2);
131: rnorm = PetscAbsReal(r1-r2)/r2;
132: if (rnorm > tol && rank == 0) {
133: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
134: }
135: MatNorm(A,NORM_INFINITY,&r1);
136: MatNorm(sA,NORM_INFINITY,&r2);
137: rnorm = PetscAbsReal(r1-r2)/r2;
138: if (rnorm > tol && rank == 0) {
139: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
140: }
141: MatNorm(A,NORM_1,&r1);
142: MatNorm(sA,NORM_1,&r2);
143: rnorm = PetscAbsReal(r1-r2)/r2;
144: if (rnorm > tol && rank == 0) {
145: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs);
146: }
148: /* Test MatGetOwnershipRange() */
149: MatGetOwnershipRange(sA,&rstart,&rend);
150: MatGetOwnershipRange(A,&i2,&j2);
151: i2 -= rstart; j2 -= rend;
152: if (i2 || j2) {
153: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
154: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
155: }
157: /* Test MatDiagonalScale() */
158: MatDiagonalScale(A,x,x);
159: MatDiagonalScale(sA,x,x);
160: MatMultEqual(A,sA,10,&flg);
163: /* Test MatGetDiagonal(), MatScale() */
164: MatGetDiagonal(A,s1);
165: MatGetDiagonal(sA,s2);
166: VecNorm(s1,NORM_1,&r1);
167: VecNorm(s2,NORM_1,&r2);
168: r1 -= r2;
169: if (r1<-tol || r1>tol) {
170: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
171: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
172: }
174: MatScale(A,alpha);
175: MatScale(sA,alpha);
177: /* Test MatGetRowMaxAbs() */
178: MatGetRowMaxAbs(A,s1,NULL);
179: MatGetRowMaxAbs(sA,s2,NULL);
181: VecNorm(s1,NORM_1,&r1);
182: VecNorm(s2,NORM_1,&r2);
183: r1 -= r2;
184: if (r1<-tol || r1>tol) {
185: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
186: }
188: /* Test MatMult(), MatMultAdd() */
189: MatMultEqual(A,sA,10,&flg);
190: if (!flg) {
191: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
192: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
193: }
195: MatMultAddEqual(A,sA,10,&flg);
196: if (!flg) {
197: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
198: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
199: }
201: /* Test MatMultTranspose(), MatMultTransposeAdd() */
202: for (i=0; i<10; i++) {
203: VecSetRandom(x,rctx);
204: MatMultTranspose(A,x,s1);
205: MatMultTranspose(sA,x,s2);
206: VecNorm(s1,NORM_1,&r1);
207: VecNorm(s2,NORM_1,&r2);
208: r1 -= r2;
209: if (r1<-tol || r1>tol) {
210: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
211: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
212: }
213: }
214: for (i=0; i<10; i++) {
215: VecSetRandom(x,rctx);
216: VecSetRandom(y,rctx);
217: MatMultTransposeAdd(A,x,y,s1);
218: MatMultTransposeAdd(sA,x,y,s2);
219: VecNorm(s1,NORM_1,&r1);
220: VecNorm(s2,NORM_1,&r2);
221: r1 -= r2;
222: if (r1<-tol || r1>tol) {
223: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
224: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
225: }
226: }
228: /* Test MatDuplicate() */
229: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
230: MatEqual(sA,sB,&flg);
231: if (!flg) {
232: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
233: }
234: MatMultEqual(sA,sB,5,&flg);
235: if (!flg) {
236: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
237: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
238: }
239: MatMultAddEqual(sA,sB,5,&flg);
240: if (!flg) {
241: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
242: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
243: }
244: MatDestroy(&sB);
245: VecDestroy(&u);
246: VecDestroy(&x);
247: VecDestroy(&y);
248: VecDestroy(&s1);
249: VecDestroy(&s2);
250: MatDestroy(&sA);
251: MatDestroy(&A);
252: PetscRandomDestroy(&rctx);
253: PetscFinalize();
254: return 0;
255: }
257: /*TEST
259: test:
260: nsize: {{1 3}}
261: args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
263: TEST*/