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

From: karl.meerbergen_at_[hidden]
Date: 2008-06-10 07:09:46

Author: karlmeerbergen
Date: 2008-06-10 07:09:46 EDT (Tue, 10 Jun 2008)
New Revision: 46290
URL: http://svn.boost.org/trac/boost/changeset/46290

Log:
added hseqr

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

Added: sandbox/libs/numeric/bindings/lapack/test/hseqr.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/hseqr.cpp	2008-06-10 07:09:46 EDT (Tue, 10 Jun 2008)
@@ -0,0 +1,122 @@
+/*hseqr.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>
+#include <boost/numeric/bindings/lapack/trevc.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 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;
+    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', H, values);
+    cout << "\nH:\n" << H << endl;
+    cout << "\nvalues: " << values << endl;
+
+    cout << "\nHSEQR for eigenvalues and Schur vectors." << endl;
+    Hessenberg(H);
+    cout << "H:\n" << H << endl;
+    lapack::hseqr('S', 'I', H, values, Z);
+    cout << "\nH: " << H << endl;
+    cout << "\nvalues: " << values << endl;
+    cout << "\nZ: " << Z << endl;
+
+    cout << "\n==================================" << endl;
+    cout << "Recalculating original matrix..." << endl;
+    ublas::matrix<double, ublas::column_major> cH(n,n);
+    cH = ublas::prod( H, ublas::herm(Z) );
+    cH = ublas::prod( Z, cH );
+    cout << "'New' original matrix:\n" << cH << endl;
+    cout << "==================================" << 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 << "\nG:\n" << G << endl;
+    cout << "\nvalues: " << values << endl;
+
+    cout << "\nHSEQR for eigenvalues and Schur vectors." << endl;
+    Hessenberg(G);
+    cout << "G:\n" << G << endl;
+    ublas::matrix<complex<double>, ublas::column_major> cZ(n,n);
+    lapack::hseqr('S', 'I', G, values, cZ);
+    cout << "\nG:\n " << G << endl;
+    cout << "\nvalues: " << values << endl;
+    cout << "\nZ:\n " << Z << endl;
+
+    cout << "\n==================================" << endl;
+    cout << "Recalculating original matrix..." << endl;
+    ublas::matrix<complex<double>, ublas::column_major> origG(G);
+    origG = ublas::prod( G, ublas::herm(cZ) );
+    origG = ublas::prod( cZ, origG );
+    cout << "'New' original matrix:\n" << origG << endl;
+    cout << "==================================" << 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);
+            }
+        }
+    }
+}
+