Program Listing for File Hamiltonian_CSR.hpp

↰ Return to documentation for file (src_ed/Hamiltonian/Hamiltonian_CSR.hpp)

/*
 Implementation of the Hamiltonian in CSR format (own implementation)
*/

#ifndef Hamiltonian_csr
#define Hamiltonian_csr

#include "Hamiltonian_base.hpp"
#include "CSR_hermitian.hpp"


template<typename HilbertField>
class Hamiltonian_CSR : public Hamiltonian<HilbertField>
{
    public:

        CSR_hermitian<HilbertField> H_csr;

        Hamiltonian_CSR(
            shared_ptr<model> the_model,
            const map<string, double> &value,
            sector _sec
        );
        void mult_add(vector<HilbertField> &x, vector<HilbertField> &y);
        void diag(vector<double> &d);
        //vector<shared_ptr<state<HilbertField>>> states(double& GS_energy);

    private:
        map<shared_ptr<HS_Hermitian_operator>, double> sparse_ops;
        void HS_ops_map(const map<string, double> &value);

};


template<typename HilbertField>
Hamiltonian_CSR<HilbertField>::Hamiltonian_CSR(
    shared_ptr<model> _the_model,
    const map<string, double> &value,
    sector _sec
) {
    std::lock_guard<std::mutex> ctor_lock(Hamiltonian_ctor_mutex());
    this->the_model = _the_model;
    this->sec = _sec;
    this->B = _the_model->provide_basis(_sec);
    this->dim = this->B->dim;
    if(this->dim == 0) return;

    HS_ops_map(value);

    int num=1;
    #ifdef _OPENMP
    num = omp_get_max_threads()/omp_get_num_threads();
    #endif
    map<index_pair,HilbertField> E;
    H_csr.diag.assign(this->dim, 0.0);
    if(global_bool("verb_Hilbert"))cout << "constructing the CSR Hamiltonian..." << endl;
    for(auto& h : sparse_ops){
        h.first->CSR_map(E, H_csr.diag, h.second);
    }
    size_t row = 0;
    size_t count=0;
    H_csr.Iptr.reserve(this->dim/2);
    H_csr.J.reserve(E.size());
    H_csr.v.reserve(E.size());
    H_csr.Iptr.push_back(0);
    for(auto &x : E) {
        if(x.first.r != row) {
            for(size_t i=row; i<x.first.r; i++) H_csr.Iptr.push_back(count);
            row = x.first.r;
        }
        H_csr.J.push_back(x.first.c);
        H_csr.v.push_back(x.second);
        count++;
    }
    H_csr.Iptr.push_back(count);
    this->H_ptr = &H_csr;
}


template<typename HilbertField>
void Hamiltonian_CSR<HilbertField>::mult_add(
    vector<HilbertField> &x,
    vector<HilbertField> &y
) {
    H_csr.apply(x, y); // applies the CSR matrix
}


template<typename HilbertField>
void Hamiltonian_CSR<HilbertField>::diag(vector<double> &d){
    for(auto& h : sparse_ops) h.first->diag(d, h.second);
}

template<typename HilbertField>
void Hamiltonian_CSR<HilbertField>::HS_ops_map(const map<string, double> &value)
{
    bool is_complex = false;
    if(typeid(HilbertField) == typeid(Complex)) is_complex = true;
    //create a vector of values keys to have random access iterator
    vector<string> keys;
    keys.reserve(value.size());
    for (auto& x : value) {
        keys.push_back(x.first);
    }

    //construct the Hamiltonian in parallel
    #pragma omp parallel for schedule(dynamic,1)
    for (auto& x : keys) {
        Hermitian_operator& op = *this->the_model->term.at(x);
        {
            std::lock_guard<std::mutex> lock(op.hs_op_mutex);
            if(op.HS_operator.find(this->sec) == op.HS_operator.end()){
                // cout << this->the_model->name+" : building "+op.name+" in "+to_string<sector>(this->sec)+'\n' << std::flush;
                op.HS_operator[this->sec] = op.build_HS_operator(this->sec, is_complex); // ***TEMPO***
            }
        }
    }

    keys.resize(0);
    //then add it to sparse_ops
    for(const auto& x : value){
        Hermitian_operator& op = *this->the_model->term.at(x.first);
        std::lock_guard<std::mutex> lock(op.hs_op_mutex);
        sparse_ops[op.HS_operator.at(this->sec)] = x.second;
    }
}


#endif