From 9e76055a0721f9ce3d96b13d3460fe11cd7e00f1 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 19 Oct 2010 10:17:12 +0000
Subject: [PATCH] Solver : support for exact preallocation of linearSystems
(only used in linearSystemPETSc for now)
---
Solver/dofManager.h | 47 +++++++++++++++++++++++++++
Solver/linearSystem.h | 1 +
Solver/linearSystemPETSc.cpp | 13 +++++---
Solver/linearSystemPETSc.h | 62 ++++++++++++++++++++++++++++++++----
Solver/sparsityPattern.cpp | 57 +++++++++++++++++++++++----------
Solver/sparsityPattern.h | 10 +++---
6 files changed, 159 insertions(+), 31 deletions(-)
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index baed02f0e0..542613c2a5 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -274,6 +274,53 @@ class dofManager{
{
getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
}
+
+ inline void insertInSparsityPatternLinConst(const Dof &R, const Dof &C)
+ {
+ std::map<Dof, int>::iterator itR = unknown.find(R);
+ if (itR != unknown.end())
+ {
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(C);
+ if (itConstraint != constraints.end()){
+ for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
+ insertInSparsityPattern(R,(itConstraint->second).linear[i].first);
+ }
+ }
+ }
+ else{ // test function ; (no shift ?)
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(R);
+ if (itConstraint != constraints.end()){
+ for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
+ insertInSparsityPattern((itConstraint->second).linear[i].first, C);
+ }
+ }
+ }
+ }
+
+ inline void insertInSparsityPattern(const Dof &R, const Dof &C)
+ {
+ if (_isParallel && !_parallelFinalized) _parallelFinalize();
+ if (!_current->isAllocated()) _current->allocate(sizeOfR());
+ std::map<Dof, int>::iterator itR = unknown.find(R);
+ if (itR != unknown.end()){
+ std::map<Dof, int>::iterator itC = unknown.find(C);
+ if (itC != unknown.end()){
+ _current->insertInSparsityPattern(itR->second, itC->second);
+ }
+ else{
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
+ if (itFixed != fixed.end()) {
+ }else insertInSparsityPatternLinConst(R, C);
+ }
+ }
+ if (itR == unknown.end())
+ {
+ insertInSparsityPatternLinConst(R, C);
+ }
+ }
+
inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
{
if (_isParallel && !_parallelFinalized) _parallelFinalize();
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index 687aa14d5c..f120581da0 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -22,6 +22,7 @@ class linearSystemBase {
virtual void zeroRightHandSide() = 0;
virtual int systemSolve() = 0;
void setParameter (std::string key, std::string value);
+ virtual void insertInSparsityPattern(int _row, int _col){};
};
template <class scalar>
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index b2f6a800dc..ada1607eb8 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -9,6 +9,8 @@
template <>
void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
{
+ if (!_entriesPreAllocated)
+ preAllocateEntries();
fullMatrix<PetscScalar> &modval = *const_cast<fullMatrix<PetscScalar> *>(&val);
for (int ii = 0; ii < val.size1(); ii++)
for (int jj = 0; jj < ii; jj++) {
@@ -93,13 +95,16 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
_try(MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE));
if (Msg::GetCommSize() > 1) {
_try(MatSetType(_a, MATMPIBAIJ));
- _try(MatSetFromOptions(_a));
- _try(MatMPIBAIJSetPreallocation(_a, _blockSize, 5, NULL, 0, NULL));
} else {
_try(MatSetType(_a, MATSEQBAIJ));
- _try(MatSetFromOptions(_a));
- _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, NULL)); //todo preAllocate off-diagonal part
}
+ _try(MatSetFromOptions(_a));
+ _try(MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd));
+ _try(MatGetSize(_a, &_globalSize, &_localSize));
+ _globalSize /= _blockSize;
+ _localSize /= _blockSize;
+ _localRowStart /= _blockSize;
+ _localRowEnd /= _blockSize;
// override the default options with the ones from the option
// database (if any)
_try(VecCreate(PETSC_COMM_WORLD, &_x));
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 7bca997fb1..761a3c3a4f 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -37,6 +37,8 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "linearSystem.h"
+#include "sparsityPattern.h"
+#include <vector>
#if defined(HAVE_PETSC)
#include <petsc.h>
@@ -46,17 +48,19 @@ template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
protected:
int _blockSize; // for block Matrix
- bool _isAllocated, _kspAllocated;
+ bool _isAllocated, _kspAllocated, _entriesPreAllocated;
Mat _a;
Vec _b, _x;
KSP _ksp;
+ int _localRowStart, _localRowEnd, _localSize, _globalSize;
+ sparsityPattern _sparsity;
static void _try(int ierr) { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
void _kspCreate() {
_try(KSPCreate(PETSC_COMM_WORLD, &_ksp));
PC pc;
_try(KSPGetPC(_ksp, &pc));
// set some default options
- _try(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect solver
+ //_try(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect solver
_try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
_try(PCFactorSetLevels(pc, 1));
_try(KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
@@ -72,7 +76,47 @@ class linearSystemPETSc : public linearSystem<scalar> {
if (_kspAllocated)
_try (KSPDestroy (_ksp));
}
+ void insertInSparsityPattern (int i, int j) {
+ i -= _localRowStart;
+ if (i<0 || i>= _localSize) return;
+ _sparsity.insertEntry (i,j);
+ }
virtual bool isAllocated() const { return _isAllocated; }
+ virtual void preAllocateEntries() {
+ if (_entriesPreAllocated) return;
+ if (!_isAllocated) Msg::Fatal("system must be allocated first");
+ if (_sparsity.getNbRows() == 0) {
+ PetscInt prealloc = 300;
+ PetscTruth set;
+ PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
+ if (_blockSize == 0) {
+ _try(MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL));
+ } else {
+ _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL));
+ }
+ } else {
+ std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
+ for (int i = 0; i < _localSize; i++) {
+ int n;
+ const int *r = _sparsity.getRow(i, n);
+ for (int j = 0; j < n; j++) {
+ if (r[j] >= _localRowStart && r[j] < _localRowEnd)
+ nByRowDiag[i] ++;
+ else
+ nByRowOffDiag[i] ++;
+ }
+ }
+ if (_blockSize == 0) {
+ _try(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]));
+ _try(MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
+ } else {
+ _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]));
+ _try(MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]));
+ }
+ _sparsity.clear();
+ }
+ _entriesPreAllocated = true;
+ }
virtual void allocate(int nbRows)
{
clear();
@@ -81,11 +125,9 @@ class linearSystemPETSc : public linearSystem<scalar> {
// override the default options with the ones from the option
// database (if any)
_try(MatSetFromOptions(_a));
+ _try(MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd));
+ _try(MatGetSize(_a, &_globalSize, &_localSize));
// preallocation option must be set after other options
- PetscInt prealloc = 300;
- PetscTruth set;
- PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
- _try(MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL));
_try(VecCreate(PETSC_COMM_WORLD, &_x));
_try(VecSetSizes(_x, nbRows, PETSC_DETERMINE));
// override the default options with the ones from the option
@@ -93,6 +135,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
_try(VecSetFromOptions(_x));
_try(VecDuplicate(_x, &_b));
_isAllocated = true;
+ _entriesPreAllocated = false;
}
void print() {
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
@@ -149,6 +192,8 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
virtual void addToMatrix(int row, int col, const scalar &val)
{
+ if (!_entriesPreAllocated)
+ preAllocateEntries();
PetscInt i = row, j = col;
PetscScalar s = val;
_try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
@@ -167,7 +212,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
virtual void zeroMatrix()
{
- if (_isAllocated) {
+ if (_isAllocated && _entriesPreAllocated) {
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
_try(MatZeroEntries(_a));
@@ -193,6 +238,9 @@ class linearSystemPETSc : public linearSystem<scalar> {
_try(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
_try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
+ /*MatInfo info;
+ MatGetInfo(_a, MAT_LOCAL, &info);
+ printf("mallocs %.0f nz_allocated %.0f nz_used %.0f nz_unneeded %.0f\n", info.mallocs, info.nz_allocated, info.nz_used, info.nz_unneeded);*/
_try(VecAssemblyBegin(_b));
_try(VecAssemblyEnd(_b));
_try(KSPSolve(_ksp, _b, _x));
diff --git a/Solver/sparsityPattern.cpp b/Solver/sparsityPattern.cpp
index f7cc6662a7..a879201be3 100644
--- a/Solver/sparsityPattern.cpp
+++ b/Solver/sparsityPattern.cpp
@@ -17,17 +17,42 @@
// and the impact on the memory usage for this operation
sparsityPattern::~sparsityPattern()
+{
+ clear();
+}
+
+void sparsityPattern::clear()
{
for(int i = 0; i <_nRows; i++) {
if (_rowsj[i]) free(_rowsj[i]);
}
- free(_nByRow);
- free(_nAllocByRows);
- free(_rowsj);
+ if (_nByRow) free(_nByRow);
+ if (_nAllocByRow) free(_nAllocByRow);
+ if (_rowsj) free(_rowsj);
+ _nByRow = NULL;
+ _rowsj = NULL;
+ _nAllocByRow = NULL;
+ _nRows = 0;
+ _nRowsAlloc = 0;
}
-void sparsityPattern::addEntry (int i, int j)
+void sparsityPattern::insertEntry (int i, int j)
{
+ if (i >= _nRows) {
+ if (i >= _nRowsAlloc) {
+ _nRowsAlloc = (i+1)*3/2;
+ _rowsj = (int**)realloc (_rowsj, sizeof(int*)*_nRowsAlloc);
+ _nByRow = (int*)realloc (_nByRow, sizeof(int)*_nRowsAlloc);
+ _nAllocByRow = (int*)realloc(_nAllocByRow, sizeof(int)*_nRowsAlloc);
+ }
+ for (int k = _nRows; k <= i; k++) {
+ _nByRow[k] = 0;
+ _nAllocByRow[k] = 0;
+ _rowsj[k] = NULL;
+ }
+ _nRows = i+1;
+ }
+
int n = _nByRow[i];
int *rowj = _rowsj[i];
int k = 0;
@@ -57,30 +82,30 @@ void sparsityPattern::addEntry (int i, int j)
}
}
_nByRow[i] = n+1;
- if (_nByRow[i]> _nAllocByRows[i]) {
+ if (_nByRow[i]> _nAllocByRow[i]) {
int na = (n+1)*3/2;
_rowsj[i] = (int*)realloc(_rowsj[i], (na*sizeof(int)));
- _nAllocByRows[i] = na;
+ _nAllocByRow[i] = na;
}
memmove(&_rowsj[i][k+1], &_rowsj[i][k], (n-k)*sizeof (int));
_rowsj[i][k] = j;
}
-sparsityPattern::sparsityPattern (int ni)
+sparsityPattern::sparsityPattern ()
{
- _nRows = ni;
- _rowsj = (int**)malloc (sizeof(int*)*ni);
- _nByRow = (int*)malloc(sizeof(int)*ni);
- _nAllocByRows = (int*)malloc(sizeof(int)*ni);
- for (int i = 0; i < ni; i++) {
- _nByRow[i] = 0;
- _nAllocByRows[i] = 0;
- _rowsj[i] = NULL;
- }
+ _nRows = 0;
+ _nRowsAlloc = 0;
+ _rowsj = NULL;
+ _nByRow = NULL;
+ _nAllocByRow = NULL;
}
const int *sparsityPattern::getRow (int i, int &size) const
{
+ if (i >= _nRows){
+ size = 0;
+ return NULL;
+ }
size = _nByRow[i];
return _rowsj[i];
}
diff --git a/Solver/sparsityPattern.h b/Solver/sparsityPattern.h
index ab1284c969..4f826208db 100644
--- a/Solver/sparsityPattern.h
+++ b/Solver/sparsityPattern.h
@@ -11,15 +11,17 @@
// - the impact on the memory for this operation
class sparsityPattern {
- int *_nByRow, *_nAllocByRows;
+ int *_nByRow, *_nAllocByRow;
int **_rowsj;
- int _nRows;
+ int _nRows, _nRowsAlloc;
public :
- void addEntry (int i, int j);
+ void insertEntry (int i, int j);
const int* getRow (int line, int &size) const;
- sparsityPattern (int nRows);
+ void clear();
+ sparsityPattern ();
~sparsityPattern();
+ inline int getNbRows() {return _nRows;}
};
#endif
--
GitLab