GetFEM  5.4.2
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 
130  gmm::clear(H);
131 
132  outer.set_rhsnorm(gmm::vect_norm2(b));
133  if (outer.get_rhsnorm() == 0.0) { clear(x); return; }
134 
135  mult(A, scaled(x, -1.0), b, w);
136  mult(M, w, r);
137  beta = gmm::vect_norm2(r);
138 
139  iteration inner = outer;
140  inner.reduce_noisy();
141  inner.set_maxiter(m);
142  inner.set_name("GMRes inner iter");
143 
144  while (! outer.finished(beta)) {
145 
146  gmm::copy(gmm::scaled(r, 1.0/beta), KS[0]);
147  gmm::clear(s);
148  s[0] = beta;
149  gmm::copy(s, gamma);
150 
151  inner.set_maxiter(m - st.tb_deb + 1);
152  size_type i = st.tb_deb - 1; inner.init();
153 
154  do {
155  mult(A, KS[i], u);
156  mult(M, u, KS[i+1]);
157  orthogonalize_with_refinment(KS, mat_col(H, i), i);
158  H(i+1, i) = a = gmm::vect_norm2(KS[i+1]);
159  gmm::scale(KS[i+1], R(1) / a);
160 
161  gmm::copy(mat_col(H, i), mat_col(Hess, i));
162  gmm::copy(mat_col(H, i), mat_col(Hobl, i));
163 
164 
165  for (size_type l = 0; l < i; ++l)
166  Apply_Givens_rotation_left(H(l,i), H(l+1,i), c_rot[l], s_rot[l]);
167 
168  Givens_rotation(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
169  Apply_Givens_rotation_left(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
170  H(i+1, i) = T(0);
171  Apply_Givens_rotation_left(s[i], s[i+1], c_rot[i], s_rot[i]);
172 
173  ++inner, ++outer, ++i;
174  } while (! inner.finished(gmm::abs(s[i])));
175 
176  if (inner.converged()) {
177  gmm::copy(s, y);
178  upper_tri_solve(H, y, i, false);
179  combine(KS, y, x, i);
180  mult(A, gmm::scaled(x, T(-1)), b, w);
181  mult(M, w, r);
182  beta = gmm::vect_norm2(r); // + verif sur beta ... à faire
183  break;
184  }
185 
186  gmm::clear(gam); gam[m] = s[i];
187  for (size_type l = m; l > 0; --l)
188  Apply_Givens_rotation_left(gam[l-1], gam[l], gmm::conj(c_rot[l-1]),
189  -s_rot[l-1]);
190 
191  mult(KS.mat(), gam, r);
192  beta = gmm::vect_norm2(r);
193 
194  mult(Hess, scaled(y, T(-1)), gamma, ztest);
195  // En fait, d'après Caroline qui s'y connait ztest et gam devrait
196  // être confondus
197  // Quand on aura vérifié que ça marche, il faudra utiliser gam à la
198  // place de ztest.
199  if (st.tb_def < p) {
200  T nss = H(m,m-1) / ztest[m];
201  nss /= gmm::abs(nss); // ns à calculer plus tard aussi
202  gmm::copy(KS.mat(), W); gmm::copy(scaled(r, nss /beta), mat_col(W, m));
203 
204  // Computation of the oblique matrix
205  sub_interval SUBI(0, m);
206  add(scaled(sub_vector(ztest, SUBI), -Hobl(m, m-1) / ztest[m]),
207  sub_vector(mat_col(Hobl, m-1), SUBI));
208  Hobl(m, m-1) *= nss * beta / ztest[m];
209 
210  /* **************************************************************** */
211  /* Locking */
212  /* **************************************************************** */
213 
214  // Computation of the Ritz eigenpairs.
215  std::vector<std::complex<R> > eval(m);
216  dense_matrix<T> YB(m-st.tb_def, m-st.tb_def);
217  std::vector<char> pure(m-st.tb_def, 0);
218  gmm::clear(YB);
219 
220  select_eval(Hobl, eval, YB, pure, st);
221 
222  if (st.conv != 0) {
223  // DEFLATION using the QR Factorization of YB
224 
225  T alpha = Lock(W, Hobl,
226  sub_matrix(YB, sub_interval(0, m-st.tb_def)),
227  sub_interval(st.tb_def, m-st.tb_def),
228  (st.tb_defwant < p));
229  // ns *= alpha; // à calculer plus tard ??
230  // V(:,m+1) = alpha*V(:, m+1); ça devait servir à qlq chose ...
231 
232 
233  // Clean the portions below the diagonal corresponding
234  // to the lock Schur vectors
235 
236  for (size_type j = st.tb_def; j < st.tb_deftot; ++j) {
237  if ( pure[j-st.tb_def] == 0)
238  gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
239  else if (pure[j-st.tb_def] == 1) {
240  gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
241  sub_interval(j, 2)));
242  ++j;
243  }
244  else GMM_ASSERT3(false, "internal error");
245  }
246 
247  if (!st.ok) {
248 
249  // attention si m = 0;
250  size_type mm = std::min(k+st.nb_unwant+st.nb_nolong, m-1);
251 
252  if (eval_sort[m-mm-1].second != R(0)
253  && eval_sort[m-mm-1].second == -eval_sort[m-mm].second) ++mm;
254 
255  std::vector<complex<R> > shifts(m-mm);
256  for (size_type i = 0; i < m-mm; ++i)
257  shifts[i] = eval_sort[i].second;
258 
259  apply_shift_to_Arnoldi_factorization(W, Hobl, shifts, mm,
260  m-mm, true);
261 
262  st.fin = mm;
263  }
264  else
265  st.fin = st.tb_deftot;
266 
267 
268  /* ************************************************************** */
269  /* Purge */
270  /* ************************************************************** */
271 
272  if (st.nb_nolong + st.nb_unwant > 0) {
273 
274  std::vector<std::complex<R> > eval(m);
275  dense_matrix<T> YB(st.fin, st.tb_deftot);
276  std::vector<char> pure(st.tb_deftot, 0);
277  gmm::clear(YB);
278  st.nb_un = st.nb_nolong + st.nb_unwant;
279 
280  select_eval_for_purging(Hobl, eval, YB, pure, st);
281 
282  T alpha = Lock(W, Hobl, YB, sub_interval(0, st.fin), ok);
283 
284  // Clean the portions below the diagonal corresponding
285  // to the unwanted lock Schur vectors
286 
287  for (size_type j = 0; j < st.tb_deftot; ++j) {
288  if ( pure[j] == 0)
289  gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
290  else if (pure[j] == 1) {
291  gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
292  sub_interval(j, 2)));
293  ++j;
294  }
295  else GMM_ASSERT3(false, "internal error");
296  }
297 
298  gmm::dense_matrix<T> z(st.nb_un, st.fin - st.nb_un);
299  sub_interval SUBI(0, st.nb_un), SUBJ(st.nb_un, st.fin - st.nb_un);
300  sylvester(sub_matrix(Hobl, SUBI),
301  sub_matrix(Hobl, SUBJ),
302  sub_matrix(gmm::scaled(Hobl, -T(1)), SUBI, SUBJ), z);
303 
304  }
305 
306  }
307 
308  }
309  }
310  }
311 
312 
313  template < typename Mat, typename Vec, typename VecB, typename Precond >
314  void idgmres(const Mat &A, Vec &x, const VecB &b,
315  const Precond &M, size_type m, iteration& outer) {
316  typedef typename linalg_traits<Mat>::value_type T;
317  modified_gram_schmidt<T> orth(m, vect_size(x));
318  gmres(A, x, b, M, m, outer, orth);
319  }
320 
321 
322  // Lock stage of an implicit restarted Arnoldi process.
323  // 1- QR factorization of YB through Householder matrices
324  // Q(Rl) = YB
325  // (0 )
326  // 2- Update of the Arnoldi factorization.
327  // H <- Q*HQ, W <- WQ
328  // 3- Restore the Hessemberg form of H.
329 
330  template <typename T, typename MATYB>
331  void Lock(gmm::dense_matrix<T> &W, gmm::dense_matrix<T> &H,
332  const MATYB &YB, const sub_interval SUB,
333  bool restore, T &ns) {
334 
335  size_type n = mat_nrows(W), m = mat_ncols(W) - 1;
336  size_type ncols = mat_ncols(YB), nrows = mat_nrows(YB);
337  size_type begin = min(SUB); end = max(SUB) - 1;
338  sub_interval SUBR(0, nrows), SUBC(0, ncols);
339  T alpha(1);
340 
341  GMM_ASSERT2(((end-begin) == ncols) && (m == mat_nrows(H))
342  && (m+1 == mat_ncols(H)), "dimensions mismatch");
343 
344  // DEFLATION using the QR Factorization of YB
345 
346  dense_matrix<T> QR(n_rows, n_rows);
347  gmmm::copy(YB, sub_matrix(QR, SUBR, SUBC));
348  gmm::clear(submatrix(QR, SUBR, sub_interval(ncols, nrows-ncols)));
349  qr_factor(QR);
350 
351 
352  apply_house_left(QR, sub_matrix(H, SUB));
353  apply_house_right(QR, sub_matrix(H, SUBR, SUB));
354  apply_house_right(QR, sub_matrix(W, sub_interval(0, n), SUB));
355 
356  // Restore to the initial block hessenberg form
357 
358  if (restore) {
359 
360  // verifier quand m = 0 ...
361  gmm::dense_matrix tab_p(end - st.tb_deftot, end - st.tb_deftot);
362  gmm::copy(identity_matrix(), tab_p);
363 
364  for (size_type j = end-1; j >= st.tb_deftot+2; --j) {
365 
366  size_type jm = j-1;
367  std::vector<T> v(jm - st.tb_deftot);
368  sub_interval SUBtot(st.tb_deftot, jm - st.tb_deftot);
369  sub_interval SUBtot2(st.tb_deftot, end - st.tb_deftot);
370  gmm::copy(sub_vector(mat_row(H, j), SUBtot), v);
371  house_vector_last(v);
372  w.resize(end);
373  col_house_update(sub_matrix(H, SUBI, SUBtot), v, w);
374  w.resize(end - st.tb_deftot);
375  row_house_update(sub_matrix(H, SUBtot, SUBtot2), v, w);
376  gmm::clear(sub_vector(mat_row(H, j),
377  sub_interval(st.tb_deftot, j-1-st.tb_deftot)));
378  w.resize(end - st.tb_deftot);
379  col_house_update(sub_matrix(tab_p, sub_interval(0, end-st.tb_deftot),
380  sub_interval(0, jm-st.tb_deftot)), v, w);
381  w.resize(n);
382  col_house_update(sub_matrix(W, sub_interval(0, n), SUBtot), v, w);
383  }
384 
385  // restore positive subdiagonal elements
386 
387  std::vector<T> d(fin-st.tb_deftot); d[0] = T(1);
388 
389  // We compute d[i+1] in order
390  // (d[i+1] * H(st.tb_deftot+i+1,st.tb_deftoti)) / d[i]
391  // be equal to |H(st.tb_deftot+i+1,st.tb_deftot+i))|.
392  for (size_type j = 0; j+1 < end-st.tb_deftot; ++j) {
393  T e = H(st.tb_deftot+j, st.tb_deftot+j-1);
394  d[j+1] = (e == T(0)) ? T(1) : d[j] * gmm::abs(e) / e;
395  scale(sub_vector(mat_row(H, st.tb_deftot+j+1),
396  sub_interval(st.tb_deftot, m-st.tb_deftot)), d[j+1]);
397  scale(mat_col(H, st.tb_deftot+j+1), T(1) / d[j+1]);
398  scale(mat_col(W, st.tb_deftot+j+1), T(1) / d[j+1]);
399  }
400 
401  alpha = tab_p(end-st.tb_deftot-1, end-st.tb_deftot-1) / d[end-st.tb_deftot-1];
402  alpha /= gmm::abs(alpha);
403  scale(mat_col(W, m), alpha);
404 
405  }
406 
407  return alpha;
408  }
409 
410 
411 
412 
413 
414 
415 
416 
417  // Apply p implicit shifts to the Arnoldi factorization
418  // AV = VH+H(k+p+1,k+p) V(:,k+p+1) e_{k+p}*
419  // and produces the following new Arnoldi factorization
420  // A(VQ) = (VQ)(Q*HQ)+H(k+p+1,k+p) V(:,k+p+1) e_{k+p}* Q
421  // where only the first k columns are relevant.
422  //
423  // Dan Sorensen and Richard J. Radke, 11/95
424  template<typename T, typename C>
425  apply_shift_to_Arnoldi_factorization(dense_matrix<T> V, dense_matrix<T> H,
426  std::vector<C> Lambda, size_type &k,
427  size_type p, bool true_shift = false) {
428 
429 
430  size_type k1 = 0, num = 0, kend = k+p, kp1 = k + 1;
431  bool mark = false;
432  T c, s, x, y, z;
433 
434  dense_matrix<T> q(1, kend);
435  gmm::clear(q); q(0,kend-1) = T(1);
436  std::vector<T> hv(3), w(std::max(kend, mat_nrows(V)));
437 
438  for(size_type jj = 0; jj < p; ++jj) {
439  // compute and apply a bulge chase sweep initiated by the
440  // implicit shift held in w(jj)
441 
442  if (abs(Lambda[jj].real()) == 0.0) {
443  // apply a real shift using 2 by 2 Givens rotations
444 
445  for (size_type k1 = 0, k2 = 0; k2 != kend-1; k1 = k2+1) {
446  k2 = k1;
447  while (h(k2+1, k2) != T(0) && k2 < kend-1) ++k2;
448 
449  Givens_rotation(H(k1, k1) - Lambda[jj], H(k1+1, k1), c, s);
450 
451  for (i = k1; i <= k2; ++i) {
452  if (i > k1) Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
453 
454  // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
455  // Vérifier qu'au final H(i+1,i) est bien un réel positif.
456 
457  // apply rotation from left to rows of H
458  row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
459  c, s, 0, 0);
460 
461  // apply rotation from right to columns of H
462  size_type ip2 = std::min(i+2, kend);
463  col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
464  c, s, 0, 0);
465 
466  // apply rotation from right to columns of V
467  col_rot(V, c, s, i, i+1);
468 
469  // accumulate e' Q so residual can be updated k+p
470  Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
471  // peut être que nous utilisons G au lieu de G* et que
472  // nous allons trop loin en k2.
473  }
474  }
475 
476  num = num + 1;
477  }
478  else {
479 
480  // Apply a double complex shift using 3 by 3 Householder
481  // transformations
482 
483  if (jj == p || mark)
484  mark = false; // skip application of conjugate shift
485  else {
486  num = num + 2; // mark that a complex conjugate
487  mark = true; // pair has been applied
488 
489  // Indices de fin de boucle à surveiller... de près !
490  for (size_type k1 = 0, k3 = 0; k3 != kend-2; k1 = k3+1) {
491  k3 = k1;
492  while (h(k3+1, k3) != T(0) && k3 < kend-2) ++k3;
493  size_type k2 = k1+1;
494 
495 
496  x = H(k1,k1) * H(k1,k1) + H(k1,k2) * H(k2,k1)
497  - 2.0*Lambda[jj].real() * H(k1,k1) + gmm::abs_sqr(Lambda[jj]);
498  y = H(k2,k1) * (H(k1,k1) + H(k2,k2) - 2.0*Lambda[jj].real());
499  z = H(k2+1,k2) * H(k2,k1);
500 
501  for (size_type i = k1; i <= k3; ++i) {
502  if (i > k1) {
503  x = H(i, i-1);
504  y = H(i+1, i-1);
505  z = H(i+2, i-1);
506  // Ne pas oublier de nettoyer H(i+1,i-1) et H(i+2,i-1)
507  // (les mettre à zéro).
508  }
509 
510  hv[0] = x; hv[1] = y; hv[2] = z;
511  house_vector(v);
512 
513  // Vérifier qu'au final H(i+1,i) est bien un réel positif
514 
515  // apply transformation from left to rows of H
516  w.resize(kend-i);
517  row_house_update(sub_matrix(H, sub_interval(i, 2),
518  sub_interval(i, kend-i)), v, w);
519 
520  // apply transformation from right to columns of H
521 
522  size_type ip3 = std::min(kend, i + 3);
523  w.resize(ip3);
524  col_house_update(sub_matrix(H, sub_interval(0, ip3),
525  sub_interval(i, 2)), v, w);
526 
527  // apply transformation from right to columns of V
528 
529  w.resize(mat_nrows(V));
530  col_house_update(sub_matrix(V, sub_interval(0, mat_nrows(V)),
531  sub_interval(i, 2)), v, w);
532 
533  // accumulate e' Q so residual can be updated k+p
534 
535  w.resize(1);
536  col_house_update(sub_matrix(q, sub_interval(0,1),
537  sub_interval(i,2)), v, w);
538 
539  }
540  }
541 
542  // clean up step with Givens rotation
543 
544  i = kend-2;
545  c = x; s = y;
546  if (i > k1) Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
547 
548  // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
549  // Vérifier qu'au final H(i+1,i) est bien un réel positif.
550 
551  // apply rotation from left to rows of H
552  row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
553  c, s, 0, 0);
554 
555  // apply rotation from right to columns of H
556  size_type ip2 = std::min(i+2, kend);
557  col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
558  c, s, 0, 0);
559 
560  // apply rotation from right to columns of V
561  col_rot(V, c, s, i, i+1);
562 
563  // accumulate e' Q so residual can be updated k+p
564  Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
565 
566  }
567  }
568  }
569 
570  // update residual and store in the k+1 -st column of v
571 
572  k = kend - num;
573  scale(mat_col(V, kend), q(0, k));
574 
575  if (k < mat_nrows(H)) {
576  if (true_shift)
577  gmm::copy(mat_col(V, kend), mat_col(V, k));
578  else
579  // v(:,k+1) = v(:,kend+1) + v(:,k+1)*h(k+1,k);
580  // v(:,k+1) = v(:,kend+1) ;
581  gmm::add(scaled(mat_col(V, kend), H(kend, kend-1)),
582  scaled(mat_col(V, k), H(k, k-1)), mat_col(V, k));
583  }
584 
585  H(k, k-1) = vect_norm2(mat_col(V, k));
586  scale(mat_col(V, kend), T(1) / H(k, k-1));
587  }
588 
589 
590 
591  template<typename MAT, typename EVAL, typename PURE>
592  void select_eval(const MAT &Hobl, EVAL &eval, MAT &YB, PURE &pure,
593  idgmres_state &st) {
594 
595  typedef typename linalg_traits<MAT>::value_type T;
596  typedef typename number_traits<T>::magnitude_type R;
597  size_type m = st.m;
598 
599  // Computation of the Ritz eigenpairs.
600 
601  col_matrix< std::vector<T> > evect(m-st.tb_def, m-st.tb_def);
602  // std::vector<std::complex<R> > eval(m);
603  std::vector<R> ritznew(m, T(-1));
604 
605  // dense_matrix<T> evect_lock(st.tb_def, st.tb_def);
606 
607  sub_interval SUB1(st.tb_def, m-st.tb_def);
608  implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
609  sub_vector(eval, SUB1), evect);
610  sub_interval SUB2(0, st.tb_def);
611  implicit_qr_algorithm(sub_matrix(Hobl, SUB2),
612  sub_vector(eval, SUB2), /* evect_lock */);
613 
614  for (size_type l = st.tb_def; l < m; ++l)
615  ritznew[l] = gmm::abs(evect(m-st.tb_def-1, l-st.tb_def) * Hobl(m, m-1));
616 
617  std::vector< std::pair<T, size_type> > eval_sort(m);
618  for (size_type l = 0; l < m; ++l)
619  eval_sort[l] = std::pair<T, size_type>(eval[l], l);
620  std::sort(eval_sort.begin(), eval_sort.end(), compare_vp());
621 
622  std::vector<size_type> index(m);
623  for (size_type l = 0; l < m; ++l) index[l] = eval_sort[l].second;
624 
625  std::vector<bool> kept(m, false);
626  std::fill(kept.begin(), kept.begin()+st.tb_def, true);
627 
628  apply_permutation(eval, index);
629  apply_permutation(evect, index);
630  apply_permutation(ritznew, index);
631  apply_permutation(kept, index);
632 
633  // Which are the eigenvalues that converged ?
634  //
635  // nb_want is the number of eigenvalues of
636  // Hess(tb_def+1:n,tb_def+1:n) that converged and are WANTED
637  //
638  // nb_unwant is the number of eigenvalues of
639  // Hess(tb_def+1:n,tb_def+1:n) that converged and are UNWANTED
640  //
641  // nb_nolong is the number of eigenvalues of
642  // Hess(1:tb_def,1:tb_def) that are NO LONGER WANTED.
643  //
644  // tb_deftot is the number of the deflated eigenvalues
645  // that is tb_def + nb_want + nb_unwant
646  //
647  // tb_defwant is the number of the wanted deflated eigenvalues
648  // that is tb_def + nb_want - nb_nolong
649 
650  st.nb_want = 0, st.nb_unwant = 0, st.nb_nolong = 0;
651  size_type j, ind;
652 
653  for (j = 0, ind = 0; j < m-p; ++j) {
654  if (ritznew[j] == R(-1)) {
655  if (std::imag(eval[j]) != R(0)) {
656  st.nb_nolong += 2; ++j; // à adapter dans le cas complexe ...
657  }
658  else st.nb_nolong++;
659  }
660  else {
661  if (ritznew[j]
662  < tol_vp * gmm::abs(eval[j])) {
663 
664  for (size_type l = 0, l < m-st.tb_def; ++l)
665  YB(l, ind) = std::real(evect(l, j));
666  kept[j] = true;
667  ++j; ++st.nb_unwant; ind++;
668 
669  if (std::imag(eval[j]) != R(0)) {
670  for (size_type l = 0, l < m-st.tb_def; ++l)
671  YB(l, ind) = std::imag(evect(l, j));
672  pure[ind-1] = 1;
673  pure[ind] = 2;
674 
675  kept[j] = true;
676 
677  st.nb_unwant++;
678  ++ind;
679  }
680  }
681  }
682  }
683 
684 
685  for (; j < m; ++j) {
686  if (ritznew[j] != R(-1)) {
687 
688  for (size_type l = 0, l < m-st.tb_def; ++l)
689  YB(l, ind) = std::real(evect(l, j));
690  pure[ind] = 1;
691  ++ind;
692  kept[j] = true;
693  ++st.nb_want;
694 
695  if (ritznew[j]
696  < tol_vp * gmm::abs(eval[j])) {
697  for (size_type l = 0, l < m-st.tb_def; ++l)
698  YB(l, ind) = std::imag(evect(l, j));
699  pure[ind] = 2;
700 
701  j++;
702  kept[j] = true;
703 
704  st.nb_want++;
705  ++ind;
706  }
707  }
708  }
709 
710  std::vector<T> shift(m - st.tb_def - st.nb_want - st.nb_unwant);
711  for (size_type j = 0, i = 0; j < m; ++j)
712  if (!kept[j]) shift[i++] = eval[j];
713 
714  // st.conv (st.nb_want+st.nb_unwant) is the number of eigenpairs that
715  // have just converged.
716  // st.tb_deftot is the total number of eigenpairs that have converged.
717 
718  size_type st.conv = ind;
719  size_type st.tb_deftot = st.tb_def + st.conv;
720  size_type st.tb_defwant = st.tb_def + st.nb_want - st.nb_nolong;
721 
722  sub_interval SUBYB(0, st.conv);
723 
724  if ( st.tb_defwant >= p ) { // An invariant subspace has been found.
725 
726  st.nb_unwant = 0;
727  st.nb_want = p + st.nb_nolong - st.tb_def;
728  st.tb_defwant = p;
729 
730  if ( pure[st.conv - st.nb_want + 1] == 2 ) {
731  ++st.nb_want; st.tb_defwant = ++p;// il faudrait que ce soit un p local
732  }
733 
734  SUBYB = sub_interval(st.conv - st.nb_want, st.nb_want);
735  // YB = YB(:, st.conv-st.nb_want+1 : st.conv); // On laisse en suspend ..
736  // pure = pure(st.conv-st.nb_want+1 : st.conv,1); // On laisse suspend ..
737  st.conv = st.nb_want;
738  st.tb_deftot = st.tb_def + st.conv;
739  st.ok = true;
740  }
741 
742  }
743 
744 
745 
746  template<typename MAT, typename EVAL, typename PURE>
747  void select_eval_for_purging(const MAT &Hobl, EVAL &eval, MAT &YB,
748  PURE &pure, idgmres_state &st) {
749 
750  typedef typename linalg_traits<MAT>::value_type T;
751  typedef typename number_traits<T>::magnitude_type R;
752  size_type m = st.m;
753 
754  // Computation of the Ritz eigenpairs.
755 
756  col_matrix< std::vector<T> > evect(st.tb_deftot, st.tb_deftot);
757 
758  sub_interval SUB1(0, st.tb_deftot);
759  implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
760  sub_vector(eval, SUB1), evect);
761  std::fill(eval.begin() + st.tb_deftot, eval.end(), std::complex<R>(0));
762 
763  std::vector< std::pair<T, size_type> > eval_sort(m);
764  for (size_type l = 0; l < m; ++l)
765  eval_sort[l] = std::pair<T, size_type>(eval[l], l);
766  std::sort(eval_sort.begin(), eval_sort.end(), compare_vp());
767 
768  std::vector<bool> sorted(m);
769  std::fill(sorted.begin(), sorted.end(), false);
770 
771  std::vector<size_type> ind(m);
772  for (size_type l = 0; l < m; ++l) ind[l] = eval_sort[l].second;
773 
774  std::vector<bool> kept(m, false);
775  std::fill(kept.begin(), kept.begin()+st.tb_def, true);
776 
777  apply_permutation(eval, ind);
778  apply_permutation(evect, ind);
779 
780  size_type j;
781  for (j = 0; j < st.tb_deftot; ++j) {
782 
783  for (size_type l = 0, l < st.tb_deftot; ++l)
784  YB(l, j) = std::real(evect(l, j));
785 
786  if (std::imag(eval[j]) != R(0)) {
787  for (size_type l = 0, l < m-st.tb_def; ++l)
788  YB(l, j+1) = std::imag(evect(l, j));
789  pure[j] = 1;
790  pure[j+1] = 2;
791 
792  j += 2;
793  }
794  else ++j;
795  }
796  }
797 
798 
799 
800 
801 
802 
803 }
804 
805 #endif
gmm_dense_sylvester.h
Sylvester equation solver.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
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
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::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
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
gmm_kernel.h
Include the base gmm files.
gmm_iter.h
Iteration object.
gmm::add
void add(const L1 &l1, L2 &l2)
*â€/
Definition: gmm_blas.h:1268