Commit 3da3ef08 authored by Jonathan Wunderlich's avatar Jonathan Wunderlich
Browse files

Merge branch '326-refactor-eigensolver' into 'feature'

[326-refactor-eigensolver] use unique_ptr in Vectors, added new constructor to...

Closes #326

See merge request !309
parents d8f4b327 e358a3b5
Pipeline #168325 failed with stages
in 23 minutes and 36 seconds
......@@ -71,6 +71,10 @@ public:
z[i] = v[i];
}
/// selects part of v (indices included)
RVectorT(const RVectorT &v, int startIdx, int endIdx)
: z(v.begin() + startIdx, v.begin() + endIdx), N(z.size()) {}
template<typename REAL1, typename Constraint = mpp_ba::enable_if_constructible<REAL,
REAL1,
REAL>>
......
......@@ -55,8 +55,7 @@ Vector::Vector(const Vector &u) : VectorMatrixBase(u),
Vector::Vector(Scalar b, const Vector &u) : VectorMatrixBase(u),
data(b, u.size()), dirichletFlags(u.dirichletFlags) {
if (mpp_ba::isNearZero(b)) (*this).accumulateFlag = false;
else (*this).accumulateFlag = true;
(*this).accumulateFlag = true;
}
Vector::Vector(const constAB<Operator, Vector> &Ov) :
......@@ -86,8 +85,9 @@ Vector &Vector::operator=(const Vector &u) {
(*this)(r, i) = u(rr, i);
}
return *this;
} else if (u.size() < size()) THROW(
"Not Implemented - since not all entries of '*this' are defined. Use '+=' ");
} else if (u.size() < size())
THROW(
"Not Implemented - since not all entries of '*this' are defined. Use '+=' ");
return *this;
}
......@@ -111,11 +111,12 @@ Vector &Vector::operator/=(Scalar b) {
data /= b;
return *this;
}
Vector &Vector::operator+=(const Vector &u) {
// data += u.data;
// return *this;
if(u.GetAccumulateFlag() && !(*this).GetAccumulateFlag()){
Warning("Adding accumulated and not accumulated Vector");
if (u.GetAccumulateFlag() && !(*this).GetAccumulateFlag()) {
Warning("Adding accumulated and not accumulated Vector");
}
if (u.size() == size()) {
data += u.data;
......@@ -136,7 +137,7 @@ Vector &Vector::operator+=(const Vector &u) {
for (int i = 0; i < r.n(); ++i)
(*this)(r, i) += u(rr, i);
}
return *this;
return *this;
}
void Vector::ConsistentAddition(const Vector &u) {
......@@ -156,18 +157,18 @@ void Vector::ConsistentAddition(const Vector &u) {
}
}
Scalar ConsistentScalarProduct(const Vector &u, const Vector &v){
Scalar ConsistentScalarProduct(const Vector &u, const Vector &v) {
double s = 0.0;
for (row r = u.rows(); r != u.rows_end(); ++r) {
procset p = u.find_procset(r());
if (p == u.procsets_end()) {
for (int i = 0; i < r.n(); ++i) {
s += u(r, i) * v(r,i);
s += u(r, i) * v(r, i);
}
} else {
if (p.master() == PPM->proc()) {
for (int i = 0; i < r.n(); ++i) {
s += u(r, i) * v(r,i);
s += u(r, i) * v(r, i);
}
}
}
......@@ -176,7 +177,7 @@ Scalar ConsistentScalarProduct(const Vector &u, const Vector &v){
}
Vector &Vector::operator-=(const Vector &u) {
if(u.GetAccumulateFlag() && !(*this).GetAccumulateFlag()){
if (u.GetAccumulateFlag() && !(*this).GetAccumulateFlag()) {
Warning("Subtracting accumulated and not accumulated Vector");
}
data -= u.data;
......@@ -190,7 +191,7 @@ Vector &Vector::operator=(const constAB<Scalar, Vector> &au) {
}
Vector &Vector::operator+=(const constAB<Scalar, Vector> &au) {
if(au.second().GetAccumulateFlag() && !(*this).GetAccumulateFlag()){
if (au.second().GetAccumulateFlag() && !(*this).GetAccumulateFlag()) {
Warning("Adding accumulated and not accumulated Vector");
}
data.MultiplyPlus(au.first(), au.second().data);
......@@ -198,7 +199,7 @@ Vector &Vector::operator+=(const constAB<Scalar, Vector> &au) {
}
Vector &Vector::operator-=(const constAB<Scalar, Vector> &au) {
if(au.second().GetAccumulateFlag() && !(*this).GetAccumulateFlag()){
if (au.second().GetAccumulateFlag() && !(*this).GetAccumulateFlag()) {
Warning("Subtracting accumulated and not accumulated Vector");
}
data.MultiplyMinus(au.first(), au.second().data);
......@@ -299,11 +300,11 @@ void Vector::DirichletConsistent() {
}
void Vector::Accumulate() {
if(accumulateFlag == true){
if (accumulateFlag == true) {
Warning("Accumulating already accumulated Vector!");
} else {
accumulateFlag = true;
}
accumulateFlag = true;
}
if (identify()) AccumulateIdentify();
if (parallel()) AccumulateParallel();
if (!identify()) return;
......@@ -584,79 +585,103 @@ constAB<int, Vectors> operator*(const int &a, const Vectors &v) {
Vectors::Vectors(int n, const IDiscretizationT<> &disc, int spaceLevel, int timeLevel)
: VectorMatrixBase(disc, spaceLevel, timeLevel), V(n) {
if (n < 1) return;
V[0] = new Vector(disc, spaceLevel, timeLevel);
V[0] = std::make_unique<Vector>(disc, spaceLevel, timeLevel);
for (int i = 1; i < n; ++i) {
V[i] = new Vector(*V[0]);
V[i] = std::make_unique<Vector>(*V[0]);
}
}
Vectors::Vectors(int n, Scalar b, const IDiscretizationT<> &disc, int spaceLevel, int timeLevel)
: VectorMatrixBase(disc, spaceLevel, timeLevel), V(n) {
if (n < 1) return;
V[0] = new Vector(b, disc, spaceLevel, timeLevel);
V[0] = std::make_unique<Vector>(b, disc, spaceLevel, timeLevel);
for (int i = 1; i < n; ++i) {
V[i] = new Vector(b, *V[0]);
V[i] = std::make_unique<Vector>(b, *V[0]);
}
}
Vectors::Vectors(int n, const Vector &u) : VectorMatrixBase(u), V(n) {
for (int i = 0; i < n; ++i) {
V[i] = new Vector(u);
V[i] = std::make_unique<Vector>(u);
}
}
Vectors::Vectors(int n, const Vectors &u) : VectorMatrixBase(u), V(n) {
for (int i = 0; i < std::min(n, u.size()); ++i) {
V[i] = new Vector(u[i]);
V[i] = std::make_unique<Vector>(u[i]);
}
for (int i = std::min(n, u.size()); i < n; ++i) {
V[i] = new Vector(0.0, u[0]);
V[i] = std::make_unique<Vector>(0.0, u[0]);
}
}
Vectors::Vectors(const Vectors &u) : VectorMatrixBase(u[0]), V(u.size()) {
for (int i = 0; i < V.size(); ++i) {
V[i] = new Vector(*(u.V[i]));
V[i] = std::make_unique<Vector>(*(u.V[i]));
*(V[i]) = *(u.V[i]);
}
}
Vectors::Vectors(const constAB<Operator, Vectors> &Ov) :
VectorMatrixBase(Ov.second()), V(Ov.second().size()) {
for (int i = 0; i < V.size(); ++i) {
V[i] = new Vector(Ov.second()[i]);
int computeLength(int length, int startIdx, int endIdx) {
if (endIdx < 0 || endIdx >= length) endIdx = length - 1;
return endIdx - startIdx + 1;
}
Vectors::Vectors(const Vectors &u, int startIdx, int endIdx)
: VectorMatrixBase(u), V(computeLength(u.size(), startIdx, endIdx)) {
for (int i = 0; i < V.size(); ++i, ++startIdx) {
V[i] = std::make_unique<Vector>(*u.V[startIdx]);
}
Ov.first().multiply(*this, Ov.second());
}
Vectors::~Vectors() {
Vectors::Vectors(const constAB<Operator, Vectors> &Ov) :
VectorMatrixBase(Ov.second()), V(Ov.second().size()) {
for (int i = 0; i < V.size(); ++i) {
delete V[i];
V[i] = std::make_unique<Vector>(Ov.second()[i]);
}
Ov.first().multiply(*this, Ov.second());
}
void Vectors::resize(int N) {
for (int i = 0; i < size(); ++i) delete V[i];
V.resize(N);
for (int i = 0; i < size(); ++i) V[i] = new Vector(disc, spaceLevel, timeLevel);
if (N == V.size()) {
for (int i = 0; i < V.size(); ++i)
*V[i] = 0.0;
} else if (N < V.size()) {
V.erase(V.begin() + N, V.end());
for (int i = 0; i < V.size(); ++i)
*V[i] = 0.0;
} else {
if (V.size() == 0)
V.push_back(std::make_unique<Vector>(0.0, disc, spaceLevel, timeLevel));
else {
for (int i = 0; i < V.size(); ++i)
*V[i] = 0.0;
}
for (int i = V.size(); i < N; ++i)
V.push_back(std::make_unique<Vector>(0.0, *V[0]));
}
}
void Vectors::resizeData(int N) {
int N_old = size();
for (int i = N; i < size(); ++i) delete V[i];
V.resize(N);
for (int i = N_old; i < size(); ++i) V[i] = new Vector(disc, spaceLevel, timeLevel);
if (N == V.size()) return;
if (N < V.size()) {
V.erase(V.begin() + N, V.end());
} else {
if (V.size() == 0)
V.push_back(std::make_unique<Vector>(0.0, disc, spaceLevel, timeLevel));
for (int i = V.size(); i < N; ++i)
V.push_back(std::make_unique<Vector>(0.0, *V[0]));
}
}
void Vectors::set(const Vectors &u) {
for (int i = 0; i < min(size(), u.size()); ++i)
*(V[i]) = u[i];
for (int i = min(size(), u.size()); i < size(); ++i)
*(V[i]) = 0;
*(V[i]) = 0.0;
}
void Vectors::erase(int i) {
mout << OUT("TODO") << endl;
V.erase(V.begin() + i);
}
......@@ -868,7 +893,7 @@ Loader &Vectors::load(Loader &loader) {
}
Vector operator+(const Vector &u, const Vector &v) {
if(u.GetAccumulateFlag() && !v.GetAccumulateFlag()){
if (u.GetAccumulateFlag() && !v.GetAccumulateFlag()) {
Warning("Adding accumulated and not accumulated Vector");
}
Vector w = u;
......@@ -876,7 +901,7 @@ Vector operator+(const Vector &u, const Vector &v) {
}
Vector operator-(const Vector &u, const Vector &v) {
if(u.GetAccumulateFlag() && !v.GetAccumulateFlag()){
if (u.GetAccumulateFlag() && !v.GetAccumulateFlag()) {
Warning("Subtracting accumulated and not accumulated Vector");
}
Vector w = u;
......@@ -885,7 +910,7 @@ Vector operator-(const Vector &u, const Vector &v) {
Scalar operator*(const Vector &u, const Vector &v) {
// Todo: Is there a case where both Vectors are on different communicators?
if(u.GetAccumulateFlag() && v.GetAccumulateFlag()){
if (u.GetAccumulateFlag() && v.GetAccumulateFlag()) {
Warning("Both vectors accumulated in scalar product!");
}
return PPM->SumOnCommSplit(u.data * v.data, u.CommSplit());
......
......@@ -250,15 +250,12 @@ double norm(double &u);
class Vectors : public VectorMatrixBase {
protected:
vector<Vector *> V;
vector<std::unique_ptr<Vector>> V;
public:
Vectors(int n, const IDiscretizationT<> &disc, int spaceLevel = -1, int timeLevel = -1);
Vectors(int n,
Scalar b,
const IDiscretizationT<> &disc,
int spaceLevel = -1,
int timeLevel = -1);
Vectors(int n, Scalar b, const IDiscretizationT<> &disc,
int spaceLevel = -1, int timeLevel = -1);
Vectors(int n, const Vector &u);
......@@ -266,9 +263,10 @@ public:
Vectors(const Vectors &u);
Vectors(const constAB<Operator, Vectors> &Ov);
/// selects part of u (indices included)
Vectors(const Vectors &u, int startIdx, int endIdx = -1);
~Vectors();
Vectors(const constAB<Operator, Vectors> &Ov);
int size() const { return V.size(); }
......
......@@ -2,7 +2,7 @@
Eigenpair::Eigenpair(const IDiscretizationT<> &disc, int spaceLevel, int timeLevel)
: VectorMatrixBase(disc, spaceLevel, timeLevel), U(disc, spaceLevel, timeLevel), lambda(0) {
: VectorMatrixBase(disc, spaceLevel, timeLevel), U(disc, spaceLevel, timeLevel), lambda(0) {
U = 0; //TODO: Can be removed if Vector(s) are initialized with 0 in constructors anyway
}
......@@ -34,7 +34,8 @@ Eigenpair &Eigenpair::operator-=(const double &shift) {
Eigenpairs::Eigenpairs(int N, const IDiscretizationT<> &disc, int spaceLevel, int timeLevel)
: VectorMatrixBase(disc, spaceLevel, timeLevel), U(N, disc, spaceLevel, timeLevel), lambda(N) {
for (int i = 0; i < N; ++i) //TODO: Can be removed if Vector(s) are initialized with 0 in constructors anyway
for (int i = 0;
i < N; ++i) //TODO: Can be removed if Vector(s) are initialized with 0 in constructors anyway
U[i] = 0;
}
......@@ -42,6 +43,10 @@ Eigenpairs::Eigenpairs(int N, const Eigenpairs &EP) : VectorMatrixBase(EP), U(N,
Eigenpairs::Eigenpairs(const Eigenpairs &EP) : VectorMatrixBase(EP), U(EP.U), lambda(EP.lambda) {}
Eigenpairs::Eigenpairs(const Eigenpairs &EP, int startIdx, int endIdx)
: VectorMatrixBase(EP), U(EP.U, startIdx, endIdx),
lambda(EP.lambda, startIdx, startIdx + U.size()) {}
void Eigenpairs::set(const Eigenpairs &EP) {
for (int i = 0; i < min(size(), EP.size()); ++i) {
U[i] = EP.U[i];
......
......@@ -33,12 +33,8 @@ public:
Eigenpair &operator-=(const double &shift);
Eigenfct &getEigenfct() { return U; }
const Eigenfct &getEigenfct() const { return U; }
Eigenvalue &getEigenvalue() { return lambda; }
const Eigenvalue &getEigenvalue() const { return lambda; }
friend inline Logging &operator<<(Logging &s, const Eigenpair &ep) {
......@@ -69,6 +65,9 @@ public:
Eigenpairs(const Eigenpairs &EP);
/// selects part of EP (indices included)
Eigenpairs(const Eigenpairs &EP, int startIdx, int endIdx = -1);
void set(const Eigenpairs &EP);
#ifdef BUILD_IA
......@@ -84,24 +83,14 @@ public:
Eigenpairs &operator-=(const double shift);
Eigenpair getEigenpair(int n) { return Eigenpair(U[n], lambda[n]); }
const Eigenpair getEigenpair(int n) const { return Eigenpair(U[n], lambda[n]); }
Eigenfcts &getEigenfcts() { return U; }
const Eigenfcts &getEigenfcts() const { return U; }
Eigenvalues &getEigenvalues() { return lambda; }
const Eigenvalues &getEigenvalues() const { return lambda; }
Eigenfct &operator()(int n) { return U[n]; }
const Eigenfct &operator()(int n) const { return U[n]; }
Eigenvalue &operator[](unsigned int n) { return lambda[n]; }
const Eigenvalue &operator[](unsigned int n) const { return lambda[n]; }
int size() const { return U.size(); }
......@@ -123,6 +112,12 @@ public:
friend inline Loader &operator>>(Loader &l, Eigenpairs &EP) {
return l >> EP.lambda >> EP.U;
}
template<bool VERIFIED>
friend class IAEigenvalueMethod;
template<bool VERIFIED>
friend class IAHomotopyMethod;
};
#endif // EIGENPAIR_H
......@@ -30,7 +30,7 @@ template<bool VERIFIED>
void IAEigenvalueMethod<VERIFIED>::computeEigenpairs(IAEigenvalueAssemble<VERIFIED> &assemble,
Eigenpairs &EP, Matrix &A, Matrix &B,
const double &t, bool init) {
computeEigenpairs(assemble, EP.getEigenfcts(), EP.getEigenvalues(), A, B, t, init);
computeEigenpairs(assemble, EP.U, EP.lambda, A, B, t, init);
}
......@@ -39,7 +39,7 @@ void IAEigenvalueMethod<VERIFIED>::computeEigenpairs(IAEigenvalueAssemble<VERIFI
Eigenpairs &EP, const double &t) {
Matrix A(EP(0));
Matrix B(EP(0));
computeEigenpairs(assemble, EP.getEigenfcts(), EP.getEigenvalues(), A, B, t, true);
computeEigenpairs(assemble, EP.U, EP.lambda, A, B, t, true);
}
template<bool VERIFIED>
......@@ -298,7 +298,8 @@ double IAHomotopyMethod<VERIFIED>::operator()(IAEigenvalueAssemble<VERIFIED> &as
outputStep(rho, EP.size());
if (nearZero >= 0.0 && EP.size() > 0 && EP[0] > nearZero) break;
if (EP.size() == 0) break;
if (nearZero >= 0.0 && EP[0] > nearZero) break;
}
checkHomotopy(assemble, EP, rho);
......@@ -420,12 +421,8 @@ void IAHomotopyMethod<VERIFIED>::computeNewRho(IAEigenvalueAssemble<VERIFIED> &a
if (clusterSize == n) break;
}
Eigenpairs GEP(clusterSize, EP);
Eigenpairs GEP(EP, n - clusterSize);
IALowerEigenvalueBounds lambda;
for (int i = 0; i < clusterSize; ++i) {
GEP(i) = EP(n - clusterSize + i);
GEP[i] = EP[n - clusterSize + i];
}
GM(assemble, levelGoerisch, GEP, lambda, rho, false, step);
rho = lambda[0];
......@@ -439,7 +436,7 @@ void IAHomotopyMethod<VERIFIED>::checkHomotopy(IAEigenvalueAssemble<VERIFIED> &a
const double &rho) {
if (EP.size() > 0) {
IAUpperEigenvalueBounds Lambda;
RM(assemble, EP.getEigenfcts(), Lambda, step);
RM(assemble, EP.U, Lambda, false, step);
if (Lambda[Lambda.size() - 1] > rho) ERROR("Homotopy failed")
}
mout << "exactly " << EP.size() << " eigenvalue"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment