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