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