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