"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp" between
armadillo-10.8.2.tar.xz and armadillo-11.0.0.tar.xz

About: Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use.

newarp_UpperHessenbergEigen_meat.hpp  (armadillo-10.8.2.tar.xz):newarp_UpperHessenbergEigen_meat.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 24 skipping to change at line 24
// See the License for the specific language governing permissions and // See the License for the specific language governing permissions and
// limitations under the License. // limitations under the License.
// ------------------------------------------------------------------------ // ------------------------------------------------------------------------
namespace newarp namespace newarp
{ {
template<typename eT> template<typename eT>
inline inline
UpperHessenbergEigen<eT>::UpperHessenbergEigen() UpperHessenbergEigen<eT>::UpperHessenbergEigen()
: n(0) : n_rows(0)
, computed(false) , computed(false)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
} }
template<typename eT> template<typename eT>
inline inline
UpperHessenbergEigen<eT>::UpperHessenbergEigen(const Mat<eT>& mat_obj) UpperHessenbergEigen<eT>::UpperHessenbergEigen(const Mat<eT>& mat_obj)
: n(mat_obj.n_rows) : n_rows(mat_obj.n_rows)
, computed(false) , computed(false)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
compute(mat_obj); compute(mat_obj);
} }
template<typename eT> template<typename eT>
inline inline
void void
UpperHessenbergEigen<eT>::compute(const Mat<eT>& mat_obj) UpperHessenbergEigen<eT>::compute(const Mat<eT>& mat_obj)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
arma_debug_check( (mat_obj.is_square() == false), "newarp::UpperHessenbergEige n::compute(): matrix must be square" ); arma_debug_check( (mat_obj.is_square() == false), "newarp::UpperHessenbergEige n::compute(): matrix must be square" );
n = blas_int(mat_obj.n_rows); n_rows = mat_obj.n_rows;
mat_Z.set_size(n, n); mat_Z.set_size(n_rows, n_rows);
mat_T.set_size(n, n); mat_T.set_size(n_rows, n_rows);
evals.set_size(n); evals.set_size(n_rows);
mat_Z.eye(); mat_Z.eye();
mat_T = mat_obj; mat_T = mat_obj;
blas_int want_T = blas_int(1); blas_int want_T = blas_int(1);
blas_int want_Z = blas_int(1); blas_int want_Z = blas_int(1);
blas_int n = blas_int(n_rows);
blas_int ilo = blas_int(1); blas_int ilo = blas_int(1);
blas_int ihi = blas_int(n); blas_int ihi = blas_int(n_rows);
blas_int iloz = blas_int(1); blas_int iloz = blas_int(1);
blas_int ihiz = blas_int(n); blas_int ihiz = blas_int(n_rows);
blas_int info = blas_int(0); blas_int info = blas_int(0);
podarray<eT> wr(static_cast<uword>(n)); podarray<eT> wr(n_rows);
podarray<eT> wi(static_cast<uword>(n)); podarray<eT> wi(n_rows);
arma_extra_debug_print("lapack::lahqr()");
lapack::lahqr(&want_T, &want_Z, &n, &ilo, &ihi, mat_T.memptr(), &n, wr.memptr( ), wi.memptr(), &iloz, &ihiz, mat_Z.memptr(), &n, &info); lapack::lahqr(&want_T, &want_Z, &n, &ilo, &ihi, mat_T.memptr(), &n, wr.memptr( ), wi.memptr(), &iloz, &ihiz, mat_Z.memptr(), &n, &info);
for(blas_int i = 0; i < n; i++) if(info != 0) { arma_stop_runtime_error("lapack::lahqr(): failed to compute a
{ ll eigenvalues"); return; }
evals(i) = std::complex<eT>(wr[i], wi[i]);
}
if(info > 0) { arma_stop_runtime_error("lapack::lahqr(): failed to compute al l eigenvalues"); return; } for(uword i=0; i < n_rows; i++) { evals(i) = std::complex<eT>(wr[i], wi[i]); }
char side = 'R'; char side = 'R';
char howmny = 'B'; char howmny = 'B';
blas_int m = blas_int(0); blas_int m = blas_int(0);
podarray<eT> work(static_cast<uword>(3 * n)); podarray<eT> work(3*n);
arma_extra_debug_print("lapack::trevc()");
lapack::trevc(&side, &howmny, (blas_int*) NULL, &n, mat_T.memptr(), &n, (eT*) NULL, &n, mat_Z.memptr(), &n, &n, &m, work.memptr(), &info); lapack::trevc(&side, &howmny, (blas_int*) NULL, &n, mat_T.memptr(), &n, (eT*) NULL, &n, mat_Z.memptr(), &n, &n, &m, work.memptr(), &info);
if(info < 0) { arma_stop_logic_error("lapack::trevc(): illegal value"); retur n; } if(info != 0) { arma_stop_runtime_error("lapack::trevc(): illegal value"); re turn; }
computed = true; computed = true;
} }
template<typename eT> template<typename eT>
inline inline
Col< std::complex<eT> > Col< std::complex<eT> >
UpperHessenbergEigen<eT>::eigenvalues() UpperHessenbergEigen<eT>::eigenvalues()
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 116 skipping to change at line 116
template<typename eT> template<typename eT>
inline inline
Mat< std::complex<eT> > Mat< std::complex<eT> >
UpperHessenbergEigen<eT>::eigenvectors() UpperHessenbergEigen<eT>::eigenvectors()
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvec tors(): need to call compute() first" ); arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvec tors(): need to call compute() first" );
// Lapack will set the imaginary parts of real eigenvalues to be exact zero // Lapack will set the imaginary parts of real eigenvalues to be exact zero
Mat< std::complex<eT> > evecs(n, n, arma_zeros_indicator()); Mat< std::complex<eT> > evecs(n_rows, n_rows, arma_zeros_indicator());
std::complex<eT>* col_ptr = evecs.memptr(); std::complex<eT>* col_ptr = evecs.memptr();
for(blas_int i = 0; i < n; i++) for(uword i=0; i < n_rows; i++)
{ {
if(cx_attrib::is_real(evals(i), eT(0))) if(cx_attrib::is_real(evals(i), eT(0)))
{ {
// for real eigenvector, normalise and copy // for real eigenvector, normalise and copy
eT z_norm = norm(mat_Z.col(i)); const eT z_norm = norm(mat_Z.col(i));
for(blas_int j = 0; j < n; j++) for(uword j=0; j < n_rows; j++)
{ {
col_ptr[j] = std::complex<eT>(mat_Z(j, i) / z_norm, eT(0)); col_ptr[j] = std::complex<eT>(mat_Z(j, i) / z_norm, eT(0));
} }
col_ptr += n; col_ptr += n_rows;
} }
else else
{ {
// complex eigenvectors are stored in consecutive columns // complex eigenvectors are stored in consecutive columns
eT r2 = dot(mat_Z.col(i), mat_Z.col(i)); const eT r2 = dot(mat_Z.col(i ), mat_Z.col(i ));
eT i2 = dot(mat_Z.col(i + 1), mat_Z.col(i + 1)); const eT i2 = dot(mat_Z.col(i+1), mat_Z.col(i+1));
eT z_norm = std::sqrt(r2 + i2); const eT z_norm = std::sqrt(r2 + i2);
eT* z_ptr = mat_Z.colptr(i); const eT* z_ptr = mat_Z.colptr(i);
for(blas_int j = 0; j < n; j++) for(uword j=0; j < n_rows; j++)
{ {
col_ptr[j ] = std::complex<eT>(z_ptr[j] / z_norm, z_ptr[j + n] / z_no col_ptr[j ] = std::complex<eT>(z_ptr[j] / z_norm, z_ptr[j + n_ro
rm); ws] / z_norm);
col_ptr[j + n] = std::conj(col_ptr[j]); col_ptr[j + n_rows] = std::conj(col_ptr[j]);
} }
i++; i++;
col_ptr += 2 * n; col_ptr += 2 * n_rows;
} }
} }
return evecs; return evecs;
} }
} // namespace newarp } // namespace newarp
 End of changes. 24 change blocks. 
31 lines changed or deleted 32 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)