Commit for #6

Added early exit for internal BiCGStab(l). Added test example based on
#6. Resolved #6 for internal solvers.

.gitignore

@@ -5,4 +5,3 @@ compile64_vs2013/
 examples/Data/
 examples/Matrices/
 examples/Grids/
-tests/

CMakeLists.txt

@@ -226,6 +226,8 @@ add_subdirectory(projects)
 endif(COMPILE_PROJECTS)
 
 if(COMPILE_TESTS)
+enable_testing()
+include(CTest)
 add_subdirectory(tests)
 endif(COMPILE_TESTS)
 
@@ -266,9 +268,9 @@ install(FILES
 
 install(TARGETS inmost EXPORT inmost-targets
         ARCHIVE DESTINATION lib
-		LIBRARY DESTINATION lib
-		RUNTIME DESTINATION bin
-		PUBLIC_HEADER DESTINATION include)
+	LIBRARY DESTINATION lib
+	RUNTIME DESTINATION bin
+	PUBLIC_HEADER DESTINATION include)
 
 
 install(EXPORT inmost-targets DESTINATION cmake)

Tests/CMakeLists.txt

+add_subdirectory(solver_test000)
\ No newline at end of file

Tests/solver_test000/CMakeLists.txt

+project(solver_test000)
+set(SOURCE main.cpp)
+
+add_executable(solver_test000 ${SOURCE})
+target_link_libraries(solver_test000 inmost)
+
+if(USE_MPI)
+  message("linking solver_test000 with MPI")
+  target_link_libraries(solver_test000 ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(solver_test000 PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+
+
+if(USE_SOLVER_ANI)
+  message("linking solver_test000 with ani3d")
+  target_link_libraries(solver_test000 ani3d)
+endif()
+if(USE_SOLVER_PETSC)
+  message("linking solver_test000 with PETSc")
+  target_link_libraries(solver_test000 ${PETSC_LIBRARIES})
+endif()
+if(USE_SOLVER_TRILINOS)
+  message("linking solver_test000 with Trilinos")
+  target_link_libraries(solver_test000 ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
+endif()
+
+
+add_test(NAME solver_test000_serial_inner_ilu2                 COMMAND $<TARGET_FILE:solver_test000> 0 0)
+add_test(NAME solver_test000_serial_inner_mliluc               COMMAND $<TARGET_FILE:solver_test000> 0 1)
+if(USE_SOLVER_PETSC)
+add_test(NAME solver_test000_serial_petsc                      COMMAND $<TARGET_FILE:solver_test000> 0 2)
+endif()
+
+if(USE_SOLVER_TIRLINOS)
+add_test(NAME solver_test000_serial_trilinos_aztec             COMMAND $<TARGET_FILE:solver_test000> 0 3)
+add_test(NAME solver_test000_serial_trilinos_ifpack            COMMAND $<TARGET_FILE:solver_test000> 0 4)
+add_test(NAME solver_test000_serial_trilinos_ml                COMMAND $<TARGET_FILE:solver_test000> 0 5)
+add_test(NAME solver_test000_serial_trilinos_belos             COMMAND $<TARGET_FILE:solver_test000> 0 6)
+endif()
+
+if(USE_SOLVER_ANI)
+add_test(NAME solver_test000_serial_ani COMMAND $<TARGET_FILE:solver_test000> 2 0 7)
+endif()
+
+if( USE_MPI )
+if( EXISTS ${MPIEXEC} )
+add_test(NAME solver_test000_parallel_normal_inner_ilu2        COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 0)
+add_test(NAME solver_test000_parallel_permute1_inner_ilu2      COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 0)
+add_test(NAME solver_test000_parallel_permute2_inner_ilu2      COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 0)
+
+add_test(NAME solver_test000_parallel_normal_inner_mliluc      COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 1)
+add_test(NAME solver_test000_parallel_permute1_inner_mliluc    COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 1)
+add_test(NAME solver_test000_parallel_permute2_inner_mliluc    COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 1)
+
+if(USE_SOLVER_PETSC)
+add_test(NAME solver_test000_parallel_normal_petsc             COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 2)
+#add_test(NAME solver_test000_parallel_permute1_petsc           COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 2)
+#add_test(NAME solver_test000_parallel_permute2_petsc           COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 2)
+endif()
+
+if(USE_SOLVER_TRILINOS)
+add_test(NAME solver_test000_parallel_normal_trilinos_aztec    COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 3)
+add_test(NAME solver_test000_parallel_permute1_trilinos_aztec  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 3)
+add_test(NAME solver_test000_parallel_permute2_trilinos_aztec  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 3)
+
+add_test(NAME solver_test000_parallel_normal_trilinos_ifpack   COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 4)
+add_test(NAME solver_test000_parallel_permute1_trilinos_ifpack COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 4)
+add_test(NAME solver_test000_parallel_permute2_trilinos_ifpack COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 4)
+
+add_test(NAME solver_test000_parallel_normal_trilinos_ml       COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 5)
+add_test(NAME solver_test000_parallel_permute1_trilinos_ml     COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 5)
+add_test(NAME solver_test000_parallel_permute2_trilinos_ml     COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 5)
+
+#add_test(NAME solver_test000_parallel_normal_trilinos_belos    COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 0 6)
+#add_test(NAME solver_test000_parallel_permute1_trilinos_belos  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 1 6)
+#add_test(NAME solver_test000_parallel_permute2_trilinos_belos  COMMAND ${MPIEXEC} -np 4 $<TARGET_FILE:solver_test000> 2 6)
+endif()
+endif()
+endif()

Tests/solver_test000/main.cpp

+#include <cstdio>
+#include <cmath>
+
+#include "../../inmost.h"
+using namespace INMOST;
+
+int main(int argc,char ** argv)
+{
+	int permut = 0;
+	int solver = 0; // 0 - INNER_ILU2, 1 - INNER_MLILUC, 2 - PETSc
+	// 3 - Trilinos_Aztec, 4 - Trilinos_Ifpack,
+	// 5 - Trilinos_ML, 6 - Trilinos_Belos, 7 - ANI
+	int rank,procs,newrank;
+
+	Solver::Initialize(&argc,&argv,"database.txt"); // Initialize the solver and MPI activity
+#if defined(USE_MPI)
+	MPI_Comm_rank(MPI_COMM_WORLD,&rank);  // Get the rank of the current process
+	MPI_Comm_size(MPI_COMM_WORLD,&procs); // Get the total number of processors used
+#else
+	rank = 0;
+	procs = 1;
+#endif
+
+	if (argc > 1)  permut = atoi(argv[1]);
+	if (argc > 2)  solver = atoi(argv[2]);
+
+	if (permut < procs)  newrank = (rank + permut) % procs;
+	else newrank = (permut - rank) % procs;
+
+	std::cout << rank <<  " -> " << newrank << std::endl;
+
+	Solver::Type type;
+	switch(solver)
+	{
+	case 0: type = Solver::INNER_ILU2; break;
+	case 1: type = Solver::INNER_MLILUC; break;
+	case 2: type = Solver::PETSc; break;
+	case 3: type = Solver::Trilinos_Aztec; break;
+	case 4: type = Solver::Trilinos_Ifpack; break;
+	case 5: type = Solver::Trilinos_ML; break;
+	case 6: type = Solver::Trilinos_Belos; break;
+	case 7: type = Solver::ANI; break;
+	}
+
+	{
+		Solver S(type); // Specify the linear solver
+
+
+		Solver::Matrix A; // Declare the matrix of the linear system to be solved
+		Solver::Vector x,b; // Declare the solution and the right-hand side vectors
+
+
+		INMOST_DATA_ENUM_TYPE mbeg, mend;
+
+		mbeg = newrank * 10;
+		mend = (newrank+1) * 10;
+		A.SetInterval(mbeg,mend);
+		x.SetInterval(mbeg,mend);
+		b.SetInterval(mbeg,mend);
+
+		for( INMOST_DATA_ENUM_TYPE i = mbeg; i !=  mend; i++ )
+		{
+			A[i][i] = 10.0;
+			b[i] = i;
+			x[i] = 0.0;
+		}
+		if (rank==0)  std::cout << "next call S.SetMatrix(A);" << std::endl;
+		S.SetMatrix(A); // Compute the preconditioner for the original matrix
+		if (rank==0)  std::cout << "next call S.Solve(b,x);" << std::endl;
+		if( !S.Solve(b,x) )   // Solve the linear system with the previously computted preconditioner
+		{
+			if( rank == 0 )
+				std::cout << S.GetReason() << std::endl;
+			MPI_Abort(MPI_COMM_WORLD,-1);
+
+		}
+		
+
+		double err = 0;
+		for( INMOST_DATA_ENUM_TYPE i = mbeg; i !=  mend; i++ )
+		{
+			err += (x[i] - static_cast<double>(i)/10.0);
+		}
+#if defined(USE_MPI)
+		double tmp;
+		MPI_Allreduce(&err,&tmp,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+		err = tmp;
+#endif
+
+		if (rank==0)  
+		{
+			std::cout << "done, error " << err << "\t\t\t" << std::endl;
+		}
+
+		if( fabs(err) > 1.0e-5 || err!=err) MPI_Abort(MPI_COMM_WORLD,-1);
+
+		x.Save("output.rhs");
+		b.Save("b.rhs");
+		A.Save("A.mtx");
+	}
+
+	Solver::Finalize(); // Finalize solver and close MPI activity
+	return 0;
+}
\ No newline at end of file

examples/Solver/main.cpp

@@ -151,6 +151,7 @@ int main(int argc,char ** argv)
 #if defined(USE_PARTITIONER)
 	Partitioner::Initialize(&argc,&argv);
 #endif
+	/*
 	{
 		Mesh m;
 		//m.SetFileOption("VERBOSITY","2");
@@ -161,7 +162,7 @@ int main(int argc,char ** argv)
 	Partitioner::Finalize();
 	Mesh::Finalize();
 	return 0;
-
+	*/
 
 	if( argc > 1 )
 	{
@@ -520,5 +521,6 @@ int main(int argc,char ** argv)
 	Partitioner::Finalize();
 #endif
 	Solver::Finalize();
+	Mesh::Finalize();
 	return 0;
 }

inmost_solver.h

@@ -94,9 +94,9 @@ namespace INMOST
 			/// Restore initial nonparallel state of the Matrix with no overlap.
 			void RestoreMatrix(Matrix & m);
 			/// Prepare parallel state of the Vector.
-			void PrepareVector(Vector & v);
+			void PrepareVector(Vector & v) const;
 			/// Restore initial nonparallel state of the Vector.
-			void RestoreVector(Vector & v); //Restore initial nonparallel state of the vector
+			void RestoreVector(Vector & v) const;
 			/// Retrieve the processor number by binary search for the specified global index.
 			INMOST_DATA_ENUM_TYPE GetProcessor(INMOST_DATA_ENUM_TYPE gind) const; //retrive processor by binary search in global_to_proc
 			void      GetOverlapRegion(INMOST_DATA_ENUM_TYPE proc, INMOST_DATA_ENUM_TYPE & mbeg, INMOST_DATA_ENUM_TYPE & mend) const;
@@ -113,9 +113,9 @@ namespace INMOST
 			/// Sum shared values in parallel vector.
 			void                  Accumulate(Vector & x); // sum shared values in parallel vector
 			/// Get the sum of num elements of real array on all processes.
-			void                  Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num);
+			void                  Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num) const;
 			/// Get the communicator which the solver is associated with.
-			INMOST_MPI_Comm         GetComm() {return comm;}
+			INMOST_MPI_Comm         GetComm() const {return comm;}
 			// Access to arrays below allows to organize manual exchange
 			INMOST_MPI_Request    * GetSendRequests() {assert(!send_requests.empty()); return &send_requests[0];}
 			INMOST_MPI_Request    * GetRecvRequests() {assert(!recv_requests.empty()); return &recv_requests[0];}
@@ -128,6 +128,9 @@ namespace INMOST
 			//for debug
 			//~ void                 BeginSequentialCode() {for(int i = 0; i < rank; i++) MPI_Barrier(comm);}
 			//~ void                   EndSequentialCode() {for(int i = rank; i < size; i++) MPI_Barrier(comm);}
+
+			/// Get the scalar product of the specified interval of the distributed vector.
+			INMOST_DATA_REAL_TYPE ScalarProd(Vector const & left, Vector const & right, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const;
 		};
 
 		/// Distributed vector class.
@@ -176,8 +179,7 @@ namespace INMOST
 			/// Get the communicator which the vector is associated with.
 			INMOST_MPI_Comm        GetCommunicator() const {return comm;}
 
-			/// Get the scalar product of the specified interval of the distributed vector.
-			INMOST_DATA_REAL_TYPE ScalarProd(Vector const & other, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const;
+			
 
 
 			/// Save the distributed vector to a single data file using parallel MPI I/O.

solver.cpp

@@ -77,14 +77,14 @@ namespace INMOST
 		}
 		return size;
 	}
-	void Solver::OrderInfo::Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num)
+	void Solver::OrderInfo::Integrate(INMOST_DATA_REAL_TYPE * inout, INMOST_DATA_ENUM_TYPE num) const
 	{
 #if defined(USE_MPI)
+		if( GetSize() == 1 ) return;
 		int ierr = 0;
-		INMOST_DATA_REAL_TYPE * temp = new INMOST_DATA_REAL_TYPE [num];
-		memcpy(temp,inout,sizeof(INMOST_DATA_REAL_TYPE)*num);
-		GUARD_MPI(MPI_Allreduce(temp,inout,num,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
-		delete [] temp;
+		dynarray<INMOST_DATA_REAL_TYPE,64> temp(num);
+		memcpy(temp.data(),inout,sizeof(INMOST_DATA_REAL_TYPE)*num);
+		GUARD_MPI(MPI_Allreduce(temp.data(),inout,num,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
 #else
 		(void) inout;
 		(void) num;
@@ -147,7 +147,7 @@ namespace INMOST
 			{
 				storage_type temp(size*2);
 				//assemble array that includes rank
-				for(int k = 0; k < size; k++)
+				for(int k = 0; k < size; ++k)
 				{
 					temp[2*k+0] = global_overlap[2*k];
 					temp[2*k+1] = k;
@@ -158,6 +158,8 @@ namespace INMOST
 				MPI_Group oldg, newg;
 				MPI_Comm newcomm;
 				std::vector<int> ranks(size);
+				for(int k = 0; k < size; ++k)
+					ranks[k] = temp[2*k+1];
 				GUARD_MPI(MPI_Comm_group(comm,&oldg));
 				GUARD_MPI(MPI_Group_incl(oldg,size,&ranks[0],&newg));
 				GUARD_MPI(MPI_Comm_create(comm,newg,&newcomm));
@@ -553,14 +555,14 @@ namespace INMOST
 		have_matrix = false;
 	}
 	
-	void Solver::OrderInfo::PrepareVector(Vector & v)
+	void Solver::OrderInfo::PrepareVector(Vector & v) const
 	{
 		if( !have_matrix ) throw PrepareMatrixFirst;
 		v.SetInterval(local_vector_begin,local_vector_end);
 		v.isParallel() = true;
 	}
 	
-	void Solver::OrderInfo::RestoreVector(Vector & v)
+	void Solver::OrderInfo::RestoreVector(Vector & v) const
 	{
 		assert(have_matrix);
 		if( v.isParallel() )
@@ -703,6 +705,7 @@ namespace INMOST
 	{
 		//std::cout << __FUNCTION__ << " start" << std::endl;
 #if defined(USE_MPI)
+		if( GetSize() == 1 ) return;
 		//use MPI_Put/MPI_Get to update vector
 		assert(x.isParallel()); //the vector was prepared
 		INMOST_DATA_ENUM_TYPE i, j = 1, k, l = 0;
@@ -749,6 +752,7 @@ namespace INMOST
 	{
 		//std::cout << __FUNCTION__ << " start" << std::endl;
 #if defined(USE_MPI)
+		if( GetSize() == 1 ) return;
 		//use MPI_Put/MPI_Get to update vector
 		assert(x.isParallel()); //the vector was prepared
 		INMOST_DATA_ENUM_TYPE i, j = 1, k, l = 0;
@@ -794,17 +798,12 @@ namespace INMOST
 		//std::cout << __FUNCTION__ << " end" << std::endl;
 	}
 
-	INMOST_DATA_REAL_TYPE Solver::Vector::ScalarProd(Vector const & other, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const
+	INMOST_DATA_REAL_TYPE Solver::OrderInfo::ScalarProd(Vector const & left, Vector const & right, INMOST_DATA_ENUM_TYPE index_begin, INMOST_DATA_ENUM_TYPE index_end) const
 	{
 		INMOST_DATA_REAL_TYPE ret = 0;
-		for(INMOST_DATA_ENUM_TYPE i = index_begin; i < index_end; i++)
-			ret += data[i]*other[i];
-#if defined(USE_MPI)
-		INMOST_DATA_REAL_TYPE temp;
-		int ierr = 0;
-		GUARD_MPI(MPI_Allreduce(&ret,&temp,1,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm));
-		ret = temp;
-#endif
+		for(INMOST_DATA_ENUM_TYPE i = index_begin; i < index_end; ++i)
+			ret += left[i]*right[i];
+		Integrate(&ret,1);
 		return ret;
 	}
 
@@ -1833,7 +1832,8 @@ namespace INMOST
 			}
 			else sol->EnumParameter(":fill") = static_cast<INMOST_DATA_ENUM_TYPE>(preconditioner_fill_level);
 
-			sol->EnumParameter("levels") = solver_gmres_substeps;
+			if( sizeof(KSOLVER) == sizeof(BCGSL_solver) )
+				sol->EnumParameter("levels") = solver_gmres_substeps;
 
 			if (!sol->isInitialized()) 
 			{
@@ -2158,9 +2158,9 @@ namespace INMOST
 			{
 				List->set("Num Blocks",100);
 				List->set("Block Size",1);
-				List->set("Maximum Iterations",maximum_iterations);
+				List->set("Maximum Iterations",static_cast<int>(maximum_iterations));
 				List->set("Maximum Restarts",20);
-				List->set("Convergence Tolerance",relative_tolerance);
+				List->set("Convergence Tolerance",static_cast<double>(relative_tolerance));
 			}
 			//int verbosity = Belos::Warnings + Belos::Errors + Belos::StatusTestDetails + Belos::Debug + Belos::FinalSummary + Belos::TimingDetails;
 			//List->set("Verbosity",verbosity);
@@ -2187,8 +2187,8 @@ namespace INMOST
 			
 			Teuchos::RCP< Belos::SolverManager<double,Epetra_MultiVector,Epetra_Operator> > BelosSolver =
 			//	Teuchos::rcp(new Belos::BlockCgSolMgr<double,Epetra_MultiVector,Epetra_Operator >);
-			//	Teuchos::rcp(new Belos::BlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator >);
-				Teuchos::rcp(new Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator >);
+				Teuchos::rcp(new Belos::BlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator >);
+			//	Teuchos::rcp(new Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator >);
 			BelosSolver->setParameters(List);
 			BelosSolver->setProblem(Problem);
 			Belos::ReturnType ret = BelosSolver->solve();

solver_bcgsl.hpp

@@ -20,7 +20,7 @@ namespace INMOST
 		INMOST_DATA_ENUM_TYPE iters, maxits, l, last_it;
 		INMOST_DATA_REAL_TYPE resid;
 		INMOST_DATA_REAL_TYPE * tau, * sigma, * gamma, *theta1, * theta2, * theta3;
-		Solver::Vector r0, t, * u, * r;
+		Solver::Vector r_tilde, x0, t, * u, * r;
 		Solver::Matrix * Alink;
 		Method * prec;
 		std::string reason;
@@ -66,7 +66,8 @@ namespace INMOST
 		{
 			if (isInitialized()) Finalize();
 			if (prec != NULL && !prec->isInitialized()) prec->Initialize();
-			info->PrepareVector(r0);
+			info->PrepareVector(r_tilde);
+			info->PrepareVector(x0);
 			info->PrepareVector(t);
 			tau = new INMOST_DATA_REAL_TYPE[l * 5 + l*l];
 			sigma = tau + l*l;
@@ -134,35 +135,51 @@ namespace INMOST
 			if (!isFinalized()) Finalize();
 			if (prec != NULL) delete prec;
 		}
+		void ApplyOperator(Solver::Vector & Input, Solver::Vector & Output)
+		{
+			if (prec != NULL) //right preconditioning here! for left preconditioner have to reverse order
+			{
+				prec->Solve(Input, t); 
+				info->Update(t);
+				Alink->MatVec(1.0,t,0,Output);
+				info->Update(Output);
+			}
+			else
+			{
+				Alink->MatVec(1.0,t,0,Output);
+				info->Update(Output);
+			}
+		}
 		bool Solve(Solver::Vector & RHS, Solver::Vector & SOL)
 		{
 			assert(isInitialized());
 			INMOST_DATA_ENUM_TYPE vbeg,vend, vlocbeg, vlocend;
-			INMOST_DATA_REAL_TYPE rho0 = 1, rho1, alpha = 0, beta, omega = 1;
+			INMOST_DATA_REAL_TYPE rho0 = 1, rho1, alpha = 0, beta, omega = 1, eta;
 			INMOST_DATA_REAL_TYPE resid0, resid, temp[2];
-			bool is_parallel = info->GetSize() > 1;
 			iters = 0;
 			info->PrepareVector(SOL);
 			info->PrepareVector(RHS);
-			if( is_parallel ) info->Update(SOL);
-			if( is_parallel ) info->Update(RHS);
+			info->Update(SOL);
+			info->Update(RHS);
 			if( prec != NULL ) prec->ReplaceSOL(SOL);
 			if( prec != NULL ) prec->ReplaceRHS(RHS);
 			info->GetLocalRegion(info->GetRank(),vlocbeg,vlocend);
 			info->GetVectorRegion(vbeg,vend);
+			//r[0] = b
 			std::copy(RHS.Begin(),RHS.End(),r[0].Begin());
 			{
+				// r[0] = r[0] - A x
 				Alink->MatVec(-1,SOL,1,r[0]); //global multiplication, r probably needs an update
-				if( is_parallel ) info->Update(r[0]); // r is good
-			}
-			std::copy(r[0].Begin(),r[0].End(),r0.Begin());
-			std::fill(u[0].Begin(),u[0].End(),0);
-			{
-				resid = 0;
-				for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k != vlocend; k++) resid += r[0][k]*r[0][k];
-				if( is_parallel ) info->Integrate(&resid,1);
+				info->Update(r[0]); // r is good
+				std::copy(x0.Begin(),x0.End(),SOL.Begin()); //x0 = x
+				std::fill(SOL.Begin(),SOL.End(),0.0); //x = 0
 			}
-			last_resid = resid = resid0 = sqrt(resid);
+			std::copy(r[0].Begin(),r[0].End(),r_tilde.Begin()); // r_tilde = r[0]
+			std::fill(u[0].Begin(),u[0].End(),0); // u[0] = 0
+			resid = info->ScalarProd(r[0],r[0],vlocbeg,vlocend); //resid = dot(r[0],r[0])
+			for(INMOST_DATA_ENUM_TYPE k = vbeg; k != vend; k++) // r_tilde = r[0] / dot(r[0],r[0])
+				r_tilde[k] /= resid;
+			last_resid = resid = resid0 = sqrt(resid); //resid = sqrt(dot(r[0],r[0])
 			last_it = 0;
 #if defined(REPORT_RESIDUAL)
 			if( info->GetRank() == 0 ) 
@@ -174,6 +191,13 @@ namespace INMOST
 				fflush(stdout);
 			}
 #endif
+
+			if( last_resid < atol || last_resid < rtol*resid0 ) 
+			{
+				reason = "initial solution satisfy tolerances";
+				goto exit;
+			}
+
 			long double tt, ts, tp, ttt;
 			INMOST_DATA_ENUM_TYPE i = 0;
 			while( true )
@@ -184,57 +208,49 @@ namespace INMOST
 				tt = Timer();
 				for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++)
 				{
-					rho1 = 0.0;
-					for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k < vlocend; ++k)
-						rho1 += r[j][k]*r0[k];
-					if( is_parallel ) info->Integrate(&rho1,1);
+					rho1 = info->ScalarProd(r[j],r_tilde,vlocbeg,vlocend); // rho1 = dot(r[j],r_tilde)
 					if( fabs(rho0) < 1.0e-54 ) 
 					{
 						reason = "denominator(1) is zero";
 						goto exit;
 					}
-					beta = alpha * rho1/rho0;
+					beta = alpha * (rho1/rho0);
 					rho0 = rho1;
 					for(INMOST_DATA_ENUM_TYPE i = 0; i < j+1; i++)
 						for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
 							u[i][k] = r[i][k] - beta*u[i][k];
+
+					ApplyOperator(u[j],u[j+1]); // u[j+1] = A*R*u[j]
+					eta = info->ScalarProd(u[j+1],r_tilde,vlocbeg,vlocend); //eta = dot(u[j+1],r_tilde)
+					if( fabs(eta) < 1.0e-54 ) 
 					{
-						ttt = Timer();
-						if (prec != NULL) prec->Solve(u[j], t);
-						tp += Timer() - ttt;
-						if( is_parallel ) info->Update(t);
-						ttt = Timer();
-						Alink->MatVec(1.0,t,0,u[j+1]);
-						ts += Timer() - ttt;
-						if( is_parallel ) info->Update(u[j+1]);
-					}
-					{
-						temp[0] = 0;
-						for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k < vlocend; ++k)
-							temp[0] += u[j+1][k]*r0[k];
-						if( is_parallel ) info->Integrate(temp,1);
-						if( fabs(temp[0]) < 1.0e-54 ) 
-						{
-							reason = "denominator(2) is zero";
-							goto exit;
-						}
-						alpha = rho0 / temp[0];
+						reason = "denominator(2) is zero";
+						goto exit;
 					}
+					alpha = rho0 / eta;
+
+					for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
+						SOL[k] += alpha*u[0][k];
+
 					for(INMOST_DATA_ENUM_TYPE i = 0; i < j+1; i++)
-						for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
+						for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k) //r[i] = r[i] - alpha * u[i+1]
 							r[i][k] -= alpha*u[i+1][k];
+
+					
+					resid = info->ScalarProd(r[j],r[j],vlocbeg,vlocend); // resid = dot(r[j],r[j])
+					resid = sqrt(resid); // resid = sqrt(dot(r[j],r[j]))
+
+					if( resid < atol || resid < rtol*resid0 ) 
 					{
-						ttt = Timer();
-						if (prec != NULL) prec->Solve(r[j], t);
-						tp += Timer() - ttt;
-						if( is_parallel ) info->Update(t);
-						ttt = Timer();
-						Alink->MatVec(1.0,t,0,r[j+1]);
-						ts += Timer() - ttt;
-						if( is_parallel ) info->Update(r[j+1]);
+						reason = "early exit in bi-cg block";
+						last_resid = resid;
+						goto exit;
 					}
-					for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
-						SOL[k] += alpha*u[0][k];
+					
+
+					ApplyOperator(r[j],r[j+1]); // r[j+1] = A*R*r[j]
+
+					
 				}
 
 				for(INMOST_DATA_ENUM_TYPE j = 1; j < l+1; j++)
@@ -244,7 +260,7 @@ namespace INMOST
 						tau[i-1 + (j-1)*l] = 0;
 						for(INMOST_DATA_ENUM_TYPE k = vlocbeg; k < vlocend; ++k)
 							tau[i-1 + (j-1)*l] += r[j][k]*r[i][k];
-						if(is_parallel) info->Integrate(&tau[i-1 + (j-1)*l],1);
+						info->Integrate(&tau[i-1 + (j-1)*l],1);
 						tau[i-1 + (j-1)*l] /= sigma[i-1];
 						for(INMOST_DATA_ENUM_TYPE k = vbeg; k < vend; ++k)
 							r[j][k] -= tau[i-1 + (j-1)*l]*r[i][k];
@@ -256,24 +272,24 @@ namespace INMOST
 						temp[0] += r[j][k]*r[j][k];
 						temp[1] += r[0][k]*r[j][k];
 					}
-					if(is_parallel) info->Integrate(temp,2);
+					info->Integrate(temp,2);
 					sigma[j-1] = temp[0]+1.0e-35; //REVIEW
 					theta2[j-1] = temp[1]/sigma[j-1];
 				}
 				omega =