SuperMatrix.h

Go to the documentation of this file.

#ifndef SUPERMATRIX
#define SUPERMATRIX
 
#include <set>
#include <boost/multi_array.hpp>
 
#include "MemCalc.h" // from TOOLS
//include "DmrgTypedefs.h"
//include "DmrgHamiltonianTerms.h"
#include "tensors/SiteOperator.h"
 
 
template<typename Symmetry, typename Scalar=double>
class SuperMatrix
{
typedef SparseMatrix<double,ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> MatrixType;
typedef SiteOperator<Symmetry,Scalar> OperatorType;
public:
    
    OperatorType &operator() (size_t i, size_t j)       {return data[i][j];} // write
    OperatorType  operator() (size_t i, size_t j) const {return data[i][j];} // read
    
    void set (size_t Daux1, size_t Daux2, size_t D)
    {
        N_rows = Daux1;
        N_cols = Daux2;
        data.resize(boost::extents[Daux1][Daux2]);
        innerResize(D);
        setZero();
    }
    
    void setRowVector (size_t Daux, size_t D)
    {
        set(1,Daux,D);
    }
    
    void setColVector (size_t Daux, size_t D)
    {
        set(Daux,1,D);
    }
    
    void setMatrix (size_t Daux, size_t D)
    {
        set(Daux,Daux,D);
    }
    
    SuperMatrix<Symmetry,Scalar> row (size_t i)
    {
        SuperMatrix<Symmetry,Scalar> Mout;
        Mout.setRowVector(N_cols,D());        // instead of auxdim
        for (size_t j=0; j<N_cols; ++j)        //   Bug? N_rows
        {
            Mout(0,j) = data[i][j];
        }
        return Mout;
    }
    
    SuperMatrix<Symmetry,Scalar> col (size_t i)
    {
        SuperMatrix<Symmetry,Scalar> Mout;
        Mout.setColVector(N_rows,D());          // instead of auxdim
        for (size_t j=0; j<N_rows; ++j)         //  Bug? N_cols
        {
            Mout(j,0) = data[j][i];
        }
        return Mout;
    }
    
    void setZero()
    {
        for (size_t i=0; i<N_rows; ++i)
        for (size_t j=0; j<N_cols; ++j)
        {
            data[i][j].data.setZero();
        }
    }
    
    inline size_t rows() const {return N_rows;}
    inline size_t cols() const {return N_cols;}
    
    inline size_t auxdim() const
    {
        assert(N_rows == N_cols and "auxdim called although SuperMatrix is not quadratic");
        return std::max(N_rows,N_cols);
    }
    
    double memory (MEMUNIT memunit=GB) const
    {
        double out = 0.;
        
        for (size_t i=0; i<N_rows; ++i)
        for (size_t j=0; j<N_cols; ++j)
        {
            return out += calc_memory<Scalar>(data[i][j].data, memunit);
        }
        return out;
    }
    
    size_t D() const
    {
        size_t Dres = data[0][0].data.rows();
        for (size_t i=0; i<N_rows; ++i)
        for (size_t j=0; j<N_cols; ++j)
        {
            assert(data[i][j].data.rows() == Dres);
            assert(data[i][j].data.cols() == Dres);
        }
        return Dres;
    }
    
    std::vector<typename Symmetry::qType> calc_qOp() const
    {
        std::set<typename Symmetry::qType> qOps;
        for (size_t i=0; i<N_rows; ++i)
        for (size_t j=0; j<N_cols; ++j)
        {
            qOps.insert(data[i][j].Q);
        }
        
        std::vector<typename Symmetry::qType> out(qOps.size());
        copy(qOps.begin(), qOps.end(), out.begin());
        
        return out;
    }
    
private:
    
    size_t N_rows;
    size_t N_cols;
    
    boost::multi_array<OperatorType,2> data;
    
    void innerResize (size_t D)
    {
        for (size_t i=0; i<N_rows; ++i)
        for (size_t j=0; j<N_cols; ++j)
        {
            data[i][j].data.resize(D,D);
        }
    }
};
 
template<typename Symmetry, typename Scalar>
SuperMatrix<Symmetry,Scalar> tensor_product (const SuperMatrix<Symmetry,Scalar> &M1, const SuperMatrix<Symmetry,Scalar> &M2)
{
    assert(M1.D() == M2.D());
    
    SuperMatrix<Symmetry,Scalar> Mout;
    
    if (M1.rows() == 1)
    {
        Mout.setRowVector(M1.cols()*M2.cols(), M1.D());             // instead of auxdim
    }
    else if (M1.cols() == 1)
    {
        Mout.setColVector(M1.rows()*M2.rows(), M1.D());             // instead of auxdim
    }
    else
    {
        Mout.set(M1.rows()*M2.rows(), M1.cols()*M2.cols(), M1.D()); // instead of auxdim
    }
    
    for (size_t r1=0; r1<M1.rows(); ++r1)
    for (size_t c1=0; c1<M1.cols(); ++c1)
    for (size_t r2=0; r2<M2.rows(); ++r2)
    for (size_t c2=0; c2<M2.cols(); ++c2)
    {
        Mout(r1*M2.rows()+r2, c1*M2.cols()+c2).data = M1(r1,c1).data * M2(r2,c2).data;
        auto Qsum = Symmetry::reduceSilent(M1(r1,c1).Q, M2(r2,c2).Q);
        // should be only one term, non-Abelian symmetries have a different code
        assert(Qsum.size() == 1 and "tensor_product of SuperMatrices called with wrong symmetry!");
        Mout(r1*M2.rows()+r2, c1*M2.cols()+c2).Q = Qsum[0];
    }
    
    return Mout;
}
 
