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 {
operator()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>
apply_permutation(CONT & cont,const IND & 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 >
idgmres(const Mat & A,Vec & x,const VecB & b,const Precond & M,size_type m,size_type p,size_type k,double tol_vp,iteration & outer,Basis & KS)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 >
idgmres(const Mat & A,Vec & x,const VecB & b,const Precond & M,size_type m,iteration & outer)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>
Lock(gmm::dense_matrix<T> & W,gmm::dense_matrix<T> & H,const MATYB & YB,const sub_interval SUB,bool restore,T & ns)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>
select_eval(const MAT & Hobl,EVAL & eval,MAT & YB,PURE & pure,idgmres_state & st)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>
select_eval_for_purging(const MAT & Hobl,EVAL & eval,MAT & YB,PURE & pure,idgmres_state & st)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
806