Commit c9109b2d authored by Kirill Terekhov's avatar Kirill Terekhov

Features

Fixed NaN output into vtk files.

Added support for SuperLU (serial) solver.
parent ff8a089d
...@@ -34,6 +34,7 @@ option(USE_SOLVER_METIS "Use METIS for matrix reordering" OFF) ...@@ -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_MONDRIAAN "Use Mondriaan for matrix reordering" OFF)
option(USE_SOLVER_PETSC "Use PETSc solvers" OFF) option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
option(USE_SOLVER_TRILINOS "Use Trilinos 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_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
#option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF) #option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
#option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates 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) ...@@ -183,6 +184,22 @@ if(USE_SOLVER_TRILINOS)
endif() 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) if(USE_AUTODIFF_OPENCL)
find_package(OpenCL) find_package(OpenCL)
if(OPENCL_FOUND) if(OPENCL_FOUND)
......
...@@ -30,6 +30,10 @@ if(USE_SOLVER) ...@@ -30,6 +30,10 @@ if(USE_SOLVER)
message("linking ADFVDiscr with Mondriaan") message("linking ADFVDiscr with Mondriaan")
target_link_libraries(ADFVDiscr ${MONDRIAAN_LIBRARIES}) target_link_libraries(ADFVDiscr ${MONDRIAAN_LIBRARIES})
endif() endif()
if(USE_SOLVER_SUPERLU)
message("linking ADFVDiscr with SuperLU")
target_link_libraries(ADFVDiscr ${SUPERLU_LIBRARIES})
endif()
endif() endif()
......
...@@ -29,6 +29,10 @@ if(USE_SOLVER) ...@@ -29,6 +29,10 @@ if(USE_SOLVER)
message("linking ADMFD with Mondriaan") message("linking ADMFD with Mondriaan")
target_link_libraries(ADMFD ${MONDRIAAN_LIBRARIES}) target_link_libraries(ADMFD ${MONDRIAAN_LIBRARIES})
endif() endif()
if(USE_SOLVER_SUPERLU)
message("linking ADMFD with SuperLU")
target_link_libraries(ADMFD ${SUPERLU_LIBRARIES})
endif()
endif() endif()
......
...@@ -449,8 +449,8 @@ int main(int argc,char ** argv) ...@@ -449,8 +449,8 @@ int main(int argc,char ** argv)
//Change matrix by nonlinear correction for DMP //Change matrix by nonlinear correction for DMP
++total; ++total;
if( !isDMP ) //if( !isDMP )
//if( false ) if( false )
{ {
const real var = 0;//0.25; const real var = 0;//0.25;
vMatrix vnKGRAD = nKGRAD; vMatrix vnKGRAD = nKGRAD;
...@@ -557,11 +557,12 @@ int main(int argc,char ** argv) ...@@ -557,11 +557,12 @@ int main(int argc,char ** argv)
if( R.Norm() < 1.0e-4 ) break; if( R.Norm() < 1.0e-4 ) break;
Solver S(Solver::INNER_MPTILUC); Solver S(Solver::INNER_MPTILUC);
S.SetMatrix(R.GetJacobian()); //Solver S(Solver::SUPERLU);
S.SetParameterReal("relative_tolerance", 1.0e-14); S.SetParameterReal("relative_tolerance", 1.0e-14);
S.SetParameterReal("absolute_tolerance", 1.0e-12); S.SetParameterReal("absolute_tolerance", 1.0e-12);
S.SetParameterReal("drop_tolerance", 1.0e-2); S.SetParameterReal("drop_tolerance", 1.0e-2);
S.SetParameterReal("reuse_tolerance", 1.0e-4); S.SetParameterReal("reuse_tolerance", 1.0e-4);
S.SetMatrix(R.GetJacobian());
//std::fill(Update.Begin(),Update.End(),0.0); //std::fill(Update.Begin(),Update.End(),0.0);
if( S.Solve(R.GetResidual(),Update) ) if( S.Solve(R.GetResidual(),Update) )
{ {
......
...@@ -16,29 +16,7 @@ if(OPENGL_FOUND) ...@@ -16,29 +16,7 @@ if(OPENGL_FOUND)
set_target_properties(DrawMatrix PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") set_target_properties(DrawMatrix PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif() endif()
endif(USE_MPI) 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) install(TARGETS DrawMatrix EXPORT inmost-targets RUNTIME DESTINATION bin)
......
...@@ -29,6 +29,10 @@ if(USE_SOLVER) ...@@ -29,6 +29,10 @@ if(USE_SOLVER)
message("linking FVDiscr with Mondriaan") message("linking FVDiscr with Mondriaan")
target_link_libraries(FVDiscr ${MONDRIAAN_LIBRARIES}) target_link_libraries(FVDiscr ${MONDRIAAN_LIBRARIES})
endif() endif()
if(USE_SOLVER_SUPERLU)
message("linking FVDiscr with SuperLU")
target_link_libraries(FVDiscr ${SUPERLU_LIBRARIES})
endif()
endif() endif()
......
...@@ -37,6 +37,10 @@ if(USE_SOLVER) ...@@ -37,6 +37,10 @@ if(USE_SOLVER)
message("linking MatSolve with Mondriaan") message("linking MatSolve with Mondriaan")
target_link_libraries(MatSolve ${MONDRIAAN_LIBRARIES}) target_link_libraries(MatSolve ${MONDRIAAN_LIBRARIES})
endif() endif()
if(USE_SOLVER_SUPERLU)
message("linking MatSolve with SuperLU")
target_link_libraries(MatSolve ${SUPERLU_LIBRARIES})
endif()
endif() endif()
......
...@@ -29,6 +29,10 @@ if(USE_SOLVER) ...@@ -29,6 +29,10 @@ if(USE_SOLVER)
message("linking Solver with Mondriaan") message("linking Solver with Mondriaan")
target_link_libraries(Solver ${MONDRIAAN_LIBRARIES}) target_link_libraries(Solver ${MONDRIAAN_LIBRARIES})
endif() endif()
if(USE_SOLVER_SUPERLU)
message("linking Solver with SuperLU")
target_link_libraries(Solver ${SUPERLU_LIBRARIES})
endif()
endif(USE_SOLVER) endif(USE_SOLVER)
if(USE_PARTITIONER) if(USE_PARTITIONER)
......
...@@ -20,6 +20,7 @@ ...@@ -20,6 +20,7 @@
#cmakedefine USE_SOLVER_PETSC #cmakedefine USE_SOLVER_PETSC
#cmakedefine USE_SOLVER_TRILINOS #cmakedefine USE_SOLVER_TRILINOS
#cmakedefine USE_SOLVER_ANI #cmakedefine USE_SOLVER_ANI
#cmakedefine USE_SOLVER_SUPERLU
#cmakedefine USE_NONLINEAR #cmakedefine USE_NONLINEAR
#cmakedefine USE_NONLINEAR_TRILINOS #cmakedefine USE_NONLINEAR_TRILINOS
......
...@@ -54,7 +54,8 @@ namespace INMOST ...@@ -54,7 +54,8 @@ namespace INMOST
PETSc, ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/ 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/ ANI, ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/
FCBIILU2, ///< external FCBIILU2 Solver (BIILU2 parallel F2C version). 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); static std::string TypeName(Type t);
......
...@@ -273,7 +273,11 @@ safe_output: ...@@ -273,7 +273,11 @@ safe_output:
case DATA_REAL: case DATA_REAL:
{ {
Storage::real_array arr = it->RealArray(tags[i]); 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<double>(arr[m]);
fprintf(f,"%14e ",val != val ? -0.9999E30 : val);
}
fprintf(f,"\n"); fprintf(f,"\n");
} }
break; break;
...@@ -288,7 +292,11 @@ safe_output: ...@@ -288,7 +292,11 @@ safe_output:
case DATA_VARIABLE: case DATA_VARIABLE:
{ {
Storage::var_array arr = it->VariableArray(tags[i]); 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<double>(arr[m].GetValue());
fprintf(f,"%14e ",val != val ? -0.9999E30 : val);
}
fprintf(f,"\n"); fprintf(f,"\n");
} }
break; break;
...@@ -346,7 +354,11 @@ safe_output: ...@@ -346,7 +354,11 @@ safe_output:
case DATA_REAL: case DATA_REAL:
{ {
Storage::real_array arr = it->RealArray(tags[i]); 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<double>(arr[m]);
fprintf(f,"%14e ",(val != val ? -0.9999E30 : val));
}
fprintf(f,"\n"); fprintf(f,"\n");
} }
break; break;
...@@ -361,7 +373,11 @@ safe_output: ...@@ -361,7 +373,11 @@ safe_output:
case DATA_VARIABLE: case DATA_VARIABLE:
{ {
Storage::var_array arr = it->VariableArray(tags[i]); 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<double>(arr[m].GetValue());
fprintf(f,"%14e ",(val != val ? -0.9999E30 : val));
}
fprintf(f,"\n"); fprintf(f,"\n");
} }
break; break;
......
...@@ -5,12 +5,14 @@ set(SOURCE ...@@ -5,12 +5,14 @@ set(SOURCE
${CMAKE_CURRENT_SOURCE_DIR}/solver_ani.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_ani.cpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.cpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_petsc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_petsc.cpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_superlu.cpp
) )
set(HEADER set(HEADER
${HEADER} ${HEADER}
${CMAKE_CURRENT_SOURCE_DIR}/solver_prototypes.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_prototypes.hpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_petsc.h ${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_ilu2.hpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_ddpqiluc2.hpp
${CMAKE_CURRENT_SOURCE_DIR}/solver_bcgsl.hpp ${CMAKE_CURRENT_SOURCE_DIR}/solver_bcgsl.hpp
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#include "inmost_solver.h" #include "inmost_solver.h"
#if defined(USE_SOLVER) #if defined(USE_SOLVER)
#include "solver_petsc.h" #include "solver_petsc.h"
#include "solver_superlu.h"
#include "solver_ani.h" #include "solver_ani.h"
#include "solver_ilu2.hpp" #include "solver_ilu2.hpp"
#include "solver_ddpqiluc2.hpp" #include "solver_ddpqiluc2.hpp"
...@@ -43,11 +44,10 @@ ...@@ -43,11 +44,10 @@
#include "Teuchos_TestForException.hpp" #include "Teuchos_TestForException.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp" #include "Teuchos_XMLParameterListHelpers.hpp"
#endif #endif
//#define USE_OMP //#define USE_OMP
//#define KSOLVER BCGSL_solver #define KSOLVER BCGSL_solver
#define KSOLVER BCGS_solver //#define KSOLVER BCGS_solver
//#define ACCELERATED_CONDEST //#define ACCELERATED_CONDEST
//#define PRINT_CONDEST //#define PRINT_CONDEST
...@@ -61,6 +61,7 @@ namespace INMOST ...@@ -61,6 +61,7 @@ namespace INMOST
static std::string ani_database_file = ""; static std::string ani_database_file = "";
static std::string fcbiilu2_database_file = ""; static std::string fcbiilu2_database_file = "";
static std::string k3biilu2_database_file = ""; static std::string k3biilu2_database_file = "";
static std::string superlu_database_file = "";
...@@ -86,6 +87,7 @@ namespace INMOST ...@@ -86,6 +87,7 @@ namespace INMOST
case ANI: return "ANI"; case ANI: return "ANI";
case FCBIILU2: return "FCBIILU2"; case FCBIILU2: return "FCBIILU2";
case K3BIILU2: return "K3BIILU2"; case K3BIILU2: return "K3BIILU2";
case SUPERLU: return "SUPERLU";
} }
return "Unknown"; return "Unknown";
} }
...@@ -852,53 +854,58 @@ namespace INMOST ...@@ -852,53 +854,58 @@ namespace INMOST
(void)argv; (void)argv;
if( database != NULL ) if( database != NULL )
{ {
FILE * f = fopen(database,"r"); FILE * f = fopen(database,"r");
if( f ) if( f )
{ {
//std::fstream file(database,std::ios::in); //std::fstream file(database,std::ios::in);
char str[4096]; char str[4096];
//while( !file.eof() && file.getline(str,4096) ) //while( !file.eof() && file.getline(str,4096) )
while( !feof(f) && fgets(str,4096,f) ) while( !feof(f) && fgets(str,4096,f) )
{ {
int k = 0, l; int k = 0, l;
for(k = 0; k < (int)strlen(str); ++k) for(k = 0; k < (int)strlen(str); ++k)
{ {
if( str[k] == ':' ) break; if( str[k] == ':' ) break;
} }
if( k == strlen(str) ) continue; //invalid line if( k == strlen(str) ) continue; //invalid line
for(l = 0; l < k; ++l) str[l] = tolower(str[l]); for(l = 0; l < k; ++l) str[l] = tolower(str[l]);
l = (int)strlen(str)-1; // Right-trim string l = (int)strlen(str)-1; // Right-trim string
while(l > 0 && isspace(str[l]) ) --l; while(l > 0 && isspace(str[l]) ) --l;
str[l+1] = 0; str[l+1] = 0;
l = k+1; l = k+1;
while(l < (int)strlen(str) && isspace(str[l]) ) ++l; while(l < (int)strlen(str) && isspace(str[l]) ) ++l;
if( l == strlen(str) ) continue; //skip empty entry if( l == strlen(str) ) continue; //skip empty entry
if( !strncmp(str,"petsc",k) ) if( !strncmp(str,"petsc",k) )
petsc_database_file = std::string(str+l); petsc_database_file = std::string(str+l);
else if( !strncmp(str,"trilinos_ifpack",k) ) else if( !strncmp(str,"trilinos_ifpack",k) )
trilinos_ifpack_database_file = std::string(str+l); trilinos_ifpack_database_file = std::string(str+l);
else if( !strncmp(str,"trilinos_aztec",k) ) else if( !strncmp(str,"trilinos_aztec",k) )
trilinos_aztec_database_file = std::string(str+l); trilinos_aztec_database_file = std::string(str+l);
else if( !strncmp(str,"trilinos_ml",k) ) else if( !strncmp(str,"trilinos_ml",k) )
trilinos_ml_database_file = std::string(str+l); trilinos_ml_database_file = std::string(str+l);
else if( !strncmp(str,"trilinos_belos",k) ) else if( !strncmp(str,"trilinos_belos",k) )
trilinos_belos_database_file = std::string(str+l); trilinos_belos_database_file = std::string(str+l);
else if( !strncmp(str,"ani",k) ) else if( !strncmp(str,"ani",k) )
ani_database_file = std::string(str+l); ani_database_file = std::string(str+l);
else if( !strncmp(str,"fcbiilu2",k) ) else if( !strncmp(str,"fcbiilu2",k) )
fcbiilu2_database_file = std::string(str+l); fcbiilu2_database_file = std::string(str+l);
else if( !strncmp(str,"k3biilu2",k) ) else if( !strncmp(str,"k3biilu2",k) )
k3biilu2_database_file = std::string(str+l); k3biilu2_database_file = std::string(str+l);
} else if( !strncmp(str,"superlu",k) )
//file.close(); superlu_database_file = std::string(str+l);
fclose(f); }
} //file.close();
fclose(f);
}
} }
//std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl; //std::cout << "PETSc \"" << petsc_database_file << "\"" << std::endl;
//std::cout << "Trilinos_Ifpack \"" << trilinos_ifpack_database_file << "\"" << std::endl; //std::cout << "Trilinos_Ifpack \"" << trilinos_ifpack_database_file << "\"" << std::endl;
#if defined(USE_SOLVER_PETSC) #if defined(USE_SOLVER_PETSC)
SolverInitializePetsc(argc,argv,petsc_database_file.c_str()); SolverInitializePetsc(argc,argv,petsc_database_file.c_str());
#endif #endif
#if defined(USE_SOLVER_SUPERLU)
SolverInitializeSuperLU(argc,argv,superlu_database_file.c_str());
#endif
#if defined(USE_SOLVER_ANI) #if defined(USE_SOLVER_ANI)
SolverInitializeAni(argc,argv,ani_database_file.c_str()); SolverInitializeAni(argc,argv,ani_database_file.c_str());
#endif #endif
...@@ -922,9 +929,12 @@ namespace INMOST ...@@ -922,9 +929,12 @@ namespace INMOST
} }
} }
} }
#endif
#if defined(USE_SOLVER_SUPERLU)
#endif #endif
Sparse::CreateRowEntryType(); Sparse::CreateRowEntryType();
is_initialized = true; is_initialized = true;
} }
void Solver::Finalize() void Solver::Finalize()
...@@ -1070,6 +1080,12 @@ namespace INMOST ...@@ -1070,6 +1080,12 @@ namespace INMOST
SolverInitDataK3biilu2(&solver_data,_comm,name.c_str()); SolverInitDataK3biilu2(&solver_data,_comm,name.c_str());
} }
#endif #endif
#if defined(USE_SOLVER_SUPERLU)
if( _pack == SUPERLU )
{
SolverInitDataSuperLU(&solver_data);
}
#endif
#if defined(USE_SOLVER_TRILINOS) #if defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -1271,6 +1287,13 @@ namespace INMOST ...@@ -1271,6 +1287,13 @@ namespace INMOST
solver_data = NULL; solver_data = NULL;
} }
} }
#endif
#if defined(USE_SOLVER_SUPERLU)
if(_pack == SUPERLU )
{
SolverDestroyDataSuperLU(&solver_data);
matrix_data = NULL;
}
#endif #endif
if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2) if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2)
{ {
...@@ -1721,6 +1744,40 @@ namespace INMOST ...@@ -1721,6 +1744,40 @@ namespace INMOST
ok = true; ok = true;
} }
#endif #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 defined(USE_SOLVER_TRILINOS)
if( _pack == Trilinos_Aztec || if( _pack == Trilinos_Aztec ||
_pack == Trilinos_Belos || _pack == Trilinos_Belos ||
...@@ -2493,6 +2550,16 @@ namespace INMOST ...@@ -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 #endif
if(_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2) if(_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC || _pack == INNER_MPTILU2)
{ {
......
...@@ -2,8 +2,6 @@ ...@@ -2,8 +2,6 @@
#define SOLVER_ANI_H_INCLUDED #define SOLVER_ANI_H_INCLUDED
#include "inmost_solver.h" #include "inmost_solver.h"
#if defined(USE_SOLVER_ANI)
void MatrixCopyDataAni(void ** ppA, void * pB); void MatrixCopyDataAni(void ** ppA, void * pB);
void MatrixAssignDataAni(void * pA, void * pB); void MatrixAssignDataAni(void * pA, void * pB);
void MatrixInitDataAni(void ** ppA, INMOST_MPI_Comm comm, const char * name); 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 ...@@ -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); bool SolverSolveAni(void * data, void * rhs_data, void * sol_data);
int SolverIterationNumberAni(void * data); int SolverIterationNumberAni(void * data);
double SolverResidualNormAni(void * data); double SolverResidualNormAni(void * data);
#endif //USE_SOLVER_ANI
#endif //SOLVER_ANI_H_INCLUDED #endif //SOLVER_ANI_H_INCLUDED
#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<SuperLU_data *>(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<SuperLU_data *>(*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<void *>(new SuperLU_data);
SuperLU_data * m = static_cast<SuperLU_data *>(*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<SuperLU_data *>(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<SuperLU_data *>(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 )