TBTK
ArnoldiIterator.h
1 /* Copyright 2016 Kristofer Bj√∂rnson
2  *
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
23 #ifndef COM_DAFER45_SOLVER_ARNOLDI_ITERATOR
24 #define COM_DAFER45_SOLVER_ARNOLDI_ITERATOR
25 
26 #include "TBTK/Model.h"
27 #include "TBTK/Solver/LUSolver.h"
28 #include "TBTK/Solver/Solver.h"
29 
30 #include <complex>
31 
32 namespace TBTK{
33 namespace Solver{
34 
40 class ArnoldiIterator : public Solver, public Communicator{
41 public:
44 
46  virtual ~ArnoldiIterator();
47 
57  enum class Mode {Normal, ShiftAndInvert};
58 
62  void setMode(Mode mode);
63 
67  Mode getMode() const;
68 
72  void setNumEigenValues(int numEigenValues);
73 
77  int getNumEigenValues() const;
78 
83  void setCalculateEigenVectors(bool calculateEigenVectors);
84 
88  bool getCalculateEigenVectors() const;
89 
94  void setNumLanczosVectors(int numLanczosVectors);
95 
100  int getNumLanczosVectors() const;
101 
107  void setTolerance(double tolerance);
108 
114  void setMaxIterations(int maxIterations);
115 
129  void setCentralValue(double centralValue);
130 
132  void run();
133 
135  const std::complex<double>* getEigenValues() const;
136 
138  const double getEigenValue(int state) const;
139 
144  const std::complex<double> getAmplitude(int state, const Index &index);
145 private:
147  Mode mode;
148 
150  int numEigenValues;
151 
154  bool calculateEigenVectors;
155 
158  int numLanczosVectors;
159 
162  double shift;
163 
166  double tolerance;
167 
170  int maxIterations;
171 
173  std::complex<double> *residuals;
174 
176  std::complex<double> *eigenValues;
177 
179  std::complex<double> *eigenVectors;
180 
182  LUSolver luSolver;
183 
186  void initNormal();
187 
190  void initShiftAndInvert();
191 
193  void arnoldiLoop();
194 
196  void checkZnaupdInfo(int info) const;
197 
199  bool executeReverseCommunicationMessage(
200  int ido,
201  int basisSize,
202  double *workd,
203  int *ipntr,
204  Matrix<double> &b
205  );
206 
208  bool executeReverseCommunicationMessage(
209  int ido,
210  int basisSize,
211  std::complex<double> *workd,
212  int *ipntr,
213  Matrix<std::complex<double>> &b
214  );
215 
216  void checkZneupdIerr(int ierr) const;
217 
220  void sort();
221 
223  void mergeSortSplit(
224  std::complex<double> *dataIn,
225  std::complex<double> *dataOut,
226  int *orderIn,
227  int *orderOut,
228  int first,
229  int end
230  );
231 
233  void mergeSortMerge(
234  std::complex<double> *dataIn,
235  std::complex<double> *dataOut,
236  int *orderIn,
237  int *orderOut,
238  int first,
239  int middle,
240  int end
241  );
242 };
243 
244 inline void ArnoldiIterator::setMode(Mode mode){
245  this->mode = mode;
246 }
247 
249  return mode;
250 }
251 
252 inline void ArnoldiIterator::setNumEigenValues(int numEigenValues){
253  this->numEigenValues = numEigenValues;
254 }
255 
257  return numEigenValues;
258 }
259 
260 inline void ArnoldiIterator::setCalculateEigenVectors(bool calculateEigenVectors){
261  this->calculateEigenVectors = calculateEigenVectors;
262 }
263 
265  return calculateEigenVectors;
266 }
267 
268 inline void ArnoldiIterator::setNumLanczosVectors(int numLanczosVectors){
269  this->numLanczosVectors = numLanczosVectors;
270 }
271 
273  return numLanczosVectors;
274 }
275 
276 inline void ArnoldiIterator::setTolerance(double tolerance){
277  this->tolerance = tolerance;
278 }
279 
280 inline void ArnoldiIterator::setMaxIterations(int maxIterations){
281  this->maxIterations = maxIterations;
282 }
283 
284 inline void ArnoldiIterator::setCentralValue(double centralValue){
285  shift = centralValue;
286 }
287 
288 inline const std::complex<double>* ArnoldiIterator::getEigenValues() const{
289  return eigenValues;
290 }
291 
292 inline const double ArnoldiIterator::getEigenValue(int state) const{
293  return real(eigenValues[state]);
294 }
295 
296 inline const std::complex<double> ArnoldiIterator::getAmplitude(
297  int state,
298  const Index &index
299 ){
300  const Model &model = getModel();
301  return eigenVectors[model.getBasisSize()*state + model.getBasisIndex(index)];
302 }
303 
304 }; //End of namespace Solver
305 }; //End of namesapce TBTK
306 
307 #endif
Base class for Solvers.
Definition: Solver.h:32
Model & getModel()
Definition: Solver.h:63
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: ArnoldiIterator.h:296
bool getCalculateEigenVectors() const
Definition: ArnoldiIterator.h:264
void setMode(Mode mode)
Definition: ArnoldiIterator.h:244
Container of Model related information.
int getBasisIndex(const Index &index) const
Definition: Model.h:296
Solves a Model using Arnoldi iteration.
Definition: ArnoldiIterator.h:40
void setCalculateEigenVectors(bool calculateEigenVectors)
Definition: ArnoldiIterator.h:260
void setTolerance(double tolerance)
Definition: ArnoldiIterator.h:276
const std::complex< double > * getEigenValues() const
Definition: ArnoldiIterator.h:288
int getBasisSize() const
Definition: Model.h:292
void setMaxIterations(int maxIterations)
Definition: ArnoldiIterator.h:280
Mode
Definition: ArnoldiIterator.h:57
Base class for Solvers.
int getNumLanczosVectors() const
Definition: ArnoldiIterator.h:272
Definition: Matrix.h:33
Flexible physical index.
Definition: Index.h:69
Definition: ModelFactory.h:35
Solves Mx = b for x, where M is a SparseMatrix.
void setNumLanczosVectors(int numLanczosVectors)
Definition: ArnoldiIterator.h:268
Mode getMode() const
Definition: ArnoldiIterator.h:248
Solves Mx = b for x, where M is a SparseMatrix.
Definition: LUSolver.h:37
void setNumEigenValues(int numEigenValues)
Definition: ArnoldiIterator.h:252
Definition: Communicator.h:28
Container of Model related information.
Definition: Model.h:52
void setCentralValue(double centralValue)
Definition: ArnoldiIterator.h:284
const double getEigenValue(int state) const
Definition: ArnoldiIterator.h:292
int getNumEigenValues() const
Definition: ArnoldiIterator.h:256