//References:
//  http://abel.ee.ucla.edu/cvxopt/userguide/lapack.html
//
#include <boost/numeric/ublas/io.hpp>
//{ublas_gbsv
//  2011-11-16.
//    The following code, up to the comment starting with }ublas_gbsv
//    was copied from:
//       http://svn.boost.org/svn/boost/sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gbsv.cpp
#include <boost/numeric/bindings/std/vector.hpp>
#include <boost/numeric/bindings/ublas/banded.hpp>
#include <boost/numeric/bindings/ublas/vector.hpp>
#include <boost/numeric/bindings/bandwidth.hpp>
#include <boost/numeric/bindings/lapack/computational/gbtrf.hpp>
#include <boost/numeric/bindings/lapack/computational/gbtrs.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)
{
  // if the matrix has kl lower and ku upper diagonals, then we should have
  // allocated kl lower and kl+ku upper diagonals
  fortran_int_t const kl = bindings::bandwidth_lower (a);
  fortran_int_t const ku = bindings::bandwidth_upper (a) - kl;
  std::vector<fortran_int_t> piv(a.size1());
  int ret = lapack::gbtrf(/*kl, ku, */a, piv);
  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 DGBTRF 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 DGBTRF");
  }

  ret = lapack::gbtrs(/*kl, ku, */a, piv, b);
  if (ret < 0) {
    //CStdString err;
    //err.Format("banded::Solve: argument %d in DGBTRS had an illegal value", -ret);
    //throw RuntimeError(err);
    throw std::runtime_error("banded::Solve: argument %d in DGBTRS had an illegal value");
  }
}

//}ublas_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();
        boost::numeric::ublas::vector<value_t>
      x(n)
        ;
      x(0)=  7.0;
      x(1)=  6.0;
      x(2)=- 4.0;
      x(3)=-15.0;
      std::cout<<"b=\n"<<x<<"\n";
    #if 1
      {
          boost::numeric::ublas::banded_matrix<value_t>&au=a;
          InPlaceSolve(au, x);
          std::cout<<"x=\n"<<x<<"\n";
      }
    #endif      
  }   
  
int main()
{
    test_tridiag();
    return 0;
}