SolverK3BIILU2 refactored using SolverInterface

Source/Headers/CMakeLists.txt

@@ -13,7 +13,7 @@ set(HEADER
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_mesh.h
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver.h
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_interface.h
-        ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_factory
+        ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver_master.h
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_partitioner.h
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_autodiff.h
         ${CMAKE_CURRENT_SOURCE_DIR}/inmost_expression.h

Source/Headers/inmost.h

@@ -7,7 +7,7 @@
 #include "inmost_dense.h"
 #include "inmost_solver.h"
 #include "inmost_solver_interface.h"
-#include "inmost_solver_factory.h"
+#include "inmost_solver_master.h"
 #include "inmost_partitioner.h"
 #include "inmost_variable.h"
 #include "inmost_nonlinear.h"

Source/Headers/inmost_solver.h

-//
-// Created by Dmitri Bagaev on 22.09.16.
-//
-
 #ifndef INMOST_SOLVER_INCLUDED
 #define INMOST_SOLVER_INCLUDED
 

Source/Solver/CMakeLists.txt

@@ -8,6 +8,8 @@ set(SOURCE
         )
 
 add_subdirectory(solver_inner)
+#add_subdirectory(solver_fcbiilu2)
+add_subdirectory(solver_k3biilu2)
 
 if (USE_SOLVER_PETSC)
     add_subdirectory(solver_petsc)
@@ -25,28 +27,6 @@ if (USE_SOLVER_SUPERLU)
     add_subdirectory(solver_superlu)
 endif()
 
-if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h" AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp" AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/fcbiilu2.cpp")
-    set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_FCBIILU2)
-    set(HAVE_SOLVER_FCBIILU2 TRUE PARENT_SCOPE)
-    list(APPEND HEADER ${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.h)
-    list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_fcbiilu2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/fcbiilu2.cpp)
-else ()
-    set(HAVE_SOLVER_FCBIILU2 FALSE)
-endif ()
-if (USE_MPI)
-    if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h" AND EXISTS
-            "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp" AND EXISTS
-            "${CMAKE_CURRENT_SOURCE_DIR}/k3d.h" AND EXISTS
-            "${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp")
-        set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_K3BIILU2)
-        set(HAVE_SOLVER_K3BIILU2 TRUE PARENT_SCOPE)
-        list(APPEND HEADER ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h ${CMAKE_CURRENT_SOURCE_DIR}/k3d.h)
-        list(APPEND SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp)
-    else ()
-        set(HAVE_SOLVER_K3BIILU2 FALSE PARENT_SCOPE)
-    endif ()
-endif ()
-
 
 set(HEADER ${HEADER} PARENT_SCOPE)
 set(SOURCE ${SOURCE} PARENT_SCOPE)

Source/Solver/solver.cpp

@@ -25,6 +25,10 @@
 #include "solver_superlu/SolverSUPERLU.h"
 #endif
 
+#if defined (HAVE_SOLVER_K3BIILU2)
+#include "solver_k3biilu2/SolverK3BIILU2.h"
+#endif
+
 namespace INMOST {
 
     int *Solver::argc = NULL;
@@ -109,6 +113,9 @@ namespace INMOST {
 #endif
 #if defined(USE_SOLVER_SUPERLU)
         SolverMaster::registerSolver<SolverSUPERLU>("superlu");
+#endif
+#if defined(HAVE_SOLVER_K3BIILU2)
+        SolverMaster::registerSolver<SolverK3BIILU2>("k3biilu2");
 #endif
         Sparse::CreateRowEntryType();
     }

Source/Solver/solver_k3biilu2/CMakeLists.txt

+if (USE_MPI)
+    if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h" AND EXISTS
+            "${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp" AND EXISTS
+            "${CMAKE_CURRENT_SOURCE_DIR}/k3d.h" AND EXISTS
+            "${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp")
+        set(SOLVER_DEFINITIONS ${SOLVER_DEFINITIONS} -DHAVE_SOLVER_K3BIILU2 PARENT_SCOPE)
+        set(HAVE_SOLVER_K3BIILU2 TRUE PARENT_SCOPE)
+        set(HEADER ${HEADER} 
+            ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.h 
+            ${CMAKE_CURRENT_SOURCE_DIR}/k3d.h
+            ${CMAKE_CURRENT_SOURCE_DIR}/SolverK3BIILU2.h PARENT_SCOPE)
+
+        set(SOURCE ${SOURCE} 
+            ${CMAKE_CURRENT_SOURCE_DIR}/solver_k3biilu2.cpp 
+            ${CMAKE_CURRENT_SOURCE_DIR}/k3d.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/SolverK3BIILU2.cpp PARENT_SCOPE)
+    else ()
+        set(HAVE_SOLVER_K3BIILU2 FALSE PARENT_SCOPE)
+    endif ()
+endif()
\ No newline at end of file

