#include <algorithm>
#include <boost/numeric/bindings/lapack/computational/gbcon.hpp>
#include <boost/numeric/bindings/lapack/computational/gbtrf.hpp>
#include <boost/numeric/bindings/lapack/auxiliary/lange.hpp>
#include <boost/numeric/bindings/ublas.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
#include <boost/numeric/ublas/operation/num_columns.hpp>
#include <boost/numeric/ublas/operation/num_rows.hpp>
#include <boost/numeric/ublas/traits.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <cstddef>
#include <iostream>
#include <string>


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


int main()
{
	typedef double value_type;
	typedef ublas::banded_matrix<value_type, ublas::column_major> matrix_type;
	typedef ublas::matrix_traits<matrix_type>::size_type size_type; 
	typedef matrix_type::array_type array_type; 
	typedef ublas::type_traits<value_type>::real_type result_type;
	typedef ublas::banded_matrix<value_type,ublas::row_major,array_type> auxiliary_matrix_type;

	const size_type nr = 4;
	const size_type nc = 4;
	const size_type kl = 1;
	const size_type ku = 2;
	const size_type n = std::min(nr, nc);

    matrix_type A(nr,nc,kl,ku);
    A(0,0) = -0.23; A(0,1) = 2.54; A(0,2) = -3.66;         /* 0 */ // 1st row
    A(1,0) = -6.98; A(1,1) = 2.46; A(1,2) = -2.73; A(1,3) = -2.13; // 2nd row
            /* 0 */ A(2,1) = 2.56; A(2,2) =  2.46; A(2,3) =  4.07; // 3rd row
            /* 0 */        /* 0 */ A(3,2) = -4.78; A(3,3) = -3.82; // 4th row


	// Compute LUP decomposition

	auxiliary_matrix_type AB(A, kl, kl+ku);

	std::cout << "AB = " << AB << std::endl;

	ublas::vector< ::fortran_int_t > ipiv(n);

	bindings::lapack::gbtrf(AB, ipiv);

	std::cout << "GBTRF -> AB = " << AB << std::endl;
	std::cout << "GBTRF -> IPIV = " << ipiv << std::endl;

	// Compute matrix 1-norm

	result_type norm(0);

	//NOTE: bindings::lapack::lange is actually broken

	norm = ublas::norm_1(A);

	// Estimate reciprocal matrix condition number 

	result_type rcond(0);

	bindings::lapack::gbcon('O', AB, ipiv, norm, rcond);

	std::cout << "GBCON -> RCOND = " << rcond << std::endl;
}