1 // SPDX-License-Identifier: Apache-2.0 2 // 3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) 4 // Copyright 2008-2016 National ICT Australia (NICTA) 5 // 6 // Licensed under the Apache License, Version 2.0 (the "License"); 7 // you may not use this file except in compliance with the License. 8 // You may obtain a copy of the License at 9 // http://www.apache.org/licenses/LICENSE-2.0 10 // 11 // Unless required by applicable law or agreed to in writing, software 12 // distributed under the License is distributed on an "AS IS" BASIS, 13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 // See the License for the specific language governing permissions and 15 // limitations under the License. 16 // ------------------------------------------------------------------------ 17 18 19 namespace newarp 20 { 21 22 23 template<typename eT> 24 class DoubleShiftQR 25 { 26 private: 27 28 uword n; // Dimension of the matrix 29 Mat<eT> mat_H; // A copy of the matrix to be factorised 30 eT shift_s; // Shift constant 31 eT shift_t; // Shift constant 32 Mat<eT> ref_u; // Householder reflectors 33 Col<unsigned short> ref_nr; // How many rows does each reflector affects 34 // 3 - A general reflector 35 // 2 - A Givens rotation 36 // 1 - An identity transformation 37 const eT prec; // Approximately zero 38 const eT eps_rel; 39 const eT eps_abs; 40 bool computed; // Whether matrix has been factorised 41 42 inline void compute_reflector(const eT& x1, const eT& x2, const eT& x3, uword ind); 43 arma_inline void compute_reflector(const eT* x, uword ind); 44 45 // Update the block X = H(il:iu, il:iu) 46 inline void update_block(uword il, uword iu); 47 48 // P = I - 2 * u * u' = P' 49 // PX = X - 2 * u * (u'X) 50 inline void apply_PX(Mat<eT>& X, uword oi, uword oj, uword nrow, uword ncol, uword u_ind); 51 52 // x is a pointer to a vector 53 // Px = x - 2 * dot(x, u) * u 54 inline void apply_PX(eT* x, uword u_ind); 55 56 // XP = X - 2 * (X * u) * u' 57 inline void apply_XP(Mat<eT>& X, uword oi, uword oj, uword nrow, uword ncol, uword u_ind); 58 59 60 public: 61 62 inline DoubleShiftQR(uword size); 63 64 inline DoubleShiftQR(const Mat<eT>& mat_obj, eT s, eT t); 65 66 inline void compute(const Mat<eT>& mat_obj, eT s, eT t); 67 68 inline Mat<eT> matrix_QtHQ(); 69 70 inline void apply_QtY(Col<eT>& y); 71 72 inline void apply_YQ(Mat<eT>& Y); 73 }; 74 75 76 } // namespace newarp 77