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{
41 public:
44 
46  virtual ~ArnoldiIterator();
47 
57  enum class Mode {Normal, ShiftAndInvert};
58 
60  void setMode(Mode mode);
61 
65  void setNumEigenValues(int numEigenValues);
66 
68  int getNumEigenValues() const;
69 
71  void setCalculateEigenVectors(bool calculateEigenVectors);
72 
74  bool getCalculateEigenVectors() const;
75 
78  void setNumLanczosVectors(int numLanczosVectors);
79 
81  void setTolerance(double tolerance);
82 
85  void setMaxIterations(int maxIterations);
86 
88  void setCentralValue(double centralValue);
89 
91  void run();
92 
94  const std::complex<double>* getEigenValues() const;
95 
97  const double getEigenValue(int state) const;
98 
103  const std::complex<double> getAmplitude(int state, const Index &index);
104 private:
106  Mode mode;
107 
109  int numEigenValues;
110 
113  bool calculateEigenVectors;
114 
117  int numLanczosVectors;
118 
121  double shift;
122 
125  double tolerance;
126 
129  int maxIterations;
130 
132  std::complex<double> *residuals;
133 
135  std::complex<double> *eigenValues;
136 
138  std::complex<double> *eigenVectors;
139 
141  LUSolver luSolver;
142 
145  void initNormal();
146 
149  void initShiftAndInvert();
150 
152  void arnoldiLoop();
153 
155  void checkZnaupdInfo(int info) const;
156 
158  bool executeReverseCommunicationMessage(
159  int ido,
160  int basisSize,
161  double *workd,
162  int *ipntr,
163  Matrix<double> &b
164  );
165 
167  bool executeReverseCommunicationMessage(
168  int ido,
169  int basisSize,
170  std::complex<double> *workd,
171  int *ipntr,
172  Matrix<std::complex<double>> &b
173  );
174 
175  void checkZneupdIerr(int ierr) const;
176 
179  void sort();
180 
182  void mergeSortSplit(
183  std::complex<double> *dataIn,
184  std::complex<double> *dataOut,
185  int *orderIn,
186  int *orderOut,
187  int first,
188  int end
189  );
190 
192  void mergeSortMerge(
193  std::complex<double> *dataIn,
194  std::complex<double> *dataOut,
195  int *orderIn,
196  int *orderOut,
197  int first,
198  int middle,
199  int end
200  );
201 };
202 
203 inline void ArnoldiIterator::setMode(Mode mode){
204  this->mode = mode;
205 }
206 
207 inline void ArnoldiIterator::setNumEigenValues(int numEigenValues){
208  this->numEigenValues = numEigenValues;
209 }
210 
212  return numEigenValues;
213 }
214 
215 inline void ArnoldiIterator::setCalculateEigenVectors(bool calculateEigenVectors){
216  this->calculateEigenVectors = calculateEigenVectors;
217 }
218 
220  return calculateEigenVectors;
221 }
222 
223 inline void ArnoldiIterator::setNumLanczosVectors(int numLanczosVectors){
224  this->numLanczosVectors = numLanczosVectors;
225 }
226 
227 inline void ArnoldiIterator::setTolerance(double tolerance){
228  this->tolerance = tolerance;
229 }
230 
231 inline void ArnoldiIterator::setMaxIterations(int maxIterations){
232  this->maxIterations = maxIterations;
233 }
234 
235 inline void ArnoldiIterator::setCentralValue(double centralValue){
236  shift = centralValue;
237 }
238 
239 inline const std::complex<double>* ArnoldiIterator::getEigenValues() const{
240  return eigenValues;
241 }
242 
243 inline const double ArnoldiIterator::getEigenValue(int state) const{
244  return real(eigenValues[state]);
245 }
246 
247 inline const std::complex<double> ArnoldiIterator::getAmplitude(
248  int state,
249  const Index &index
250 ){
251  const Model &model = getModel();
252  return eigenVectors[model.getBasisSize()*state + model.getBasisIndex(index)];
253 }
254 
255 }; //End of namespace Solver
256 }; //End of namesapce TBTK
257 
258 #endif
Base class for Solvers.
Definition: Solver.h:32
Model & getModel()
Definition: Solver.h:57
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: ArnoldiIterator.h:247
bool getCalculateEigenVectors() const
Definition: ArnoldiIterator.h:219
void setMode(Mode mode)
Definition: ArnoldiIterator.h:203
Container of Model related information.
int getBasisIndex(const Index &index) const
Definition: Model.h:214
Solves a Model using Arnoldi iteration.
Definition: ArnoldiIterator.h:40
void setCalculateEigenVectors(bool calculateEigenVectors)
Definition: ArnoldiIterator.h:215
void setTolerance(double tolerance)
Definition: ArnoldiIterator.h:227
const std::complex< double > * getEigenValues() const
Definition: ArnoldiIterator.h:239
int getBasisSize() const
Definition: Model.h:210
void setMaxIterations(int maxIterations)
Definition: ArnoldiIterator.h:231
Mode
Definition: ArnoldiIterator.h:57
Base class for Solvers.
Definition: Matrix.h:33
Flexible physical index.
Definition: Index.h:70
Definition: ModelFactory.h:35
Solves Hx = y for x, where H is a SparseMatrix.
void setNumLanczosVectors(int numLanczosVectors)
Definition: ArnoldiIterator.h:223
Solves Hx = y for x, where H is a SparseMatrix.
Definition: LUSolver.h:37
void setNumEigenValues(int numEigenValues)
Definition: ArnoldiIterator.h:207
Container of Model related information.
Definition: Model.h:52
void setCentralValue(double centralValue)
Definition: ArnoldiIterator.h:235
const double getEigenValue(int state) const
Definition: ArnoldiIterator.h:243
int getNumEigenValues() const
Definition: ArnoldiIterator.h:211