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 !PetscDefined(HAVE_MUMPS)
20: PetscMPIInt size;
21: #endif
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));
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_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
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 PetscDefined(HAVE_MUMPS)
116: PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_WORLD, PETSC_ERR_SUP, "Spectrum slicing with MUMPS is not available for complex scalars");
117: PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
118: /*
119: must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for
120: matrix inertia), currently there is no better way of setting this in program
121: */
122: PetscCall(PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1"));
123: #else
124: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
125: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Configure with MUMPS if you want to run this example in parallel");
126: #endif
128: if (!setfromoptions_first) {
129: /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
130: PetscCall(KSPSetFromOptions(ksp));
131: }
133: /* get inertia */
134: PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve", &solve, NULL));
135: PetscCall(PetscOptionsGetBool(NULL, NULL, "-sameA", &sameA, NULL));
136: PetscCall(KSPSetUp(ksp));
137: PetscCall(PCFactorGetMatrix(pc, &F));
138: PetscCall(MatGetInertia(F, &in, NULL, NULL));
139: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
140: if (solve) {
141: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving the intermediate KSP\n"));
142: PetscCall(KSPSolve(ksp, b, x));
143: } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "NOT Solving the intermediate KSP\n"));
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Solve the linear system
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
149: if (sameA) {
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting A\n"));
151: PetscCall(MatAXPY(A, 1.1, B, DIFFERENT_NONZERO_PATTERN));
152: PetscCall(KSPSetOperators(ksp, A, A));
153: } else {
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Setting B\n"));
155: PetscCall(MatAXPY(B, 1.1, A, DIFFERENT_NONZERO_PATTERN));
156: PetscCall(KSPSetOperators(ksp, B, B));
157: }
158: PetscCall(KSPSetUp(ksp));
159: PetscCall(PCFactorGetMatrix(pc, &F));
160: PetscCall(MatGetInertia(F, &in, NULL, NULL));
161: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in));
162: PetscCall(KSPSolve(ksp, b, x));
163: PetscCall(MatDestroy(&B));
165: /* Free work space.*/
166: PetscCall(KSPDestroy(&ksp));
167: PetscCall(VecDestroy(&u));
168: PetscCall(VecDestroy(&x));
169: PetscCall(VecDestroy(&b));
170: PetscCall(MatDestroy(&A));
172: PetscCall(PetscFinalize());
173: return 0;
174: }
176: /*TEST
178: build:
179: requires: !complex
181: test:
182: args:
184: test:
185: suffix: 2
186: args: -sameA
188: test:
189: suffix: 3
190: args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}
192: TEST*/