$search
00001 // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) 00002 // Copyright (C) 2008-2011 Conrad Sanderson 00003 // Copyright (C) 2009 Edmund Highcock 00004 // Copyright (C) 2011 Stanislav Funiak 00005 // 00006 // This file is part of the Armadillo C++ library. 00007 // It is provided without any warranty of fitness 00008 // for any purpose. You can redistribute this file 00009 // and/or modify it under the terms of the GNU 00010 // Lesser General Public License (LGPL) as published 00011 // by the Free Software Foundation, either version 3 00012 // of the License or (at your option) any later version. 00013 // (see http://www.opensource.org/licenses for more info) 00014 00015 00018 00019 00020 // 00021 // symmetric/hermitian matrices 00022 // 00023 00024 00026 template<typename T1> 00027 inline 00028 bool 00029 eig_sym 00030 ( 00031 Col<typename T1::pod_type>& eigval, 00032 const Base<typename T1::elem_type,T1>& X, 00033 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00034 ) 00035 { 00036 arma_extra_debug_sigprint(); 00037 arma_ignore(junk); 00038 00039 // unwrap_check not used as T1::elem_type and T1::pod_type may not be the same. 00040 // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::eig_sym() makes a copy of X 00041 00042 const bool status = auxlib::eig_sym(eigval, X); 00043 00044 if(status == false) 00045 { 00046 eigval.reset(); 00047 arma_bad("eig_sym(): failed to converge", false); 00048 } 00049 00050 return status; 00051 } 00052 00053 00054 00056 template<typename T1> 00057 inline 00058 Col<typename T1::pod_type> 00059 eig_sym 00060 ( 00061 const Base<typename T1::elem_type,T1>& X, 00062 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00063 ) 00064 { 00065 arma_extra_debug_sigprint(); 00066 arma_ignore(junk); 00067 00068 Col<typename T1::pod_type> out; 00069 const bool status = auxlib::eig_sym(out, X); 00070 00071 if(status == false) 00072 { 00073 out.reset(); 00074 arma_bad("eig_sym(): failed to converge"); 00075 } 00076 00077 return out; 00078 } 00079 00080 00081 00083 template<typename T1> 00084 inline 00085 bool 00086 eig_sym 00087 ( 00088 Col<typename T1::pod_type>& eigval, 00089 Mat<typename T1::elem_type>& eigvec, 00090 const Base<typename T1::elem_type,T1>& X, 00091 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00092 ) 00093 { 00094 arma_extra_debug_sigprint(); 00095 arma_ignore(junk); 00096 00097 arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_sym(): eigval is an alias of eigvec" ); 00098 00099 const bool status = auxlib::eig_sym(eigval, eigvec, X); 00100 00101 if(status == false) 00102 { 00103 eigval.reset(); 00104 eigvec.reset(); 00105 arma_bad("eig_sym(): failed to converge", false); 00106 } 00107 00108 return status; 00109 } 00110 00111 00112 00113 // 00114 // general matrices 00115 // 00116 00117 00118 00120 template<typename T1> 00121 inline 00122 bool 00123 eig_gen 00124 ( 00125 Col< std::complex<typename T1::pod_type> >& eigval, 00126 Mat<typename T1::elem_type>& l_eigvec, 00127 Mat<typename T1::elem_type>& r_eigvec, 00128 const Base<typename T1::elem_type,T1>& X, 00129 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00130 ) 00131 { 00132 arma_extra_debug_sigprint(); 00133 arma_ignore(junk); 00134 00135 arma_debug_check 00136 ( 00137 ((&l_eigvec) == (&r_eigvec)), 00138 "eig_gen(): l_eigvec is an alias of r_eigvec" 00139 ); 00140 00141 arma_debug_check 00142 ( 00143 ( 00144 (((void*)(&eigval)) == ((void*)(&l_eigvec))) 00145 || 00146 (((void*)(&eigval)) == ((void*)(&r_eigvec))) 00147 ), 00148 "eig_gen(): eigval is an alias of l_eigvec or r_eigvec" 00149 ); 00150 00151 const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b'); 00152 00153 if(status == false) 00154 { 00155 eigval.reset(); 00156 l_eigvec.reset(); 00157 r_eigvec.reset(); 00158 arma_bad("eig_gen(): failed to converge", false); 00159 } 00160 00161 return status; 00162 } 00163 00164 00165 00169 template<typename eT, typename T1> 00170 inline 00171 bool 00172 eig_gen 00173 ( 00174 Col< std::complex<eT> >& eigval, 00175 Mat< std::complex<eT> >& eigvec, 00176 const Base<eT, T1>& X, 00177 const char side = 'r', 00178 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00179 ) 00180 { 00181 arma_extra_debug_sigprint(); 00182 arma_ignore(junk); 00183 00184 //std::cout << "real" << std::endl; 00185 00186 arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" ); 00187 00188 Mat<eT> dummy_eigvec; 00189 Mat<eT> tmp_eigvec; 00190 00191 bool status; 00192 00193 switch(side) 00194 { 00195 case 'r': 00196 status = auxlib::eig_gen(eigval, dummy_eigvec, tmp_eigvec, X, side); 00197 break; 00198 00199 case 'l': 00200 status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side); 00201 break; 00202 00203 default: 00204 arma_stop("eig_gen(): parameter 'side' is invalid"); 00205 status = false; 00206 } 00207 00208 if(status == false) 00209 { 00210 eigval.reset(); 00211 eigvec.reset(); 00212 arma_bad("eig_gen(): failed to converge", false); 00213 } 00214 else 00215 { 00216 const uword n = eigval.n_elem; 00217 00218 if(n > 0) 00219 { 00220 eigvec.set_size(n,n); 00221 00222 for(uword j=0; j<n; ++j) 00223 { 00224 if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) ) 00225 { 00226 // eigvec.col(j) = Mat< std::complex<eT> >( tmp_eigvec.col(j), tmp_eigvec.col(j+1) ); 00227 // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -tmp_eigvec.col(j+1) ); 00228 00229 for(uword i=0; i<n; ++i) 00230 { 00231 eigvec.at(i,j) = std::complex<eT>( tmp_eigvec.at(i,j), tmp_eigvec.at(i,j+1) ); 00232 eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eigvec.at(i,j+1) ); 00233 } 00234 00235 ++j; 00236 } 00237 else 00238 { 00239 // eigvec.col(i) = tmp_eigvec.col(i); 00240 00241 for(uword i=0; i<n; ++i) 00242 { 00243 eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0)); 00244 } 00245 00246 } 00247 } 00248 } 00249 } 00250 00251 return status; 00252 } 00253 00254 00255 00259 template<typename T, typename T1> 00260 inline 00261 bool 00262 eig_gen 00263 ( 00264 Col<std::complex<T> >& eigval, 00265 Mat<std::complex<T> >& eigvec, 00266 const Base<std::complex<T>, T1>& X, 00267 const char side = 'r', 00268 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0 00269 ) 00270 { 00271 arma_extra_debug_sigprint(); 00272 arma_ignore(junk); 00273 00274 //std::cout << "complex" << std::endl; 00275 00276 arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(): eigval is an alias of eigvec" ); 00277 00278 Mat< std::complex<T> > dummy_eigvec; 00279 00280 bool status; 00281 00282 switch(side) 00283 { 00284 case 'r': 00285 status = auxlib::eig_gen(eigval, dummy_eigvec, eigvec, X, side); 00286 break; 00287 00288 case 'l': 00289 status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side); 00290 break; 00291 00292 default: 00293 arma_stop("eig_gen(): parameter 'side' is invalid"); 00294 status = false; 00295 } 00296 00297 if(status == false) 00298 { 00299 eigval.reset(); 00300 eigvec.reset(); 00301 arma_bad("eig_gen(): failed to converge", false); 00302 } 00303 00304 return status; 00305 } 00306 00307 00308 00310