Program Listing for File block_matrix.hpp

↰ Return to documentation for file (src_util/block_matrix.hpp)

#ifndef block_matrix_h
#define block_matrix_h

#include <iostream>
#include <vector>

#include "matrix.hpp"

template<typename T>
struct block_matrix
{
    size_t r;
    size_t c;
    vector<matrix<T>> block;
    vector<size_t> cumul_r;
    vector<size_t> cumul_c;

    block_matrix(): r(0),c(0)
    {
        set_size();
    }

    block_matrix(size_t n)
    {
        block.assign(n, matrix<T>());
    }

    block_matrix(vector<size_t> dim)
    {
    block.assign(dim.size(), matrix<T>());
    for(size_t i=0; i<dim.size(); i++) block[i].set_size(dim[i]);
    set_size();
    }

    block_matrix(vector<int> dim)
    {
    block.assign(dim.size(), matrix<T>());
    for(size_t i=0; i<dim.size(); i++) block[i].set_size(dim[i]);
    set_size();
    }

  block_matrix(vector<matrix<T> > &_block) : block(_block)
    {
    set_size();
    }

    block_matrix& operator=(const block_matrix<T>& A)
    {
        r = A.r;
        c = A.c;
        cumul_r = A.cumul_r;
        cumul_c = A.cumul_c;
        block = A.block;
        return(*this);
    }


    void set_size(){
        cumul_r.clear();
        cumul_r.reserve(block.size());
        r = 0;
        for(size_t i=0; i<block.size(); ++i){
            cumul_r.push_back(r);
            r += block[i].r;
        }
        cumul_c.clear();
        cumul_c.reserve(block.size());
        c = 0;
        for(size_t i=0; i<block.size(); ++i){
            cumul_c.push_back(c);
            c += block[i].c;
        }
    }

    inline T operator()(const size_t &i, const size_t &j)const{
        size_t I= block.size()-1; while(cumul_r[I]>i) I--;
        size_t J= block.size()-1; while(cumul_c[J]>i) J--;
        if(I != J) return(0.0);
        else return (block[I])(i-cumul_r[I],j-cumul_c[I]);
    }

    matrix<T> build_matrix(){
    matrix<T> A(r,c);
        for(size_t i=0; i<block.size(); i++){
            for(size_t j=0; j<block[i].r; j++){
                for(size_t k=0; k<block[i].c; k++){
                    A(cumul_r[i]+j,cumul_c[i]+k) = block[i](j,k);
                }
            }
        }
    return A;
    }


    template <typename S>
    block_matrix<T>& operator-=(const block_matrix<S> &A){
        for(size_t i=0; i<block.size(); ++i) block[i].v -= A.block[i].v;
        return *this;
    }

    template <typename S>
    block_matrix<T>& operator+=(const block_matrix<S> &A){
        for(size_t i=0; i<block.size(); ++i) block[i].v += A.block[i].v;
        return *this;
    }

    template <typename S>
    void add(matrix<S> &Y, T z) const
    {
        for(size_t ib=0; ib<block.size(); ++ib){ // loop over sets of rows
            size_t br = block[ib].r;
            size_t bc = block[ib].c;
            size_t R = cumul_r[ib];
            size_t C = cumul_c[ib];
            for(size_t i=0; i<br; ++i){
                for(size_t j=0; j<bc; ++j){
                    Y(R+i,C+j) += z*block[ib](i,j);
                }
            }
        }
    }

    void add(T d){
        for(size_t i=0; i<block.size(); ++i) block[i].add(d);
    }


    void clear(){
        for(size_t i=0; i<block.size(); ++i) block[i].clear();
    }

    double diff_sq(const block_matrix<T> &A)
    {
        double z=0.0;
        for(size_t i = 0; i < block.size(); ++i) z += block[i].diff_sq(*A.block[i]);
        return(z);
    }

    void inverse()
    {
        for(size_t i = 0; i < block.size(); ++i) block[i].inverse();
    }


    T trace(){
        T z(0.0);
        for(size_t i = 0; i < block.size(); ++i) z += block[i].trace();
        return z;
    }

