diff --git a/MathLib/LinAlg/Solvers/CGParallel.cpp b/MathLib/LinAlg/Solvers/CGParallel.cpp
index 77210867b02f7bf06b321d1723f7a633e2cd6033..acd28b0a5d2813440cc2249cadfbfc02515bf893 100644
--- a/MathLib/LinAlg/Solvers/CGParallel.cpp
+++ b/MathLib/LinAlg/Solvers/CGParallel.cpp
@@ -65,18 +65,23 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const
 		return 0;
 	}
 
+	unsigned k;
 	for (unsigned l = 1; l <= nsteps; ++l) {
 #ifndef NDEBUG
 		std::cout << "Step " << l << ", resid=" << resid / nrmb << std::endl;
 #endif
 		// r^ = C r
-		blas::copy(N, r, rhat);
+		// 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, num_threads);
 
-		unsigned k;
 		if (l > 1) {
 			double beta = rho / rho1;
 			// p = r^ + beta * p
@@ -84,10 +89,20 @@ unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const
 			for (k = 0; k < N; k++) {
 				p[k] = rhat[k] + beta * p[k];
 			}
-		} else blas::copy(N, rhat, p);
+		} else {
+//			blas::copy(N, rhat, p);
+			#pragma omp parallel for
+			for (k = 0; k < N; k++) {
+				p[k] = rhat[k];
+			}
+		}
 
 		// q = Ap
-		blas::setzero(N, q);
+//		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
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 06b01cd6d0a88e248b92f4d078f4d7a7f84c418b..ccad16867335050e3840b2711fa74c7fd07d7b33 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -15,7 +15,6 @@ double scpr(double const * const v, double const * const w, unsigned n,
 {
         long double res (0.0);
         unsigned k;
-        omp_set_num_threads (num_of_threads);
         #pragma omp parallel
         {
                 #pragma omp parallel for reduction (+:res)