/*test_uBlas.cpp
 */

#include<iostream>
#include<complex>
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/matrix_proxy.hpp"
#include "boost/numeric/ublas/io.hpp"
#include "boost/numeric/bindings/traits/ublas_matrix.hpp"
#include "boost/numeric/bindings/traits/ublas_vector.hpp"
#include "boost/numeric/bindings/lapack/geev.hpp"
#include "boost/numeric/bindings/lapack/hseqr.hpp"

using std::cout;
using std::endl;
using std::vector;
using std::complex;

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

void geev(int);
void hseqr(int);
template <typename T>
void Hessenberg(ublas::matrix<T, ublas::column_major>& );

int main(){
    cout << "I'm testing uBlas." << endl;

    int n = 5;
//  geev(n);
    hseqr(n);

}

void hseqr(int n){
    cout << "\nCalculating eigenvalues using LAPACK's hseqr." << endl;
    ublas::matrix<double, ublas::column_major> H(n,n);
    Hessenberg(H);
    cout << "\nUpper Hessenberg matrix H:\n" << H << endl;

    ublas::vector<complex<double> > values(n);
    ublas::matrix<double, ublas::column_major> Z(n,n);

    cout << "\nHSEQR for only eigenvalues." << endl;
    lapack::hseqr('E', 'N', H, values, Z);
    cout << "Done with hseqr." << endl;
    cout << "\nH: " << H << endl;
    cout << "\nvalues: " << values << endl;
    cout << "\nZ: " << Z << endl;

    cout << "\nHSEQR for eigenvalues and Schur vectors." << endl;
    Hessenberg(H);
    cout << "H:\n" << H << endl;
    cout << "\nCalling hseqr." << endl;
    lapack::hseqr('S', 'N', H, values, Z);
    cout << "\nH: " << H << endl;
    cout << "\nvalues: " << values << endl;
    cout << "\nZ: " << Z << endl;


}

template <typename T>
void Hessenberg(ublas::matrix<T, ublas::column_major>& H){
    T k = 1;
    for( unsigned int i = 0; i < H.size1(); ++i ){
        for( unsigned int j = i; j <= H.size2(); ++j ){
            if( j > 0 ){
                H(i,j-1) = k;
                k += 1;
            }
        }
    }
}



void geev(int n){
    cout << "\nCalculating eigenvalues using LAPACK's geev." << endl;
    ublas::matrix<double, ublas::column_major> A(n,n);
    for( unsigned int i = 0; i < n; ++i ) A(i,i) = i+1;
    A(2,1) = 1;
    cout << "A:\n" << A << endl;

    ublas::vector<complex<double> > values(n);
    ublas::matrix<double, ublas::column_major> Vectors_left(n,n);
    ublas::matrix<double, ublas::column_major> Vectors_right(n,n);

    lapack::geev(A, values, &Vectors_left, &Vectors_right, lapack::optimal_workspace());
    cout << "eigenvalues: " << values << endl;
    cout << "eigenvectors_right: " << Vectors_right<< endl;
    cout << "eigenvectors_left: " << Vectors_left<< endl;

}