No commit message

No commit message

Solver/linearSystemCSR.cpp

@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <complex>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "linearSystemCSR.h"
@@ -135,6 +136,49 @@ void linearSystemCSR<double>::preAllocateEntries ()
   }
 }
 
+template<>
+void linearSystemCSR<std::complex<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_strict (_ai, nnz);
+  CSRList_Resize_strict (_ptr, nnz);
+  INDEX_TYPE *jptr = (INDEX_TYPE*) _jptr->array;
+  INDEX_TYPE *ai = (INDEX_TYPE*) _ai->array;
+  INDEX_TYPE *ptr = (INDEX_TYPE*) _ptr->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];
+      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();
+  // we do this after _sparsity.clear so that the peak memory usage is reduced
+  CSRList_Resize_strict (_a, nnz);
+  std::complex<double> *a = ( std::complex<double> * ) _a->array;
+  for (int i = 0; i < nnz; i++) {
+    a[i] = std::complex<double>();
+  }
+}
+
 template<>
 void linearSystemCSR<double>::allocate(int nbRows)
 {
@@ -173,6 +217,44 @@ void linearSystemCSR<double>::allocate(int nbRows)
   _x = new std::vector<double>(nbRows);
 }
 
+template<>
+void linearSystemCSR<std::complex<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;
+    sorted = false;
+    something = 0;
+    return;
+  }
+
+  _a    = CSRList_Create(nbRows, nbRows, sizeof(std::complex<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<std::complex<double> >(nbRows);
+  _x = new std::vector<std::complex<double> >(nbRows);
+}
+
 const int NSTACK = 50;
 const unsigned int M_sort2 = 7;
 
@@ -339,6 +421,18 @@ void linearSystemCSR<double>::getMatrix(INDEX_TYPE*& jptr,INDEX_TYPE*& ai,double
   sorted = true;
 }
 
+template<>
+void linearSystemCSR<std::complex<double> >::getMatrix(INDEX_TYPE*& jptr,INDEX_TYPE*& ai,double *& a)
+{
+  jptr = (INDEX_TYPE*) _jptr->array;
+  ai = (INDEX_TYPE*) _ai->array;
+  a = ( double  * ) _a->array;
+  if (!sorted)
+    sortColumns_(_b->size(), CSRList_Nbr(_a), (INDEX_TYPE *) _ptr->array, jptr,
+                 ai, (std::complex<double> *)a);
+  sorted = true;
+}
+
 #if defined(HAVE_GMM)
 
 #include "gmm.h"
@@ -417,4 +511,43 @@ int linearSystemCSRTaucs<double>::systemSolve()
   return 1;
 }
 
+template<>
+int linearSystemCSRTaucs<std::complex<double> >::systemSolve()
+{
+  if (!_a)return 1;
+  if(!sorted){
+    sortColumns_(_b->size(),
+                CSRList_Nbr(_a),
+                (INDEX_TYPE *) _ptr->array,
+                (INDEX_TYPE *) _jptr->array,
+                (INDEX_TYPE *) _ai->array,
+                (std::complex<double>*) _a->array);
+  }
+  sorted = true;
+
+  taucs_ccs_matrix myVeryCuteTaucsMatrix;
+  myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
+  //myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ptr->array;
+  //myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_ai->array;
+  myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ai->array;
+  myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_jptr->array;
+  myVeryCuteTaucsMatrix.values.z = (taucs_dcomplex*)_a->array;
+  myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DCOMPLEX;
+  char* options[] = { (char*)"taucs.factor.mf=true", NULL };
+  double t1 = Cpu();
+  int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
+                              NULL,
+                              1,
+                              &(*_x)[0],
+                              &(*_b)[0],
+                              options,
+                              NULL);
+  double t2 = Cpu();
+  Msg::Debug("TAUCS has solved %d unknowns in %8.3f seconds", _b->size(), t2 - t1);
+  if (result != TAUCS_SUCCESS){
+    Msg::Error("Taucs Was Not Successfull %d", result);
+  }
+  return 1;
+}
+
 #endif

Solver/linearSystemCSR.h

@@ -12,6 +12,7 @@
 #include "linearSystem.h"
 #include "sparsityPattern.h"
 
+
 typedef int INDEX_TYPE ;
 typedef struct {
   int nmax;
@@ -53,6 +54,7 @@ class linearSystemCSR : public linearSystem<scalar> {
   virtual void preAllocateEntries ();
   virtual void addToMatrix(int il, int ic, const scalar &val) 
   {
+    
     if (!_entriesPreAllocated)
       preAllocateEntries();
     INDEX_TYPE  *jptr  = (INDEX_TYPE*) _jptr->array;
@@ -123,11 +125,11 @@ class linearSystemCSR : public linearSystem<scalar> {
   }
   virtual void addToRightHandSide(int row, const scalar &val)
   {
-    if(val != 0.0) (*_b)[row] += val;
+    if(val != scalar()) (*_b)[row] += val;
   }
   virtual void addToSolution(int row, const scalar &val)
   {
-    if(val != 0.0) (*_x)[row] += val;
+    if(val != scalar()) (*_x)[row] += val;
   }
   virtual void getFromRightHandSide(int row, scalar &val) const
   {
@@ -142,17 +144,17 @@ class linearSystemCSR : public linearSystem<scalar> {
     if (!_a) return;
     int N = CSRList_Nbr(_a);
     scalar *a = (scalar*) _a->array;
-    for (int i = 0; i < N; i++) a[i] = 0;
+    for (int i = 0; i < N; i++) a[i] = scalar();
   }
   virtual void zeroRightHandSide()
   {
     if (!_b) return;
-    for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
+    for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = scalar();
   }
   virtual void zeroSolution()
   {
     if (!_x) return;
-    for(unsigned int i = 0; i < _x->size(); i++) (*_x)[i] = 0.;
+    for(unsigned int i = 0; i < _x->size(); i++) (*_x)[i] = scalar();
   }
 
   virtual double normInfRightHandSide() const
@@ -160,8 +162,7 @@ class linearSystemCSR : public linearSystem<scalar> {
     double nor = 0.;
     double temp;
     for(unsigned int i = 0; i < _b->size(); i++){
-      temp = (*_b)[i];
-      if(temp < 0) temp = -temp;
+      temp = std::abs((*_b)[i]);
       if(nor < temp) nor = temp;
     }
     return nor;
@@ -194,7 +195,7 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
  public:
   linearSystemCSRTaucs(){}
   virtual ~linearSystemCSRTaucs(){}
-  virtual void addToMatrix(int il, int ic, const double &val)
+  virtual void addToMatrix(int il, int ic, const scalar &val)
   {
     if (il <= ic) {
       linearSystemCSR<scalar>::addToMatrix(il, ic, val);