//Purpose:
//  Solve tridiagonal system with matrix RHS using
//  gbsv.
//References:
//  http://abel.ee.ucla.edu/cvxopt/userguide/lapack.html
//  http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/lle/functn_gbsv.html
//
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/bindings/lapack/driver/gbsv.hpp>
#include <boost/numeric/bindings/std/vector.hpp>
#include <boost/numeric/bindings/ublas/banded.hpp>
#include <boost/numeric/bindings/ublas/matrix.hpp>
#include <boost/numeric/bindings/bandwidth.hpp>
#include <vector>
#include <stdexcept>

namespace ublas = boost::numeric::ublas;
namespace lapack = boost::numeric::bindings::lapack;
namespace bindings = boost::numeric::bindings;

// solves the equation Ax = B, and puts the solution in B
// A is mutated by this routine
template <typename MatrA, typename MatrB>
void InPlaceSolve(MatrA& a, MatrB& b)
{
  std::vector<fortran_int_t> piv(a.size1());
  int ret = lapack::gbsv(a,piv,b);
  if (ret < 0) {
    //CStdString err;
    //err.Format("banded::Solve: argument %d in DGBTRF had an illegal value", -ret);
    //throw RuntimeError(err);
    throw std::runtime_error("banded::Solve: argument %d in GBSV had an illegal value");
  }
  if (ret > 0) {
    //CStdString err;
    //err.Format("banded::Solve: the (%d,%d) diagonal element is 0 after DGBTRF", ret, ret);
    //throw RuntimeError(err);
    throw std::runtime_error("banded::Solve: the (%d,%d) diagonal element is 0 after GBSV");
  }
}

  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:
          static
        unsigned const
      kl=1
      ;
          static
        unsigned const
      ku=1
      ;
      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()
        , kl
        , kl+ku
          //this, although counterintuitive,
          //is what do_typename function in
          //the source code referenced in the 
          //  {ublas_gbsv
          //comment does.
        )
      , 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  ]; 
          }
      }
};

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();
        unsigned const
      m=2;
        ublas::matrix
        < value_t
        , ublas::column_major
        >
      x(n,m)
        ;
      for( unsigned j=0; j<m; ++j)
      {
          value_t pow10=std::pow(10,value_t(j));
          x(0,j)=  7.0*pow10;
          x(1,j)=  6.0*pow10;
          x(2,j)=- 4.0*pow10;
          x(3,j)=-15.0*pow10;
      }
      std::cout<<"b=\n"<<x<<"\n";
    #if 1
      {
          ublas::banded_matrix<value_t>&au=a;
          InPlaceSolve(au, x);
          std::cout<<"x=\n"<<x<<"\n";
          std::cout<<"a(mutated)=\n"<<au<<"\n";
      }
    #endif      
  }   
  
int main()
{
    test_tridiag();
    return 0;
}