Source/Solver/solver_k3biilu2/SolverK3BIILU2.cpp

+#include "SolverK3BIILU2.h"
+
+namespace INMOST {
+
+	SolverK3BIILU2::SolverK3BIILU2() {
+		
+	}
+
+    SolverK3BIILU2::SolverK3BIILU2(const SolverInterface *other) {
+    	const SolverK3BIILU2 *k3other = static_cast<const SolverK3BIILU2 *>(other);
+    	SolverCopyDataK3biilu2(&solver_data, k3other->solver_data, k3other->communicator);
+		if (k3other->matrix_data != NULL) {
+			MatrixCopyDataK3biilu2(&matrix_data, k3other->matrix_data);
+			SolverSetMatrixK3biilu2(solver_data, matrix_data, false, false);
+		}
+    }
+
+	void SolverK3BIILU2::Assign(const SolverInterface *other) {
+		const SolverK3BIILU2 *k3other = static_cast<const SolverK3BIILU2 *>(other);
+		SolverAssignDataK3biilu2(solver_data, k3other->solver_data);
+		if (k3other->matrix_data != NULL) {
+			if (matrix_data != NULL ) {
+				MatrixAssignDataK3biilu2(matrix_data, k3other->matrix_data);
+			} else {
+				MatrixCopyDataK3biilu2(&matrix_data, k3other->matrix_data);
+			}
+			SolverSetMatrixK3biilu2(solver_data, matrix_data, false, false);
+		}
+	}
+
+	void SolverK3BIILU2::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {
+		SolverInitDataK3biilu2(&solver_data, communicator, prefix.c_str());
+		SolverInitializeK3biilu2(argc, argv, parameters_file);
+	}
+
+	void SolverK3BIILU2::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
+		bool modified_pattern = ModifiedPattern;
+			//~ if( A.comm != comm ) throw DifferentCommunicatorInSolver;
+		if (matrix_data == NULL) {
+			MatrixInitDataK3biilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
+				modified_pattern = true;
+		}
+		if (modified_pattern) {
+			global_size = local_size = A.Size();
+		
+			MatrixDestroyDataK3biilu2(&matrix_data);
+			MatrixInitDataK3biilu2(&matrix_data, A.GetCommunicator(), A.GetName().c_str());
+			INMOST_DATA_ENUM_TYPE nnz = 0, k = 0, q = 1, shift = 0;
+			int nproc, myid;
+#if defined(USE_MPI)
+			MPI_Comm_size(A.GetCommunicator(), &nproc);
+			MPI_Comm_rank(A.GetCommunicator(), &myid);
+#else
+			nproc = 1;
+			myid = 0;
+#endif
+			int * ibl = (int *) malloc(sizeof(int) * (nproc + 1));
+			ibl[0] = 0;
+			int n = A.Size();
+#if defined(USE_MPI)
+			INMOST_DATA_ENUM_TYPE mbeg,mend;
+			A.GetInterval(mbeg, mend);
+			n = mend - mbeg;
+			int block_end = mend;
+			MPI_Allgather(&block_end, 1, MPI_INT, &ibl[1], 1, MPI_INT, A.GetCommunicator());
+#else
+			ibl[1] = n;
+#endif
+			int * ia = (int *) malloc(sizeof(int) * (n + 1));
+			for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
+				nnz += it->Size();
+			}
+			int * ja = (int *) malloc(sizeof(int) * nnz);
+			double * values = (double *) malloc(sizeof(double) * nnz);
+			ia[0] = shift;
+			for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
+				for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++) {
+					ja[k] = jt->first + 0;
+					values[k] = jt->second;
+		                                //std::cout<<"# q="<<q<<" k="<<k<<" ja="<<ja[k]<<" a="<<values[k]<<std::endl;//db!
+					k++;
+				}
+				shift += it->Size();
+				ia[q++] = shift;
+			}
+			MatrixFillK3biilu2(matrix_data,n,nproc,ibl,ia,ja,values);
+		} else {
+			INMOST_DATA_ENUM_TYPE nnz = 0, k = 0;
+			for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
+				nnz += it->Size();
+			}
+			double * values = (double *) malloc(sizeof(double) * nnz);
+			k = 0;
+			for(Sparse::Matrix::iterator it = A.Begin(); it != A.End(); it++) {
+				for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); jt++) {
+					values[k++] = jt->second;
+				}
+			}
+			MatrixFillValuesK3biilu2(matrix_data,values);
+		}
+		MatrixFinalizeK3biilu2(matrix_data);
+		SolverSetMatrixK3biilu2(solver_data, matrix_data, modified_pattern, OldPreconditioner);
+	}
+
+	bool SolverK3BIILU2::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
+		INMOST_DATA_ENUM_TYPE vbeg, vend;
+        RHS.GetInterval(vbeg, vend);
+		void *rhs_data = NULL;
+		void *solution_data = NULL;
+		
+		VectorInitDataK3biilu2(&rhs_data, RHS.GetCommunicator(), RHS.GetName().c_str());
+		VectorPreallocateK3biilu2(rhs_data, local_size);
+			
+		VectorInitDataK3biilu2(&solution_data, SOL.GetCommunicator(), SOL.GetName().c_str());
+		VectorPreallocateK3biilu2(solution_data,local_size);	
+			
+		VectorFillK3biilu2(rhs_data, &RHS[vbeg]);
+		VectorFinalizeK3biilu2(rhs_data);
+
+		VectorFillK3biilu2(solution_data, &SOL[vbeg]);
+		VectorFinalizeK3biilu2(solution_data);
+
+		bool result = SolverSolveK3biilu2(solver_data, rhs_data, solution_data);
+		if (result) VectorLoadK3biilu2(solution_data, &SOL[vbeg]);
+		return result;
+	}
+
+	bool SolverK3BIILU2::Clear() {
+		if (matrix_data != NULL) {
+			MatrixDestroyDataK3biilu2(&matrix_data);
+		}
+		SolverDestroyDataK3biilu2(&solver_data);
+		return true;
+	}
+
+	bool SolverK3BIILU2::isMatrixSet() {
+		return matrix_data != NULL;
+	}
+
+	void SolverK3BIILU2::SetDefaultParameters() {
+
+	}
+
+	SolverParameter SolverK3BIILU2::GetParameter(std::string name) const {
+		throw INMOST::SolverUnsupportedOperation;
+	}
+
+	void SolverK3BIILU2::SetParameter(std::string name, std::string value) {
+		throw INMOST::SolverUnsupportedOperation;
+	}
+
+	const INMOST_DATA_ENUM_TYPE SolverK3BIILU2::Iterations() const {
+		return SolverIterationNumberK3biilu2(solver_data);
+	}
+
+	const INMOST_DATA_REAL_TYPE SolverK3BIILU2::Residual() const {
+		return SolverResidualNormK3biilu2(solver_data);
+	}
+
+	const std::string SolverK3BIILU2::ReturnReason() const {
+		return "Unspecified for K3BIILU2";
+	}
+
+	const std::string SolverK3BIILU2::SolverName() const {
+		return "k3biilu2";
+	}
+
+	void SolverK3BIILU2::Finalize() {
+		SolverFinalizeK3biilu2();
+	}
+
+	SolverK3BIILU2::~SolverK3BIILU2() {
+		this->Clear();
+	}
+
+}
\ No newline at end of file

