Actual source code: ex76.c
petsc-3.3-p6 2013-02-11
2: static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n";
4: #include <petscmat.h>
8: int main(int argc,char **args)
9: {
10: Vec x,y,b;
11: Mat A; /* linear system matrix */
12: Mat sA,sC; /* symmetric part of the matrices */
13: PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,col[3],block, row,Ii,J,n1,lvl;
15: PetscMPIInt size;
16: PetscReal norm2,tol=1.e-10,err[10];
17: PetscScalar neg_one = -1.0,four=4.0,value[3];
18: IS perm,cperm;
19: PetscRandom rdm;
20: PetscInt reorder=0,displ=0;
21: MatFactorInfo factinfo;
22: PetscBool equal;
23: PetscBool TestAIJ=PETSC_FALSE,TestBAIJ=PETSC_TRUE;
24: PetscInt TestShift=0;
26: PetscInitialize(&argc,&args,(char *)0,help);
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
28: if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
29: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
30: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
31: PetscOptionsGetInt(PETSC_NULL,"-reorder",&reorder,PETSC_NULL);
32: PetscOptionsGetBool(PETSC_NULL,"-testaij",&TestAIJ,PETSC_NULL);
33: PetscOptionsGetInt(PETSC_NULL,"-testShift",&TestShift,PETSC_NULL);
34: PetscOptionsGetInt(PETSC_NULL,"-displ",&displ,PETSC_NULL);
36: n = mbs*bs;
37: if (TestAIJ){ /* A is in aij format */
38: MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,PETSC_NULL,&A);
39: TestBAIJ = PETSC_FALSE;
40: } else { /* A is in baij format */
41: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL,&A);
42: TestAIJ = PETSC_FALSE;
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: }
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);
58: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
59: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
60: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
61: } else if (prob ==2){ /* matrix for the five point stencil */
62: n1 = (int) (PetscSqrtReal((PetscReal)n) + 0.001);
63: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,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: Ii = j + n1*i;
67: if (i>0) {
68: J = Ii - n1;
69: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
70: }
71: if (i<n1-1) {
72: J = Ii + n1;
73: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
74: }
75: if (j>0) {
76: J = Ii - 1;
77: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
78: }
79: if (j<n1-1) {
80: J = Ii + 1;
81: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
82: }
83: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
84: }
85: }
86: }
87: } else { /* bs > 1 */
88: for (block=0; block<n/bs; block++){
89: /* diagonal blocks */
90: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
91: for (i=1+block*bs; i<bs-1+block*bs; i++) {
92: col[0] = i-1; col[1] = i; col[2] = i+1;
93: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
94: }
95: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
96: value[0]=-1.0; value[1]=4.0;
97: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
99: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
100: value[0]=4.0; value[1] = -1.0;
101: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
102: }
103: /* off-diagonal blocks */
104: value[0]=-1.0;
105: for (i=0; i<(n/bs-1)*bs; i++){
106: col[0]=i+bs;
107: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
108: col[0]=i; row=i+bs;
109: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
110: }
111: }
113: if (TestShift){
114: /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
115: for (i=0; i<bs; i++){
116: row = i; col[0] = i; value[0] = 0.0;
117: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
118: }
119: }
121: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
122: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
124: /* Test MatConvert */
125: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
126: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
127: MatMultEqual(A,sA,20,&equal);
128: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA");
130: /* Test MatGetOwnershipRange() */
131: MatGetOwnershipRange(A,&Ii,&J);
132: MatGetOwnershipRange(sA,&i,&j);
133: if (i-Ii || j-J){
134: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
135: }
137: /* Vectors */
138: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
139: PetscRandomSetFromOptions(rdm);
140: VecCreateSeq(PETSC_COMM_SELF,n,&x);
141: VecDuplicate(x,&b);
142: VecDuplicate(x,&y);
143: VecSetRandom(x,rdm);
145: /* Test MatReordering() - not work on sbaij matrix */
146: if (reorder){
147: MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);
148: } else {
149: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);
150: }
151: ISDestroy(&cperm);
153: /* initialize factinfo */
154: MatFactorInfoInitialize(&factinfo);
155: if (TestShift == 1){
156: factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
157: factinfo.shiftamount = 0.1;
158: } else if (TestShift == 2){
159: factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
160: }
161:
162: /* Test MatCholeskyFactor(), MatICCFactor() */
163: /*------------------------------------------*/
164: /* Test aij matrix A */
165: if (TestAIJ){
166: if (displ>0){PetscPrintf(PETSC_COMM_SELF,"AIJ: \n");}
167: i = 0;
168: for (lvl=-1; lvl<10; lvl++){
169: if (lvl==-1) { /* Cholesky factor */
170: factinfo.fill = 5.0;
171: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
172: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
173: } else { /* incomplete Cholesky factor */
174: factinfo.fill = 5.0;
175: factinfo.levels = lvl;
176: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
177: MatICCFactorSymbolic(sC,A,perm,&factinfo);
178: }
179: MatCholeskyFactorNumeric(sC,A,&factinfo);
181: MatMult(A,x,b);
182: MatSolve(sC,b,y);
183: MatDestroy(&sC);
185: /* Check the error */
186: VecAXPY(y,neg_one,x);
187: VecNorm(y,NORM_2,&norm2);
188: if (displ>0){PetscPrintf(PETSC_COMM_SELF," lvl: %d, error: %G\n", lvl,norm2);}
189: err[i++] = norm2;
190: }
191: }
192:
193: /* Test baij matrix A */
194: if (TestBAIJ){
195: if (displ>0){PetscPrintf(PETSC_COMM_SELF,"BAIJ: \n");}
196: i = 0;
197: for (lvl=-1; lvl<10; lvl++){
198: if (lvl==-1) { /* Cholesky factor */
199: factinfo.fill = 5.0;
200: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
201: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
202: } else { /* incomplete Cholesky factor */
203: factinfo.fill = 5.0;
204: factinfo.levels = lvl;
205: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
206: MatICCFactorSymbolic(sC,A,perm,&factinfo);
207: }
208: MatCholeskyFactorNumeric(sC,A,&factinfo);
210: MatMult(A,x,b);
211: MatSolve(sC,b,y);
212: MatDestroy(&sC);
214: /* Check the error */
215: VecAXPY(y,neg_one,x);
216: VecNorm(y,NORM_2,&norm2);
217: if (displ>0){PetscPrintf(PETSC_COMM_SELF," lvl: %d, error: %G\n", lvl,norm2);}
218: err[i++] = norm2;
219: }
220: }
222: /* Test sbaij matrix sA */
223: if (displ>0){PetscPrintf(PETSC_COMM_SELF,"SBAIJ: \n");}
224: i = 0;
225: for (lvl=-1; lvl<10; lvl++){
226: if (lvl==-1) { /* Cholesky factor */
227: factinfo.fill = 5.0;
228: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
229: MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
230: } else { /* incomplete Cholesky factor */
231: factinfo.fill = 5.0;
232: factinfo.levels = lvl;
233: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
234: MatICCFactorSymbolic(sC,sA,perm,&factinfo);
235: }
236: MatCholeskyFactorNumeric(sC,sA,&factinfo);
238: if (lvl==0 && bs==1){ /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
239: /*
240: Mat B;
241: MatDuplicate(sA,MAT_COPY_VALUES,&B);
242: MatICCFactor(B,perm,&factinfo);
243: MatEqual(sC,B,&equal);
244: if (!equal){
245: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
246: }
247: MatDestroy(&B);
248: */
249: }
252: MatMult(sA,x,b);
253: MatSolve(sC,b,y);
255: /* Test MatSolves() */
256: if (bs == 1) {
257: Vecs xx,bb;
258: VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
259: VecsDuplicate(xx,&bb);
260: MatSolves(sC,bb,xx);
261: VecsDestroy(xx);
262: VecsDestroy(bb);
263: }
264: MatDestroy(&sC);
266: /* Check the error */
267: VecAXPY(y,neg_one,x);
268: VecNorm(y,NORM_2,&norm2);
269: if (displ>0){PetscPrintf(PETSC_COMM_SELF," lvl: %d, error: %G\n", lvl,norm2); }
270: err[i] -= norm2;
271: if (err[i] > tol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER," level: %d, err: %G\n", lvl,err[i]);
272: }
274: ISDestroy(&perm);
275: MatDestroy(&A);
276: MatDestroy(&sA);
277: VecDestroy(&x);
278: VecDestroy(&y);
279: VecDestroy(&b);
280: PetscRandomDestroy(&rdm);
282: PetscFinalize();
283: return 0;
284: }