    block_matrix<T>& operator+=(const double a){for(size_t i=0; i<block.size(); ++i) block[i].add(a); return *this;}
    block_matrix<T>& operator+=(const Complex a){for(size_t i=0; i<block.size(); ++i) block[i].add(a); return *this;}
    block_matrix<T>& operator-=(const double a){for(size_t i=0; i<block.size(); ++i) block[i].add(-a); return *this;}
    block_matrix<T>& operator-=(const Complex a){for(size_t i=0; i<block.size(); ++i) block[i].add(-a); return *this;}



    template <typename S, typename U>
    void mult_left(const matrix<S> &X, matrix<U> &Y) const
    {
        for(size_t ib=0; ib<block.size(); ++ib){ // loop over sets of rows of Y
            for(size_t i=0; i<block[ib].r; ++i){
                for(size_t j=0; j<X.c; ++j){
                    U z(0.0);
                    for(size_t k=0; k<block[ib].c; ++k){
                        z += block[ib](i,k)*X(k+cumul_c[ib],j);
                    }
                    Y(i+cumul_r[ib],j) = z;
                }
            }
        }
    }

    template <typename S, typename U>
    void mult_right(const matrix<S> &X, matrix<U> &Y) const
    {
        for(size_t ib=0; ib<block.size(); ++ib){ // loop over sets of rows of Y
            for(size_t i=0; i<block[ib].c; ++i){
                for(size_t j=0; j<X.r; ++j){
                    U z(0.0);
                    for(size_t k=0; k<block[ib].r; ++k){
                        z += block[ib](k,i)*X(j,k+cumul_r[ib]);
                    }
                    Y(j,i+cumul_c[ib]) = z;
                }
            }
        }
    }




    template<typename S, typename U>
    void simil(matrix<S>  &A, const matrix<U> &B)
    {
        A.zero();
        for(size_t ij=0; ij<block.size(); ij++){
            size_t J_r = (int)cumul_r[ij];
            size_t J_c = (int)cumul_c[ij];
            for(size_t j = 0; j < block[ij].c; j++){
                for(size_t ii=0; ii<block.size(); ii++){
                    size_t I_r = (int)cumul_r[ii];
                    size_t I_c = (int)cumul_c[ii];
                    for(size_t i = 0; i < block[ii].c; ++i){
                        S z = 0.0;
                        for(size_t k = 0; k < block[ii].r; ++k){
                            S zz = 0.0;
                            for(size_t l = 0; l < block[ij].r; ++l){
                                zz += block[ij](l,j)*B(k+I_r,l+J_r);
                            }
                            z += zz*conjugate(block[ii](k,i));
                        }
                        A(i+I_c,j+J_c) += z;
                    }
                }
            }
        }
    }




    template<typename S, typename U>
    void simil_conjugate(matrix<S>  &A, const matrix<U> &B)
    {
        A.zero();
        for(size_t ij=0; ij<block.size(); ij++)
        {
            size_t J_r = (int)cumul_r[ij];
            size_t J_c = (int)cumul_c[ij];
            for(size_t j = 0; j < block[ij].c; j ++)
            {
                for(size_t ii=0; ii<block.size(); ii++)
                {
                    size_t I_r = (int)cumul_r[ii];
                    size_t I_c = (int)cumul_c[ii];
                    for(size_t i = 0; i < block[ii].c; ++i)
                    {
                        S z = 0.0;
                        for(size_t k = 0; k < block[ii].r; ++k)
                        {
                            S zz = 0.0;
                            for(size_t l = 0; l < block[ij].r; ++l)
                            {
                                zz += conjugate(block[ij](j,l))*B(k+I_r,l+J_r);
                            }
                            z += zz*conjugate(block[ii](i,k));
                        }
                        A(i+I_c,j+J_c) += z;
                    }
                }
            }
        }
    }



};

template <typename T>
std::ostream & operator<<(std::ostream &flux, block_matrix<T> &A){
    for(size_t i=0; i<A.block.size(); ++i) flux << A.block[i];
    return flux;
}


#endif