Source/Solver/solver_k3biilu2/SolverK3BIILU2.h

+#ifndef INMOST_SOLVERK3BIILU2_H
+#define INMOST_SOLVERK3BIILU2_H
+
+#include "inmost_solver_interface.h"
+#include "solver_k3biilu2.h"
+
+namespace INMOST {
+
+    class SolverK3BIILU2 : public SolverInterface {
+    private:
+        void *solver_data;
+        void *matrix_data;
+        INMOST_DATA_ENUM_TYPE local_size, global_size;
+    public:
+        SolverK3BIILU2();
+
+        SolverK3BIILU2(const SolverInterface *other);
+
+        virtual void Assign(const SolverInterface *other);
+
+        virtual void Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix);
+
+        virtual void SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner);
+
+        virtual bool Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL);
+
+        virtual bool Clear();
+
+        virtual bool isMatrixSet();
+
+        virtual void SetDefaultParameters();
+
+        virtual SolverParameter GetParameter(std::string name) const;
+
+        virtual void SetParameter(std::string name, std::string value);
+
+        virtual const INMOST_DATA_ENUM_TYPE Iterations() const;
+
+        virtual const INMOST_DATA_REAL_TYPE Residual() const;
+
+        virtual const std::string ReturnReason() const;
+
+        virtual const std::string SolverName() const;
+
+        virtual void Finalize();
+
+        virtual ~SolverK3BIILU2();
+
+    };
+
+}
+
+
+#endif //INMOST_SOLVERK3BIILU2_H

