"Fossies" - the Fresh Open Source Software Archive

Member "armadillo-9.800.3/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp" (16 Jun 2016, 4191 Bytes) of package /linux/misc/armadillo-9.800.3.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "newarp_UpperHessenbergEigen_meat.hpp" see the Fossies "Dox" file reference documentation.

    1 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
    2 // Copyright 2008-2016 National ICT Australia (NICTA)
    3 // 
    4 // Licensed under the Apache License, Version 2.0 (the "License");
    5 // you may not use this file except in compliance with the License.
    6 // You may obtain a copy of the License at
    7 // http://www.apache.org/licenses/LICENSE-2.0
    8 // 
    9 // Unless required by applicable law or agreed to in writing, software
   10 // distributed under the License is distributed on an "AS IS" BASIS,
   11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   12 // See the License for the specific language governing permissions and
   13 // limitations under the License.
   14 // ------------------------------------------------------------------------
   15 
   16 
   17 namespace newarp
   18 {
   19 
   20 
   21 template<typename eT>
   22 inline
   23 UpperHessenbergEigen<eT>::UpperHessenbergEigen()
   24   : n(0)
   25   , computed(false)
   26   {
   27   arma_extra_debug_sigprint();
   28   }
   29 
   30 
   31 
   32 template<typename eT>
   33 inline
   34 UpperHessenbergEigen<eT>::UpperHessenbergEigen(const Mat<eT>& mat_obj)
   35   : n(mat_obj.n_rows)
   36   , computed(false)
   37   {
   38   arma_extra_debug_sigprint();
   39   
   40   compute(mat_obj);
   41   }
   42 
   43 
   44 
   45 template<typename eT>
   46 inline
   47 void
   48 UpperHessenbergEigen<eT>::compute(const Mat<eT>& mat_obj)
   49   {
   50   arma_extra_debug_sigprint();
   51   
   52   arma_debug_check( (mat_obj.is_square() == false), "newarp::UpperHessenbergEigen::compute(): matrix must be square" );
   53   
   54   n = blas_int(mat_obj.n_rows);
   55   
   56   mat_Z.set_size(n, n);
   57   mat_T.set_size(n, n);
   58   evals.set_size(n);
   59   
   60   mat_Z.eye();
   61   mat_T = mat_obj;
   62   
   63   blas_int want_T = blas_int(1);
   64   blas_int want_Z = blas_int(1);
   65   
   66   blas_int ilo  = blas_int(1);
   67   blas_int ihi  = blas_int(n);
   68   blas_int iloz = blas_int(1);
   69   blas_int ihiz = blas_int(n);
   70   
   71   blas_int info = blas_int(0);
   72   
   73   podarray<eT> wr(static_cast<uword>(n));
   74   podarray<eT> wi(static_cast<uword>(n));
   75   
   76   
   77   lapack::lahqr(&want_T, &want_Z, &n, &ilo, &ihi, mat_T.memptr(), &n, wr.memptr(), wi.memptr(), &iloz, &ihiz, mat_Z.memptr(), &n, &info);
   78   
   79   for(blas_int i = 0; i < n; i++)
   80     {
   81     evals(i) = std::complex<eT>(wr[i], wi[i]);
   82     }
   83   
   84   if(info > 0)  { arma_stop_runtime_error("lapack::lahqr(): failed to compute all eigenvalues"); return; }
   85   
   86   char     side   = 'R';
   87   char     howmny = 'B';
   88   blas_int m      = blas_int(0);
   89   
   90   podarray<eT> work(static_cast<uword>(3 * n));
   91   
   92   lapack::trevc(&side, &howmny, (blas_int*) NULL, &n, mat_T.memptr(), &n, (eT*) NULL, &n, mat_Z.memptr(), &n, &n, &m, work.memptr(), &info);
   93   
   94   if(info < 0)  { arma_stop_logic_error("lapack::trevc(): illegal value"); return; }
   95   
   96   computed = true;
   97   }
   98 
   99 
  100 
  101 template<typename eT>
  102 inline
  103 Col< std::complex<eT> >
  104 UpperHessenbergEigen<eT>::eigenvalues()
  105   {
  106   arma_extra_debug_sigprint();
  107   
  108   arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvalues(): need to call compute() first" );
  109 
  110   return evals;
  111   }
  112 
  113 
  114 
  115 template<typename eT>
  116 inline
  117 Mat< std::complex<eT> >
  118 UpperHessenbergEigen<eT>::eigenvectors()
  119   {
  120   arma_extra_debug_sigprint();
  121   
  122   arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvectors(): need to call compute() first" );
  123 
  124   // Lapack will set the imaginary parts of real eigenvalues to be exact zero
  125   Mat< std::complex<eT> > evecs(n, n);
  126   
  127   std::complex<eT>* col_ptr = evecs.memptr();
  128   
  129   for(blas_int i = 0; i < n; i++)
  130     {
  131     if(cx_attrib::is_real(evals(i), eT(0)))
  132       {
  133       // for real eigenvector, normalise and copy
  134       eT z_norm = norm(mat_Z.col(i));
  135       
  136       for(blas_int j = 0; j < n; j++)
  137         {
  138         col_ptr[j] = std::complex<eT>(mat_Z(j, i) / z_norm, eT(0));
  139         }
  140 
  141       col_ptr += n;
  142       }
  143     else
  144       {
  145       // complex eigenvectors are stored in consecutive columns
  146       eT r2 = dot(mat_Z.col(i), mat_Z.col(i));
  147       eT i2 = dot(mat_Z.col(i + 1), mat_Z.col(i + 1));
  148       
  149       eT  z_norm = std::sqrt(r2 + i2);
  150       eT* z_ptr  = mat_Z.colptr(i);
  151       
  152       for(blas_int j = 0; j < n; j++)
  153         {
  154         col_ptr[j    ] = std::complex<eT>(z_ptr[j] / z_norm, z_ptr[j + n] / z_norm);
  155         col_ptr[j + n] = std::conj(col_ptr[j]);
  156         }
  157 
  158       i++;
  159       col_ptr += 2 * n;
  160       }
  161     }
  162 
  163   return evecs;
  164   }
  165 
  166 
  167 }  // namespace newarp