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;
24: PetscInitialize(&argc, &args, (char *)0, help);
25: MPI_Comm_size(PETSC_COMM_WORLD, &size);
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: 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;
48: value[1] = 2.0;
49: value[2] = -1.0;
50: for (i = 1; i < n - 1; i++) {
51: col[0] = i - 1;
52: col[1] = i;
53: col[2] = i + 1;
54: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
55: }
56: i = n - 1;
57: col[0] = 0;
58: col[1] = n - 2;
59: col[2] = n - 1;
61: value[0] = 0.1;
62: value[1] = -1;
63: value[2] = 2;
64: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
66: i = 0;
67: col[0] = 0;
68: col[1] = 1;
69: col[2] = n - 1;
71: value[0] = 2.0;
72: value[1] = -1.0;
73: value[2] = 0.1;
74: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
75: } else if (prob == 2) { /* matrix for the five point stencil */
76: n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
78: for (i = 0; i < n1; i++) {
79: for (j = 0; j < n1; j++) {
80: Ii = j + n1 * i;
81: if (i > 0) {
82: J = Ii - n1;
83: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
84: }
85: if (i < n1 - 1) {
86: J = Ii + n1;
87: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
88: }
89: if (j > 0) {
90: J = Ii - 1;
91: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
92: }
93: if (j < n1 - 1) {
94: J = Ii + 1;
95: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
96: }
97: MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES);
98: }
99: }
100: }
101: } else { /* bs > 1 */
102: for (block = 0; block < n / bs; block++) {
103: /* diagonal blocks */
104: value[0] = -1.0;
105: value[1] = 4.0;
106: value[2] = -1.0;
107: for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
108: col[0] = i - 1;
109: col[1] = i;
110: col[2] = i + 1;
111: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
112: }
113: i = bs - 1 + block * bs;
114: col[0] = bs - 2 + block * bs;
115: col[1] = bs - 1 + block * bs;
117: value[0] = -1.0;
118: value[1] = 4.0;
119: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
121: i = 0 + block * bs;
122: col[0] = 0 + block * bs;
123: col[1] = 1 + block * bs;
125: value[0] = 4.0;
126: value[1] = -1.0;
127: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
128: }
129: /* off-diagonal blocks */
130: value[0] = -1.0;
131: for (i = 0; i < (n / bs - 1) * bs; i++) {
132: col[0] = i + bs;
133: MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
134: col[0] = i;
135: row = i + bs;
136: MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
137: }
138: }
140: if (TestShift) {
141: /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
142: for (i = 0; i < bs; i++) {
143: row = i;
144: col[0] = i;
145: value[0] = 0.0;
146: MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
147: }
148: }
150: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
153: /* Test MatConvert */
154: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
155: MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA);
156: MatMultEqual(A, sA, 20, &equal);
159: /* Test MatGetOwnershipRange() */
160: MatGetOwnershipRange(A, &Ii, &J);
161: MatGetOwnershipRange(sA, &i, &j);
164: /* Vectors */
165: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
166: PetscRandomSetFromOptions(rdm);
167: VecCreateSeq(PETSC_COMM_SELF, n, &x);
168: VecDuplicate(x, &b);
169: VecDuplicate(x, &y);
170: VecSetRandom(x, rdm);
172: /* Test MatReordering() - not work on sbaij matrix */
173: if (reorder) {
174: MatGetOrdering(A, MATORDERINGRCM, &perm, &cperm);
175: } else {
176: MatGetOrdering(A, MATORDERINGNATURAL, &perm, &cperm);
177: }
178: ISDestroy(&cperm);
180: /* initialize factinfo */
181: MatFactorInfoInitialize(&factinfo);
182: if (TestShift == 1) {
183: factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
184: factinfo.shiftamount = 0.1;
185: } else if (TestShift == 2) {
186: factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
187: }
189: /* Test MatCholeskyFactor(), MatICCFactor() */
190: /*------------------------------------------*/
191: /* Test aij matrix A */
192: if (TestAIJ) {
193: if (displ) PetscPrintf(PETSC_COMM_WORLD, "AIJ: \n");
194: i = 0;
195: for (lvl = -1; lvl < 10; lvl++) {
196: if (lvl == -1) { /* Cholesky factor */
197: factinfo.fill = 5.0;
199: MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
200: MatCholeskyFactorSymbolic(sC, A, perm, &factinfo);
201: } else { /* incomplete Cholesky factor */
202: factinfo.fill = 5.0;
203: 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 residual */
215: VecAXPY(y, neg_one, x);
216: VecNorm(y, NORM_2, &norm2);
218: if (displ) PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2);
219: }
220: }
222: /* Test baij matrix A */
223: if (TestBAIJ) {
224: if (displ) PetscPrintf(PETSC_COMM_WORLD, "BAIJ: \n");
225: i = 0;
226: for (lvl = -1; lvl < 10; lvl++) {
227: if (lvl == -1) { /* Cholesky factor */
228: factinfo.fill = 5.0;
230: MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
231: MatCholeskyFactorSymbolic(sC, A, perm, &factinfo);
232: } else { /* incomplete Cholesky factor */
233: factinfo.fill = 5.0;
234: factinfo.levels = lvl;
236: MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC);
237: MatICCFactorSymbolic(sC, A, perm, &factinfo);
238: }
239: MatCholeskyFactorNumeric(sC, A, &factinfo);
241: MatMult(A, x, b);
242: MatSolve(sC, b, y);
243: MatDestroy(&sC);
245: /* Check the residual */
246: VecAXPY(y, neg_one, x);
247: VecNorm(y, NORM_2, &norm2);
248: if (displ) PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2);
249: }
250: }
252: /* Test sbaij matrix sA */
253: if (displ) PetscPrintf(PETSC_COMM_WORLD, "SBAIJ: \n");
254: i = 0;
255: for (lvl = -1; lvl < 10; lvl++) {
256: if (lvl == -1) { /* Cholesky factor */
257: factinfo.fill = 5.0;
259: MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC);
260: MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo);
261: } else { /* incomplete Cholesky factor */
262: factinfo.fill = 5.0;
263: factinfo.levels = lvl;
265: MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC);
266: MatICCFactorSymbolic(sC, sA, perm, &factinfo);
267: }
268: MatCholeskyFactorNumeric(sC, sA, &factinfo);
270: if (lvl == 0 && bs == 1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
271: /*
272: Mat B;
273: MatDuplicate(sA,MAT_COPY_VALUES,&B);
274: MatICCFactor(B,perm,&factinfo);
275: MatEqual(sC,B,&equal);
276: if (!equal) {
277: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
278: }
279: MatDestroy(&B);
280: */
281: }
283: MatMult(sA, x, b);
284: MatSolve(sC, b, y);
286: /* Test MatSolves() */
287: if (bs == 1) {
288: Vecs xx, bb;
289: VecsCreateSeq(PETSC_COMM_SELF, n, 4, &xx);
290: VecsDuplicate(xx, &bb);
291: MatSolves(sC, bb, xx);
292: VecsDestroy(xx);
293: VecsDestroy(bb);
294: }
295: MatDestroy(&sC);
297: /* Check the residual */
298: VecAXPY(y, neg_one, x);
299: VecNorm(y, NORM_2, &norm2);
300: if (displ) PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2);
301: }
303: ISDestroy(&perm);
304: MatDestroy(&A);
305: MatDestroy(&sA);
306: VecDestroy(&x);
307: VecDestroy(&y);
308: VecDestroy(&b);
309: PetscRandomDestroy(&rdm);
311: PetscFinalize();
312: return 0;
313: }
315: /*TEST
317: test:
318: args: -bs {{1 2 3 4 5 6 7 8}}
320: test:
321: suffix: 3
322: args: -testaij
323: output_file: output/ex76_1.out
325: TEST*/