From eda5bc26f44ee9a6f83dcf8c91f17296d7fc509d Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Mon, 12 Feb 2024 14:52:43 +0100 Subject: Move into version control --- .../newarp_UpperHessenbergEigen_meat.hpp | 168 +++++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 src/armadillo/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp (limited to 'src/armadillo/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp') diff --git a/src/armadillo/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp b/src/armadillo/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp new file mode 100644 index 0000000..a758205 --- /dev/null +++ b/src/armadillo/include/armadillo_bits/newarp_UpperHessenbergEigen_meat.hpp @@ -0,0 +1,168 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + +namespace newarp +{ + + +template +inline +UpperHessenbergEigen::UpperHessenbergEigen() + : n_rows(0) + , computed(false) + { + arma_extra_debug_sigprint(); + } + + + +template +inline +UpperHessenbergEigen::UpperHessenbergEigen(const Mat& mat_obj) + : n_rows(mat_obj.n_rows) + , computed(false) + { + arma_extra_debug_sigprint(); + + compute(mat_obj); + } + + + +template +inline +void +UpperHessenbergEigen::compute(const Mat& mat_obj) + { + arma_extra_debug_sigprint(); + + arma_debug_check( (mat_obj.is_square() == false), "newarp::UpperHessenbergEigen::compute(): matrix must be square" ); + + n_rows = mat_obj.n_rows; + + mat_Z.set_size(n_rows, n_rows); + mat_T.set_size(n_rows, n_rows); + evals.set_size(n_rows); + + mat_Z.eye(); + mat_T = mat_obj; + + blas_int want_T = 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 ihi = blas_int(n_rows); + blas_int iloz = blas_int(1); + blas_int ihiz = blas_int(n_rows); + + blas_int info = blas_int(0); + + podarray wr(n_rows); + podarray 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); + + if(info != 0) { arma_stop_runtime_error("lapack::lahqr(): failed to compute all eigenvalues"); return; } + + for(uword i=0; i < n_rows; i++) { evals(i) = std::complex(wr[i], wi[i]); } + + char side = 'R'; + char howmny = 'B'; + blas_int m = blas_int(0); + + podarray 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); + + if(info != 0) { arma_stop_runtime_error("lapack::trevc(): illegal value"); return; } + + computed = true; + } + + + +template +inline +Col< std::complex > +UpperHessenbergEigen::eigenvalues() + { + arma_extra_debug_sigprint(); + + arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvalues(): need to call compute() first" ); + + return evals; + } + + + +template +inline +Mat< std::complex > +UpperHessenbergEigen::eigenvectors() + { + arma_extra_debug_sigprint(); + + arma_debug_check( (computed == false), "newarp::UpperHessenbergEigen::eigenvectors(): need to call compute() first" ); + + // Lapack will set the imaginary parts of real eigenvalues to be exact zero + Mat< std::complex > evecs(n_rows, n_rows, arma_zeros_indicator()); + + std::complex* col_ptr = evecs.memptr(); + + for(uword i=0; i < n_rows; i++) + { + if(cx_attrib::is_real(evals(i), eT(0))) + { + // for real eigenvector, normalise and copy + const eT z_norm = norm(mat_Z.col(i)); + + for(uword j=0; j < n_rows; j++) + { + col_ptr[j] = std::complex(mat_Z(j, i) / z_norm, eT(0)); + } + + col_ptr += n_rows; + } + else + { + // complex eigenvectors are stored in consecutive columns + const eT r2 = dot(mat_Z.col(i ), mat_Z.col(i )); + const eT i2 = dot(mat_Z.col(i+1), mat_Z.col(i+1)); + + const eT z_norm = std::sqrt(r2 + i2); + const eT* z_ptr = mat_Z.colptr(i); + + for(uword j=0; j < n_rows; j++) + { + col_ptr[j ] = std::complex(z_ptr[j] / z_norm, z_ptr[j + n_rows] / z_norm); + col_ptr[j + n_rows] = std::conj(col_ptr[j]); + } + + i++; + col_ptr += 2 * n_rows; + } + } + + return evecs; + } + + +} // namespace newarp -- cgit v1.2.1