Actual source code: ex76.c
2: static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Vec x,y,b;
9: Mat A; /* linear system matrix */
10: Mat sA,sC; /* symmetric part of the matrices */
11: PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,col[3],block, row,Ii,J,n1,lvl;
12: PetscMPIInt size;
13: PetscReal norm2;
14: PetscScalar neg_one = -1.0,four=4.0,value[3];
15: IS perm,cperm;
16: PetscRandom rdm;
17: PetscBool reorder = PETSC_FALSE,displ = PETSC_FALSE;
18: MatFactorInfo factinfo;
19: PetscBool equal;
20: PetscBool TestAIJ = PETSC_FALSE,TestBAIJ = PETSC_TRUE;
21: PetscInt TestShift=0;
23: PetscInitialize(&argc,&args,(char*)0,help);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
28: PetscOptionsGetBool(NULL,NULL,"-reorder",&reorder,NULL);
29: PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL);
31: PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);
33: n = mbs*bs;
34: if (TestAIJ) { /* A is in aij format */
35: MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A);
36: TestBAIJ = PETSC_FALSE;
37: } else { /* A is in baij format */
38: MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A);
39: TestAIJ = PETSC_FALSE;
40: }
42: /* Assemble matrix */
43: if (bs == 1) {
44: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
45: if (prob == 1) { /* tridiagonal matrix */
46: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
47: for (i=1; i<n-1; i++) {
48: col[0] = i-1; col[1] = i; col[2] = i+1;
49: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
50: }
51: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
53: value[0]= 0.1; value[1]=-1; value[2]=2;
54: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
56: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
58: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
59: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
60: } else if (prob ==2) { /* matrix for the five point stencil */
61: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
63: for (i=0; i<n1; i++) {
64: for (j=0; j<n1; j++) {
65: Ii = j + n1*i;
66: if (i>0) {
67: J = Ii - n1;
68: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
69: }
70: if (i<n1-1) {
71: J = Ii + n1;
72: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
73: }
74: if (j>0) {
75: J = Ii - 1;
76: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
77: }
78: if (j<n1-1) {
79: J = Ii + 1;
80: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
81: }
82: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
83: }
84: }
85: }
86: } else { /* bs > 1 */
87: for (block=0; block<n/bs; block++) {
88: /* diagonal blocks */
89: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
90: for (i=1+block*bs; i<bs-1+block*bs; i++) {
91: col[0] = i-1; col[1] = i; col[2] = i+1;
92: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
93: }
94: 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;
101: value[0]=4.0; value[1] = -1.0;
102: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
103: }
104: /* off-diagonal blocks */
105: value[0]=-1.0;
106: for (i=0; i<(n/bs-1)*bs; i++) {
107: col[0]=i+bs;
108: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
109: col[0]=i; row=i+bs;
110: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
111: }
112: }
114: if (TestShift) {
115: /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
116: for (i=0; i<bs; i++) {
117: row = i; col[0] = i; value[0] = 0.0;
118: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
119: }
120: }
122: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
125: /* Test MatConvert */
126: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
127: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
128: MatMultEqual(A,sA,20,&equal);
131: /* Test MatGetOwnershipRange() */
132: MatGetOwnershipRange(A,&Ii,&J);
133: MatGetOwnershipRange(sA,&i,&j);
136: /* Vectors */
137: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
138: PetscRandomSetFromOptions(rdm);
139: VecCreateSeq(PETSC_COMM_SELF,n,&x);
140: VecDuplicate(x,&b);
141: VecDuplicate(x,&y);
142: VecSetRandom(x,rdm);
144: /* Test MatReordering() - not work on sbaij matrix */
145: if (reorder) {
146: MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);
147: } else {
148: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);
149: }
150: ISDestroy(&cperm);
152: /* initialize factinfo */
153: MatFactorInfoInitialize(&factinfo);
154: if (TestShift == 1) {
155: factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
156: factinfo.shiftamount = 0.1;
157: } else if (TestShift == 2) {
158: factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
159: }
161: /* Test MatCholeskyFactor(), MatICCFactor() */
162: /*------------------------------------------*/
163: /* Test aij matrix A */
164: if (TestAIJ) {
165: if (displ) {
166: PetscPrintf(PETSC_COMM_WORLD,"AIJ: \n");
167: }
168: i = 0;
169: for (lvl=-1; lvl<10; lvl++) {
170: if (lvl==-1) { /* Cholesky factor */
171: factinfo.fill = 5.0;
173: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
174: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
175: } else { /* incomplete Cholesky factor */
176: factinfo.fill = 5.0;
177: factinfo.levels = lvl;
179: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
180: MatICCFactorSymbolic(sC,A,perm,&factinfo);
181: }
182: MatCholeskyFactorNumeric(sC,A,&factinfo);
184: MatMult(A,x,b);
185: MatSolve(sC,b,y);
186: MatDestroy(&sC);
188: /* Check the residual */
189: VecAXPY(y,neg_one,x);
190: VecNorm(y,NORM_2,&norm2);
192: if (displ) {
193: PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
194: }
195: }
196: }
198: /* Test baij matrix A */
199: if (TestBAIJ) {
200: if (displ) {
201: PetscPrintf(PETSC_COMM_WORLD,"BAIJ: \n");
202: }
203: i = 0;
204: for (lvl=-1; lvl<10; lvl++) {
205: if (lvl==-1) { /* Cholesky factor */
206: factinfo.fill = 5.0;
208: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
209: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
210: } else { /* incomplete Cholesky factor */
211: factinfo.fill = 5.0;
212: factinfo.levels = lvl;
214: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
215: MatICCFactorSymbolic(sC,A,perm,&factinfo);
216: }
217: MatCholeskyFactorNumeric(sC,A,&factinfo);
219: MatMult(A,x,b);
220: MatSolve(sC,b,y);
221: MatDestroy(&sC);
223: /* Check the residual */
224: VecAXPY(y,neg_one,x);
225: VecNorm(y,NORM_2,&norm2);
226: if (displ) {
227: PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
228: }
229: }
230: }
232: /* Test sbaij matrix sA */
233: if (displ) {
234: PetscPrintf(PETSC_COMM_WORLD,"SBAIJ: \n");
235: }
236: i = 0;
237: for (lvl=-1; lvl<10; lvl++) {
238: if (lvl==-1) { /* Cholesky factor */
239: factinfo.fill = 5.0;
241: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
242: MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
243: } else { /* incomplete Cholesky factor */
244: factinfo.fill = 5.0;
245: factinfo.levels = lvl;
247: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
248: MatICCFactorSymbolic(sC,sA,perm,&factinfo);
249: }
250: MatCholeskyFactorNumeric(sC,sA,&factinfo);
252: if (lvl==0 && bs==1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
253: /*
254: Mat B;
255: MatDuplicate(sA,MAT_COPY_VALUES,&B);
256: MatICCFactor(B,perm,&factinfo);
257: MatEqual(sC,B,&equal);
258: if (!equal) {
259: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
260: }
261: MatDestroy(&B);
262: */
263: }
265: MatMult(sA,x,b);
266: MatSolve(sC,b,y);
268: /* Test MatSolves() */
269: if (bs == 1) {
270: Vecs xx,bb;
271: VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
272: VecsDuplicate(xx,&bb);
273: MatSolves(sC,bb,xx);
274: VecsDestroy(xx);
275: VecsDestroy(bb);
276: }
277: MatDestroy(&sC);
279: /* Check the residual */
280: VecAXPY(y,neg_one,x);
281: VecNorm(y,NORM_2,&norm2);
282: if (displ) {
283: PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
284: }
285: }
287: ISDestroy(&perm);
288: MatDestroy(&A);
289: MatDestroy(&sA);
290: VecDestroy(&x);
291: VecDestroy(&y);
292: VecDestroy(&b);
293: PetscRandomDestroy(&rdm);
295: PetscFinalize();
296: return 0;
297: }
299: /*TEST
301: test:
302: args: -bs {{1 2 3 4 5 6 7 8}}
304: test:
305: suffix: 3
306: args: -testaij
307: output_file: output/ex76_1.out
309: TEST*/