GetFEM  5.4.3
gmm_dense_Householder.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2020 Yves Renard, Caroline Lecalvez
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_dense_Householder.h
33  @author Caroline Lecalvez <Caroline.Lecalvez@gmm.insa-toulouse.fr>
34  @author Yves Renard <Yves.Renard@insa-lyon.fr>
35  @date June 5, 2003.
36  @brief Householder for dense matrices.
37 */
38 
39 #ifndef GMM_DENSE_HOUSEHOLDER_H
40 #define GMM_DENSE_HOUSEHOLDER_H
41 
42 #include "gmm_kernel.h"
43 
44 namespace gmm {
45  ///@cond DOXY_SHOW_ALL_FUNCTIONS
46 
47  /* ********************************************************************* */
48  /* Rank one update (complex and real version) */
49  /* ********************************************************************* */
50 
51  template <typename Matrix, typename VecX, typename VecY>
52  inline void rank_one_update(Matrix &A, const VecX& x,
53  const VecY& y, row_major) {
54  typedef typename linalg_traits<Matrix>::value_type T;
55  size_type N = mat_nrows(A);
56  GMM_ASSERT2(N <= vect_size(x) && mat_ncols(A) <= vect_size(y),
57  "dimensions mismatch");
58  typename linalg_traits<VecX>::const_iterator itx = vect_const_begin(x);
59  for (size_type i = 0; i < N; ++i, ++itx) {
60  typedef typename linalg_traits<Matrix>::sub_row_type row_type;
61  row_type row = mat_row(A, i);
62  typename linalg_traits<typename org_type<row_type>::t>::iterator
63  it = vect_begin(row), ite = vect_end(row);
64  typename linalg_traits<VecY>::const_iterator ity = vect_const_begin(y);
65  T tx = *itx;
66  for (; it != ite; ++it, ++ity) *it += conj_product(*ity, tx);
67  }
68  }
69 
70  template <typename Matrix, typename VecX, typename VecY>
71  inline void rank_one_update(Matrix &A, const VecX& x,
72  const VecY& y, col_major) {
73  typedef typename linalg_traits<Matrix>::value_type T;
74  size_type M = mat_ncols(A);
75  GMM_ASSERT2(mat_nrows(A) <= vect_size(x) && M <= vect_size(y),
76  "dimensions mismatch");
77  typename linalg_traits<VecY>::const_iterator ity = vect_const_begin(y);
78  for (size_type i = 0; i < M; ++i, ++ity) {
79  typedef typename linalg_traits<Matrix>::sub_col_type col_type;
80  col_type col = mat_col(A, i);
81  typename linalg_traits<typename org_type<col_type>::t>::iterator
82  it = vect_begin(col), ite = vect_end(col);
83  typename linalg_traits<VecX>::const_iterator itx = vect_const_begin(x);
84  T ty = *ity;
85  for (; it != ite; ++it, ++itx) *it += conj_product(ty, *itx);
86  }
87  }
88 
89  ///@endcond
90  template <typename Matrix, typename VecX, typename VecY>
91  inline void rank_one_update(const Matrix &AA, const VecX& x,
92  const VecY& y) {
93  Matrix& A = const_cast<Matrix&>(AA);
94  rank_one_update(A, x, y, typename principal_orientation_type<typename
95  linalg_traits<Matrix>::sub_orientation>::potype());
96  }
97  ///@cond DOXY_SHOW_ALL_FUNCTIONS
98 
99  /* ********************************************************************* */
100  /* Rank two update (complex and real version) */
101  /* ********************************************************************* */
102 
103  template <typename Matrix, typename VecX, typename VecY>
104  inline void rank_two_update(Matrix &A, const VecX& x,
105  const VecY& y, row_major) {
106  typedef typename linalg_traits<Matrix>::value_type value_type;
107  size_type N = mat_nrows(A);
108  GMM_ASSERT2(N <= vect_size(x) && mat_ncols(A) <= vect_size(y),
109  "dimensions mismatch");
110  typename linalg_traits<VecX>::const_iterator itx1 = vect_const_begin(x);
111  typename linalg_traits<VecY>::const_iterator ity2 = vect_const_begin(y);
112  for (size_type i = 0; i < N; ++i, ++itx1, ++ity2) {
113  typedef typename linalg_traits<Matrix>::sub_row_type row_type;
114  row_type row = mat_row(A, i);
115  typename linalg_traits<typename org_type<row_type>::t>::iterator
116  it = vect_begin(row), ite = vect_end(row);
117  typename linalg_traits<VecX>::const_iterator itx2 = vect_const_begin(x);
118  typename linalg_traits<VecY>::const_iterator ity1 = vect_const_begin(y);
119  value_type tx = *itx1, ty = *ity2;
120  for (; it != ite; ++it, ++ity1, ++itx2)
121  *it += conj_product(*ity1, tx) + conj_product(*itx2, ty);
122  }
123  }
124 
125  template <typename Matrix, typename VecX, typename VecY>
126  inline void rank_two_update(Matrix &A, const VecX& x,
127  const VecY& y, col_major) {
128  typedef typename linalg_traits<Matrix>::value_type value_type;
129  size_type M = mat_ncols(A);
130  GMM_ASSERT2(mat_nrows(A) <= vect_size(x) && M <= vect_size(y),
131  "dimensions mismatch");
132  typename linalg_traits<VecX>::const_iterator itx2 = vect_const_begin(x);
133  typename linalg_traits<VecY>::const_iterator ity1 = vect_const_begin(y);
134  for (size_type i = 0; i < M; ++i, ++ity1, ++itx2) {
135  typedef typename linalg_traits<Matrix>::sub_col_type col_type;
136  col_type col = mat_col(A, i);
137  typename linalg_traits<typename org_type<col_type>::t>::iterator
138  it = vect_begin(col), ite = vect_end(col);
139  typename linalg_traits<VecX>::const_iterator itx1 = vect_const_begin(x);
140  typename linalg_traits<VecY>::const_iterator ity2 = vect_const_begin(y);
141  value_type ty = *ity1, tx = *itx2;
142  for (; it != ite; ++it, ++itx1, ++ity2)
143  *it += conj_product(ty, *itx1) + conj_product(tx, *ity2);
144  }
145  }
146 
147  ///@endcond
148  template <typename Matrix, typename VecX, typename VecY>
149  inline void rank_two_update(const Matrix &AA, const VecX& x,
150  const VecY& y) {
151  Matrix& A = const_cast<Matrix&>(AA);
152  rank_two_update(A, x, y, typename principal_orientation_type<typename
153  linalg_traits<Matrix>::sub_orientation>::potype());
154  }
155  ///@cond DOXY_SHOW_ALL_FUNCTIONS
156 
157  /* ********************************************************************* */
158  /* Householder vector computation (complex and real version) */
159  /* ********************************************************************* */
160 
161  template <typename VECT> void house_vector(const VECT &VV) {
162  VECT &V = const_cast<VECT &>(VV);
163  typedef typename linalg_traits<VECT>::value_type T;
164  typedef typename number_traits<T>::magnitude_type R;
165 
166  R mu = vect_norm2(V), abs_v0 = gmm::abs(V[0]);
167  if (mu != R(0))
168  gmm::scale(V, (abs_v0 == R(0)) ? T(R(1) / mu)
169  : (safe_divide(T(abs_v0), V[0]) / (abs_v0 + mu)));
170  if (gmm::real(V[vect_size(V)-1]) * R(0) != R(0))
171  gmm::clear(V);
172  V[0] = T(1);
173  }
174 
175  template <typename VECT> void house_vector_last(const VECT &VV) {
176  VECT &V = const_cast<VECT &>(VV);
177  typedef typename linalg_traits<VECT>::value_type T;
178  typedef typename number_traits<T>::magnitude_type R;
179 
180  size_type m = vect_size(V);
181  R mu = vect_norm2(V), abs_v0 = gmm::abs(V[m-1]);
182  if (mu != R(0))
183  gmm::scale(V, (abs_v0 == R(0)) ? T(R(1) / mu)
184  : ((abs_v0 / V[m-1]) / (abs_v0 + mu)));
185  if (gmm::real(V[0]) * R(0) != R(0))
186  gmm::clear(V);
187  V[m-1] = T(1);
188  }
189 
190  /* ********************************************************************* */
191  /* Householder updates (complex and real version) */
192  /* ********************************************************************* */
193 
194  // multiply A to the left by the reflector stored in V. W is a temporary.
195  template <typename MAT, typename VECT1, typename VECT2> inline
196  void row_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
197  VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
198  typedef typename linalg_traits<MAT>::value_type value_type;
199  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
200 
201  gmm::mult(conjugated(A),
202  scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
203  rank_one_update(A, V, W);
204  }
205 
206  // multiply A to the right by the reflector stored in V. W is a temporary.
207  template <typename MAT, typename VECT1, typename VECT2> inline
208  void col_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
209  VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
210  typedef typename linalg_traits<MAT>::value_type value_type;
211  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
212 
213  gmm::mult(A,
214  scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
215  rank_one_update(A, W, V);
216  }
217 
218  ///@endcond
219 
220  /* ********************************************************************* */
221  /* Hessenberg reduction with Householder. */
222  /* ********************************************************************* */
223 
224  template <typename MAT1, typename MAT2>
225  void Hessenberg_reduction(const MAT1& AA, const MAT2 &QQ, bool compute_Q){
226  MAT1& A = const_cast<MAT1&>(AA); MAT2& Q = const_cast<MAT2&>(QQ);
227  typedef typename linalg_traits<MAT1>::value_type value_type;
228  if (compute_Q) gmm::copy(identity_matrix(), Q);
229  size_type n = mat_nrows(A); if (n < 2) return;
230  std::vector<value_type> v(n), w(n);
231  sub_interval SUBK(0,n);
232  for (size_type k = 1; k+1 < n; ++k) {
233  sub_interval SUBI(k, n-k), SUBJ(k-1,n-k+1);
234  v.resize(n-k);
235  for (size_type j = k; j < n; ++j) v[j-k] = A(j, k-1);
236  house_vector(v);
237  row_house_update(sub_matrix(A, SUBI, SUBJ), v, sub_vector(w, SUBJ));
238  col_house_update(sub_matrix(A, SUBK, SUBI), v, w);
239  // is it possible to "unify" the two on the common part of the matrix?
240  if (compute_Q) col_house_update(sub_matrix(Q, SUBK, SUBI), v, w);
241  }
242  }
243 
244  /* ********************************************************************* */
245  /* Householder tridiagonalization for symmetric matrices */
246  /* ********************************************************************* */
247 
248  template <typename MAT1, typename MAT2>
249  void Householder_tridiagonalization(const MAT1 &AA, const MAT2 &QQ,
250  bool compute_Q) {
251  MAT1 &A = const_cast<MAT1 &>(AA); MAT2 &Q = const_cast<MAT2 &>(QQ);
252  typedef typename linalg_traits<MAT1>::value_type T;
253  typedef typename number_traits<T>::magnitude_type R;
254 
255  size_type n = mat_nrows(A); if (n < 2) return;
256  std::vector<T> v(n), p(n), w(n), ww(n);
257  sub_interval SUBK(0,n);
258 
259  for (size_type k = 1; k+1 < n; ++k) { // not optimized ...
260  sub_interval SUBI(k, n-k);
261  v.resize(n-k); p.resize(n-k); w.resize(n-k);
262  for (size_type l = k; l < n; ++l)
263  { v[l-k] = w[l-k] = A(l, k-1); A(l, k-1) = A(k-1, l) = T(0); }
264  house_vector(v);
265  R norm = vect_norm2_sqr(v);
266  A(k-1, k) = gmm::conj(A(k, k-1) = w[0] - T(2)*v[0]*vect_hp(w, v)/norm);
267 
268  gmm::mult(sub_matrix(A, SUBI), gmm::scaled(v, T(-2) / norm), p);
269  gmm::add(p, gmm::scaled(v, -vect_hp(v, p) / norm), w);
270  rank_two_update(sub_matrix(A, SUBI), v, w);
271  // it should be possible to compute only the upper or lower part
272  if (compute_Q)
273  col_house_update(sub_matrix(Q, SUBK, SUBI), v, ww);
274  }
275  }
276 
277  /* ********************************************************************* */
278  /* Real and complex Givens rotations */
279  /* ********************************************************************* */
280 
281  template <typename T> void Givens_rotation(T a, T b, T &c, T &s) {
282  typedef typename number_traits<T>::magnitude_type R;
283  R aa = gmm::abs(a), bb = gmm::abs(b);
284  if (bb == R(0)) { c = T(1); s = T(0); return; }
285  if (aa == R(0)) { c = T(0); s = b / bb; return; }
286  if (bb > aa)
287  { T t = -safe_divide(a,b); s = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); c = s * t; }
288  else
289  { T t = -safe_divide(b,a); c = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); s = c * t; }
290  }
291 
292  // Apply Q* v
293  template <typename T> inline
294  void Apply_Givens_rotation_left(T &x, T &y, T c, T s)
295  { T t1=x, t2=y; x = gmm::conj(c)*t1 - gmm::conj(s)*t2; y = c*t2 + s*t1; }
296 
297  // Apply v^T Q
298  template <typename T> inline
299  void Apply_Givens_rotation_right(T &x, T &y, T c, T s)
300  { T t1=x, t2=y; x = c*t1 - s*t2; y = gmm::conj(c)*t2 + gmm::conj(s)*t1; }
301 
302  template <typename MAT, typename T>
303  void row_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
304  MAT &A = const_cast<MAT &>(AA); // can be specialized for row matrices
305  for (size_type j = 0; j < mat_ncols(A); ++j)
306  Apply_Givens_rotation_left(A(i,j), A(k,j), c, s);
307  }
308 
309  template <typename MAT, typename T>
310  void col_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
311  MAT &A = const_cast<MAT &>(AA); // can be specialized for column matrices
312  for (size_type j = 0; j < mat_nrows(A); ++j)
313  Apply_Givens_rotation_right(A(j,i), A(j,k), c, s);
314  }
315 
316 }
317 
318 #endif
319 
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:558
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:512
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:545
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
gmm::conjugated
conjugated_return< L >::return_type conjugated(const L &v)
return a conjugated view of the input matrix or vector.
Definition: gmm_conjugated.h:302
gmm_kernel.h
Include the base gmm files.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49