TBTK
Need a break? Support the development by playing Polarity Puzzles
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/CArray.h"
27 #include "TBTK/Model.h"
28 #include "TBTK/Solver/LUSolver.h"
29 #include "TBTK/Solver/Solver.h"
30 
31 #include <complex>
32 
33 namespace TBTK{
34 namespace Solver{
35 
60 class ArnoldiIterator : public Solver, public Communicator{
61 public:
64 
74  enum class Mode {Normal, ShiftAndInvert};
75 
79  void setMode(Mode mode);
80 
84  Mode getMode() const;
85 
89  void setNumEigenValues(int numEigenValues);
90 
94  int getNumEigenValues() const;
95 
100  void setCalculateEigenVectors(bool calculateEigenVectors);
101 
105  bool getCalculateEigenVectors() const;
106 
111  void setNumLanczosVectors(int numLanczosVectors);
112 
117  int getNumLanczosVectors() const;
118 
124  void setTolerance(double tolerance);
125 
131  void setMaxIterations(int maxIterations);
132 
146  void setCentralValue(double centralValue);
147 
149  void run();
150 
153 
155  const double getEigenValue(int state) const;
156 
161  const std::complex<double> getAmplitude(int state, const Index &index);
162 private:
164  Mode mode;
165 
167  int numEigenValues;
168 
171  bool calculateEigenVectors;
172 
175  int numLanczosVectors;
176 
179  double shift;
180 
183  double tolerance;
184 
187  int maxIterations;
188 
190  CArray<std::complex<double>> residuals;
191 
193  CArray<std::complex<double>> eigenValues;
194 
196  CArray<std::complex<double>> eigenVectors;
197 
199 
201  LUSolver luSolver;
202 
205  void initNormal();
206 
209  void initShiftAndInvert();
210 
212  void arnoldiLoop();
213 
215  void checkZnaupdInfo(int info) const;
216 
218  bool executeReverseCommunicationMessage(
219  int ido,
220  int basisSize,
221  double *workd,
222  int *ipntr,
223  Matrix<double> &b
224  );
225 
227  bool executeReverseCommunicationMessage(
228  int ido,
229  int basisSize,
230  std::complex<double> *workd,
231  int *ipntr,
232  Matrix<std::complex<double>> &b
233  );
234 
235  void checkZneupdIerr(int ierr) const;
236 
239  void sort();
240 
242  void mergeSortSplit(
243  std::complex<double> *dataIn,
244  std::complex<double> *dataOut,
245  int *orderIn,
246  int *orderOut,
247  int first,
248  int end
249  );
250 
252  void mergeSortMerge(
253  std::complex<double> *dataIn,
254  std::complex<double> *dataOut,
255  int *orderIn,
256  int *orderOut,
257  int first,
258  int middle,
259  int end
260  );
261 };
262 
263 inline void ArnoldiIterator::setMode(Mode mode){
264  this->mode = mode;
265 }
266 
268  return mode;
269 }
270 
271 inline void ArnoldiIterator::setNumEigenValues(int numEigenValues){
272  this->numEigenValues = numEigenValues;
273 }
274 
276  return numEigenValues;
277 }
278 
279 inline void ArnoldiIterator::setCalculateEigenVectors(bool calculateEigenVectors){
280  this->calculateEigenVectors = calculateEigenVectors;
281 }
282 
284  return calculateEigenVectors;
285 }
286 
287 inline void ArnoldiIterator::setNumLanczosVectors(int numLanczosVectors){
288  this->numLanczosVectors = numLanczosVectors;
289 }
290 
292  return numLanczosVectors;
293 }
294 
295 inline void ArnoldiIterator::setTolerance(double tolerance){
296  this->tolerance = tolerance;
297 }
298 
299 inline void ArnoldiIterator::setMaxIterations(int maxIterations){
300  this->maxIterations = maxIterations;
301 }
302 
303 inline void ArnoldiIterator::setCentralValue(double centralValue){
304  shift = centralValue;
305 }
306 
308  return eigenValues;
309 }
310 
311 inline const double ArnoldiIterator::getEigenValue(int state) const{
312  return real(eigenValues[state]);
313 }
314 
315 inline const std::complex<double> ArnoldiIterator::getAmplitude(
316  int state,
317  const Index &index
318 ){
319  const Model &model = getModel();
320  return eigenVectors[model.getBasisSize()*state + model.getBasisIndex(index)];
321 }
322 
323 }; //End of namespace Solver
324 }; //End of namesapce TBTK
325 
326 #endif
Base class for Solvers.
Definition: Solver.h:42
Model & getModel()
Definition: Solver.h:73
int getBasisIndex(const Index &index) const
Definition: Model.h:331
const std::complex< double > getAmplitude(int state, const Index &index)
Definition: ArnoldiIterator.h:315
void setMode(Mode mode)
Definition: ArnoldiIterator.h:263
Container of Model related information.
Mode getMode() const
Definition: ArnoldiIterator.h:267
Container for a C style array.
bool getCalculateEigenVectors() const
Definition: ArnoldiIterator.h:283
int getNumEigenValues() const
Definition: ArnoldiIterator.h:275
Solves a Model using Arnoldi iteration.
Definition: ArnoldiIterator.h:60
void setCalculateEigenVectors(bool calculateEigenVectors)
Definition: ArnoldiIterator.h:279
void setTolerance(double tolerance)
Definition: ArnoldiIterator.h:295
void setMaxIterations(int maxIterations)
Definition: ArnoldiIterator.h:299
const CArray< std::complex< double > > & getEigenValues() const
Definition: ArnoldiIterator.h:307
Mode
Definition: ArnoldiIterator.h:74
Base class for Solvers.
Definition: Matrix.h:33
Physical index.
Definition: Index.h:44
Container for a C style array.
Definition: CArray.h:44
Definition: Boolean.h:32
Definition: SparseMatrix.h:35
Solves Mx = b for x, where M is a SparseMatrix.
void setNumLanczosVectors(int numLanczosVectors)
Definition: ArnoldiIterator.h:287
Solves Mx = b for x, where M is a SparseMatrix.
Definition: LUSolver.h:37
int getBasisSize() const
Definition: Model.h:311
void setNumEigenValues(int numEigenValues)
Definition: ArnoldiIterator.h:271
Base class for classes that can communicate their status during execution.
Definition: Communicator.h:56
Container of Model related information.
Definition: Model.h:57
int getNumLanczosVectors() const
Definition: ArnoldiIterator.h:291
void setCentralValue(double centralValue)
Definition: ArnoldiIterator.h:303
const double getEigenValue(int state) const
Definition: ArnoldiIterator.h:311