/*test_uBlas.cpp
 $Author: jlconlin $
 $Id: test_Arnoldi.cpp 309 2008-05-19 15:22:12Z jlconlin $
 $Revision: 309 $
 $Date: 2008-05-19 11:22:12 -0400 (Mon, 19 May 2008) $

 This file is used for testing/learning about Boost.org uBlas library.

 */

#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 basic(int);
void geev(int);
void hseqr(int);
template <typename T>
void Hessenberg(ublas::matrix<T, ublas::column_major>& );
template <typename T>
void Hessenberg(ublas::matrix<complex<T>, ublas::column_major>& );

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

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

}
void basic(int n){
    int m = 6;
    cout << "\n'Regular' matrix." << endl;
    ublas::matrix<double> A(m,n);
    for (unsigned i = 0; i < A.size1 (); ++ i)
        for (unsigned j = 0; j < A.size2 (); ++ j)
            A (i, j) = n * i + j;
    cout << "Here is my matrix A:\n" << A << endl;
    cout << "Here is a 2x2 submatrix of A:\n";
    ublas::matrix_range<ublas::matrix<double> > mr(A,
            ublas::range(1,3),ublas::range(2,4) );
    cout << mr << endl;

    cout << "\nColumn Major matrix." << endl;
    ublas::matrix<double, ublas::column_major> cm(m,n);
    for (unsigned i = 0; i < cm.size1 (); ++ i)
        for (unsigned j = 0; j < cm.size2 (); ++ j)
            cm (i, j) = n * i + j;
    cout << "Here is my matrix cm:\n" << cm << endl;

    cout << "\n\nRow Major rm." << endl;
    ublas::matrix<double, ublas::row_major> rm(m,n);
    for (unsigned i = 0; i < rm.size1 (); ++ i)
        for (unsigned j = 0; j < rm.size2 (); ++ j)
            rm (i, j) = n * i + j;
    cout << "Here is my matrix rm:\n" << rm << endl;

}

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<complex<double>, ublas::column_major> Z(n,n);

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

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

//  cout << "\nHSEQR for only eigenvalues.  Complex version" << endl;
//  ublas::matrix<complex<double>, ublas::column_major> G(n,n);
//  Hessenberg(G);
//  cout << "\nG:\n" << G << endl;
//  lapack::hseqr('E', G, values);
//  cout << "Done with hseqr." << endl;
//  cout << "\nG:\n" << G << endl;
//  cout << "\nvalues: " << values << endl;
//  cout << "\nZ:\n" << Z << endl;

//  cout << "\nGSEQR for eigenvalues and Schur vectors." << endl;
//  Hessenberg(G);
//  cout << "G:\n" << G << endl;
//  cout << "\nCalling hseqr." << endl;
//  lapack::hseqr('S', 'I', G, values, Z);
//  cout << "\nG:\n " << G << endl;
//  cout << "\nvalues: " << values << endl;
//  cout << "\nZ:\n " << 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;
            }
        }
    }
}

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

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;

}