tri_diag_mat.h

#pragma once
#include <cgv/math/vec.h>
#include <cgv/math/mat.h>
 
namespace cgv {
    namespace math {
 
 
template<typename T>
class tri_diag_mat
{   
protected:
    vec<T> _data;
 
    unsigned _dim;
 
public:
    tri_diag_mat()
    {
        _dim=0;
    }
 
    tri_diag_mat(unsigned n)
    {
        resize(n);
    }
 
 
    tri_diag_mat(unsigned n,const T& c)
    {
        resize(n);
        fill(c);
    }
 
    tri_diag_mat(const tri_diag_mat<T>& m)  
    {
        _data = m._data;
        _dim = m._dim;
    }
 
    tri_diag_mat<T>& operator = (const tri_diag_mat<T>& m) 
    {
        _data = m._data;
        _dim  = m._dim;
        return *this;
    }
 
    operator const mat<T>() const
    {
        mat<T> m;
        m.zeros(_dim,_dim);
 
        for(int b = -1;b <= 1;b++)
            for(unsigned i =0; i < _dim;i++)
            {
                if(b== -1 && i==0)
                    continue;
                if(b== 1 && i==_dim-1)
                    continue;
 
                m(i,b+i)=operator()(i,b);
            }
        
        return m;
    }
 
    void identity()
    {
        
        for(unsigned i = 0; i < 3*_dim; i++)
        {
            if( i < _dim || i > 2*_dim)
                _data(i)=0;
            else
                _data(i)=1;
        }
 
    }
 
    operator T*()
    {
        return (T*)_data;
    }
 
    operator const T*() const
    {
        return (const T*)_data;
    }
 
    void fill(const T& c)
    {
        _data.fill(c);
        _data(0) = 0;       
        _data(3*_dim-1) =0;
    }
 
    bool is_square()
    {
        return true;
    }
 
    //transpose
    void transpose()
    {
 
    }
 
    void zeros()
    {
        _data.zeros(); 
    }
 
    void zeros(unsigned n)
    {
        resize(n);
        zeros();
    }
 
 
    void resize(unsigned n)
    {
        _dim=n;
        _data.resize(3*n);
        _data(0) = 0;       
        _data(3*n-1) =0;
    }
 
    unsigned nrows()
    {
        return _dim;
    }
        
    unsigned ncols()
    {
        return _dim;
    }
 
    unsigned size()
    {
        return 3*_dim;
    }
 
    T& operator()(unsigned r,int  b)
    {
        assert(b >= -1 && b <=1);   
        return _data[ (b+1)*_dim + r];
    } 
 
    const T& operator()(unsigned r,int  b) const
    {
        assert(b >= -1 && b <=1);   
        return _data[ (b+1)*_dim + r];
    } 
 
    
    vec<T> band(int b) const
    {   
        return _data.sub_vec((b+1)*_dim,_dim);
    }
    
 
    vec<T> operator*(const vec<T>& v) const
    {
        vec<T> r;
        r.zeros(v.dim());
                
        for(unsigned i = 0; i < v.dim(); i++)
            for(int b = -1; b <= 1;b++)
            {
                if((i+b >= 0) && (i + b<  v.dim()))
                    r(i) += operator()(i,b)*v(i+b);
            }
 
        return r;
    }
 
 
 
};
 
    }
}