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