BlockDiagonalizer.h

/* Copyright 2017 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_TBTK_SOLVER_BLOCK_DIAGONALIZER
#define COM_DAFER45_TBTK_SOLVER_BLOCK_DIAGONALIZER
 
#include "TBTK/BlockStructureDescriptor.h"
#include "TBTK/CArray.h"
#include "TBTK/Communicator.h"
#include "TBTK/Model.h"
#include "TBTK/Solver/Solver.h"
#include "TBTK/Timer.h"
 
#include <complex>
 
namespace TBTK{
namespace Solver{
 
class BlockDiagonalizer : public Solver, public Communicator{
public:
    class SelfConsistencyCallback{
    public:
        virtual bool selfConsistencyCallback(
            BlockDiagonalizer &blockDiagonalizer
        ) = 0;
    };
 
    BlockDiagonalizer();
 
    void setSelfConsistencyCallback(
        SelfConsistencyCallback &selfConsistencyCallback
    );
 
    void setMaxIterations(int maxIterations);
 
    void run();
 
    const double getEigenValue(int state);
 
    const double getEigenValue(const Index &blockIndex, int state);
 
    const std::complex<double> getAmplitude(int state, const Index &index);
 
    const std::complex<double> getAmplitude(
        const Index &blockIndex,
        int state,
        const Index &intraBlockIndex
    );
 
    unsigned int getFirstStateInBlock(const Index &index) const;
 
    unsigned int getLastStateInBlock(const Index &index) const;
 
    void setParallelExecution(bool parallelExecution);
private:
    CArray<std::complex<double>> hamiltonian;
 
    CArray<double> eigenValues;
 
    CArray<std::complex<double>> eigenVectors;
 
    BlockStructureDescriptor blockStructureDescriptor;
 
    std::vector<unsigned int> blockSizes;
 
    std::vector<unsigned int> blockOffsets;
 
    std::vector<unsigned int> eigenVectorSizes;
 
    std::vector<unsigned int> eigenVectorOffsets;
 
    int maxIterations;
 
    bool parallelExecution;
 
    SelfConsistencyCallback *selfConsistencyCallback;
 
    void init();
 
    void update();
 
    void solve();
};
 
inline void BlockDiagonalizer::setSelfConsistencyCallback(
    SelfConsistencyCallback &selfConsistencyCallback
){
    this->selfConsistencyCallback = &selfConsistencyCallback;
}
 
inline void BlockDiagonalizer::setMaxIterations(int maxIterations){
    this->maxIterations = maxIterations;
}
 
inline const std::complex<double> BlockDiagonalizer::getAmplitude(
    int state,
    const Index &index
){
    const Model &model = getModel();
    unsigned int block = blockStructureDescriptor.getBlockIndex(state);
    unsigned int offset = eigenVectorOffsets.at(block);
    unsigned int linearIndex = model.getBasisIndex(index);
    unsigned int firstStateInBlock
        = blockStructureDescriptor.getFirstStateInBlock(block);
    unsigned int lastStateInBlock = firstStateInBlock
        + blockStructureDescriptor.getNumStatesInBlock(block)-1;
    offset += (
        state - firstStateInBlock
    )*blockStructureDescriptor.getNumStatesInBlock(block);
    if(
        linearIndex >= firstStateInBlock
        && linearIndex <= lastStateInBlock
    ){
        return eigenVectors[
            offset + (linearIndex - firstStateInBlock)
        ];
    }
    else{
        return 0;
    }
}
 
inline const std::complex<double> BlockDiagonalizer::getAmplitude(
    const Index &blockIndex,
    int state,
    const Index &intraBlockIndex
){
    int firstStateInBlock = getModel().getHoppingAmplitudeSet(
    ).getFirstIndexInBlock(blockIndex);
    unsigned int block = blockStructureDescriptor.getBlockIndex(
        firstStateInBlock
    );
    TBTKAssert(
        state >= 0
        && state < (int)blockStructureDescriptor.getNumStatesInBlock(
            block
        ),
        "BlockDiagonalizer::getAmplitude()",
        "Out of bound error. The block with block Index "
        << blockIndex.toString() << " has "
        << blockStructureDescriptor.getNumStatesInBlock(block)
        << " states, but state " << state << " was requested.",
        ""
    );
    unsigned int offset = eigenVectorOffsets.at(block)
        + state*blockStructureDescriptor.getNumStatesInBlock(block);
    unsigned int linearIndex = getModel().getBasisIndex(
        Index(blockIndex, intraBlockIndex)
    );
 
    return eigenVectors[offset + (linearIndex - firstStateInBlock)];
}
 
inline const double BlockDiagonalizer::getEigenValue(int state){
    return eigenValues[state];
}
 
inline const double BlockDiagonalizer::getEigenValue(
    const Index &blockIndex,
    int state
){
    int offset = getModel().getHoppingAmplitudeSet().getFirstIndexInBlock(
        blockIndex
    );
 
    return eigenValues[offset + state];
}
 
inline unsigned int BlockDiagonalizer::getFirstStateInBlock(
    const Index &index
) const{
    unsigned int linearIndex = getModel().getBasisIndex(index);
    unsigned int block
        = blockStructureDescriptor.getBlockIndex(linearIndex);
 
    return blockStructureDescriptor.getFirstStateInBlock(block);
}
 
inline unsigned int BlockDiagonalizer::getLastStateInBlock(
    const Index &index
) const{
    unsigned int linearIndex = getModel().getBasisIndex(index);
    unsigned int block
        = blockStructureDescriptor.getBlockIndex(linearIndex);
 
    return getFirstStateInBlock(index)
        + blockStructureDescriptor.getNumStatesInBlock(block) - 1;
}
 
inline void BlockDiagonalizer::setParallelExecution(
    bool parallelExecution
){
    this->parallelExecution = parallelExecution;
}
 
};  //End of namespace Solver
};  //End of namespace TBTK
 
#endif