diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt index 0a4100756c2cbb7a1d0c3c78b4c3f1e18e750c24..b936f23343264acdb8672d0fb083da46fea26e89 100644 --- a/Solver/CMakeLists.txt +++ b/Solver/CMakeLists.txt @@ -18,6 +18,7 @@ set(SRC luaFunction.cpp functionSpace.cpp filters.cpp + sparsityPattern.cpp ) file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) diff --git a/Solver/sparsityPattern.cpp b/Solver/sparsityPattern.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f7cc6662a7bf11b3939370e71f66fe2104a72712 --- /dev/null +++ b/Solver/sparsityPattern.cpp @@ -0,0 +1,86 @@ +// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. +// +// Contributor(s): +// Jonathan Lambrechts +// + +#include "sparsityPattern.h" + +#include<stdlib.h> +#include<string.h> + +// this class has been optimized, please before changing anything, check twice : +// the impact on the performance to assemble typical High Order FE problems +// and the impact on the memory usage for this operation + +sparsityPattern::~sparsityPattern() +{ + for(int i = 0; i <_nRows; i++) { + if (_rowsj[i]) free(_rowsj[i]); + } + free(_nByRow); + free(_nAllocByRows); + free(_rowsj); +} + +void sparsityPattern::addEntry (int i, int j) +{ + int n = _nByRow[i]; + int *rowj = _rowsj[i]; + int k = 0; + int k0 = 0, k1 = n; + if (n>20) { + while (k1-k0 > 20) { + int k2 = ((k0+k1)/2); + if (rowj[k2] > j) + k1 = k2; + else if (rowj[k2] < j) + k0 = k2+1; + else + return; + } + for (k = k0; k < k1; k++) { + if (rowj[k] >= j) { + if (rowj[k] == j) return; + break; + } + } + } else { // this "if() else" is not strictly necessary but with it, gcc unroll the for(k) loop + for (k = 0; k < n; k++) { + if (rowj[k] >= j) { + if (rowj[k] == j) return; + break; + } + } + } + _nByRow[i] = n+1; + if (_nByRow[i]> _nAllocByRows[i]) { + int na = (n+1)*3/2; + _rowsj[i] = (int*)realloc(_rowsj[i], (na*sizeof(int))); + _nAllocByRows[i] = na; + } + memmove(&_rowsj[i][k+1], &_rowsj[i][k], (n-k)*sizeof (int)); + _rowsj[i][k] = j; +} + +sparsityPattern::sparsityPattern (int ni) +{ + _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; + } +} + +const int *sparsityPattern::getRow (int i, int &size) const +{ + size = _nByRow[i]; + return _rowsj[i]; +} diff --git a/Solver/sparsityPattern.h b/Solver/sparsityPattern.h new file mode 100644 index 0000000000000000000000000000000000000000..ab1284c96908aae690e26e2d5abd3c015cadaf87 --- /dev/null +++ b/Solver/sparsityPattern.h @@ -0,0 +1,25 @@ +// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. + +#ifndef _SPARSITY_PATTERN_H_ +#define _SPARSITY_PATTERN_H_ + +// this class has been optimized, please before changing anything, check twice : +// - the impact on the performance to assemble typical High Order FE problems +// - the impact on the memory for this operation + +class sparsityPattern { + int *_nByRow, *_nAllocByRows; + int **_rowsj; + int _nRows; + + public : + void addEntry (int i, int j); + const int* getRow (int line, int &size) const; + sparsityPattern (int nRows); + ~sparsityPattern(); +}; + +#endif