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: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
25: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-reorder", &reorder, NULL));
30: PetscCall(PetscOptionsGetBool(NULL, NULL, "-testaij", &TestAIJ, NULL));
31: PetscCall(PetscOptionsGetInt(NULL, NULL, "-testShift", &TestShift, NULL));
32: PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));
34: n = mbs * bs;
35: if (TestAIJ) { /* A is in aij format */
36: PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, n, n, nz, NULL, &A));
37: TestBAIJ = PETSC_FALSE;
38: } else { /* A is in baij format */
39: PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A));
40: TestAIJ = PETSC_FALSE;
41: }
43: /* Assemble matrix */
44: if (bs == 1) {
45: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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);
77: PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
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: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
84: }
85: if (i < n1 - 1) {
86: J = Ii + n1;
87: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
88: }
89: if (j > 0) {
90: J = Ii - 1;
91: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
92: }
93: if (j < n1 - 1) {
94: J = Ii + 1;
95: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
96: }
97: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
134: col[0] = i;
135: row = i + bs;
136: PetscCall(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: PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
147: }
148: }
150: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
151: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
153: /* Test MatConvert */
154: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
155: PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));
156: PetscCall(MatMultEqual(A, sA, 20, &equal));
157: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_USER, "A != sA");
159: /* Test MatGetOwnershipRange() */
160: PetscCall(MatGetOwnershipRange(A, &Ii, &J));
161: PetscCall(MatGetOwnershipRange(sA, &i, &j));
162: PetscCheck(i == Ii && j == J, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetOwnershipRange() in MatSBAIJ format");
164: /* Vectors */
165: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
166: PetscCall(PetscRandomSetFromOptions(rdm));
167: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
168: PetscCall(VecDuplicate(x, &b));
169: PetscCall(VecDuplicate(x, &y));
170: PetscCall(VecSetRandom(x, rdm));
172: /* Test MatReordering() - not work on sbaij matrix */
173: if (reorder) {
174: PetscCall(MatGetOrdering(A, MATORDERINGRCM, &perm, &cperm));
175: } else {
176: PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &cperm));
177: }
178: PetscCall(ISDestroy(&cperm));
180: /* initialize factinfo */
181: PetscCall(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) PetscCall(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: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
200: PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
201: } else { /* incomplete Cholesky factor */
202: factinfo.fill = 5.0;
203: factinfo.levels = lvl;
205: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
206: PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
207: }
208: PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));
210: PetscCall(MatMult(A, x, b));
211: PetscCall(MatSolve(sC, b, y));
212: PetscCall(MatDestroy(&sC));
214: /* Check the residual */
215: PetscCall(VecAXPY(y, neg_one, x));
216: PetscCall(VecNorm(y, NORM_2, &norm2));
218: if (displ) PetscCall(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) PetscCall(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: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
231: PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
232: } else { /* incomplete Cholesky factor */
233: factinfo.fill = 5.0;
234: factinfo.levels = lvl;
236: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
237: PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
238: }
239: PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));
241: PetscCall(MatMult(A, x, b));
242: PetscCall(MatSolve(sC, b, y));
243: PetscCall(MatDestroy(&sC));
245: /* Check the residual */
246: PetscCall(VecAXPY(y, neg_one, x));
247: PetscCall(VecNorm(y, NORM_2, &norm2));
248: if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
249: }
250: }
252: /* Test sbaij matrix sA */
253: if (displ) PetscCall(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: PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
260: PetscCall(MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo));
261: } else { /* incomplete Cholesky factor */
262: factinfo.fill = 5.0;
263: factinfo.levels = lvl;
265: PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
266: PetscCall(MatICCFactorSymbolic(sC, sA, perm, &factinfo));
267: }
268: PetscCall(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: PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&B));
274: PetscCall(MatICCFactor(B,perm,&factinfo));
275: PetscCall(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: PetscCall(MatDestroy(&B));
280: */
281: }
283: PetscCall(MatMult(sA, x, b));
284: PetscCall(MatSolve(sC, b, y));
286: /* Test MatSolves() */
287: if (bs == 1) {
288: Vecs xx, bb;
289: PetscCall(VecsCreateSeq(PETSC_COMM_SELF, n, 4, &xx));
290: PetscCall(VecsDuplicate(xx, &bb));
291: PetscCall(MatSolves(sC, bb, xx));
292: PetscCall(VecsDestroy(xx));
293: PetscCall(VecsDestroy(bb));
294: }
295: PetscCall(MatDestroy(&sC));
297: /* Check the residual */
298: PetscCall(VecAXPY(y, neg_one, x));
299: PetscCall(VecNorm(y, NORM_2, &norm2));
300: if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
301: }
303: PetscCall(ISDestroy(&perm));
304: PetscCall(MatDestroy(&A));
305: PetscCall(MatDestroy(&sA));
306: PetscCall(VecDestroy(&x));
307: PetscCall(VecDestroy(&y));
308: PetscCall(VecDestroy(&b));
309: PetscCall(PetscRandomDestroy(&rdm));
311: PetscCall(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*/