diff --git a/Numeric/gmshLinearSystemCSR.cpp b/Numeric/gmshLinearSystemCSR.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5defc796c031548d7d23a04fceb76db1ebc6bc0
--- /dev/null
+++ b/Numeric/gmshLinearSystemCSR.cpp
@@ -0,0 +1,345 @@
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "gmshLinearSystemCSR.h"
+
+#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
+#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
+
+// #if defined(HAVE_TAUCS)
+// #include "taucs.h"
+// #endif
+
+// #if defined(HAVE_GMM)
+// #include <gmm.h>
+// #include <gmm_interface.h>
+// #endif
+
+
+void *CSRMalloc(size_t size)
+{
+ void *ptr;
+ if (!size) return(NULL);
+ ptr = malloc(size);
+ return(ptr);
+}
+
+void *CSRRealloc(void *ptr, size_t size)
+{
+ if (!size) return(NULL);
+ ptr = realloc(ptr,size);
+ return(ptr);
+}
+
+void CSRList_Realloc(CSRList_T *liste,int n)
+{
+ char* temp;
+ if (n <= 0) return;
+ if (liste->array == NULL) {
+ liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+ liste->array = (char *)CSRMalloc(liste->nmax * liste->size);
+ }
+ else {
+ if (n > liste->nmax) {
+ liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+ temp = (char *)CSRRealloc(liste->array, liste->nmax * liste->size);
+ liste->array = temp;
+ }
+ }
+}
+
+
+CSRList_T *CSRList_Create(int n, int incr, int size)
+{
+ CSRList_T *liste;
+
+ if (n <= 0) n = 1 ;
+ if (incr <= 0) incr = 1;
+
+ liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
+
+ liste->nmax = 0;
+ liste->incr = incr;
+ liste->size = size;
+ liste->n = 0;
+ liste->isorder = 0;
+ liste->array = NULL;
+
+ CSRList_Realloc(liste,n);
+ return(liste);
+}
+
+void CSRList_Delete(CSRList_T *liste)
+{
+ if (liste != 0) {
+ free(liste->array);
+ free(liste);
+ }
+}
+
+
+
+void CSRList_Add(CSRList_T *liste, void *data)
+{
+ liste->n++;
+
+ CSRList_Realloc(liste,liste->n);
+ liste->isorder = 0;
+ memcpy(&liste->array[(liste->n - 1) * liste->size],data,liste->size);
+}
+
+int CSRList_Nbr(CSRList_T *liste)
+{
+ return(liste->n);
+}
+
+template<>
+void gmshLinearSystemCSR<double> :: allocate(int _nbRows)
+{
+ if(a_) {
+ CSRList_Delete(a_);
+ CSRList_Delete(ai_);
+ CSRList_Delete(ptr_);
+ CSRList_Delete(jptr_);
+ delete _x;
+ delete _b;
+ delete something;
+ }
+
+ if (_nbRows == 0){
+ a_ = 0;
+ ai_ = 0;
+ ptr_ = 0;
+ jptr_ = 0;
+ _b = 0;
+ _x = 0;
+ something = 0;
+ return;
+ }
+
+ a_ = CSRList_Create (_nbRows, _nbRows, sizeof(double));
+ ai_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+ ptr_ = CSRList_Create (_nbRows, _nbRows, sizeof(INDEX_TYPE));
+ jptr_ = CSRList_Create (_nbRows+1, _nbRows, sizeof(INDEX_TYPE));
+
+ something = new char[_nbRows];
+
+ for (int i=0;i<_nbRows;i++)something[i] = 0;
+
+ _b = new std::vector<double>(_nbRows);
+ _x = new std::vector<double>(_nbRows);
+}
+
+const int NSTACK = 50;
+const int M_sort2 = 7;
+
+static void free_ivector(int *v, long nl, long nh){
+ // free an int vector allocated with ivector()
+ free((char*)(v+nl-1));
+}
+
+static int * ivector(long nl, long nh) {
+ // allocate an int vector with subscript range v[nl..nh]
+ int *v;
+ v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
+ if (!v) fprintf(stderr, "allocation failure in ivector()\n");
+ return v-nl+1;
+}
+
+static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj){
+ if(ai<bi)return -1;
+ if(ai>bi)return 1;
+ if(aj<bj)return -1;
+ if(aj>bj)return 1;
+ return 0;
+}
+template <class scalar>
+void _sort2_xkws (unsigned long n, double arr[], INDEX_TYPE ai[] , INDEX_TYPE aj [] ) {
+
+ unsigned long i,ir=n,j,k,l=1;
+ int *istack,jstack=0;
+ INDEX_TYPE tempi;
+ scalar a,temp;
+ int b,c;
+
+ istack=ivector(1,NSTACK);
+ for (;;) {
+ if (ir-l < M_sort2) {
+ for (j=l+1;j<=ir;j++) {
+ a=arr[j -1];
+ b=ai[j -1];
+ c=aj[j -1];
+ for (i=j-1;i>=1;i--) {
+ if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
+ arr[i+1 -1]=arr[i -1];
+ ai[i+1 -1]=ai[i -1];
+ aj[i+1 -1]=aj[i -1];
+ }
+ arr[i+1 -1]=a;
+ ai[i+1 -1]=b;
+ aj[i+1 -1]=c;
+ }
+ if (!jstack) {
+ free_ivector(istack,1,NSTACK);
+ return;
+ }
+ ir=istack[jstack];
+ l=istack[jstack-1];
+ jstack -= 2;
+ }
+ else {
+ k=(l+ir) >> 1;
+ SWAP(arr[k -1],arr[l+1 -1])
+ SWAPI(ai[k -1],ai[l+1 -1])
+ SWAPI(aj[k -1],aj[l+1 -1])
+ if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
+ SWAP(arr[l+1 -1],arr[ir -1])
+ SWAPI(ai[l+1 -1],ai[ir -1])
+ SWAPI(aj[l+1 -1],aj[ir -1])
+ }
+ if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
+ SWAP(arr[l -1],arr[ir -1])
+ SWAPI(ai[l -1],ai[ir -1])
+ SWAPI(aj[l -1],aj[ir -1])
+ }
+ if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
+ SWAP(arr[l+1 -1],arr[l -1])
+ SWAPI(ai[l+1 -1],ai[l -1])
+ SWAPI(aj[l+1 -1],aj[l -1])
+ }
+ i=l+1;
+ j=ir;
+ a=arr[l -1];
+ b=ai[l -1];
+ c=aj[l -1];
+ for (;;) {
+ do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
+ do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
+ if (j < i) break;
+ SWAP(arr[i -1],arr[j -1])
+ SWAPI(ai[i -1],ai[j -1])
+ SWAPI(aj[i -1],aj[j -1])
+ }
+ arr[l -1]=arr[j -1];
+ arr[j -1]=a;
+ ai[l -1]=ai[j -1];
+ ai[j -1]=b;
+ aj[l -1]=aj[j -1];
+ aj[j -1]=c;
+ jstack += 2;
+ if (jstack > NSTACK) {
+ Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
+ throw;
+ }
+ if (ir-i+1 >= j-l) {
+ istack[jstack]=ir;
+ istack[jstack-1]=i;
+ ir=j-1;
+ }
+ else {
+ istack[jstack]=j-1;
+ istack[jstack-1]=l;
+ l=i;
+ }
+ }
+ }
+}
+
+template <class scalar>
+void sortColumns (int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a) {
+ // replace pointers by lines
+ int *count = new int [NbLines];
+
+ for(int i=0;i<NbLines;i++){
+ count[i] = 0;
+ INDEX_TYPE _position = jptr[i];
+ while(1){
+ count[i]++;
+ INDEX_TYPE _position_temp = _position;
+ _position = ptr[_position];
+ ptr[_position_temp] = i;
+ if (_position == 0) break;
+ }
+ }
+ _sort2_xkws<double>(nnz,a,ptr,ai);
+ jptr[0] = 0;
+ for(int i=1;i<=NbLines;i++){
+ jptr[i] = jptr[i-1] + count[i-1];
+ }
+
+ for(int i=0;i<NbLines;i++){
+ for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
+ ptr[j] = j+1;
+ }
+ ptr[jptr[i+1]] = 0;
+ }
+
+
+ delete[] count;
+}
+
+#if defined(HAVE_GMM)
+#include "gmm.h"
+template<>
+int gmshLinearSystemCSRGmm<double> :: systemSolve()
+{
+ sortColumns(_b->size(),
+ CSRList_Nbr(a_),
+ (INDEX_TYPE *) ptr_->array,
+ (INDEX_TYPE *) jptr_->array,
+ (INDEX_TYPE *) ai_->array,
+ (double*) a_->array);
+ sorted = true;
+
+ // for (int i=0;i<_b->size();i++)
+ // printf("%d ",((INDEX_TYPE *) jptr_->array)[i]);
+ // printf("\n");
+
+ gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0> ref((double*) a_->array,
+ (INDEX_TYPE *) ai_->array,
+ (INDEX_TYPE *) jptr_->array,
+ _b->size(), _b->size());
+ gmm::csr_matrix<double,0> M;
+ M.init_with(ref);
+
+ gmm::ildltt_precond<gmm::csr_matrix<double,0> > P(M, 10, 1.e-10);
+ gmm::iteration iter(_prec);
+ iter.set_noisy(_noisy);
+ if(_gmres) gmm::gmres(M, *_x, *_b, P, 100, iter);
+ else gmm::cg(M, *_x, *_b, P, iter);
+ return 1;
+}
+#endif
+
+#if defined(HAVE_TAUCSw)
+#include "taucs.h"
+template<>
+int gmshLinearSystemCSRTaucs<double> :: systemSolve()
+{
+ sortColumns(_b->size(),
+ CSRList_Nbr(a_),
+ (INDEX_TYPE *) ptr_->array,
+ (INDEX_TYPE *) jptr_->array,
+ (INDEX_TYPE *) ai_->array,
+ (double*) a_->array);
+ sorted = true;
+
+ taucs_ccs_matrix myVeryCuteTaucsMatrix;
+ myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
+ myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)jptr_->array;
+ myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)ai_->array;
+ myVeryCuteTaucsMatrix.values.d = (double*) a_->array;
+ char* options[] = { "taucs.factor.LLT=true", NULL };
+ int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
+ NULL, 1, &(*_x)[0],&(*_b)[0],
+ options,NULL);
+ if (result != TAUCS_SUCCESS){
+ Msg::Error("Taucs Was Not Successfull");
+ }
+ return 1;
+}
+#endif
diff --git a/Numeric/gmshLinearSystemCSR.h b/Numeric/gmshLinearSystemCSR.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b8a1fb5438c5713d3c091c6b6e876f1f8f33105
--- /dev/null
+++ b/Numeric/gmshLinearSystemCSR.h
@@ -0,0 +1,169 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _GMSH_LINEAR_SYSTEM_TAUCS_H_
+#define _GMSH_LINEAR_SYSTEM_TAUCS_H_
+
+// Interface to TAUCS
+
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "gmshLinearSystem.h"
+
+#include <vector>
+
+typedef int INDEX_TYPE ;
+
+typedef struct {
+ int nmax;
+ int size;
+ int incr;
+ int n;
+ int isorder;
+ char *array;
+} CSRList_T;
+
+void CSRList_Add(CSRList_T *liste, void *data);
+int CSRList_Nbr(CSRList_T *liste);
+
+template <class scalar>
+class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
+ protected:
+ bool sorted;
+ char *something;
+ CSRList_T *a_,*ai_,*ptr_,*jptr_;
+ std::vector<scalar> *_b, *_x;
+ public:
+ gmshLinearSystemCSR()
+ : a_(0),sorted(false) {}
+ virtual bool isAllocated() const { return a_ != 0; }
+ virtual void allocate(int) ;
+ virtual ~gmshLinearSystemCSR()
+ {
+ allocate(0);
+ }
+ void addToMatrix ( int il, int ic, double val)
+ {
+ // if (sorted)throw;
+
+ INDEX_TYPE *jptr = (INDEX_TYPE*) jptr_->array;
+ INDEX_TYPE *ptr = (INDEX_TYPE*) ptr_->array;
+ INDEX_TYPE *ai = (INDEX_TYPE*) ai_->array;
+ scalar *a = ( scalar * ) a_->array;
+
+ INDEX_TYPE position_ = jptr[il];
+
+ if(something[il]) {
+ while(1){
+ if(ai[position_] == ic){
+ a[position_] += val;
+ // if (il == 0) printf("FOUND %d %d %d\n",il,ic,position_);
+ return;
+ }
+ if (ptr[position_] == 0)break;
+ position_ = ptr[position_];
+ }
+ }
+
+ INDEX_TYPE zero = 0;
+ CSRList_Add (a_, &val);
+ CSRList_Add (ai_, &ic);
+ CSRList_Add (ptr_, &zero);
+ // The pointers may have been modified
+ // if there has been a reallocation in CSRList_Add
+
+ ptr = (INDEX_TYPE*) ptr_->array;
+ ai = (INDEX_TYPE*) ai_->array;
+ a = (scalar*) a_->array;
+
+ INDEX_TYPE n = CSRList_Nbr(a_) - 1;
+
+ // if (il == 0) printf("NOT FOUND %d %d %d\n",il,ic,n);
+
+ if(!something[il]) {
+ jptr[il] = n;
+ something[il] = 1;
+ }
+ else ptr[position_] = n;
+
+/* for (int i=0;i<_b->size();i++)printf("%d ",something[i]?jptr[i]:-1); */
+/* printf("\n"); */
+/* for (int i=0;i<CSRList_Nbr(ai_);i++)printf("(%d %d %g)",ai[i],ptr[i],a[i]); */
+/* printf("\n"); */
+
+ }
+
+ virtual scalar getFromMatrix (int _row, int _col) const
+ {
+ throw;
+ }
+ virtual void addToRightHandSide(int _row, scalar _val)
+ {
+ if(_val != 0.0) (*_b)[_row] += _val;
+ }
+ virtual scalar getFromRightHandSide(int _row) const
+ {
+ return (*_b)[_row];
+ }
+ virtual scalar getFromSolution(int _row) const
+ {
+ return (*_x)[_row];
+ }
+ virtual void zeroMatrix()
+ {
+ int N=CSRList_Nbr(a_);
+ scalar *a = (scalar*) a_->array;
+ for (int i=0;i<N;i++)a[i]=0;
+ }
+ virtual void zeroRightHandSide()
+ {
+ for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+ }
+};
+
+
+template <class scalar>
+class gmshLinearSystemCSRGmm : public gmshLinearSystemCSR<scalar> {
+ private:
+ double _prec;
+ int _noisy, _gmres;
+ public:
+ gmshLinearSystemCSRGmm()
+ : _prec(1.e-8), _noisy(0), _gmres(0) {}
+ virtual ~gmshLinearSystemCSRGmm()
+ {}
+ void setPrec(double p){ _prec = p; }
+ void setNoisy(int n){ _noisy = n; }
+ void setGmres(int n){ _gmres = n; }
+ virtual int systemSolve()
+#if defined(HAVE_GMM)
+ ;
+#else
+ {
+ Msg::Error("Gmm++ is not available in this version of Gmsh");
+ }
+#endif
+};
+
+template <class scalar>
+class gmshLinearSystemCSRTaucs : public gmshLinearSystemCSR<scalar> {
+ private:
+ public:
+ gmshLinearSystemCSRTaucs()
+ {}
+ virtual ~gmshLinearSystemCSRTaucs()
+ {}
+ virtual int systemSolve()
+#if defined(HAVE_TAUCS)
+ ;
+#else
+ {
+ Msg::Error("Taucs is not available in this version of Gmsh");
+ }
+#endif
+};
+
+
+#endif