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