Actual source code: ex76.c
petsc-3.7.3 2016-08-01
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(NULL,NULL,"-bs",&bs,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-reorder",&reorder,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL);
34: PetscOptionsGetInt(NULL,NULL,"-displ",&displ,NULL);
36: n = mbs*bs;
37: if (TestAIJ) { /* A is in aij format */
38: MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A);
39: TestBAIJ = PETSC_FALSE;
40: } else { /* A is in baij format */
41: ierr =MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A);
42: TestAIJ = PETSC_FALSE;
43: }
45: /* Assemble matrix */
46: if (bs == 1) {
47: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,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;
56: value[0]= 0.1; value[1]=-1; value[2]=2;
57: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
59: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
61: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
62: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
63: } else if (prob ==2) { /* matrix for the five point stencil */
64: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
65: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
66: for (i=0; i<n1; i++) {
67: for (j=0; j<n1; j++) {
68: Ii = j + n1*i;
69: if (i>0) {
70: J = Ii - n1;
71: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
72: }
73: if (i<n1-1) {
74: J = Ii + n1;
75: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
76: }
77: if (j>0) {
78: J = Ii - 1;
79: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
80: }
81: if (j<n1-1) {
82: J = Ii + 1;
83: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
84: }
85: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
86: }
87: }
88: }
89: } else { /* bs > 1 */
90: for (block=0; block<n/bs; block++) {
91: /* diagonal blocks */
92: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
93: for (i=1+block*bs; i<bs-1+block*bs; i++) {
94: col[0] = i-1; col[1] = i; col[2] = i+1;
95: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
96: }
97: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
99: value[0]=-1.0; value[1]=4.0;
100: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
102: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
104: value[0]=4.0; value[1] = -1.0;
105: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
106: }
107: /* off-diagonal blocks */
108: value[0]=-1.0;
109: for (i=0; i<(n/bs-1)*bs; i++) {
110: col[0]=i+bs;
111: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
112: col[0]=i; row=i+bs;
113: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
114: }
115: }
117: if (TestShift) {
118: /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
119: for (i=0; i<bs; i++) {
120: row = i; col[0] = i; value[0] = 0.0;
121: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
122: }
123: }
125: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
128: /* Test MatConvert */
129: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
130: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
131: MatMultEqual(A,sA,20,&equal);
132: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA");
134: /* Test MatGetOwnershipRange() */
135: MatGetOwnershipRange(A,&Ii,&J);
136: MatGetOwnershipRange(sA,&i,&j);
137: if (i-Ii || j-J) {
138: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
139: }
141: /* Vectors */
142: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
143: PetscRandomSetFromOptions(rdm);
144: VecCreateSeq(PETSC_COMM_SELF,n,&x);
145: VecDuplicate(x,&b);
146: VecDuplicate(x,&y);
147: VecSetRandom(x,rdm);
149: /* Test MatReordering() - not work on sbaij matrix */
150: if (reorder) {
151: MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);
152: } else {
153: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);
154: }
155: ISDestroy(&cperm);
157: /* initialize factinfo */
158: MatFactorInfoInitialize(&factinfo);
159: if (TestShift == 1) {
160: factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
161: factinfo.shiftamount = 0.1;
162: } else if (TestShift == 2) {
163: factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
164: }
166: /* Test MatCholeskyFactor(), MatICCFactor() */
167: /*------------------------------------------*/
168: /* Test aij matrix A */
169: if (TestAIJ) {
170: if (displ>0) {
171: PetscPrintf(PETSC_COMM_SELF,"AIJ: \n");
172: }
173: i = 0;
174: for (lvl=-1; lvl<10; lvl++) {
175: if (lvl==-1) { /* Cholesky factor */
176: factinfo.fill = 5.0;
178: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
179: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
180: } else { /* incomplete Cholesky factor */
181: factinfo.fill = 5.0;
182: factinfo.levels = lvl;
184: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
185: MatICCFactorSymbolic(sC,A,perm,&factinfo);
186: }
187: MatCholeskyFactorNumeric(sC,A,&factinfo);
189: MatMult(A,x,b);
190: MatSolve(sC,b,y);
191: MatDestroy(&sC);
193: /* Check the error */
194: VecAXPY(y,neg_one,x);
195: VecNorm(y,NORM_2,&norm2);
197: if (displ>0) {
198: PetscPrintf(PETSC_COMM_SELF," lvl: %D, error: %g\n", lvl,(double)norm2);
199: }
200: err[i++] = norm2;
201: }
202: }
204: /* Test baij matrix A */
205: if (TestBAIJ) {
206: if (displ>0) {
207: PetscPrintf(PETSC_COMM_SELF,"BAIJ: \n");
208: }
209: i = 0;
210: for (lvl=-1; lvl<10; lvl++) {
211: if (lvl==-1) { /* Cholesky factor */
212: factinfo.fill = 5.0;
214: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
215: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
216: } else { /* incomplete Cholesky factor */
217: factinfo.fill = 5.0;
218: factinfo.levels = lvl;
220: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
221: MatICCFactorSymbolic(sC,A,perm,&factinfo);
222: }
223: MatCholeskyFactorNumeric(sC,A,&factinfo);
225: MatMult(A,x,b);
226: MatSolve(sC,b,y);
227: MatDestroy(&sC);
229: /* Check the error */
230: VecAXPY(y,neg_one,x);
231: VecNorm(y,NORM_2,&norm2);
232: if (displ>0) {
233: PetscPrintf(PETSC_COMM_SELF," lvl: %D, error: %g\n", lvl,(double)norm2);
234: }
235: err[i++] = norm2;
236: }
237: }
239: /* Test sbaij matrix sA */
240: if (displ>0) {
241: PetscPrintf(PETSC_COMM_SELF,"SBAIJ: \n");
242: }
243: i = 0;
244: for (lvl=-1; lvl<10; lvl++) {
245: if (lvl==-1) { /* Cholesky factor */
246: factinfo.fill = 5.0;
248: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
249: MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
250: } else { /* incomplete Cholesky factor */
251: factinfo.fill = 5.0;
252: factinfo.levels = lvl;
254: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
255: MatICCFactorSymbolic(sC,sA,perm,&factinfo);
256: }
257: MatCholeskyFactorNumeric(sC,sA,&factinfo);
259: if (lvl==0 && bs==1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
260: /*
261: Mat B;
262: MatDuplicate(sA,MAT_COPY_VALUES,&B);
263: MatICCFactor(B,perm,&factinfo);
264: MatEqual(sC,B,&equal);
265: if (!equal) {
266: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
267: }
268: MatDestroy(&B);
269: */
270: }
273: MatMult(sA,x,b);
274: MatSolve(sC,b,y);
276: /* Test MatSolves() */
277: if (bs == 1) {
278: Vecs xx,bb;
279: VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
280: VecsDuplicate(xx,&bb);
281: MatSolves(sC,bb,xx);
282: VecsDestroy(xx);
283: VecsDestroy(bb);
284: }
285: MatDestroy(&sC);
287: /* Check the error */
288: VecAXPY(y,neg_one,x);
289: VecNorm(y,NORM_2,&norm2);
290: if (displ>0) {
291: PetscPrintf(PETSC_COMM_SELF," lvl: %D, error: %g\n", lvl,(double)norm2);
292: }
293: err[i] -= norm2;
294: if (err[i] > tol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER," level: %d, err: %g\n", lvl,(double)err[i]);
295: }
297: ISDestroy(&perm);
298: MatDestroy(&A);
299: MatDestroy(&sA);
300: VecDestroy(&x);
301: VecDestroy(&y);
302: VecDestroy(&b);
303: PetscRandomDestroy(&rdm);
305: PetscFinalize();
306: return 0;
307: }