GetFEM  5.4.2
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)) 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)) gmm::clear(V);
185  V[m-1] = T(1);
186  }
187 
188  /* ********************************************************************* */
189  /* Householder updates (complex and real version) */
190  /* ********************************************************************* */
191 
192  // multiply A to the left by the reflector stored in V. W is a temporary.
193  template <typename MAT, typename VECT1, typename VECT2> inline
194  void row_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
195  VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
196  typedef typename linalg_traits<MAT>::value_type value_type;
197  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
198 
199  gmm::mult(conjugated(A),
200  scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
201  rank_one_update(A, V, W);
202  }
203 
204  // multiply A to the right by the reflector stored in V. W is a temporary.
205  template <typename MAT, typename VECT1, typename VECT2> inline
206  void col_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
207  VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
208  typedef typename linalg_traits<MAT>::value_type value_type;
209  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
210 
211  gmm::mult(A,
212  scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
213  rank_one_update(A, W, V);
214  }
215 
216  ///@endcond
217 
218  /* ********************************************************************* */
219  /* Hessenberg reduction with Householder. */
220  /* ********************************************************************* */
221 
222  template <typename MAT1, typename MAT2>
223  void Hessenberg_reduction(const MAT1& AA, const MAT2 &QQ, bool compute_Q){
224  MAT1& A = const_cast<MAT1&>(AA); MAT2& Q = const_cast<MAT2&>(QQ);
225  typedef typename linalg_traits<MAT1>::value_type value_type;
226  if (compute_Q) gmm::copy(identity_matrix(), Q);
227  size_type n = mat_nrows(A); if (n < 2) return;
228  std::vector<value_type> v(n), w(n);
229  sub_interval SUBK(0,n);
230  for (size_type k = 1; k+1 < n; ++k) {
231  sub_interval SUBI(k, n-k), SUBJ(k-1,n-k+1);
232  v.resize(n-k);
233  for (size_type j = k; j < n; ++j) v[j-k] = A(j, k-1);
234  house_vector(v);
235  row_house_update(sub_matrix(A, SUBI, SUBJ), v, sub_vector(w, SUBJ));
236  col_house_update(sub_matrix(A, SUBK, SUBI), v, w);
237  // is it possible to "unify" the two on the common part of the matrix?
238  if (compute_Q) col_house_update(sub_matrix(Q, SUBK, SUBI), v, w);
239  }
240  }
241 
242  /* ********************************************************************* */
243  /* Householder tridiagonalization for symmetric matrices */
244  /* ********************************************************************* */
245 
246  template <typename MAT1, typename MAT2>
247  void Householder_tridiagonalization(const MAT1 &AA, const MAT2 &QQ,
248  bool compute_q) {
249  MAT1 &A = const_cast<MAT1 &>(AA); MAT2 &Q = const_cast<MAT2 &>(QQ);
250  typedef typename linalg_traits<MAT1>::value_type T;
251  typedef typename number_traits<T>::magnitude_type R;
252 
253  size_type n = mat_nrows(A); if (n < 2) return;
254  std::vector<T> v(n), p(n), w(n), ww(n);
255  sub_interval SUBK(0,n);
256 
257  for (size_type k = 1; k+1 < n; ++k) { // not optimized ...
258  sub_interval SUBI(k, n-k);
259  v.resize(n-k); p.resize(n-k); w.resize(n-k);
260  for (size_type l = k; l < n; ++l)
261  { v[l-k] = w[l-k] = A(l, k-1); A(l, k-1) = A(k-1, l) = T(0); }
262  house_vector(v);
263  R norm = vect_norm2_sqr(v);
264  A(k-1, k) = gmm::conj(A(k, k-1) = w[0] - T(2)*v[0]*vect_hp(w, v)/norm);
265 
266  gmm::mult(sub_matrix(A, SUBI), gmm::scaled(v, T(-2) / norm), p);
267  gmm::add(p, gmm::scaled(v, -vect_hp(v, p) / norm), w);
268  rank_two_update(sub_matrix(A, SUBI), v, w);
269  // it should be possible to compute only the upper or lower part
270 
271  if (compute_q) col_house_update(sub_matrix(Q, SUBK, SUBI), v, ww);
272  }
273  }
274 
275  /* ********************************************************************* */
276  /* Real and complex Givens rotations */
277  /* ********************************************************************* */
278 
279  template <typename T> void Givens_rotation(T a, T b, T &c, T &s) {
280  typedef typename number_traits<T>::magnitude_type R;
281  R aa = gmm::abs(a), bb = gmm::abs(b);
282  if (bb == R(0)) { c = T(1); s = T(0); return; }
283  if (aa == R(0)) { c = T(0); s = b / bb; return; }
284  if (bb > aa)
285  { T t = -safe_divide(a,b); s = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); c = s * t; }
286  else
287  { T t = -safe_divide(b,a); c = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); s = c * t; }
288  }
289 
290  // Apply Q* v
291  template <typename T> inline
292  void Apply_Givens_rotation_left(T &x, T &y, T c, T s)
293  { T t1=x, t2=y; x = gmm::conj(c)*t1 - gmm::conj(s)*t2; y = c*t2 + s*t1; }
294 
295  // Apply v^T Q
296  template <typename T> inline
297  void Apply_Givens_rotation_right(T &x, T &y, T c, T s)
298  { T t1=x, t2=y; x = c*t1 - s*t2; y = gmm::conj(c)*t2 + gmm::conj(s)*t1; }
299 
300  template <typename MAT, typename T>
301  void row_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
302  MAT &A = const_cast<MAT &>(AA); // can be specialized for row matrices
303  for (size_type j = 0; j < mat_ncols(A); ++j)
304  Apply_Givens_rotation_left(A(i,j), A(k,j), c, s);
305  }
306 
307  template <typename MAT, typename T>
308  void col_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
309  MAT &A = const_cast<MAT &>(AA); // can be specialized for column matrices
310  for (size_type j = 0; j < mat_nrows(A); ++j)
311  Apply_Givens_rotation_right(A(j,i), A(j,k), c, s);
312  }
313 
314 }
315 
316 #endif
317 
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:511
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:544
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:557
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:294
gmm_kernel.h
Include the base gmm files.