ArnoldiIterator.h

/* Copyright 2016 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef COM_DAFER45_SOLVER_ARNOLDI_ITERATOR
#define COM_DAFER45_SOLVER_ARNOLDI_ITERATOR
 
#include "TBTK/CArray.h"
#include "TBTK/Model.h"
#include "TBTK/Solver/LUSolver.h"
#include "TBTK/Solver/Solver.h"
 
#include <complex>
 
namespace TBTK{
namespace Solver{
 
class ArnoldiIterator : public Solver, public Communicator{
public:
    ArnoldiIterator();
 
    enum class Mode {Normal, ShiftAndInvert};
 
    void setMode(Mode mode);
 
    Mode getMode() const;
 
    void setNumEigenValues(int numEigenValues);
 
    int getNumEigenValues() const;
 
    void setCalculateEigenVectors(bool calculateEigenVectors);
 
    bool getCalculateEigenVectors() const;
 
    void setNumLanczosVectors(int numLanczosVectors);
 
    int getNumLanczosVectors() const;
 
    void setTolerance(double tolerance);
 
    void setMaxIterations(int maxIterations);
 
    void setCentralValue(double centralValue);
 
    void run();
 
    const CArray<std::complex<double>>& getEigenValues() const;
 
    const double getEigenValue(int state) const;
 
    const std::complex<double> getAmplitude(int state, const Index &index);
private:
    Mode mode;
 
    int numEigenValues;
 
    bool calculateEigenVectors;
 
    int numLanczosVectors;
 
    double shift;
 
    double tolerance;
 
    int maxIterations;
 
    CArray<std::complex<double>> residuals;
 
    CArray<std::complex<double>> eigenValues;
 
    CArray<std::complex<double>> eigenVectors;
 
    SparseMatrix<std::complex<double>> matrix;
 
    LUSolver luSolver;
 
    void initNormal();
 
    void initShiftAndInvert();
 
    void arnoldiLoop();
 
    void checkZnaupdInfo(int info) const;
 
    bool executeReverseCommunicationMessage(
        int ido,
        int basisSize,
        double *workd,
        int *ipntr,
        Matrix<double> &b
    );
 
    bool executeReverseCommunicationMessage(
        int ido,
        int basisSize,
        std::complex<double> *workd,
        int *ipntr,
        Matrix<std::complex<double>> &b
    );
 
    void checkZneupdIerr(int ierr) const;
 
    void sort();
 
    void mergeSortSplit(
        std::complex<double> *dataIn,
        std::complex<double> *dataOut,
        int *orderIn,
        int *orderOut,
        int first,
        int end
    );
 
    void mergeSortMerge(
        std::complex<double> *dataIn,
        std::complex<double> *dataOut,
        int *orderIn,
        int *orderOut,
        int first,
        int middle,
        int end
    );
};
 
inline void ArnoldiIterator::setMode(Mode mode){
    this->mode = mode;
}
 
inline ArnoldiIterator::Mode ArnoldiIterator::getMode() const{
    return mode;
}
 
inline void ArnoldiIterator::setNumEigenValues(int numEigenValues){
    this->numEigenValues = numEigenValues;
}
 
inline int ArnoldiIterator::getNumEigenValues() const{
    return numEigenValues;
}
 
inline void ArnoldiIterator::setCalculateEigenVectors(bool calculateEigenVectors){
    this->calculateEigenVectors = calculateEigenVectors;
}
 
inline bool ArnoldiIterator::getCalculateEigenVectors() const{
    return calculateEigenVectors;
}
 
inline void ArnoldiIterator::setNumLanczosVectors(int numLanczosVectors){
    this->numLanczosVectors = numLanczosVectors;
}
 
inline int ArnoldiIterator::getNumLanczosVectors() const{
    return numLanczosVectors;
}
 
inline void ArnoldiIterator::setTolerance(double tolerance){
    this->tolerance = tolerance;
}
 
inline void ArnoldiIterator::setMaxIterations(int maxIterations){
    this->maxIterations = maxIterations;
}
 
inline void ArnoldiIterator::setCentralValue(double centralValue){
    shift = centralValue;
}
 
inline const CArray<std::complex<double>>& ArnoldiIterator::getEigenValues() const{
    return eigenValues;
}
 
inline const double ArnoldiIterator::getEigenValue(int state) const{
    return real(eigenValues[state]);
}
 
inline const std::complex<double> ArnoldiIterator::getAmplitude(
    int state,
    const Index &index
){
    const Model &model = getModel();
    return eigenVectors[model.getBasisSize()*state + model.getBasisIndex(index)];
}
 
};  //End of namespace Solver
};  //End of namesapce TBTK
 
#endif