fn_eig.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // - Edmund Highcock (edmund dot highcock at merton dot ox dot ac dot uk)
00007 // 
00008 // This file is part of the Armadillo C++ library.
00009 // It is provided without any warranty of fitness
00010 // for any purpose. You can redistribute this file
00011 // and/or modify it under the terms of the GNU
00012 // Lesser General Public License (LGPL) as published
00013 // by the Free Software Foundation, either version 3
00014 // of the License or (at your option) any later version.
00015 // (see http://www.opensource.org/licenses for more info)
00016 
00017 
00018 //! \addtogroup fn_eig
00019 //! @{
00020 
00021 
00022 //
00023 // symmetric/hermitian matrices
00024 //
00025 
00026 
00027 //! Eigenvalues of real/complex symmetric/hermitian matrix X
00028 template<typename T1>
00029 inline
00030 void
00031 eig_sym(Col<typename T1::pod_type>& eigval, const Base<typename T1::elem_type,T1>& X, const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0)
00032   {
00033   arma_extra_debug_sigprint();
00034   
00035   typedef typename T1::elem_type eT;
00036   
00037   // unwrap_check not used as T1::elem_type and T1::pod_type may not be the same.
00038   // furthermore, it doesn't matter if A is an alias of S, as auxlib::eig() makes a copy of A
00039 
00040   const unwrap<T1> tmp(X.get_ref());
00041   const Mat<eT>& A = tmp.M;
00042 
00043   auxlib::eig_sym(eigval, A);
00044   }
00045 
00046 
00047 
00048 //! Eigenvalues of real/complex symmetric/hermitian matrix X
00049 template<typename T1>
00050 inline
00051 Col<typename T1::pod_type>
00052 eig_sym(const Base<typename T1::elem_type,T1>& X, const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0)
00053   {
00054   arma_extra_debug_sigprint();
00055   
00056   Col<typename T1::pod_type> out;
00057   eig_sym(out, X);
00058   
00059   return out;
00060   }
00061 
00062 
00063 //! Eigenvalues and eigenvectors of real/complex symmetric/hermitian matrix X
00064 template<typename T1> 
00065 inline
00066 void
00067 eig_sym
00068   (
00069   Col<typename T1::pod_type>& eigval,
00070   Mat<typename T1::elem_type>& eigvec,
00071   const Base<typename T1::elem_type,T1>& X,
00072   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00073   )
00074   {
00075   arma_extra_debug_sigprint();
00076 
00077   typedef typename T1::elem_type eT;
00078   
00079   const unwrap<T1> tmp(X.get_ref());
00080   const Mat<eT>& A = tmp.M;
00081   
00082   auxlib::eig_sym(eigval, eigvec, A);
00083   }
00084 
00085 
00086 
00087 //
00088 // general matrices
00089 //
00090 
00091 
00092 
00093 //! Eigenvalues and eigenvectors (both left and right) of general real/complex square matrix X
00094 template<typename T1>
00095 inline
00096 void
00097 eig_gen
00098   (
00099   Col< std::complex<typename T1::pod_type> >& eigval, 
00100   Mat<typename T1::elem_type>&                l_eigvec,
00101   Mat<typename T1::elem_type>&                r_eigvec,
00102   const Base<typename T1::elem_type,T1>&      X,
00103   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00104   )
00105   {
00106   arma_extra_debug_sigprint();
00107 
00108   typedef typename T1::elem_type eT;
00109   
00110   const unwrap<T1> tmp(X.get_ref());
00111   const Mat<eT>& A = tmp.M;
00112 
00113   auxlib::eig_gen(eigval, l_eigvec, r_eigvec, A, 'b');
00114   }
00115 
00116 
00117 
00118 //! Eigenvalues and eigenvectors of general real square matrix X.
00119 //! Optional argument 'side' specifies which eigenvectors should be computed:
00120 //! 'r' for right (default) and 'l' for left.
00121 template<typename eT, typename T1>
00122 inline
00123 void
00124 eig_gen
00125   (
00126   Col< std::complex<eT> >& eigval, 
00127   Mat< std::complex<eT> >& eigvec,
00128   const Base<eT, T1>& X, 
00129   const char side = 'r',
00130   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00131   )
00132   {
00133   arma_extra_debug_sigprint();
00134 
00135   //std::cout << "real" << std::endl;
00136 
00137   const unwrap<T1> tmp(X.get_ref());
00138   const Mat<eT>& A = tmp.M;
00139 
00140   Mat<eT> dummy_eigvec;
00141   Mat<eT> tmp_eigvec;
00142   
00143   switch(side)
00144     {
00145     case 'r':
00146       auxlib::eig_gen(eigval, dummy_eigvec, tmp_eigvec, A, side);
00147       break;
00148 
00149     case 'l':
00150       auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, A, side);
00151       break;
00152       
00153     default:
00154       arma_stop("eig_gen(): parameter 'side' is invalid");
00155       return;
00156     }
00157 
00158 
00159   const u32 n = A.n_rows;
00160 
00161   if(n > 0)
00162     {
00163     eigvec.set_size(n,n);
00164 
00165     for(u32 j=0; j<n; ++j)
00166       {
00167       if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) )
00168         {
00169         // eigvec.col(j)   = Mat< std::complex<eT> >( tmp_eigvec.col(j),  tmp_eigvec.col(j+1) );
00170         // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -tmp_eigvec.col(j+1) );
00171 
00172         for(u32 i=0; i<n; ++i)
00173           {
00174           eigvec.at(i,j)   = std::complex<eT>( tmp_eigvec.at(i,j),  tmp_eigvec.at(i,j+1) );
00175           eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eigvec.at(i,j+1) );
00176           }
00177 
00178         ++j;
00179         }
00180       else
00181         {
00182         // eigvec.col(i) = tmp_eigvec.col(i);
00183 
00184         for(u32 i=0; i<n; ++i)
00185           {
00186           eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0));
00187           }
00188 
00189         }
00190       }
00191     }
00192 
00193   }
00194 
00195 
00196 
00197 //! Eigenvalues and eigenvectors of general complex square matrix X
00198 //! Optional argument 'side' specifies which eigenvectors should be computed:
00199 //! 'r' for right (default) and 'l' for left.
00200 template<typename T, typename T1>
00201 inline
00202 void
00203 eig_gen
00204   (
00205   Col< std::complex<T> >& eigval, 
00206   Mat< std::complex<T> >& eigvec,
00207   const Base<std::complex<T>, T1>& X, 
00208   const char side = 'r',
00209   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00210   )
00211   {
00212   arma_extra_debug_sigprint();
00213   //std::cout << "complex" << std::endl;
00214 
00215   typedef typename std::complex<T> eT;
00216 
00217   const unwrap<T1> tmp(X.get_ref());
00218   const Mat<eT>& A = tmp.M;
00219 
00220   Mat<eT> dummy_eigvec;
00221   
00222   switch(side)
00223     {
00224     case 'r':
00225       auxlib::eig_gen(eigval, dummy_eigvec, eigvec, A, side);
00226       break;
00227 
00228     case 'l':
00229       auxlib::eig_gen(eigval, eigvec, dummy_eigvec, A, side);
00230       break;
00231       
00232     default:
00233       arma_stop("eig_gen(): parameter 'side' is invalid");
00234     }
00235   }
00236 
00237 
00238 
00239 //! @}
00240