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