1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 3 4 // This library is free software; you can redistribute it and/or 5 // modify it under the terms of the GNU Lesser General Public 6 // License as published by the Free Software Foundation; either 7 // version 2.1 of the License, or (at your option) any later version. 8 9 // This library is distributed in the hope that it will be useful, 10 // but WITHOUT ANY WARRANTY; without even the implied warranty of 11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12 // Lesser General Public License for more details. 13 14 // You should have received a copy of the GNU Lesser General Public 15 // License along with this library; if not, write to the Free Software 16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 18 19 20 #ifndef LIBMESH_SLEPC_EIGEN_SOLVER_H 21 #define LIBMESH_SLEPC_EIGEN_SOLVER_H 22 23 #include "libmesh/libmesh_config.h" 24 25 #ifdef LIBMESH_HAVE_SLEPC 26 27 // Local includes 28 #include "libmesh/eigen_solver.h" 29 #include "libmesh/slepc_macro.h" 30 31 // SLEPc include files. 32 EXTERN_C_FOR_SLEPC_BEGIN 33 # include "libmesh/ignore_warnings.h" 34 # include <slepceps.h> 35 # include "libmesh/restore_warnings.h" 36 EXTERN_C_FOR_SLEPC_END 37 38 namespace libMesh 39 { 40 template <typename T> class PetscVector; 41 42 /** 43 * This class provides an interface to the SLEPc 44 * eigenvalue solver library from http://slepc.upv.es/. 45 * 46 * \author Steffen Peterson 47 * \date 2005 48 * \brief EigenSolver implementation based on SLEPc. 49 */ 50 template <typename T> 51 class SlepcEigenSolver : public EigenSolver<T> 52 { 53 54 public: 55 56 /** 57 * Constructor. Initializes Petsc data structures 58 */ 59 SlepcEigenSolver(const Parallel::Communicator & comm_in); 60 61 62 /** 63 * Destructor. 64 */ 65 ~SlepcEigenSolver(); 66 67 68 /** 69 * Release all memory and clear data structures. 70 */ 71 virtual void clear() override; 72 73 74 /** 75 * Initialize data structures if not done so already. 76 */ 77 virtual void init() override; 78 79 80 /** 81 * This function calls the SLEPc solver to compute 82 * the eigenpairs of the SparseMatrix matrix_A. \p nev is 83 * the number of eigenpairs to be computed and 84 * \p ncv is the number of basis vectors to be 85 * used in the solution procedure. Return values 86 * are the number of converged eigen values and the 87 * number of the iterations carried out by the eigen 88 * solver. 89 */ 90 virtual std::pair<unsigned int, unsigned int> 91 solve_standard (SparseMatrix<T> & matrix_A, 92 int nev, 93 int ncv, 94 const double tol, 95 const unsigned int m_its) override; 96 97 /** 98 * Same as above except that matrix_A is a ShellMatrix 99 * in this case. 100 */ 101 virtual std::pair<unsigned int, unsigned int> 102 solve_standard (ShellMatrix<T> & shell_matrix, 103 int nev, 104 int ncv, 105 const double tol, 106 const unsigned int m_its) override; 107 108 /** 109 * Same as above except that matrix_A is a ShellMatrix 110 * in this case. 111 */ 112 virtual std::pair<unsigned int, unsigned int> 113 solve_standard (ShellMatrix<T> & shell_matrix, 114 SparseMatrix<T> & precond, 115 int nev, 116 int ncv, 117 const double tol, 118 const unsigned int m_its) override; 119 120 /** 121 * Same as above except that precond is a ShellMatrix 122 * in this case. 123 */ 124 virtual std::pair<unsigned int, unsigned int> 125 solve_standard (ShellMatrix<T> & shell_matrix, 126 ShellMatrix<T> & precond, 127 int nev, 128 int ncv, 129 const double tol, 130 const unsigned int m_its) override; 131 132 133 /** 134 * This function calls the SLEPc solver to compute 135 * the eigenpairs for the generalized eigenproblem 136 * defined by the matrix_A and matrix_B, 137 * which are of type SparseMatrix. The argument 138 * \p nev is the number of eigenpairs to be computed 139 * and \p ncv is the number of basis vectors to be 140 * used in the solution procedure. Return values 141 * are the number of converged eigen values and the 142 * number of the iterations carried out by the eigen 143 * solver. 144 */ 145 virtual std::pair<unsigned int, unsigned int> 146 solve_generalized(SparseMatrix<T> & matrix_A, 147 SparseMatrix<T> & matrix_B, 148 int nev, 149 int ncv, 150 const double tol, 151 const unsigned int m_its) override; 152 153 /** 154 * Solve generalized eigenproblem when matrix_A is of 155 * type ShellMatrix, matrix_B is of type SparseMatrix. 156 */ 157 virtual std::pair<unsigned int, unsigned int> 158 solve_generalized(ShellMatrix<T> & matrix_A, 159 SparseMatrix<T> & matrix_B, 160 int nev, 161 int ncv, 162 const double tol, 163 const unsigned int m_its) override; 164 165 /** 166 * Solve generalized eigenproblem when matrix_A is of 167 * type SparseMatrix, matrix_B is of type ShellMatrix. 168 * When using this function, one should use the 169 * command line options: 170 * -st_ksp_type gmres -st_pc_type none 171 * or 172 * -st_ksp_type gmres -st_pc_type jacobi 173 * or similar. 174 */ 175 virtual std::pair<unsigned int, unsigned int> 176 solve_generalized(SparseMatrix<T> & matrix_A, 177 ShellMatrix<T> & matrix_B, 178 int nev, 179 int ncv, 180 const double tol, 181 const unsigned int m_its) override; 182 183 /** 184 * Solve generalized eigenproblem when both matrix_A and 185 * matrix_B are of type ShellMatrix. 186 * When using this function, one should use the 187 * command line options: 188 * -st_ksp_type gmres -st_pc_type none 189 * or 190 * -st_ksp_type gmres -st_pc_type jacobi 191 * or similar. 192 */ 193 virtual std::pair<unsigned int, unsigned int> 194 solve_generalized(ShellMatrix<T> & matrix_A, 195 ShellMatrix<T> & matrix_B, 196 int nev, 197 int ncv, 198 const double tol, 199 const unsigned int m_its) override; 200 201 202 virtual std::pair<unsigned int, unsigned int> 203 solve_generalized(ShellMatrix<T> & matrix_A, 204 ShellMatrix<T> & matrix_B, 205 SparseMatrix<T> & precond, 206 int nev, 207 int ncv, 208 const double tol, 209 const unsigned int m_its) override; 210 211 virtual std::pair<unsigned int, unsigned int> 212 solve_generalized(ShellMatrix<T> & matrix_A, 213 ShellMatrix<T> & matrix_B, 214 ShellMatrix<T> & precond, 215 int nev, 216 int ncv, 217 const double tol, 218 const unsigned int m_its) override; 219 220 /** 221 * \returns The real and imaginary part of the ith eigenvalue and 222 * copies the respective eigenvector to the solution vector. 223 * 224 * \note The eigenpair may be complex even for real-valued matrices. 225 */ 226 virtual std::pair<Real, Real> 227 get_eigenpair (dof_id_type i, 228 NumericVector<T> & solution_in) override; 229 230 /** 231 * Same as above, but does not copy the eigenvector. 232 */ 233 virtual std::pair<Real, Real> 234 get_eigenvalue (dof_id_type i) override; 235 236 /** 237 * \returns The relative error \f$ ||A x - \lambda x|| / |\lambda x| \f$ 238 * of the ith eigenpair (or the equivalent for a general eigenvalue problem). 239 */ 240 Real get_relative_error (unsigned int i); 241 242 /** 243 * Attach a deflation space defined by a single vector. 244 */ 245 virtual void attach_deflation_space(NumericVector<T> & deflation_vector) override; 246 247 /** 248 * Use \p initial_space_in as the initial guess. 249 */ 250 virtual void 251 set_initial_space(NumericVector<T> & initial_space_in) override; 252 253 /** 254 * \returns The raw SLEPc \p EPS pointer. 255 */ eps()256 EPS eps() { this->init(); return _eps; } 257 258 private: 259 260 /** 261 * Helper function that actually performs the standard eigensolve. 262 */ 263 std::pair<unsigned int, unsigned int> _solve_standard_helper (Mat mat, 264 Mat precond, 265 int nev, 266 int ncv, 267 const double tol, 268 const unsigned int m_its); 269 270 /** 271 * Helper function that actually performs the generalized eigensolve. 272 */ 273 std::pair<unsigned int, unsigned int> _solve_generalized_helper (Mat mat_A, 274 Mat mat_B, 275 Mat precond, 276 int nev, 277 int ncv, 278 const double tol, 279 const unsigned int m_its); 280 281 /** 282 * Tells Slepc to use the user-specified solver stored in 283 * \p _eigen_solver_type 284 */ 285 void set_slepc_solver_type (); 286 287 /** 288 * Tells Slepc to deal with the type of problem stored in 289 * \p _eigen_problem_type 290 */ 291 void set_slepc_problem_type (); 292 293 /** 294 * Tells Slepc to compute the spectrum at the position 295 * stored in \p _position_of_spectrum 296 */ 297 void set_slepc_position_of_spectrum(); 298 299 /** 300 * Internal function if shell matrix mode is used, this just 301 * calls the shell matrix's matrix multiplication function. 302 * See PetscLinearSolver for a similar implementation. 303 */ 304 static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest); 305 306 /** 307 * Internal function if shell matrix mode is used, this just 308 * calls the shell matrix's get_diagonal function. 309 * Required in order to use Jacobi preconditioning. 310 */ 311 static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest); 312 313 /** 314 * Eigenproblem solver context 315 */ 316 EPS _eps; 317 318 /** 319 * A vector used for initial space. The vector will be used as the basis for EPS. 320 */ 321 PetscVector<T>* _initial_space; 322 }; 323 324 } // namespace libMesh 325 326 327 #endif // #ifdef LIBMESH_HAVE_SLEPC 328 #endif // LIBMESH_SLEPC_EIGEN_SOLVER_H 329