DMRGBlock.hpp

#ifndef __DMRG_BLOCK_HPP__
#define __DMRG_BLOCK_HPP__

#include <petscmat.h>
#include "QuantumNumbers.hpp"
#include "MiscTools.hpp"
#include "kron.hpp"
#include "linalg_tools.hpp"
#include <string>
#include <map>

typedef enum
{
    OpSm = -1,  
    OpSz = 0,   
    OpSp = +1,  
    OpEye= +2   
} Op_t;

static const std::map<Op_t, std::string> OpString =
{
    {OpSm, "Sm"},
    {OpSz, "Sz"},
    {OpSp, "Sp"}
};

#define OpToCStr(OP)    ((OpString.find(OP)->second).c_str())

#define OpToStr(OP)     (OpString.find(OP)->second)

typedef enum
{
    SideLeft = 0,   
    SideRight = 1   
} Side_t;

typedef enum
{
    SpinOneHalf  = 102, 
    SpinOne      = 101, 
    SpinNull     = -2   
} Spin_t;

static const std::map<std::string,Spin_t> SpinTypes = {
    {"1/2",         SpinOneHalf },
    {"1",           SpinOne     }
};

static const std::map<Spin_t,std::string> SpinTypesString = {
    {SpinOneHalf,   "1/2"       },
    {SpinOne,       "1"         }
};

static const std::vector<Op_t> BasicOpTypes = { OpSz, OpSp };
static const std::vector<Side_t> SideTypes = { SideLeft, SideRight };

namespace Block {

    class SpinBase
    {

    private:

        MPI_Comm        mpi_comm = PETSC_COMM_SELF;

        PetscMPIInt     mpi_rank;

        PetscMPIInt     mpi_size;

        Spin_t          spin_type = SpinOneHalf;

        std::string     spin_type_str = "1/2";

        PetscInt        _loc_dim = 2;

        std::vector<PetscScalar>    _loc_qn_list = {+0.5, -0.5};

        std::vector<PetscInt>       _loc_qn_size = {1, 1};

        PetscBool init = PETSC_FALSE;

        PetscBool init_once = PETSC_FALSE;

        PetscBool mpi_init = PETSC_FALSE;

        PetscBool verbose = PETSC_FALSE;

        PetscBool log_io = PETSC_FALSE;

        PetscInt num_sites;

        PetscInt num_states;

        PetscBool init_Sm = PETSC_FALSE;

        std::vector<Mat> SzData;

        std::vector<Mat> SpData;

        std::vector<Mat> SmData;

        PetscBool init_save = PETSC_FALSE;

        PetscBool disk_set = PETSC_FALSE;

        PetscBool saved = PETSC_FALSE;

        PetscBool retrieved = PETSC_FALSE;

        std::string save_dir;

        std::string read_dir;

        std::string write_dir;

        PetscInt num_reads = 0;

        PetscErrorCode SaveOperator(
            const std::string& OpName,
            const size_t& isite,
            Mat& Op,
            const MPI_Comm& comm_in
            );

        PetscErrorCode Retrieve_NoChecks();

        /* Number of subcommunicators to be used when performing the rotation. */
        PetscInt nsubcomm = 1;

        /* Tells the subcommunicator to use for rotating this site */
        std::vector< PetscMPIInt > site_color;

        enum RotMethod { mmmmult=0, matptap=1 };

        RotMethod rot_method = mmmmult;

    protected:

        PetscBool MPIInitialized() const { return mpi_init; }

    public:

        virtual PetscInt loc_dim() const
        {
            return _loc_dim;
        }

        virtual std::vector<PetscScalar> loc_qn_list() const
        {
            return _loc_qn_list;
        }

        virtual std::vector<PetscInt> loc_qn_size() const
        {
            return _loc_qn_size;
        }

        virtual PetscErrorCode MatSpinSzCreate(Mat& Sz);

        virtual PetscErrorCode MatSpinSpCreate(Mat& Sp);

        PetscErrorCode Initialize(
            const MPI_Comm& comm_in       
        );

