Commit 73291432 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Add support for SuperLU_dist, automatic download and install via cmake; other minor updates

parent 0524b955
......@@ -309,7 +309,6 @@ if(USE_SOLVER_SUPERLU)
endif( WIN32 )
message("-- SUPERLU_INCLUDE_DIR: ${SUPERLU_INCLUDES}")
message("-- SUPERLU_LIBRARIES: ${SUPERLU_LIBRARIES}")
include_directories(${SUPERLU_INCLUDES})
else(USE_SOLVER_SUPERLU_DOWNLOAD)
set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
endif(USE_SOLVER_SUPERLU_DOWNLOAD)
......@@ -317,9 +316,91 @@ if(USE_SOLVER_SUPERLU)
MESSAGE("\nFound SuperLU! Here are the details: ")
MESSAGE("INCLUDES: ${SUPERLU_INCLUDES}")
MESSAGE("LIBRARIES: ${SUPERLU_LIBRARIES}")
include_directories(${SUPERLU_INCLUDES})
endif()
find_package(SUPERLU_DIST)
if(NOT SUPERLU_DIST_FOUND)
#set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
MESSAGE("SUPERLU_DIST NOT FOUND")
if( USE_SOLVER_SUPERLU_DOWNLOAD )
if( USE_MPI AND USE_PARTITIONER_PARMETIS )
message("-- Download SuperLU_DIST")
set(C_STANDARD_FLAG "-std=c99")
string (REPLACE ";" "|" PASS_PARMETIS_LIBRARIES "${PARMETIS_LIBRARIES}")
message("-- Pass parmetis libraries ${PASS_PARMETIS_LIBRARIES} ")
#do we need blas?
find_package(LAPACK)
if( NOT LAPACK_FOUND )
ExternalProject_Add(LAPACK
URL "https://www.netlib.org/clapack/clapack-3.2.1-CMAKE.tgz"
UPDATE_DISCONNECTED 1
PREFIX "${LIB_DOWNLOAD_PATH}")
if( WIN32 )
set(LAPACK_LIBRARIES "${LIB_DOWNLOAD_PATH}/lib/lapack.lib" CACHE PATH "LAPACK_LIBRARIES is set" FORCE)
else( WIN32 )
set(LAPACK_LIBRARIES "${LIB_DOWNLOAD_PATH}/lib/liblapack.a" CACHE PATH "LAPACK_LIBRARIES is set" FORCE)
endif( WIN32 )
set(LAPACK_LIBRARIES "${SUPERLU_LIBRARIES};${LAPACK_LIBRARIES};${PARMETIS_LIBRARIES};${LIB_DOWNLOAD_PATH}/lib/superlu_dist.lib" CACHE PATH "LAPACK_LIBRARIES is set" FORCE)
endif()
message("-- Pass lapack libraries ${LAPACK_LIBRARIES} ")
ExternalProject_Add(SuperLU_DIST
GIT_REPOSITORY "https://github.com/xiaoyeli/superlu_dist.git"
GIT_TAG "master"
UPDATE_DISCONNECTED 1
PREFIX "${LIB_DOWNLOAD_PATH}"
#this command is responsible for update of repository from git
#UPDATE_COMMAND ""
LIST_SEPARATOR | # Use the alternate list separator
CMAKE_ARGS "-DTPL_PARMETIS_INCLUDE_DIRS=${PARMETIS_INCLUDE_DIR} "
"-DTPL_PARMETIS_LIBRARIES=${PASS_PARMETIS_LIBRARIES} "
"-DTPL_LAPACK_LIBRARIES=${LAPACK_LIBRARIES}"
"-Denable_examples=OFF "
"-Denable_tests=OFF "
"-Denable_doc=OFF "
"-DCMAKE_INSTALL_PREFIX=${LIB_DOWNLOAD_PATH} "
"-DUSE_XSDK_DEFAULTS=TRUE "
"-DXSDK_ENABLE_Fortran=FALSE "
"-DNOFORTRAN=TRUE "
"-DCMAKE_C_FLAGS='${C_STANDARD_FLAG} -DMETIS_EXPORT=' "
"-DCMAKE_CXX_FLAGS='-DMETIS_EXPORT=' "
"-DCMAKE_BUILD_TYPE=Release "
"-DMPI_C_COMPILER=${MPI_C_COMPILER} "
"-DMPI_CXX_COMPILER=${MPI_CXX_COMPILER} "
"-DCMAKE_DISABLE_FIND_PACKAGE_OpenMP=TRUE "
"-DCMAKE_POSITION_INDEPENDENT_CODE=ON "
)
message("-- Linking SuperLU_dist and ParMETIS")
add_dependencies(SuperLU_DIST ParMETIS)
add_definitions(-DUSE_SOLVER_SUPERLU_DIST)
#same as above
#set(SUPERLU_INCLUDES "${LIB_DOWNLOAD_PATH}/include" CACHE PATH "SUPERLU_INCLUDE_DIR is set" FORCE)
#extended
#TODO: properly install lapack
if( WIN32 )
set(SUPERLU_LIBRARIES "${SUPERLU_LIBRARIES};${LIB_DOWNLOAD_PATH}/lib/superlu_dist.lib;${LAPACK_LIBRARIES};${PARMETIS_LIBRARIES}" CACHE PATH "SUPERLU_LIBRARIES is set" FORCE)
else( WIN32 )
set(SUPERLU_LIBRARIES "${SUPERLU_LIBRARIES};${METIS_LIBRARIES};${LIB_DOWNLOAD_PATH}/lib/libsuperlu_dist.a;${LAPACK_LIBRARIES};${PARMETIS_LIBRARIES}" CACHE PATH "SUPERLU_LIBRARIES is set" FORCE)
endif( WIN32 )
message("-- SUPERLU_INCLUDE_DIR: ${SUPERLU_INCLUDES}")
message("-- SUPERLU_LIBRARIES: ${SUPERLU_LIBRARIES}")
else()
message("-- Skip SuperLU_DIST as either USE_MPI or USE_PARTITIONER_PARMETIS are not checked for compilation")
endif()
else(USE_SOLVER_SUPERLU_DOWNLOAD)
set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
endif(USE_SOLVER_SUPERLU_DOWNLOAD)
else()
MESSAGE("\nFound SuperLU_DIST! Here are the details: ")
set(SUPERLU_INCLUDES "${SUPERLU_INCLUDES};${SUPERLU_DIST_INCLUDES}")
set(SUPERLU_LIBRARIES "${SUPERLU_LIBRARIES};${SUPERLU_DIST_LIBRARIES}")
MESSAGE("INCLUDES: ${SUPERLU_INCLUDES}")
MESSAGE("LIBRARIES: ${SUPERLU_LIBRARIES}")
endif()
include_directories(${SUPERLU_INCLUDES})
endif()
if(USE_OPENCL)
......
......@@ -36,7 +36,7 @@ int main(int argc, char ** argv)
if( argc < 2 )
{
std::cout << "Usage: " << argv[0] << " mesh [axis=2 (0:x,1:y,2:z)] [control=0 (0:shear rate,1:pressure drop,2:maximal velocity)] [control_value=10] [visc=1.0e-5] [mesh_out=grid_out.pmf] [addsym=0]" << std::endl;
std::cout << "Usage: " << argv[0] << " mesh [axis=2 (0:x,1:y,2:z)] [control=0 (0:shear rate,1:pressure drop,2:maximal velocity)] [control_value=10] [visc=1.0e-5] [Diam=-1(autodetect)] [mesh_out=grid_out.pmf] [addsym=0]" << std::endl;
return 0;
}
......@@ -58,7 +58,7 @@ int main(int argc, char ** argv)
int control = 0;
double control_value = 10;
double Len = 8;
double Diam = 1;
double Diam = -1;
double shear = 10;//0; //shear rate
double dp = 0;//36e-6;
double vmax = 0;
......@@ -68,8 +68,9 @@ int main(int argc, char ** argv)
if( argc > 3 ) control = atoi(argv[3]);
if( argc > 4 ) control_value = atof(argv[4]);
if( argc > 5 ) visc = atof(argv[5]);
if( argc > 6 ) fout = std::string(argv[6]);
if( argc > 7 ) addsym = atof(argv[7]);
if( argc > 6 ) Diam = atof(argv[6]);
if( argc > 7 ) fout = std::string(argv[7]);
if( argc > 8 ) addsym = atof(argv[8]);
......@@ -119,7 +120,41 @@ int main(int argc, char ** argv)
}
}
Len = cmax[axis] - cmin[axis];
Diam = std::max(cmax[(axis+1)%3]-cmin[(axis+1)%3],cmax[(axis+2)%3]-cmin[(axis+2)%3]);
if( Diam == -1 ) Diam = std::max(cmax[(axis+1)%3]-cmin[(axis+1)%3],cmax[(axis+2)%3]-cmin[(axis+2)%3]);
if( true ) //project nodes to diameter
{
MarkerType bnode = m->CreateMarker();
double n[3];
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
it->UnitNormal(n);
if( fabs(n[axis]) < 1.0e-2 ) it->getNodes().SetMarker(bnode);
}
int ix = (axis+1)%3, iy = (axis+2)%3;
//double cx = 0.5*(cmax[ix]+cmin[ix]), cy = 0.5*(cmax[iy]+cmin[iy]), R = 0.5*Diam;
double cx = 0.5, cy = 0.5, R = 0.5*Diam;
std::cout << "center of cylinder " << cx << " " << cy << std::endl;
double tot_corr = 0;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) if( it->GetMarker(bnode) )
{
double & x = it->Coords()[ix];
double & y = it->Coords()[iy];
double r = sqrt(pow(x-cx,2)+pow(y-cy,2));
if( r )
{
double dx = (R/r-1.0)*(x-cx);
double dy = (R/r-1.0)*(y-cy);
//std::cout << "at " << x << "," << y << " r " << r << " dx " << dx << " dy " << dy << std::endl;
x += dx;
y += dy;
tot_corr += sqrt(dx*dx+dy*dy);
}
else std::cout << "node at center of cylinder: " << x << "," << y << std::endl;
}
m->ReleaseMarker(bnode,NODE);
std::cout << "projection done tot_corr " << tot_corr << std::endl;
}
if( shear )
......
......@@ -50,7 +50,7 @@ int main(int argc, char ** argv)
std::string fout = "grid_out.pmf";
double Umax = 2.25;
bool fix_cylinder = false;
int fix_cylinder = 0;
if( argc > 2 ) fout = std::string(argv[2]);
if( argc > 3 ) Umax = atof(argv[3]);
if( argc > 4 ) fix_cylinder = atoi(argv[4]);
......
......@@ -74,7 +74,7 @@ void rotate_tensor(double nrm[3], const rMatrix & Kin, rMatrix & Kout, double ph
double qmat[16], qrmat[16], sclmat[16], qnrm;
struct quat q, qr;
rMatrix N = rMatrix::FromVector(nrm,3);
rMatrix NX = rMatrix::CrossProduct(nrm);
rMatrix NX = rMatrix::CrossProductMatrix(nrm);
rMatrix U(3,1), NXU;
U(0,0) = 0;
U(1,0) = 0;
......
......@@ -35,7 +35,7 @@ namespace INMOST
AbstractEntry::Access<hessian_variable>(const Storage& e) const {return Unknown(e);}
#if defined(USE_MESH)
#if defined(USE_MESH) //Automatizator class does not exist without mesh
Automatizator * Automatizator::CurrentAutomatizator = NULL;
bool print_ad_ctor = false;
bool GetAutodiffPrint() {return print_ad_ctor;}
......
......@@ -1640,7 +1640,7 @@ namespace INMOST
/// and vector. For a x b equivalent is CrossProduct(a)*b.
/// @param vec Array of elements representing a vector.
/// @return A matrix representing cross product.
static Matrix<Var, pool_array_t<Var> > CrossProduct(const Var vec[3])
static Matrix<Var, pool_array_t<Var> > CrossProductMatrix(const Var vec[3])
{
// | 0 -z y |
// | z 0 -x |
......
......@@ -58,8 +58,8 @@ namespace INMOST
shell_expression() {}// if( GetAutodiffPrint() ) std::cout << this << " Shell Created for " << dynamic_cast<basic_expression *>(this) << std::endl;}
shell_expression(const shell_expression & other) {}//std::cout << this << " Shell Created from " << &other << std::endl;}
__INLINE virtual INMOST_DATA_REAL_TYPE GetValue() const {return static_cast<const Derived *>(this)->GetValue(); }
__INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const { return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
__INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const { return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
__INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const { if( mult ) return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
__INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const { if( mult ) return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
__INLINE virtual void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {return static_cast<const Derived *>(this)->GetHessian(multJ,J,multH,H); }
operator Derived & () {return *static_cast<Derived *>(this);}
operator const Derived & () const {return *static_cast<const Derived *>(this);}
......@@ -99,8 +99,8 @@ namespace INMOST
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
__INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
__INLINE INMOST_DATA_ENUM_TYPE GetIndex() const { return index; }
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const {if( index != ENUMUNDEF ) r[index] += mult;}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const {if( index != ENUMUNDEF ) r[index] += mult;}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const {if( mult && index != ENUMUNDEF ) r[index] += mult;}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const {if( mult && index != ENUMUNDEF ) r[index] += mult;}
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {if( index != ENUMUNDEF ) J.Push(index,multJ); (void)multH; (void)H;}
__INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE i) const {return index == i? 1.0 : 0.0;}
__INLINE var_expression & operator =(var_expression const & other)
......@@ -157,8 +157,9 @@ namespace INMOST
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
r[it->first] += it->second*mult;
r.AddRow(mult, entries);
//~ for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
//~ r[it->first] += it->second*mult;
}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
......@@ -424,8 +425,9 @@ namespace INMOST
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
r[it->first] += it->second*mult;
r.AddRow(mult,entries);
//for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
// r[it->first] += it->second*mult;
}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
......@@ -635,8 +637,9 @@ namespace INMOST
{
if( entries )
{
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
r[it->first] += it->second*mult;
r.AddRow(mult,*entries);
//~ for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
//~ r[it->first] += it->second*mult;
}
}
/// Retrive derivatives with multiplier into Sparse::Row structure.
......@@ -821,8 +824,9 @@ namespace INMOST
/// Retrive derivatives with multiplier into Sparse::RowMerger structure.
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
r[it->first] += it->second*mult;
r.AddRow(mult,*entries);
//~ for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
//~ r[it->first] += it->second*mult;
}
/// Retrive derivatives with multiplier into Sparse::Row structure.
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
......
......@@ -6,6 +6,7 @@
#include "solver_prototypes.hpp"
#include "solver_bcgsl.hpp"
#define KSOLVER BCGS_solver
//#define KSOLVER BCGSL_solver
//#define ACCELERATED_CONDEST
namespace INMOST {
......
......@@ -19,7 +19,7 @@ using namespace INMOST;
#define REORDER_RCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
#define REORDER_METIS_ND
//~ #define REORDER_METIS_ND
#endif
#if defined(USE_SOLVER_MONDRIAAN)
//#define REORDER_MONDRIAAN
......@@ -1129,7 +1129,7 @@ static bool allow_pivot = true;
i = wend;
#elif defined(REORDER_RCM)
{
if( verbosity )
if( verbosity > 1 )
printf("Reordering with RCM\n");
//create a symmetric graph of the matrix A + A^T
std::vector<INMOST_DATA_ENUM_TYPE> xadj(wend-wbeg+1), adjncy;
......
......@@ -8,8 +8,8 @@
//#define LFILL //control, that factorization is not less then fill for ilu2
//select one of the two rescaling techniques
#define RESCALE_EQUALIZE_1NORM //equalize 1-norms of each row and each column to 1
//#define RESCALE_EQUALIZE_2NORM
//#define RESCALE_EQUALIZE_1NORM //equalize 1-norms of each row and each column to 1
#define RESCALE_EQUALIZE_2NORM
//#define RESCALE_MAXIMUM_TRANSVERSAL //use rescaling produced by maximum product transversal algorithm that bounds all values between -1 and 1
#define REORDER_MAXIMUM_TRANSVERSAL
......
#include "SolverSUPERLU.h"
namespace INMOST {
namespace INMOST
{
SolverSUPERLU::SolverSUPERLU() {
SolverSUPERLU::SolverSUPERLU()
{
#if defined(USE_SOLVER_SUPERLU_DIST)
set_default_options_dist(&options_);
options_.PrintStat = NO;
options_.IterRefine = SLU_DOUBLE;
options_.ReplaceTinyPivot = YES;
options_.Equil = YES;
options_.Fact = DOFACT;
PStatInit(&stat_);
#else // USE_SOLVER_SUPERLU_DIST
set_default_options(&options);
StatInit(&stat);
perm_c = NULL;
perm_r = NULL;
#endif //USE_SOLVER_SUPERLU_DIST
a_size = 0;
}
SolverInterface *SolverSUPERLU::Copy(const SolverInterface *other) {
SolverInterface *SolverSUPERLU::Copy(const SolverInterface *other)
{
throw INMOST::SolverUnsupportedOperation;
}
void SolverSUPERLU::Assign(const SolverInterface *other) {
void SolverSUPERLU::Assign(const SolverInterface *other)
{
throw INMOST::SolverUnsupportedOperation;
}
void SolverSUPERLU::Setup(int *argc, char ***argv, SolverParameters &p) {
void SolverSUPERLU::Setup(int *argc, char ***argv, SolverParameters &p)
{
//read options from file and arguments
}
void SolverSUPERLU::SetMatrix(Sparse::Matrix &A, bool ModifiedPattern, bool OldPreconditioner) {
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();
int mend = A.GetLastIndex();
int size = 0;
remap = new int[mend - mbeg];
#if defined(USE_SOLVER_SUPERLU_DIST)
int mpisize = 1;
MPI_Comm_size(GetCommunicator(),&mpisize);
nnz = 0;
size = mend-mbeg;
for (int k = 0; k < size; ++k)
nnz += A[k + mbeg].Size();
ia = (int *) malloc(sizeof(int) * (size + 1));
ja = (int *) malloc(sizeof(int) * nnz);
a = (double *) malloc(sizeof(double) * nnz);
int q = 0, f = 0;
ia[0] = 0;
for (int k = 0; k < size; ++k)
{
Sparse::Row &r = A[k + mbeg];
for (INMOST_DATA_ENUM_TYPE l = 0; l < r.Size(); ++l)
{
ja[q] = r.GetIndex(l);
a[q] = r.GetValue(l);
++q;
}
ia[f + 1] = q;
f++;
}
int dims[2] = {0,0};
MPI_Dims_create(mpisize,2,dims);
//~ std::cout << "Size: " << mpisize << " dims: " << dims[0] << "," << dims[1] << std::endl;
superlu_gridinit(GetCommunicator(),dims[0],dims[1],&grid);
int offset = 0;
int local_size = size;
int global_size = size;
MPI_Exscan(&local_size,&offset,1,MPI_INT,MPI_SUM,GetCommunicator());
MPI_Allreduce(&local_size,&global_size,1,MPI_INT,MPI_SUM,GetCommunicator());
//std::cout << "M = " << m << " N = " << n << " local offset " << offset << std::endl;
//~ dCreate_CompRowLoc_Matrix_dist(&A_,m, n, pCSRMatrixFormat->nonzeros,local_size[0],offset,pCSRMatrixFormat->pData,pCSRMatrixFormat->pCols,pCSRMatrixFormat->pRows,SLU_NR_loc, SLU_D, SLU_GE);
dCreate_CompRowLoc_Matrix_dist(&(this->A),global_size, global_size, nnz,local_size,offset,a,ja,ia,SLU_NR_loc, SLU_D, SLU_GE);
ScalePermstructInit(global_size,global_size,&ScalePermstruct);
LUstructInit(global_size,&LUstruct);
a_size = size;
#else //USE_SOLVER_SUPERLU_DIST
remap = new int[mend - mbeg];
for (int k = 0; k < mend - mbeg; ++k) {
Sparse::Row &r = A[k + mbeg];
if (r.Size()) {
......@@ -70,50 +127,83 @@ namespace INMOST {
f++;
}
}
dCreate_CompCol_Matrix(&(this->A), size, size, nnz, a, ja, ia, SLU_NR, SLU_D, SLU_GE);
a_size = size;
dCreate_CompCol_Matrix(&(this->A), size, size, nnz, a, ja, ia, SLU_NR, SLU_D, SLU_GE);
perm_c = new int[size];
perm_r = new int[size];
#endif //USE_SOLVER_SUPERLU_DIST
}
bool SolverSUPERLU::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL) {
double *inout = new double[a_size];
int mbeg = RHS.GetFirstIndex(), mend = RHS.GetLastIndex();
for (int k = 0; k < mend - mbeg; ++k) {
if (remap[k] != -1) {
inout[remap[k]] = RHS[k + mbeg];
}
}
SuperMatrix B;
dCreate_Dense_Matrix(&B, a_size, 1, inout, a_size, SLU_DN, SLU_D, SLU_GE);
dgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);
Destroy_SuperMatrix_Store(&B);
bool ret = (info == 0);
for (int k = 0; k < mend - mbeg; ++k) {
if (remap[k] != -1) {
SOL[k + mbeg] = inout[remap[k]];
}
}
delete[] inout;
bool SolverSUPERLU::Solve(INMOST::Sparse::Vector &RHS, INMOST::Sparse::Vector &SOL)
{
bool ret = false;
#if defined(USE_SOLVER_SUPERLU_DIST)
//~ std::cout << __FILE__ << ":" << __LINE__ << " AStype " << (A.Stype == SLU_NR_loc) << " ADtype " << (A.Dtype == SLU_D? 1 : 0) << " AGEtype " << (A.Mtype == SLU_GE? 1 : 0) << std::endl;
double pBerr_[2];
int mbeg = RHS.GetFirstIndex(), mend = RHS.GetLastIndex();
SOL.SetInterval(mbeg,mend);
memcpy(&SOL[mbeg],&RHS[mbeg], (mend-mbeg)*sizeof(double));
//set_default_options_dist(&options_);
pdgssvx(&options_, &A, &ScalePermstruct, &SOL[mbeg], a_size, 1, &grid, &LUstruct, &SOLVEstruct, pBerr_, &stat_, &info);
//~ PStatPrint(&options_, &stat_, &grid);
options_.Fact = FACTORED;
#else //USE_SOLVER_SUPERLU_DIST
//~ std::cout << __FILE__ << ":" << __LINE__ << std::endl;
double *inout = new double[a_size];
int mbeg = RHS.GetFirstIndex(), mend = RHS.GetLastIndex();
for (int k = 0; k < mend - mbeg; ++k)
if (remap[k] != -1)
inout[remap[k]] = RHS[k + mbeg];
SuperMatrix B;
dCreate_Dense_Matrix(&B, a_size, 1, inout, a_size, SLU_DN, SLU_D, SLU_GE);
dgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);
Destroy_SuperMatrix_Store(&B);
for (int k = 0; k < mend - mbeg; ++k)
if (remap[k] != -1)
SOL[k + mbeg] = inout[remap[k]];
delete[] inout;
#endif //USE_SOLVER_SUPERLU_DIST
ret = (info == 0);
return ret;
}
bool SolverSUPERLU::Clear() {
Destroy_CompCol_Matrix(&A);
Destroy_SuperNode_Matrix(&L);
Destroy_CompCol_Matrix(&U);
StatFree(&stat);
if (perm_c != NULL) delete[] perm_c;
if (perm_r != NULL) delete[] perm_r;
if (remap != NULL) delete[] remap; //allocated outside
bool SolverSUPERLU::Clear()
{
//~ std::cout << "call clear!" << std::endl;
#if defined(USE_SOLVER_SUPERLU_DIST)
PStatFree(&stat_);
Destroy_CompRowLoc_Matrix_dist(&A);
ScalePermstructFree(&ScalePermstruct);
Destroy_LU(a_size,&grid,&LUstruct);
LUstructFree(&LUstruct);
if( options_.SolveInitialized )
dSolveFinalize(&options_,&SOLVEstruct);
superlu_gridexit(&grid);
#else //USE_SOLVER_SUPERLU_DIST
Destroy_CompCol_Matrix(&A);
Destroy_SuperNode_Matrix(&L);
Destroy_CompCol_Matrix(&U);
StatFree(&stat);
if (perm_c != NULL) delete[] perm_c;
if (perm_r != NULL) delete[] perm_r;
if (remap != NULL) delete[] remap; //allocated outside
#endif
return true;
}
bool SolverSUPERLU::isMatrixSet() {
bool SolverSUPERLU::isMatrixSet()
{
return a_size != 0;
}
std::string SolverSUPERLU::GetParameter(std::string name) const {
std::string SolverSUPERLU::GetParameter(std::string name) const
{
#if !defined(SILENCE_SET_PARAMETER)
std::cout << "SolverSUPERLU::GetParameter unsupported operation" << std::endl;
#endif
......@@ -121,46 +211,52 @@ namespace INMOST {
return "";
}
void SolverSUPERLU::SetParameter(std::string name, std::string value) {
void SolverSUPERLU::SetParameter(std::string name, std::string value)
{
#if !defined(SILENCE_SET_PARAMETER)
std::cout << "SolverSUPERLU::SetParameter unsupported operation" << std::endl;
#endif
//throw INMOST::SolverUnsupportedOperation;
}
INMOST_DATA_ENUM_TYPE SolverSUPERLU::Iterations() const {
INMOST_DATA_ENUM_TYPE SolverSUPERLU::Iterations() const
{
return 1;
}
INMOST_DATA_REAL_TYPE SolverSUPERLU::Residual() const {
INMOST_DATA_REAL_TYPE SolverSUPERLU::Residual() const
{
return 0;
}
const std::string SolverSUPERLU::ReturnReason() const {
const std::string SolverSUPERLU::ReturnReason() const
{
char reason_str[256];
if (info <= a_size) {
if (info <= a_size)
sprintf(reason_str, "diagonal entry of U-factor is exactly singular at %d/%d", info, a_size);
} else if (info > a_size) {
else if (info > a_size)
sprintf(reason_str, "memory allocation failed after %d bytes were allocated", info - a_size);
} else if (info == 0) {
else if (info == 0)
strcpy(reason_str, "factorization was successfull");
} else {
else
sprintf(reason_str, "unknown exit code %d", info);
}
return std::string(reason_str);
}
const std::string SolverSUPERLU::SolverName() const {
const std::string SolverSUPERLU::SolverName() const
{
return "superlu";
}