[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Whitespace fixe
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Whitespace fixes |
|
Date: |
Wed, 20 Dec 2023 16:47:25 -0500 |
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch master
in repository getfem.
The following commit(s) were added to refs/heads/master by this push:
new dc63475e Whitespace fixes
dc63475e is described below
commit dc63475e720030ddbed8f967951d999fd61f9fe7
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Wed Dec 20 22:47:16 2023 +0100
Whitespace fixes
---
src/gmm/gmm_solver_idgmres.h | 869 +++++++++++++++++++++----------------------
1 file changed, 426 insertions(+), 443 deletions(-)
diff --git a/src/gmm/gmm_solver_idgmres.h b/src/gmm/gmm_solver_idgmres.h
index d675dc81..ba987095 100644
--- a/src/gmm/gmm_solver_idgmres.h
+++ b/src/gmm/gmm_solver_idgmres.h
@@ -46,7 +46,7 @@ namespace gmm {
template <typename T> compare_vp {
bool operator()(const std::pair<T, size_type> &a,
- const std::pair<T, size_type> &b) const
+ const std::pair<T, size_type> &b) const
{ return (gmm::abs(a.first) > gmm::abs(b.first)); }
}
@@ -57,35 +57,35 @@ namespace gmm {
idgmres_state(size_type mm, size_type pp, size_type kk)
: m(mm), tb_deb(1), tb_def(0), p(pp), k(kk), nb_want(0),
- nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
- conv(0), nb_un(0), fin(0), ok(false); {}
+ nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
+ conv(0), nb_un(0), fin(0), ok(false); {}
}
- idgmres_state(size_type mm, size_type pp, size_type kk)
- : m(mm), tb_deb(1), tb_def(0), p(pp), k(kk), nb_want(0),
- nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
- conv(0), nb_un(0), fin(0), ok(false); {}
-
+ idgmres_state(size_type mm, size_type pp, size_type kk)
+ : m(mm), tb_deb(1), tb_def(0), p(pp), k(kk), nb_want(0),
+ nb_unwant(0), nb_nolong(0), tb_deftot(0), tb_defwant(0),
+ conv(0), nb_un(0), fin(0), ok(false); {}
+
template <typename CONT, typename IND>
apply_permutation(CONT &cont, const IND &ind) {
size_type m = ind.end() - ind.begin();
std::vector<bool> sorted(m, false);
-
+
for (size_type l = 0; l < m; ++l)
if (!sorted[l] && ind[l] != l) {
- typeid(cont[0]) aux = cont[l];
- k = ind[l];
- cont[l] = cont[k];
- sorted[l] = true;
-
- for(k2 = ind[k]; k2 != l; k2 = ind[k]) {
- cont[k] = cont[k2];
- sorted[k] = true;
- k = k2;
- }
- cont[k] = aux;
+ typeid(cont[0]) aux = cont[l];
+ k = ind[l];
+ cont[l] = cont[k];
+ sorted[l] = true;
+
+ for(k2 = ind[k]; k2 != l; k2 = ind[k]) {
+ cont[k] = cont[k2];
+ sorted[k] = true;
+ k = k2;
+ }
+ cont[k] = aux;
}
}
@@ -95,7 +95,7 @@ namespace gmm {
See: C. Le Calvez, B. Molina, Implicitly restarted and deflated
FOM and GMRES, numerical applied mathematics,
(30) 2-3 (1999) pp191-212.
-
+
@param A Real or complex unsymmetric matrix.
@param x initial guess vector and final result.
@param b right hand side
@@ -109,14 +109,14 @@ namespace gmm {
@param KS
*/
template < typename Mat, typename Vec, typename VecB, typename Precond,
- typename Basis >
+ typename Basis >
void 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) {
+ size_type m, size_type p, size_type k, double tol_vp,
+ iteration &outer, Basis& KS) {
typedef typename linalg_traits<Mat>::value_type T;
typedef typename number_traits<T>::magnitude_type R;
-
+
R a, beta;
idgmres_state st(m, p, k);
@@ -125,13 +125,12 @@ namespace gmm {
std::vector<T> y(m+1), ztest(m+1), gam(m+1);
std::vector<T> gamma(m+1);
gmm::dense_matrix<T> H(m+1, m), Hess(m+1, m),
- Hobl(m+1, m), W(vect_size(x), m+1);
-
+ Hobl(m+1, m), W(vect_size(x), m+1);
gmm::clear(H);
outer.set_rhsnorm(gmm::vect_norm2(b));
if (outer.get_rhsnorm() == 0.0) { clear(x); return; }
-
+
mult(A, scaled(x, -1.0), b, w);
mult(M, w, r);
beta = gmm::vect_norm2(r);
@@ -140,9 +139,9 @@ namespace gmm {
inner.reduce_noisy();
inner.set_maxiter(m);
inner.set_name("GMRes inner iter");
-
+
while (! outer.finished(beta)) {
-
+
gmm::copy(gmm::scaled(r, 1.0/beta), KS[0]);
gmm::clear(s);
s[0] = beta;
@@ -150,47 +149,46 @@ namespace gmm {
inner.set_maxiter(m - st.tb_deb + 1);
size_type i = st.tb_deb - 1; inner.init();
-
+
do {
- mult(A, KS[i], u);
- mult(M, u, KS[i+1]);
- orthogonalize_with_refinment(KS, mat_col(H, i), i);
- H(i+1, i) = a = gmm::vect_norm2(KS[i+1]);
- gmm::scale(KS[i+1], R(1) / a);
-
- gmm::copy(mat_col(H, i), mat_col(Hess, i));
- gmm::copy(mat_col(H, i), mat_col(Hobl, i));
-
-
- for (size_type l = 0; l < i; ++l)
- Apply_Givens_rotation_left(H(l,i), H(l+1,i), c_rot[l], s_rot[l]);
-
- Givens_rotation(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
- Apply_Givens_rotation_left(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
- H(i+1, i) = T(0);
- Apply_Givens_rotation_left(s[i], s[i+1], c_rot[i], s_rot[i]);
-
- ++inner, ++outer, ++i;
+ mult(A, KS[i], u);
+ mult(M, u, KS[i+1]);
+ orthogonalize_with_refinment(KS, mat_col(H, i), i);
+ H(i+1, i) = a = gmm::vect_norm2(KS[i+1]);
+ gmm::scale(KS[i+1], R(1) / a);
+
+ gmm::copy(mat_col(H, i), mat_col(Hess, i));
+ gmm::copy(mat_col(H, i), mat_col(Hobl, i));
+
+ for (size_type l = 0; l < i; ++l)
+ Apply_Givens_rotation_left(H(l,i), H(l+1,i), c_rot[l], s_rot[l]);
+
+ Givens_rotation(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
+ Apply_Givens_rotation_left(H(i,i), H(i+1,i), c_rot[i], s_rot[i]);
+ H(i+1, i) = T(0);
+ Apply_Givens_rotation_left(s[i], s[i+1], c_rot[i], s_rot[i]);
+
+ ++inner, ++outer, ++i;
} while (! inner.finished(gmm::abs(s[i])));
if (inner.converged()) {
- gmm::copy(s, y);
- upper_tri_solve(H, y, i, false);
- combine(KS, y, x, i);
- mult(A, gmm::scaled(x, T(-1)), b, w);
- mult(M, w, r);
- beta = gmm::vect_norm2(r); // + verif sur beta ... à faire
- break;
+ gmm::copy(s, y);
+ upper_tri_solve(H, y, i, false);
+ combine(KS, y, x, i);
+ mult(A, gmm::scaled(x, T(-1)), b, w);
+ mult(M, w, r);
+ beta = gmm::vect_norm2(r); // + verif sur beta ... à faire
+ break;
}
gmm::clear(gam); gam[m] = s[i];
for (size_type l = m; l > 0; --l)
- Apply_Givens_rotation_left(gam[l-1], gam[l], gmm::conj(c_rot[l-1]),
- -s_rot[l-1]);
+ Apply_Givens_rotation_left(gam[l-1], gam[l], gmm::conj(c_rot[l-1]),
+ -s_rot[l-1]);
mult(KS.mat(), gam, r);
beta = gmm::vect_norm2(r);
-
+
mult(Hess, scaled(y, T(-1)), gamma, ztest);
// En fait, d'après Caroline qui s'y connait ztest et gam devrait
// être confondus
@@ -198,124 +196,119 @@ namespace gmm {
// place de ztest.
if (st.tb_def < p) {
T nss = H(m,m-1) / ztest[m];
- nss /= gmm::abs(nss); // ns à calculer plus tard aussi
- gmm::copy(KS.mat(), W); gmm::copy(scaled(r, nss /beta), mat_col(W, m));
-
- // Computation of the oblique matrix
- sub_interval SUBI(0, m);
- add(scaled(sub_vector(ztest, SUBI), -Hobl(m, m-1) / ztest[m]),
- sub_vector(mat_col(Hobl, m-1), SUBI));
- Hobl(m, m-1) *= nss * beta / ztest[m];
-
- /* **************************************************************** */
- /* Locking */
- /* **************************************************************** */
-
- // Computation of the Ritz eigenpairs.
- std::vector<std::complex<R> > eval(m);
- dense_matrix<T> YB(m-st.tb_def, m-st.tb_def);
- std::vector<char> pure(m-st.tb_def, 0);
- gmm::clear(YB);
-
- select_eval(Hobl, eval, YB, pure, st);
-
- if (st.conv != 0) {
- // DEFLATION using the QR Factorization of YB
-
- T alpha = Lock(W, Hobl,
- sub_matrix(YB, sub_interval(0, m-st.tb_def)),
- sub_interval(st.tb_def, m-st.tb_def),
- (st.tb_defwant < p));
- // ns *= alpha; // à calculer plus tard ??
- // V(:,m+1) = alpha*V(:, m+1); ça devait servir à qlq chose ...
-
-
- // Clean the portions below the diagonal corresponding
- // to the lock Schur vectors
-
- for (size_type j = st.tb_def; j < st.tb_deftot; ++j) {
- if ( pure[j-st.tb_def] == 0)
- gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
- else if (pure[j-st.tb_def] == 1) {
- gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
- sub_interval(j, 2)));
- ++j;
- }
- else GMM_ASSERT3(false, "internal error");
- }
-
- if (!st.ok) {
-
- // attention si m = 0;
- size_type mm = std::min(k+st.nb_unwant+st.nb_nolong, m-1);
-
- if (eval_sort[m-mm-1].second != R(0)
- && eval_sort[m-mm-1].second == -eval_sort[m-mm].second) ++mm;
-
- std::vector<complex<R> > shifts(m-mm);
- for (size_type i = 0; i < m-mm; ++i)
- shifts[i] = eval_sort[i].second;
-
- apply_shift_to_Arnoldi_factorization(W, Hobl, shifts, mm,
- m-mm, true);
-
- st.fin = mm;
- }
- else
- st.fin = st.tb_deftot;
-
-
- /* ************************************************************** */
- /* Purge */
- /* ************************************************************** */
-
- if (st.nb_nolong + st.nb_unwant > 0) {
-
- std::vector<std::complex<R> > eval(m);
- dense_matrix<T> YB(st.fin, st.tb_deftot);
- std::vector<char> pure(st.tb_deftot, 0);
- gmm::clear(YB);
- st.nb_un = st.nb_nolong + st.nb_unwant;
-
- select_eval_for_purging(Hobl, eval, YB, pure, st);
-
- T alpha = Lock(W, Hobl, YB, sub_interval(0, st.fin), ok);
-
- // Clean the portions below the diagonal corresponding
- // to the unwanted lock Schur vectors
-
- for (size_type j = 0; j < st.tb_deftot; ++j) {
- if ( pure[j] == 0)
- gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
- else if (pure[j] == 1) {
- gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
- sub_interval(j, 2)));
- ++j;
- }
- else GMM_ASSERT3(false, "internal error");
- }
-
- gmm::dense_matrix<T> z(st.nb_un, st.fin - st.nb_un);
- sub_interval SUBI(0, st.nb_un), SUBJ(st.nb_un, st.fin - st.nb_un);
- sylvester(sub_matrix(Hobl, SUBI),
- sub_matrix(Hobl, SUBJ),
- sub_matrix(gmm::scaled(Hobl, -T(1)), SUBI, SUBJ), z);
-
- }
-
- }
-
+ nss /= gmm::abs(nss); // ns à calculer plus tard aussi
+ gmm::copy(KS.mat(), W); gmm::copy(scaled(r, nss /beta), mat_col(W, m));
+
+ // Computation of the oblique matrix
+ sub_interval SUBI(0, m);
+ add(scaled(sub_vector(ztest, SUBI), -Hobl(m, m-1) / ztest[m]),
+ sub_vector(mat_col(Hobl, m-1), SUBI));
+ Hobl(m, m-1) *= nss * beta / ztest[m];
+
+ /* **************************************************************** */
+ /* Locking */
+ /* **************************************************************** */
+
+ // Computation of the Ritz eigenpairs.
+ std::vector<std::complex<R> > eval(m);
+ dense_matrix<T> YB(m-st.tb_def, m-st.tb_def);
+ std::vector<char> pure(m-st.tb_def, 0);
+ gmm::clear(YB);
+
+ select_eval(Hobl, eval, YB, pure, st);
+
+ if (st.conv != 0) {
+ // DEFLATION using the QR Factorization of YB
+
+ T alpha = Lock(W, Hobl,
+ sub_matrix(YB, sub_interval(0, m-st.tb_def)),
+ sub_interval(st.tb_def, m-st.tb_def),
+ (st.tb_defwant < p));
+ // ns *= alpha; // à calculer plus tard ??
+ // V(:,m+1) = alpha*V(:, m+1); ça devait servir à qlq chose ...
+
+
+ // Clean the portions below the diagonal corresponding
+ // to the lock Schur vectors
+ for (size_type j = st.tb_def; j < st.tb_deftot; ++j) {
+ if ( pure[j-st.tb_def] == 0)
+ gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
+ else if (pure[j-st.tb_def] == 1) {
+ gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
+ sub_interval(j, 2)));
+ ++j;
+ }
+ else GMM_ASSERT3(false, "internal error");
+ }
+
+ if (!st.ok) {
+ // attention si m = 0;
+ size_type mm = std::min(k+st.nb_unwant+st.nb_nolong, m-1);
+
+ if (eval_sort[m-mm-1].second != R(0)
+ && eval_sort[m-mm-1].second == -eval_sort[m-mm].second)
+ ++mm;
+
+ std::vector<complex<R> > shifts(m-mm);
+ for (size_type i = 0; i < m-mm; ++i)
+ shifts[i] = eval_sort[i].second;
+
+ apply_shift_to_Arnoldi_factorization(W, Hobl, shifts, mm,
+ m-mm, true);
+
+ st.fin = mm;
+ }
+ else
+ st.fin = st.tb_deftot;
+
+
+ /* ************************************************************** */
+ /* Purge */
+ /* ************************************************************** */
+
+ if (st.nb_nolong + st.nb_unwant > 0) {
+
+ std::vector<std::complex<R> > eval(m);
+ dense_matrix<T> YB(st.fin, st.tb_deftot);
+ std::vector<char> pure(st.tb_deftot, 0);
+ gmm::clear(YB);
+ st.nb_un = st.nb_nolong + st.nb_unwant;
+
+ select_eval_for_purging(Hobl, eval, YB, pure, st);
+
+ T alpha = Lock(W, Hobl, YB, sub_interval(0, st.fin), ok);
+
+ // Clean the portions below the diagonal corresponding
+ // to the unwanted lock Schur vectors
+ for (size_type j = 0; j < st.tb_deftot; ++j) {
+ if ( pure[j] == 0)
+ gmm::clear(sub_vector(mat_col(Hobl,j), sub_interval(j+1,m-j)));
+ else if (pure[j] == 1) {
+ gmm::clear(sub_matrix(Hobl, sub_interval(j+2,m-j-1),
+ sub_interval(j, 2)));
+ ++j;
+ }
+ else GMM_ASSERT3(false, "internal error");
+ }
+
+ gmm::dense_matrix<T> z(st.nb_un, st.fin - st.nb_un);
+ sub_interval SUBI(0, st.nb_un), SUBJ(st.nb_un, st.fin - st.nb_un);
+ sylvester(sub_matrix(Hobl, SUBI),
+ sub_matrix(Hobl, SUBJ),
+ sub_matrix(gmm::scaled(Hobl, -T(1)), SUBI, SUBJ), z);
+ }
+ }
}
}
}
-
+
template < typename Mat, typename Vec, typename VecB, typename Precond >
void idgmres(const Mat &A, Vec &x, const VecB &b,
- const Precond &M, size_type m, iteration& outer) {
+ const Precond &M, size_type m, iteration& outer) {
typedef typename linalg_traits<Mat>::value_type T;
modified_gram_schmidt<T> orth(m, vect_size(x));
- gmres(A, x, b, M, m, outer, orth);
+ gmres(A, x, b, M, m, outer, orth);
}
@@ -328,9 +321,9 @@ namespace gmm {
// 3- Restore the Hessemberg form of H.
template <typename T, typename MATYB>
- void Lock(gmm::dense_matrix<T> &W, gmm::dense_matrix<T> &H,
- const MATYB &YB, const sub_interval SUB,
- bool restore, T &ns) {
+ void Lock(gmm::dense_matrix<T> &W, gmm::dense_matrix<T> &H,
+ const MATYB &YB, const sub_interval SUB,
+ bool restore, T &ns) {
size_type n = mat_nrows(W), m = mat_ncols(W) - 1;
size_type ncols = mat_ncols(YB), nrows = mat_nrows(YB);
@@ -338,82 +331,74 @@ namespace gmm {
sub_interval SUBR(0, nrows), SUBC(0, ncols);
T alpha(1);
- GMM_ASSERT2(((end-begin) == ncols) && (m == mat_nrows(H))
- && (m+1 == mat_ncols(H)), "dimensions mismatch");
-
+ GMM_ASSERT2(((end-begin) == ncols) && (m == mat_nrows(H))
+ && (m+1 == mat_ncols(H)), "dimensions mismatch");
+
// DEFLATION using the QR Factorization of YB
-
+
dense_matrix<T> QR(n_rows, n_rows);
gmmm::copy(YB, sub_matrix(QR, SUBR, SUBC));
gmm::clear(submatrix(QR, SUBR, sub_interval(ncols, nrows-ncols)));
- qr_factor(QR);
-
+ qr_factor(QR);
apply_house_left(QR, sub_matrix(H, SUB));
apply_house_right(QR, sub_matrix(H, SUBR, SUB));
apply_house_right(QR, sub_matrix(W, sub_interval(0, n), SUB));
-
+
// Restore to the initial block hessenberg form
-
+
if (restore) {
-
+
// verifier quand m = 0 ...
gmm::dense_matrix tab_p(end - st.tb_deftot, end - st.tb_deftot);
gmm::copy(identity_matrix(), tab_p);
-
+
for (size_type j = end-1; j >= st.tb_deftot+2; --j) {
-
- size_type jm = j-1;
- std::vector<T> v(jm - st.tb_deftot);
- sub_interval SUBtot(st.tb_deftot, jm - st.tb_deftot);
- sub_interval SUBtot2(st.tb_deftot, end - st.tb_deftot);
- gmm::copy(sub_vector(mat_row(H, j), SUBtot), v);
- house_vector_last(v);
- w.resize(end);
- col_house_update(sub_matrix(H, SUBI, SUBtot), v, w);
- w.resize(end - st.tb_deftot);
- row_house_update(sub_matrix(H, SUBtot, SUBtot2), v, w);
- gmm::clear(sub_vector(mat_row(H, j),
- sub_interval(st.tb_deftot, j-1-st.tb_deftot)));
- w.resize(end - st.tb_deftot);
- col_house_update(sub_matrix(tab_p, sub_interval(0, end-st.tb_deftot),
- sub_interval(0, jm-st.tb_deftot)), v, w);
- w.resize(n);
- col_house_update(sub_matrix(W, sub_interval(0, n), SUBtot), v, w);
+
+ size_type jm = j-1;
+ std::vector<T> v(jm - st.tb_deftot);
+ sub_interval SUBtot(st.tb_deftot, jm - st.tb_deftot);
+ sub_interval SUBtot2(st.tb_deftot, end - st.tb_deftot);
+ gmm::copy(sub_vector(mat_row(H, j), SUBtot), v);
+ house_vector_last(v);
+ w.resize(end);
+ col_house_update(sub_matrix(H, SUBI, SUBtot), v, w);
+ w.resize(end - st.tb_deftot);
+ row_house_update(sub_matrix(H, SUBtot, SUBtot2), v, w);
+ gmm::clear(sub_vector(mat_row(H, j),
+ sub_interval(st.tb_deftot, j-1-st.tb_deftot)));
+ w.resize(end - st.tb_deftot);
+ col_house_update(sub_matrix(tab_p, sub_interval(0, end-st.tb_deftot),
+ sub_interval(0, jm-st.tb_deftot)), v, w);
+ w.resize(n);
+ col_house_update(sub_matrix(W, sub_interval(0, n), SUBtot), v, w);
}
-
+
// restore positive subdiagonal elements
-
+
std::vector<T> d(fin-st.tb_deftot); d[0] = T(1);
-
- // We compute d[i+1] in order
- // (d[i+1] * H(st.tb_deftot+i+1,st.tb_deftoti)) / d[i]
+
+ // We compute d[i+1] in order
+ // (d[i+1] * H(st.tb_deftot+i+1,st.tb_deftoti)) / d[i]
// be equal to |H(st.tb_deftot+i+1,st.tb_deftot+i))|.
for (size_type j = 0; j+1 < end-st.tb_deftot; ++j) {
- T e = H(st.tb_deftot+j, st.tb_deftot+j-1);
- d[j+1] = (e == T(0)) ? T(1) : d[j] * gmm::abs(e) / e;
- scale(sub_vector(mat_row(H, st.tb_deftot+j+1),
- sub_interval(st.tb_deftot, m-st.tb_deftot)), d[j+1]);
- scale(mat_col(H, st.tb_deftot+j+1), T(1) / d[j+1]);
- scale(mat_col(W, st.tb_deftot+j+1), T(1) / d[j+1]);
+ T e = H(st.tb_deftot+j, st.tb_deftot+j-1);
+ d[j+1] = (e == T(0)) ? T(1) : d[j] * gmm::abs(e) / e;
+ scale(sub_vector(mat_row(H, st.tb_deftot+j+1),
+ sub_interval(st.tb_deftot, m-st.tb_deftot)), d[j+1]);
+ scale(mat_col(H, st.tb_deftot+j+1), T(1) / d[j+1]);
+ scale(mat_col(W, st.tb_deftot+j+1), T(1) / d[j+1]);
}
alpha = tab_p(end-st.tb_deftot-1, end-st.tb_deftot-1) /
d[end-st.tb_deftot-1];
alpha /= gmm::abs(alpha);
scale(mat_col(W, m), alpha);
-
}
-
+
return alpha;
}
-
-
-
-
-
-
// Apply p implicit shifts to the Arnoldi factorization
// AV = VH+H(k+p+1,k+p) V(:,k+p+1) e_{k+p}*
// and produces the following new Arnoldi factorization
@@ -422,10 +407,9 @@ namespace gmm {
//
// Dan Sorensen and Richard J. Radke, 11/95
template<typename T, typename C>
- apply_shift_to_Arnoldi_factorization(dense_matrix<T> V, dense_matrix<T> H,
- std::vector<C> Lambda, size_type &k,
- size_type p, bool true_shift = false) {
-
+ apply_shift_to_Arnoldi_factorization(dense_matrix<T> V, dense_matrix<T> H,
+ std::vector<C> Lambda, size_type &k,
+ size_type p, bool true_shift = false) {
size_type k1 = 0, num = 0, kend = k+p, kp1 = k + 1;
bool mark = false;
@@ -438,132 +422,138 @@ namespace gmm {
for(size_type jj = 0; jj < p; ++jj) {
// compute and apply a bulge chase sweep initiated by the
// implicit shift held in w(jj)
-
+
if (abs(Lambda[jj].real()) == 0.0) {
- // apply a real shift using 2 by 2 Givens rotations
+ // apply a real shift using 2 by 2 Givens rotations
- for (size_type k1 = 0, k2 = 0; k2 != kend-1; k1 = k2+1) {
- k2 = k1;
- while (h(k2+1, k2) != T(0) && k2 < kend-1) ++k2;
+ for (size_type k1 = 0, k2 = 0; k2 != kend-1; k1 = k2+1) {
+ k2 = k1;
+ while (h(k2+1, k2) != T(0) && k2 < kend-1)
+ ++k2;
- Givens_rotation(H(k1, k1) - Lambda[jj], H(k1+1, k1), c, s);
-
- for (i = k1; i <= k2; ++i) {
- if (i > k1) Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
-
- // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
- // Vérifier qu'au final H(i+1,i) est bien un réel positif.
+ Givens_rotation(H(k1, k1) - Lambda[jj], H(k1+1, k1), c, s);
+
+ for (i = k1; i <= k2; ++i) {
+ if (i > k1)
+ Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
+
+ // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
+ // Vérifier qu'au final H(i+1,i) est bien un réel positif.
// apply rotation from left to rows of H
- row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
- c, s, 0, 0);
-
- // apply rotation from right to columns of H
+ row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
+ c, s, 0, 0);
+
+ // apply rotation from right to columns of H
size_type ip2 = std::min(i+2, kend);
col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
- c, s, 0, 0);
-
+ c, s, 0, 0);
+
// apply rotation from right to columns of V
- col_rot(V, c, s, i, i+1);
-
+ col_rot(V, c, s, i, i+1);
+
// accumulate e' Q so residual can be updated k+p
- Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
- // peut être que nous utilisons G au lieu de G* et que
- // nous allons trop loin en k2.
- }
- }
-
- num = num + 1;
+ Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
+ // peut être que nous utilisons G au lieu de G* et que
+ // nous allons trop loin en k2.
+ }
+ }
+
+ num = num + 1;
}
else {
-
- // Apply a double complex shift using 3 by 3 Householder
- // transformations
-
- if (jj == p || mark)
- mark = false; // skip application of conjugate shift
- else {
- num = num + 2; // mark that a complex conjugate
- mark = true; // pair has been applied
-
- // Indices de fin de boucle à surveiller... de près !
- for (size_type k1 = 0, k3 = 0; k3 != kend-2; k1 = k3+1) {
- k3 = k1;
- while (h(k3+1, k3) != T(0) && k3 < kend-2) ++k3;
- size_type k2 = k1+1;
+
+ // Apply a double complex shift using 3 by 3 Householder
+ // transformations
+
+ if (jj == p || mark)
+ mark = false; // skip application of conjugate shift
+ else {
+ num = num + 2; // mark that a complex conjugate
+ mark = true; // pair has been applied
+
+ // Indices de fin de boucle à surveiller... de près !
+ for (size_type k1 = 0, k3 = 0; k3 != kend-2; k1 = k3+1) {
+ k3 = k1;
+ while (h(k3+1, k3) != T(0) && k3 < kend-2)
+ ++k3;
+ size_type k2 = k1+1;
x = H(k1,k1) * H(k1,k1) + H(k1,k2) * H(k2,k1)
- - 2.0*Lambda[jj].real() * H(k1,k1) + gmm::abs_sqr(Lambda[jj]);
- y = H(k2,k1) * (H(k1,k1) + H(k2,k2) - 2.0*Lambda[jj].real());
- z = H(k2+1,k2) * H(k2,k1);
-
- for (size_type i = k1; i <= k3; ++i) {
- if (i > k1) {
- x = H(i, i-1);
- y = H(i+1, i-1);
- z = H(i+2, i-1);
- // Ne pas oublier de nettoyer H(i+1,i-1) et H(i+2,i-1)
- // (les mettre à zéro).
- }
-
- hv[0] = x; hv[1] = y; hv[2] = z;
- house_vector(v);
-
- // Vérifier qu'au final H(i+1,i) est bien un réel positif
-
- // apply transformation from left to rows of H
- w.resize(kend-i);
- row_house_update(sub_matrix(H, sub_interval(i, 2),
- sub_interval(i, kend-i)), v, w);
-
- // apply transformation from right to columns of H
-
- size_type ip3 = std::min(kend, i + 3);
- w.resize(ip3);
+ - 2.0*Lambda[jj].real() * H(k1,k1) + gmm::abs_sqr(Lambda[jj]);
+ y = H(k2,k1) * (H(k1,k1) + H(k2,k2) - 2.0*Lambda[jj].real());
+ z = H(k2+1,k2) * H(k2,k1);
+
+ for (size_type i = k1; i <= k3; ++i) {
+ if (i > k1) {
+ x = H(i, i-1);
+ y = H(i+1, i-1);
+ z = H(i+2, i-1);
+ // Ne pas oublier de nettoyer H(i+1,i-1) et H(i+2,i-1)
+ // (les mettre à zéro).
+ }
+
+ hv[0] = x; hv[1] = y; hv[2] = z;
+ house_vector(v);
+
+ // Vérifier qu'au final H(i+1,i) est bien un réel positif
+
+ // apply transformation from left to rows of H
+ w.resize(kend-i);
+ row_house_update(sub_matrix(H, sub_interval(i, 2),
+ sub_interval(i, kend-i)),
+ v, w);
+
+ // apply transformation from right to columns of H
+
+ size_type ip3 = std::min(kend, i + 3);
+ w.resize(ip3);
col_house_update(sub_matrix(H, sub_interval(0, ip3),
- sub_interval(i, 2)), v, w);
-
- // apply transformation from right to columns of V
-
- w.resize(mat_nrows(V));
- col_house_update(sub_matrix(V, sub_interval(0, mat_nrows(V)),
- sub_interval(i, 2)), v, w);
-
- // accumulate e' Q so residual can be updated k+p
-
- w.resize(1);
- col_house_update(sub_matrix(q, sub_interval(0,1),
- sub_interval(i,2)), v, w);
-
- }
- }
-
- // clean up step with Givens rotation
-
- i = kend-2;
- c = x; s = y;
- if (i > k1) Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
-
- // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
- // Vérifier qu'au final H(i+1,i) est bien un réel positif.
-
- // apply rotation from left to rows of H
- row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
- c, s, 0, 0);
-
- // apply rotation from right to columns of H
- size_type ip2 = std::min(i+2, kend);
- col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
- c, s, 0, 0);
-
- // apply rotation from right to columns of V
- col_rot(V, c, s, i, i+1);
-
- // accumulate e' Q so residual can be updated k+p
- Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
-
- }
+ sub_interval(i, 2)),
+ v, w);
+
+ // apply transformation from right to columns of V
+
+ w.resize(mat_nrows(V));
+ col_house_update(sub_matrix(V, sub_interval(0, mat_nrows(V)),
+ sub_interval(i, 2)),
+ v, w);
+
+ // accumulate e' Q so residual can be updated k+p
+ w.resize(1);
+ col_house_update(sub_matrix(q, sub_interval(0,1),
+ sub_interval(i,2)),
+ v, w);
+ }
+ }
+
+ // clean up step with Givens rotation
+
+ i = kend-2;
+ c = x;
+ s = y;
+ if (i > k1)
+ Givens_rotation(H(i, i-1), H(i+1, i-1), c, s);
+
+ // Ne pas oublier de nettoyer H(i+1,i-1) (le mettre à zéro).
+ // Vérifier qu'au final H(i+1,i) est bien un réel positif.
+
+ // apply rotation from left to rows of H
+ row_rot(sub_matrix(H, sub_interval(i,2), sub_interval(i, kend-i)),
+ c, s, 0, 0);
+
+ // apply rotation from right to columns of H
+ size_type ip2 = std::min(i+2, kend);
+ col_rot(sub_matrix(H, sub_interval(0, ip2), sub_interval(i, 2))
+ c, s, 0, 0);
+
+ // apply rotation from right to columns of V
+ col_rot(V, c, s, i, i+1);
+
+ // accumulate e' Q so residual can be updated k+p
+ Apply_Givens_rotation_left(q(0,i), q(0,i+1), c, s);
+ }
}
}
@@ -571,15 +561,15 @@ namespace gmm {
k = kend - num;
scale(mat_col(V, kend), q(0, k));
-
+
if (k < mat_nrows(H)) {
if (true_shift)
- gmm::copy(mat_col(V, kend), mat_col(V, k));
+ gmm::copy(mat_col(V, kend), mat_col(V, k));
else
- // v(:,k+1) = v(:,kend+1) + v(:,k+1)*h(k+1,k);
- // v(:,k+1) = v(:,kend+1) ;
- gmm::add(scaled(mat_col(V, kend), H(kend, kend-1)),
- scaled(mat_col(V, k), H(k, k-1)), mat_col(V, k));
+ // v(:,k+1) = v(:,kend+1) + v(:,k+1)*h(k+1,k);
+ // v(:,k+1) = v(:,kend+1) ;
+ gmm::add(scaled(mat_col(V, kend), H(kend, kend-1)),
+ scaled(mat_col(V, k), H(k, k-1)), mat_col(V, k));
}
H(k, k-1) = vect_norm2(mat_col(V, k));
@@ -590,38 +580,38 @@ namespace gmm {
template<typename MAT, typename EVAL, typename PURE>
void select_eval(const MAT &Hobl, EVAL &eval, MAT &YB, PURE &pure,
- idgmres_state &st) {
+ idgmres_state &st) {
typedef typename linalg_traits<MAT>::value_type T;
typedef typename number_traits<T>::magnitude_type R;
size_type m = st.m;
// Computation of the Ritz eigenpairs.
-
+
col_matrix< std::vector<T> > evect(m-st.tb_def, m-st.tb_def);
// std::vector<std::complex<R> > eval(m);
std::vector<R> ritznew(m, T(-1));
-
+
// dense_matrix<T> evect_lock(st.tb_def, st.tb_def);
-
+
sub_interval SUB1(st.tb_def, m-st.tb_def);
implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
- sub_vector(eval, SUB1), evect);
+ sub_vector(eval, SUB1), evect);
sub_interval SUB2(0, st.tb_def);
implicit_qr_algorithm(sub_matrix(Hobl, SUB2),
- sub_vector(eval, SUB2), /* evect_lock */);
-
+ sub_vector(eval, SUB2), /* evect_lock */);
+
for (size_type l = st.tb_def; l < m; ++l)
ritznew[l] = gmm::abs(evect(m-st.tb_def-1, l-st.tb_def) * Hobl(m, m-1));
-
+
std::vector< std::pair<T, size_type> > eval_sort(m);
for (size_type l = 0; l < m; ++l)
eval_sort[l] = std::pair<T, size_type>(eval[l], l);
std::sort(eval_sort.begin(), eval_sort.end(), compare_vp());
-
+
std::vector<size_type> index(m);
for (size_type l = 0; l < m; ++l) index[l] = eval_sort[l].second;
-
+
std::vector<bool> kept(m, false);
std::fill(kept.begin(), kept.begin()+st.tb_def, true);
@@ -630,107 +620,105 @@ namespace gmm {
apply_permutation(ritznew, index);
apply_permutation(kept, index);
- // Which are the eigenvalues that converged ?
+ // Which are the eigenvalues that converged ?
//
- // nb_want is the number of eigenvalues of
- // Hess(tb_def+1:n,tb_def+1:n) that converged and are WANTED
+ // nb_want is the number of eigenvalues of
+ // Hess(tb_def+1:n,tb_def+1:n) that converged and are WANTED
//
- // nb_unwant is the number of eigenvalues of
- // Hess(tb_def+1:n,tb_def+1:n) that converged and are UNWANTED
+ // nb_unwant is the number of eigenvalues of
+ // Hess(tb_def+1:n,tb_def+1:n) that converged and are UNWANTED
//
- // nb_nolong is the number of eigenvalues of
- // Hess(1:tb_def,1:tb_def) that are NO LONGER WANTED.
+ // nb_nolong is the number of eigenvalues of
+ // Hess(1:tb_def,1:tb_def) that are NO LONGER WANTED.
//
- // tb_deftot is the number of the deflated eigenvalues
- // that is tb_def + nb_want + nb_unwant
+ // tb_deftot is the number of the deflated eigenvalues
+ // that is tb_def + nb_want + nb_unwant
//
- // tb_defwant is the number of the wanted deflated eigenvalues
- // that is tb_def + nb_want - nb_nolong
-
+ // tb_defwant is the number of the wanted deflated eigenvalues
+ // that is tb_def + nb_want - nb_nolong
+
st.nb_want = 0, st.nb_unwant = 0, st.nb_nolong = 0;
size_type j, ind;
-
+
for (j = 0, ind = 0; j < m-p; ++j) {
if (ritznew[j] == R(-1)) {
- if (std::imag(eval[j]) != R(0)) {
- st.nb_nolong += 2; ++j; // à adapter dans le cas complexe ...
- }
- else st.nb_nolong++;
- }
- else {
- if (ritznew[j]
- < tol_vp * gmm::abs(eval[j])) {
-
- for (size_type l = 0, l < m-st.tb_def; ++l)
- YB(l, ind) = std::real(evect(l, j));
- kept[j] = true;
- ++j; ++st.nb_unwant; ind++;
-
- if (std::imag(eval[j]) != R(0)) {
- for (size_type l = 0, l < m-st.tb_def; ++l)
- YB(l, ind) = std::imag(evect(l, j));
- pure[ind-1] = 1;
- pure[ind] = 2;
-
- kept[j] = true;
-
- st.nb_unwant++;
- ++ind;
- }
- }
+ if (std::imag(eval[j]) != R(0)) {
+ st.nb_nolong += 2; ++j; // à adapter dans le cas complexe ...
+ } else
+ st.nb_nolong++;
+ } else {
+ if (ritznew[j] < tol_vp * gmm::abs(eval[j])) {
+
+ for (size_type l = 0, l < m-st.tb_def; ++l)
+ YB(l, ind) = std::real(evect(l, j));
+ kept[j] = true;
+ ++j; ++st.nb_unwant; ind++;
+
+ if (std::imag(eval[j]) != R(0)) {
+ for (size_type l = 0, l < m-st.tb_def; ++l)
+ YB(l, ind) = std::imag(evect(l, j));
+ pure[ind-1] = 1;
+ pure[ind] = 2;
+
+ kept[j] = true;
+
+ st.nb_unwant++;
+ ++ind;
+ }
+ }
}
}
-
-
+
+
for (; j < m; ++j) {
if (ritznew[j] != R(-1)) {
- for (size_type l = 0, l < m-st.tb_def; ++l)
- YB(l, ind) = std::real(evect(l, j));
- pure[ind] = 1;
- ++ind;
- kept[j] = true;
- ++st.nb_want;
-
- if (ritznew[j]
- < tol_vp * gmm::abs(eval[j])) {
- for (size_type l = 0, l < m-st.tb_def; ++l)
- YB(l, ind) = std::imag(evect(l, j));
- pure[ind] = 2;
-
- j++;
- kept[j] = true;
-
- st.nb_want++;
- ++ind;
- }
+ for (size_type l = 0, l < m-st.tb_def; ++l)
+ YB(l, ind) = std::real(evect(l, j));
+ pure[ind] = 1;
+ ++ind;
+ kept[j] = true;
+ ++st.nb_want;
+
+ if (ritznew[j] < tol_vp * gmm::abs(eval[j])) {
+ for (size_type l = 0, l < m-st.tb_def; ++l)
+ YB(l, ind) = std::imag(evect(l, j));
+ pure[ind] = 2;
+
+ j++;
+ kept[j] = true;
+
+ st.nb_want++;
+ ++ind;
+ }
}
}
-
+
std::vector<T> shift(m - st.tb_def - st.nb_want - st.nb_unwant);
for (size_type j = 0, i = 0; j < m; ++j)
- if (!kept[j]) shift[i++] = eval[j];
-
+ if (!kept[j])
+ shift[i++] = eval[j];
+
// st.conv (st.nb_want+st.nb_unwant) is the number of eigenpairs that
// have just converged.
// st.tb_deftot is the total number of eigenpairs that have converged.
-
+
size_type st.conv = ind;
size_type st.tb_deftot = st.tb_def + st.conv;
size_type st.tb_defwant = st.tb_def + st.nb_want - st.nb_nolong;
-
+
sub_interval SUBYB(0, st.conv);
-
+
if ( st.tb_defwant >= p ) { // An invariant subspace has been found.
-
+
st.nb_unwant = 0;
st.nb_want = p + st.nb_nolong - st.tb_def;
st.tb_defwant = p;
-
+
if ( pure[st.conv - st.nb_want + 1] == 2 ) {
- ++st.nb_want; st.tb_defwant = ++p;// il faudrait que ce soit un p local
+ ++st.nb_want; st.tb_defwant = ++p;// il faudrait que ce soit un p local
}
-
+
SUBYB = sub_interval(st.conv - st.nb_want, st.nb_want);
// YB = YB(:, st.conv-st.nb_want+1 : st.conv); // On laisse en suspend ..
// pure = pure(st.conv-st.nb_want+1 : st.conv,1); // On laisse suspend ..
@@ -738,28 +726,27 @@ namespace gmm {
st.tb_deftot = st.tb_def + st.conv;
st.ok = true;
}
-
}
template<typename MAT, typename EVAL, typename PURE>
void select_eval_for_purging(const MAT &Hobl, EVAL &eval, MAT &YB,
- PURE &pure, idgmres_state &st) {
+ PURE &pure, idgmres_state &st) {
typedef typename linalg_traits<MAT>::value_type T;
typedef typename number_traits<T>::magnitude_type R;
size_type m = st.m;
// Computation of the Ritz eigenpairs.
-
+
col_matrix< std::vector<T> > evect(st.tb_deftot, st.tb_deftot);
-
+
sub_interval SUB1(0, st.tb_deftot);
implicit_qr_algorithm(sub_matrix(Hobl, SUB1),
- sub_vector(eval, SUB1), evect);
+ sub_vector(eval, SUB1), evect);
std::fill(eval.begin() + st.tb_deftot, eval.end(), std::complex<R>(0));
-
+
std::vector< std::pair<T, size_type> > eval_sort(m);
for (size_type l = 0; l < m; ++l)
eval_sort[l] = std::pair<T, size_type>(eval[l], l);
@@ -767,37 +754,33 @@ namespace gmm {
std::vector<bool> sorted(m);
std::fill(sorted.begin(), sorted.end(), false);
-
+
std::vector<size_type> ind(m);
for (size_type l = 0; l < m; ++l) ind[l] = eval_sort[l].second;
-
+
std::vector<bool> kept(m, false);
std::fill(kept.begin(), kept.begin()+st.tb_def, true);
apply_permutation(eval, ind);
apply_permutation(evect, ind);
-
+
size_type j;
for (j = 0; j < st.tb_deftot; ++j) {
-
+
for (size_type l = 0, l < st.tb_deftot; ++l)
- YB(l, j) = std::real(evect(l, j));
-
+ YB(l, j) = std::real(evect(l, j));
+
if (std::imag(eval[j]) != R(0)) {
- for (size_type l = 0, l < m-st.tb_def; ++l)
- YB(l, j+1) = std::imag(evect(l, j));
- pure[j] = 1;
- pure[j+1] = 2;
-
- j += 2;
+ for (size_type l = 0, l < m-st.tb_def; ++l)
+ YB(l, j+1) = std::imag(evect(l, j));
+ pure[j] = 1;
+ pure[j+1] = 2;
+
+ j += 2;
}
else ++j;
}
}
-
-
-
-
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Whitespace fixes,
Konstantinos Poulios <=