diff --git a/MathLib/LinAlg/Solvers/CGParallel.cpp b/MathLib/LinAlg/Solvers/CGParallel.cpp index 79f1cf26e710ecaef8c4ff6524b7d0b0fc1046a5..fa2b8153c0725556789d5aefd73f2565c6c15b4a 100644 --- a/MathLib/LinAlg/Solvers/CGParallel.cpp +++ b/MathLib/LinAlg/Solvers/CGParallel.cpp @@ -84,7 +84,7 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const mat->precondApply(rhat); // rho = r * r^; - rho = scpr(r, rhat, N, num_threads); + rho = scpr(r, rhat, N); if (l > 1) { double beta = rho / rho1; @@ -102,29 +102,29 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const } // q = Ap -// blas::setzero(N, q); - #pragma omp parallel for - for (k = 0; k < N; k++) { - q[k] = 0.0; - } mat->amux(D_ONE, p, q); // alpha = rho / p*q - double alpha = rho / scpr(p, q, N, num_threads); + double alpha = rho / scpr(p, q, N); - // x += alpha * p - #pragma omp parallel for - for (k = 0; k < N; k++) { - x[k] += alpha * p[k]; - } + #pragma omp parallel + { + // x += alpha * p + #pragma omp for nowait + for (k = 0; k < N; k++) { + x[k] += alpha * p[k]; + } - // r -= alpha * q - #pragma omp parallel for - for (k = 0; k < N; k++) { - r[k] -= alpha * q[k]; + // r -= alpha * q + #pragma omp for nowait + for (k = 0; k < N; k++) { + r[k] -= alpha * q[k]; + } + + #pragma omp barrier } - resid = sqrt(scpr(r, r, N, num_threads)); + resid = sqrt(scpr(r, r, N)); if (resid <= eps * nrmb) { eps = resid / nrmb; diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp index ccad16867335050e3840b2711fa74c7fd07d7b33..35bd2612dc18b2b2ae3e1384638e4bd7f989443a 100644 --- a/MathLib/MathTools.cpp +++ b/MathLib/MathTools.cpp @@ -10,15 +10,14 @@ namespace MathLib { #ifdef _OPENMP -double scpr(double const * const v, double const * const w, unsigned n, - unsigned num_of_threads) +double scpr(double const * const v, double const * const w, unsigned n) { - long double res (0.0); + long double res (v[0]*w[0]); unsigned k; #pragma omp parallel { #pragma omp parallel for reduction (+:res) - for (k = 0; k<n; k++) { + for (k = 1; k<n; k++) { res += v[k] * w[k]; } } diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h index 1d88785af6d6bbe991890b0e1349f0759d996870..4a9d2b804f9272a0c531be13973088508f023d34 100644 --- a/MathLib/MathTools.h +++ b/MathLib/MathTools.h @@ -35,8 +35,7 @@ double scpr(const T* v0, const T* v1, size_t n) } #ifdef _OPENMP -double scpr(double const * const v, double const * const w, unsigned n, - unsigned num_of_threads); +double scpr(double const * const v, double const * const w, unsigned n); #endif /**