        PetscErrorCode Initialize(
            const MPI_Comm& comm_in,      
            const PetscInt& num_sites_in, 
            const PetscInt& num_states_in,
            const PetscBool& init_ops = PETSC_TRUE
            );

        PetscErrorCode Initialize(
            const MPI_Comm& comm_in,                    
            const PetscInt& num_sites_in,               
            const std::vector<PetscReal>& qn_list_in,   
            const std::vector<PetscInt>& qn_size_in,    
            const PetscBool& init_ops = PETSC_TRUE
            );

        PetscErrorCode Initialize(
            const PetscInt& num_sites_in,   
            const QuantumNumbers& qn_in     
            );

        PetscErrorCode InitializeFromDisk(
            const MPI_Comm& comm_in,        
            const std::string& block_path   
            );

        PetscErrorCode InitializeSave(
            const std::string& save_dir_in  
            );

        PetscErrorCode SetDiskStorage(
            const std::string& read_dir_in,     
            const std::string& write_dir_in     
            );

        std::string SaveDir() const { return save_dir; }

        PetscBool SaveInitialized() const { return init_save; }

        PetscErrorCode SaveAndDestroy();

        PetscErrorCode SaveBlockInfo();

        PetscErrorCode RetrieveOperator(
            const std::string& OpName,
            const size_t& isite,
            Mat& Op,
            const MPI_Comm& comm_in=MPI_COMM_NULL
            );

        PetscErrorCode Retrieve();

        PetscErrorCode EnsureSaved();

        PetscErrorCode EnsureRetrieved();

        PetscErrorCode Destroy();

        QuantumNumbers Magnetization;

        PetscErrorCode CheckOperatorArray(
            const Op_t& OpType  
            ) const;

        PetscBool Initialized() const { return init; }

        PetscBool Saved() const { return saved; }

        MPI_Comm MPIComm() const { return mpi_comm; }

        PetscInt NumSites() const {return num_sites; }

        PetscInt NumStates() const {return num_states; }

        Mat     H = nullptr;

        Mat Sz(const PetscInt& Isite) const {
            if(Isite >= num_sites) throw std::runtime_error("Attempted to access non-existent site.");
            return SzData[Isite];
        }

        Mat Sp(const PetscInt& Isite) const {
            if(Isite >= num_sites) throw std::runtime_error("Attempted to access non-existent site.");
            return SpData[Isite];
        }

        Mat Sm(const PetscInt& Isite) const {
            if(Isite >= num_sites) throw std::runtime_error("Attempted to access non-existent site.");
            if(!init_Sm) throw std::runtime_error("Sm matrices were not initialized on this block.");
            return SmData[Isite];
        }

        const std::vector<Mat>& Sz() const { return SzData; }

        const std::vector<Mat>& Sp() const { return SpData; }

        const std::vector<Mat>& Sm() const { return SmData; }

        PetscErrorCode CheckOperators() const;

        PetscErrorCode CheckSectors() const;

        PetscErrorCode MatOpGetNNZs(
            const Op_t& OpType,
            const PetscInt& isite,
            std::vector<PetscInt>& nnzs
            ) const;

        PetscErrorCode MatGetNNZs(
            const Mat& matin,
            std::vector<PetscInt>& nnzs
            ) const;

        PetscErrorCode MatOpCheckOperatorBlocks(
            const Op_t& op_t,       
            const PetscInt& isite   
            ) const;

        PetscErrorCode MatCheckOperatorBlocks(
            const Op_t& op_t,       
            const Mat& matin        
            ) const;

        PetscErrorCode CheckOperatorBlocks() const;

        PetscErrorCode GetOperatorBlocks(Op_t Operator);

        PetscErrorCode CreateSm();

        PetscErrorCode DestroySm();

        PetscErrorCode RotateOperators(
            SpinBase& Source,           
            const Mat& RotMatT          
            );

        PetscErrorCode AssembleOperators();
    };

}

#endif // __DMRG_BLOCK_HPP__