#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/blas.hpp>
#include <iomanip>

  template
  < typename Value
  >
  struct
tridiag_ublas
: boost::numeric::ublas::banded_matrix<Value>
{
        typedef
      boost::numeric::ublas::banded_matrix<Value>
    super_t
      ;
        typedef
      std::vector<Value>
    diag_t
      ;
 private:      
    tridiag_ublas()
      ;
      unsigned
    my_size
      ;
 public:
      unsigned
    size()const
      {
          return my_size;
      }
    tridiag_ublas
      ( diag_t const&a_u
      , diag_t const&a_d
      , diag_t const&a_l
      )
      : super_t( a_d.size(), a_d.size(), 1, 1)
      , my_size(a_d.size())
      {
          super_t&s=*this;
          for(unsigned i=0; i<my_size; ++i)
          {
              if(0<i        )s( i, i-1)=a_l[i-1];
                             s( i, i  )=a_d[i  ];
              if(i<my_size-1)s( i, i+1)=a_u[i  ]; 
          }
      }
};

  template
  < typename Value
  >
  std::ostream&
operator<<
  ( std::ostream& sout
  , tridiag_ublas<Value>const& a
  )
  {
      std::size_t const size=a.size();
      sout<<"{ ";
      for(unsigned i=0; i<size; ++i)
      {
          if(0<i) sout<<", ";
          sout<<"{ ";
          for(unsigned j=0; j<size; ++j)
          {
              if(0<j) sout<<", ";
              unsigned const w=5;
              int const del_index=i-j;
              //sout<<"del_index="<<del_index<<"\n";
              if(std::abs(del_index)<=1) sout<<std::setw(w)<<a( i, j);
              else                       sout<<std::setw(w)<<"_";
          }
          sout<<"}\n";
      }
      sout<<"}\n";
      return sout;
  }          

  template
  < typename Value
  >
  std::ostream&
operator<<
  ( std::ostream& sout
  , boost::numeric::ublas::matrix<Value>const& a
  )
  {
      std::size_t const size1=a.size1();
      std::size_t const size2=a.size2();
      sout<<"{ ";
      for(unsigned i1=0; i1<size1; ++i1)
      {
          if(0<i1) sout<<", ";
          sout<<"{ ";
          for(unsigned i2=0; i2<size2; ++i2)
          {
              if(0<i2) sout<<", ";
              unsigned const w=5;
              sout<<std::setw(w)<<a( i1, i2);
          }
          sout<<"}\n";
      }
      sout<<"}\n";
      return sout;
  }          

void test_tridiag()
 /**@brief
  *  Demonstrate ublas tridiagonal solve with
  *  input, originally, from an array_dyn.
  *  The test problem is that
  *  from Example 3.2E on page 136 of [MarLop1991].
  **@references
  *  [MarLop1991]
  *    Maron, Melvin J.
  *    Lopez, Robert J.
  *    _Numerical Analysis: A Practical Approach_
  *    Wadsworth(1991)
  */
  {
      typedef double value_t;
      typedef tridiag_ublas<value_t> trid_t;
      typedef std::vector<value_t> vv_t;
        trid_t 
      a ( vv_t( { -3.0, -2.0, 1.0} ) //upper
        , vv_t( {  5.0,  4.0, 3.0, 5.0} ) //diagonal
        , vv_t( {  1.0, -1.0, 2.0} ) //lower
        );
      std::cout<<"a.size()="<<a.size()<<"\na=\n"<<a<<"\n";
        unsigned const 
      n=a.size();
          typedef
        boost::numeric::ublas::matrix<value_t>
      array_t
        ;
        array_t
      r0(n,1)
        ;
      r0(0,0)=  7.0;
      r0(1,0)=  6.0;
      r0(2,0)=- 4.0;
      r0(3,0)=-15.0;
      std::cout<<"r0=\n"<<r0<<"\n";
    #if 1
      {
          array_t x(r0.size1(), r0.size2());
          value_t t=1.0;
          boost::numeric::ublas::blas_3::tsm(x,t,r0,t);
          std::cout<<"x=\n"<<x<<"\n";
      }
    #endif      
  }   
  
int main()
{
    test_tridiag();
    return 0;
}