diff --git a/CMakeLists.txt b/CMakeLists.txt index fb03317eaf1f030bdbb2e0ad2cfdc965e55819ba..e7d014f5fbdbab7e4676fb727a6262663568dbf6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,6 +34,7 @@ option(USE_SOLVER_METIS "Use METIS for matrix reordering" OFF) option(USE_SOLVER_MONDRIAAN "Use Mondriaan for matrix reordering" OFF) option(USE_SOLVER_PETSC "Use PETSc solvers" OFF) option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF) +option(USE_SOLVER_SUPERLU "Use SuperLU solver" OFF) #option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF) #option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF) #option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF) @@ -183,6 +184,22 @@ if(USE_SOLVER_TRILINOS) endif() +if(USE_SOLVER_SUPERLU) + find_package(SUPERLU) + if(NOT SUPERLU_FOUND) + set(USE_SOLVER_SUPERLU OFF) + MESSAGE("SUPERLU NOT FOUND") + else() + MESSAGE("\nFound SuperLU! Here are the details: ") + MESSAGE("INCLUDES: ${SUPERLU_INCLUDES}") + MESSAGE("LIBRARIES: ${SUPERLU_LIBRARIES}") + + include_directories(${SUPERLU_INCLUDES}) + set(USE_SOLVER_SUPERLU ON) + + endif() +endif() + if(USE_AUTODIFF_OPENCL) find_package(OpenCL) if(OPENCL_FOUND) diff --git a/Examples/ADFVDiscr/CMakeLists.txt b/Examples/ADFVDiscr/CMakeLists.txt index 24df7cb5ed9cf0080d57daf4cf6aaa586fd7d403..f81abafe1aa4dd631303adaa804e3ef7b11d4cbd 100644 --- a/Examples/ADFVDiscr/CMakeLists.txt +++ b/Examples/ADFVDiscr/CMakeLists.txt @@ -30,6 +30,10 @@ if(USE_SOLVER) message("linking ADFVDiscr with Mondriaan") target_link_libraries(ADFVDiscr ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking ADFVDiscr with SuperLU") + target_link_libraries(ADFVDiscr ${SUPERLU_LIBRARIES}) + endif() endif() diff --git a/Examples/ADMFD/CMakeLists.txt b/Examples/ADMFD/CMakeLists.txt index 1397294925769f6a15985a12bf14c5ef685ffc6e..ceb33316da4d38136161481381c6ddc28b79e8d8 100644 --- a/Examples/ADMFD/CMakeLists.txt +++ b/Examples/ADMFD/CMakeLists.txt @@ -29,6 +29,10 @@ if(USE_SOLVER) message("linking ADMFD with Mondriaan") target_link_libraries(ADMFD ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking ADMFD with SuperLU") + target_link_libraries(ADMFD ${SUPERLU_LIBRARIES}) + endif() endif() diff --git a/Examples/ADMFD/main.cpp b/Examples/ADMFD/main.cpp index 14308d02b67508342397caa66e78a104f0a3ae72..198b629521c7ce7903a58be6fdda2ca3183be7de 100644 --- a/Examples/ADMFD/main.cpp +++ b/Examples/ADMFD/main.cpp @@ -449,8 +449,8 @@ int main(int argc,char ** argv) //Change matrix by nonlinear correction for DMP ++total; - if( !isDMP ) - //if( false ) + //if( !isDMP ) + if( false ) { const real var = 0;//0.25; vMatrix vnKGRAD = nKGRAD; @@ -557,11 +557,12 @@ int main(int argc,char ** argv) if( R.Norm() < 1.0e-4 ) break; Solver S(Solver::INNER_MPTILUC); - S.SetMatrix(R.GetJacobian()); - S.SetParameterReal("relative_tolerance", 1.0e-14); + //Solver S(Solver::SUPERLU); + S.SetParameterReal("relative_tolerance", 1.0e-14); S.SetParameterReal("absolute_tolerance", 1.0e-12); S.SetParameterReal("drop_tolerance", 1.0e-2); S.SetParameterReal("reuse_tolerance", 1.0e-4); + S.SetMatrix(R.GetJacobian()); //std::fill(Update.Begin(),Update.End(),0.0); if( S.Solve(R.GetResidual(),Update) ) { diff --git a/Examples/DrawMatrix/CMakeLists.txt b/Examples/DrawMatrix/CMakeLists.txt index c53110192528849e9e0e3f0682eca387b2e1b129..5773fcc3fee6204486cdf19be5827504f8746b7a 100644 --- a/Examples/DrawMatrix/CMakeLists.txt +++ b/Examples/DrawMatrix/CMakeLists.txt @@ -16,29 +16,7 @@ if(OPENGL_FOUND) set_target_properties(DrawMatrix PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") endif() endif(USE_MPI) - if(USE_SOLVER) - if(USE_SOLVER_ANI3D) - message("linking DrawMatrix with ANI3d") - target_link_libraries(DrawMatrix ani3d) - endif() - if(USE_SOLVER_PETSC) - message("linking DrawMatrix with PETSc") - add_definitions(${PETSC_DEFINITIONS}) - target_link_libraries(DrawMatrix ${PETSC_LIBRARIES}) - endif() - if(USE_SOLVER_TRILINOS) - message("linking DrawMatrix with Trilinos") - target_link_libraries(DrawMatrix ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES}) - endif() - if(USE_SOLVER_METIS) - message("linking DrawMatrix with Metis") - target_link_libraries(DrawMatrix ${METIS_LIBRARIES}) - endif() - if(USE_SOLVER_MONDRIAAN) - message("linking DrawMatrix with Mondriaan") - target_link_libraries(DrawMatrix ${MONDRIAAN_LIBRARIES}) - endif() - endif() + install(TARGETS DrawMatrix EXPORT inmost-targets RUNTIME DESTINATION bin) diff --git a/Examples/FVDiscr/CMakeLists.txt b/Examples/FVDiscr/CMakeLists.txt index 373a1701ef737f3fb915672655e647d75ee13b04..2012e75e7b81f4977270e0dbeecb63b930ef582c 100644 --- a/Examples/FVDiscr/CMakeLists.txt +++ b/Examples/FVDiscr/CMakeLists.txt @@ -29,6 +29,10 @@ if(USE_SOLVER) message("linking FVDiscr with Mondriaan") target_link_libraries(FVDiscr ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking FVDiscr with SuperLU") + target_link_libraries(FVDiscr ${SUPERLU_LIBRARIES}) + endif() endif() diff --git a/Examples/MatSolve/CMakeLists.txt b/Examples/MatSolve/CMakeLists.txt index b1333dfe32d996b653007329410febc1918be5dd..50b25dae718652181afe5d2be48904c103e59676 100644 --- a/Examples/MatSolve/CMakeLists.txt +++ b/Examples/MatSolve/CMakeLists.txt @@ -37,6 +37,10 @@ if(USE_SOLVER) message("linking MatSolve with Mondriaan") target_link_libraries(MatSolve ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking MatSolve with SuperLU") + target_link_libraries(MatSolve ${SUPERLU_LIBRARIES}) + endif() endif() diff --git a/Examples/Solver/CMakeLists.txt b/Examples/Solver/CMakeLists.txt index 4463e40575effe99c555254f35878b5dbf47084d..f98254da70a52b3d178156c5ecfa4786d8d4a707 100644 --- a/Examples/Solver/CMakeLists.txt +++ b/Examples/Solver/CMakeLists.txt @@ -29,6 +29,10 @@ if(USE_SOLVER) message("linking Solver with Mondriaan") target_link_libraries(Solver ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking Solver with SuperLU") + target_link_libraries(Solver ${SUPERLU_LIBRARIES}) + endif() endif(USE_SOLVER) if(USE_PARTITIONER) diff --git a/Source/Headers/inmost_options_cmake.h b/Source/Headers/inmost_options_cmake.h index 79d7a5b0fab620c9e82332404c80c08454cb2b65..47ccdda26b3a559c681fd86f9ac6d50d00b684a3 100644 --- a/Source/Headers/inmost_options_cmake.h +++ b/Source/Headers/inmost_options_cmake.h @@ -20,6 +20,7 @@ #cmakedefine USE_SOLVER_PETSC #cmakedefine USE_SOLVER_TRILINOS #cmakedefine USE_SOLVER_ANI +#cmakedefine USE_SOLVER_SUPERLU #cmakedefine USE_NONLINEAR #cmakedefine USE_NONLINEAR_TRILINOS diff --git a/Source/Headers/inmost_solver.h b/Source/Headers/inmost_solver.h index d0aa2183db27f650819a72e8c8fb2d2c33f511a4..93230b66de6bf0d914ed01e799e6307ab868fd0c 100644 --- a/Source/Headers/inmost_solver.h +++ b/Source/Headers/inmost_solver.h @@ -54,7 +54,8 @@ namespace INMOST PETSc, ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/ ANI, ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/ FCBIILU2, ///< external FCBIILU2 Solver (BIILU2 parallel F2C version). - K3BIILU2 ///< inner K3BIILU2 Solver (BIILU2 parallel version). + K3BIILU2, ///< inner K3BIILU2 Solver (BIILU2 parallel version). + SUPERLU ///< external Solver SuperLU @see https://github.com/starseeker/SuperLU }; static std::string TypeName(Type t); diff --git a/Source/IO/mesh_vtk_file.cpp b/Source/IO/mesh_vtk_file.cpp index 340bb4f7bcd47a251a5da53c7721d426eff768ac..b4a2db46db2b502816dafc1d2a6c2c05f2857213 100644 --- a/Source/IO/mesh_vtk_file.cpp +++ b/Source/IO/mesh_vtk_file.cpp @@ -273,7 +273,11 @@ safe_output: case DATA_REAL: { Storage::real_array arr = it->RealArray(tags[i]); - for(unsigned int m = 0; m < comps; m++) fprintf(f,"%14e ",arr[m]); + for(unsigned int m = 0; m < comps; m++) + { + double val = static_cast(arr[m]); + fprintf(f,"%14e ",val != val ? -0.9999E30 : val); + } fprintf(f,"\n"); } break; @@ -288,7 +292,11 @@ safe_output: case DATA_VARIABLE: { Storage::var_array arr = it->VariableArray(tags[i]); - for(unsigned int m = 0; m < comps; m++) fprintf(f,"%14e ",arr[m].GetValue()); + for(unsigned int m = 0; m < comps; m++) + { + double val = static_cast(arr[m].GetValue()); + fprintf(f,"%14e ",val != val ? -0.9999E30 : val); + } fprintf(f,"\n"); } break; @@ -346,7 +354,11 @@ safe_output: case DATA_REAL: { Storage::real_array arr = it->RealArray(tags[i]); - for(unsigned int m = 0; m < comps; m++) fprintf(f,"%14e ",arr[m]); + for(unsigned int m = 0; m < comps; m++) + { + double val = static_cast(arr[m]); + fprintf(f,"%14e ",(val != val ? -0.9999E30 : val)); + } fprintf(f,"\n"); } break; @@ -361,7 +373,11 @@ safe_output: case DATA_VARIABLE: { Storage::var_array arr = it->VariableArray(tags[i]); - for(unsigned int m = 0; m < comps; m++) fprintf(f,"%14e ",arr[m].GetValue()); + for(unsigned int m = 0; m < comps; m++) + { + double val = static_cast(arr[m].GetValue()); + fprintf(f,"%14e ",(val != val ? -0.9999E30 : val)); + } fprintf(f,"\n"); } break; diff --git a/Source/Solver/CMakeLists.txt b/Source/Solver/CMakeLists.txt index ede096437b6ef8e9b7c9f53c0c9a86bf16dad4a1..8d10349e93383036dabc85adb9a74efc14e0aeff 100644 --- a/Source/Solver/CMakeLists.txt +++ b/Source/Solver/CMakeLists.txt @@ -5,12 +5,14 @@ set(SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/solver_ani.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_petsc.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/solver_superlu.cpp ) set(HEADER ${HEADER} ${CMAKE_CURRENT_SOURCE_DIR}/solver_prototypes.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_petsc.h + ${CMAKE_CURRENT_SOURCE_DIR}/solver_superlu.h ${CMAKE_CURRENT_SOURCE_DIR}/solver_ilu2.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_bcgsl.hpp diff --git a/Source/Solver/solver.cpp b/Source/Solver/solver.cpp index 352254aa4073095a2f992b58e43b35a8f495ade4..edc973c884949eec1f6df7e1755bd267d61e2564 100644 --- a/Source/Solver/solver.cpp +++ b/Source/Solver/solver.cpp @@ -2,6 +2,7 @@ #include "inmost_solver.h" #if defined(USE_SOLVER) #include "solver_petsc.h" +#include "solver_superlu.h" #include "solver_ani.h" #include "solver_ilu2.hpp" #include "solver_ddpqiluc2.hpp" @@ -43,11 +44,10 @@ #include "Teuchos_TestForException.hpp" #include "Teuchos_XMLParameterListHelpers.hpp" #endif - //#define USE_OMP -//#define KSOLVER BCGSL_solver -#define KSOLVER BCGS_solver +#define KSOLVER BCGSL_solver +//#define KSOLVER BCGS_solver //#define ACCELERATED_CONDEST //#define PRINT_CONDEST @@ -61,6 +61,7 @@ namespace INMOST static std::string ani_database_file = ""; static std::string fcbiilu2_database_file = ""; static std::string k3biilu2_database_file = ""; + static std::string superlu_database_file = ""; @@ -86,6 +87,7 @@ namespace INMOST case ANI: return "ANI"; case FCBIILU2: return "FCBIILU2"; case K3BIILU2: return "K3BIILU2"; + case SUPERLU: return "SUPERLU"; } return "Unknown"; } @@ -852,53 +854,58 @@ namespace INMOST (void)argv; if( database != NULL ) { - FILE * f = fopen(database,"r"); - if( f ) - { - //std::fstream file(database,std::ios::in); - char str[4096]; - //while( !file.eof() && file.getline(str,4096) ) - while( !feof(f) && fgets(str,4096,f) ) - { - int k = 0, l; - for(k = 0; k < (int)strlen(str); ++k) - { - if( str[k] == ':' ) break; - } - if( k == strlen(str) ) continue; //invalid line - for(l = 0; l < k; ++l) str[l] = tolower(str[l]); - l = (int)strlen(str)-1; // Right-trim string - while(l > 0 && isspace(str[l]) ) --l; - str[l+1] = 0; - l = k+1; - while(l < (int)strlen(str) && isspace(str[l]) ) ++l; - if( l == strlen(str) ) continue; //skip empty entry - if( !strncmp(str,"petsc",k) ) - petsc_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_ifpack",k) ) - trilinos_ifpack_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_aztec",k) ) - trilinos_aztec_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_ml",k) ) - trilinos_ml_database_file = std::string(str+l); - else if( !strncmp(str,"trilinos_belos",k) ) - trilinos_belos_database_file = std::string(str+l); - else if( !strncmp(str,"ani",k) ) - ani_database_file = std::string(str+l); - else if( !strncmp(str,"fcbiilu2",k) ) - fcbiilu2_database_file = std::string(str+l); - else if( !strncmp(str,"k3biilu2",k) ) - k3biilu2_database_file = std::string(str+l); - } - //file.close(); - fclose(f); - } + FILE * f = fopen(database,"r"); + if( f ) + { + //std::fstream file(database,std::ios::in); + char str[4096]; + //while( !file.eof() && file.getline(str,4096) ) + while( !feof(f) && fgets(str,4096,f) ) + { + int k = 0, l; + for(k = 0; k < (int)strlen(str); ++k) + { + if( str[k] == ':' ) break; + } + if( k == strlen(str) ) continue; //invalid line + for(l = 0; l < k; ++l) str[l] = tolower(str[l]); + l = (int)strlen(str)-1; // Right-trim string + while(l > 0 && isspace(str[l]) ) --l; + str[l+1] = 0; + l = k+1; + while(l < (int)strlen(str) && isspace(str[l]) ) ++l; + if( l == strlen(str) ) continue; //skip empty entry + if( !strncmp(str,"petsc",k) ) + petsc_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_ifpack",k) ) + trilinos_ifpack_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_aztec",k) ) + trilinos_aztec_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_ml",k) ) + trilinos_ml_database_file = std::string(str+l); + else if( !strncmp(str,"trilinos_belos",k) ) + trilinos_belos_database_file = std::string(str+l); + else if( !strncmp(str,"ani",k) ) + ani_database_file = std::string(str+l); + else if( !strncmp(str,"fcbiilu2",k) ) + fcbiilu2_database_file = std::string(str+l); + else if( !strncmp(str,"k3biilu2",k) ) + k3biilu2_database_file = std::string(str+l); + else if( !strncmp(str,"superlu",k) ) + superlu_database_file = std::string(str+l); + } + //file.close(); + fclose(f); + } } //std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl; //std::cout << "Trilinos_Ifpack \"" << trilinos_ifpack_database_file << "\"" << std::endl; #if defined(USE_SOLVER_PETSC) SolverInitializePetsc(argc,argv,petsc_database_file.c_str()); #endif +#if defined(USE_SOLVER_SUPERLU) + SolverInitializeSuperLU(argc,argv,superlu_database_file.c_str()); +#endif #if defined(USE_SOLVER_ANI) SolverInitializeAni(argc,argv,ani_database_file.c_str()); #endif @@ -922,9 +929,12 @@ namespace INMOST } } } +#endif +#if defined(USE_SOLVER_SUPERLU) + #endif Sparse::CreateRowEntryType(); - is_initialized = true; + is_initialized = true; } void Solver::Finalize() @@ -1070,6 +1080,12 @@ namespace INMOST SolverInitDataK3biilu2(&solver_data,_comm,name.c_str()); } #endif +#if defined(USE_SOLVER_SUPERLU) + if( _pack == SUPERLU ) + { + SolverInitDataSuperLU(&solver_data); + } +#endif #if defined(USE_SOLVER_TRILINOS) if( _pack == Trilinos_Aztec || _pack == Trilinos_Belos || @@ -1271,6 +1287,13 @@ namespace INMOST solver_data = NULL; } } +#endif +#if defined(USE_SOLVER_SUPERLU) + if(_pack == SUPERLU ) + { + SolverDestroyDataSuperLU(&solver_data); + matrix_data = NULL; + } #endif if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2) { @@ -1721,6 +1744,40 @@ namespace INMOST ok = true; } #endif +#if defined(USE_SOLVER_SUPERLU) + if( _pack == SUPERLU ) //just serial SuperLU + { + //check that the run is serial! + int * ia, * ja, nnz = 0; + double * a; + int mbeg = A.GetFirstIndex(); + int mend = A.GetLastIndex(); + for(int k = 0; k < mend-mbeg; ++k) nnz += A[k].Size(); + ia = (int *)malloc(sizeof(int)*(A.Size()+1)); + ja = (int *)malloc(sizeof(int)*nnz); + a = (double *)malloc(sizeof(double)*nnz); + int q = 0; + ia[0] = 0; + for(int k = 0; k < mend-mbeg; ++k) + { + Sparse::Row & r = A[k+mbeg]; + for(int l = 0; l < r.Size(); ++l) + { + ja[q] = r.GetIndex(l)-mbeg; + a[q] = r.GetValue(l); + ++q; + } + ia[k+1] = q; + } + MatrixFillSuperLU(solver_data,mend-mbeg,nnz,ia,ja,a); + //arrays are freed inside SuperLU + //delete [] ia; + //delete [] ja; + //delete [] a; + matrix_data = solver_data; //to avoid exception in Solve() + ok = true; + } +#endif #if defined(USE_SOLVER_TRILINOS) if( _pack == Trilinos_Aztec || _pack == Trilinos_Belos || @@ -2493,6 +2550,16 @@ namespace INMOST } } +#endif +#if defined(USE_SOLVER_SUPERLU) + if(_pack == SUPERLU ) + { + bool ret = SolverSolveSuperLU(solver_data,&*RHS.Begin(),&*SOL.Begin()); + last_it = 1; + last_resid = 0; + return_reason = SolverConvergedReasonSuperLU(solver_data); + return ret; + } #endif if(_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2) { diff --git a/Source/Solver/solver_ani.h b/Source/Solver/solver_ani.h index 2be0862b80d5a48e0f313c7b3eeae41f8946e60c..71b15abaafae68b37d4d96cc78ae1c6ad07dcc86 100644 --- a/Source/Solver/solver_ani.h +++ b/Source/Solver/solver_ani.h @@ -2,8 +2,6 @@ #define SOLVER_ANI_H_INCLUDED #include "inmost_solver.h" - -#if defined(USE_SOLVER_ANI) void MatrixCopyDataAni(void ** ppA, void * pB); void MatrixAssignDataAni(void * pA, void * pB); void MatrixInitDataAni(void ** ppA, INMOST_MPI_Comm comm, const char * name); @@ -32,7 +30,6 @@ void SolverSetMatrixAni(void * data, void * matrix_data, bool same_pattern, bool bool SolverSolveAni(void * data, void * rhs_data, void * sol_data); int SolverIterationNumberAni(void * data); double SolverResidualNormAni(void * data); -#endif //USE_SOLVER_ANI #endif //SOLVER_ANI_H_INCLUDED diff --git a/Source/Solver/solver_superlu.cpp b/Source/Solver/solver_superlu.cpp new file mode 100644 index 0000000000000000000000000000000000000000..af32e026df882f6540ecb504270194ab1e2cf0b1 --- /dev/null +++ b/Source/Solver/solver_superlu.cpp @@ -0,0 +1,90 @@ +#define _CRT_SECURE_NO_WARNINGS +#include "inmost_solver.h" +#if defined(USE_SOLVER) +#if defined(USE_SOLVER_SUPERLU) +#include "superlu/slu_ddefs.h" +#include "solver_superlu.h" + +struct SuperLU_data +{ + SuperMatrix A, L, U; + int * perm_r; + int * perm_c; + int size, info; + superlu_options_t options; + SuperLUStat_t stat; +}; + +void SolverInitializeSuperLU(int * argc,char *** argv, const char * file_options) +{ + //read options from file and arguments +} + + +void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values) +{ + SuperLU_data * m = static_cast(data); + dCreate_CompCol_Matrix(&m->A,size,size,nnz,col_values,col_positions,col_ranges,SLU_NR,SLU_D,SLU_GE); + m->size = size; + m->perm_c = new int [size]; + m->perm_r = new int [size]; +} + + +void SolverDestroyDataSuperLU(void ** data) +{ + SuperLU_data * m = static_cast(*data); + Destroy_CompCol_Matrix(&m->A); + Destroy_SuperNode_Matrix(&m->L); + Destroy_CompCol_Matrix(&m->U); + StatFree(&m->stat); + if( m->perm_c != NULL ) delete [] m->perm_c; + if( m->perm_r != NULL ) delete [] m->perm_r; + delete m; + *data = NULL; +} + +void SolverInitDataSuperLU(void ** data) +{ + if( data == NULL ) throw INMOST::DataCorruptedInSolver; + if( *data != NULL ) + SolverDestroyDataSuperLU(data); + *data = static_cast(new SuperLU_data); + SuperLU_data * m = static_cast(*data); + set_default_options(&(m->options)); + StatInit(&(m->stat)); + m->perm_c = NULL; + m->perm_r = NULL; + m->size = 0; +} + + +bool SolverSolveSuperLU(void * data, void * rhs_data, void * sol_data) +{ + SuperLU_data * m = static_cast(data); + assert(m->size != 0); + SuperMatrix B; + memcpy(sol_data,rhs_data,sizeof(double)*m->size); + dCreate_Dense_Matrix(&B,m->size,1,(double *)sol_data,m->size,SLU_DN,SLU_D,SLU_GE); + dgssv(&m->options,&m->A,m->perm_c,m->perm_r,&m->L,&m->U,&B,&m->stat,&m->info); + Destroy_SuperMatrix_Store(&B); + return m->info == 0; +} + +const char * SolverConvergedReasonSuperLU(void * data) +{ + SuperLU_data * m = static_cast(data); + static char reason_str[4096]; + if( m->info <= m->size ) + sprintf(reason_str,"diagonal entry of U-factor is exactly singular at %d/%d",m->info,m->size); + else if( m->info > m->size ) + sprintf(reason_str,"memory allocation failed after %d bytes were allocated",m->info-m->size); + else if( m->info == 0 ) + strcpy(reason_str,"factorization was successfull"); + else + sprintf(reason_str,"unknown exit code %d",m->info); + return reason_str; +} + +#endif +#endif diff --git a/Source/Solver/solver_superlu.h b/Source/Solver/solver_superlu.h new file mode 100644 index 0000000000000000000000000000000000000000..0df01ea808843f2978327eebbd3a6ba95e069e10 --- /dev/null +++ b/Source/Solver/solver_superlu.h @@ -0,0 +1,13 @@ +#ifndef _SOLVER_SUPERLU +#define _SOLVER_SUPERLU + +//IF CERTAIN FUNCTIONALITY IS NOT AVAILABLE, PLEASE THROW INMOST::NotImplemented EXCEPTION + +void SolverInitializeSuperLU(int * argc,char *** argv, const char * file_options); +void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values); +void SolverDestroyDataSuperLU(void ** data); +void SolverInitDataSuperLU(void ** data); +bool SolverSolveSuperLU(void * data, void * rhs_data, void * sol_data); +const char * SolverConvergedReasonSuperLU(void * data); + +#endif diff --git a/Tests/solver_test000/CMakeLists.txt b/Tests/solver_test000/CMakeLists.txt index f2868618ef54ffc87ebe908a2fa751e57b390f35..526be7f416ff55024712da7aff0e129140637c10 100644 --- a/Tests/solver_test000/CMakeLists.txt +++ b/Tests/solver_test000/CMakeLists.txt @@ -36,6 +36,10 @@ if(USE_SOLVER_MONDRIAAN) message("linking solver_test000 with Mondriaan") target_link_libraries(solver_test000 ${MONDRIAAN_LIBRARIES}) endif() +if(USE_SOLVER_SUPERLU) + message("linking solver_test000 with SuperLU") + target_link_libraries(solver_test000 ${SUPERLU_LIBRARIES}) +endif() add_test(NAME solver_test000_serial_inner_ilu2 COMMAND $ 0 0) diff --git a/Tests/solver_test001/CMakeLists.txt b/Tests/solver_test001/CMakeLists.txt index 903092c5440ac71e1bf4e19b272c352a0fb2f27f..5b89a6b4b5149f225b261b94f71b07cff5b1a255 100644 --- a/Tests/solver_test001/CMakeLists.txt +++ b/Tests/solver_test001/CMakeLists.txt @@ -34,7 +34,10 @@ if(USE_SOLVER_MONDRIAAN) message("linking solver_test001 with Mondriaan") target_link_libraries(solver_test001 ${MONDRIAAN_LIBRARIES}) endif() - +if(USE_SOLVER_SUPERLU) + message("linking solver_test001 with SuperLU") + target_link_libraries(solver_test001 ${SUPERLU_LIBRARIES}) +endif() set(TOKAMAK_TESTS utm300 utm1700a diff --git a/Tests/solver_test002/CMakeLists.txt b/Tests/solver_test002/CMakeLists.txt index 905767670250a9150777b2ddbf4500c8aeebbce1..26538f20d45c194d78b87f3657c12b01308bf145 100644 --- a/Tests/solver_test002/CMakeLists.txt +++ b/Tests/solver_test002/CMakeLists.txt @@ -37,6 +37,10 @@ if(USE_SOLVER) message("linking solver_test002 with Mondriaan") target_link_libraries(solver_test002 ${MONDRIAAN_LIBRARIES}) endif() + if(USE_SOLVER_SUPERLU) + message("linking solver_test002 with SuperLU") + target_link_libraries(solver_test002 ${SUPERLU_LIBRARIES}) + endif() endif() add_test(NAME solver_test002_serial_inner_ilu2 COMMAND $ 0 20) diff --git a/cmake_find/FindSuperLU.cmake b/cmake_find/FindSuperLU.cmake new file mode 100644 index 0000000000000000000000000000000000000000..94b8fa49aedff1ea671b242bbe8f15845bbe9bf4 --- /dev/null +++ b/cmake_find/FindSuperLU.cmake @@ -0,0 +1,25 @@ + +# Umfpack lib usually requires linking to a blas library. +# It is up to the user of this module to find a BLAS and link to it. + +if (SUPERLU_INCLUDES AND SUPERLU_LIBRARIES) + set(SUPERLU_FIND_QUIETLY TRUE) +endif (SUPERLU_INCLUDES AND SUPERLU_LIBRARIES) + + +find_path(SUPERLU_INCLUDES + NAMES + superlu/supermatrix.h + PATHS + $ENV{SUPERLUDIR}/include + ${SUPERLUDIR}/include + ${INCLUDE_INSTALL_DIR} +) + +find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR}/lib ${SUPERLUDIR}/lib ${LIB_INSTALL_DIR}) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(SUPERLU DEFAULT_MSG + SUPERLU_INCLUDES SUPERLU_LIBRARIES) + +mark_as_advanced(SUPERLU_INCLUDES SUPERLU_LIBRARIES) diff --git a/inmost-config.cmake.in b/inmost-config.cmake.in index 5ae985026bcde9e2dfd4119a829551c7e0d12324..12fab576c0e47801c767c2dfb5daa13166bdc044 100644 --- a/inmost-config.cmake.in +++ b/inmost-config.cmake.in @@ -76,6 +76,11 @@ if( USE_SOLVER_TRILINOS ) endif(MSVC) endif( USE_SOLVER_TRILINOS ) +if( USE_SOLVER_SUPERLU ) + list(APPEND INMOST_LIBRARIES "@SUPERLU_LIBRARIES") + list(APPEND INMOST_INCLUDE_DIRS "@SUPERLU_INCLUDES") +endif( USE_SOLVER_SUPERLU ) + if( USE_SOLVER_PETSC ) list(APPEND INMOST_LIBRARIES "@PETSC_LIBRARIES@")