Program Listing for File mcf_set.hpp

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

#ifndef mcf_set_h
#define mcf_set_h

#include "Green_function_set.hpp"
#include "matrix_continued_fraction.hpp"
#include "global_parameter.hpp"

template<typename T>
struct mcf_set : Green_function_set
{
    vector<matrix_continued_fraction<T>> e;
    vector<matrix_continued_fraction<T>> h;
    vector<matrix_continued_fraction<T>> combined;

    mcf_set(shared_ptr<symmetry_group> _group, int mixing)
        : Green_function_set(_group, mixing)
    {
        e.resize(group->g);
        h.resize(group->g);
        combined.resize(group->g);
    }


    void build_combined()
    {
        if(!global_bool("combine_mcf")) return;
        for(size_t r = 0; r < group->g; ++r){
            bool has_e = e[r].floors() > 0;
            bool has_h = h[r].floors() > 0;
            if(has_e && has_h){
                int M0 = e[r].floors() + h[r].floors();
                combined[r] = combine_via_lanczos(e[r], h[r], M0);
            }
            // single-sector blocks: combined stays empty; Green_function handles them.
        }
    }

    // Virtual method implementations
    void Green_function(const Complex &z, block_matrix<Complex> &G) override;
    void integrated_Green_function(block_matrix<Complex> &G) override;
    void write_hdf5(H5::Group& grp) override;
    void read_hdf5(H5::Group& grp) override;
};


//==============================================================================
// Inline implementations


template<typename T>
inline void mcf_set<T>::Green_function(const Complex &z, block_matrix<Complex> &G)
{
    if(global_bool("combine_mcf")){
        for(size_t r = 0; r < group->g; ++r)
            if(combined[r].floors() > 0)
                G.block[r] += combined[r].evaluate(z);
        return;
    }
    for(size_t r = 0; r < group->g; ++r){
        if(e[r].floors() > 0)
            G.block[r] += e[r].evaluate(z);
        if(h[r].floors() > 0){
            matrix<Complex> Gh = h[r].evaluate(z);
            Gh.transpose();
            G.block[r] += Gh;
        }
    }
}


template<typename T>
inline void mcf_set<T>::integrated_Green_function(block_matrix<Complex> &G)
{
    for(size_t r = 0; r < group->g; ++r) {
        if(h[r].floors() == 0 || h[r].p == 0) continue;
        int p = h[r].p;
        matrix<Complex> Wh = h[r].W;
        Wh.hermitian_conjugate();
        matrix<Complex> M(p);
        M.product(Wh, h[r].W);  // M = W^H * W
        M.transpose();           // (W^H W)^T: cf. G(b,a) += M(a,b) convention
        G.block[r] += M;
    }
}


template<typename T>
inline void mcf_set<T>::write_hdf5(H5::Group& grp)
{
    h5_write_attr(grp, "nblocks", (int)group->g);
    for(size_t r = 0; r < group->g; ++r){
        H5::Group bg = grp.createGroup("block_" + to_string(r));
        H5::Group eg = bg.createGroup("e");
        h5_write_mcf(eg, e[r]);
        H5::Group hg = bg.createGroup("h");
        h5_write_mcf(hg, h[r]);
    }
}


template<typename T>
inline void mcf_set<T>::read_hdf5(H5::Group& grp)
{
    int nblocks = h5_read_attr_int(grp, "nblocks");
    e.resize(nblocks);
    h.resize(nblocks);
    for(int r = 0; r < nblocks; ++r){
        H5::Group bg = grp.openGroup("block_" + to_string(r));
        H5::Group eg = bg.openGroup("e");
        h5_read_mcf(eg, e[r]);
        H5::Group hg = bg.openGroup("h");
        h5_read_mcf(hg, h[r]);
    }
}


#endif /* mcf_set_h */