Solver : support for exact preallocation of linearSystems (only

used in linearSystemPETSc for now)

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();

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>

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));

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));

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];
 }

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