Skip to content
Snippets Groups Projects
Commit 203ff04d authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

prepare petsc

parent efa4f073
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 158 additions and 136 deletions
CMakeLists.txt
+ 34
− 6
View file @ 203ff04d
......@@ -36,6 +36,7 @@ option(ENABLE_OCC "Enable OpenCASCADE geometrical models" ON)
option(ENABLE_OOFELIE "Set compiler options to match oofelie settings" OFF)
option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
option(ENABLE_PARSER "Build the GEO file parser" ON)
option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ON)
option(ENABLE_POST "Build the post-processing module" ON)
option(ENABLE_QT "Build QT GUI" OFF)
option(ENABLE_TAUCS "Enable Taucs linear algebra solver" ON)
......@@ -120,11 +121,11 @@ macro(append_gmsh_src DIRNAME FILES)
set(GMSH_DIRS ${GMSH_DIRS};${DIRNAME} PARENT_SCOPE)
endmacro(append_gmsh_src)
macro(find_all_libraries VARNAME LISTNAME SUFFIX)
macro(find_all_libraries VARNAME LISTNAME PATH SUFFIX)
set(${VARNAME})
list(LENGTH ${LISTNAME} NUM_LIST)
foreach(LIB ${${LISTNAME}})
find_library(FOUND_LIB ${LIB} PATH_SUFFIXES ${SUFFIX})
find_library(FOUND_LIB ${LIB} PATHS ${PATH} PATH_SUFFIXES ${SUFFIX})
if(FOUND_LIB)
list(APPEND ${VARNAME} ${FOUND_LIB})
endif(FOUND_LIB)
......@@ -161,14 +162,14 @@ if(ENABLE_BLAS_LAPACK)
set(MKL_PATH ia32/lib)
endif(HAVE_64BIT_SIZE_T)
set(MKL_LIBS_REQUIRED libguide40 mkl_intel_c mkl_intel_thread mkl_core)
find_all_libraries(LAPACK_LIBRARIES MKL_LIBS_REQUIRED ${MKL_PATH})
find_all_libraries(LAPACK_LIBRARIES MKL_LIBS_REQUIRED "" ${MKL_PATH})
if(LAPACK_LIBRARIES)
set(HAVE_BLAS TRUE)
set(HAVE_LAPACK TRUE)
list(APPEND CONFIG_OPTIONS "IntelMKL")
else(LAPACK_LIBRARIES)
set(REFLAPACK_LIBS_REQUIRED lapack blas g2c gcc)
find_all_libraries(LAPACK_LIBRARIES REFLAPACK_LIBS_REQUIRED "")
find_all_libraries(LAPACK_LIBRARIES REFLAPACK_LIBS_REQUIRED "" "")
if(LAPACK_LIBRARIES)
set(HAVE_BLAS TRUE)
set(HAVE_LAPACK TRUE)
......@@ -184,7 +185,7 @@ if(ENABLE_BLAS_LAPACK)
set(MKL_PATH lib/32)
endif(HAVE_64BIT_SIZE_T)
set(MKL_LIBS_REQUIRED mkl_gf_lp64 iomp5 mkl_gnu_thread mkl_core)
find_all_libraries(LAPACK_LIBRARIES MKL_LIBS_REQUIRED ${MKL_PATH})
find_all_libraries(LAPACK_LIBRARIES MKL_LIBS_REQUIRED "" ${MKL_PATH})
if(LAPACK_LIBRARIES)
set(HAVE_BLAS TRUE)
set(HAVE_LAPACK TRUE)
......@@ -498,6 +499,33 @@ if(ENABLE_TAUCS)
endif(TAUCS_LIB)
endif(ENABLE_TAUCS)
if(ENABLE_PETSC)
set(ENV_PETSC_DIR $ENV{PETSC_DIR})
set(ENV_PETSC_ARCH $ENV{PETSC_ARCH})
if(EXISTS ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/conf/petscvariables)
message(STATUS "Found PETSc 3")
set(HAVE_PETSC TRUE)
list(APPEND CONFIG_OPTIONS "PETSc")
file(STRINGS ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/conf/petscvariables VARS NEWLINE_CONSUME)
# append include directories
list(APPEND EXTERNAL_INCLUDES ${ENV_PETSC_DIR}/include)
list(APPEND EXTERNAL_INCLUDES ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/include)
string(REGEX MATCH "PACKAGES_INCLUDES = [-_/\\.a-zA-Z0-9 ]*" PETSC_PACKAGES_INC ${VARS})
string(REPLACE "PACKAGES_INCLUDES = " "" PETSC_PACKAGES_INC ${PETSC_PACKAGES_INC})
string(REPLACE "-I" "" PETSC_PACKAGES_INC ${PETSC_PACKAGES_INC})
string(REPLACE " " ";" PETSC_PACKAGES_INC ${PETSC_PACKAGES_INC})
foreach(VAR ${PETSC_PACKAGES_INC})
list(APPEND EXTERNAL_INCLUDES ${VAR})
endforeach(VAR)
# append libraries
set(PETSC_LIBS_REQUIRED petscksp petscdm petscmat petscvec petsc)
find_all_libraries(PETSC_LIBS PETSC_LIBS_REQUIRED ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/lib "")
string(REGEX MATCH "PACKAGES_LIBS = [-_/\\.a-zA-Z0-9 ]*" PETSC_PACKAGES_LIBS ${VARS})
string(REPLACE "PACKAGES_LIBS = " "" PETSC_PACKAGES_LIBS ${PETSC_PACKAGES_LIBS})
list(APPEND EXTERNAL_LIBRARIES ${PETSC_LIBS} "${PETSC_PACKAGES_LIBS}")
endif(EXISTS ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/conf/petscvariables)
endif(ENABLE_PETSC)
if(ENABLE_OCC)
if(WIN32)
set(OCC_SYS_NAME win32)
......@@ -649,7 +677,7 @@ set_target_properties(lib PROPERTIES OUTPUT_NAME Gmsh)
if(MSVC)
SET_TARGET_PROPERTIES(lib PROPERTIES DEBUG_POSTFIX d)
if(ENABLE_OOFELIE)
INCLUDE(${CMAKE_SOURCE_DIR}/utils/misc/withOOFELIE.cmake)
include(${CMAKE_SOURCE_DIR}/utils/misc/withOOFELIE.cmake)
endif(ENABLE_OOFELIE)
endif(MSVC)
......
Common/GmshConfig.h.in
+ 1
− 0
View file @ 203ff04d
......@@ -35,6 +35,7 @@
#cmakedefine HAVE_OCC_MESH_CONSTRAINTS
#cmakedefine HAVE_OPENGL
#cmakedefine HAVE_OSMESA
#cmakedefine HAVE_PETSC
#cmakedefine HAVE_QT
#cmakedefine HAVE_SOLVER
#cmakedefine HAVE_TAUCS
......
Post/PViewDataRemote.h
+ 1
− 1
View file @ 203ff04d
......@@ -39,7 +39,7 @@ class PViewDataRemote : public PViewData {
int getNumElements(int step=-1, int ent=-1){ return -1; }
void setMin(double min){ _min = min; }
void setMax(double max){ _max = max; }
void setBoundingBox(SBoundingBox3d bbox){ _bbox = bbox; }
void setBoundingBox(SBoundingBox3d &bbox){ _bbox = bbox; }
void setTime(double time){ _time = time; }
bool isRemote(){ return true; }
int fillRemoteVertexArrays(std::string &options)
......
Post/PViewVertexArrays.cpp
+ 3
− 2
View file @ 203ff04d
......@@ -1179,10 +1179,11 @@ void PView::fillVertexArray(ConnectionManager *remote, int length, const char *b
data->setMin(min);
data->setMax(max);
data->setTime(time);
data->setBoundingBox(bbox);
}
}
p->getOptions()->tmpBBox += bbox;
// not perfect (does not take transformations into account)
p->getOptions()->tmpBBox = bbox;
switch(type){
case 1:
......
Solver/linearSystemCSR.cpp
+ 31
− 34
View file @ 203ff04d
......@@ -45,7 +45,6 @@ static void CSRList_Realloc(CSRList_T *liste,int n)
}
}
static CSRList_T *CSRList_Create(int n, int incr, int size)
{
CSRList_T *liste;
......@@ -77,7 +76,6 @@ static void CSRList_Delete(CSRList_T *liste)
void CSRList_Add(CSRList_T *liste, void *data)
{
liste->n++;
CSRList_Realloc(liste,liste->n);
liste->isorder = 0;
memcpy(&liste->array[(liste->n - 1) * liste->size], data, liste->size);
......@@ -89,40 +87,40 @@ int CSRList_Nbr(CSRList_T *liste)
}
template<>
void linearSystemCSR<double>::allocate(int _nbRows)
void linearSystemCSR<double>::allocate(int nbRows)
{
if(a_) {
CSRList_Delete(a_);
CSRList_Delete(ai_);
CSRList_Delete(ptr_);
CSRList_Delete(jptr_);
if(_a) {
CSRList_Delete(_a);
CSRList_Delete(_ai);
CSRList_Delete(_ptr);
CSRList_Delete(_jptr);
delete _x;
delete _b;
delete something;
}
if(_nbRows == 0){
a_ = 0;
ai_ = 0;
ptr_ = 0;
jptr_ = 0;
if(nbRows == 0){
_a = 0;
_ai = 0;
_ptr = 0;
_jptr = 0;
_b = 0;
_x = 0;
something = 0;
return;
}
a_ = CSRList_Create (_nbRows, _nbRows, sizeof(double));
ai_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
ptr_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
_a = CSRList_Create(nbRows, nbRows, sizeof(double));
_ai = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
_ptr = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
_jptr = CSRList_Create(nbRows + 1, nbRows, sizeof(INDEX_TYPE));
something = new char[_nbRows];
something = new char[nbRows];
for (int i=0;i<_nbRows;i++)something[i] = 0;
for (int i = 0; i < nbRows; i++) something[i] = 0;
_b = new std::vector<double>(_nbRows);
_x = new std::vector<double>(_nbRows);
_b = new std::vector<double>(nbRows);
_x = new std::vector<double>(nbRows);
}
const int NSTACK = 50;
......@@ -291,16 +289,16 @@ int linearSystemCSRGmm<double>::systemSolve()
{
if (!sorted)
sortColumns_(_b->size(),
CSRList_Nbr(a_),
(INDEX_TYPE *) ptr_->array,
(INDEX_TYPE *) jptr_->array,
(INDEX_TYPE *) ai_->array,
(double*) a_->array);
CSRList_Nbr(_a),
(INDEX_TYPE *) _ptr->array,
(INDEX_TYPE *) _jptr->array,
(INDEX_TYPE *) _ai->array,
(double*) _a->array);
sorted = true;
gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>
ref((double*)a_->array, (INDEX_TYPE *) ai_->array,
(INDEX_TYPE *)jptr_->array, _b->size(), _b->size());
ref((double*)_a->array, (INDEX_TYPE *) _ai->array,
(INDEX_TYPE *)_jptr->array, _b->size(), _b->size());
gmm::csr_matrix<double,0> M;
M.init_with(ref);
......@@ -313,4 +311,3 @@ int linearSystemCSRGmm<double>::systemSolve()
}
#endif
Solver/linearSystemCSR.h
+ 37
− 43
View file @ 203ff04d
......@@ -29,12 +29,12 @@ class linearSystemCSR : public linearSystem<scalar> {
protected:
bool sorted;
char *something;
CSRList_T *a_,*ai_,*ptr_,*jptr_;
CSRList_T *_a, *_ai, *_ptr, *_jptr;
std::vector<scalar> *_b, *_x;
public:
linearSystemCSR()
: sorted(false), a_(0) {}
virtual bool isAllocated() const { return a_ != 0; }
: sorted(false), _a(0) {}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int) ;
virtual void clear()
{
......@@ -46,67 +46,64 @@ class linearSystemCSR : public linearSystem<scalar> {
}
virtual void addToMatrix(int il, int ic, double val)
{
// if (sorted)throw;
INDEX_TYPE *jptr = (INDEX_TYPE*) _jptr->array;
INDEX_TYPE *ptr = (INDEX_TYPE*) _ptr->array;
INDEX_TYPE *ai = (INDEX_TYPE*) _ai->array;
scalar *a = ( scalar * ) _a->array;
INDEX_TYPE *jptr = (INDEX_TYPE*) jptr_->array;
INDEX_TYPE *ptr = (INDEX_TYPE*) ptr_->array;
INDEX_TYPE *ai = (INDEX_TYPE*) ai_->array;
scalar *a = ( scalar * ) a_->array;
INDEX_TYPE position_ = jptr[il];
INDEX_TYPE position = jptr[il];
if(something[il]) {
while(1){
if(ai[position_] == ic){
a[position_] += val;
// if (il == 0) printf("FOUND %d %d %d\n",il,ic,position_);
if(ai[position] == ic){
a[position] += val;
return;
}
if (ptr[position_] == 0)break;
position_ = ptr[position_];
if (ptr[position] == 0) break;
position = ptr[position];
}
}
INDEX_TYPE zero = 0;
CSRList_Add (a_, &val);
CSRList_Add (ai_, &ic);
CSRList_Add (ptr_, &zero);
CSRList_Add(_a, &val);
CSRList_Add(_ai, &ic);
CSRList_Add(_ptr, &zero);
// The pointers may have been modified
// if there has been a reallocation in CSRList_Add
ptr = (INDEX_TYPE*) ptr_->array;
ai = (INDEX_TYPE*) ai_->array;
a = (scalar*) a_->array;
ptr = (INDEX_TYPE*) _ptr->array;
ai = (INDEX_TYPE*) _ai->array;
a = (scalar*) _a->array;
INDEX_TYPE n = CSRList_Nbr(a_) - 1;
INDEX_TYPE n = CSRList_Nbr(_a) - 1;
if(!something[il]) {
jptr[il] = n;
something[il] = 1;
}
else ptr[position_] = n;
else ptr[position] = n;
}
virtual scalar getFromMatrix (int _row, int _col) const
virtual scalar getFromMatrix (int row, int col) const
{
throw;
Msg::Error("getFromMatrix not implemented for CSR");
return 0;
}
virtual void addToRightHandSide(int _row, scalar _val)
virtual void addToRightHandSide(int row, scalar val)
{
if(_val != 0.0) (*_b)[_row] += _val;
if(val != 0.0) (*_b)[row] += val;
}
virtual scalar getFromRightHandSide(int _row) const
virtual scalar getFromRightHandSide(int row) const
{
return (*_b)[_row];
return (*_b)[row];
}
virtual scalar getFromSolution(int _row) const
virtual scalar getFromSolution(int row) const
{
return (*_x)[_row];
return (*_x)[row];
}
virtual void zeroMatrix()
{
int N=CSRList_Nbr(a_);
scalar *a = (scalar*) a_->array;
int N = CSRList_Nbr(_a);
scalar *a = (scalar*) _a->array;
for (int i = 0; i < N; i++) a[i] = 0;
}
virtual void zeroRightHandSide()
......@@ -121,22 +118,19 @@ class linearSystemCSRGmm : public linearSystemCSR<scalar> {
double _prec;
int _noisy, _gmres;
public:
linearSystemCSRGmm()
: _prec(1.e-8), _noisy(0), _gmres(0) {}
virtual ~linearSystemCSRGmm()
{}
linearSystemCSRGmm() : _prec(1.e-8), _noisy(0), _gmres(0) {}
virtual ~linearSystemCSRGmm(){}
void setPrec(double p){ _prec = p; }
void setNoisy(int n){ _noisy = n; }
void setGmres(int n){ _gmres = n; }
virtual int systemSolve()
#if defined(HAVE_GMM)
;
#else
#if !defined(HAVE_GMM)
{
Msg::Error("Gmm++ is not available in this version of Gmsh");
return 0;
}
#endif
;
};
#endif
Solver/linearSystemFull.h
+ 14
− 15
View file @ 203ff04d
......@@ -22,12 +22,12 @@ class linearSystemFull : public linearSystem<scalar> {
public :
linearSystemFull() : _a(0), _b(0), _x(0){}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
virtual void allocate(int nbRows)
{
clear();
_a = new fullMatrix<scalar>(_nbRows, _nbRows);
_b = new fullVector<scalar>(_nbRows);
_x = new fullVector<scalar>(_nbRows);
_a = new fullMatrix<scalar>(nbRows, nbRows);
_b = new fullVector<scalar>(nbRows);
_x = new fullVector<scalar>(nbRows);
}
virtual ~linearSystemFull()
{
......@@ -42,25 +42,25 @@ class linearSystemFull : public linearSystem<scalar> {
}
_a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
virtual void addToMatrix(int row, int col, scalar val)
{
if(_val != 0.0) (*_a)(_row, _col) += _val;
if(val != 0.0) (*_a)(row, col) += val;
}
virtual scalar getFromMatrix(int _row, int _col) const
virtual scalar getFromMatrix(int row, int col) const
{
return (*_a)(_row, _col);
return (*_a)(row, col);
}
virtual void addToRightHandSide(int _row, scalar _val)
virtual void addToRightHandSide(int row, scalar val)
{
if(_val != 0.0) (*_b)(_row) += _val;
if(val != 0.0) (*_b)(row) += val;
}
virtual scalar getFromRightHandSide(int _row) const
virtual scalar getFromRightHandSide(int row) const
{
return (*_b)(_row);
return (*_b)(row);
}
virtual scalar getFromSolution(int _row) const
virtual scalar getFromSolution(int row) const
{
return (*_x)(_row);
return (*_x)(row);
}
virtual void zeroMatrix()
{
......@@ -74,7 +74,6 @@ class linearSystemFull : public linearSystem<scalar> {
{
if (_b->size())
_a->lu_solve(*_b, *_x);
// _x->print("********* mySol");
return 1;
}
};
......
Solver/linearSystemGMM.h
+ 19
− 19
View file @ 203ff04d
......@@ -26,12 +26,12 @@ class linearSystemGmm : public linearSystem<scalar> {
linearSystemGmm()
: _a(0), _b(0), _x(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
virtual void allocate(int nbRows)
{
clear();
_a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
_b = new std::vector<scalar>(_nbRows);
_x = new std::vector<scalar>(_nbRows);
_a = new gmm::row_matrix< gmm::wsvector<scalar> >(nbRows, nbRows);
_b = new std::vector<scalar>(nbRows);
_x = new std::vector<scalar>(nbRows);
}
virtual ~linearSystemGmm()
{
......@@ -46,25 +46,25 @@ class linearSystemGmm : public linearSystem<scalar> {
}
_a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
virtual void addToMatrix(int row, int col, scalar val)
{
if(_val != 0.0) (*_a)(_row, _col) += _val;
if(val != 0.0) (*_a)(row, col) += val;
}
virtual scalar getFromMatrix (int _row, int _col) const
virtual scalar getFromMatrix (int row, int col) const
{
return (*_a)(_row, _col);
return (*_a)(row, col);
}
virtual void addToRightHandSide(int _row, scalar _val)
virtual void addToRightHandSide(int row, scalar val)
{
if(_val != 0.0) (*_b)[_row] += _val;
if(val != 0.0) (*_b)[row] += val;
}
virtual scalar getFromRightHandSide(int _row) const
virtual scalar getFromRightHandSide(int row) const
{
return (*_b)[_row];
return (*_b)[row];
}
virtual scalar getFromSolution(int _row) const
virtual scalar getFromSolution(int row) const
{
return (*_x)[_row];
return (*_x)[row];
}
virtual void zeroMatrix()
{
......@@ -100,11 +100,11 @@ class linearSystemGmm : public linearSystem<scalar> {
}
virtual bool isAllocated() const { return false; }
virtual void allocate(int nbRows) {}
virtual void addToMatrix(int _row, int _col, scalar val) {}
virtual scalar getFromMatrix(int _row, int _col) const { return 0.; }
virtual void addToRightHandSide(int _row, scalar val) {}
virtual scalar getFromRightHandSide(int _row) const { return 0.; }
virtual scalar getFromSolution(int _row) const { return 0.; }
virtual void addToMatrix(int row, int col, scalar val) {}
virtual scalar getFromMatrix(int row, int col) const { return 0.; }
virtual void addToRightHandSide(int row, scalar val) {}
virtual scalar getFromRightHandSide(int row) const { return 0.; }
virtual scalar getFromSolution(int row) const { return 0.; }
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual int systemSolve() { return 0; }
......
Solver/linearSystemTAUCS.cpp
+ 11
− 11
View file @ 203ff04d
......@@ -26,21 +26,21 @@ int linearSystemCSRTaucs<double>::systemSolve()
{
if(!sorted){
sortColumns_(_b->size(),
CSRList_Nbr(a_),
(INDEX_TYPE *) ptr_->array,
(INDEX_TYPE *) jptr_->array,
(INDEX_TYPE *) ai_->array,
(double*) a_->array);
CSRList_Nbr(_a),
(INDEX_TYPE *) _ptr->array,
(INDEX_TYPE *) _jptr->array,
(INDEX_TYPE *) _ai->array,
(double*) _a->array);
}
sorted = true;
taucs_ccs_matrix myVeryCuteTaucsMatrix;
myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
//myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ptr_->array;
//myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)ai_->array;
myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)ai_->array;
myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)jptr_->array;
myVeryCuteTaucsMatrix.values.d = (double*) a_->array;
//myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ptr->array;
//myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_ai->array;
myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ai->array;
myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_jptr->array;
myVeryCuteTaucsMatrix.values.d = (double*)_a->array;
myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
char* options[] = { "taucs.factor.LLT=true", NULL };
......
Solver/linearSystemTAUCS.h
+ 7
− 5
View file @ 203ff04d
......@@ -19,11 +19,13 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
linearSystemCSR<scalar>::addToMatrix(il, ic, val);
}
virtual int systemSolve()
#if defined(HAVE_TAUCS)
;
#else
{ throw; }
#if !defined(HAVE_TAUCS)
{
Msg::Error("TAUCS is not available in this version of Gmsh");
return 0;
}
#endif
;
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment