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