max@0: // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
max@0: // Copyright (C) 2008-2011 Conrad Sanderson
max@0: // Copyright (C) 2009 Edmund Highcock
max@0: // Copyright (C) 2011 Stanislav Funiak
max@0: // 
max@0: // This file is part of the Armadillo C++ library.
max@0: // It is provided without any warranty of fitness
max@0: // for any purpose. You can redistribute this file
max@0: // and/or modify it under the terms of the GNU
max@0: // Lesser General Public License (LGPL) as published
max@0: // by the Free Software Foundation, either version 3
max@0: // of the License or (at your option) any later version.
max@0: // (see http://www.opensource.org/licenses for more info)
max@0: 
max@0: 
max@0: //! \addtogroup fn_eig
max@0: //! @{
max@0: 
max@0: 
max@0: //
max@0: // symmetric/hermitian matrices
max@0: //
max@0: 
max@0: 
max@0: //! Eigenvalues of real/complex symmetric/hermitian matrix X
max@0: template<typename T1>
max@0: inline
max@0: bool
max@0: eig_sym
max@0:   (
max@0:          Col<typename T1::pod_type>&     eigval,
max@0:   const Base<typename T1::elem_type,T1>& X,
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   // unwrap_check not used as T1::elem_type and T1::pod_type may not be the same.
max@0:   // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::eig_sym() makes a copy of X
max@0:   
max@0:   const bool status = auxlib::eig_sym(eigval, X);
max@0:   
max@0:   if(status == false)
max@0:     {
max@0:     eigval.reset();
max@0:     arma_bad("eig_sym(): failed to converge", false);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues of real/complex symmetric/hermitian matrix X
max@0: template<typename T1>
max@0: inline
max@0: Col<typename T1::pod_type>
max@0: eig_sym
max@0:   (
max@0:   const Base<typename T1::elem_type,T1>& X,
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   Col<typename T1::pod_type> out;
max@0:   const bool status = auxlib::eig_sym(out, X);
max@0: 
max@0:   if(status == false)
max@0:     {
max@0:     out.reset();
max@0:     arma_bad("eig_sym(): failed to converge");
max@0:     }
max@0:   
max@0:   return out;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors of real/complex symmetric/hermitian matrix X
max@0: template<typename T1> 
max@0: inline
max@0: bool
max@0: eig_sym
max@0:   (
max@0:          Col<typename T1::pod_type>&     eigval,
max@0:          Mat<typename T1::elem_type>&    eigvec,
max@0:   const Base<typename T1::elem_type,T1>& X,
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_sym(): eigval is an alias of eigvec" );
max@0:   
max@0:   const bool status = auxlib::eig_sym(eigval, eigvec, X);
max@0:   
max@0:   if(status == false)
max@0:     {
max@0:     eigval.reset();
max@0:     eigvec.reset();
max@0:     arma_bad("eig_sym(): failed to converge", false);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // general matrices
max@0: //
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors (both left and right) of general real/complex square matrix X
max@0: template<typename T1>
max@0: inline
max@0: bool
max@0: eig_gen
max@0:   (
max@0:          Col< std::complex<typename T1::pod_type> >& eigval, 
max@0:          Mat<typename T1::elem_type>&                l_eigvec,
max@0:          Mat<typename T1::elem_type>&                r_eigvec,
max@0:   const Base<typename T1::elem_type,T1>&             X,
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   arma_debug_check
max@0:     (
max@0:     ((&l_eigvec) == (&r_eigvec)),
max@0:     "eig_gen(): l_eigvec is an alias of r_eigvec"
max@0:     );
max@0:   
max@0:   arma_debug_check
max@0:     (
max@0:       (
max@0:       (((void*)(&eigval)) == ((void*)(&l_eigvec)))
max@0:       ||
max@0:       (((void*)(&eigval)) == ((void*)(&r_eigvec)))
max@0:       ),
max@0:     "eig_gen(): eigval is an alias of l_eigvec or r_eigvec"
max@0:     );
max@0:   
max@0:   const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b');
max@0:   
max@0:   if(status == false)
max@0:     {
max@0:     eigval.reset();
max@0:     l_eigvec.reset();
max@0:     r_eigvec.reset();
max@0:     arma_bad("eig_gen(): failed to converge", false);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors of general real square matrix X.
max@0: //! Optional argument 'side' specifies which eigenvectors should be computed:
max@0: //! 'r' for right (default) and 'l' for left.
max@0: template<typename eT, typename T1>
max@0: inline
max@0: bool
max@0: eig_gen
max@0:   (
max@0:         Col< std::complex<eT> >& eigval, 
max@0:         Mat< std::complex<eT> >& eigvec,
max@0:   const Base<eT, T1>&            X, 
max@0:   const char                     side = 'r',
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   //std::cout << "real" << std::endl;
max@0:   
max@0:   arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
max@0:   
max@0:   Mat<eT> dummy_eigvec;
max@0:   Mat<eT> tmp_eigvec;
max@0:   
max@0:   bool status;
max@0:   
max@0:   switch(side)
max@0:     {
max@0:     case 'r':
max@0:       status = auxlib::eig_gen(eigval, dummy_eigvec, tmp_eigvec, X, side);
max@0:       break;
max@0:     
max@0:     case 'l':
max@0:       status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side);
max@0:       break;
max@0:       
max@0:     default:
max@0:       arma_stop("eig_gen(): parameter 'side' is invalid");
max@0:       status = false;
max@0:     }
max@0:   
max@0:   if(status == false)
max@0:     {
max@0:     eigval.reset();
max@0:     eigvec.reset();
max@0:     arma_bad("eig_gen(): failed to converge", false);
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword n = eigval.n_elem;
max@0:     
max@0:     if(n > 0)
max@0:       {
max@0:       eigvec.set_size(n,n);
max@0:       
max@0:       for(uword j=0; j<n; ++j)
max@0:         {
max@0:         if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) )
max@0:           {
max@0:           // eigvec.col(j)   = Mat< std::complex<eT> >( tmp_eigvec.col(j),  tmp_eigvec.col(j+1) );
max@0:           // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -tmp_eigvec.col(j+1) );
max@0:           
max@0:           for(uword i=0; i<n; ++i)
max@0:             {
max@0:             eigvec.at(i,j)   = std::complex<eT>( tmp_eigvec.at(i,j),  tmp_eigvec.at(i,j+1) );
max@0:             eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eigvec.at(i,j+1) );
max@0:             }
max@0:           
max@0:           ++j;
max@0:           }
max@0:         else
max@0:           {
max@0:           // eigvec.col(i) = tmp_eigvec.col(i);
max@0:           
max@0:           for(uword i=0; i<n; ++i)
max@0:             {
max@0:             eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0));
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       }
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! Eigenvalues and eigenvectors of general complex square matrix X
max@0: //! Optional argument 'side' specifies which eigenvectors should be computed:
max@0: //! 'r' for right (default) and 'l' for left.
max@0: template<typename T, typename T1>
max@0: inline
max@0: bool
max@0: eig_gen
max@0:   (
max@0:          Col<std::complex<T> >&    eigval, 
max@0:          Mat<std::complex<T> >&    eigvec,
max@0:   const Base<std::complex<T>, T1>& X, 
max@0:   const char                       side = 'r',
max@0:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   //std::cout << "complex" << std::endl;
max@0:   
max@0:   arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" );
max@0:   
max@0:   Mat< std::complex<T> > dummy_eigvec;
max@0:   
max@0:   bool status;
max@0:   
max@0:   switch(side)
max@0:     {
max@0:     case 'r':
max@0:       status = auxlib::eig_gen(eigval, dummy_eigvec, eigvec, X, side);
max@0:       break;
max@0:     
max@0:     case 'l':
max@0:       status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side);
max@0:       break;
max@0:       
max@0:     default:
max@0:       arma_stop("eig_gen(): parameter 'side' is invalid");
max@0:       status = false;
max@0:     }
max@0:   
max@0:   if(status == false)
max@0:     {
max@0:     eigval.reset();
max@0:     eigvec.reset();
max@0:     arma_bad("eig_gen(): failed to converge", false);
max@0:     }
max@0:   
max@0:   return status;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! @}
max@0: