Actual source code: ex76.c
1: /*$Id: ex76.c,v 1.18 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests matrix permutation for factorization and solve on matrix with MatSBAIJ format. Modified from ex74.c\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Vec x,y,b;
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;
16: int lf; /* level of fill for icc */
17: PetscReal norm1,norm2,tol=1.e-10;
18: PetscScalar neg_one = -1.0,four=4.0,value[3];
19: IS perm;
20: PetscRandom rand;
21: PetscTruth reorder=PETSC_TRUE;
22: MatFactorInfo factinfo;
24: PetscInitialize(&argc,&args,(char *)0,help);
25: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
27: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
28: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
30: n = mbs*bs;
31: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
32: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);
34: /* Test MatGetOwnershipRange() */
35: MatGetOwnershipRange(A,&I,&J);
36: MatGetOwnershipRange(sA,&i,&j);
37: if (i-I || j-J){
38: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
39: }
41: /* Assemble matrix */
42: if (bs == 1){
43: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
44: if (prob == 1){ /* tridiagonal matrix */
45: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
46: for (i=1; i<n-1; i++) {
47: col[0] = i-1; col[1] = i; col[2] = i+1;
48: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
49: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
50: }
51: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
52: value[0]= 0.1; value[1]=-1; value[2]=2;
53: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
54: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
56: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
57: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
58: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
59: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
60: }
61: else if (prob ==2){ /* matrix for the five point stencil */
62: n1 = (int) (sqrt((PetscReal)n) + 0.001);
63: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
64: for (i=0; i<n1; i++) {
65: for (j=0; j<n1; j++) {
66: I = j + n1*i;
67: if (i>0) {
68: J = I - n1;
69: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
70: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
71: }
72: if (i<n1-1) {
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 (j>0) {
78: J = I - 1;
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<n1-1) {
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: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
88: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
89: }
90: }
91: }
92: }
93: else { /* bs > 1 */
94: for (block=0; block<n/bs; block++){
95: /* diagonal blocks */
96: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
97: for (i=1+block*bs; i<bs-1+block*bs; i++) {
98: col[0] = i-1; col[1] = i; col[2] = i+1;
99: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
100: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
101: }
102: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
103: value[0]=-1.0; value[1]=4.0;
104: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
105: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
107: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
108: value[0]=4.0; value[1] = -1.0;
109: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
110: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
111: }
112: /* off-diagonal blocks */
113: value[0]=-1.0;
114: for (i=0; i<(n/bs-1)*bs; i++){
115: col[0]=i+bs;
116: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
117: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
118: col[0]=i; row=i+bs;
119: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
120: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
121: }
122: if (bs == 2){
123: /* insert a value to off-diag blocks */
124: row = 2; col[0] = 5; value[0] = 0.01;
125: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
126: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
127: row = 0; col[0] = 3; value[0] = 0.01;
128: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
129: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
130: }
131: }
132: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
134: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
135: MatView(A, PETSC_VIEWER_DRAW_WORLD);
136: MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
138: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
139: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
140: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n"); */
141: /* MatView(sA, PETSC_VIEWER_DRAW_WORLD); */
142: /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD); */
144: /* Vectors */
145: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
146: VecCreateSeq(PETSC_COMM_SELF,n,&x);
147: VecDuplicate(x,&b);
148: VecDuplicate(x,&y);
149: VecSetRandom(rand,x);
151: /* Test MatReordering() */
152: PetscMalloc(mbs*sizeof(int),&ip_ptr);
153: for (i=0; i<mbs; i++) ip_ptr[i] = i;
154: if(reorder){
155: i = ip_ptr[1]; ip_ptr[1] = ip_ptr[mbs-2]; ip_ptr[mbs-2] = i;
156: /* i = ip_ptr[0]; ip_ptr[0] = ip_ptr[mbs-1]; ip_ptr[mbs-1] = i; */
157: /* i = ip_ptr[2]; ip_ptr[2] = ip_ptr[mbs-3]; ip_ptr[mbs-3] = i; */
158: }
159: ISCreateGeneral(PETSC_COMM_SELF,mbs,ip_ptr,&perm);
160: ISSetPermutation(perm);
161:
162: /* Test MatCholeskyFactor(), MatICCFactor() */
163: norm1 = tol;
164: for (lf=-1; lf<10*bs; lf += bs){
165: if (lf==-1) { /* Cholesky factor */
166: factinfo.fill = 5.0;
167: MatCholeskyFactorSymbolic(sA,perm,&factinfo,&sC);
168: } else { /* incomplete Cholesky factor */
169: factinfo.fill = 5.0;
170: factinfo.levels = lf;
171: MatICCFactorSymbolic(sA,perm,&factinfo,&sC);
172: }
173: MatCholeskyFactorNumeric(sA,&sC);
174: /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */ /* view factored matrix */
175: /* MatView(sC, PETSC_VIEWER_STDOUT_WORLD); */
176:
177: MatMult(sA,x,b);
178: MatSolve(sC,b,y);
179: if (bs == 1) {
180: Vecs xx,bb;
181: VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
182: VecsDuplicate(xx,&bb);
183: MatSolves(sC,bb,xx);
184: VecsDestroy(xx);
185: VecsDestroy(bb);
186: }
187: MatDestroy(sC);
189: /* Check the error */
190: VecAXPY(&neg_one,x,y);
191: VecNorm(y,NORM_2,&norm2);
192: /* printf("lf: %d, error: %g\n", lf,norm2); */
193: if (10*norm1 < norm2 && lf-bs != -1){
194: PetscPrintf(PETSC_COMM_SELF,"lf=%d, %d, Norm of error=%g, %g\n",lf-bs,lf,norm1,norm2);
195: }
196: norm1 = norm2;
197: if (norm2 < tol && lf != -1) break;
198: }
200: ISDestroy(perm);
201: PetscFree(ip_ptr);
202: MatDestroy(A);
203: MatDestroy(sA);
204: VecDestroy(x);
205: VecDestroy(y);
206: VecDestroy(b);
207: PetscRandomDestroy(rand);
209: PetscFinalize();
210: return 0;
211: }