[Boost-commit] svn:boost r46643 - sandbox/libs/numeric/bindings/lapack/test

From: karl.meerbergen_at_[hidden]
Date: 2008-06-24 10:25:57

Author: karlmeerbergen
Date: 2008-06-24 10:25:56 EDT (Tue, 24 Jun 2008)
New Revision: 46643
URL: http://svn.boost.org/trac/boost/changeset/46643

Log:
added test for geev

Added:
   sandbox/libs/numeric/bindings/lapack/test/ublas_geev.cpp   (contents, props changed)

Added: sandbox/libs/numeric/bindings/lapack/test/ublas_geev.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_geev.cpp	2008-06-24 10:25:56 EDT (Tue, 24 Jun 2008)
@@ -0,0 +1,66 @@
+#include<iostream>
+#include<cstdlib>
+#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 "utils.h"
+
+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);
+template <typename T>
+void Hessenberg(ublas::matrix<T, ublas::column_major>& );
+
+int main(){
+    cout << "I'm testing uBlas." << endl;
+
+    int n = 5;
+    geev(n);
+
+}
+void geev(int n){
+    cout << "\nCalculating eigenvalues using LAPACK's geev." << endl;
+    ublas::matrix<double, ublas::column_major> A(n,n);
+    Hessenberg(A);
+    print_m(A);
+
+    ublas::vector<complex<double> > values(n);
+    ublas::matrix<complex<double>, ublas::column_major>* Vectors_left = 0;
+    ublas::matrix<complex<double>, ublas::column_major> Vectors_right(n,n);
+
+    lapack::geev(A, values, Vectors_left, &Vectors_right, lapack::optimal_workspace());
+    print_v(values, "values"); cout << endl;
+    print_m(Vectors_right, "Vectors_right"); cout << endl;
+
+    Hessenberg(A);
+    cout << "A*x = l*x." << endl;
+    for( int i = 0; i < Vectors_right.size2(); ++i ){
+        ublas::vector<complex<double> > tmp(n);
+        tmp = ublas::prod( A, column(Vectors_right, i) );
+        cout << tmp - values(i)*column(Vectors_right,i) << 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;
+            }
+        }
+    }
+}
+