Chris@49: // Copyright (C) 2008-2012 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2008-2012 Conrad Sanderson
Chris@49: // Copyright (C) 2009 Edmund Highcock
Chris@49: // Copyright (C) 2011 Stanislav Funiak
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup fn_eig
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: //
Chris@49: // symmetric/hermitian matrices
Chris@49: //
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues of real/complex symmetric/hermitian matrix X
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: bool
Chris@49: eig_sym
Chris@49:   (
Chris@49:          Col<typename T1::pod_type>&     eigval,
Chris@49:   const Base<typename T1::elem_type,T1>& X,
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   // unwrap_check not used as T1::elem_type and T1::pod_type may not be the same.
Chris@49:   // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::eig_sym() makes a copy of X
Chris@49:   
Chris@49:   const bool status = auxlib::eig_sym(eigval, X);
Chris@49:   
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     eigval.reset();
Chris@49:     arma_bad("eig_sym(): failed to converge", false);
Chris@49:     }
Chris@49:   
Chris@49:   return status;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues of real/complex symmetric/hermitian matrix X
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: Col<typename T1::pod_type>
Chris@49: eig_sym
Chris@49:   (
Chris@49:   const Base<typename T1::elem_type,T1>& X,
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   Col<typename T1::pod_type> out;
Chris@49:   const bool status = auxlib::eig_sym(out, X);
Chris@49: 
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     out.reset();
Chris@49:     arma_bad("eig_sym(): failed to converge");
Chris@49:     }
Chris@49:   
Chris@49:   return out;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues and eigenvectors of real/complex symmetric/hermitian matrix X
Chris@49: template<typename T1> 
Chris@49: inline
Chris@49: bool
Chris@49: eig_sym
Chris@49:   (
Chris@49:          Col<typename T1::pod_type>&     eigval,
Chris@49:          Mat<typename T1::elem_type>&    eigvec,
Chris@49:   const Base<typename T1::elem_type,T1>& X,
Chris@49:   const char* method =                   "",
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   arma_debug_check( void_ptr(&eigval) == void_ptr(&eigvec), "eig_sym(): eigval is an alias of eigvec" );
Chris@49:   
Chris@49:   bool use_divide_and_conquer = false;
Chris@49:   
Chris@49:   const char sig = method[0];
Chris@49:   
Chris@49:   switch(sig)
Chris@49:     {
Chris@49:     case '\0':
Chris@49:     case 's':
Chris@49:       break;
Chris@49:       
Chris@49:     case 'd':
Chris@49:       use_divide_and_conquer = true;
Chris@49:       break;
Chris@49:     
Chris@49:     default:
Chris@49:       {
Chris@49:       arma_stop("eig_sym(): unknown method specified");
Chris@49:       return false;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   const bool status = (use_divide_and_conquer == false) ? auxlib::eig_sym(eigval, eigvec, X) : auxlib::eig_sym_dc(eigval, eigvec, X);
Chris@49:   
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     eigval.reset();
Chris@49:     eigvec.reset();
Chris@49:     arma_bad("eig_sym(): failed to converge", false);
Chris@49:     }
Chris@49:   
Chris@49:   return status;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //
Chris@49: // general matrices
Chris@49: //
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues and eigenvectors (both left and right) of general real/complex square matrix X
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: bool
Chris@49: eig_gen
Chris@49:   (
Chris@49:          Col< std::complex<typename T1::pod_type> >& eigval, 
Chris@49:          Mat<typename T1::elem_type>&                l_eigvec,
Chris@49:          Mat<typename T1::elem_type>&                r_eigvec,
Chris@49:   const Base<typename T1::elem_type,T1>&             X,
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     ((&l_eigvec) == (&r_eigvec)),
Chris@49:     "eig_gen(): l_eigvec is an alias of r_eigvec"
Chris@49:     );
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:       (
Chris@49:       (((void*)(&eigval)) == ((void*)(&l_eigvec)))
Chris@49:       ||
Chris@49:       (((void*)(&eigval)) == ((void*)(&r_eigvec)))
Chris@49:       ),
Chris@49:     "eig_gen(): eigval is an alias of l_eigvec or r_eigvec"
Chris@49:     );
Chris@49:   
Chris@49:   const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b');
Chris@49:   
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     eigval.reset();
Chris@49:     l_eigvec.reset();
Chris@49:     r_eigvec.reset();
Chris@49:     arma_bad("eig_gen(): failed to converge", false);
Chris@49:     }
Chris@49:   
Chris@49:   return status;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues and eigenvectors of general real square matrix X.
Chris@49: //! Optional argument 'side' specifies which eigenvectors should be computed:
Chris@49: //! 'r' for right (default) and 'l' for left.
Chris@49: template<typename eT, typename T1>
Chris@49: inline
Chris@49: bool
Chris@49: eig_gen
Chris@49:   (
Chris@49:         Col< std::complex<eT> >& eigval, 
Chris@49:         Mat< std::complex<eT> >& eigvec,
Chris@49:   const Base<eT, T1>&            X, 
Chris@49:   const char                     side = 'r',
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   //std::cout << "real" << std::endl;
Chris@49:   
Chris@49:   arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
Chris@49:   
Chris@49:   Mat<eT> dummy_eigvec;
Chris@49:   Mat<eT> tmp_eigvec;
Chris@49:   
Chris@49:   bool status;
Chris@49:   
Chris@49:   switch(side)
Chris@49:     {
Chris@49:     case 'r':
Chris@49:       status = auxlib::eig_gen(eigval, dummy_eigvec, tmp_eigvec, X, side);
Chris@49:       break;
Chris@49:     
Chris@49:     case 'l':
Chris@49:       status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side);
Chris@49:       break;
Chris@49:       
Chris@49:     default:
Chris@49:       arma_stop("eig_gen(): parameter 'side' is invalid");
Chris@49:       status = false;
Chris@49:     }
Chris@49:   
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     eigval.reset();
Chris@49:     eigvec.reset();
Chris@49:     arma_bad("eig_gen(): failed to converge", false);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword n = eigval.n_elem;
Chris@49:     
Chris@49:     if(n > 0)
Chris@49:       {
Chris@49:       eigvec.set_size(n,n);
Chris@49:       
Chris@49:       for(uword j=0; j<n; ++j)
Chris@49:         {
Chris@49:         if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) )
Chris@49:           {
Chris@49:           // eigvec.col(j)   = Mat< std::complex<eT> >( tmp_eigvec.col(j),  tmp_eigvec.col(j+1) );
Chris@49:           // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -tmp_eigvec.col(j+1) );
Chris@49:           
Chris@49:           for(uword i=0; i<n; ++i)
Chris@49:             {
Chris@49:             eigvec.at(i,j)   = std::complex<eT>( tmp_eigvec.at(i,j),  tmp_eigvec.at(i,j+1) );
Chris@49:             eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eigvec.at(i,j+1) );
Chris@49:             }
Chris@49:           
Chris@49:           ++j;
Chris@49:           }
Chris@49:         else
Chris@49:           {
Chris@49:           // eigvec.col(i) = tmp_eigvec.col(i);
Chris@49:           
Chris@49:           for(uword i=0; i<n; ++i)
Chris@49:             {
Chris@49:             eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0));
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return status;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Eigenvalues and eigenvectors of general complex square matrix X
Chris@49: //! Optional argument 'side' specifies which eigenvectors should be computed:
Chris@49: //! 'r' for right (default) and 'l' for left.
Chris@49: template<typename T, typename T1>
Chris@49: inline
Chris@49: bool
Chris@49: eig_gen
Chris@49:   (
Chris@49:          Col<std::complex<T> >&    eigval, 
Chris@49:          Mat<std::complex<T> >&    eigvec,
Chris@49:   const Base<std::complex<T>, T1>& X, 
Chris@49:   const char                       side = 'r',
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   //std::cout << "complex" << std::endl;
Chris@49:   
Chris@49:   arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
Chris@49:   
Chris@49:   Mat< std::complex<T> > dummy_eigvec;
Chris@49:   
Chris@49:   bool status;
Chris@49:   
Chris@49:   switch(side)
Chris@49:     {
Chris@49:     case 'r':
Chris@49:       status = auxlib::eig_gen(eigval, dummy_eigvec, eigvec, X, side);
Chris@49:       break;
Chris@49:     
Chris@49:     case 'l':
Chris@49:       status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side);
Chris@49:       break;
Chris@49:       
Chris@49:     default:
Chris@49:       arma_stop("eig_gen(): parameter 'side' is invalid");
Chris@49:       status = false;
Chris@49:     }
Chris@49:   
Chris@49:   if(status == false)
Chris@49:     {
Chris@49:     eigval.reset();
Chris@49:     eigvec.reset();
Chris@49:     arma_bad("eig_gen(): failed to converge", false);
Chris@49:     }
Chris@49:   
Chris@49:   return status;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}
Chris@49: