GetFEM  5.5
gmm_solver_idgmres.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2026 Yves Renard, Caroline Lecalvez
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /**@file gmm_solver_idgmres.h
32  @author Caroline Lecalvez <Caroline.Lecalvez@gmm.insa-tlse.fr>
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date October 6, 2003.
35  @brief Implicitly restarted and deflated Generalized Minimum Residual.
36 */
37 #ifndef GMM_IDGMRES_H
38 #define GMM_IDGMRES_H
39 
40 #include "gmm_kernel.h"
41 #include "gmm_iter.h"
42 #include "gmm_dense_sylvester.h"
43 
44 namespace gmm {
45 
46  template <typename T> compare_vp {
47  bool operator()(const std::pair<T, size_type> &a,
48  const std::pair<T, size_type> &b) const
49  { return (gmm::abs(a.first) > gmm::abs(b.first)); }
50  }
51 
52  struct idgmres_state {
53  size_type m, tb_deb, tb_def, p, k, nb_want, nb_unwant;
54  size_type nb_nolong, tb_deftot, tb_defwant, conv, nb_un, fin;
55  bool ok;
56 
57  idgmres_state(size_type mm, size_type pp, size_type kk)
58  : m(mm), tb_deb(1), tb_def(0), p(pp), k(kk), nb_want(0),
59  nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
60  conv(0), nb_un(0), fin(0), ok(false); {}
61  }
62 
63  idgmres_state(size_type mm, size_type pp, size_type kk)
64  : m(mm), tb_deb(1), tb_def(0), p(pp), k(kk), nb_want(0),
65  nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
66  conv(0), nb_un(0), fin(0), ok(false); {}
67 
68 
69  template <typename CONT, typename IND>
70  apply_permutation(CONT &cont, const IND &ind) {
71  size_type m = ind.end() - ind.begin();
72  std::vector<bool> sorted(m, false);
73 
74  for (size_type l = 0; l < m; ++l)
75  if (!sorted[l] && ind[l] != l) {
76 
77  typeid(cont[0]) aux = cont[l];
78  k = ind[l];
79  cont[l] = cont[k];
80  sorted[l] = true;
81 
82  for(k2 = ind[k]; k2 != l; k2 = ind[k]) {
83  cont[k] = cont[k2];
84  sorted[k] = true;
85  k = k2;
86  }
87  cont[k] = aux;
88  }
89  }
90 
91 
92  /** Implicitly restarted and deflated Generalized Minimum Residual
93 
94  See: C. Le Calvez, B. Molina, Implicitly restarted and deflated
95  FOM and GMRES, numerical applied mathematics,
96  (30) 2-3 (1999) pp191-212.
97 
98  @param A Real or complex unsymmetric matrix.
99  @param x initial guess vector and final result.
100  @param b right hand side
101  @param M preconditionner
102  @param m size of the subspace between two restarts
103  @param p number of converged ritz values seeked
104  @param k size of the remaining Krylov subspace when the p ritz values
105  have not yet converged 0 <= p <= k < m.
106  @param tol_vp : tolerance on the ritz values.
107  @param outer
108  @param KS
109  */
110  template < typename Mat, typename Vec, typename VecB, typename Precond,
111  typename Basis >
112  void idgmres(const Mat &A, Vec &x, const VecB &b, const Precond &M,
113  size_type m, size_type p, size_type k, double tol_vp,
114  iteration &outer, Basis& KS) {
115 
116  typedef typename linalg_traits<Mat>::value_type T;
117  typedef typename number_traits<T>::magnitude_type R;
118 
119  R a, beta;
120  idgmres_state st(m, p, k);
121 
122  std::vector<T> w(vect_size(x)), r(vect_size(x)), u(vect_size(x));
123  std::vector<T> c_rot(m+1), s_rot(m+1), s(m+1);
124  std::vector<T> y(m+1), ztest(m+1), gam(m+1);
125  std::vector<T> gamma(m+1);
126  gmm::dense_matrix<T> H(m+1, m), Hess(m+1, m),
127  Hobl(m+1, m), W(vect_size(x), m+1);
128  gmm::clear(H);
129 
130  outer.set_rhsnorm(gmm::vect_norm2(b));
131  if (outer.get_rhsnorm() == 0.0) { clear(x); return; }
132 
133  mult(A, scaled(x, -1.0), b, w);
134  mult(M, w, r);
135  beta = gmm::vect_norm2(r);
136 
137  iteration inner = outer;
138  inner.reduce_noisy();
139  inner.set_maxiter(m);
140  inner.set_name("GMRes inner iter");
141 
142  while (! outer.finished(beta)) {
143 
144  gmm::copy(gmm::scaled(r, 1.0/beta), KS[0]);
145  gmm::clear(s);
146  s[0] = beta;
147  gmm::copy(s, gamma);
148 
149  inner.set_maxiter(m - st.tb_deb + 1);
150  size_type i = st.tb_deb - 1; inner.init();
151 
152  do {
153  mult(A, KS[i], u);
154  mult(M, u, KS[i+1]);
155  orthogonalize_with_refinment(KS, mat_col(H, i), i);
156  H(i+1, i) = a = gmm::vect_norm2(KS[i+1]);
157  gmm::scale(KS[i+1], R(1) / a);
158 
159  gmm::copy(mat_col(H, i), mat_col(Hess, i));
160  gmm::copy(mat_col(H, i), mat_col(Hobl, i));
161 
162  for (size_type l = 0; l < i; ++l)
163  Apply_Givens_rotation_left(H(l,i), H(l+1,i), c_rot[l], s_rot[l]);
164 
165  Givens_rotation(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
166  Apply_Givens_rotation_left(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
167  H(i+1, i) = T(0);
168  Apply_Givens_rotation_left(s[i], s[i+1], c_rot[i], s_rot[i]);
169 
170  ++inner, ++outer, ++i;
171  } while (! inner.finished(gmm::abs(s[i])));
172 
173  if (inner.converged()) {
174  gmm::copy(s, y);
175  upper_tri_solve(H, y, i, false);
176  combine(KS, y, x, i);
177  mult(A, gmm::scaled(x, T(-1)), b, w);
178  mult(M, w, r);
179  beta = gmm::vect_norm2(r); // + verif sur beta ... à faire
180  break;
181  }
182 
183  gmm::clear(gam); gam[m] = s[i];
184  for (size_type l = m; l > 0; --l)
185  Apply_Givens_rotation_left(gam[l-1], gam[l], gmm::conj(c_rot[l-1]),
186  -s_rot[l-1]);
187 
188  mult(KS.mat(), gam, r);
189  beta = gmm::vect_norm2(r);
190 
191  mult(Hess, scaled(y, T(-1)), gamma, ztest);
192  // En fait, d'après Caroline qui s'y connait ztest et gam devrait
193  // être confondus
194  // Quand on aura vérifié que ça marche, il faudra utiliser gam à la
195  // place de ztest.
196  if (st.tb_def < p) {
197  T nss = H(m,m-1) / ztest[m];
198  nss /= gmm::abs(nss); // ns à calculer plus tard aussi
199  gmm::copy(KS.mat(), W); gmm::copy(scaled(r, nss /beta), mat_col(W, m));
200 
201  // Computation of the oblique matrix
202  sub_interval SUBI(0, m);
203  add(scaled(sub_vector(ztest, SUBI), -Hobl(m, m-1) / ztest[m]),
204  sub_vector(mat_col(Hobl, m-1), SUBI));
205  Hobl(m, m-1) *= nss * beta / ztest[m];
206 
207  /* **************************************************************** */
208  /* Locking */
209  /* **************************************************************** */
210 
211  // Computation of the Ritz eigenpairs.
212  std::vector<std::complex<R> > eval(m);
213  dense_matrix<T> YB(m-st.tb_def, m-st.tb_def);
214  std::vector<char> pure(m-st.tb_def, 0);
215  gmm::clear(YB);
216 
217  select_eval(Hobl, eval, YB, pure, st);
218 
219  if (st.conv != 0) {
220  // DEFLATION using the QR Factorization of YB
221 
222  T alpha = Lock(W, Hobl,
223  sub_matrix(YB, sub_interval(0, m-st.tb_def)),
224  sub_interval(st.tb_def, m-st.tb_def),
225  (st.tb_defwant < p));
226  // ns *= alpha; // à calculer plus tard ??
227  // V(:,m+1) = alpha*V(:, m+1); ça devait servir à qlq chose ...
228 
229 
230  // Clean the portions below the diagonal corresponding
231  // to the lock Schur vectors
232  for (size_type j = st.tb_def; j < st.tb_deftot; ++j) {
233  if ( pure[j-st.tb_def] == 0)
234  gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
235  else if (pure[j-st.tb_def] == 1) {
236  gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
237  sub_interval(j, 2)));
238  ++j;
239  }
240  else GMM_ASSERT3(false, "internal error");
241  }
242 
243  if (!st.ok) {
244  // attention si m = 0;
245  size_type mm = std::min(k+st.nb_unwant+st.nb_nolong, m-1);
246 
247  if (eval_sort[m-mm-1].second != R(0)
248  && eval_sort[m-mm-1].second == -eval_sort[m-mm].second)
249  ++mm;
250 
251  std::vector<complex<R> > shifts(m-mm);
252  for (size_type i = 0; i < m-mm; ++i)
253  shifts[i] = eval_sort[i].second;
254 
255  apply_shift_to_Arnoldi_factorization(W, Hobl, shifts, mm,
256  m-mm, true);
257 
258  st.fin = mm;
259  }
260  else
261  st.fin = st.tb_deftot;
262 
263 
264  /* ************************************************************** */
265  /* Purge */
266  /* ************************************************************** */
267 
268  if (st.nb_nolong + st.nb_unwant > 0) {
269 
270  std::vector<std::complex<R> > eval(m);
271  dense_matrix<T> YB(st.fin, st.tb_deftot);
272  std::vector<char> pure(st.tb_deftot, 0);
273  gmm::clear(YB);
274  st.nb_un = st.nb_nolong + st.nb_unwant;
275 
276  select_eval_for_purging(Hobl, eval, YB, pure, st);
277 
278  T alpha = Lock(W, Hobl, YB, sub_interval(0, st.fin), ok);
279 
280  // Clean the portions below the diagonal corresponding
281  // to the unwanted lock Schur vectors
282  for (size_type j = 0; j < st.tb_deftot; ++j) {
283  if ( pure[j] == 0)
284  gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
285  else if (pure[j] == 1) {
286  gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
287  sub_interval(j, 2)));
288  ++j;
289  }
290  else GMM_ASSERT3(false, "internal error");
291  }
292 
293  gmm::dense_matrix<T> z(st.nb_un, st.fin - st.nb_un);
294  sub_interval SUBI(0, st.nb_un), SUBJ(st.nb_un, st.fin - st.nb_un);
295  sylvester(sub_matrix(Hobl, SUBI),
296  sub_matrix(Hobl, SUBJ),
297  sub_matrix(gmm::scaled(Hobl, -T(1)), SUBI, SUBJ), z);
298  }
299  }
300  }
301  }
302  }
303 
304 
305  template < typename Mat, typename Vec, typename VecB, typename Precond >
306  void idgmres(const Mat &A, Vec &x, const VecB &b,
307  const Precond &M, size_type m, iteration& outer) {
308  typedef typename linalg_traits<Mat>::value_type T;
309  modified_gram_schmidt<T> orth(m, vect_size(x));
310  gmres(A, x, b, M, m, outer, orth);
311  }
312 
313 
314  // Lock stage of an implicit restarted Arnoldi process.
315  // 1- QR factorization of YB through Householder matrices
316  // Q(Rl) = YB
317  // (0 )
318  // 2- Update of the Arnoldi factorization.
319  // H <- Q*HQ, W <- WQ
320  // 3- Restore the Hessemberg form of H.
321 
322  template <typename T, typename MATYB>
323  void Lock(gmm::dense_matrix<T> &W, gmm::dense_matrix<T> &H,
324  const MATYB &YB, const sub_interval SUB,
325  bool restore, T &ns) {
326 
327  size_type n = mat_nrows(W), m = mat_ncols(W) - 1;
328  size_type ncols = mat_ncols(YB), nrows = mat_nrows(YB);
329  size_type begin = min(SUB); end = max(SUB) - 1;
330  sub_interval SUBR(0, nrows), SUBC(0, ncols);
331  T alpha(1);
332 
333  GMM_ASSERT2(((end-begin) == ncols) && (m == mat_nrows(H))
334  && (m+1 == mat_ncols(H)), "dimensions mismatch");
335 
336  // DEFLATION using the QR Factorization of YB
337 
338  dense_matrix<T> QR(n_rows, n_rows);
339  gmmm::copy(YB, sub_matrix(QR, SUBR, SUBC));
340  gmm::clear(submatrix(QR, SUBR, sub_interval(ncols, nrows-ncols)));
341  qr_factor(QR);
342 
343  apply_house_left(QR, sub_matrix(H, SUB));
344  apply_house_right(QR, sub_matrix(H, SUBR, SUB));
345  apply_house_right(QR, sub_matrix(W, sub_interval(0, n), SUB));
346 
347  // Restore to the initial block hessenberg form
348 
349  if (restore) {
350 
351  // verifier quand m = 0 ...
352  gmm::dense_matrix tab_p(end - st.tb_deftot, end - st.tb_deftot);
353  gmm::copy(identity_matrix(), tab_p);
354 
355  for (size_type j = end-1; j >= st.tb_deftot+2; --j) {
356 
357  size_type jm = j-1;
358  std::vector<T> v(jm - st.tb_deftot);
359  sub_interval SUBtot(st.tb_deftot, jm - st.tb_deftot);
360  sub_interval SUBtot2(st.tb_deftot, end - st.tb_deftot);
361  gmm::copy(sub_vector(mat_row(H, j), SUBtot), v);
362  house_vector_last(v);
363  w.resize(end);
364  col_house_update(sub_matrix(H, SUBI, SUBtot), v, w);
365  w.resize(end - st.tb_deftot);
366  row_house_update(sub_matrix(H, SUBtot, SUBtot2), v, w);
367  gmm::clear(sub_vector(mat_row(H, j),
368  sub_interval(st.tb_deftot, j-1-st.tb_deftot)));
369  w.resize(end - st.tb_deftot);
370  col_house_update(sub_matrix(tab_p, sub_interval(0, end-st.tb_deftot),
371  sub_interval(0, jm-st.tb_deftot)), v, w);
372  w.resize(n);
373  col_house_update(sub_matrix(W, sub_interval(0, n), SUBtot), v, w);
374  }
375 
376  // restore positive subdiagonal elements
377 
378  std::vector<T> d(fin-st.tb_deftot); d[0] = T(1);
379 
380  // We compute d[i+1] in order
381  // (d[i+1] * H(st.tb_deftot+i+1,st.tb_deftoti)) / d[i]
382  // be equal to |H(st.tb_deftot+i+1,st.tb_deftot+i))|.
383  for (size_type j = 0; j+1 < end-st.tb_deftot; ++j) {
384  T e = H(st.tb_deftot+j, st.tb_deftot+j-1);
385  d[j+1] = (e == T(0)) ? T(1) : d[j] * gmm::abs(e) / e;
386  scale(sub_vector(mat_row(H, st.tb_deftot+j+1),
387  sub_interval(st.tb_deftot, m-st.tb_deftot)), d[j+1]);
388  scale(mat_col(H, st.tb_deftot+j+1), T(1) / d[j+1]);
389  scale(mat_col(W, st.tb_deftot+j+1), T(1) / d[j+1]);
390  }
391 
392  alpha = tab_p(end-st.tb_deftot-1, end-st.tb_deftot-1) / d[end-st.tb_deftot-1];
393  alpha /= gmm::abs(alpha);
394  scale(mat_col(W, m), alpha);
395  }
396 
397  return alpha;
398  }
399 
400 
401  // Apply p implicit shifts to the Arnoldi factorization
402  // AV = VH+H(k+p+1,k+p) V(:,k+p+1) e_{k+p}*
403  // and produces the following new Arnoldi factorization
404  // A(VQ) = (VQ)(Q*HQ)+H(k+p+1,k+p) V(:,k+p+1) e_{k+p}* Q
405  // where only the first k columns are relevant.
406  //
407  // Dan Sorensen and Richard J. Radke, 11/95
408  template<typename T, typename C>
409  apply_shift_to_Arnoldi_factorization(dense_matrix<T> V, dense_matrix<T> H,
410  std::vector<C> Lambda, size_type &k,
411  size_type p, bool true_shift = false) {
412 
413  size_type k1 = 0, num = 0, kend = k+p, kp1 = k + 1;
414  bool mark = false;
415  T c, s, x, y, z;
416 
417  dense_matrix<T> q(1, kend);
418  gmm::clear(q); q(0,kend-1) = T(1);
419  std::vector<T> hv(3), w(std::max(kend, mat_nrows(V)));
420 
421  for(size_type jj = 0; jj < p; ++jj) {
422  // compute and apply a bulge chase sweep initiated by the
423  // implicit shift held in w(jj)
424 
425  if (abs(Lambda[jj].real()) == 0.0) {
426  // apply a real shift using 2 by 2 Givens rotations
427 
428  for (size_type k1 = 0, k2 = 0; k2 != kend-1; k1 = k2+1) {
429  k2 = k1;
430  while (h(k2+1, k2) != T(0) && k2 < kend-1)
431  ++k2;
432 
433  Givens_rotation(H(k1, k1) - Lambda[jj], H(k1+1, k1), c, s);
434 
435  for (i = k1; i <= k2; ++i) {
436  if (i > k1)
437  Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
438 
439  // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
440  // Vérifier qu'au final H(i+1,i) est bien un réel positif.
441 
442  // apply rotation from left to rows of H
443  row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
444  c, s, 0, 0);
445 
446  // apply rotation from right to columns of H
447  size_type ip2 = std::min(i+2, kend);
448  col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
449  c, s, 0, 0);
450 
451  // apply rotation from right to columns of V
452  col_rot(V, c, s, i, i+1);
453 
454  // accumulate e' Q so residual can be updated k+p
455  Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
456  // peut être que nous utilisons G au lieu de G* et que
457  // nous allons trop loin en k2.
458  }
459  }
460 
461  num = num + 1;
462  }
463  else {
464 
465  // Apply a double complex shift using 3 by 3 Householder
466  // transformations
467 
468  if (jj == p || mark)
469  mark = false; // skip application of conjugate shift
470  else {
471  num = num + 2; // mark that a complex conjugate
472  mark = true; // pair has been applied
473 
474  // Indices de fin de boucle à surveiller... de près !
475  for (size_type k1 = 0, k3 = 0; k3 != kend-2; k1 = k3+1) {
476  k3 = k1;
477  while (h(k3+1, k3) != T(0) && k3 < kend-2)
478  ++k3;
479  size_type k2 = k1+1;
480 
481 
482  x = H(k1,k1) * H(k1,k1) + H(k1,k2) * H(k2,k1)
483  - 2.0*Lambda[jj].real() * H(k1,k1) + gmm::abs_sqr(Lambda[jj]);
484  y = H(k2,k1) * (H(k1,k1) + H(k2,k2) - 2.0*Lambda[jj].real());
485  z = H(k2+1,k2) * H(k2,k1);
486 
487  for (size_type i = k1; i <= k3; ++i) {
488  if (i > k1) {
489  x = H(i, i-1);
490  y = H(i+1, i-1);
491  z = H(i+2, i-1);
492  // Ne pas oublier de nettoyer H(i+1,i-1) et H(i+2,i-1)
493  // (les mettre à zéro).
494  }
495 
496  hv[0] = x; hv[1] = y; hv[2] = z;
497  house_vector(v);
498 
499  // Vérifier qu'au final H(i+1,i) est bien un réel positif
500 
501  // apply transformation from left to rows of H
502  w.resize(kend-i);
503  row_house_update(sub_matrix(H, sub_interval(i, 2),
504  sub_interval(i, kend-i)),
505  v, w);
506 
507  // apply transformation from right to columns of H
508 
509  size_type ip3 = std::min(kend, i + 3);
510  w.resize(ip3);
511  col_house_update(sub_matrix(H, sub_interval(0, ip3),
512  sub_interval(i, 2)),
513  v, w);
514 
515  // apply transformation from right to columns of V
516 
517  w.resize(mat_nrows(V));
518  col_house_update(sub_matrix(V, sub_interval(0, mat_nrows(V)),
519  sub_interval(i, 2)),
520  v, w);
521 
522  // accumulate e' Q so residual can be updated k+p
523  w.resize(1);
524  col_house_update(sub_matrix(q, sub_interval(0,1),
525  sub_interval(i,2)),
526  v, w);
527  }
528  }
529 
530  // clean up step with Givens rotation
531 
532  i = kend-2;
533  c = x;
534  s = y;
535  if (i > k1)
536  Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
537 
538  // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
539  // Vérifier qu'au final H(i+1,i) est bien un réel positif.
540 
541  // apply rotation from left to rows of H
542  row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
543  c, s, 0, 0);
544 
545  // apply rotation from right to columns of H
546  size_type ip2 = std::min(i+2, kend);
547  col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
548  c, s, 0, 0);
549 
550  // apply rotation from right to columns of V
551  col_rot(V, c, s, i, i+1);
552 
553  // accumulate e' Q so residual can be updated k+p
554  Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
555  }
556  }
557  }
558 
559  // update residual and store in the k+1 -st column of v
560 
561  k = kend - num;
562  scale(mat_col(V, kend), q(0, k));
563 
564  if (k < mat_nrows(H)) {
565  if (true_shift)
566  gmm::copy(mat_col(V, kend), mat_col(V, k));
567  else
568  // v(:,k+1) = v(:,kend+1) + v(:,k+1)*h(k+1,k);
569  // v(:,k+1) = v(:,kend+1) ;
570  gmm::add(scaled(mat_col(V, kend), H(kend, kend-1)),
571  scaled(mat_col(V, k), H(k, k-1)), mat_col(V, k));
572  }
573 
574  H(k, k-1) = vect_norm2(mat_col(V, k));
575  scale(mat_col(V, kend), T(1) / H(k, k-1));
576  }
577 
578 
579 
580  template<typename MAT, typename EVAL, typename PURE>
581  void select_eval(const MAT &Hobl, EVAL &eval, MAT &YB, PURE &pure,
582  idgmres_state &st) {
583 
584  typedef typename linalg_traits<MAT>::value_type T;
585  typedef typename number_traits<T>::magnitude_type R;
586  size_type m = st.m;
587 
588  // Computation of the Ritz eigenpairs.
589 
590  col_matrix< std::vector<T> > evect(m-st.tb_def, m-st.tb_def);
591  // std::vector<std::complex<R> > eval(m);
592  std::vector<R> ritznew(m, T(-1));
593 
594  // dense_matrix<T> evect_lock(st.tb_def, st.tb_def);
595 
596  sub_interval SUB1(st.tb_def, m-st.tb_def);
597  implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
598  sub_vector(eval, SUB1), evect);
599  sub_interval SUB2(0, st.tb_def);
600  implicit_qr_algorithm(sub_matrix(Hobl, SUB2),
601  sub_vector(eval, SUB2), /* evect_lock */);
602 
603  for (size_type l = st.tb_def; l < m; ++l)
604  ritznew[l] = gmm::abs(evect(m-st.tb_def-1, l-st.tb_def) * Hobl(m, m-1));
605 
606  std::vector< std::pair<T, size_type> > eval_sort(m);
607  for (size_type l = 0; l < m; ++l)
608  eval_sort[l] = std::pair<T, size_type>(eval[l], l);
609  std::sort(eval_sort.begin(), eval_sort.end(), compare_vp());
610 
611  std::vector<size_type> index(m);
612  for (size_type l = 0; l < m; ++l) index[l] = eval_sort[l].second;
613 
614  std::vector<bool> kept(m, false);
615  std::fill(kept.begin(), kept.begin()+st.tb_def, true);
616 
617  apply_permutation(eval, index);
618  apply_permutation(evect, index);
619  apply_permutation(ritznew, index);
620  apply_permutation(kept, index);
621 
622  // Which are the eigenvalues that converged ?
623  //
624  // nb_want is the number of eigenvalues of
625  // Hess(tb_def+1:n,tb_def+1:n) that converged and are WANTED
626  //
627  // nb_unwant is the number of eigenvalues of
628  // Hess(tb_def+1:n,tb_def+1:n) that converged and are UNWANTED
629  //
630  // nb_nolong is the number of eigenvalues of
631  // Hess(1:tb_def,1:tb_def) that are NO LONGER WANTED.
632  //
633  // tb_deftot is the number of the deflated eigenvalues
634  // that is tb_def + nb_want + nb_unwant
635  //
636  // tb_defwant is the number of the wanted deflated eigenvalues
637  // that is tb_def + nb_want - nb_nolong
638 
639  st.nb_want = 0, st.nb_unwant = 0, st.nb_nolong = 0;
640  size_type j, ind;
641 
642  for (j = 0, ind = 0; j < m-p; ++j) {
643  if (ritznew[j] == R(-1)) {
644  if (std::imag(eval[j]) != R(0)) {
645  st.nb_nolong += 2; ++j; // à adapter dans le cas complexe ...
646  } else
647  st.nb_nolong++;
648  } else {
649  if (ritznew[j] < tol_vp * gmm::abs(eval[j])) {
650 
651  for (size_type l = 0, l < m-st.tb_def; ++l)
652  YB(l, ind) = std::real(evect(l, j));
653  kept[j] = true;
654  ++j; ++st.nb_unwant; ind++;
655 
656  if (std::imag(eval[j]) != R(0)) {
657  for (size_type l = 0, l < m-st.tb_def; ++l)
658  YB(l, ind) = std::imag(evect(l, j));
659  pure[ind-1] = 1;
660  pure[ind] = 2;
661 
662  kept[j] = true;
663 
664  st.nb_unwant++;
665  ++ind;
666  }
667  }
668  }
669  }
670 
671 
672  for (; j < m; ++j) {
673  if (ritznew[j] != R(-1)) {
674 
675  for (size_type l = 0, l < m-st.tb_def; ++l)
676  YB(l, ind) = std::real(evect(l, j));
677  pure[ind] = 1;
678  ++ind;
679  kept[j] = true;
680  ++st.nb_want;
681 
682  if (ritznew[j] < tol_vp * gmm::abs(eval[j])) {
683  for (size_type l = 0, l < m-st.tb_def; ++l)
684  YB(l, ind) = std::imag(evect(l, j));
685  pure[ind] = 2;
686 
687  j++;
688  kept[j] = true;
689 
690  st.nb_want++;
691  ++ind;
692  }
693  }
694  }
695 
696  std::vector<T> shift(m - st.tb_def - st.nb_want - st.nb_unwant);
697  for (size_type j = 0, i = 0; j < m; ++j)
698  if (!kept[j])
699  shift[i++] = eval[j];
700 
701  // st.conv (st.nb_want+st.nb_unwant) is the number of eigenpairs that
702  // have just converged.
703  // st.tb_deftot is the total number of eigenpairs that have converged.
704 
705  size_type st.conv = ind;
706  size_type st.tb_deftot = st.tb_def + st.conv;
707  size_type st.tb_defwant = st.tb_def + st.nb_want - st.nb_nolong;
708 
709  sub_interval SUBYB(0, st.conv);
710 
711  if ( st.tb_defwant >= p ) { // An invariant subspace has been found.
712 
713  st.nb_unwant = 0;
714  st.nb_want = p + st.nb_nolong - st.tb_def;
715  st.tb_defwant = p;
716 
717  if ( pure[st.conv - st.nb_want + 1] == 2 ) {
718  ++st.nb_want; st.tb_defwant = ++p;// il faudrait que ce soit un p local
719  }
720 
721  SUBYB = sub_interval(st.conv - st.nb_want, st.nb_want);
722  // YB = YB(:, st.conv-st.nb_want+1 : st.conv); // On laisse en suspend ..
723  // pure = pure(st.conv-st.nb_want+1 : st.conv,1); // On laisse suspend ..
724  st.conv = st.nb_want;
725  st.tb_deftot = st.tb_def + st.conv;
726  st.ok = true;
727  }
728  }
729 
730 
731 
732  template<typename MAT, typename EVAL, typename PURE>
733  void select_eval_for_purging(const MAT &Hobl, EVAL &eval, MAT &YB,
734  PURE &pure, idgmres_state &st) {
735 
736  typedef typename linalg_traits<MAT>::value_type T;
737  typedef typename number_traits<T>::magnitude_type R;
738  size_type m = st.m;
739 
740  // Computation of the Ritz eigenpairs.
741 
742  col_matrix< std::vector<T> > evect(st.tb_deftot, st.tb_deftot);
743 
744  sub_interval SUB1(0, st.tb_deftot);
745  implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
746  sub_vector(eval, SUB1), evect);
747  std::fill(eval.begin() + st.tb_deftot, eval.end(), std::complex<R>(0));
748 
749  std::vector< std::pair<T, size_type> > eval_sort(m);
750  for (size_type l = 0; l < m; ++l)
751  eval_sort[l] = std::pair<T, size_type>(eval[l], l);
752  std::sort(eval_sort.begin(), eval_sort.end(), compare_vp());
753 
754  std::vector<bool> sorted(m);
755  std::fill(sorted.begin(), sorted.end(), false);
756 
757  std::vector<size_type> ind(m);
758  for (size_type l = 0; l < m; ++l) ind[l] = eval_sort[l].second;
759 
760  std::vector<bool> kept(m, false);
761  std::fill(kept.begin(), kept.begin()+st.tb_def, true);
762 
763  apply_permutation(eval, ind);
764  apply_permutation(evect, ind);
765 
766  size_type j;
767  for (j = 0; j < st.tb_deftot; ++j) {
768 
769  for (size_type l = 0, l < st.tb_deftot; ++l)
770  YB(l, j) = std::real(evect(l, j));
771 
772  if (std::imag(eval[j]) != R(0)) {
773  for (size_type l = 0, l < m-st.tb_def; ++l)
774  YB(l, j+1) = std::imag(evect(l, j));
775  pure[j] = 1;
776  pure[j+1] = 2;
777 
778  j += 2;
779  }
780  else ++j;
781  }
782  }
783 
784 
785 }
786 
787 #endif
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:556
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
gmm::qr_factor
void qr_factor(const MAT &A_)
QR factorization using Householder method (complex and real version).
Definition: gmm_dense_qr.h:48
gmm_dense_sylvester.h
Sylvester equation solver.
gmm_iter.h
Iteration object.
gmm_kernel.h
Include the base gmm files.
gmm::gmres
void gmres(const Mat &A, Vec &x, const VecB &b, const Precond &M, int restart, iteration &outer, Basis &KS)
Generalized Minimum Residual.
Definition: gmm_solver_gmres.h:89
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:46