1 static char help[] = "Estimate eigenvalues with KSP.\n\n";
2
3 /*
4 Test example that demonstrates how KSP can estimate eigenvalues.
5
6 Contributed by: Pablo Brubeck <brubeck@protonmail.com>
7 */
8 #include <petscksp.h>
9
main(int argc,char ** args)10 int main(int argc, char **args)
11 {
12 Vec x, b; /* approx solution, RHS */
13 Mat A; /* linear system matrix */
14 KSP ksp; /* linear solver context */
15 PC pc; /* preconditioner context */
16 PetscInt n = 6, i, j, col[2];
17 PetscScalar value[4];
18 PetscMPIInt size;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &args, NULL, help));
22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
24
25 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26 Compute the matrix and right-hand-side vector that define
27 the linear system, Ax = b.
28 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
29
30 /*
31 Create vectors. Note that we form 1 vector from scratch and
32 then duplicate as needed.
33 */
34 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
35 PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
36 PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
37 PetscCall(VecSetFromOptions(x));
38 PetscCall(VecDuplicate(x, &b));
39
40 /*
41 Create matrix. When using MatCreate(), the matrix format can
42 be specified at runtime.
43
44 Performance tuning note: For problems of substantial size,
45 preallocation of matrix memory is crucial for attaining good
46 performance. See the matrix chapter of the users manual for details.
47 */
48 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
49 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
50 PetscCall(MatSetFromOptions(A));
51 PetscCall(MatSetUp(A));
52
53 /*
54 Assemble matrix
55 */
56 value[0] = 2.0;
57 value[1] = -1.0;
58 value[2] = -1.0;
59 value[3] = 2.0;
60 for (i = 0; 2 * i < n; i++) {
61 col[0] = 2 * i;
62 col[1] = 2 * i + 1;
63 PetscCall(MatSetValues(A, 2, col, 2, col, value, INSERT_VALUES));
64 for (j = 0; j < 4; j++) value[j] *= 3.0;
65 }
66 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
67 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
68
69 /*
70 Set random right-hand-side vector.
71 */
72 PetscCall(VecSetRandom(b, NULL));
73
74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75 Create the linear solver and set various options
76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77 /*
78 Create linear solver context
79 */
80 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
81
82 /*
83 Set operators. Here the matrix that defines the linear system
84 also serves as the matrix from which the preconditioner is constructed.
85 */
86 PetscCall(KSPSetOperators(ksp, A, A));
87
88 /*
89 Set linear solver defaults for this problem (optional).
90 - By extracting the KSP and PC contexts from the KSP context,
91 we can then directly call any KSP and PC routines to set
92 various options.
93 - The following four statements are optional; all of these
94 parameters could alternatively be specified at runtime via
95 KSPSetFromOptions();
96 */
97 PetscCall(KSPGetPC(ksp, &pc));
98 PetscCall(PCSetType(pc, PCJACOBI));
99 PetscCall(KSPSetTolerances(ksp, 1.e-5, 1.e-5, PETSC_CURRENT, PETSC_CURRENT));
100
101 /*
102 Set runtime options, e.g.,
103 -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
104 These options will override those specified above as long as
105 KSPSetFromOptions() is called _after_ any other customization
106 routines.
107 */
108 PetscCall(KSPSetFromOptions(ksp));
109
110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111 Solve the linear system
112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113 /*
114 Solve linear system
115 */
116 PetscCall(KSPSolve(ksp, b, x));
117
118 /*
119 View solver info; we could instead use the option -ksp_view to
120 print this info to the screen at the conclusion of KSPSolve().
121 */
122 PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
123
124 /*
125 Free work space. All PETSc objects should be destroyed when they
126 are no longer needed.
127 */
128 PetscCall(VecDestroy(&x));
129 PetscCall(VecDestroy(&b));
130 PetscCall(MatDestroy(&A));
131 PetscCall(KSPDestroy(&ksp));
132
133 /*
134 Always call PetscFinalize() before exiting a program. This routine
135 - finalizes the PETSc libraries as well as MPI
136 - provides summary and diagnostic information if certain runtime
137 options are chosen (e.g., -log_view).
138 */
139 PetscCall(PetscFinalize());
140 return 0;
141 }
142
143 /*TEST
144
145 test:
146 suffix: cg_none
147 filter: grep -v "variant HERMITIAN"
148 args: -ksp_type cg -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
149
150 test:
151 suffix: cg_jacobi
152 filter: grep -v "variant HERMITIAN"
153 args: -ksp_type cg -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
154
155 test:
156 suffix: fcg_none
157 args: -ksp_type fcg -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
158
159 test:
160 suffix: fcg_jacobi
161 args: -ksp_type fcg -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
162
163 test:
164 suffix: minres_none
165 args: -ksp_type minres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
166
167 test:
168 suffix: minres_jacobi
169 args: -ksp_type minres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
170
171 test:
172 suffix: gmres_none
173 args: -ksp_type gmres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
174
175 test:
176 suffix: gmres_jacobi_left
177 args: -ksp_type gmres -ksp_pc_side left -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
178
179 test:
180 suffix: gmres_jacobi_right
181 args: -ksp_type gmres -ksp_pc_side right -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
182
183 test:
184 suffix: fgmres_none
185 args: -ksp_type fgmres -pc_type none -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
186
187 test:
188 suffix: fgmres_jacobi
189 args: -ksp_type fgmres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -ksp_converged_reason
190
191 TEST*/
192