Commit c7cbb878 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Performance optimization

parent 64b92aed
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include "DOFVector.h" #include "DOFVector.h"
#include "BasisFunction.h" #include "BasisFunction.h"
#include "Lagrange.h" #include "Lagrange.h"
#include "OpenMP.h"
namespace AMDiS { namespace AMDiS {
...@@ -14,19 +15,34 @@ namespace AMDiS { ...@@ -14,19 +15,34 @@ namespace AMDiS {
/* are those corresponding to these barycentric coordinates. */ /* are those corresponding to these barycentric coordinates. */
/****************************************************************************/ /****************************************************************************/
BasisFunction::~BasisFunction()
{
DELETE nDOF;
}
BasisFunction::BasisFunction(const ::std::string& name_, int dim_, int degree_) BasisFunction::BasisFunction(const ::std::string& name_, int dim_, int degree_)
: name(name_), degree(degree_), dim(dim_) : name(name_),
degree(degree_),
dim(dim_)
{ {
FUNCNAME("BasisFunction::BasisFunction()"); FUNCNAME("BasisFunction::BasisFunction()");
nDOF = NEW DimVec<int>(dim, DEFAULT_VALUE, -1); nDOF = NEW DimVec<int>(dim, DEFAULT_VALUE, -1);
dow = Global::getGeo(WORLD);
grdTmpVec1.resize(omp_get_max_threads());
grdTmpVec2.resize(omp_get_max_threads());
for (int i = 0; i < omp_get_max_threads(); i++) {
grdTmpVec1[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
grdTmpVec2[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
}
}; };
BasisFunction::~BasisFunction()
{
DELETE nDOF;
for (int i = 0; i < grdTmpVec1.size(); i++) {
DELETE grdTmpVec1[i];
DELETE grdTmpVec2[i];
}
}
/****************************************************************************/ /****************************************************************************/
/* some routines for evaluation of a finite element function, its gradient */ /* some routines for evaluation of a finite element function, its gradient */
...@@ -50,9 +66,8 @@ namespace AMDiS { ...@@ -50,9 +66,8 @@ namespace AMDiS {
const WorldVector<double> *uh_loc, const WorldVector<double> *uh_loc,
WorldVector<double>* values) const WorldVector<double>* values) const
{ {
static WorldVector<double> Values(DEFAULT_VALUE, 0.); static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
WorldVector<double> *val = (NULL != values) ? values : &Values; WorldVector<double> *val = (NULL != values) ? values : &Values;
int dow = Global::getGeo(WORLD);
for (int n = 0; n < dow; n++) for (int n = 0; n < dow; n++)
(*val)[n] = 0; (*val)[n] = 0;
...@@ -70,32 +85,27 @@ namespace AMDiS { ...@@ -70,32 +85,27 @@ namespace AMDiS {
const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda, const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda, const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, const double *uh_loc,
WorldVector<double>* grd_uh) const WorldVector<double>* val) const
{ {
static WorldVector<double> grd; TEST_EXIT_DBG(val)("return value is NULL\n");
DimVec<double> grd_b(dim, DEFAULT_VALUE, 0.0);
WorldVector<double> *val;
DimVec<double> grd1(dim, DEFAULT_VALUE, 0.);
val = grd_uh ? grd_uh : &grd; DimVec<double> *grdTmp1 = grdTmpVec1[omp_get_thread_num()];
DimVec<double> *grdTmp2 = grdTmpVec2[omp_get_thread_num()];
for (int j = 0; j < dim + 1; j++) for (int j = 0; j < dim + 1; j++)
grd1[j] = 0.0; (*grdTmp2)[j] = 0.0;
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
(*(*grdPhi)[i])(lambda, grd_b); (*(*grdPhi)[i])(lambda, *grdTmp1);
for (int j = 0; j < dim + 1; j++) for (int j = 0; j < dim + 1; j++)
grd1[j] += uh_loc[i] * grd_b[j]; (*grdTmp2)[j] += uh_loc[i] * (*grdTmp1)[j];
} }
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++) { for (int i = 0; i < dow; i++) {
(*val)[i] = 0; (*val)[i] = 0.0;
for (int j = 0; j < dim + 1; j++) for (int j = 0; j < dim + 1; j++)
(*val)[i] += grd_lambda[j][i] * grd1[j]; (*val)[i] += grd_lambda[j][i] * (*grdTmp2)[j];
} }
return ((*val)); return ((*val));
...@@ -105,7 +115,7 @@ namespace AMDiS { ...@@ -105,7 +115,7 @@ namespace AMDiS {
const DimVec<WorldVector<double> >& grd_lambda, const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, WorldMatrix<double>* D2_uh) const const double *uh_loc, WorldMatrix<double>* D2_uh) const
{ {
static WorldMatrix<double> D2(DEFAULT_VALUE, 0.); static WorldMatrix<double> D2(DEFAULT_VALUE, 0.0);
DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0); DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0); DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
WorldMatrix<double> *val = D2_uh ? D2_uh : &D2; WorldMatrix<double> *val = D2_uh ? D2_uh : &D2;
...@@ -117,8 +127,6 @@ namespace AMDiS { ...@@ -117,8 +127,6 @@ namespace AMDiS {
D2_tmp[k][l] += uh_loc[i] * D2_b[k][l]; D2_tmp[k][l] += uh_loc[i] * D2_b[k][l];
} }
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++) for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) { for (int j = 0; j < dow; j++) {
(*val)[i][j] = 0.0; (*val)[i][j] = 0.0;
......
...@@ -432,6 +432,11 @@ namespace AMDiS { ...@@ -432,6 +432,11 @@ namespace AMDiS {
*/ */
int dim; int dim;
/** \brief
* Dimension of the world.
*/
int dow;
/** \brief /** \brief
* Number of DOFs at the different positions * Number of DOFs at the different positions
*/ */
...@@ -451,6 +456,19 @@ namespace AMDiS { ...@@ -451,6 +456,19 @@ namespace AMDiS {
* Vector of second derivatives * Vector of second derivatives
*/ */
::std::vector<D2BasFctType*> *d2Phi; ::std::vector<D2BasFctType*> *d2Phi;
/** \brief
* Is used by function evalGrdUh. To make it thread safe, we need a
* temporary DimVec for every thread.
*/
::std::vector<DimVec<double>* > grdTmpVec1;
/** \brief
* Is used by function evalGrdUh. To make it thread safe, we need a
* temporary DimVec for every thread.
*/
::std::vector<DimVec<double>* > grdTmpVec2;
}; };
} }
......
...@@ -53,7 +53,7 @@ namespace AMDiS { ...@@ -53,7 +53,7 @@ namespace AMDiS {
// define result vector // define result vector
static DOFVector<WorldVector<double> > *result = NULL; static DOFVector<WorldVector<double> > *result = NULL;
if(grad) { if (grad) {
result = grad; result = grad;
} else { } else {
if(result && result->getFESpace() != feSpace) { if(result && result->getFESpace() != feSpace) {
...@@ -62,8 +62,6 @@ namespace AMDiS { ...@@ -62,8 +62,6 @@ namespace AMDiS {
} }
} }
int i, j;
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim(); int dim = mesh->getDim();
...@@ -80,11 +78,11 @@ namespace AMDiS { ...@@ -80,11 +78,11 @@ namespace AMDiS {
int numNodes = 0; int numNodes = 0;
int numDOFs = 0; int numDOFs = 0;
for(i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
GeoIndex geoIndex = INDEX_OF_DIM(i, dim); GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
int numPositionNodes = mesh->getGeo(geoIndex); int numPositionNodes = mesh->getGeo(geoIndex);
int numPreDOFs = admin->getNumberOfPreDOFs(i); int numPreDOFs = admin->getNumberOfPreDOFs(i);
for(j = 0; j < numPositionNodes; j++) { for (int j = 0; j < numPositionNodes; j++) {
int dofs = basFcts->getNumberOfDOFs(geoIndex); int dofs = basFcts->getNumberOfDOFs(geoIndex);
numNodeDOFs.push_back(dofs); numNodeDOFs.push_back(dofs);
numDOFs += dofs; numDOFs += dofs;
...@@ -99,7 +97,7 @@ namespace AMDiS { ...@@ -99,7 +97,7 @@ namespace AMDiS {
TEST_EXIT_DBG(numDOFs == basFcts->getNumber()) TEST_EXIT_DBG(numDOFs == basFcts->getNumber())
("number of dofs != number of basis functions\n"); ("number of dofs != number of basis functions\n");
for(i = 0; i < numDOFs; i++) { for (int i = 0; i < numDOFs; i++) {
bary.push_back(basFcts->getCoords(i)); bary.push_back(basFcts->getCoords(i));
} }
...@@ -112,22 +110,19 @@ namespace AMDiS { ...@@ -112,22 +110,19 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag); ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while(elInfo) { while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL); const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
int localDOFNr = 0; int localDOFNr = 0;
for(i = 0; i < numNodes; i++) { // for all nodes for (int i = 0; i < numNodes; i++) { // for all nodes
for(j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node for (int j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j]; DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
if(!visited[dofIndex]) { if (!visited[dofIndex]) {
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
result[dofIndex] = basFcts->evalGrdUh(*(bary[localDOFNr]), localUh, &((*result)[dofIndex]));
grdLambda,
localUh,
NULL);
visited[dofIndex] = true; visited[dofIndex] = true;
} }
...@@ -168,8 +163,6 @@ namespace AMDiS { ...@@ -168,8 +163,6 @@ namespace AMDiS {
DOFVector<double> volume(feSpace, "volume"); DOFVector<double> volume(feSpace, "volume");
volume.set(0.0); volume.set(0.0);
int i;
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim(); int dim = mesh->getDim();
...@@ -190,17 +183,15 @@ namespace AMDiS { ...@@ -190,17 +183,15 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag); ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while(elInfo) { while (elInfo) {
double det = elInfo->getDet(); double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL); const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
const WorldVector<double> &grd = basFcts->evalGrdUh(bary, WorldVector<double> grd;
grdLambda, basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
localUh,
NULL);
for(i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
DegreeOfFreedom dofIndex = dof[i][numPreDOFs]; DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
(*result)[dofIndex] += grd * det; (*result)[dofIndex] += grd * det;
volume[dofIndex] += det; volume[dofIndex] += det;
...@@ -212,8 +203,8 @@ namespace AMDiS { ...@@ -212,8 +203,8 @@ namespace AMDiS {
DOFVector<double>::Iterator volIt(&volume, USED_DOFS); DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS); DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
for(volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) { for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
if(*volIt != 0.0) { if (*volIt != 0.0) {
*grdIt *= 1.0/(*volIt); *grdIt *= 1.0/(*volIt);
} }
} }
...@@ -248,7 +239,6 @@ namespace AMDiS { ...@@ -248,7 +239,6 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints(); int numPoints = quadrature->getNumPoints();
int nBasFcts = basFcts->getNumber();
static WorldVector<double> *grd = NULL; static WorldVector<double> *grd = NULL;
...@@ -349,7 +339,6 @@ namespace AMDiS { ...@@ -349,7 +339,6 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints(); int numPoints = quadrature->getNumPoints();
int nBasFcts = basFcts->getNumber();
int i, j, k, l, iq; int i, j, k, l, iq;
static WorldMatrix<double> *vec = NULL; static WorldMatrix<double> *vec = NULL;
...@@ -403,7 +392,6 @@ namespace AMDiS { ...@@ -403,7 +392,6 @@ namespace AMDiS {
for (i = 0; i < nBasFcts; i++) { for (i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n"); WARNING("not tested after index correction\n");
//(*(basFcts->getD2Phi(j)))(quad->getLambda(i), D2Phi);
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi); (*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi);
for (k = 0; k < parts; k++) for (k = 0; k < parts; k++)
...@@ -441,7 +429,7 @@ namespace AMDiS { ...@@ -441,7 +429,7 @@ namespace AMDiS {
this->set(0.0); this->set(0.0);
DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, nBasisFcts); DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts); double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts);
if (feSpace->getMesh() == sourceFeSpace->getMesh()) { if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
...@@ -454,13 +442,13 @@ namespace AMDiS { ...@@ -454,13 +442,13 @@ namespace AMDiS {
while (elInfo) { while (elInfo) {
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); basisFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
source->getLocalVector(el, sourceLocalCoeffs); source->getLocalVector(el, sourceLocalCoeffs);
for (int i = 0; i < nBasisFcts; i++) { for (int i = 0; i < nBasisFcts; i++) {
if (vec[localIndices[i]] == 0.0) { if (vec[myLocalIndices[i]] == 0.0) {
coords = basisFcts->getCoords(i); coords = basisFcts->getCoords(i);
vec[localIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor; vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
} }
} }
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
...@@ -483,12 +471,12 @@ namespace AMDiS { ...@@ -483,12 +471,12 @@ namespace AMDiS {
while (nextTraverse) { while (nextTraverse) {
basisFcts->getLocalIndices(elInfo1->getElement(), basisFcts->getLocalIndices(elInfo1->getElement(),
feSpace->getAdmin(), feSpace->getAdmin(),
localIndices); myLocalIndices);
source->getLocalVector(elInfo2->getElement(), source->getLocalVector(elInfo2->getElement(),
sourceLocalCoeffs); sourceLocalCoeffs);
for (int i = 0; i < nBasisFcts; i++) { for (int i = 0; i < nBasisFcts; i++) {
if (vec[localIndices[i]] == 0.0) { if (vec[myLocalIndices[i]] == 0.0) {
coords1 = basisFcts->getCoords(i); coords1 = basisFcts->getCoords(i);
elInfo1->coordToWorld(*coords1, &worldVec); elInfo1->coordToWorld(*coords1, &worldVec);
elInfo2->worldToCoord(worldVec, &coords2); elInfo2->worldToCoord(worldVec, &coords2);
...@@ -502,7 +490,7 @@ namespace AMDiS { ...@@ -502,7 +490,7 @@ namespace AMDiS {
} }
if (isPositive) { if (isPositive) {
vec[localIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs); vec[myLocalIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);
} }
} }
} }
...@@ -512,7 +500,6 @@ namespace AMDiS { ...@@ -512,7 +500,6 @@ namespace AMDiS {
} }
} }
FREE_MEMORY(localIndices, DegreeOfFreedom, nBasisFcts);
FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts); FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts);
} }
...@@ -539,7 +526,7 @@ namespace AMDiS { ...@@ -539,7 +526,7 @@ namespace AMDiS {
int vNumBasFcts = vBasFcts->getNumber(); int vNumBasFcts = vBasFcts->getNumber();
if (feSpace->getMesh() == vFESpace->getMesh()) { if (feSpace->getMesh() == vFESpace->getMesh()) {
DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts); DegreeOfFreedom *myLocalIndices = localIndices[omp_get_thread_num()];
WorldVector<double> *vLocalCoeffs = NEW WorldVector<double>[vNumBasFcts]; WorldVector<double> *vLocalCoeffs = NEW WorldVector<double>[vNumBasFcts];
Mesh *mesh = feSpace->getMesh(); Mesh *mesh = feSpace->getMesh();
TraverseStack stack; TraverseStack stack;
...@@ -550,20 +537,19 @@ namespace AMDiS { ...@@ -550,20 +537,19 @@ namespace AMDiS {
while (elInfo) { while (elInfo) {
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); basFcts->getLocalIndices(el, feSpace->getAdmin(), myLocalIndices);
v->getLocalVector(el, vLocalCoeffs); v->getLocalVector(el, vLocalCoeffs);
for (int i = 0; i < numBasFcts; i++) { for (int i = 0; i < numBasFcts; i++) {
if (vec[localIndices[i]] == nul) { if (vec[myLocalIndices[i]] == nul) {
coords = basFcts->getCoords(i); coords = basFcts->getCoords(i);
vec[localIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor; vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
} }
} }
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
DELETE [] vLocalCoeffs; DELETE [] vLocalCoeffs;
} else { } else {
ERROR_EXIT("not yet for dual traverse\n"); ERROR_EXIT("not yet for dual traverse\n");
...@@ -651,10 +637,8 @@ namespace AMDiS { ...@@ -651,10 +637,8 @@ namespace AMDiS {
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j]; DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
if (!visited[dofIndex]) { if (!visited[dofIndex]) {
grd = basFcts->evalGrdUh(*(bary[localDOFNr]), basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
grdLambda, localUh, &grd);
localUh,
NULL);
for (int k = 0; k < dow; k++) { for (int k = 0; k < dow; k++) {
(*result)[k][dofIndex] = grd[k]; (*result)[k][dofIndex] = grd[k];
......
...@@ -38,6 +38,7 @@ ...@@ -38,6 +38,7 @@
#include "CreatorInterface.h" #include "CreatorInterface.h"
#include "Serializable.h" #include "Serializable.h"
#include "DOFMatrix.h" #include "DOFMatrix.h"
#include "BasisFunction.h"
#include <vector> #include <vector>
#include <memory> #include <memory>
#include <list> #include <list>
...@@ -69,18 +70,17 @@ namespace AMDiS { ...@@ -69,18 +70,17 @@ namespace AMDiS {
class DOFVectorBase : public DOFIndexed<T> class DOFVectorBase : public DOFIndexed<T>
{ {
public: public:
DOFVectorBase() DOFVectorBase()
: feSpace(NULL), : feSpace(NULL),
elementVector(NULL), elementVector(NULL),
boundaryManager(NULL) boundaryManager(NULL),
nBasFcts(0)
{}; {};
DOFVectorBase(const FiniteElemSpace *f, ::std::string n) DOFVectorBase(const FiniteElemSpace *f, ::std::string n);
: feSpace(f),
name(n), virtual ~DOFVectorBase();
elementVector(NULL),
boundaryManager(NULL)
{};
virtual const T *getLocalVector(const Element *el, T* localVec) const; virtual const T *getLocalVector(const Element *el, T* localVec) const;
...@@ -156,9 +156,14 @@ namespace AMDiS { ...@@ -156,9 +156,14 @@ namespace AMDiS {
*/ */
T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind); T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);
inline ::std::vector<Operator*>& getOperators() { return operators; }; inline ::std::vector<Operator*>& getOperators() {
return operators;
};
inline ::std::vector<double*>& getOperatorFactor() {