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