fn_eig.hpp

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