From 5eae401eca188348e6c6d87be8843f77d21e95df Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 19 Oct 2010 17:55:00 +0000
Subject: [PATCH] Solver : linearSystemCSR : support for preallocation using
sparsityPattern
---
Solver/linearSystemCSR.cpp | 47 ++++++++++++++++++++++++++++++++++++++
Solver/linearSystemCSR.h | 16 ++++++++++++-
2 files changed, 62 insertions(+), 1 deletion(-)
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 898910ba6a..20ab399b1e 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/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> >();
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 6387315df9..661afc0b6c 100644
--- a/Solver/linearSystemCSR.h
+++ b/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)
--
GitLab