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
27: PetscInitialize(&argc, &args, (char *)0, help);
28: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
29: 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: MatCreate(PETSC_COMM_WORLD, &A);
35: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
36: MatSetFromOptions(A);
37: MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
38: MatSeqAIJSetPreallocation(A, 5, NULL);
39: MatSetUp(A);
41: MatGetOwnershipRange(A, &Istart, &Iend);
43: PetscLogStageRegister("Assembly", &stage);
44: 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: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
52: }
53: if (i < m - 1) {
54: J = Ii + n;
55: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
56: }
57: if (j > 0) {
58: J = Ii - 1;
59: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
60: }
61: if (j < n - 1) {
62: J = Ii + 1;
63: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
64: }
65: v = 4.0;
66: MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
67: }
68: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
70: PetscLogStagePop();
72: /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
73: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
75: /* Create parallel vectors. */
76: VecCreate(PETSC_COMM_WORLD, &u);
77: VecSetSizes(u, PETSC_DECIDE, m * n);
78: VecSetFromOptions(u);
79: VecDuplicate(u, &b);
80: 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: PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL);
89: if (flg) {
90: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
91: PetscRandomSetFromOptions(rctx);
92: VecSetRandom(u, rctx);
93: PetscRandomDestroy(&rctx);
94: } else {
95: VecSet(u, 1.0);
96: }
97: MatMult(A, u, b);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create the linear solver and set various options
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: /* Create linear solver context */
103: KSPCreate(PETSC_COMM_WORLD, &ksp);
105: /* Set operators. */
106: KSPSetOperators(ksp, A, A);
108: KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
110: PetscOptionsGetBool(NULL, NULL, "-setfromoptions_first", &setfromoptions_first, NULL);
111: if (setfromoptions_first) {
112: /* code path for changing from KSPLSQR to KSPREONLY */
113: KSPSetFromOptions(ksp);
114: }
115: KSPSetType(ksp, KSPPREONLY);
116: KSPGetPC(ksp, &pc);
117: 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: 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: PetscOptionsInsertString(NULL, "-mat_mumps_icntl_13 1");
128: #else
129: MPI_Comm_size(PETSC_COMM_WORLD, &size);
131: #endif
133: if (!setfromoptions_first) {
134: /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
135: KSPSetFromOptions(ksp);
136: }
138: /* get inertia */
139: PetscOptionsGetBool(NULL, NULL, "-solve", &solve, NULL);
140: PetscOptionsGetBool(NULL, NULL, "-sameA", &sameA, NULL);
141: KSPSetUp(ksp);
142: PCFactorGetMatrix(pc, &F);
143: MatGetInertia(F, &in, NULL, NULL);
144: PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in);
145: if (solve) {
146: PetscPrintf(PETSC_COMM_WORLD, "Solving the intermediate KSP\n");
147: KSPSolve(ksp, b, x);
148: } else PetscPrintf(PETSC_COMM_WORLD, "NOT Solving the intermediate KSP\n");
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Solve the linear system
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: MatDuplicate(A, MAT_COPY_VALUES, &B);
154: if (sameA) {
155: PetscPrintf(PETSC_COMM_WORLD, "Setting A\n");
156: MatAXPY(A, 1.1, B, DIFFERENT_NONZERO_PATTERN);
157: KSPSetOperators(ksp, A, A);
158: } else {
159: PetscPrintf(PETSC_COMM_WORLD, "Setting B\n");
160: MatAXPY(B, 1.1, A, DIFFERENT_NONZERO_PATTERN);
161: KSPSetOperators(ksp, B, B);
162: }
163: KSPSetUp(ksp);
164: PCFactorGetMatrix(pc, &F);
165: MatGetInertia(F, &in, NULL, NULL);
166: PetscPrintf(PETSC_COMM_WORLD, "INERTIA=%" PetscInt_FMT "\n", in);
167: KSPSolve(ksp, b, x);
168: MatDestroy(&B);
170: /* Free work space.*/
171: KSPDestroy(&ksp);
172: VecDestroy(&u);
173: VecDestroy(&x);
174: VecDestroy(&b);
175: MatDestroy(&A);
177: 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*/