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 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