Actual source code: ex8.c
1: static char help[] = "Solves a linear system in parallel with KSP. \n\
2: Contributed by Jose E. Roman, SLEPc developer, for testing repeated call of KSPSetOperators(), 2014 \n\n";
4: #include <petscksp.h>
5: int main(int argc, char **args)
6: {
7: Vec x, b, u; /* approx solution, RHS, exact solution */
8: Mat A; /* linear system matrix */
9: KSP ksp; /* linear solver context */
10: PetscRandom rctx; /* random number generator context */
11: PetscInt i, j, Ii, J, Istart, Iend, m = 8, n = 7;
12: PetscBool flg = PETSC_FALSE;
13: PetscScalar v;
14: PC pc;
15: PetscInt in;
16: Mat F, B;
17: PetscBool solve = PETSC_FALSE, sameA = PETSC_FALSE, setfromoptions_first = PETSC_FALSE;
18: PetscLogStage stage;
19: #if !defined(PETSC_HAVE_MUMPS)
20: PetscMPIInt size;
21: #endif
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
26: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
27: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28: Compute the matrix and right-hand-side vector that define
29: the linear system, Ax = b.
30: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
33: PetscCall(MatSetFromOptions(A));
34: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
35: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
36: PetscCall(MatSetUp(A));
38: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
40: PetscCall(PetscLogStageRegister("Assembly", &stage));
41: PetscCall(PetscLogStagePush(stage));
42: for (Ii = Istart; Ii < Iend; Ii++) {
43: v = -1.0;
44: i = Ii / n;
45: j = Ii - i * n;
46: if (i > 0) {
47: J = Ii - n;
48: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
49: }
50: if (i < m - 1) {
51: J = Ii + n;
52: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
53: }
54: if (j > 0) {
55: J = Ii - 1;
56: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
57: }
58: if (j < n - 1) {
59: J = Ii + 1;
60: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
61: }
62: v = 4.0;
63: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
64: }
65: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
67: PetscCall(PetscLogStagePop());
69: /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
70: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
72: /* Create parallel vectors. */
73: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
74: PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
75: PetscCall(VecSetFromOptions(u));
76: PetscCall(VecDuplicate(u, &b));
77: PetscCall(VecDuplicate(b, &x));
79: /*
80: Set exact solution; then compute right-hand-side vector.
81: By default we use an exact solution of a vector with all
82: elements of 1.0; Alternatively, using the runtime option
83: -random_sol forms a solution vector with random components.
84: */
85: PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL));
86: if (flg) {
87: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
88: PetscCall(PetscRandomSetFromOptions(rctx));
89: PetscCall(VecSetRandom(u, rctx));
90: PetscCall(PetscRandomDestroy(&rctx));
91: } else {
92: PetscCall(VecSet(u, 1.0));
93: }
94: PetscCall(MatMult(A, u, b));
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create the linear solver and set various options
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: /* Create linear solver context */
100: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
102: /* Set operators. */
103: PetscCall(KSPSetOperators(ksp, A, A));
105: PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
107: PetscCall(PetscOptionsGetBool(NULL, NULL, "-setfromoptions_first", &setfromoptions_first, NULL));
108: if (setfromoptions_first) {
109: /* code path for changing from KSPLSQR to KSPREONLY */
110: PetscCall(KSPSetFromOptions(ksp));
111: }
112: PetscCall(KSPSetType(ksp, KSPPREONLY));
113: PetscCall(KSPGetPC(ksp, &pc));
114: PetscCall(PCSetType(pc, PCCHOLESKY));
115: #if defined(PETSC_HAVE_MUMPS)
116: #if defined(PETSC_USE_COMPLEX)
117: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Spectrum slicing with MUMPS is not available for complex scalars");
118: #endif
119: PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
120: /*
121: must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for
122: matrix inertia), currently there is no better way of setting this in program
123: */
124: PetscCall(PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1"));
125: #else
126: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
127: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Configure with MUMPS if you want to run this example in parallel");
128: #endif
130: if (!setfromoptions_first) {
131: /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
132: PetscCall(KSPSetFromOptions(ksp));
133: }
135: /* get inertia */
136: PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve", &solve, NULL));
137: PetscCall(PetscOptionsGetBool(NULL, NULL, "-sameA", &sameA, NULL));
138: PetscCall(KSPSetUp(ksp));
139: PetscCall(PCFactorGetMatrix(pc, &F));
140: PetscCall(MatGetInertia(F, &in, NULL, NULL));
141: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
142: if (solve) {
143: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving the intermediate KSP\n"));
144: PetscCall(KSPSolve(ksp, b, x));
145: } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NOT Solving the intermediate KSP\n"));
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Solve the linear system
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
151: if (sameA) {
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting A\n"));
153: PetscCall(MatAXPY(A, 1.1, B, DIFFERENT_NONZERO_PATTERN));
154: PetscCall(KSPSetOperators(ksp, A, A));
155: } else {
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting B\n"));
157: PetscCall(MatAXPY(B, 1.1, A, DIFFERENT_NONZERO_PATTERN));
158: PetscCall(KSPSetOperators(ksp, B, B));
159: }
160: PetscCall(KSPSetUp(ksp));
161: PetscCall(PCFactorGetMatrix(pc, &F));
162: PetscCall(MatGetInertia(F, &in, NULL, NULL));
163: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
164: PetscCall(KSPSolve(ksp, b, x));
165: PetscCall(MatDestroy(&B));
167: /* Free work space.*/
168: PetscCall(KSPDestroy(&ksp));
169: PetscCall(VecDestroy(&u));
170: PetscCall(VecDestroy(&x));
171: PetscCall(VecDestroy(&b));
172: PetscCall(MatDestroy(&A));
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: /*TEST
180: build:
181: requires: !complex
183: test:
184: args:
186: test:
187: suffix: 2
188: args: -sameA
190: test:
191: suffix: 3
192: args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}
194: TEST*/