Actual source code: ex76.c
petsc-3.8.4 2018-03-24
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;
13: PetscMPIInt size;
14: PetscReal norm2,tol=1.e-10,err[10];
15: PetscScalar neg_one = -1.0,four=4.0,value[3];
16: IS perm,cperm;
17: PetscRandom rdm;
18: PetscBool reorder = PETSC_FALSE,displ = PETSC_FALSE;
19: MatFactorInfo factinfo;
20: PetscBool equal;
21: PetscBool TestAIJ = PETSC_FALSE,TestBAIJ = PETSC_TRUE;
22: PetscInt TestShift=0;
24: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
27: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
29: PetscOptionsGetBool(NULL,NULL,"-reorder",&reorder,NULL);
30: PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);
34: n = mbs*bs;
35: if (TestAIJ) { /* A is in aij format */
36: MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A);
37: TestBAIJ = PETSC_FALSE;
38: } else { /* A is in baij format */
39: ierr =MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A);
40: TestAIJ = PETSC_FALSE;
41: }
43: /* Assemble matrix */
44: if (bs == 1) {
45: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
46: if (prob == 1) { /* tridiagonal matrix */
47: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
48: for (i=1; i<n-1; i++) {
49: col[0] = i-1; col[1] = i; col[2] = i+1;
50: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
51: }
52: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
54: value[0]= 0.1; value[1]=-1; value[2]=2;
55: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
57: 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 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
63: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!");
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;
97: value[0]=-1.0; value[1]=4.0;
98: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
100: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
102: value[0]=4.0; value[1] = -1.0;
103: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
104: }
105: /* off-diagonal blocks */
106: value[0]=-1.0;
107: for (i=0; i<(n/bs-1)*bs; i++) {
108: col[0]=i+bs;
109: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
110: col[0]=i; row=i+bs;
111: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
112: }
113: }
115: if (TestShift) {
116: /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
117: for (i=0; i<bs; i++) {
118: row = i; col[0] = i; value[0] = 0.0;
119: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
120: }
121: }
123: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
126: /* Test MatConvert */
127: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
128: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
129: MatMultEqual(A,sA,20,&equal);
130: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA");
132: /* Test MatGetOwnershipRange() */
133: MatGetOwnershipRange(A,&Ii,&J);
134: MatGetOwnershipRange(sA,&i,&j);
135: if (i-Ii || j-J) {
136: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
137: }
139: /* Vectors */
140: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
141: PetscRandomSetFromOptions(rdm);
142: VecCreateSeq(PETSC_COMM_SELF,n,&x);
143: VecDuplicate(x,&b);
144: VecDuplicate(x,&y);
145: VecSetRandom(x,rdm);
147: /* Test MatReordering() - not work on sbaij matrix */
148: if (reorder) {
149: MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);
150: } else {
151: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);
152: }
153: ISDestroy(&cperm);
155: /* initialize factinfo */
156: MatFactorInfoInitialize(&factinfo);
157: if (TestShift == 1) {
158: factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
159: factinfo.shiftamount = 0.1;
160: } else if (TestShift == 2) {
161: factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
162: }
164: /* Test MatCholeskyFactor(), MatICCFactor() */
165: /*------------------------------------------*/
166: /* Test aij matrix A */
167: if (TestAIJ) {
168: if (displ) {
169: PetscPrintf(PETSC_COMM_WORLD,"AIJ: \n");
170: }
171: i = 0;
172: for (lvl=-1; lvl<10; lvl++) {
173: if (lvl==-1) { /* Cholesky factor */
174: factinfo.fill = 5.0;
176: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
177: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
178: } else { /* incomplete Cholesky factor */
179: factinfo.fill = 5.0;
180: factinfo.levels = lvl;
182: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
183: MatICCFactorSymbolic(sC,A,perm,&factinfo);
184: }
185: MatCholeskyFactorNumeric(sC,A,&factinfo);
187: MatMult(A,x,b);
188: MatSolve(sC,b,y);
189: MatDestroy(&sC);
191: /* Check the error */
192: VecAXPY(y,neg_one,x);
193: VecNorm(y,NORM_2,&norm2);
195: if (displ) {
196: PetscPrintf(PETSC_COMM_WORLD," lvl: %D, error: %g\n", lvl,(double)norm2);
197: }
198: err[i++] = norm2;
199: }
200: }
202: /* Test baij matrix A */
203: if (TestBAIJ) {
204: if (displ) {
205: PetscPrintf(PETSC_COMM_WORLD,"BAIJ: \n");
206: }
207: i = 0;
208: for (lvl=-1; lvl<10; lvl++) {
209: if (lvl==-1) { /* Cholesky factor */
210: factinfo.fill = 5.0;
212: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
213: MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
214: } else { /* incomplete Cholesky factor */
215: factinfo.fill = 5.0;
216: factinfo.levels = lvl;
218: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
219: MatICCFactorSymbolic(sC,A,perm,&factinfo);
220: }
221: MatCholeskyFactorNumeric(sC,A,&factinfo);
223: MatMult(A,x,b);
224: MatSolve(sC,b,y);
225: MatDestroy(&sC);
227: /* Check the error */
228: VecAXPY(y,neg_one,x);
229: VecNorm(y,NORM_2,&norm2);
230: if (displ) {
231: PetscPrintf(PETSC_COMM_WORLD," lvl: %D, error: %g\n", lvl,(double)norm2);
232: }
233: err[i++] = norm2;
234: }
235: }
237: /* Test sbaij matrix sA */
238: if (displ) {
239: PetscPrintf(PETSC_COMM_WORLD,"SBAIJ: \n");
240: }
241: i = 0;
242: for (lvl=-1; lvl<10; lvl++) {
243: if (lvl==-1) { /* Cholesky factor */
244: factinfo.fill = 5.0;
246: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
247: MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
248: } else { /* incomplete Cholesky factor */
249: factinfo.fill = 5.0;
250: factinfo.levels = lvl;
252: MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
253: MatICCFactorSymbolic(sC,sA,perm,&factinfo);
254: }
255: MatCholeskyFactorNumeric(sC,sA,&factinfo);
257: if (lvl==0 && bs==1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
258: /*
259: Mat B;
260: MatDuplicate(sA,MAT_COPY_VALUES,&B);
261: MatICCFactor(B,perm,&factinfo);
262: MatEqual(sC,B,&equal);
263: if (!equal) {
264: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
265: }
266: MatDestroy(&B);
267: */
268: }
271: MatMult(sA,x,b);
272: MatSolve(sC,b,y);
274: /* Test MatSolves() */
275: if (bs == 1) {
276: Vecs xx,bb;
277: VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
278: VecsDuplicate(xx,&bb);
279: MatSolves(sC,bb,xx);
280: VecsDestroy(xx);
281: VecsDestroy(bb);
282: }
283: MatDestroy(&sC);
285: /* Check the error */
286: VecAXPY(y,neg_one,x);
287: VecNorm(y,NORM_2,&norm2);
288: if (displ) {
289: PetscPrintf(PETSC_COMM_WORLD," lvl: %D, error: %g\n", lvl,(double)norm2);
290: }
291: err[i] -= norm2;
292: if (err[i] > tol) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_USER," level: %d, err: %g\n", lvl,(double)err[i]);
293: }
295: ISDestroy(&perm);
296: MatDestroy(&A);
297: MatDestroy(&sA);
298: VecDestroy(&x);
299: VecDestroy(&y);
300: VecDestroy(&b);
301: PetscRandomDestroy(&rdm);
303: PetscFinalize();
304: return ierr;
305: }