GetFEM  5.5
gmm_dense_qr.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2026 Yves Renard
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_qr.h
32  @author Caroline Lecalvez, Caroline.Lecalvez@gmm.insa-tlse.fr, Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date September 12, 2003.
34  @brief Dense QR factorization.
35 */
36 #ifndef GMM_DENSE_QR_H
37 #define GMM_DENSE_QR_H
38 
39 #include "gmm_dense_Householder.h"
40 
41 namespace gmm {
42 
43 
44  /**
45  QR factorization using Householder method (complex and real version).
46  */
47  template <typename MAT>
48  void qr_factor(const MAT &A_) {
49  MAT &A = const_cast<MAT &>(A_);
50  typedef typename linalg_traits<MAT>::value_type T;
51  typedef typename number_traits<T>::magnitude_type R;
52 
53  const size_type m = mat_nrows(A), n = mat_ncols(A);
54  GMM_ASSERT2(m >= n, "dimensions mismatch");
55 
56  std::vector<T> W(m), V(m);
57  // skip the last reflection for *real* square matrices (m-1 limit)
58  // lapack does the same in {sd}larfg.f but not in {cz}larfg.f
59  const size_type jmax = (m==n && !is_complex(T())) ? m-1
60  : n;
61  for (size_type j = 0; j < jmax; ++j) {
62  sub_interval SUBI(j, m-j), SUBJ(j, n-j);
63  V.resize(m-j); W.resize(n-j);
64 
65  for (size_type i = j; i < m; ++i) V[i-j] = A(i, j);
66  house_vector(V);
67 
68  R Vnorm2 = vect_norm2_sqr(V);
69  if (Vnorm2 > R(1) || gmm::imag(A(j, j)) != R(0)) { // compatible with {sdcz}larfg.f
70  auto subA = sub_matrix(A, SUBI, SUBJ);
71  gmm::mult(conjugated(subA),
72  scaled(V, T(R(-2)/Vnorm2)), W);
73  rank_one_update(subA, V, W);
74  for (size_type i = j+1; i < m; ++i) A(i, j) = V[i-j];
75  }
76  }
77  }
78 
79 
80  // QR comes from QR_factor(QR) where the upper triangular part stands for R
81  // and the lower part contains the Householder reflectors.
82  // A <- AQ
83  template <typename MAT1, typename MAT2>
84  void apply_house_right(const MAT1 &QR, const MAT2 &A_) {
85  MAT2 &A = const_cast<MAT2 &>(A_);
86  typedef typename linalg_traits<MAT1>::value_type T;
87  typedef typename number_traits<T>::magnitude_type R;
88 
89  const size_type m = mat_nrows(QR), n = mat_ncols(QR);
90  GMM_ASSERT2(m == mat_ncols(A), "dimensions mismatch");
91  if (m == 0) return;
92  std::vector<T> V(m), W(mat_nrows(A));
93  V[0] = T(1);
94  // assume that the last reflection was skipped for *real* square matrices
95  const size_type jmax = (m==n && !is_complex(T())) ? m-1
96  : n;
97  for (size_type j = 0; j < jmax; ++j) {
98  V.resize(m-j);
99  for (size_type i = j+1; i < m; ++i)
100  V[i-j] = QR(i, j);
101  R Vnorm2 = vect_norm2_sqr(V);
102  if (Vnorm2 > R(1) || gmm::imag(QR(j, j)) != R(0)) { // compatible with {sdcz}larfg.f
103  auto subA = sub_matrix(A, sub_interval(0, mat_nrows(A)),
104  sub_interval(j, m-j));
105  gmm::mult(subA, scaled(V, T(R(-2)/Vnorm2)), W);
106  rank_one_update(subA, W, V);
107  }
108  }
109  }
110 
111  // QR comes from QR_factor(QR) where the upper triangular part stands for R
112  // and the lower part contains the Householder reflectors.
113  // A <- Q*A
114  template <typename MAT1, typename MAT2>
115  void apply_house_left(const MAT1 &QR, const MAT2 &A_) {
116  MAT2 &A = const_cast<MAT2 &>(A_);
117  typedef typename linalg_traits<MAT1>::value_type T;
118  typedef typename number_traits<T>::magnitude_type R;
119 
120  const size_type m = mat_nrows(QR), n = mat_ncols(QR);
121  GMM_ASSERT2(m == mat_nrows(A), "dimensions mismatch");
122  if (m == 0) return;
123  std::vector<T> V(m), W(mat_ncols(A));
124  V[0] = T(1);
125  // assume that the last reflection was skipped for *real* square matrices
126  const size_type jmax = (m==n && !is_complex(T())) ? m-1
127  : n;
128  for (size_type j = 0; j < jmax; ++j) {
129  V.resize(m-j);
130  for (size_type i = j+1; i < m; ++i)
131  V[i-j] = QR(i, j);
132  R Vnorm2 = vect_norm2_sqr(V);
133  if (Vnorm2 > R(1) || gmm::imag(QR(j, j)) != R(0)) { // compatible with {sdcz}larfg.f
134  auto subA = sub_matrix(A, sub_interval(j, m-j),
135  sub_interval(0, mat_ncols(A)));
136  gmm::mult(conjugated(subA),
137  scaled(V, T(R(-2)/Vnorm2)), W);
138  rank_one_update(subA, V, W);
139  }
140  }
141  }
142 
143  /** Compute the QR factorization, where Q is assembled. */
144  template <typename MAT1, typename MAT2, typename MAT3>
145  void qr_factor(const MAT1 &A, const MAT2 &QQ, const MAT3 &RR) {
146  MAT2 &Q = const_cast<MAT2 &>(QQ); MAT3 &R = const_cast<MAT3 &>(RR);
147  typedef typename linalg_traits<MAT1>::value_type T;
148 
149  const size_type m = mat_nrows(A), n = mat_ncols(A);
150  GMM_ASSERT2(m >= n, "dimensions mismatch");
151  gmm::copy(A, Q);
152 
153  std::vector<T> W(m);
154  dense_matrix<T> VV(m, n);
155  // skip the last reflection for *real* square matrices (m-1 limit)
156  // lapack does the same in dlarfg.f but not in zlarfg.f
157  const size_type jmax = (m==n && !is_complex(T())) ? m-1
158  : n;
159  for (size_type j = 0; j < jmax; ++j) {
160  sub_interval SUBI(j, m-j), SUBJ(j, n-j);
161 
162  for (size_type i = j; i < m; ++i) VV(i,j) = Q(i, j);
163  house_vector(sub_vector(mat_col(VV,j), SUBI));
164 
165  row_house_update(sub_matrix(Q, SUBI, SUBJ),
166  sub_vector(mat_col(VV,j), SUBI), sub_vector(W, SUBJ));
167  }
168 
169  gmm::copy(sub_matrix(Q, sub_interval(0, n), sub_interval(0, n)), R);
170  gmm::copy(identity_matrix(), Q);
171  // assume that the last reflection was skipped for *real* square matrices
172  const size_type jstart = (m==n && !is_complex(T())) ? m-2
173  : n-1;
174  for (size_type j = jstart; j != size_type(-1); --j) {
175  sub_interval SUBI(j, m-j), SUBJ(j, n-j);
176  row_house_update(sub_matrix(Q, SUBI, SUBJ),
177  sub_vector(mat_col(VV,j), SUBI), sub_vector(W, SUBJ));
178  }
179  }
180 
181  ///@cond DOXY_SHOW_ALL_FUNCTIONS
182  template <typename TA, typename TV, typename Ttol,
183  typename MAT, typename VECT>
184  void extract_eig(const MAT &A, VECT &V, Ttol tol, TA, TV) {
185  size_type n = mat_nrows(A);
186  if (n == 0) return;
187  tol *= Ttol(2);
188  Ttol tol_i = tol * gmm::abs(A(0,0)), tol_cplx = tol_i;
189  for (size_type i = 0; i < n; ++i) {
190  if (i < n-1) {
191  tol_i = (gmm::abs(A(i,i))+gmm::abs(A(i+1,i+1)))*tol;
192  tol_cplx = std::max(tol_cplx, tol_i);
193  }
194  if ((i < n-1) && gmm::abs(A(i+1,i)) >= tol_i) {
195  TA tr = A(i,i) + A(i+1, i+1);
196  TA det = A(i,i)*A(i+1, i+1) - A(i,i+1)*A(i+1, i);
197  TA delta = tr*tr - TA(4) * det;
198  if (delta < -tol_cplx) {
199  GMM_WARNING1("A complex eigenvalue has been detected : "
200  << std::complex<TA>(tr/TA(2), gmm::sqrt(-delta)/TA(2)));
201  V[i] = V[i+1] = tr / TA(2);
202  }
203  else {
204  delta = std::max(TA(0), delta);
205  V[i ] = TA(tr + gmm::sqrt(delta))/ TA(2);
206  V[i+1] = TA(tr - gmm::sqrt(delta))/ TA(2);
207  }
208  ++i;
209  }
210  else
211  V[i] = TV(A(i,i));
212  }
213  }
214 
215  template <typename TA, typename TV, typename Ttol,
216  typename MAT, typename VECT>
217  void extract_eig(const MAT &A, VECT &V, Ttol tol, TA, std::complex<TV>) {
218  size_type n = mat_nrows(A);
219  tol *= Ttol(2);
220  for (size_type i = 0; i < n; ++i)
221  if ((i == n-1) ||
222  gmm::abs(A(i+1,i)) < (gmm::abs(A(i,i))+gmm::abs(A(i+1,i+1)))*tol)
223  V[i] = std::complex<TV>(A(i,i));
224  else {
225  TA tr = A(i,i) + A(i+1, i+1);
226  TA det = A(i,i)*A(i+1, i+1) - A(i,i+1)*A(i+1, i);
227  TA delta = tr*tr - TA(4) * det;
228  if (delta < TA(0)) {
229  V[i] = std::complex<TV>(tr / TA(2), gmm::sqrt(-delta) / TA(2));
230  V[i+1] = std::complex<TV>(tr / TA(2), -gmm::sqrt(-delta)/ TA(2));
231  }
232  else {
233  V[i ] = TA(tr + gmm::sqrt(delta)) / TA(2);
234  V[i+1] = TA(tr - gmm::sqrt(delta)) / TA(2);
235  }
236  ++i;
237  }
238  }
239 
240  template <typename TA, typename TV, typename Ttol,
241  typename MAT, typename VECT>
242  void extract_eig(const MAT &A, VECT &V, Ttol tol, std::complex<TA>, TV) {
243  typedef std::complex<TA> T;
244  size_type n = mat_nrows(A);
245  if (n == 0) return;
246  tol *= Ttol(2);
247  Ttol tol_i = tol * gmm::abs(A(0,0)), tol_cplx = tol_i;
248  for (size_type i = 0; i < n; ++i) {
249  if (i < n-1) {
250  tol_i = (gmm::abs(A(i,i))+gmm::abs(A(i+1,i+1)))*tol;
251  tol_cplx = std::max(tol_cplx, tol_i);
252  }
253  if ((i == n-1) || gmm::abs(A(i+1,i)) < tol_i) {
254  if (gmm::abs(std::imag(A(i,i))) > tol_cplx)
255  GMM_WARNING1("A complex eigenvalue has been detected : "
256  << T(A(i,i)) << " : " << gmm::abs(std::imag(A(i,i)))
257  / gmm::abs(std::real(A(i,i))) << " : " << tol_cplx);
258  V[i] = std::real(A(i,i));
259  }
260  else {
261  T tr = A(i,i) + A(i+1, i+1);
262  T det = A(i,i)*A(i+1, i+1) - A(i,i+1)*A(i+1, i);
263  T delta = tr*tr - TA(4) * det;
264  T a1 = (tr + gmm::sqrt(delta)) / TA(2);
265  T a2 = (tr - gmm::sqrt(delta)) / TA(2);
266  if (gmm::abs(std::imag(a1)) > tol_cplx)
267  GMM_WARNING1("A complex eigenvalue has been detected : " << a1);
268  if (gmm::abs(std::imag(a2)) > tol_cplx)
269  GMM_WARNING1("A complex eigenvalue has been detected : " << a2);
270 
271  V[i] = std::real(a1); V[i+1] = std::real(a2);
272  ++i;
273  }
274  }
275  }
276 
277  template <typename TA, typename TV, typename Ttol,
278  typename MAT, typename VECT>
279  void extract_eig(const MAT &A, VECT &V, Ttol tol,
280  std::complex<TA>, std::complex<TV>) {
281  size_type n = mat_nrows(A);
282  tol *= Ttol(2);
283  for (size_type i = 0; i < n; ++i)
284  if ((i == n-1) ||
285  gmm::abs(A(i+1,i)) < (gmm::abs(A(i,i))+gmm::abs(A(i+1,i+1)))*tol)
286  V[i] = std::complex<TV>(A(i,i));
287  else {
288  std::complex<TA> tr = A(i,i) + A(i+1, i+1);
289  std::complex<TA> det = A(i,i)*A(i+1, i+1) - A(i,i+1)*A(i+1, i);
290  std::complex<TA> delta = tr*tr - TA(4) * det;
291  V[i] = (tr + gmm::sqrt(delta)) / TA(2);
292  V[i+1] = (tr - gmm::sqrt(delta)) / TA(2);
293  ++i;
294  }
295  }
296 
297  ///@endcond
298  /**
299  Compute eigenvalue vector.
300  */
301  template <typename MAT, typename Ttol, typename VECT> inline
302  void extract_eig(const MAT &A, const VECT &V, Ttol tol) {
303  extract_eig(A, const_cast<VECT&>(V), tol,
304  typename linalg_traits<MAT>::value_type(),
305  typename linalg_traits<VECT>::value_type());
306  }
307 
308  /* ********************************************************************* */
309  /* Stop criterion for QR algorithms */
310  /* ********************************************************************* */
311 
312  template <typename MAT, typename Ttol>
313  void qr_stop_criterion(MAT &A, size_type &p, size_type &q, Ttol tol) {
314  typedef typename linalg_traits<MAT>::value_type T;
315  typedef typename number_traits<T>::magnitude_type R;
316  R rmin = default_min(R()) * R(2);
317  size_type n = mat_nrows(A);
318  if (n <= 2) { q = n; p = 0; }
319  else {
320  for (size_type i = 1; i < n-q; ++i)
321  if (gmm::abs(A(i,i-1)) < (gmm::abs(A(i,i))+ gmm::abs(A(i-1,i-1)))*tol
322  || gmm::abs(A(i,i-1)) < rmin)
323  A(i,i-1) = T(0);
324 
325  while ((q < n-1 && A(n-1-q, n-2-q) == T(0)) ||
326  (q < n-2 && A(n-2-q, n-3-q) == T(0))) ++q;
327  if (q >= n-2) q = n;
328  p = n-q; if (p) --p; if (p) --p;
329  while (p > 0 && A(p,p-1) != T(0)) --p;
330  }
331  }
332 
333  template <typename MAT, typename Ttol> inline
334  void symmetric_qr_stop_criterion(const MAT &AA, size_type &p, size_type &q,
335  Ttol tol) {
336  typedef typename linalg_traits<MAT>::value_type T;
337  typedef typename number_traits<T>::magnitude_type R;
338  R rmin = default_min(R()) * R(2);
339  MAT& A = const_cast<MAT&>(AA);
340  size_type n = mat_nrows(A);
341  if (n <= 1) { q = n; p = 0; }
342  else {
343  for (size_type i = 1; i < n-q; ++i)
344  if (gmm::abs(A(i,i-1)) < (gmm::abs(A(i,i))+ gmm::abs(A(i-1,i-1)))*tol
345  || gmm::abs(A(i,i-1)) < rmin)
346  A(i,i-1) = T(0);
347 
348  while (q < n-1 && A(n-1-q, n-2-q) == T(0)) ++q;
349  if (q >= n-1) q = n;
350  p = n-q; if (p) --p; if (p) --p;
351  while (p > 0 && A(p,p-1) != T(0)) --p;
352  }
353  }
354 
355  template <typename VECT1, typename VECT2, typename Ttol> inline
356  void symmetric_qr_stop_criterion(const VECT1 &diag, const VECT2 &sdiag_,
357  size_type &p, size_type &q, Ttol tol) {
358  typedef typename linalg_traits<VECT2>::value_type T;
359  typedef typename number_traits<T>::magnitude_type R;
360  R rmin = default_min(R()) * R(2);
361  VECT2 &sdiag = const_cast<VECT2 &>(sdiag_);
362  size_type n = vect_size(diag);
363  if (n <= 1) { q = n; p = 0; return; }
364  for (size_type i = 1; i < n-q; ++i)
365  if (gmm::abs(sdiag[i-1]) < (gmm::abs(diag[i])+ gmm::abs(diag[i-1]))*tol
366  || gmm::abs(sdiag[i-1]) < rmin)
367  sdiag[i-1] = T(0);
368  while (q < n-1 && sdiag[n-2-q] == T(0)) ++q;
369  if (q >= n-1) q = n;
370  p = n-q; if (p) --p; if (p) --p;
371  while (p > 0 && sdiag[p-1] != T(0)) --p;
372  }
373 
374  /* ********************************************************************* */
375  /* 2x2 blocks reduction for Schur vectors */
376  /* ********************************************************************* */
377 
378  template <typename MATH, typename MATQ, typename Ttol>
379  void block2x2_reduction(MATH &H, MATQ &Q, Ttol tol) {
380  typedef typename linalg_traits<MATH>::value_type T;
381  typedef typename number_traits<T>::magnitude_type R;
382 
383  size_type n = mat_nrows(H), nq = mat_nrows(Q);
384  if (n < 2) return;
385  sub_interval SUBQ(0, nq), SUBL(0, 2);
386  std::vector<T> v(2), w(std::max(n, nq)); v[0] = T(1);
387  tol *= Ttol(2);
388  Ttol tol_i = tol * gmm::abs(H(0,0)), tol_cplx = tol_i;
389  for (size_type i = 0; i < n-1; ++i) {
390  tol_i = (gmm::abs(H(i,i))+gmm::abs(H(i+1,i+1)))*tol;
391  tol_cplx = std::max(tol_cplx, tol_i);
392 
393  if (gmm::abs(H(i+1,i)) > tol_i) { // 2x2 block detected
394  T tr = (H(i+1, i+1) - H(i,i)) / T(2);
395  T delta = tr*tr + H(i,i+1)*H(i+1, i);
396 
397  if (is_complex(T()) || gmm::real(delta) >= R(0)) {
398  sub_interval SUBI(i, 2);
399  T theta = (tr - gmm::sqrt(delta)) / H(i+1,i);
400  R a = gmm::abs(theta);
401  v[1] = (a == R(0)) ? T(-1)
402  : gmm::conj(theta) * (R(1) - gmm::sqrt(a*a + R(1)) / a);
403  row_house_update(sub_matrix(H, SUBI), v, sub_vector(w, SUBL));
404  col_house_update(sub_matrix(H, SUBI), v, sub_vector(w, SUBL));
405  col_house_update(sub_matrix(Q, SUBQ, SUBI), v, sub_vector(w, SUBQ));
406  }
407  ++i;
408  }
409  }
410  }
411 
412  /* ********************************************************************* */
413  /* Basic qr algorithm. */
414  /* ********************************************************************* */
415 
416  #define tol_type_for_qr typename number_traits<typename \
417  linalg_traits<MAT1>::value_type>::magnitude_type
418  #define default_tol_for_qr \
419  (gmm::default_tol(tol_type_for_qr()) * tol_type_for_qr(3))
420 
421  // QR method for real or complex square matrices based on QR factorisation.
422  // eigval has to be a complex vector if A has complex eigeinvalues.
423  // Very slow method. Use implicit_qr_method instead.
424  template <typename MAT1, typename VECT, typename MAT2>
425  void rudimentary_qr_algorithm(const MAT1 &A, const VECT &eigval_,
426  const MAT2 &eigvect_,
427  tol_type_for_qr tol = default_tol_for_qr,
428  bool compvect = true) {
429  VECT &eigval = const_cast<VECT &>(eigval_);
430  MAT2 &eigvect = const_cast<MAT2 &>(eigvect_);
431 
432  typedef typename linalg_traits<MAT1>::value_type value_type;
433 
434  size_type n = mat_nrows(A), p, q = 0, ite = 0;
435  dense_matrix<value_type> Q(n, n), R(n,n), A1(n,n);
436  gmm::copy(A, A1);
437 
438  Hessenberg_reduction(A1, eigvect, compvect);
439  qr_stop_criterion(A1, p, q, tol);
440 
441  while (q < n) {
442  qr_factor(A1, Q, R);
443  gmm::mult(R, Q, A1);
444  if (compvect) { gmm::mult(eigvect, Q, R); gmm::copy(R, eigvect); }
445 
446  qr_stop_criterion(A1, p, q, tol);
447  ++ite;
448  GMM_ASSERT1(ite < n*1000, "QR algorithm failed");
449  }
450  if (compvect) block2x2_reduction(A1, Q, tol);
451  extract_eig(A1, eigval, tol);
452  }
453 
454  template <typename MAT1, typename VECT>
455  void rudimentary_qr_algorithm(const MAT1 &a, VECT &eigval,
456  tol_type_for_qr tol = default_tol_for_qr) {
457  dense_matrix<typename linalg_traits<MAT1>::value_type> m(0,0);
458  rudimentary_qr_algorithm(a, eigval, m, tol, false);
459  }
460 
461  /* ********************************************************************* */
462  /* Francis QR step. */
463  /* ********************************************************************* */
464 
465  template <typename MAT1, typename MAT2>
466  void Francis_qr_step(const MAT1& HH, const MAT2 &QQ, bool compute_Q) {
467  MAT1& H = const_cast<MAT1&>(HH); MAT2& Q = const_cast<MAT2&>(QQ);
468  typedef typename linalg_traits<MAT1>::value_type value_type;
469  size_type n = mat_nrows(H), nq = mat_nrows(Q);
470 
471  std::vector<value_type> v(3), w(std::max(n, nq));
472 
473  value_type s = H(n-2, n-2) + H(n-1, n-1);
474  value_type t = H(n-2, n-2) * H(n-1, n-1) - H(n-2, n-1) * H(n-1, n-2);
475  value_type x = H(0, 0) * H(0, 0) + H(0,1) * H(1, 0) - s * H(0,0) + t;
476  value_type y = H(1, 0) * (H(0,0) + H(1,1) - s);
477  value_type z = H(1, 0) * H(2, 1);
478 
479  sub_interval SUBQ(0, nq);
480 
481  for (size_type k = 0; k < n - 2; ++k) {
482  v[0] = x; v[1] = y; v[2] = z;
483  house_vector(v);
484  size_type r = std::min(k+4, n), q = (k==0) ? 0 : k-1;
485  sub_interval SUBI(k, 3), SUBJ(0, r), SUBK(q, n-q);
486 
487  row_house_update(sub_matrix(H, SUBI, SUBK), v, sub_vector(w, SUBK));
488  col_house_update(sub_matrix(H, SUBJ, SUBI), v, sub_vector(w, SUBJ));
489 
490  if (compute_Q)
491  col_house_update(sub_matrix(Q, SUBQ, SUBI), v, sub_vector(w, SUBQ));
492 
493  x = H(k+1, k); y = H(k+2, k);
494  if (k < n-3) z = H(k+3, k);
495  }
496  sub_interval SUBI(n-2,2), SUBJ(0, n), SUBK(n-3,3), SUBL(0, 3);
497  v.resize(2);
498  v[0] = x; v[1] = y;
499  house_vector(v);
500  row_house_update(sub_matrix(H, SUBI, SUBK), v, sub_vector(w, SUBL));
501  col_house_update(sub_matrix(H, SUBJ, SUBI), v, sub_vector(w, SUBJ));
502  if (compute_Q)
503  col_house_update(sub_matrix(Q, SUBQ, SUBI), v, sub_vector(w, SUBQ));
504  }
505 
506  /* ********************************************************************* */
507  /* Wilkinson Double shift QR step (from Lapack). */
508  /* ********************************************************************* */
509 
510  template <typename MAT1, typename MAT2, typename Ttol>
511  void Wilkinson_double_shift_qr_step(const MAT1& HH, const MAT2 &QQ,
512  Ttol tol, bool exc, bool compute_Q) {
513  MAT1& H = const_cast<MAT1&>(HH); MAT2& Q = const_cast<MAT2&>(QQ);
514  typedef typename linalg_traits<MAT1>::value_type T;
515  typedef typename number_traits<T>::magnitude_type R;
516 
517  size_type n = mat_nrows(H), nq = mat_nrows(Q), m;
518  std::vector<T> v(3), w(std::max(n, nq));
519  const R dat1(0.75), dat2(-0.4375);
520  T h33, h44, h43h34, v1(0), v2(0), v3(0);
521 
522  if (exc) { /* Exceptional shift. */
523  R s = gmm::abs(H(n-1, n-2)) + gmm::abs(H(n-2, n-3));
524  h33 = h44 = dat1 * s;
525  h43h34 = dat2*s*s;
526  }
527  else { /* Wilkinson double shift. */
528  h44 = H(n-1,n-1); h33 = H(n-2, n-2);
529  h43h34 = H(n-1, n-2) * H(n-2, n-1);
530  }
531 
532  /* Look for two consecutive small subdiagonal elements. */
533  /* Determine the effect of starting the double-shift QR iteration at */
534  /* row m, and see if this would make H(m-1, m-2) negligible. */
535  for (m = n-2; m != 0; --m) {
536  T h11 = H(m-1, m-1), h22 = H(m, m);
537  T h21 = H(m, m-1), h12 = H(m-1, m);
538  T h44s = h44 - h11, h33s = h33 - h11;
539  v1 = (h33s*h44s-h43h34) / h21 + h12;
540  v2 = h22 - h11 - h33s - h44s;
541  v3 = H(m+1, m);
542  R s = gmm::abs(v1) + gmm::abs(v2) + gmm::abs(v3);
543  v1 /= s; v2 /= s; v3 /= s;
544  if (m == 1) break;
545  T h00 = H(m-2, m-2);
546  T h10 = H(m-1, m-2);
547  R tst1 = gmm::abs(v1)*(gmm::abs(h00)+gmm::abs(h11)+gmm::abs(h22));
548  if (gmm::abs(h10)*(gmm::abs(v2)+gmm::abs(v3)) <= tol * tst1) break;
549  }
550 
551  /* Double shift QR step. */
552  sub_interval SUBQ(0, nq);
553  for (size_type k = (m == 0) ? 0 : m-1; k < n-2; ++k) {
554  v[0] = v1; v[1] = v2; v[2] = v3;
555  house_vector(v);
556  size_type r = std::min(k+4, n), q = (k==0) ? 0 : k-1;
557  sub_interval SUBI(k, 3), SUBJ(0, r), SUBK(q, n-q);
558 
559  row_house_update(sub_matrix(H, SUBI, SUBK), v, sub_vector(w, SUBK));
560  col_house_update(sub_matrix(H, SUBJ, SUBI), v, sub_vector(w, SUBJ));
561  if (k > m-1) { H(k+1, k-1) = T(0); if (k < n-3) H(k+2, k-1) = T(0); }
562 
563  if (compute_Q)
564  col_house_update(sub_matrix(Q, SUBQ, SUBI), v, sub_vector(w, SUBQ));
565 
566  v1 = H(k+1, k); v2 = H(k+2, k);
567  if (k < n-3) v3 = H(k+3, k);
568  }
569  sub_interval SUBI(n-2,2), SUBJ(0, n), SUBK(n-3,3), SUBL(0, 3);
570  v.resize(2); v[0] = v1; v[1] = v2;
571  house_vector(v);
572  row_house_update(sub_matrix(H, SUBI, SUBK), v, sub_vector(w, SUBL));
573  col_house_update(sub_matrix(H, SUBJ, SUBI), v, sub_vector(w, SUBJ));
574  if (compute_Q)
575  col_house_update(sub_matrix(Q, SUBQ, SUBI), v, sub_vector(w, SUBQ));
576  }
577 
578  /* ********************************************************************* */
579  /* Implicit QR algorithm. */
580  /* ********************************************************************* */
581 
582  // QR method for real or complex square matrices based on an
583  // implicit QR factorisation. eigval has to be a complex vector
584  // if A has complex eigenvalues. Complexity about 10n^3, 25n^3 if
585  // eigenvectors are computed
586  template <typename MAT1, typename VECT, typename MAT2>
587  void implicit_qr_algorithm(const MAT1 &A, const VECT &eigval_,
588  const MAT2 &Q_,
589  tol_type_for_qr tol = default_tol_for_qr,
590  bool compvect = true) {
591  VECT &eigval = const_cast<VECT &>(eigval_);
592  MAT2 &Q = const_cast<MAT2 &>(Q_);
593  typedef typename linalg_traits<MAT1>::value_type T;
594  typedef typename number_traits<T>::magnitude_type R;
595 
596  size_type n(mat_nrows(A)), q(0), p(0);
597  dense_matrix<T> H(n,n);
598  gmm::copy(A, H);
599 
600  Hessenberg_reduction(H, Q, compvect);
601  qr_stop_criterion(H, p, q, tol);
602 
603  size_type ite(0), its(0);
604  sub_interval SUBK(0,0);
605  while (q < n) {
606  sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(Q));
607  if (compvect) SUBK = SUBI;
608 // Francis_qr_step(sub_matrix(H, SUBI),
609 // sub_matrix(Q, SUBJ, SUBK), compvect);
610  Wilkinson_double_shift_qr_step(sub_matrix(H, SUBI),
611  sub_matrix(Q, SUBJ, SUBK),
612  tol, (its == 10 || its == 20), compvect);
613  size_type q_old(q);
614  qr_stop_criterion(H, p, q, tol*R(2));
615  if (q != q_old) its = 0;
616  ++its; ++ite;
617  GMM_ASSERT1(ite < n*100, "QR algorithm failed");
618  }
619  if (compvect) block2x2_reduction(H, Q, tol);
620  extract_eig(H, eigval, tol);
621  }
622 
623  template <typename MAT1, typename VECT>
624  void implicit_qr_algorithm(const MAT1 &a, VECT &eigval,
625  tol_type_for_qr tol = default_tol_for_qr) {
626  dense_matrix<typename linalg_traits<MAT1>::value_type> m(1,1);
627  implicit_qr_algorithm(a, eigval, m, tol, false);
628  }
629 
630  /* ********************************************************************* */
631  /* Implicit symmetric QR step with Wilkinson Shift. */
632  /* ********************************************************************* */
633 
634  template <typename MAT1, typename MAT2>
635  void symmetric_Wilkinson_qr_step(const MAT1& MM, const MAT2 &ZZ,
636  bool compute_z) {
637  MAT1& M = const_cast<MAT1&>(MM); MAT2& Z = const_cast<MAT2&>(ZZ);
638  typedef typename linalg_traits<MAT1>::value_type T;
639  typedef typename number_traits<T>::magnitude_type R;
640  size_type n = mat_nrows(M);
641 
642  for (size_type i = 0; i < n; ++i) {
643  M(i, i) = T(gmm::real(M(i, i)));
644  if (i > 0) {
645  T a = (M(i, i-1) + gmm::conj(M(i-1, i)))/R(2);
646  M(i, i-1) = a; M(i-1, i) = gmm::conj(a);
647  }
648  }
649 
650  R d = gmm::real(M(n-2, n-2) - M(n-1, n-1)) / R(2);
651  R e = gmm::abs_sqr(M(n-1, n-2));
652  R nu = d + gmm::sgn(d)*gmm::sqrt(d*d+e);
653  if (nu == R(0)) { M(n-1, n-2) = T(0); return; }
654  R mu = gmm::real(M(n-1, n-1)) - e / nu;
655  T x = M(0,0) - T(mu), z = M(1, 0), c, s;
656 
657  for (size_type k = 1; k < n; ++k) {
658  Givens_rotation(x, z, c, s);
659 
660  if (k > 1) Apply_Givens_rotation_left(M(k-1,k-2), M(k,k-2), c, s);
661  Apply_Givens_rotation_left(M(k-1,k-1), M(k,k-1), c, s);
662  Apply_Givens_rotation_left(M(k-1,k ), M(k,k ), c, s);
663  if (k < n-1) Apply_Givens_rotation_left(M(k-1,k+1), M(k,k+1), c, s);
664  if (k > 1) Apply_Givens_rotation_right(M(k-2,k-1), M(k-2,k), c, s);
665  Apply_Givens_rotation_right(M(k-1,k-1), M(k-1,k), c, s);
666  Apply_Givens_rotation_right(M(k ,k-1), M(k,k) , c, s);
667  if (k < n-1) Apply_Givens_rotation_right(M(k+1,k-1), M(k+1,k), c, s);
668 
669  if (compute_z) col_rot(Z, c, s, k-1, k);
670  if (k < n-1) { x = M(k, k-1); z = M(k+1, k-1); }
671  }
672  }
673 
674  template <typename VECT1, typename VECT2, typename MAT>
675  void symmetric_Wilkinson_qr_step(const VECT1& diag_, const VECT2& sdiag_,
676  const MAT &ZZ, bool compute_z) {
677  VECT1& diag = const_cast<VECT1&>(diag_);
678  VECT2& sdiag = const_cast<VECT2&>(sdiag_);
679  MAT& Z = const_cast<MAT&>(ZZ);
680  typedef typename linalg_traits<VECT2>::value_type T;
681  typedef typename number_traits<T>::magnitude_type R;
682 
683  size_type n = vect_size(diag);
684  R d = (diag[n-2] - diag[n-1]) / R(2);
685  R e = gmm::abs_sqr(sdiag[n-2]);
686  R nu = d + gmm::sgn(d)*gmm::sqrt(d*d+e);
687  if (nu == R(0)) { sdiag[n-2] = T(0); return; }
688  R mu = diag[n-1] - e / nu;
689  T x = diag[0] - T(mu), z = sdiag[0], c, s;
690 
691  T a01(0), a02(0);
692  T a10(0), a11(diag[0]), a12(gmm::conj(sdiag[0])), a13(0);
693  T a20(0), a21(sdiag[0]), a22(diag[1]), a23(gmm::conj(sdiag[1]));
694  T a31(0), a32(sdiag[1]);
695 
696  for (size_type k = 1; k < n; ++k) {
697  Givens_rotation(x, z, c, s);
698 
699  if (k > 1) Apply_Givens_rotation_left(a10, a20, c, s);
700  Apply_Givens_rotation_left(a11, a21, c, s);
701  Apply_Givens_rotation_left(a12, a22, c, s);
702  if (k < n-1) Apply_Givens_rotation_left(a13, a23, c, s);
703 
704  if (k > 1) Apply_Givens_rotation_right(a01, a02, c, s);
705  Apply_Givens_rotation_right(a11, a12, c, s);
706  Apply_Givens_rotation_right(a21, a22, c, s);
707  if (k < n-1) Apply_Givens_rotation_right(a31, a32, c, s);
708 
709  if (compute_z) col_rot(Z, c, s, k-1, k);
710 
711  diag[k-1] = gmm::real(a11);
712  diag[k] = gmm::real(a22);
713  if (k > 1) sdiag[k-2] = (gmm::conj(a01) + a10) / R(2);
714  sdiag[k-1] = (gmm::conj(a12) + a21) / R(2);
715 
716  x = sdiag[k-1]; z = (gmm::conj(a13) + a31) / R(2);
717 
718  a01 = a12; a02 = a13;
719  a10 = a21; a11 = a22; a12 = a23; a13 = T(0);
720  a20 = a31; a21 = a32; a31 = T(0);
721 
722  if (k < n-1) {
723  sdiag[k] = (gmm::conj(a23) + a32) / R(2);
724  a22 = T(diag[k+1]); a32 = sdiag[k+1]; a23 = gmm::conj(a32);
725  }
726  }
727  }
728 
729  /* ********************************************************************* */
730  /* Implicit QR algorithm for symmetric or hermitian matrices. */
731  /* ********************************************************************* */
732 
733  // implicit QR method for real square symmetric matrices or complex
734  // hermitian matrices.
735  // eigval has to be a complex vector if A has complex eigeinvalues.
736  // complexity about 4n^3/3, 9n^3 if eigenvectors are computed
737  template <typename MAT1, typename VECT, typename MAT2>
738  void symmetric_qr_algorithm_old(const MAT1 &A, const VECT &eigval_,
739  const MAT2 &eigvect_,
740  tol_type_for_qr tol = default_tol_for_qr,
741  bool compvect = true) {
742  VECT &eigval = const_cast<VECT &>(eigval_);
743  MAT2 &eigvect = const_cast<MAT2 &>(eigvect_);
744  typedef typename linalg_traits<MAT1>::value_type T;
745  typedef typename number_traits<T>::magnitude_type R;
746 
747  size_type n(mat_nrows(A)), q(0), p(0), ite(0);
748  dense_matrix<T> Tri(n, n);
749  gmm::copy(A, Tri);
750 
751  if (compvect) gmm::copy(identity_matrix(), eigvect);
752  Householder_tridiagonalization(Tri, eigvect, compvect);
753  symmetric_qr_stop_criterion(Tri, p, q, tol);
754 
755  while (q < n) {
756  sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(eigvect)), SUBK(p, n-p-q);
757  if (!compvect) SUBK = sub_interval(0,0);
758  symmetric_Wilkinson_qr_step(sub_matrix(Tri, SUBI),
759  sub_matrix(eigvect, SUBJ, SUBK), compvect);
760 
761  symmetric_qr_stop_criterion(Tri, p, q, tol*R(2));
762  ++ite;
763  GMM_ASSERT1(ite < n*100, "QR algorithm failed. Probably, your matrix"
764  " is not real symmetric or complex hermitian");
765  }
766 
767  extract_eig(Tri, eigval, tol);
768  }
769 
770  template <typename MAT1, typename VECT, typename MAT2>
771  void symmetric_qr_algorithm(const MAT1 &A, const VECT &eigval_,
772  const MAT2 &eigvect_,
773  tol_type_for_qr tol = default_tol_for_qr,
774  bool compvect = true) {
775  VECT &eigval = const_cast<VECT &>(eigval_);
776  MAT2 &eigvect = const_cast<MAT2 &>(eigvect_);
777  typedef typename linalg_traits<MAT1>::value_type T;
778  typedef typename number_traits<T>::magnitude_type R;
779 
780  size_type n = mat_nrows(A), q = 0, p, ite = 0;
781  if (compvect) gmm::copy(identity_matrix(), eigvect);
782  if (n == 0) return;
783  if (n == 1) { eigval[0]=gmm::real(A(0,0)); return; }
784  dense_matrix<T> Tri(n, n);
785  gmm::copy(A, Tri);
786 
787  Householder_tridiagonalization(Tri, eigvect, compvect);
788 
789  std::vector<R> diag(n);
790  std::vector<T> sdiag(n);
791  for (size_type i = 0; i < n; ++i)
792  { diag[i] = gmm::real(Tri(i, i)); if (i+1 < n) sdiag[i] = Tri(i+1, i); }
793 
794  symmetric_qr_stop_criterion(diag, sdiag, p, q, tol);
795 
796  while (q < n) {
797  sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(eigvect)), SUBK(p, n-p-q);
798  if (!compvect) SUBK = sub_interval(0,0);
799 
800  symmetric_Wilkinson_qr_step(sub_vector(diag, SUBI),
801  sub_vector(sdiag, SUBI),
802  sub_matrix(eigvect, SUBJ, SUBK), compvect);
803 
804  symmetric_qr_stop_criterion(diag, sdiag, p, q, tol*R(3));
805  ++ite;
806  GMM_ASSERT1(ite < n*100, "QR algorithm failed.");
807  }
808 
809  gmm::copy(diag, eigval);
810  }
811 
812 
813  template <typename MAT1, typename VECT>
814  void symmetric_qr_algorithm(const MAT1 &a, VECT &eigval,
815  tol_type_for_qr tol = default_tol_for_qr) {
816  dense_matrix<typename linalg_traits<MAT1>::value_type> m(0,0);
817  symmetric_qr_algorithm(a, eigval, m, tol, false);
818  }
819 
820 
821 }
822 
823 #endif
824 
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
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::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
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_dense_Householder.h
Householder for dense matrices.
gmm::qr_factor
void qr_factor(const MAT &A_)
QR factorization using Householder method (complex and real version).
Definition: gmm_dense_qr.h:48
gmm::extract_eig
void extract_eig(const MAT &A, const VECT &V, Ttol tol)
Compute eigenvalue vector.
Definition: gmm_dense_qr.h:302
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48