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