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