diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index eccdc6bf93a18c5287aac8c55b7ab39bf017fe02..460d0099b0e13c8f9239f45b15fbcb25609e3bb6 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1101,7 +1101,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value;
- myAssembler.getDofValue(value, v, 0, 1);
+ myAssembler.getDofValue(v, 0, 1, value);
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
@@ -1350,8 +1350,8 @@ bool GFaceCompound::parametrize_conformal() const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value1, value2;
- myAssembler.getDofValue(value1, v,0,1);
- myAssembler.getDofValue(value2, v,0,2);
+ myAssembler.getDofValue(v, 0, 1, value1);
+ myAssembler.getDofValue(v, 0, 2, value2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 7822b83d7ec5b8c598ebc8b526a24a7a3b76b538..8f371dccefc3c5790481a18bd10fd353a6b98d25 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -371,7 +371,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
+ myAssembler.getDofValue(v_3D, 0, 1, _sizes[v_2D]);
}
delete _lsys;
#endif
@@ -472,7 +472,7 @@ void backgroundMesh::updateSizes(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
+ myAssembler.getDofValue(v_3D, 0, 1, _sizes[v_2D]);
}
delete _lsys;
#endif
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 281cde696c9e9b889d18c708d75074e97ffd82e5..15ef36313e750fcd2ae1726bcc147b12f3f2190c 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -600,8 +600,8 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
SPoint2 param;
reparamMeshVertexOnFace((*it), gf, param);
SPoint2 dparam;
- myAssembler.getDofValue(dparam[0], (*it), 0, getTag());
- myAssembler.getDofValue(dparam[1], (*it), 1, getTag());
+ myAssembler.getDofValue((*it), 0, getTag(), dparam[0]);
+ myAssembler.getDofValue((*it), 1, getTag(), dparam[1]);
SPoint2 newp = param+dparam;
dx += newp.x() * newp.x() + newp.y() * newp.y();
(*it)->setParameter(0, newp.x());
@@ -717,9 +717,9 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
double ax, ay, az;
- myAssembler.getDofValue(ax, it->first, 0, getTag());
- myAssembler.getDofValue(ay, it->first, 1, getTag());
- myAssembler.getDofValue(az, it->first, 2, getTag());
+ myAssembler.getDofValue(it->first, 0, getTag(), ax);
+ myAssembler.getDofValue(it->first, 1, getTag(), ay);
+ myAssembler.getDofValue(it->first, 2, getTag(), az);
it->first->x() += ax;
it->first->y() += ay;
it->first->z() += az;
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 04630c7ae53d15b0cb8022530662ef3844a7934e..01ac94f1a846f656b7fb42d8e735048942455fc7 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -22,6 +22,9 @@ class fullVector
scalar *_data;
friend class fullMatrix<scalar>;
public:
+ inline const scalar* getDataPtr() const{
+ return _data;
+ }
fullVector(int r) : _r(r)
{
_data = new scalar[_r];
@@ -204,6 +207,12 @@ class fullMatrix
}
return *this;
}
+ void operator += (const fullMatrix<scalar> &other)
+ {
+ if(_r != other._r || _c!= other._c)
+ Msg::Error("sum matrices of different sizes\n");
+ for(int i = 0; i < _r * _c; ++i) _data[i] += other._data[i];
+ }
inline scalar operator () (int i, int j) const
{
return _data[i + _r * j];
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index bb5634416319a46bdf538dbb9149f025e03760b3..cf67e2addafc2f706a55ca8c0cdcdab321fdb265 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -189,7 +189,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
itv !=distance_map2.end() ; ++itv){
MVertex *v = itv->first;
double value;
- myAssembler.getDofValue(value, v, 0, 1);
+ myAssembler.getDofValue(v, 0, 1, value);
value = std::min(0.9999, value);
double dist = -mu * log(1. - value);
itv->second = dist;
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 942a375092395e1364b4ce628c27d497ad06661f..2c1c0126f5103ca2efe5524df47999fd9130f3b5 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -50,7 +50,7 @@ template<class T> struct dofTraits
{
typedef T VecType;
typedef T MatType;
- inline static void mult (VecType &r, const MatType &m, const VecType &v) { r = m*v;}
+ inline static void gemm (VecType &r, const MatType &m, const VecType &v) { r += m*v;}
inline static void neg (VecType &r) { r = -r;}
inline static void setToZero (VecType &r) { r = 0;}
};
@@ -59,7 +59,7 @@ template<class T> struct dofTraits<fullMatrix<T> >
{
typedef fullMatrix<T> VecType;
typedef fullMatrix<T> MatType;
- inline static void mult (VecType &r, const MatType &m, const VecType &v) { m.mult(v, r);}
+ inline static void gemm (VecType &r, const MatType &m, const VecType &v) { r.gemm(m, v);}
inline static void neg (VecType &r) { r.scale(-1.);}
inline static void setToZero (VecType &r) { r.scale(0);}
};
@@ -158,39 +158,47 @@ class dofManager{
inline void getDofValue(std::vector<Dof> &keys,std::vector<dataVec> &Vals)
{
int ndofs=keys.size();
- Vals.reserve(Vals.size()+ndofs);
- for (int i=0;i<ndofs;++i) Vals.push_back(getDofValue(keys[i]));
+ size_t originalSize = Vals.size();
+ Vals.resize(originalSize+ndofs);
+ for (int i=0;i<ndofs;++i) getDofValue(keys[i], Vals[originalSize+i]);
}
- inline dataVec getDofValue(Dof key) const
+ inline void getDofValue(Dof key, dataVec &val) const
{
{
typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
- if (it != fixed.end()) return it->second;
+ if (it != fixed.end()) {
+ val = it->second;
+ return ;
+ }
}
{
std::map<Dof, int>::const_iterator it = unknown.find(key);
- if (it != unknown.end())
- return _current->getFromSolution(it->second);
+ if (it != unknown.end()) {
+ _current->getFromSolution(it->second,val );
+ return;
+ }
}
{
typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it = constraints.find(key);
if (it != constraints.end())
{
- dataVec somme(0.0);
+ dataVec tmp(val);
+ dofTraits<T>::setToZero(val);
for (int i=0;i<(it->second).linear.size();i++)
{
std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
- somme = somme + getDofValue(((it->second).linear[i]).first)*((it->second).linear[i]).second;
+ getDofValue(((it->second).linear[i]).first, tmp);
+ dofTraits<T>::gemm(val,((it->second).linear[i]).second, tmp);
}
- somme = somme + (it->second).shift;
- return somme;
+ val += (it->second).shift;
+ return ;
}
}
- return dataVec(0.0);
+ dofTraits<T>::setToZero(val);
}
- inline void getDofValue(dataVec &v, int ent, int type) const
+ inline void getDofValue(int ent, int type, dataVec &v) const
{
Dof key(ent, type);
{
@@ -203,7 +211,7 @@ class dofManager{
{
std::map<Dof, int>::const_iterator it = unknown.find(key);
if (it != unknown.end()){
- v = _current->getFromSolution(it->second);
+ _current->getFromSolution(it->second, v);
return;
}
}
@@ -212,22 +220,23 @@ class dofManager{
if (it != constraints.end())
{
dofTraits<T>::setToZero(v);
+ dataVec tmp(v);
for (int i=0;i<(it->second).linear.size();i++)
{
std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
- v = v + getDofValue(((it->second).linear[i]).first)*((it->second).linear[i]).second;
+ getDofValue(((it->second).linear[i]).first, tmp);
+ dofTraits<T>::gemm(v, ((it->second).linear[i]).second, tmp);
}
- v = v + (it->second).shift;
+ v += (it->second).shift;
return ;
}
}
dofTraits<T>::setToZero(v);
}
-
- inline dataVec getDofValue(dataVec &value, MVertex *v, int iComp, int iField) const
+ inline void getDofValue(MVertex *v, int iComp, int iField, dataVec &value) const
{
- getDofValue(value, v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
+ getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
}
inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
@@ -244,7 +253,8 @@ class dofManager{
if (itFixed != fixed.end()) {
// tmp = -value * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::mult(tmp , value, itFixed->second);
+ dofTraits<T>::setToZero(tmp);
+ dofTraits<T>::gemm(tmp , value, itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(itR->second, tmp);
}else assembleLinConst(R, C, value);
@@ -255,9 +265,7 @@ class dofManager{
assembleLinConst(R, C, value);
}
}
-
-
- inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C, fullMatrix<dataMat> &m)
+ inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C, const fullMatrix<dataMat> &m)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
@@ -284,7 +292,8 @@ class dofManager{
if (itFixed != fixed.end()){
// tmp = -m(i,j) * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::setToZero(tmp);
+ dofTraits<T>::gemm(tmp , m(i,j), itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(NR[i], tmp);
}else assembleLinConst(R[i], C[j], m(i,j));
@@ -301,7 +310,7 @@ class dofManager{
}
}
- inline void assemble(std::vector<Dof> &R, fullVector<dataMat> &m) // for linear forms
+ inline void assemble(std::vector<Dof> &R, const fullVector<dataMat> &m) // for linear forms
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
std::vector<int> NR(R.size());
@@ -334,7 +343,7 @@ class dofManager{
}
- inline void assemble(std::vector<Dof> &R, fullMatrix<dataMat> &m)
+ inline void assemble(std::vector<Dof> &R, const fullMatrix<dataMat> &m)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
std::vector<int> NR(R.size());
@@ -360,7 +369,8 @@ class dofManager{
if (itFixed != fixed.end()){
// tmp = -m(i,j) * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::setToZero(tmp);
+ dofTraits<T>::gemm(tmp , m(i,j), itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(NR[i], tmp);
} else assembleLinConst(R[i],R[j],m(i,j));
@@ -401,7 +411,7 @@ class dofManager{
assemble(Dof(entR, typeR), value);
}
inline void assemble(MVertex *vR, int iCompR, int iFieldR,
- const dataMat &value)
+ const dataMat &value)
{
assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR), value);
}
@@ -446,11 +456,18 @@ class dofManager{
itConstraint = constraints.find(C);
if (itConstraint != constraints.end())
{
+ dataMat tmp(value);
for (int i=0;i<(itConstraint->second).linear.size();i++)
{
- assemble(R,(itConstraint->second).linear[i].first,(itConstraint->second).linear[i].second*value);
+ dofTraits<T>::setToZero(tmp);
+ dofTraits<T>::gemm(tmp,(itConstraint->second).linear[i].second,value);
+ assemble(R,(itConstraint->second).linear[i].first,tmp);
}
- _current->addToRightHandSide(itR->second, -value * (itConstraint->second).shift);
+ dataMat tmp2(value);
+ dofTraits<T>::setToZero(tmp2);
+ dofTraits<T>::gemm(tmp2, value, itConstraint->second.shift);
+ dofTraits<T>::neg(tmp2);
+ _current->addToRightHandSide(itR->second, tmp2);
}
}
else // test function ; (no shift ?)
@@ -459,9 +476,12 @@ class dofManager{
itConstraint = constraints.find(R);
if (itConstraint != constraints.end())
{
+ dataMat tmp(value);
for (int i=0;i<(itConstraint->second).linear.size();i++)
{
- assemble((itConstraint->second).linear[i].first,C,(itConstraint->second).linear[i].second*value);
+ dofTraits<T>::setToZero(tmp);
+ dofTraits<T>::gemm(tmp,itConstraint->second.linear[i].second,value);
+ assemble((itConstraint->second).linear[i].first,C,tmp);
}
}
}
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 65540b388c94fbea3910aa7c8b355ceeafd2d324..a4a70c44bcd8d52ec4af757ec3298416230d9b46 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -225,9 +225,9 @@ static double vonMises(dofManager<double> *a, MElement *e,
double valy[256];
double valz[256];
for (int k = 0; k < e->getNumVertices(); k++){
- a->getDofValue(valx[k], e->getVertex(k), 0, _tag);
- a->getDofValue(valy[k], e->getVertex(k), 1, _tag);
- a->getDofValue(valz[k], e->getVertex(k), 2, _tag);
+ a->getDofValue(e->getVertex(k), 0, _tag, valx[k]);
+ a->getDofValue(e->getVertex(k), 1, _tag, valy[k]);
+ a->getDofValue(e->getVertex(k), 2, _tag, valz[k]);
}
double gradux[3];
double graduy[3];
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index 3b6f9982ed784ff7ab7fc4acbe4e51ec9d26ef6c..5ff9c10a148ecfab0d5a6a24a8ce31a27382839c 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -18,11 +18,11 @@ class linearSystem {
virtual bool isAllocated() const = 0;
virtual void allocate(int nbRows) = 0;
virtual void clear() = 0;
- virtual void addToMatrix(int _row, int _col, scalar val) = 0;
- virtual scalar getFromMatrix(int _row, int _col) const = 0;
- virtual void addToRightHandSide(int _row, scalar val) = 0;
- virtual scalar getFromRightHandSide(int _row) const = 0;
- virtual scalar getFromSolution(int _row) const = 0;
+ virtual void addToMatrix(int _row, int _col, const scalar &val) = 0;
+ virtual void getFromMatrix(int _row, int _col, scalar &val) const = 0;
+ virtual void addToRightHandSide(int _row, const scalar &val) = 0;
+ virtual void getFromRightHandSide(int _row, scalar &val) const = 0;
+ virtual void getFromSolution(int _row, scalar &val) const = 0;
virtual void zeroMatrix() = 0;
virtual void zeroRightHandSide() = 0;
virtual int systemSolve() = 0;
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 207a1060ac08c92320656b9f70da9f3825f08690..a86ef09e1cad4d0c96dd70004eaaa94cac7871f9 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -74,7 +74,7 @@ static void CSRList_Delete(CSRList_T *liste)
}
}
-void CSRList_Add(CSRList_T *liste, void *data)
+void CSRList_Add(CSRList_T *liste,const void *data)
{
liste->n++;
CSRList_Realloc(liste, liste->n);
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index ed7e369596d353822b57c8e73ade13c6efcab3cf..4dc0f47ceb931333a73881a7a0155f15438e9bfa 100644
--- a/Solver/linearSystemCSR.h
+++ b/Solver/linearSystemCSR.h
@@ -23,7 +23,7 @@ typedef struct {
char *array;
} CSRList_T;
-void CSRList_Add(CSRList_T *liste, void *data);
+void CSRList_Add(CSRList_T *liste, const void *data);
int CSRList_Nbr(CSRList_T *liste);
template <class scalar>
@@ -46,7 +46,7 @@ class linearSystemCSR : public linearSystem<scalar> {
{
allocate(0);
}
- virtual void addToMatrix(int il, int ic, double val)
+ virtual void addToMatrix(int il, int ic, const double &val)
{
INDEX_TYPE *jptr = (INDEX_TYPE*) _jptr->array;
INDEX_TYPE *ptr = (INDEX_TYPE*) _ptr->array;
@@ -87,22 +87,21 @@ class linearSystemCSR : public linearSystem<scalar> {
}
virtual void getMatrix(INDEX_TYPE*& jptr,INDEX_TYPE*& ai,double*& a);
- virtual scalar getFromMatrix (int row, int col) const
+ virtual void getFromMatrix (int row, int col, scalar &val) const
{
Msg::Error("getFromMatrix not implemented for CSR");
- return 0;
}
- virtual void addToRightHandSide(int row, scalar val)
+ virtual void addToRightHandSide(int row, const scalar &val)
{
if(val != 0.0) (*_b)[row] += val;
}
- virtual scalar getFromRightHandSide(int row) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
- return (*_b)[row];
+ val = (*_b)[row];
}
- virtual scalar getFromSolution(int row) const
+ virtual void getFromSolution(int row, scalar &val) const
{
- return (*_x)[row];
+ val = (*_x)[row];
}
virtual void zeroMatrix()
{
@@ -143,7 +142,7 @@ class linearSystemCSRTaucs : public linearSystemCSR<scalar> {
public:
linearSystemCSRTaucs(){}
virtual ~linearSystemCSRTaucs(){}
- virtual void addToMatrix(int il, int ic, double val)
+ virtual void addToMatrix(int il, int ic, const double &val)
{
if (il <= ic)
linearSystemCSR<scalar>::addToMatrix(il, ic, val);
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index ca0fe74a61f49d7950d5fc259559f9b3c45a6eac..586fa52aaa3fa01ba339efea3456d506e72cc065 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -42,25 +42,25 @@ class linearSystemFull : public linearSystem<scalar> {
}
_a = 0;
}
- virtual void addToMatrix(int row, int col, scalar val)
+ virtual void addToMatrix(int row, int col, const scalar &val)
{
if(val != 0.0) (*_a)(row, col) += val;
}
- virtual scalar getFromMatrix(int row, int col) const
+ virtual void getFromMatrix(int row, int col, scalar &val) const
{
- return (*_a)(row, col);
+ val = (*_a)(row, col);
}
- virtual void addToRightHandSide(int row, scalar val)
+ virtual void addToRightHandSide(int row, const scalar &val)
{
if(val != 0.0) (*_b)(row) += val;
}
- virtual scalar getFromRightHandSide(int row) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
- return (*_b)(row);
+ val = (*_b)(row);
}
- virtual scalar getFromSolution(int row) const
+ virtual void getFromSolution(int row, scalar &val) const
{
- return (*_x)(row);
+ val = (*_x)(row);
}
virtual void zeroMatrix()
{
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index d79b0ba5b1975a8c74262ff6516441cf332ef7a1..f67e4c7e99af21c15e807b3d579ec7251f556d68 100644
--- a/Solver/linearSystemGMM.h
+++ b/Solver/linearSystemGMM.h
@@ -46,25 +46,25 @@ class linearSystemGmm : public linearSystem<scalar> {
}
_a = 0;
}
- virtual void addToMatrix(int row, int col, scalar val)
+ virtual void addToMatrix(int row, int col, scalar &val)
{
if(val != 0.0) (*_a)(row, col) += val;
}
- virtual scalar getFromMatrix (int row, int col) const
+ virtual void getFromMatrix (int row, int col, scalar &val) const
{
- return (*_a)(row, col);
+ val = (*_a)(row, col);
}
- virtual void addToRightHandSide(int row, scalar val)
+ virtual void addToRightHandSide(int row, scalar &val)
{
if(val != 0.0) (*_b)[row] += val;
}
- virtual scalar getFromRightHandSide(int row) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
- return (*_b)[row];
+ val = (*_b)[row];
}
- virtual scalar getFromSolution(int row) const
+ virtual void getFromSolution(int row, scalar &val) const
{
- return (*_x)[row];
+ val = (*_x)[row];
}
virtual void zeroMatrix()
{
@@ -100,11 +100,11 @@ class linearSystemGmm : public linearSystem<scalar> {
}
virtual bool isAllocated() const { return false; }
virtual void allocate(int nbRows) {}
- virtual void addToMatrix(int row, int col, scalar val) {}
- virtual scalar getFromMatrix(int row, int col) const { return 0.; }
- virtual void addToRightHandSide(int row, scalar val) {}
- virtual scalar getFromRightHandSide(int row) const { return 0.; }
- virtual scalar getFromSolution(int row) const { return 0.; }
+ virtual void addToMatrix(int row, int col, scalar &val) {}
+ virtual void getFromMatrix(int row, int col, calar &val) const { return 0.; }
+ virtual void addToRightHandSide(int row, scalar &val) {}
+ virtual void getFromRightHandSide(int row, scalar &val) const { return 0.; }
+ virtual void getFromSolution(int row, scalar &val) const { return 0.; }
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual int systemSolve() { return 0; }
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index c77af8d464bb156544d4bc3cf5bde89ac1812911..55d908a0c13becda94e08e46bb7d1d33c66d5aec 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -3,37 +3,44 @@
#include "linearSystemPETSc.h"
#include "fullMatrix.h"
template <>
-void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col, fullMatrix<double> val)
+void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col, const fullMatrix<double> &val)
{
+ fullMatrix<double> &modval = *const_cast<fullMatrix<double> *>(&val);
for (int ii = 0; ii<val.size1(); ii++)
for (int jj = 0; jj < ii; jj ++) {
- double buff = val(ii,jj);
- val(ii,jj) = val (jj,ii);
- val(jj,ii) = buff;
+ double buff = modval(ii,jj);
+ modval(ii,jj) = modval (jj,ii);
+ modval(jj,ii) = buff;
}
PetscInt i = row, j = col;
- _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &val(0,0), ADD_VALUES));
+ _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES));
+ //transpose back so that the original matrix is not modified
+ for (int ii = 0; ii<val.size1(); ii++)
+ for (int jj = 0; jj < ii; jj ++) {
+ double buff = modval(ii,jj);
+ modval(ii,jj) = modval (jj,ii);
+ modval(jj,ii) = buff;
+ }
}
template<>
-void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row, fullMatrix<double> val)
+void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row, const fullMatrix<double> &val)
{
for (int ii = 0; ii < _blockSize; ii ++) {
PetscInt i = row * _blockSize + ii;
- _try(VecSetValues(_b, 1, &i, &val(ii,0), ADD_VALUES));
+ double v = val(ii,0);
+ _try(VecSetValues(_b, 1, &i, &v, ADD_VALUES));
}
}
template<>
-fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromMatrix(int row, int col) const
+void linearSystemPETSc<fullMatrix<double> >::getFromMatrix(int row, int col, fullMatrix<double> &val ) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
- return fullMatrix<double>(0,0);
}
template<>
-fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row) const
+void linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row, fullMatrix<double> &val) const
{
- fullMatrix<double> val(_blockSize,1);
PetscScalar *tmp;
_try(VecGetArray(_b, &tmp));
for (int i=0; i<_blockSize; i++) {
@@ -45,13 +52,11 @@ fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(
#endif
}
_try(VecRestoreArray(_b, &tmp));
- return val;
}
template<>
-fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row) const
+void linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row, fullMatrix<double> &val) const
{
- fullMatrix<double> val(_blockSize,1);
PetscScalar *tmp;
_try(VecGetArray(_x, &tmp));
for (int i=0; i<_blockSize; i++) {
@@ -63,7 +68,6 @@ fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromSolution(int r
#endif
}
_try(VecRestoreArray(_x, &tmp));
- return val;
}
template<>
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index d17a14c333e13d8dd001a04007fa240ac19f6806..fa470e20ef254d9420c7f03d68afb41982aba97e 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -83,18 +83,17 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
_isAllocated = false;
}
- virtual scalar getFromMatrix(int row, int col) const
+ virtual void getFromMatrix(int row, int col, scalar &val) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
- return 0;
}
- virtual void addToRightHandSide(int row, scalar val)
+ virtual void addToRightHandSide(int row, const scalar &val)
{
PetscInt i = row;
PetscScalar s = val;
_try(VecSetValues(_b, 1, &i, &s, ADD_VALUES));
}
- virtual scalar getFromRightHandSide(int row) const
+ virtual void getFromRightHandSide(int row, scalar &val) const
{
PetscScalar *tmp;
_try(VecGetArray(_b, &tmp));
@@ -102,27 +101,27 @@ class linearSystemPETSc : public linearSystem<scalar> {
_try(VecRestoreArray(_b, &tmp));
// FIXME specialize this routine
#if defined(PETSC_USE_COMPLEX)
- return s.real();
+ val = s.real();
#else
- return s;
+ val = s;
#endif
}
- virtual void addToMatrix(int row, int col, scalar val)
+ virtual void addToMatrix(int row, int col, const scalar &val)
{
PetscInt i = row, j = col;
PetscScalar s = val;
_try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
}
- virtual scalar getFromSolution(int row) const
+ virtual void getFromSolution(int row, scalar &val) const
{
PetscScalar *tmp;
_try(VecGetArray(_x, &tmp));
PetscScalar s = tmp[row];
_try(VecRestoreArray(_x, &tmp));
#if defined(PETSC_USE_COMPLEX)
- return s.real();
+ val = s.real();
#else
- return s;
+ val = s;
#endif
}
virtual void zeroMatrix()
@@ -152,9 +151,6 @@ class linearSystemPETSc : public linearSystem<scalar> {
// override the default options with the ones from the option
// database (if any)
_try(KSPSetFromOptions(ksp));
- PetscViewer view;
- PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat.out", &view);
- _try(MatView(_a, view));
_try(KSPSolve(ksp, _b, _x));
_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
PetscInt its;
@@ -180,11 +176,11 @@ class linearSystemPETSc : public linearSystem<scalar> {
virtual bool isAllocated() const { return false; }
virtual void allocate(int nbRows) {}
virtual void clear(){}
- virtual void addToMatrix(int row, int col, scalar val) {}
- virtual scalar getFromMatrix(int row, int col) const { return 0.; }
- virtual void addToRightHandSide(int row, scalar val) {}
- virtual scalar getFromRightHandSide(int row) const { return 0.; }
- virtual scalar getFromSolution(int row) const { return 0.; }
+ virtual void addToMatrix(int row, int col, const scalar &val) {}
+ virtual void getFromMatrix(int row, int col, scalar &val) const { return 0.; }
+ virtual void addToRightHandSide(int row, const scalar &val) {}
+ virtual void getFromRightHandSide(int row, scalar &val) const { return 0.; }
+ virtual void getFromSolution(int row, scalar &val) const { return 0.; }
virtual void zeroMatrix() {}
virtual void zeroRightHandSide() {}
virtual int systemSolve() { return 0; }
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 29acbb5a29e68499d37a2b2866448167afabd3b0..93693142e189eef344407befe9e75f500e1537f7 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1200,7 +1200,7 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value;
- myAssembler.getDofValue(value, v, 0, 1);
+ myAssembler.getDofValue(v, 0, 1, value);
if (step == 0)solution[v] = SPoint2(value,0);
else solution[v] = SPoint2(solution[v][0],value);
}