Actual source code: ex49.c
2: static char help[] = "Tests SeqSBAIJ factorizations for different block sizes\n\n";
4: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Vec x, b, u;
9: Mat A, A2;
10: KSP ksp;
11: PetscRandom rctx;
12: PetscReal norm;
13: PetscInt i, j, k, l, n = 27, its, bs = 2, Ii, J;
14: PetscBool test_hermitian = PETSC_FALSE, convert = PETSC_FALSE;
15: PetscScalar v;
18: PetscInitialize(&argc, &args, (char *)0, help);
19: PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
20: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
21: PetscOptionsGetBool(NULL, NULL, "-herm", &test_hermitian, NULL);
22: PetscOptionsGetBool(NULL, NULL, "-conv", &convert, NULL);
24: MatCreate(PETSC_COMM_SELF, &A);
25: MatSetSizes(A, n * bs, n * bs, PETSC_DETERMINE, PETSC_DETERMINE);
26: MatSetBlockSize(A, bs);
27: MatSetType(A, MATSEQSBAIJ);
28: MatSetFromOptions(A);
29: MatSeqSBAIJSetPreallocation(A, bs, n, NULL);
30: MatSeqBAIJSetPreallocation(A, bs, n, NULL);
31: MatSeqAIJSetPreallocation(A, n * bs, NULL);
32: MatMPIAIJSetPreallocation(A, n * bs, NULL, n * bs, NULL);
34: PetscRandomCreate(PETSC_COMM_SELF, &rctx);
35: for (i = 0; i < n; i++) {
36: for (j = i; j < n; j++) {
37: PetscRandomGetValue(rctx, &v);
38: if (PetscRealPart(v) < .1 || i == j) {
39: for (k = 0; k < bs; k++) {
40: for (l = 0; l < bs; l++) {
41: Ii = i * bs + k;
42: J = j * bs + l;
43: PetscRandomGetValue(rctx, &v);
44: if (Ii == J) v = PetscRealPart(v + 3 * n * bs);
45: MatSetValue(A, Ii, J, v, INSERT_VALUES);
46: if (test_hermitian) {
47: MatSetValue(A, J, Ii, PetscConj(v), INSERT_VALUES);
48: } else {
49: MatSetValue(A, J, Ii, v, INSERT_VALUES);
50: }
51: }
52: }
53: }
54: }
55: }
56: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
59: /* With complex numbers:
60: - PETSc cholesky does not support hermitian matrices
61: - CHOLMOD only supports hermitian matrices
62: - SUPERLU_DIST seems supporting both
63: */
64: if (test_hermitian) MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE);
66: {
67: Mat M;
68: MatComputeOperator(A, MATAIJ, &M);
69: MatViewFromOptions(M, NULL, "-expl_view");
70: MatDestroy(&M);
71: }
73: A2 = NULL;
74: if (convert) MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &A2);
76: VecCreate(PETSC_COMM_SELF, &u);
77: VecSetSizes(u, PETSC_DECIDE, n * bs);
78: VecSetFromOptions(u);
79: VecDuplicate(u, &b);
80: VecDuplicate(b, &x);
82: VecSet(u, 1.0);
83: MatMult(A, u, b);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Create the linear solver and set various options
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: /*
90: Create linear solver context
91: */
92: KSPCreate(PETSC_COMM_SELF, &ksp);
94: /*
95: Set operators.
96: */
97: KSPSetOperators(ksp, A2 ? A2 : A, A);
99: KSPSetFromOptions(ksp);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Solve the linear system
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: KSPSolve(ksp, b, x);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Check solution and clean up
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: /*
112: Check the error
113: */
114: VecAXPY(x, -1.0, u);
115: VecNorm(x, NORM_2, &norm);
116: KSPGetIterationNumber(ksp, &its);
118: /*
119: Print convergence information. PetscPrintf() produces a single
120: print statement from all processes that share a communicator.
121: An alternative is PetscFPrintf(), which prints to a file.
122: */
123: if (norm > 100 * PETSC_SMALL) PetscPrintf(PETSC_COMM_SELF, "Norm of residual %g iterations %" PetscInt_FMT " bs %" PetscInt_FMT "\n", (double)norm, its, bs);
125: /*
126: Free work space. All PETSc objects should be destroyed when they
127: are no longer needed.
128: */
129: KSPDestroy(&ksp);
130: VecDestroy(&u);
131: VecDestroy(&x);
132: VecDestroy(&b);
133: MatDestroy(&A);
134: MatDestroy(&A2);
135: PetscRandomDestroy(&rctx);
137: /*
138: Always call PetscFinalize() before exiting a program. This routine
139: - finalizes the PETSc libraries as well as MPI
140: - provides summary and diagnostic information if certain runtime
141: options are chosen (e.g., -log_view).
142: */
143: PetscFinalize();
144: return 0;
145: }
147: /*TEST
149: test:
150: args: -mat_type {{aij baij sbaij}} -bs {{1 2 3 4 5 6 7 8 9 10 11 12}} -pc_type cholesky -herm 0 -conv {{0 1}}
152: test:
153: nsize: {{1 4}}
154: suffix: cholmod
155: requires: suitesparse
156: args: -mat_type {{aij sbaij}} -bs 1 -pc_type cholesky -pc_factor_mat_solver_type cholmod -herm -conv {{0 1}}
158: test:
159: nsize: {{1 4}}
160: suffix: superlu_dist
161: requires: superlu_dist
162: output_file: output/ex49_cholmod.out
163: args: -mat_type mpiaij -bs 3 -pc_type cholesky -pc_factor_mat_solver_type superlu_dist -herm {{0 1}} -conv {{0 1}}
165: test:
166: suffix: mkl_pardiso
167: requires: mkl_pardiso
168: output_file: output/ex49_1.out
169: args: -bs {{1 3}} -pc_type cholesky -pc_factor_mat_solver_type mkl_pardiso
171: test:
172: suffix: cg
173: requires: complex
174: output_file: output/ex49_cg.out
175: args: -herm -ksp_cg_type hermitian -mat_type aij -ksp_type cg -pc_type jacobi -ksp_rtol 4e-07
177: test:
178: suffix: pipecg2
179: requires: complex
180: output_file: output/ex49_pipecg2.out
181: args: -herm -mat_type aij -ksp_type pipecg2 -pc_type jacobi -ksp_rtol 4e-07 -ksp_norm_type {{preconditioned unpreconditioned natural}}
183: TEST*/