diff options
Diffstat (limited to 'src/armadillo/include/armadillo_bits/newarp_DoubleShiftQR_meat.hpp')
| -rw-r--r-- | src/armadillo/include/armadillo_bits/newarp_DoubleShiftQR_meat.hpp | 399 |
1 files changed, 399 insertions, 0 deletions
diff --git a/src/armadillo/include/armadillo_bits/newarp_DoubleShiftQR_meat.hpp b/src/armadillo/include/armadillo_bits/newarp_DoubleShiftQR_meat.hpp new file mode 100644 index 0000000..1c7497d --- /dev/null +++ b/src/armadillo/include/armadillo_bits/newarp_DoubleShiftQR_meat.hpp @@ -0,0 +1,399 @@ +// 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<typename eT> +inline +void +DoubleShiftQR<eT>::compute_reflector(const eT& x1, const eT& x2, const eT& x3, uword ind) + { + arma_extra_debug_sigprint(); + + // In general case the reflector affects 3 rows + ref_nr(ind) = 3; + eT x2x3 = eT(0); + // If x3 is zero, decrease nr by 1 + if(std::abs(x3) < prec) + { + // If x2 is also zero, nr will be 1, and we can exit this function + if(std::abs(x2) < prec) + { + ref_nr(ind) = 1; + return; + } + else + { + ref_nr(ind) = 2; + } + x2x3 = std::abs(x2); + } + else + { + x2x3 = arma_hypot(x2, x3); + } + + // x1' = x1 - rho * ||x|| + // rho = -sign(x1), if x1 == 0, we choose rho = 1 + eT x1_new = x1 - ((x1 <= 0) - (x1 > 0)) * arma_hypot(x1, x2x3); + eT x_norm = arma_hypot(x1_new, x2x3); + // Double check the norm of new x + if(x_norm < prec) + { + ref_nr(ind) = 1; + return; + } + ref_u(0, ind) = x1_new / x_norm; + ref_u(1, ind) = x2 / x_norm; + ref_u(2, ind) = x3 / x_norm; + } + + +template<typename eT> +arma_inline +void +DoubleShiftQR<eT>::compute_reflector(const eT* x, uword ind) + { + arma_extra_debug_sigprint(); + + compute_reflector(x[0], x[1], x[2], ind); + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::update_block(uword il, uword iu) + { + arma_extra_debug_sigprint(); + + // Block size + uword bsize = iu - il + 1; + + // If block size == 1, there is no need to apply reflectors + if(bsize == 1) + { + ref_nr(il) = 1; + return; + } + + // For block size == 2, do a Givens rotation on M = X * X - s * X + t * I + if(bsize == 2) + { + // m00 = x00 * (x00 - s) + x01 * x10 + t + eT m00 = mat_H(il, il) * (mat_H(il, il) - shift_s) + + mat_H(il, il + 1) * mat_H(il + 1, il) + + shift_t; + // m10 = x10 * (x00 + x11 - s) + eT m10 = mat_H(il + 1, il) * (mat_H(il, il) + mat_H(il + 1, il + 1) - shift_s); + // This causes nr=2 + compute_reflector(m00, m10, 0, il); + // Apply the reflector to X + apply_PX(mat_H, il, il, 2, n - il, il); + apply_XP(mat_H, 0, il, il + 2, 2, il); + + ref_nr(il + 1) = 1; + return; + } + + // For block size >=3, use the regular strategy + eT m00 = mat_H(il, il) * (mat_H(il, il) - shift_s) + + mat_H(il, il + 1) * mat_H(il + 1, il) + + shift_t; + eT m10 = mat_H(il + 1, il) * (mat_H(il, il) + mat_H(il + 1, il + 1) - shift_s); + // m20 = x21 * x10 + eT m20 = mat_H(il + 2, il + 1) * mat_H(il + 1, il); + compute_reflector(m00, m10, m20, il); + + // Apply the first reflector + apply_PX(mat_H, il, il, 3, n - il, il); + apply_XP(mat_H, 0, il, il + (std::min)(bsize, uword(4)), 3, il); + + // Calculate the following reflectors + // If entering this loop, block size is at least 4. + for(uword i = 1; i < bsize - 2; i++) + { + compute_reflector(mat_H.colptr(il + i - 1) + il + i, il + i); + // Apply the reflector to X + apply_PX(mat_H, il + i, il + i - 1, 3, n + 1 - il - i, il + i); + apply_XP(mat_H, 0, il + i, il + (std::min)(bsize, uword(i + 4)), 3, il + i); + } + + // The last reflector + // This causes nr=2 + compute_reflector(mat_H(iu - 1, iu - 2), mat_H(iu, iu - 2), 0, iu - 1); + // Apply the reflector to X + apply_PX(mat_H, iu - 1, iu - 2, 2, n + 2 - iu, iu - 1); + apply_XP(mat_H, 0, iu - 1, il + bsize, 2, iu - 1); + + ref_nr(iu) = 1; + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::apply_PX(Mat<eT>& X, uword oi, uword oj, uword nrow, uword ncol, uword u_ind) + { + arma_extra_debug_sigprint(); + + if(ref_nr(u_ind) == 1) { return; } + + // Householder reflectors at index u_ind + Col<eT> u(ref_u.colptr(u_ind), 3, false); + + const uword stride = X.n_rows; + const eT u0_2 = 2 * u(0); + const eT u1_2 = 2 * u(1); + + eT* xptr = &X(oi, oj); + if(ref_nr(u_ind) == 2 || nrow == 2) + { + for(uword i = 0; i < ncol; i++, xptr += stride) + { + eT tmp = u0_2 * xptr[0] + u1_2 * xptr[1]; + xptr[0] -= tmp * u(0); + xptr[1] -= tmp * u(1); + } + } + else + { + const eT u2_2 = 2 * u(2); + for(uword i = 0; i < ncol; i++, xptr += stride) + { + eT tmp = u0_2 * xptr[0] + u1_2 * xptr[1] + u2_2 * xptr[2]; + xptr[0] -= tmp * u(0); + xptr[1] -= tmp * u(1); + xptr[2] -= tmp * u(2); + } + } + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::apply_PX(eT* x, uword u_ind) + { + arma_extra_debug_sigprint(); + + if(ref_nr(u_ind) == 1) { return; } + + eT u0 = ref_u(0, u_ind), + u1 = ref_u(1, u_ind), + u2 = ref_u(2, u_ind); + + // When the reflector only contains two elements, u2 has been set to zero + bool nr_is_2 = (ref_nr(u_ind) == 2); + eT dot2 = x[0] * u0 + x[1] * u1 + (nr_is_2 ? 0 : (x[2] * u2)); + dot2 *= 2; + x[0] -= dot2 * u0; + x[1] -= dot2 * u1; + if(!nr_is_2) { x[2] -= dot2 * u2; } + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::apply_XP(Mat<eT>& X, uword oi, uword oj, uword nrow, uword ncol, uword u_ind) + { + arma_extra_debug_sigprint(); + + if(ref_nr(u_ind) == 1) { return; } + + // Householder reflectors at index u_ind + Col<eT> u(ref_u.colptr(u_ind), 3, false); + uword stride = X.n_rows; + const eT u0_2 = 2 * u(0); + const eT u1_2 = 2 * u(1); + eT* X0 = &X(oi, oj); + eT* X1 = X0 + stride; // X0 => X(oi, oj), X1 => X(oi, oj + 1) + + if(ref_nr(u_ind) == 2 || ncol == 2) + { + // tmp = 2 * u0 * X0 + 2 * u1 * X1 + // X0 => X0 - u0 * tmp + // X1 => X1 - u1 * tmp + for(uword i = 0; i < nrow; i++) + { + eT tmp = u0_2 * X0[i] + u1_2 * X1[i]; + X0[i] -= tmp * u(0); + X1[i] -= tmp * u(1); + } + } + else + { + eT* X2 = X1 + stride; // X2 => X(oi, oj + 2) + const eT u2_2 = 2 * u(2); + for(uword i = 0; i < nrow; i++) + { + eT tmp = u0_2 * X0[i] + u1_2 * X1[i] + u2_2 * X2[i]; + X0[i] -= tmp * u(0); + X1[i] -= tmp * u(1); + X2[i] -= tmp * u(2); + } + } + } + + + +template<typename eT> +inline +DoubleShiftQR<eT>::DoubleShiftQR(uword size) + : n(size) + , prec(std::numeric_limits<eT>::epsilon()) + , eps_rel(prec) + , eps_abs(prec) + , computed(false) + { + arma_extra_debug_sigprint(); + } + + + +template<typename eT> +inline +DoubleShiftQR<eT>::DoubleShiftQR(const Mat<eT>& mat_obj, eT s, eT t) + : n(mat_obj.n_rows) + , mat_H(n, n) + , shift_s(s) + , shift_t(t) + , ref_u(3, n) + , ref_nr(n) + , prec(std::numeric_limits<eT>::epsilon()) + , eps_rel(prec) + , eps_abs(prec) + , computed(false) + { + arma_extra_debug_sigprint(); + + compute(mat_obj, s, t); + } + + + +template<typename eT> +void +DoubleShiftQR<eT>::compute(const Mat<eT>& mat_obj, eT s, eT t) + { + arma_extra_debug_sigprint(); + + arma_debug_check( (mat_obj.is_square() == false), "newarp::DoubleShiftQR::compute(): matrix must be square" ); + + n = mat_obj.n_rows; + mat_H.set_size(n, n); + shift_s = s; + shift_t = t; + ref_u.set_size(3, n); + ref_nr.set_size(n); + + // Make a copy of mat_obj + mat_H = mat_obj; + + // Obtain the indices of zero elements in the subdiagonal, + // so that H can be divided into several blocks + std::vector<uword> zero_ind; + zero_ind.reserve(n - 1); + zero_ind.push_back(0); + eT* Hii = mat_H.memptr(); + for(uword i = 0; i < n - 2; i++, Hii += (n + 1)) + { + // Hii[1] => mat_H(i + 1, i) + const eT h = std::abs(Hii[1]); + if(h <= eps_abs || h <= eps_rel * (std::abs(Hii[0]) + std::abs(Hii[n + 1]))) + { + Hii[1] = 0; + zero_ind.push_back(i + 1); + } + // Make sure mat_H is upper Hessenberg + // Zero the elements below mat_H(i + 1, i) + std::fill(Hii + 2, Hii + n - i, eT(0)); + } + zero_ind.push_back(n); + + for(std::vector<uword>::size_type i = 0; i < zero_ind.size() - 1; i++) + { + uword start = zero_ind[i]; + uword end = zero_ind[i + 1] - 1; + // Compute refelctors from each block X + update_block(start, end); + } + + computed = true; + } + + + +template<typename eT> +Mat<eT> +DoubleShiftQR<eT>::matrix_QtHQ() + { + arma_extra_debug_sigprint(); + + arma_debug_check( (computed == false), "newarp::DoubleShiftQR::matrix_QtHQ(): need to call compute() first" ); + + return mat_H; + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::apply_QtY(Col<eT>& y) + { + arma_extra_debug_sigprint(); + + arma_debug_check( (computed == false), "newarp::DoubleShiftQR::apply_QtY(): need to call compute() first" ); + + eT* y_ptr = y.memptr(); + for(uword i = 0; i < n - 1; i++, y_ptr++) + { + apply_PX(y_ptr, i); + } + } + + + +template<typename eT> +inline +void +DoubleShiftQR<eT>::apply_YQ(Mat<eT>& Y) + { + arma_extra_debug_sigprint(); + + arma_debug_check( (computed == false), "newarp::DoubleShiftQR::apply_YQ(): need to call compute() first" ); + + uword nrow = Y.n_rows; + for(uword i = 0; i < n - 2; i++) + { + apply_XP(Y, 0, i, nrow, 3, i); + } + + apply_XP(Y, 0, n - 2, nrow, 2, n - 2); + } + + +} // namespace newarp |