diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 9cb38e7ecd4d9f95dfece6df82994c88a642b894..26700358043419ea93222e958bcb499ad358c459 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/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
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 6c05e03643640ecbaf49b0cea18c2b024f108e0d..52f286600c49daae5e60b6774fefc4ac42dc7d6b 100644
--- a/Solver/linearSystemCSR.h
+++ b/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);