* and now the files I have forgotten in the last revision :)

AMDiS/src/BFGS_Precond.cc

+#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;
+  }
+
+}

AMDiS/src/BFGS_Precond.h

+// ============================================================================
+// ==                                                                        ==
+// == 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

AMDiS/src/QN_Precond.h

+/** \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

AMDiS/src/QN_Precond.hh

+#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;
+  }
+
+}