Source/Solver/solver_k3biilu2/solver_k3biilu2.h

+#ifndef SOLVER_K3BIILU2_H_INCLUDED
+#define SOLVER_K3BIILU2_H_INCLUDED
+#include "inmost_solver.h"
+
+
+//#define USE_SOLVER_K3BIILU2 //just for test; see flag HAVE_SOLVER_K3BIILU2
+//#if defined(USE_SOLVER_K3BIILU2)
+void MatrixCopyDataK3biilu2(void ** ppA, void * pB);
+void MatrixAssignDataK3biilu2(void * pA, void * pB);
+void MatrixInitDataK3biilu2(void ** ppA, INMOST_MPI_Comm comm, const char * name);
+void MatrixDestroyDataK3biilu2(void ** pA);
+void MatrixFillK3biilu2(void * pA, int size, int nproc, int * ibl, int * ia, int * ja, double * values);
+void MatrixFillValuesK3biilu2(void * pA, double * values);
+void MatrixFinalizeK3biilu2(void * data);
+
+void VectorInitDataK3biilu2(void ** ppA, INMOST_MPI_Comm comm, const char * name);
+void VectorCopyDataK3biilu2(void ** ppA, void * pB);
+void VectorAssignDataK3biilu2(void * pA, void * pB);
+void VectorPreallocateK3biilu2(void * pA, int size);
+void VectorFillK3biilu2(void * pA, double * values);
+void VectorLoadK3biilu2(void * pA, double * values);
+void VectorFinalizeK3biilu2(void * data);
+void VectorDestroyDataK3biilu2(void ** ppA);
+
+void SolverInitializeK3biilu2(int * argc,char *** argv, const char * file_options);
+bool SolverIsFinalizedK3biilu2();
+void SolverFinalizeK3biilu2();
+void SolverDestroyDataK3biilu2(void ** data);
+void SolverInitDataK3biilu2(void ** data, INMOST_MPI_Comm comm, const char * name);
+void SolverCopyDataK3biilu2(void ** data, void * other_data, INMOST_MPI_Comm comm);
+void SolverAssignDataK3biilu2(void * data, void * other_data);
+void SolverSetMatrixK3biilu2(void * data, void * matrix_data, bool same_pattern, bool reuse_preconditioner);
+bool SolverSolveK3biilu2(void * data, void * rhs_data, void * sol_data);
+int SolverIterationNumberK3biilu2(void * data);
+double SolverResidualNormK3biilu2(void * data);
+void SolverAddOtherStatK3biilu2(void * data, unsigned int * pivmod, double * prdens, double * t_prec, double * t_iter);
+//#endif //USE_K3SOLVER_BIILU2
+
+
+#endif //SOLVER_K3BIILU2_H_INCLUDED

Source/Solver/solver_petsc/solver_petsc.h

-//
-// Created by Dmitri Bagaev on 28.09.16.
-//
-
 #ifndef INMOST_NEW_SOLVER_PETSC_H
 #define INMOST_NEW_SOLVER_PETSC_H
 

Source/Solver/solver_superlu/SolverSUPERLU.cpp

@@ -23,6 +23,7 @@ namespace INMOST {
         }
 
         void SolverSUPERLU::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
+        	//check that the run is serial!
         	int *ia, *ja, nnz = 0;
 			double *a;
 			int mbeg = A.GetFirstIndex();