From 63a72fb392edee0e7ee7c88d1dc23e654571cf68 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 12 Sep 2022 16:51:40 +0200 Subject: [PATCH] Revise vector operations in SolverIDR to improve performance --- include/deal.II/lac/solver_idr.h | 91 +++++++++++++++++--------------- 1 file changed, 49 insertions(+), 42 deletions(-) diff --git a/include/deal.II/lac/solver_idr.h b/include/deal.II/lac/solver_idr.h index 1744585e8c..57405e21f4 100644 --- a/include/deal.II/lac/solver_idr.h +++ b/include/deal.II/lac/solver_idr.h @@ -222,7 +222,7 @@ namespace internal if (data[i] == nullptr) { data[i] = std::move(typename VectorMemory::Pointer(mem)); - data[i]->reinit(temp); + data[i]->reinit(temp, true); } return *data[i]; } @@ -324,21 +324,15 @@ SolverIDR::solve(const MatrixType & A, // Define temporary vectors which do not do not depend on s typename VectorMemory::Pointer r_pointer(this->memory); typename VectorMemory::Pointer v_pointer(this->memory); - typename VectorMemory::Pointer vhat_pointer(this->memory); typename VectorMemory::Pointer uhat_pointer(this->memory); - typename VectorMemory::Pointer ghat_pointer(this->memory); VectorType &r = *r_pointer; VectorType &v = *v_pointer; - VectorType &vhat = *vhat_pointer; VectorType &uhat = *uhat_pointer; - VectorType &ghat = *ghat_pointer; r.reinit(x, true); v.reinit(x, true); - vhat.reinit(x, true); uhat.reinit(x, true); - ghat.reinit(x, true); // Initial residual A.vmult(r, x); @@ -362,10 +356,9 @@ SolverIDR::solve(const MatrixType & A, std::normal_distribution<> normal_distribution(0.0, 1.0); for (unsigned int i = 0; i < s; ++i) { - VectorType &tmp_g = G(i, x); - VectorType &tmp_u = U(i, x); - tmp_g = 0; - tmp_u = 0; + // Initialize vectors + G(i, x); + U(i, x); // Compute random set of s orthonormalized vectors Q // Note: the first vector is chosen to be the initial @@ -436,43 +429,58 @@ SolverIDR::solve(const MatrixType & A, Mk_inv.vmult(gamma, phik); } - { - v = r; + v = r; - unsigned int j = 0; - for (unsigned int i = k; i < s; ++i, ++j) - v.add(-1.0 * gamma(j), G[i]); - preconditioner.vmult(vhat, v); + if (step > 1) + { + for (unsigned int i = k, j = 0; i < s; ++i, ++j) + v.add(-gamma(j), G[i]); + } + + preconditioner.vmult(uhat, v); - uhat = vhat; + if (step > 1) + { + uhat.sadd(omega, gamma(0), U[k]); + for (unsigned int i = k + 1, j = 1; i < s; ++i, ++j) + uhat.add(gamma(j), U[i]); + } + else uhat *= omega; - j = 0; - for (unsigned int i = k; i < s; ++i, ++j) - uhat.add(gamma(j), U[i]); - A.vmult(ghat, uhat); - } + + A.vmult(G[k], uhat); // Update G and U - // Orthogonalize ghat to Q0,..,Q_{k-1} - // and update uhat - for (unsigned int i = 0; i < k; ++i) + // Orthogonalize G[k] to Q0,..,Q_{k-1} and update uhat + if (k > 0) { - double alpha = (Q[i] * ghat) / M(i, i); - ghat.add(-alpha, G[i]); - uhat.add(-alpha, U[i]); + double alpha = Q[0] * G[k] / M(0, 0); + for (unsigned int i = 1; i < k; ++i) + { + const double 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 + if (i % 2 == 1) + uhat.add(-alpha_old, U[i - 1], -alpha, U[i]); + } + M(k, k) = G[k].add_and_dot(-alpha, G[k - 1], Q[k]); + if (k % 2 == 1) + uhat.add(-alpha, U[k - 1]); } - G[k] = ghat; - U[k] = uhat; + else + M(k, k) = G[k] * Q[k]; + + U[k].swap(uhat); // Update kth column of M - for (unsigned int i = k; i < s; ++i) + for (unsigned int i = k + 1; i < s; ++i) M(i, k) = Q[i] * G[k]; - // Orthogonalize r to Q0,...,Qk, - // update x + // Orthogonalize r to Q0,...,Qk, update x { - double beta = phi(k) / M(k, k); - r.add(-1.0 * beta, G[k]); + const double beta = phi(k) / M(k, k); + r.add(-beta, G[k]); x.add(beta, U[k]); print_vectors(step, x, r, U[k]); @@ -502,18 +510,17 @@ SolverIDR::solve(const MatrixType & A, break; // Update r and x - preconditioner.vmult(vhat, r); - A.vmult(v, vhat); + preconditioner.vmult(uhat, r); + A.vmult(v, uhat); omega = (v * r) / (v * v); - r.add(-1.0 * omega, v); - x.add(omega, vhat); + res = std::sqrt(r.add_and_dot(-1.0 * omega, v, r)); + x.add(omega, uhat); - print_vectors(step, x, r, vhat); + print_vectors(step, x, r, uhat); // Check for convergence - res = r.l2_norm(); iteration_state = this->iteration_status(step, res, x); if (iteration_state != SolverControl::iterate) break; -- 2.39.5