extrapolation.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2011 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef EXTRAPOLATION_HH
#define EXTRAPOLATION_HH
 
#include <vector>
 
namespace Kaskade
{
  template <class T>
  class ExtrapolationTableau
  {
  public:
    typedef T value_type;
 
    explicit ExtrapolationTableau(double hTarget_): hTarget(hTarget_) {}
 
 
    void clear() { data.clear(); }
 
    T const& back() const { return data.back(); }
 
    T const& operator[](int i) const {
      assert(0<=i && i<data.size());
      return data[i];
    }
 
    int size() const { return data.size(); }
 
 
    void push_back(T t, double hnew)
    {
      T tmp(t);
 
      for (int i=0; i<data.size(); ++i) {
        // at loop start: t = t(n+1,i), data[i] = t(n,i)
 
        // compute coefficients, see Deuflhard/Hohmann 7.1.2
        int n = data.size()-1;
        double a = - (hTarget-hnew)/(hnew-h[n-i]);
        double b = (hTarget-h[n-i])/(hnew-h[n-i]);
        // compute tmp = t(n+1,i+1) = a*t(n,i)+b*t(n+1,i)
        tmp = t; tmp *= b;
        data[i] *= a; tmp += data[i];
        // store data[i] = t(n+1,i)
        data[i] = t;
        // maintain loop invariant: t = t(n+1,i+1)
        t = tmp;
      }
      // store new t(n+1,n+1)
      data.push_back(t);
      h.push_back(hnew);
    }
 
 
 
  private:
    // The array data always contains a complete row of the extrapolation tableau.
    std::vector<T>      data;
    std::vector<double> h;
    double              hTarget;
  };
} // namespace Kaskade
 
#endif