diff --git a/AMDiS/src/BFGS_Precond.cc b/AMDiS/src/BFGS_Precond.cc
new file mode 100755
index 0000000000000000000000000000000000000000..0a31bb1d68b93e50f7145c1c758c9d0b75fadde7
--- /dev/null
+++ b/AMDiS/src/BFGS_Precond.cc
@@ -0,0 +1,68 @@
+#include "BFGS_Precond.h"
+#include "QN_Precond.h"
+#include "FiniteElemSpace.h"
+#include "DOFAdmin.h"
+#include "DOFVector.h"
+#include "DOFMatrix.h"
+#include "DOFIterator.h"
+
+namespace AMDiS
+{
+
+ void BFGS_Precond::addPair(DOFVector<double> *s, DOFVector<double> *y,
+ double sy)
+ {
+ FUNCNAME("BFGS_Preconditioner::addPair()");
+
+ if (!sy)
+ sy = *s * *y;
+ TEST_EXIT(sy > 0)("Vectors do not satisfy the curvature condition!\n");
+
+ BFGS_Vectors newPair(s->getFESpace());
+ newPair.setVectors(s, y, sy);
+
+ if (k < m) {
+ k++;
+ } else {
+ pairsVec.pop_back();
+ }
+
+ pairsVec.push_front(newPair);
+
+ return;
+ }
+
+ void BFGS_Precond::precon(DOFVector<double> *vec)
+ {
+ FUNCNAME("BFGS_Preconditioner::precon()");
+
+ int i = 0;
+ double beta;
+ double *alpha = GET_MEMORY(double, k);
+
+ ::std::deque<BFGS_Vectors>::iterator pairsIterator(pairsVec.begin());
+
+ while (pairsIterator != pairsVec.end()) {
+ alpha[i] = pairsIterator->s * *vec;
+ alpha[i] /= pairsIterator->sy;
+ aXpY(-alpha[i], pairsIterator->y, *vec, bound);
+ ++i;
+ ++pairsIterator;
+ }
+
+ diagPrecond(*matrix[row], vec, bound);
+
+ while (pairsIterator != pairsVec.begin()) {
+ --pairsIterator;
+ --i;
+ beta = pairsIterator->y * *vec;
+ beta /= pairsIterator->sy;
+ aXpY(alpha[i]-beta, pairsIterator->s, *vec, bound);
+ }
+
+ FREE_MEMORY(alpha, double, k);
+
+ return;
+ }
+
+}
diff --git a/AMDiS/src/BFGS_Precond.h b/AMDiS/src/BFGS_Precond.h
new file mode 100755
index 0000000000000000000000000000000000000000..5629330327dd060a81f49b6907e876b7f6a62cff
--- /dev/null
+++ b/AMDiS/src/BFGS_Precond.h
@@ -0,0 +1,168 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == crystal growth group ==
+// == ==
+// == Stiftung caesar ==
+// == Ludwig-Erhard-Allee 2 ==
+// == 53175 Bonn ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == http://www.caesar.de/cg/AMDiS ==
+// == ==
+// ============================================================================
+
+/** \file BFGS_Precond.h */
+
+#ifndef AMDIS_BFGS_PRECONDITIONER_H
+#define AMDIS_BFGS_PRECONDITIONER_H
+
+#include "Parameters.h"
+#include "Preconditioner.h"
+#include "QN_Precond.h"
+
+namespace AMDiS
+{
+
+ template<typename T> class DOFVector;
+ template<typename T> class OEMSolver;
+
+ // ============================================================================
+ // ===== class BFGS_Vectors ===================================================
+ // ============================================================================
+
+ // Auxiliar class por storing pairs of DOFVectors and their scalar products.
+ class BFGS_Vectors
+ {
+ public:
+ MEMORY_MANAGED(BFGS_Vectors);
+
+ /** \brief
+ * Constructor.
+ */
+ BFGS_Vectors(const FiniteElemSpace *fe_space)
+ : s(fe_space, "BFGS->s"), y(fe_space, "BFGS->y")
+ {};
+
+ inline void setVectors(DOFVector<double> *s_, DOFVector<double> *y_,
+ double sy_=0.0)
+ {
+ s = *s_;
+ y = *y_;
+
+ if (!sy_)
+ sy = s * y;
+ else
+ sy = sy_;
+
+ return;
+ };
+
+ private:
+ DOFVector<double> s, y;
+ double sy;
+
+ friend class BFGS_Precond;
+ };
+
+ // ============================================================================
+ // ===== class BFGS_Precond ===================================================
+ // ============================================================================
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * BFGS preconditioner.
+ */
+ class BFGS_Precond : public PreconditionerScal
+ {
+ public:
+ MEMORY_MANAGED(BFGS_Precond);
+
+ /** \brief
+ * Creator class used in the PreconditionerMap.
+ */
+ class Creator : public PreconditionerScalCreator
+ {
+ public:
+ MEMORY_MANAGED(Creator);
+
+ /** \brief
+ * Creates a new BFGS Preconditioner.
+ */
+ PreconditionerScal *create()
+ {
+ return NEW BFGS_Precond(size, row);
+ };
+ };
+
+ /** \brief
+ * Constructor.
+ */
+ BFGS_Precond(int size_ = 1, int row_ = 0) : PreconditionerScal(size_, row_) {
+ k = 0;
+
+ m = 2 * Global::getGeo(WORLD);
+ GET_PARAMETER(1, "preconditioner->max. number of pairs", "%d", &m);
+ m = std::max(0, std::min(m, 32));
+ };
+
+ /** \brief
+ * Destructor.
+ */
+ virtual ~BFGS_Precond() {};
+
+ /** \brief
+ * realisation of Preconditioner::init
+ */
+ inline void init()
+ {
+ FUNCNAME("BFGS_Precond::init()");
+ TEST_EXIT(matrix[row])("no matrix\n");
+ };
+
+ /** \brief
+ * realisation of Preconditioner::precon
+ */
+ void precon(DOFVector<double> *vec);
+
+ /** \brief
+ * realisation of Preconditioner::exit
+ */
+ inline void exit() {};
+
+
+ /* ----- Specific functions ----- */
+
+ /** \brief
+ * add pair of vectors
+ */
+ void addPair(DOFVector<double> *s, DOFVector<double> *y, double sy = 0.0);
+
+ protected:
+
+ /** \brief
+ * list of pairs
+ */
+ std::deque<BFGS_Vectors> pairsVec;
+
+ /** \brief
+ * maximal number of pairs
+ */
+ int m;
+
+ /** \brief
+ * actual number of pairs
+ */
+ int k;
+ };
+
+}
+
+#endif // AMDIS_BFGS_PRECONDITIONER_H
diff --git a/AMDiS/src/QN_Precond.h b/AMDiS/src/QN_Precond.h
new file mode 100755
index 0000000000000000000000000000000000000000..edc038d085d592a635e21735a8e5aa525cec4094
--- /dev/null
+++ b/AMDiS/src/QN_Precond.h
@@ -0,0 +1,26 @@
+/** \file QN_Precond.h */
+
+#ifndef QN_PRECOND_H
+#define QN_PRECOND_H
+
+#include "DOFVector.h"
+
+namespace AMDiS
+{
+
+ template<typename T>
+ void diagPrecond(const DOFMatrix *A, DOFVector<T>* vec,
+ DOFVector<BoundaryType> *bound = NULL);
+
+ template<typename T>
+ void aXpY(double a, const DOFVector<T>& x, DOFVector<T>& y,
+ DOFVector<BoundaryType> *bound = NULL);
+
+ template<typename T>
+ void XpaY(double a, const DOFVector<T>& x, DOFVector<T>& y,
+ DOFVector<BoundaryType> *bound = NULL);
+
+}
+
+#include "QN_Precond.hh"
+#endif // QN_PRECOND_H
diff --git a/AMDiS/src/QN_Precond.hh b/AMDiS/src/QN_Precond.hh
new file mode 100755
index 0000000000000000000000000000000000000000..5bb4a35e865a4757c603d9dd3832df5a09780b1e
--- /dev/null
+++ b/AMDiS/src/QN_Precond.hh
@@ -0,0 +1,141 @@
+#include "DOFVector.h"
+
+namespace AMDiS
+{
+
+ template<typename T>
+ void diagPrecond(DOFMatrix *A, DOFVector<T> *vec,
+ DOFVector<BoundaryType> *bound)
+ {
+ FUNCNAME("QN_Precond::diagPrecond()");
+
+ TEST(A)("No matrix; doing nothing\n");
+
+ DOFMatrix::Iterator rowIterator(A, USED_DOFS);
+ typename DOFVector<T>::Iterator vecIterator(vec, USED_DOFS);
+
+ if (bound) {
+ DOFVector<BoundaryType>::Iterator
+ bIterator(const_cast<DOFVector<BoundaryType>*>(bound), USED_DOFS);
+
+ for (vecIterator.reset(), rowIterator.reset(), bIterator.reset();
+ !vecIterator.end(); ++vecIterator, ++bIterator, ++rowIterator)
+ {
+ // only non-dirichlet nodes will be preconditioned
+ if ((*bIterator) <= 0 && rowIterator->size() != 0)
+ (*vecIterator) *= (1.0 / (*rowIterator)[0].entry);
+ }
+ } else {
+ for (vecIterator.reset(), rowIterator.reset(); !vecIterator.end();
+ ++vecIterator, ++rowIterator)
+ {
+ if (rowIterator->size() != 0)
+ (*vecIterator) *= (1.0 / (*rowIterator)[0].entry);
+ }
+ }
+
+ return;
+ }
+
+ template<typename T>
+ void aXpY(double alpha, const DOFVector<T>& x, DOFVector<T>& y,
+ DOFVector<BoundaryType> *bound)
+ {
+ FUNCNAME("QN_Precond::aXpY()");
+ const DOFAdmin *admin;
+
+ TEST_EXIT(x.getFESpace() && y.getFESpace())
+ ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
+ TEST_EXIT((admin = x.getFESpace()->getAdmin()) &&
+ (admin == y.getFESpace()->getAdmin()))
+ ("no admin or different admins: %8X, %8X\n",
+ x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
+ TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
+ ("size = %d too small: admin->size = %d\n", x.getSize(),
+ admin->getUsedSize());
+ TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
+ ("y.size = %d too small: admin->size = %d\n", y.getSize(),
+ admin->getUsedSize());
+
+ typename DOFVector<T>::Iterator
+ xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)),
+ USED_DOFS);
+ typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y),
+ USED_DOFS);
+
+ if (bound) {
+ DOFVector<BoundaryType>::Iterator
+ bIterator(const_cast<DOFVector<BoundaryType>*>(bound), USED_DOFS);
+
+ for (xIterator.reset(), yIterator.reset(), bIterator.reset();
+ !xIterator.end();
+ ++xIterator, ++yIterator, ++bIterator)
+ {
+ // only non-dirichlet nodes will be updated
+ if ((*bIterator) <= 0 )
+ *yIterator += alpha * (*xIterator);
+ }
+ } else {
+ for (xIterator.reset(), yIterator.reset(); !xIterator.end();
+ ++xIterator, ++yIterator)
+ {
+ *yIterator += alpha * (*xIterator);
+ }
+ }
+
+ return;
+ }
+
+ template<typename T>
+ void XpaY(double alpha, const DOFVector<T>& x, DOFVector<T>& y,
+ DOFVector<BoundaryType> *bound)
+ {
+ FUNCNAME("QN_Precond::aXpY()");
+ const DOFAdmin *admin;
+
+ TEST_EXIT(x.getFESpace() && y.getFESpace())
+ ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
+ TEST_EXIT((admin = x.getFESpace()->getAdmin()) &&
+ (admin == y.getFESpace()->getAdmin()))
+ ("no admin or different admins: %8X, %8X\n",
+ x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
+ TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
+ ("size = %d too small: admin->size = %d\n", x.getSize(),
+ admin->getUsedSize());
+ TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
+ ("y.size = %d too small: admin->size = %d\n", y.getSize(),
+ admin->getUsedSize());
+
+ typename DOFVector<T>::Iterator
+ xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)),
+ USED_DOFS);
+ typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y),
+ USED_DOFS);
+
+ if (bound) {
+ DOFVector<BoundaryType>::Iterator
+ bIterator(const_cast<DOFVector<BoundaryType>*>(bound), USED_DOFS);
+
+ for (xIterator.reset(), yIterator.reset(), bIterator.reset();
+ !xIterator.end();
+ ++xIterator, ++yIterator, ++bIterator)
+ {
+ // only non-dirichlet nodes will be updated
+ if ((*bIterator) <= 0 ) {
+ *yIterator *= alpha;
+ *yIterator += *xIterator;
+ }
+ }
+ } else {
+ for (xIterator.reset(), yIterator.reset(); !xIterator.end();
+ ++xIterator, ++yIterator)
+ {
+ *yIterator *= alpha;
+ *yIterator += *xIterator;
+ }
+ }
+
+ return;
+ }
+
+}