From edab03a8e0b519ea41a7ec91e200ac2d9456c4a9 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Wed, 14 Dec 2011 16:09:08 +0100
Subject: [PATCH] CG method works correctly again
---
MathLib/LinAlg/Solvers/CGParallel.cpp | 42 +++++++++++++--------------
1 file changed, 20 insertions(+), 22 deletions(-)
diff --git a/MathLib/LinAlg/Solvers/CGParallel.cpp b/MathLib/LinAlg/Solvers/CGParallel.cpp
index 742fa96b0e1..c752961e30d 100644
--- a/MathLib/LinAlg/Solvers/CGParallel.cpp
+++ b/MathLib/LinAlg/Solvers/CGParallel.cpp
@@ -80,34 +80,32 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const
}
mat->precondApply(rhat);
- #pragma omp parallel
- {
- // rho = r * r^;
- rho = scpr(r, rhat, N);
-
- if (l > 1) {
- double beta = rho / rho1;
- // p = r^ + beta * p
- #pragma omp parallel for
- for (k = 0; k < N; k++) {
- p[k] = rhat[k] + beta * p[k];
- }
- } else {
+ // rho = r * r^;
+ rho = scpr(r, rhat, N);
+
+ if (l > 1) {
+ double beta = rho / rho1;
+ // p = r^ + beta * p
+ #pragma omp parallel for
+ for (k = 0; k < N; k++) {
+ p[k] = rhat[k] + beta * p[k];
+ }
+ } else {
// blas::copy(N, rhat, p);
- #pragma omp parallel for
- for (k = 0; k < N; k++) {
- p[k] = rhat[k];
- }
+ #pragma omp parallel for
+ for (k = 0; k < N; k++) {
+ p[k] = rhat[k];
}
}
// q = Ap
mat->amux(D_ONE, p, q);
+
+ // alpha = rho / p*q
+ double alpha = rho / scpr(p, q, N);
+
#pragma omp parallel
{
- // alpha = rho / p*q
- double alpha = rho / scpr(p, q, N);
-
// x += alpha * p
#pragma omp for nowait
for (k = 0; k < N; k++) {
@@ -121,10 +119,10 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const
}
#pragma omp barrier
-
- resid = sqrt(scpr(r, r, N));
} // end #pragma omp parallel
+ resid = sqrt(scpr(r, r, N));
+
if (resid <= eps * nrmb) {
eps = resid / nrmb;
nsteps = l;
--
GitLab