partialderivative.h

Go to the documentation of this file.

#ifndef PARTIALDERIVATIVE_H
#define PARTIALDERIVATIVE_H

#include <cmath>
using namespace std;

#include <typeop.h>

template< typename EXP >
class partialderivative
{
public:

  typename typeop<EXP>::Tref exp;

  typedef typename typeop<EXP>::Tbare::FNtype FN;
  typedef typename typeop<EXP>::Tbare::Ttype T;
  
  T * dxi;

  uintc N;

  partialderivative
  (
    typename typeop<EXP>::Tref exp_
  );

  ~partialderivative();

  void operator()();

  T const squared() const
  {
    T val(0);
    for (uint i=0; i<N; ++i)
      val += dxi[i]*dxi[i];

    return val;
  }
 
  T const absvalue() const
  {
    T val(0);
    for (uint i=0; i<N; ++i)
      val += abs(dxi[i]);

    return val;
  }
};

//----------------------------------------------------
//  Implementation


template< typename EXP >
void partialderivative<EXP>::operator()()
{
  T x;
  T fnb;
  T fna;

  T h;

  // O(N).O(fn) complexity
  for (uint i=0; i<N; ++i)
  {
    h = exp.hi[i];
    x = exp.fn.xi[i];
    exp.fn.xi[i] -= h;
    exp.fn(fna);
    exp.fn.xi[i] = x + h;
    exp.fn(fnb);
    // Restore xi
    exp.fn.xi[i] = x;
    
    // Compute the partial derivative.
    assert(h!=0.0);
    dxi[i] = (fnb-fna)/h*0.5;
  }

  // Restore the current state.
  exp.fn(fna);
}

template< typename EXP >
partialderivative<EXP>::partialderivative
(
  typename typeop<EXP>::Tref exp_
)
  : exp(exp_), N(exp_.N)
{
  dxi = new T[N];
}

template< typename EXP >
partialderivative<EXP>::~partialderivative()
{
  delete[] dxi;
}




#endif