Actual source code: ex74.c
1: /*$Id: ex74.c,v 1.47 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests the various sequential routines in MatSBAIJ format.\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Vec x,y,b,s1,s2;
12: Mat A; /* linear system matrix */
13: Mat sA,sC; /* symmetric part of the matrices */
14: int n,mbs=16,bs=1,nz=3,prob=1;
15: int ierr,i,j,col[3],size,block, row,I,J,n1,*ip_ptr,inc;
16: int lf; /* level of fill for icc */
17: int *cols1,*cols2;
18: PetscReal norm1,norm2,tol=1.e-10;
19: PetscScalar neg_one = -1.0,four=4.0,value[3],alpha=0.1;
20: PetscScalar *vr1,*vr2,*vr1_wk,*vr2_wk;
21: IS perm, isrow, iscol;
22: PetscRandom rand;
23: PetscTruth getrow=PETSC_FALSE;
24: MatInfo minfo1,minfo2;
25: MatFactorInfo factinfo;
27: PetscInitialize(&argc,&args,(char *)0,help);
28: MPI_Comm_size(PETSC_COMM_WORLD,&size);
29: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
30: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
31: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
32: PetscOptionsHasName(PETSC_NULL,"-testMatGetRow",&getrow);
34: n = mbs*bs;
35: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
36: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);
38: /* Test MatGetOwnershipRange() */
39: MatGetOwnershipRange(A,&I,&J);
40: MatGetOwnershipRange(sA,&i,&j);
41: if (i-I || j-J){
42: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
43: }
45: /* Assemble matrix */
46: if (bs == 1){
47: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
48: if (prob == 1){ /* tridiagonal matrix */
49: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
50: for (i=1; i<n-1; i++) {
51: col[0] = i-1; col[1] = i; col[2] = i+1;
52: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
53: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
54: }
55: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
56: value[0]= 0.1; value[1]=-1; value[2]=2;
57: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
58: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
60: i = 0;
61: col[0] = n-1; col[1] = 1; col[2]=0;
62: value[0] = 0.1; value[1] = -1.0; value[2]=2;
63: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
64: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
65: }
66: else if (prob ==2){ /* matrix for the five point stencil */
67: n1 = (int) (sqrt((PetscReal)n) + 0.001);
68: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
69: for (i=0; i<n1; i++) {
70: for (j=0; j<n1; j++) {
71: I = j + n1*i;
72: if (i>0) {
73: J = I - n1;
74: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
75: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
76: }
77: if (i<n1-1) {
78: J = I + n1;
79: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
80: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
81: }
82: if (j>0) {
83: J = I - 1;
84: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
85: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
86: }
87: if (j<n1-1) {
88: J = I + 1;
89: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
90: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
91: }
92: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
93: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
94: }
95: }
96: }
97: }
98: else { /* bs > 1 */
99: for (block=0; block<n/bs; block++){
100: /* diagonal blocks */
101: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
102: for (i=1+block*bs; i<bs-1+block*bs; i++) {
103: col[0] = i-1; col[1] = i; col[2] = i+1;
104: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
105: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
106: }
107: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
108: value[0]=-1.0; value[1]=4.0;
109: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
110: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
112: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
113: value[0]=4.0; value[1] = -1.0;
114: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
115: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
116: }
117: /* off-diagonal blocks */
118: value[0]=-1.0;
119: for (i=0; i<(n/bs-1)*bs; i++){
120: col[0]=i+bs;
121: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
122: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
123: col[0]=i; row=i+bs;
124: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
125: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
126: }
127: }
128: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
130: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
131: MatView(A, PETSC_VIEWER_DRAW_WORLD);
132: MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
134: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
136: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
137: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
138: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
139: */
141: /* Test MatNorm(), MatDuplicate() */
142: MatNorm(A,NORM_FROBENIUS,&norm1);
143: MatDuplicate(sA,MAT_COPY_VALUES,&sC);
144: MatNorm(sC,NORM_FROBENIUS,&norm2);
145: MatDestroy(sC);
146: norm1 -= norm2;
147: if (norm1<-tol || norm1>tol){
148: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
149: }
150: MatNorm(A,NORM_INFINITY,&norm1);
151: MatNorm(sA,NORM_INFINITY,&norm2);
152: norm1 -= norm2;
153: if (norm1<-tol || norm1>tol){
154: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
155: }
156:
157: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
158: MatGetInfo(A,MAT_LOCAL,&minfo1);
159: MatGetInfo(sA,MAT_LOCAL,&minfo2);
160: /*
161: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
162: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
163: */
164: i = (int) (minfo1.nz_used - minfo2.nz_used);
165: j = (int) (minfo2.nz_allocated - minfo2.nz_used);
166: if (i<0 || j<0) {
167: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
168: }
170: MatGetSize(A,&I,&J);
171: MatGetSize(sA,&i,&j);
172: if (i-I || j-J) {
173: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
174: }
175:
176: MatGetBlockSize(A, &I);
177: MatGetBlockSize(sA, &i);
178: if (i-I){
179: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
180: }
182: /* Test MatGetRow() */
183: if (getrow){
184: row = n/2;
185: PetscMalloc(n*sizeof(PetscScalar),&vr1);
186: vr1_wk = vr1;
187: PetscMalloc(n*sizeof(PetscScalar),&vr2);
188: vr2_wk = vr2;
189: MatGetRow(A,row,&J,&cols1,&vr1);
190: vr1_wk += J-1;
191: MatGetRow(sA,row,&j,&cols2,&vr2);
192: vr2_wk += j-1;
193: VecCreateSeq(PETSC_COMM_SELF,j,&x);
194:
195: for (i=j-1; i>-1; i--){
196: VecSetValue(x,i,*vr2_wk - *vr1_wk,INSERT_VALUES);
197: vr2_wk--; vr1_wk--;
198: }
199: VecNorm(x,NORM_1,&norm2);
200: if (norm2<-tol || norm2>tol) {
201: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRow()\n");
202: }
203: VecDestroy(x);
204: MatRestoreRow(A,row,&J,&cols1,&vr1);
205: MatRestoreRow(sA,row,&j,&cols2,&vr2);
206: PetscFree(vr1);
207: PetscFree(vr2);
209: /* Test GetSubMatrix() */
210: /* get a submatrix consisting of every next block row and column of the original matrix */
211: /* for symm. matrix, iscol=isrow. */
212: PetscMalloc(n*sizeof(IS),&isrow);
213: PetscMalloc(n*sizeof(int),&ip_ptr);
214: j = 0;
215: for (n1=0; n1<mbs; n1 += 2){ /* n1: block row */
216: for (i=0; i<bs; i++) ip_ptr[j++] = n1*bs + i;
217: }
218: ISCreateGeneral(PETSC_COMM_SELF, j, ip_ptr, &isrow);
219: /* ISView(isrow, PETSC_VIEWER_STDOUT_SELF); */
220:
221: MatGetSubMatrix(sA,isrow,isrow,PETSC_DECIDE,MAT_INITIAL_MATRIX,&sC);
222: ISDestroy(isrow);
223: PetscFree(ip_ptr);
224: printf("sA =\n");
225: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
226: printf("submatrix of sA =\n");
227: MatView(sC,PETSC_VIEWER_STDOUT_WORLD);
228: MatDestroy(sC);
229: }
231: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
232: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
233: VecCreateSeq(PETSC_COMM_SELF,n,&x);
234: VecDuplicate(x,&s1);
235: VecDuplicate(x,&s2);
236: VecDuplicate(x,&y);
237: VecDuplicate(x,&b);
238: VecSetRandom(rand,x);
240: MatDiagonalScale(A,x,x);
241: MatDiagonalScale(sA,x,x);
243: MatGetDiagonal(A,s1);
244: MatGetDiagonal(sA,s2);
245:
246: VecAXPY(&neg_one,s1,s2);
247: VecNorm(s2,NORM_1,&norm1);
248: if ( norm1>tol) {
249: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%g\n",norm1);
250: }
251: /*
252: VecNorm(s1,NORM_1,&norm1);
253: VecNorm(s2,NORM_1,&norm2);
254: norm1 -= norm2;
255: if (norm1<-tol || norm1>tol) {
256: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
257: }
258: */
260: MatScale(&alpha,A);
261: MatScale(&alpha,sA);
263: /* Test MatGetRowMax() */
264: MatGetRowMax(A,s1);
265: MatGetRowMax(sA,s2);
266: VecNorm(s1,NORM_1,&norm1);
267: VecNorm(s2,NORM_1,&norm2);
268: norm1 -= norm2;
269: if (norm1<-tol || norm1>tol) {
270: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMax() \n");
271: }
273: /* Test MatMult() */
274: for (i=0; i<40; i++) {
275: VecSetRandom(rand,x);
276: MatMult(A,x,s1);
277: MatMult(sA,x,s2);
278: VecNorm(s1,NORM_1,&norm1);
279: VecNorm(s2,NORM_1,&norm2);
280: norm1 -= norm2;
281: if (norm1<-tol || norm1>tol) {
282: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %g\n",norm1);
283: }
284: }
286: /* MatMultAdd() */
287: for (i=0; i<40; i++) {
288: VecSetRandom(rand,x);
289: VecSetRandom(rand,y);
290: MatMultAdd(A,x,y,s1);
291: MatMultAdd(sA,x,y,s2);
292: VecNorm(s1,NORM_1,&norm1);
293: VecNorm(s2,NORM_1,&norm2);
294: norm1 -= norm2;
295: if (norm1<-tol || norm1>tol) {
296: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %g\n",norm1);
297: }
298: }
300: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
301: MatGetOrdering(A,MATORDERING_NATURAL,&perm,&iscol);
302: ISDestroy(iscol);
303: norm1 = tol;
304: inc = bs;
306: /* initialize factinfo */
307: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
309: for (lf=-1; lf<10; lf += inc){
310: if (lf==-1) { /* Cholesky factor */
311: factinfo.fill = 5.0;
312: MatCholeskyFactorSymbolic(sA,perm,&factinfo,&sC);
313: } else { /* incomplete Cholesky factor */
314: factinfo.fill = 5.0;
315: factinfo.levels = lf;
316: MatICCFactorSymbolic(sA,perm,&factinfo,&sC);
317: }
318: MatCholeskyFactorNumeric(sA,&sC);
319: /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */
321: /* test MatGetDiagonal on numeric factor */
322: /*
323: if (lf == -1) {
324: MatGetDiagonal(sC,s1);
325: printf(" in ex74.c, diag: \n");
326: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
327: }
328: */
330: MatMult(sA,x,b);
331: MatSolve(sC,b,y);
332: MatDestroy(sC);
333:
334: /* Check the error */
335: VecAXPY(&neg_one,x,y);
336: VecNorm(y,NORM_2,&norm2);
337: /* printf("lf: %d, error: %g\n", lf,norm2); */
338: if (10*norm1 < norm2 && lf-inc != -1){
339: PetscPrintf(PETSC_COMM_SELF,"lf=%d, %d, Norm of error=%g, %g\n",lf-inc,lf,norm1,norm2);
340: }
341: norm1 = norm2;
342: if (norm2 < tol && lf != -1) break;
343: }
345: ISDestroy(perm);
347: MatDestroy(A);
348: MatDestroy(sA);
349: VecDestroy(x);
350: VecDestroy(y);
351: VecDestroy(s1);
352: VecDestroy(s2);
353: VecDestroy(b);
354: PetscRandomDestroy(rand);
356: PetscFinalize();
357: return 0;
358: }