Skip to content

GitLab

A
amdis
Commit e551c4a2 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files
Options

* TEST_EXIT (re)introduced 

* Small changes to make execution of surfacediffusion code possible
parent 42521b59
Hide whitespace changes
Inline Side-by-side
Showing with 238 additions and 456 deletions
AMDiS/src/DiagonalPreconditioner.cc
View file @ e551c4a2
......@@ -13,7 +13,7 @@ namespace AMDiS {
TEST_EXIT_DBG(matrix[row])("no matrix\n");
TEST_EXIT_DBG(x)("no solution vector\n");
DOFMatrix::Iterator rowIterator((*matrix[row]), USED_DOFS);
DOFVector<double>::Iterator vecIterator(x, USED_DOFS);
......
AMDiS/src/Flag.h
View file @ e551c4a2
......@@ -61,57 +61,78 @@ namespace AMDiS {
/** \brief
* Compares two Flags
*/
inline bool operator==(const Flag& f) const {return (flags==f.flags);};
inline bool operator==(const Flag& f) const {
return (flags == f.flags);
};
/** \brief
* Compares two Flags
*/
inline bool operator!=(const Flag& f) const {return !(f==*this);};
inline bool operator!=(const Flag& f) const {
return !(f == *this);
};
/** \brief
* Assignment operator
*/
inline Flag& operator=(const Flag& f) {
if (this!=&f) flags=f.flags;
if (this != &f)
flags = f.flags;
return *this;
};
/** \brief
* Typecast
*/
inline operator bool() const { return isAnySet(); };
inline operator bool() const {
return isAnySet();
};
/** \brief
* Set \ref flags
*/
inline void setFlags(const unsigned long f) { flags = f; };
inline void setFlags(const unsigned long f) {
flags = f;
};
/** \brief
* Set \ref flags
*/
inline void setFlags(const Flag& f) { flags = f.flags; };
inline void setFlags(const Flag& f) {
flags = f.flags;
};
/** \brief
* Sets \ref flags to \ref flags | f
*/
inline void setFlag(const unsigned long f) { flags |= f; };
inline void setFlag(const unsigned long f) {
flags |= f;
};
/** \brief
* Sets \ref flags to \ref flags | f.flags
*/
inline void setFlag(const Flag& f) { flags |= f.flags; };
inline void setFlag(const Flag& f) {
flags |= f.flags;
};
/** \brief
* Sets \ref flags to \ref flags & ~f
*/
inline void unsetFlag(const unsigned long f) { flags &= ~f; };
inline void unsetFlag(const unsigned long f) {
flags &= ~f;
};
/** \brief
* Sets \ref flags to \ref flags & ~f.flags
*/
inline void unsetFlag(const Flag& f) { flags &= ~f.flags; };
inline void unsetFlag(const Flag& f) {
flags &= ~f.flags;
};
inline const unsigned long getFlags() const { return flags; };
inline const unsigned long getFlags() const {
return flags;
};
/** \brief
* Returns \ref flags | f.flags
......@@ -145,7 +166,7 @@ namespace AMDiS {
*/
inline Flag operator&(const Flag& f) const {
Flag r(flags);
r.flags &=f.flags;
r.flags &= f.flags;
return r;
};
......@@ -170,7 +191,9 @@ namespace AMDiS {
* Sets \ref flags to \ref flags & f.flags
*/
inline Flag& operator|=(const Flag& f) {
if (this!=&f) {flags |=f.flags;};
if (this != &f) {
flags |= f.flags;
};
return *this;
};
......@@ -193,12 +216,16 @@ namespace AMDiS {
/** \brief
* Returns !\ref isSet(f)
*/
inline bool isUnset(const Flag& f) const { return !isSet(f); };
inline bool isUnset(const Flag& f) const {
return !isSet(f);
};
/** \brief
* Returns true if \ref flags != 0
*/
inline bool isAnySet() const { return (flags != 0); };
inline bool isAnySet() const {
return (flags != 0);
};
protected:
/** \brief
......
AMDiS/src/GMResSolver.h
View file @ e551c4a2
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// ============================================================================
......@@ -22,21 +21,32 @@
#ifndef AMDIS_GMRESSOLVER_H
#define AMDIS_GMRESSOLVER_H
#include <vector>
#include "OEMSolver.h"
#include "MemoryManager.h"
namespace AMDiS {
// ============================================================================
// ===== class GMResSolver ====================================================
// ===== class GMResSolver ===================================================
// ============================================================================
/**
* \ingroup Solver
*
* \brief
* Solves a linear system by the GMRes method with restart and can be used for
* regular system matrices.
* Solves a linear system by the restarted GMRES(m) method (Generalized Minimal
* Residual) for general nonsymmetric system matrices.
*
* The implementation is based on the following book: "Numerik linearer
* Gleichungssystene", 2. Auflage, Andreas Meister. The algorithm is described
* on pages 150 and 154 (for the restarted version). The extension to the
* preconditioned GMRES(m) is based on the description in the book "Iterative
* Methods for Sparse Linear Systems", second edition, Yousef Saad, on page 268.
*
* The orthogonal system is build using the modified Gram Schmidt (MGS) method.
* The least square problem is solved using Givens transformations.
*
*/
template<typename VectorType>
class GMResSolver : public OEMSolver<VectorType>
......@@ -88,41 +98,61 @@ namespace AMDiS {
*/
void exit();
/** \brief
* One run of the GMRES(m) method. Is called from solveSystem().
*/
int gmres(MatVecMultiplier<VectorType> *mv,
VectorType *x,
VectorType *b);
private:
/** \brief
* internal used method
* Stores the tolerance boundary for numerical computations.
*/
double householderVector(double sigma,
VectorType *u,
VectorType *x,
int offset);
double TOL_;
/** \brief
* internal used method
* The parameter m of GMRES(m), i.e. the number of iterations before a restart of
* the algorithm is forced.
*/
void newBasisVector(int m, int k, VectorType **U, double **LR,
MatVecMultiplier<VectorType> *matVec,
VectorType *v, VectorType *b, double* y);
int restart_;
/** \brief
* internal used method
* Stores intermediate results in the orthogonalization step.
*/
int solve_k(MatVecMultiplier<VectorType> *mv,
VectorType *x, const VectorType *b);
VectorType *w_;
private:
double TOL;
int restart;
int dim;
VectorType *z_;
VectorType *r_;
/** \brief
* Pointers to the vectors of the orthogonal system.
*/
VectorType **v;
/** \brief
* Stores the gamma values for Givens transformations.
*/
::std::vector<double> gamma;
/** \brief
* Stores the c values for Givens transformations.
*/
::std::vector<double> c;
// pointer to memory needed for solveSystem
VectorType *r, *v;
VectorType **U;
double **LR;
double **givens, *w, *y;
/** \brief
* Stores the s values for Givens transformations.
*/
::std::vector<double> s;
double gmres_k_residual_0;
::std::vector<double> y_;
/** \brief
* H-Matrix computed in the GMRES algorithm.
*/
::std::vector< ::std::vector<double> > h;
};
......
AMDiS/src/GMResSolver.hh
View file @ e551c4a2
#include "GMResSolver.h"
#include "Preconditioner.h"
#include <algorithm>
namespace AMDiS {
template<typename VectorType>
GMResSolver<VectorType>::GMResSolver(::std::string name)
: OEMSolver<VectorType>(name),
TOL(1.e-25),
restart(0),
dim(0),
r(NULL),
v(NULL),
U(NULL),
LR(NULL),
givens(NULL),
w(NULL),
y(NULL),
gmres_k_residual_0(0.0)
TOL_(1.e-25),
restart_(0),
w_(NULL)
{
FUNCNAME("GMResSolver::GMResSolver");
GET_PARAMETER(0, name + "->restart", "%d", &restart);
TEST_EXIT(restart)("restart not set\n");
FUNCNAME("GMResSolver::GMResSolver()");
GET_PARAMETER(0, name + "->restart", "%d", &restart_);
TEST_EXIT(restart_)("restart not set\n");
}
template<typename VectorType>
......@@ -30,305 +22,162 @@ namespace AMDiS {
template<typename VectorType>
void GMResSolver<VectorType>::init()
{
r = this->vectorCreator->create();
dim = size(r);
restart = ::std::max(0, ::std::min(restart, dim));
// Create a new vector w.
w_ = this->vectorCreator->create();
z_ = this->vectorCreator->create();
r_ = this->vectorCreator->create();
// Get memory for the vector pointers and create the pointers afterwards.
v = GET_MEMORY(VectorType*, restart_ + 1);
for (int i = 0; i <= restart_; i++) {
v[i] = this->vectorCreator->create();
}
// Resize all fields to the corresponding size.
gamma.resize(restart_ + 1);
c.resize(restart_);
s.resize(restart_);
v = this->vectorCreator->create();
U = GET_MEMORY(VectorType*, restart);
LR = GET_MEMORY(double*, restart);
givens = GET_MEMORY(double*, restart);
y_.resize(restart_);
for (int i = 0; i < restart; i++) {
U[i] = this->vectorCreator->create();
LR[i] = GET_MEMORY(double, restart);
givens[i] = GET_MEMORY(double, 2);
h.resize(restart_ + 1);
for (int i = 0; i <= restart_; i++) {
h[i].resize(restart_);
}
w = GET_MEMORY(double, restart);
y = GET_MEMORY(double, restart);
}
template<typename VectorType>
void GMResSolver<VectorType>::exit()
{
if (r) {
int dim = restart;
int k = ::std::max(0,::std::min(restart, dim));
this->vectorCreator->free(r);
this->vectorCreator->free(v);
for (int i = 0; i < k; i++) {
this->vectorCreator->free(U[i]);
FREE_MEMORY(LR[i], double, k);
FREE_MEMORY(givens[i], double, 2);
// Empty used memory.
if (w_) {
this->vectorCreator->free(w_);
this->vectorCreator->free(z_);
this->vectorCreator->free(r_);
for (int i = 0; i <= restart_; i++) {
this->vectorCreator->free(v[i]);
}
FREE_MEMORY(U, VectorType*, k);
FREE_MEMORY(LR, double*, k);
FREE_MEMORY(givens, double*, k);
FREE_MEMORY(w, double, k);
FREE_MEMORY(y, double, k);
}
}
template<typename VectorType>
int GMResSolver<VectorType>::solve_k(MatVecMultiplier<VectorType> *matVec,
VectorType *x, const VectorType *b)
int GMResSolver<VectorType>::gmres(MatVecMultiplier<VectorType> *matVec,
VectorType *x,
VectorType *b)
{
FUNCNAME("GMResSolver::solve_k");
VectorType *um1 = NULL;
double c, s, wm1;
int k = ::std::max(0, ::std::min(restart, dim));
FUNCNAME("GMResSolver::gmres()");
/*--------------------------------------------------------------------------*/
/* Initialization */
/*--------------------------------------------------------------------------*/
// r_0 = b - Ax, where r_0 is already stored as the first vector in the
// matrix V.
matVec->matVec(NoTranspose, *x, *v[0]);
xpay(-1.0, *b, *v[0]);
// r = Ax
matVec->matVec(NoTranspose, *x, *r);
for (int i = 0; i < dim; i++) {
(*r)[i] -= (*b)[i];
(*r)[i] *= -1.0;
// Left preconditioning, if required.
if (this->leftPrecon) {
this->leftPrecon->precon(v[0]);
}
if (this->leftPrecon)
this->leftPrecon->precon(r);
gmres_k_residual_0 = norm(r);
if (gmres_k_residual_0 < this->tolerance) {
this->residual = gmres_k_residual_0;
return(0);
}
/*--------------------------------------------------------------------------*/
/* construct k-dimensional Krylov space */
/*--------------------------------------------------------------------------*/
wm1 = householderVector(gmres_k_residual_0, U[0], r, 0);
um1 = U[0];
double rj, sigma, maxi, val;
VectorType *uj;
int m = 0;
for (m = 0; m < k; m++) {
w[m] = wm1;
newBasisVector(m + 1, k, U, LR, matVec, r, v, y);
if (m + 1 < dim) {
double norm = 0.0;
for (int i = m + 1; i < dim; i++) {
norm += (*r)[i] * (*r)[i];
}
norm = sqrt(norm);
if (norm > TOL) {
/****************************************************************************/
/* one of the last components of r is not zero; calculate Householder */
/* vector; if m < k-1, we need the Householder vector for the calculation */
/* of the next basis vector => store it; if m == k-1 we do not need this */
/* vector anymore => do not store it! */
/****************************************************************************/
if (m < k - 1) {
um1 = U[m + 1];
(*r)[m + 1] = householderVector(norm, um1, r, m + 1);
} else {
(*r)[m + 1] = householderVector(norm, NULL, r, m + 1);
}
}
}
for (int j = 0; j < m; j++) {
rj = (*r)[j];
c = givens[j][0];
s = givens[j][1];
(*r)[j] = c * rj + s * (*r)[j + 1];
(*r)[j+1] = -s * rj + c * (*r)[j + 1];
}
if (m + 1 < dim && abs((*r)[m + 1]) > TOL) {
/****************************************************************************/
/* Find Givens rotation J_m such that, */
/* a) (J_m r)[i] = r[i], i < m, */
/* b) (J_m r)[m+1] = 0.0 */
/* => (J_m r)[m] = +- sqrt(r[m]^2 + r[m+1]^2) =: sigma */
/* */
/* */
/* |1 0 . . . . 0| */
/* |0 1 . . . . .| */
/* |. . .| c = r[m]/sigma */
/* J_m = |. . .| s = r[m+1]/sigma */
/* |. . .| */
/* |. c s| m */
/* |0 . . . . -s c| m+1 */
/* m m+1 */
/****************************************************************************/
maxi = ::std::max((*r)[m], (*r)[m+1]);
// Compute the norm of r_0 = v_0.
gamma[0] = norm(v[0]);
c = (*r)[m] / maxi;
s = (*r)[m + 1] / maxi;
sigma = maxi * sqrt(c * c + s * s);
if ((*r)[m] < 0)
sigma = -sigma;
givens[m][0] = c = (*r)[m]/sigma;
givens[m][1] = s = (*r)[m+1]/sigma;
(*r)[m] = sigma; /* r[m+1] == 0 automatically! */
wm1 = -s * w[m]; /* |wm1| is the new residual! */
w[m] *= c;
} else {
wm1 = 0.0;
}
// Normalize v_0.
*v[0] *= (1 / gamma[0]);
/****************************************************************************/
/* now, store the first m components of the column vector r in the the mth */
/* column of LR */
/****************************************************************************/
// If the norm of gamma_0 is less than the tolerance bounday, x is already
// a good solution for Ax = b.
if (gamma[0] < this->tolerance) {
this->residual = gamma[0];
return(0);
}
// Main loop of the GMRES algorithm.
for (int j = 0; j < restart_; j++) {
// w_j = A * v_j;
matVec->matVec(NoTranspose, *v[j], *w_);
for (int j = 0; j <= m; j++)
LR[j][m] = (*r)[j];
// Preconditioning of w_j
if (this->leftPrecon) {
this->leftPrecon->precon(w_);
}
/****************************************************************************/
/* test, whether tolarance is reached or not */
/****************************************************************************/
// w_j = A * v_j - \sum(h_ij * v_i)
for (int i = 0; i <= j; i++) {
// h_ij = (v_i, A * v_j)_2
h[i][j] = *w_ * (*v[i]);
if (abs(wm1) < this->tolerance) {
m++;
break;
// Update w_j
axpy(-h[i][j], *v[i], *w_);
}
/****************************************************************************/
/* tolerance not reached: calculate and store (m+1)th row of matrix L; */
/* this row is only needed for the computation of the (m+1)th column of */
/* the orthogonal matrix Q_m; this vector is the additional basis vector */
/* for the enlargement of the Krylov space; only needed for m < k-1 */
/* L[m+1][j] = <u_(m+1),u_j>_2 */
/* (m+1)th vector u = umi is allready stored at U+(m+1) */
/****************************************************************************/
if (m < k-1) {
for (int j = 0; j < m + 1; j++) {
uj = U[j];
val = 0.0;
for (int i = m + 1; i < dim; i++) {
val += (*um1)[i] * (*uj)[i];
}
LR[(m+1)][j] = 2.0 * val;
}
// h_{j + 1}{j} = \| w_j \|_2
h[j + 1][j] = norm(w_);
// Compute h_ij and h_{i + 1}{j} using parameters of the Givens rotation.
for (int i = 0; i < j; i++) {
double t1 = c[i] * h[i][j] + s[i] * h[i + 1][j];
double t2 = -s[i] * h[i][j] + c[i] * h[i + 1][j];
h[i][j] = t1;
h[i + 1][j] = t2;
}
}
/****************************************************************************/
/* and now solve the upper triangular system */
/****************************************************************************/
y[m - 1] = w[m - 1] / LR[m - 1][m - 1]; /* = LR[(m-1)*k+(m-1)]! */
for (int j = m - 2; j >= 0; j--) {
double yj = w[j];
for (int l = j + 1; l < m; l++) {
yj -= LR[j][l] * y[l];
}
// Compute new beta
double beta = sqrt(h[j][j] * h[j][j] + h[j + 1][j] * h[j + 1][j]);
// Update s_j and c_j
s[j] = h[j + 1][j] / beta;
c[j] = h[j][j] / beta;
// h_jj = beta
h[j][j] = beta;
// Update gammas
gamma[j + 1] = -s[j] * gamma[j];
gamma[j] *= c[j];
y[j] = yj / LR[j][j];
// Set v_{j + 1} to w_j and normalize it.
*v[j + 1] = *w_;
*v[j + 1] *= 1 / h[j + 1][j];
}
/****************************************************************************/
/* calculate v = 2 U_m^T [e_0,....,e_m-1] y */
/****************************************************************************/
for (int i = 0; i < m; i++) {
val = 0.0;
for (int j = i; j < m; j++) {
val += (*(U[i]))[j] * y[j];
}
(*v)[i] = 2.0 * val;
}
/****************************************************************************/
/* now solve L_m^T w = v in R^m (upper triagonal system, diagonal == 1.0 */
/****************************************************************************/
w[m - 1] = (*v)[m - 1];
for (int j = m - 2; j >= 0; j--) {
double wj = (*v)[j];
for (int l = j + 1; l < m; l++) {
wj -= LR[l][j]*w[l];
// Compute the solution.
for (int i = restart_ - 1; i >= 0; i--) {
for (int k = i + 1; k < restart_; k++) {
gamma[i] -= h[i][k] * gamma[k];
}
w[j] = wj;
gamma[i] /= h[i][i];
axpy(gamma[i], *v[i], *x);
}