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);
33: n = mbs*bs;
34: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
35: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);
37: /* Test MatGetOwnershipRange() */
38: MatGetOwnershipRange(A,&I,&J);
39: MatGetOwnershipRange(sA,&i,&j);
40: if (i-I || j-J){
41: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
42: }
44: /* Assemble matrix */
45: if (bs == 1){
46: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
47: if (prob == 1){ /* tridiagonal matrix */
48: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
49: for (i=1; i<n-1; i++) {
50: col[0] = i-1; col[1] = i; col[2] = i+1;
51: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
52: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
53: }
54: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
55: value[0]= 0.1; value[1]=-1; value[2]=2;
56: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
57: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
59: i = 0;
60: col[0] = n-1; col[1] = 1; col[2]=0;
61: value[0] = 0.1; value[1] = -1.0; value[2]=2;
62: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
63: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
64: }
65: else if (prob ==2){ /* matrix for the five point stencil */
66: n1 = (int) (sqrt((PetscReal)n) + 0.001);
67: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
68: for (i=0; i<n1; i++) {
69: for (j=0; j<n1; j++) {
70: I = j + n1*i;
71: if (i>0) {
72: J = I - n1;
73: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
74: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
75: }
76: if (i<n1-1) {
77: J = I + n1;
78: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
79: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
80: }
81: if (j>0) {
82: J = I - 1;
83: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
84: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
85: }
86: if (j<n1-1) {
87: J = I + 1;
88: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
89: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
90: }
91: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
92: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
93: }
94: }
95: }
96: }
97: else { /* bs > 1 */
98: for (block=0; block<n/bs; block++){
99: /* diagonal blocks */
100: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
101: for (i=1+block*bs; i<bs-1+block*bs; i++) {
102: col[0] = i-1; col[1] = i; col[2] = i+1;
103: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
104: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
105: }
106: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
107: value[0]=-1.0; value[1]=4.0;
108: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
109: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
111: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
112: value[0]=4.0; value[1] = -1.0;
113: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
114: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
115: }
116: /* off-diagonal blocks */
117: value[0]=-1.0;
118: for (i=0; i<(n/bs-1)*bs; i++){
119: col[0]=i+bs;
120: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
121: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
122: col[0]=i; row=i+bs;
123: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
124: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
125: }
126: }
127: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
129: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
130: MatView(A, PETSC_VIEWER_DRAW_WORLD);
131: MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
133: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
134: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
135: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
136: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
137: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
138: */
140: /* Test MatNorm(), MatDuplicate() */
141: MatNorm(A,NORM_FROBENIUS,&norm1);
142: MatDuplicate(sA,MAT_COPY_VALUES,&sC);
143: MatNorm(sC,NORM_FROBENIUS,&norm2);
144: MatDestroy(sC);
145: norm1 -= norm2;
146: if (norm1<-tol || norm1>tol){
147: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
148: }
149: MatNorm(A,NORM_INFINITY,&norm1);
150: MatNorm(sA,NORM_INFINITY,&norm2);
151: norm1 -= norm2;
152: if (norm1<-tol || norm1>tol){
153: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
154: }
155:
156: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
157: MatGetInfo(A,MAT_LOCAL,&minfo1);
158: MatGetInfo(sA,MAT_LOCAL,&minfo2);
159: /*
160: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
161: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
162: */
163: i = (int) (minfo1.nz_used - minfo2.nz_used);
164: j = (int) (minfo2.nz_allocated - minfo2.nz_used);
165: if (i<0 || j<0) {
166: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
167: }
169: MatGetSize(A,&I,&J);
170: MatGetSize(sA,&i,&j);
171: if (i-I || j-J) {
172: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
173: }
174:
175: MatGetBlockSize(A, &I);
176: MatGetBlockSize(sA, &i);
177: if (i-I){
178: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
179: }
181: /* Test MatGetRow() */
182: if (getrow){
183: row = n/2;
184: PetscMalloc(n*sizeof(PetscScalar),&vr1);
185: vr1_wk = vr1;
186: PetscMalloc(n*sizeof(PetscScalar),&vr2);
187: vr2_wk = vr2;
188: MatGetRow(A,row,&J,&cols1,&vr1);
189: vr1_wk += J-1;
190: MatGetRow(sA,row,&j,&cols2,&vr2);
191: vr2_wk += j-1;
192: VecCreateSeq(PETSC_COMM_SELF,j,&x);
193:
194: for (i=j-1; i>-1; i--){
195: VecSetValue(x,i,*vr2_wk - *vr1_wk,INSERT_VALUES);
196: vr2_wk--; vr1_wk--;
197: }
198: VecNorm(x,NORM_1,&norm2);
199: if (norm2<-tol || norm2>tol) {
200: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRow()\n");
201: }
202: VecDestroy(x);
203: MatRestoreRow(A,row,&J,&cols1,&vr1);
204: MatRestoreRow(sA,row,&j,&cols2,&vr2);
205: PetscFree(vr1);
206: PetscFree(vr2);
208: /* Test GetSubMatrix() */
209: /* get a submatrix consisting of every next block row and column of the original matrix */
210: /* for symm. matrix, iscol=isrow. */
211: PetscMalloc(n*sizeof(IS),&isrow);
212: PetscMalloc(n*sizeof(int),&ip_ptr);
213: j = 0;
214: for (n1=0; n1<mbs; n1 += 2){ /* n1: block row */
215: for (i=0; i<bs; i++) ip_ptr[j++] = n1*bs + i;
216: }
217: ISCreateGeneral(PETSC_COMM_SELF, j, ip_ptr, &isrow);
218: /* ISView(isrow, PETSC_VIEWER_STDOUT_SELF); */
219:
220: MatGetSubMatrix(sA,isrow,isrow,PETSC_DECIDE,MAT_INITIAL_MATRIX,&sC);
221: ISDestroy(isrow);
222: PetscFree(ip_ptr);
223: printf("sA =\n");
224: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
225: printf("submatrix of sA =\n");
226: MatView(sC,PETSC_VIEWER_STDOUT_WORLD);
227: MatDestroy(sC);
228: }
230: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
231: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
232: VecCreateSeq(PETSC_COMM_SELF,n,&x);
233: VecDuplicate(x,&s1);
234: VecDuplicate(x,&s2);
235: VecDuplicate(x,&y);
236: VecDuplicate(x,&b);
237: VecSetRandom(rand,x);
239: MatDiagonalScale(A,x,x);
240: MatDiagonalScale(sA,x,x);
242: MatGetDiagonal(A,s1);
243: MatGetDiagonal(sA,s2);
244: VecNorm(s1,NORM_1,&norm1);
245: VecNorm(s2,NORM_1,&norm2);
246: norm1 -= norm2;
247: if (norm1<-tol || norm1>tol) {
248: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
249: }
251: MatScale(&alpha,A);
252: MatScale(&alpha,sA);
254: /* Test MatGetRowMax() */
255: MatGetRowMax(A,s1);
256: MatGetRowMax(sA,s2);
257: VecNorm(s1,NORM_1,&norm1);
258: VecNorm(s2,NORM_1,&norm2);
259: norm1 -= norm2;
260: if (norm1<-tol || norm1>tol) {
261: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMax() \n");
262: }
264: /* Test MatMult(), MatMultAdd() */
265: for (i=0; i<40; i++) {
266: VecSetRandom(rand,x);
267: MatMult(A,x,s1);
268: MatMult(sA,x,s2);
269: VecNorm(s1,NORM_1,&norm1);
270: VecNorm(s2,NORM_1,&norm2);
271: norm1 -= norm2;
272: if (norm1<-tol || norm1>tol) {
273: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
274: }
275: }
277: for (i=0; i<40; i++) {
278: VecSetRandom(rand,x);
279: VecSetRandom(rand,y);
280: MatMultAdd(A,x,y,s1);
281: MatMultAdd(sA,x,y,s2);
282: VecNorm(s1,NORM_1,&norm1);
283: VecNorm(s2,NORM_1,&norm2);
284: norm1 -= norm2;
285: if (norm1<-tol || norm1>tol) {
286: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
287: }
288: }
290: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
291: MatGetOrdering(A,MATORDERING_NATURAL,&perm,&iscol);
292: ISDestroy(iscol);
293: norm1 = tol;
294: inc = bs;
296: /* initialize factinfo */
297: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
299: for (lf=-1; lf<10; lf += inc){
300: if (lf==-1) { /* Cholesky factor */
301: factinfo.fill = 5.0;
302: MatCholeskyFactorSymbolic(sA,perm,&factinfo,&sC);
303: } else { /* incomplete Cholesky factor */
304: factinfo.fill = 5.0;
305: factinfo.levels = lf;
306: MatICCFactorSymbolic(sA,perm,&factinfo,&sC);
307: }
308: MatCholeskyFactorNumeric(sA,&sC);
309: /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */
311: /* test MatGetDiagonal on numeric factor */
312: /*
313: if (lf == -1) {
314: MatGetDiagonal(sC,s1);
315: printf(" in ex74.c, diag: \n");
316: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
317: }
318: */
320: MatMult(sA,x,b);
321: MatSolve(sC,b,y);
322: MatDestroy(sC);
323:
324: /* Check the error */
325: VecAXPY(&neg_one,x,y);
326: VecNorm(y,NORM_2,&norm2);
327: /* printf("lf: %d, error: %g\n", lf,norm2); */
328: if (10*norm1 < norm2 && lf-inc != -1){
329: PetscPrintf(PETSC_COMM_SELF,"lf=%d, %d, Norm of error=%g, %g\n",lf-inc,lf,norm1,norm2);
330: }
331: norm1 = norm2;
332: if (norm2 < tol && lf != -1) break;
333: }
335: ISDestroy(perm);
337: MatDestroy(A);
338: MatDestroy(sA);
339: VecDestroy(x);
340: VecDestroy(y);
341: VecDestroy(s1);
342: VecDestroy(s2);
343: VecDestroy(b);
344: PetscRandomDestroy(rand);
346: PetscFinalize();
347: return 0;
348: }