diff --git a/AMDiS/src/DiagonalPreconditioner.cc b/AMDiS/src/DiagonalPreconditioner.cc
index 1678fabc84f9494515812afdea1f94e446eda05a..48a9c81f9129233d9e150a4e44d52d2e06711885 100644
--- a/AMDiS/src/DiagonalPreconditioner.cc
+++ b/AMDiS/src/DiagonalPreconditioner.cc
@@ -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);
diff --git a/AMDiS/src/Flag.h b/AMDiS/src/Flag.h
index 8d0b6a255350e1ec6e04ca79ff1709b9895ea88f..404643818eed061897a062b2bf166940b1512c6d 100644
--- a/AMDiS/src/Flag.h
+++ b/AMDiS/src/Flag.h
@@ -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
diff --git a/AMDiS/src/GMResSolver.h b/AMDiS/src/GMResSolver.h
index b4412d5feb9b87923edec6071ec08401f31ea597..bb539a6966d8d0c277a30d003bfddd7011cfdbad 100644
--- a/AMDiS/src/GMResSolver.h
+++ b/AMDiS/src/GMResSolver.h
@@ -1,5 +1,4 @@
// ============================================================================
-// == ==
// == 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;
};
diff --git a/AMDiS/src/GMResSolver.hh b/AMDiS/src/GMResSolver.hh
index fec30262d0a605ef9edc8a0fb8b5dbecee428221..9315d96ef46172298beabbba69b44dd80c7541f7 100644
--- a/AMDiS/src/GMResSolver.hh
+++ b/AMDiS/src/GMResSolver.hh
@@ -1,26 +1,18 @@
+#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);
}
- /****************************************************************************/
- /* v = [e_0,....,e_m-1] y - U_m w */
- /****************************************************************************/
-
- for (int i = 0; i < dim; i++) {
- double vi = 0.0;
- for (int j = 0; j < ::std::min(i + 1,m); j++)
- vi += (*(U[j]))[i] * w[j];
- (*v)[i] = -vi;
- }
-
- for (int j = 0; j < m; j++)
- (*v)[j] += y[j];
-
- /****************************************************************************/
- /* and now, make the update of u :-) */
- /****************************************************************************/
-
- if (this->rightPrecon)
- this->rightPrecon->precon(v);
-
- for (int i = 0; i < dim; i++)
- (*x)[i] += (*v)[i];
-
- this->residual = abs(wm1);
- return(m);
+ // Update the new residual.
+ this->residual = abs(gamma[restart_]);
+
+ return restart_;
}
+
template<typename VectorType>
int GMResSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x,
+ VectorType *b)
{
- FUNCNAME("GMResSolver::solveSystem()");
-
- double old_res = -1.0;
-
-// ::std::cout << "TEST!" << ::std::endl;
-
-// matVec->matVec(NoTranspose, *x, *r);
-// *r -= *b;
-// r->print();
+ FUNCNAME("GMResSolver::solveSystem()");
+ double old_res = -1.0;
-
- if (norm(b) < TOL) {
+ // If norm of b is smaller than the tolarance, we can assume b to be zero.
+ // Hence, x = 0 is the solution of the linear system.
+ if (norm(b) < TOL_) {
INFO(this->info, 2)("b == 0, x = 0 is the solution of the linear system\n");
setValue(*x, 0.0);
this->residual = 0.0;
return(0);
}
-
+
START_INFO();
for (int iter = 0; iter <= this->max_iter; iter++) {
- int k = solve_k(matVec, x, b);
+ // Solve Ax=b using GMRES(m).
+ int k = gmres(matVec, x, b);
+ // Check the new residual.
if (SOLVE_INFO(iter, this->residual, &old_res))
return(iter);
TEST_EXIT(k != 0)("this must not happen\n");
@@ -338,136 +187,4 @@ namespace AMDiS {
return(0);
}
- template<typename VectorType>
- double GMResSolver<VectorType>::householderVector(double sigma,
- VectorType *u,
- VectorType *x,
- int offset)
- {
- FUNCNAME("GMResSolver::householderVector()");
-
- if (dim <= 0) {
- MSG("dim <= 0\n");
- return(0.0);
- }
-
- if ((*x)[offset] >= 0) {
- if (u) {
- double beta = sqrt(0.5 / (sigma * (sigma + (*x)[offset])));
- (*u)[offset] = beta * ((*x)[offset] + sigma);
-
- for (int i = offset + 1; i < dim; i++)
- (*u)[i] = beta * (*x)[i];
- }
- return(-sigma);
-
- } else {
- if (u) {
- double beta = sqrt(0.5 / (sigma * (sigma - (*x)[offset])));
- (*u)[offset] = ((*x)[offset] - sigma) * beta;
-
- for (int i = offset + 1; i < dim; i++)
- (*u)[i] = (*x)[i] * beta;
- }
-
- return(sigma);
- }
- }
-
- template<typename VectorType>
- void GMResSolver<VectorType>::newBasisVector(int m, int k, VectorType **U,
- double **LR,
- MatVecMultiplier<VectorType> *matVec,
- VectorType *v,
- VectorType *b, double* y)
- {
- /****************************************************************************/
- /* first calculate y = 2 U_m^T e_m */
- /****************************************************************************/
-
- for (int j = 0; j < m; j++)
- (*b)[j] = 2.0*(*(U[j]))[m-1];
-
-
- /****************************************************************************/
- /* now solve L_m^T y = b in R^m (upper triagonal system, diagonal == 1.0 */
- /****************************************************************************/
-
- y[m-1] = (*b)[m-1];
- for (int j = m-2; j >= 0; j--) {
- double yj = (*b)[j];
-
- for (int l = j + 1; l < m; l++)
- yj -= LR[l][j]*y[l];
-
- y[j] = yj;
- }
-
- /****************************************************************************/
- /* b = -U_m y + e_m */
- /****************************************************************************/
-
- for (int i = 0; i < dim; i++) {
- double bi = 0.0;
-
- for (int j = 0; j < ::std::min(i + 1,m); j++)
- bi += (*(U[j]))[i] * y[j];
-
- (*b)[i] = -bi;
- }
- (*b)[m - 1] += 1.0;
-
- /****************************************************************************/
- /* v = Ab */
- /****************************************************************************/
-
- if (this->rightPrecon)
- this->rightPrecon->precon(b);
-
- matVec->matVec(NoTranspose, *b, *v);
-
- if (this->leftPrecon)
- this->leftPrecon->precon(v);
-
- /****************************************************************************/
- /* b = 2 U_m^T v in R^m */
- /****************************************************************************/
-
- for (int j = 0; j < m; j++) {
- double bj = 0.0;
-
- for (int i = j; i < dim; i++)
- bj += (*(U[j]))[i]*(*v)[i];
- (*b)[j] = 2.0*bj;
- }
-
- /****************************************************************************/
- /* now solve L_m y = b in R^m (lower triagonal system, diagonal == 1.0 */
- /****************************************************************************/
-
- y[0] = (*b)[0];
- for (int j = 1; j < m; j++) {
- double yj = (*b)[j];
-
- for (int l = 0; l < j; l++)
- yj -= LR[j][l] * y[l];
-
- y[j] = yj;
- }
-
- /****************************************************************************/
- /* v = v - U_m y */
- /****************************************************************************/
-
- for (int i = 0; i < dim; i++) {
- double vi = 0.0;
-
- for (int j = 0; j < ::std::min(i + 1,m); j++)
- vi += (*(U[j]))[i]*y[j];
-
- (*v)[i] -= vi;
- }
-
- return;
- }
}
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index 8d5358519b3f251b2996fbcf3e92fe400ad67230..1c3f7b3fbdc47d4933cba5d52efedd205e4a5394 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -346,7 +346,7 @@ namespace AMDiS {
/** \brief
* if test is false, an error message is printed and the program exits
*/
-#define TEST_EXIT(test)
+#define TEST_EXIT(test) if ((test));else ERROR_EXIT
/** \brief
* In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
diff --git a/AMDiS/src/ProblemInstat.cc b/AMDiS/src/ProblemInstat.cc
index 7a8d19569b306cc2547907302becf1abb3cfef9e..eb3d387d2d4bcc51c6f9a14d50517d53c074bc52 100644
--- a/AMDiS/src/ProblemInstat.cc
+++ b/AMDiS/src/ProblemInstat.cc
@@ -136,14 +136,14 @@ namespace AMDiS {
ProblemInstat::initialize(initFlag,adoptProblem,adoptFlag);
// === create vector for old solution ===
- if(oldSolution) {
+ if (oldSolution) {
WARNING("oldSolution already created\n");
} else {
- if(initFlag.isSet(INIT_UH_OLD)) {
+ if (initFlag.isSet(INIT_UH_OLD)) {
createUhOld();
}
- if(adoptProblem && adoptFlag.isSet(INIT_UH_OLD)) {
- ProblemInstatVec* _adoptProblem=dynamic_cast<ProblemInstatVec*>(adoptProblem);
+ if (adoptProblem && adoptFlag.isSet(INIT_UH_OLD)) {
+ ProblemInstatVec* _adoptProblem = dynamic_cast<ProblemInstatVec*>(adoptProblem);
TEST_EXIT(_adoptProblem)("can't adopt oldSolution from problem which is not instationary and vectorial");
TEST_EXIT(!oldSolution)("oldSolution already created");
oldSolution = _adoptProblem->getOldSolution();
@@ -155,15 +155,15 @@ namespace AMDiS {
}
void ProblemInstatVec::createUhOld() {
- if (oldSolution) WARNING("oldSolution already created\n");
- else {
+ if (oldSolution) {
+ WARNING("oldSolution already created\n");
+ } else {
int size = problemStat->getNumComponents();
// create oldSolution
oldSolution = NEW SystemVector("old solution",
problemStat->getFESpaces(),
size);
- int i;
- for(i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++) {
oldSolution->setDOFVector(i, NEW DOFVector<double>(problemStat->getFESpace(i),
name + "->uOld"));
oldSolution->getDOFVector(i)->refineInterpol(true);
diff --git a/AMDiS/src/ProblemInstat.h b/AMDiS/src/ProblemInstat.h
index 9d153b323fbbedbaff052181da334344742d9c86..5b3824a6c1ca5969897dce2d003c1e7c50577efc 100644
--- a/AMDiS/src/ProblemInstat.h
+++ b/AMDiS/src/ProblemInstat.h
@@ -88,15 +88,21 @@ namespace AMDiS {
/** \brief
*/
- virtual Flag markElements(AdaptInfo *adaptInfo) { return 0; };
+ virtual Flag markElements(AdaptInfo *adaptInfo) {
+ return 0;
+ };
/** \brief
*/
- virtual Flag refineMesh(AdaptInfo *adaptInfo) { return 0; };
+ virtual Flag refineMesh(AdaptInfo *adaptInfo) {
+ return 0;
+ };
/** \brief
*/
- virtual Flag coarsenMesh(AdaptInfo *adaptInfo) { return 0; };
+ virtual Flag coarsenMesh(AdaptInfo *adaptInfo) {
+ return 0;
+ };
/** \brief
* Implementation of ProblemTimeInterface::closeTimestep().
@@ -106,7 +112,9 @@ namespace AMDiS {
/** \brief
* Returns \ref name.
*/
- inline const ::std::string& getName() { return name; };
+ inline const ::std::string& getName() {
+ return name;
+ };
/** \brief
* Used by \ref problemInitial
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index bbd98a4fc535f9fb11ded84f47f19f9ed0ca5496..eb7438c8cd6e68a1c8720e255bc812410b72947f 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -775,7 +775,7 @@ namespace AMDiS {
FUNCNAME("ProblemVec::writeFiles()");
::std::list<FileWriterInterface*>::iterator it;
- for(it = fileWriters_.begin(); it != fileWriters_.end(); ++it) {
+ for (it = fileWriters_.begin(); it != fileWriters_.end(); ++it) {
(*it)->writeFiles(adaptInfo, force);
}
}