template<typename Symmetry, typename Scalar>
SuperMatrix<Symmetry,Scalar> directSum (const SuperMatrix<Symmetry,Scalar> &M1, const SuperMatrix<Symmetry,Scalar> &M2)
{
    SuperMatrix<Symmetry,Scalar> Mout;
    size_t R;
    size_t C;
    
    if (M1.rows()==1 and M2.rows()==1)
    {
        Mout.setRowVector(M1.cols()+M2.cols(), M1.D());         // instead of auxdim
        R = 0;
        C = M1.cols();
    }
    else if (M1.cols()==1 and M2.cols()==1)
    {
        Mout.setColVector(M1.rows()+M2.rows(), M1.D());     // instead of auxdim
        R = M1.rows();
        C = 0;
    }
    else
    {
        Mout.set(M1.rows()+M2.rows(), M1.cols()+M2.cols(), M1.D());            // instead of auxdim
        R = M1.rows();
        C = M1.cols();
    }
    
    for (size_t r1=0; r1<M1.rows(); ++r1)
    for (size_t c1=0; c1<M1.cols(); ++c1)
    {
        Mout(r1,c1) = M1(r1,c1);
    }
    
    for (size_t r2=0; r2<M2.rows(); ++r2)
    for (size_t c2=0; c2<M2.cols(); ++c2)
    {
        Mout(R+r2,C+c2) = M2(r2,c2);
    }
    
    return Mout;
}
 
template<typename Symmetry, typename Scalar>
std::ostream &operator<< (std::ostream& os, const SuperMatrix<Symmetry,Scalar> &M)
{
    os << std::showpos << std::setprecision(1) << std::fixed;
    for (int i=0; i<M.rows(); ++i)
    {
        for (int n=0; n<M.D(); ++n)
        {
            for (int j=0; j<M.cols(); ++j)
            {
                for (int m=0; m<M.D(); ++m)
                {
                    
                    os << Matrix<Scalar,Dynamic,Dynamic>(M(i,j).data)(n,m);
                }
                os << " ";
            }
            if (n!=M.D()-1) {os << std::endl;}
        }
        if (i!=M.rows()-1) {os << std::endl << std::endl;}
    }
    os << noshowpos;
    return os;
}
 
/*template<typename Symmetry, typename Scalar>
SuperMatrix<Symmetry, Scalar> Generator (const HamiltonianTerms<Symmetry, Scalar> &Terms)
{
    typedef SiteOperator<Symmetry,Scalar> OperatorType;
    size_t Daux = 2 + Terms.tight.size() + 2*Terms.nextn.size();
    
    std::vector<OperatorType> col;
    std::vector<OperatorType> row;
    size_t locdim;
    
    if (Terms.local.size()>0)
    {
        locdim = std::get<1>(Terms.local[0]).data.rows();
    }
    else if (Terms.tight.size()>0)
    {
        locdim = std::get<1>(Terms.tight[0]).data.rows();
    }
    else
    {
        locdim = std::get<1>(Terms.nextn[0]).data.rows();
    }
    
    OperatorType Id(Matrix<Scalar,Dynamic,Dynamic>::Identity(locdim,locdim).sparseView(),Symmetry::qvacuum());
    OperatorType Zero(SparseMatrix<Scalar>(locdim,locdim),Symmetry::qvacuum());
    
    // last row (except corner element)
    for (size_t i=0; i<Terms.nextn.size(); ++i)
    {
        row.push_back(Zero);
    }
    for (int i=0; i<Terms.tight.size(); ++i)
    {
        row.push_back(std::get<0>(Terms.tight[i]) * std::get<1>(Terms.tight[i]));
    }
    for (int i=0; i<Terms.nextn.size(); ++i)
    {
        row.push_back(std::get<0>(Terms.nextn[i]) * std::get<1>(Terms.nextn[i]));
    }
    row.push_back(Id);
    
    // first col (except corner element)
    col.push_back(Id);
    for (int i=0; i<Terms.nextn.size(); ++i)
    {
        col.push_back(std::get<2>(Terms.nextn[i]));
    }
    for (int i=0; i<Terms.tight.size(); ++i)
    {
        col.push_back(std::get<2>(Terms.tight[i]));
    }
    for (size_t i=0; i<Terms.nextn.size(); ++i)
    {
        col.push_back(Zero);
    }
    
    SuperMatrix<Symmetry,Scalar> Gout;
    Gout.setMatrix(Daux,locdim);
    Gout.setZero();
    
    for (size_t i=0; i<Daux-1; ++i)
    {
        Gout(i,0)        = col[i];
        Gout(Daux-1,i+1) = row[i];
    }
    
    // corner element : local interaction
    for (int i=0; i<Terms.local.size(); ++i)
    {
        Gout(Daux-1,0) += std::get<0>(Terms.local[i]) * std::get<1>(Terms.local[i]);
    }
    
    // nearest-neighbour transfer
    if (Terms.nextn.size() != 0)
    {
        for (size_t i=0; i<Terms.nextn.size(); ++i)
        {
            Gout(Daux-1-Terms.nextn.size()+i,1+i) = std::get<3>(Terms.nextn[i]);
        }
    }
    
    return Gout;
}*/
 
#endif