Solver : linearSystemCSR : support for preallocation using

sparsityPattern

Solver/linearSystemCSR.cpp

@@ -46,6 +46,12 @@ static void CSRList_Realloc(CSRList_T *liste, int n)
   }
 }
 
+static void CSRList_Resize(CSRList_T *liste, int n)
+{
+  CSRList_Realloc(liste, n);
+  liste->n = n;
+}
+
 static CSRList_T *CSRList_Create(int n, int incr, int size)
 {
   CSRList_T *liste;
@@ -87,6 +93,46 @@ int CSRList_Nbr(CSRList_T *liste)
   return(liste->n);
 }
 
+template<>
+void linearSystemCSR<double>::preAllocateEntries ()
+{
+  if (_entriesPreAllocated) return;
+  if (_sparsity.getNbRows() == 0) return;
+  INDEX_TYPE nnz = 0;
+  int nbRows = _b->size();
+  for (int i = 0; i < nbRows; i++){
+    int nInRow;
+    _sparsity.getRow (i, nInRow);
+    nnz += nInRow;
+  }
+  CSRList_Resize (_a, nnz);
+  CSRList_Resize (_ai, nnz);
+  CSRList_Resize (_ptr, nnz);
+  INDEX_TYPE *jptr = (INDEX_TYPE*) _jptr->array;
+  INDEX_TYPE *ai = (INDEX_TYPE*) _ai->array;
+  INDEX_TYPE *ptr = (INDEX_TYPE*) _ptr->array;
+  double *a = ( double * ) _a->array;
+  jptr[0] = 0;
+  nnz = 0;
+  for (int i = 0; i < nbRows; i++){
+    int nInRow;
+    const int *row = _sparsity.getRow (i, nInRow);
+    for (int j = 0; j < nInRow; j++) {
+      ai[nnz] = row[j];
+      a[nnz] = 0.;
+      ptr[nnz] = nnz+1;
+      nnz ++;
+    }
+    if (nInRow != 0)
+      ptr[nnz-1] = 0;
+    jptr[i+1] = nnz;
+    something[i] = (nInRow == 0 ? 0 : 1);
+  }
+  _entriesPreAllocated = true;
+  sorted = true;
+  _sparsity.clear();
+}
+
 template<>
 void linearSystemCSR<double>::allocate(int nbRows)
 {
@@ -297,6 +343,7 @@ template<>
   void linearSystemCSRGmm<double>::registerBindings(binding *b)
   {
     classBinding *cb = b->addClass< linearSystemCSRGmm<double> >("linearSystemCSRGmmdouble");
+    cb->setParentClass<linearSystem<double> >();
     cb->setDescription("Sparse matrix representation.");
     methodBinding *cm;
     cm = cb->setConstructor<linearSystemCSRGmm<double> >();

Solver/linearSystemCSR.h

@@ -10,6 +10,7 @@
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "linearSystem.h"
+#include "sparsityPattern.h"
 
 class binding;
 
@@ -29,10 +30,12 @@ int  CSRList_Nbr(CSRList_T *liste);
 template <class scalar>
 class linearSystemCSR : public linearSystem<scalar> {
  protected:
+  bool _entriesPreAllocated;
   bool sorted;
   char *something;
   CSRList_T *_a, *_ai, *_ptr, *_jptr;
   std::vector<scalar> *_b, *_x;
+  sparsityPattern _sparsity; // only used for pre-allocation, does not store the sparsity once allocated
  public:
   linearSystemCSR()
     : sorted(false), _a(0), _b(0), _x(0) {}
@@ -46,8 +49,14 @@ class linearSystemCSR : public linearSystem<scalar> {
   {
     allocate(0);
   }
+  virtual void insertInSparsityPattern (int i, int j) {
+    _sparsity.insertEntry (i,j);
+  }
+  virtual void preAllocateEntries ();
   virtual void addToMatrix(int il, int ic, const scalar &val) 
   {
+    if (!_entriesPreAllocated)
+      preAllocateEntries();
     INDEX_TYPE  *jptr  = (INDEX_TYPE*) _jptr->array;
     INDEX_TYPE  *ptr   = (INDEX_TYPE*) _ptr->array;
     INDEX_TYPE  *ai    = (INDEX_TYPE*) _ai->array;
@@ -180,8 +189,13 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
   virtual ~linearSystemCSRTaucs(){}
   virtual void addToMatrix(int il, int ic, const double &val)
   {
-    if (il <= ic)
+    if (il <= ic) {
       linearSystemCSR<scalar>::addToMatrix(il, ic, val);
+    }
+  }
+  virtual void insertInSparsityPattern(int il, int ic) {
+    if (il <= ic)
+      linearSystemCSR<scalar>::insertInSparsityPattern(il,ic);
   }
   virtual int systemSolve()
 #if !defined(HAVE_TAUCS)