MTK++ Latest version: 0.2.0

diagonalize.h
Go to the documentation of this file.
00001 
00033 #ifndef DIAGONALIZE_H
00034 #define DIAGONALIZE_H
00035 
00036 #include <boost/numeric/ublas/symmetric.hpp>
00037 #include <boost/numeric/ublas/matrix.hpp>
00038 #include <boost/numeric/ublas/io.hpp>
00039 #include "boost/numeric/bindings/lapack/syev.hpp"
00040 #include "boost/numeric/bindings/lapack/gesvd.hpp"
00041 #include "boost/numeric/bindings/lapack/gesdd.hpp"
00042 #include "boost/numeric/bindings/traits/ublas_matrix.hpp"
00043 #include "boost/numeric/bindings/traits/ublas_vector.hpp"
00044 
00045 namespace ublas  = boost::numeric::ublas;
00046 namespace lapack = boost::numeric::bindings::lapack;
00047 
00048 namespace MTKpp
00049 {
00050 
00065 inline int diagonalize(ublas::matrix<double, ublas::column_major>& eigenvectors,
00066                        ublas::vector<double>& eigenvalues) {
00067     int r = lapack::syev( 'V', 'U', eigenvectors, eigenvalues, lapack::minimal_workspace() );
00068     return r;
00069 }
00070 
00100 inline int svd(ublas::matrix<double, ublas::column_major>& H,
00101                ublas::vector<double>& s,
00102                ublas::matrix<double, ublas::column_major>& U,
00103                ublas::matrix<double, ublas::column_major>& VT) {
00104 
00105     //int r = lapack::gesvd('A', 'A', H, s, U, VT);
00106     int r = lapack::gesdd(H, s, U, VT);
00107     return r;
00108 }
00109 
00116 inline void eigenValueSort(ublas::matrix<double, ublas::column_major>& eigenvectors,
00117                            ublas::vector<double>& eigenvalues, int order) {
00118     int k;
00119     double p;
00120     int size = eigenvectors.size1();
00121 
00122     // Ascending
00123     if (!order) {
00124       for (int i = 0; i < size-1 ; ++i) {
00125         k = i;
00126         p = eigenvalues(i);
00127 
00128         for (int j = i+1; j < size; ++j) {
00129           if (eigenvalues(j) < p) {
00130             k = j;
00131             p = eigenvalues(j);
00132           }
00133         }
00134         if ( k != i ) {
00135           eigenvalues(k) = eigenvalues(i);
00136           eigenvalues(i) = p;
00137           for (int m = 0; m < size; ++m) {
00138             p = eigenvectors(m,i);
00139             eigenvectors(m,i) = eigenvectors(m,k);
00140             eigenvectors(m,k) = p;
00141           }
00142         }
00143       }
00144     }
00145     // Descending
00146     else {
00147       for (int i = 0; i < size-1 ; ++i) {
00148         k = i;
00149         p = eigenvalues(i);
00150 
00151         for (int j = i+1; j < size; ++j) {
00152           if (eigenvalues(j) > p) {
00153             k = j;
00154             p = eigenvalues(j);
00155           }
00156         }
00157         if ( k != i ) {
00158           eigenvalues(k) = eigenvalues(i);
00159           eigenvalues(i) = p;
00160           for (int m = 0; m < size; ++m) {
00161             p = eigenvectors(m,i);
00162             eigenvectors(m,i) = eigenvectors(m,k);
00163             eigenvectors(m,k) = p;
00164           }
00165         }
00166       }
00167     }
00168 }
00169 
00170 } // MTKpp namespace
00171 
00172 #endif // DIAGONALIZE_H

Generated on Fri Dec 23 2011 09:28:49 for MTK++ by Doxygen 1.7.5.1