MFEM  v4.5.2
Finite element discretization library
slepc.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 #ifdef MFEM_USE_PETSC
16 #ifdef MFEM_USE_SLEPC
17 
18 #include "linalg.hpp"
19 #include "petscinternals.hpp"
20 #include "slepc.h"
21 
22 static PetscErrorCode ierr;
23 
24 using namespace std;
25 
26 namespace mfem
27 {
28 
30 {
31  MFEMInitializeSlepc(NULL,NULL,NULL,NULL);
32 }
33 
34 void MFEMInitializeSlepc(int *argc,char*** argv)
35 {
36  MFEMInitializeSlepc(argc,argv,NULL,NULL);
37 }
38 
39 void MFEMInitializeSlepc(int *argc,char ***argv,const char rc_file[],
40  const char help[])
41 {
43  {
44  // Tell PETSc to use the same CUDA device as MFEM:
45  ierr = PetscOptionsSetValue(NULL,"-cuda_device",
46  to_string(mfem::Device::GetId()).c_str());
47  MFEM_VERIFY(!ierr,"Unable to set initial option value to PETSc");
48  }
49  ierr = SlepcInitialize(argc,argv,rc_file,help);
50  MFEM_VERIFY(!ierr,"Unable to initialize SLEPc");
51 }
52 
54 {
55  ierr = SlepcFinalize();
56  MFEM_VERIFY(!ierr,"Unable to finalize SLEPc");
57 }
58 
59 
60 SlepcEigenSolver::SlepcEigenSolver(MPI_Comm comm, const std::string &prefix)
61 {
62  clcustom = false;
63  VR = NULL;
64  VC = NULL;
65 
66  ierr = EPSCreate(comm,&eps); CCHKERRQ(comm,ierr);
67  ierr = EPSSetOptionsPrefix(eps, prefix.c_str()); PCHKERRQ(eps, ierr);
68 }
69 
70 SlepcEigenSolver::~SlepcEigenSolver()
71 {
72  delete VR;
73  delete VC;
74 
75  MPI_Comm comm;
76  ierr = PetscObjectGetComm((PetscObject)eps,&comm); PCHKERRQ(eps,ierr);
77  ierr = EPSDestroy(&eps); CCHKERRQ(comm,ierr);
78 }
79 
80 
81 void SlepcEigenSolver::SetOperator(const PetscParMatrix &op)
82 {
83  delete VR;
84  delete VC;
85  VR = VC = NULL;
86 
87  ierr = EPSSetOperators(eps,op,NULL); PCHKERRQ(eps, ierr);
88 
89  VR = new PetscParVector(op, true, false);
90  VC = new PetscParVector(op, true, false);
91 }
92 
93 void SlepcEigenSolver::SetOperators(const PetscParMatrix &op,
94  const PetscParMatrix&opB)
95 {
96  delete VR;
97  delete VC;
98  VR = VC = NULL;
99 
100  ierr = EPSSetOperators(eps,op,opB); PCHKERRQ(eps,ierr);
101 
102  VR = new PetscParVector(op, true, false);
103  VC = new PetscParVector(op, true, false);
104 }
105 
106 void SlepcEigenSolver::SetTol(double tol)
107 {
108  PetscInt max_its;
109 
110  ierr = EPSGetTolerances(eps,NULL,&max_its); PCHKERRQ(eps,ierr);
111  // Work around uninitialized maximum iterations
112  if (max_its==0) { max_its = PETSC_DECIDE; }
113  ierr = EPSSetTolerances(eps,tol,max_its); PCHKERRQ(eps,ierr);
114 }
115 
116 void SlepcEigenSolver::SetMaxIter(int max_its)
117 {
118  double tol;
119 
120  ierr = EPSGetTolerances(eps,&tol,NULL); PCHKERRQ(eps,ierr);
121  ierr = EPSSetTolerances(eps,tol,max_its); PCHKERRQ(eps,ierr);
122 }
123 
124 void SlepcEigenSolver::SetNumModes(int num_eigs)
125 {
126  ierr = EPSSetDimensions(eps,num_eigs,PETSC_DECIDE,PETSC_DECIDE);
127  PCHKERRQ(eps,ierr);
128 }
129 
130 void SlepcEigenSolver::Solve()
131 {
132  Customize();
133 
134  ierr = EPSSolve(eps); PCHKERRQ(eps,ierr);
135 }
136 
137 void SlepcEigenSolver::Customize(bool customize) const
138 {
139  if (!customize) {clcustom = true; }
140  if (!clcustom)
141  {
142  ierr = EPSSetFromOptions(eps); PCHKERRQ(eps,ierr);
143  }
144  clcustom = true;
145 }
146 
147 void SlepcEigenSolver::GetEigenvalue(unsigned int i, double & lr) const
148 {
149  ierr = EPSGetEigenvalue(eps,i,&lr,NULL); PCHKERRQ(eps,ierr);
150 }
151 
152 void SlepcEigenSolver::GetEigenvalue(unsigned int i, double & lr,
153  double & lc) const
154 {
155  ierr = EPSGetEigenvalue(eps,i,&lr,&lc); PCHKERRQ(eps,ierr);
156 }
157 
158 void SlepcEigenSolver::GetEigenvector(unsigned int i, Vector & vr) const
159 {
160  MFEM_VERIFY(VR,"Missing real vector");
161 
162  MFEM_ASSERT(vr.Size() == VR->Size(), "invalid vr.Size() = " << vr.Size()
163  << ", expected size = " << VR->Size());
164 
165  VR->PlaceMemory(vr.GetMemory());
166  ierr = EPSGetEigenvector(eps,i,*VR,NULL); PCHKERRQ(eps,ierr);
167  VR->ResetMemory();
168 
169 }
170 
171 void SlepcEigenSolver::GetEigenvector(unsigned int i, Vector & vr,
172  Vector & vc) const
173 {
174  MFEM_VERIFY(VR,"Missing real vector");
175  MFEM_VERIFY(VC,"Missing imaginary vector");
176  MFEM_ASSERT(vr.Size() == VR->Size(), "invalid vr.Size() = " << vr.Size()
177  << ", expected size = " << VR->Size());
178  MFEM_ASSERT(vc.Size() == VC->Size(), "invalid vc.Size() = " << vc.Size()
179  << ", expected size = " << VC->Size());
180 
181  VR->PlaceArray(vr.GetMemory());
182  VC->PlaceArray(vc.GetMemory());
183  ierr = EPSGetEigenvector(eps,i,*VR,*VC); PCHKERRQ(eps,ierr);
184  VR->ResetMemory();
185  VC->ResetMemory();
186 }
187 
188 int SlepcEigenSolver::GetNumConverged()
189 {
190  PetscInt num_conv;
191  ierr = EPSGetConverged(eps,&num_conv); PCHKERRQ(eps,ierr);
192  return static_cast<int>(num_conv);
193 }
194 
195 void SlepcEigenSolver::SetWhichEigenpairs(SlepcEigenSolver::Which which)
196 {
197  switch (which)
198  {
199  case SlepcEigenSolver::LARGEST_MAGNITUDE:
200  ierr = EPSSetWhichEigenpairs(eps,EPS_LARGEST_MAGNITUDE); PCHKERRQ(eps,ierr);
201  break;
202  case SlepcEigenSolver::SMALLEST_MAGNITUDE:
203  ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE); PCHKERRQ(eps,ierr);
204  break;
205  case SlepcEigenSolver::LARGEST_REAL:
206  ierr = EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL); PCHKERRQ(eps,ierr);
207  break;
208  case SlepcEigenSolver::SMALLEST_REAL:
209  ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL); PCHKERRQ(eps,ierr);
210  break;
211  case SlepcEigenSolver::LARGEST_IMAGINARY:
212  ierr = EPSSetWhichEigenpairs(eps,EPS_LARGEST_IMAGINARY); PCHKERRQ(eps,ierr);
213  break;
214  case SlepcEigenSolver::SMALLEST_IMAGINARY:
215  ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_IMAGINARY); PCHKERRQ(eps,ierr);
216  break;
217  case SlepcEigenSolver::TARGET_MAGNITUDE:
218  ierr = EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE); PCHKERRQ(eps,ierr);
219  break;
220  case SlepcEigenSolver::TARGET_REAL:
221  ierr = EPSSetWhichEigenpairs(eps,EPS_TARGET_REAL); PCHKERRQ(eps,ierr);
222  break;
223  default:
224  MFEM_ABORT("Which eigenpair not implemented!");
225  break;
226  }
227 }
228 
229 void SlepcEigenSolver::SetTarget(double target)
230 {
231  ierr = EPSSetTarget(eps,target); PCHKERRQ(eps,ierr);
232 }
233 
234 void SlepcEigenSolver::SetSpectralTransformation(
236 {
237  ST st;
238  ierr = EPSGetST(eps,&st); PCHKERRQ(eps,ierr);
239  switch (transformation)
240  {
241  case SlepcEigenSolver::SHIFT:
242  ierr = STSetType(st,STSHIFT); PCHKERRQ(eps,ierr);
243  break;
244  case SlepcEigenSolver::SHIFT_INVERT:
245  ierr = STSetType(st,STSINVERT); PCHKERRQ(eps,ierr);
246  break;
247  default:
248  MFEM_ABORT("Spectral transformation not implemented!");
249  break;
250  }
251 }
252 
253 }
254 
255 #endif // MFEM_USE_SLEPC
256 #endif // MFEM_USE_PETSC
257 #endif // MFEM_USE_MPI
static int GetId()
Get the device id of the configured device.
Definition: device.hpp:252
HYPRE_Int PetscInt
Definition: petsc.hpp:47
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:315
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
STL namespace.
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:230
void MFEMInitializeSlepc(int *argc, char ***argv, const char rc_file[], const char help[])
Definition: slepc.cpp:39
void transformation(const Vector &p, Vector &v)
void MFEMFinalizeSlepc()
Definition: slepc.cpp:53
Biwise-OR of all CUDA backends.
Definition: device.hpp:88
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Definition: device.hpp:258
struct _p_PetscObject * PetscObject
Definition: petsc.hpp:51
Vector data type.
Definition: vector.hpp:60