diff --git a/MathLib/LinAlg/Solvers/CGParallel.cpp b/MathLib/LinAlg/Solvers/CGParallel.cpp
index 8c74cb401b417a826497ae0c9ae0d4ab25cea24d..c752961e30d9b2aaebc671dfe508f6cb0aebdbcb 100644
--- a/MathLib/LinAlg/Solvers/CGParallel.cpp
+++ b/MathLib/LinAlg/Solvers/CGParallel.cpp
@@ -71,41 +71,41 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const
 		std::cout << "Step " << l << ", resid=" << resid / nrmb << std::endl;
 #endif
 
-		#pragma omp parallel
-		{
-			// r^ = C r
-			// rhat = r
-//			blas::copy(N, r, rhat);
+		// r^ = C r
+		// rhat = r
+//		blas::copy(N, r, rhat);
+		#pragma omp parallel for
+		for (k = 0; k < N; k++) {
+			rhat[k] = r[k];
+		}
+		mat->precondApply(rhat);
+
+		// 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++) {
-				rhat[k] = r[k];
+				p[k] = rhat[k] + beta * p[k];
 			}
-			mat->precondApply(rhat);
-
-			// 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 {
+		} 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);
+		// q = Ap
+		mat->amux(D_ONE, p, q);
 
-			// alpha = rho / p*q
-			double alpha = rho / scpr(p, q, N);
+		// alpha = rho / p*q
+		double alpha = rho / scpr(p, q, N);
 
+		#pragma omp parallel
+		{
 			// x += alpha * p
 			#pragma omp for nowait
 			for (k = 0; k < N; k++) {
@@ -119,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;