From: Martin Kronbichler Date: Tue, 13 Sep 2022 09:57:15 +0000 (+0200) Subject: Switch number types in SolverIDR to VectorType::value_type X-Git-Tag: v9.5.0-rc1~969^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0fc6a4d4fd33d0f02436140ec9c68ba483165299;p=dealii.git Switch number types in SolverIDR to VectorType::value_type --- diff --git a/include/deal.II/lac/solver_idr.h b/include/deal.II/lac/solver_idr.h index 57405e21f4..d3cae4e16a 100644 --- a/include/deal.II/lac/solver_idr.h +++ b/include/deal.II/lac/solver_idr.h @@ -338,8 +338,11 @@ SolverIDR::solve(const MatrixType & A, A.vmult(r, x); r.sadd(-1.0, 1.0, b); + using value_type = typename VectorType::value_type; + using real_type = typename numbers::NumberTraits::real_type; + // Check for convergent initial guess - double res = r.l2_norm(); + real_type res = r.l2_norm(); iteration_state = this->iteration_status(step, res, x); if (iteration_state == SolverControl::success) return; @@ -348,7 +351,7 @@ SolverIDR::solve(const MatrixType & A, internal::SolverIDRImplementation::TmpVectors G(s, this->memory); internal::SolverIDRImplementation::TmpVectors U(s, this->memory); internal::SolverIDRImplementation::TmpVectors Q(s, this->memory); - FullMatrix M(s, s); + FullMatrix M(s, s); // Random number generator for vector entries of // Q (normal distribution, mean=0 sigma=1) @@ -393,7 +396,7 @@ SolverIDR::solve(const MatrixType & A, M(i, i) = 1.; } - double omega = 1.; + value_type omega = 1.; bool early_exit = false; @@ -403,7 +406,7 @@ SolverIDR::solve(const MatrixType & A, ++step; // Compute phi - Vector phi(s); + Vector phi(s); for (unsigned int i = 0; i < s; ++i) phi(i) = Q[i] * r; @@ -411,10 +414,10 @@ SolverIDR::solve(const MatrixType & A, for (unsigned int k = 0; k < s; ++k) { // Solve M(k:s)*gamma = phi(k:s) - Vector gamma(s - k); + Vector gamma(s - k); { - Vector phik(s - k); - FullMatrix Mk(s - k, s - k); + Vector phik(s - k); + FullMatrix Mk(s - k, s - k); std::vector indices; unsigned int j = 0; for (unsigned int i = k; i < s; ++i, ++j) @@ -424,7 +427,7 @@ SolverIDR::solve(const MatrixType & A, } Mk.extract_submatrix_from(M, indices, indices); - FullMatrix Mk_inv(s - k, s - k); + FullMatrix Mk_inv(s - k, s - k); Mk_inv.invert(Mk); Mk_inv.vmult(gamma, phik); } @@ -454,10 +457,10 @@ SolverIDR::solve(const MatrixType & A, // Orthogonalize G[k] to Q0,..,Q_{k-1} and update uhat if (k > 0) { - double alpha = Q[0] * G[k] / M(0, 0); + value_type alpha = Q[0] * G[k] / M(0, 0); for (unsigned int i = 1; i < k; ++i) { - const double alpha_old = alpha; + const value_type alpha_old = alpha; alpha = G[k].add_and_dot(-alpha, G[i - 1], Q[i]) / M(i, i); // update uhat every other iteration to reduce vector access @@ -479,7 +482,7 @@ SolverIDR::solve(const MatrixType & A, // Orthogonalize r to Q0,...,Qk, update x { - const double beta = phi(k) / M(k, k); + const value_type beta = phi(k) / M(k, k); r.add(-beta, G[k]); x.add(beta, U[k]);