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";
3
4 #include <petscksp.h>
main(int argc,char ** args)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
22
23 PetscFunctionBeginUser;
24 PetscCall(PetscInitialize(&argc, &args, NULL, 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));
37
38 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
39
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());
68
69 /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
70 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
71
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));
78
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));
95
96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 Create the linear solver and set various options
98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 /* Create linear solver context */
100 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
101
102 /* Set operators. */
103 PetscCall(KSPSetOperators(ksp, A, A));
104
105 PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
106
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_WRONG_MPI_SIZE, "Configure with MUMPS if you want to run this example in parallel");
128 #endif
129
130 if (!setfromoptions_first) {
131 /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
132 PetscCall(KSPSetFromOptions(ksp));
133 }
134
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"));
146
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));
166
167 /* Free work space.*/
168 PetscCall(KSPDestroy(&ksp));
169 PetscCall(VecDestroy(&u));
170 PetscCall(VecDestroy(&x));
171 PetscCall(VecDestroy(&b));
172 PetscCall(MatDestroy(&A));
173
174 PetscCall(PetscFinalize());
175 return 0;
176 }
177
178 /*TEST
179
180 build:
181 requires: !complex
182
183 test:
184 args:
185
186 test:
187 suffix: 2
188 args: -sameA
189
190 test:
191 suffix: 3
192 args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}
193
194 TEST*/
195