// $Id: nonlinear.h,v 1.1.1.1 2000/11/02 08:47:17 fritzi Exp $

#ifndef __nonlinear_h
#define __nonlinear_h

#ifndef __iteration_h
#include <deal/iteration.h>
#endif
#include <math.h>

//////////
template<class MT, class VT, class MEM>
inline int nl(MT& A, VT& u, MEM& mem, IterationInfo& info)
{
  int it,reached=0;
  VT& r = mem.v(0);
  VT& d = mem.v(1);

  double res_old;
  double res;

  res     =  A.nl_residual(r,u);
  res_old =  res;
  reached =  info.check_residual(0,res);

  for (it=1;!reached;it++)
  {
    A.nl_assemble(u);      // compute Newton-Matrix
    A.nl_solve(d,r);       // search direction

    double a_old=0.;
    int    damped=0;
    for (int i=0;(i<10)&&(!damped);i++) {
      double a = pow(0.5,i);
      u.add(a-a_old,d);

      res = A.nl_residual(r,u);
      if (res<res_old)
	damped=1;
      a_old = a;
    }
    reached = info.check_residual(it,res);
    res_old = res;
  }
  if (reached<0) return 1;
  return 0;
}

//////////
template<class MT, class VT, class MEM>
inline int nonlinear(MT& A, VT& u, MEM& mem, IterationInfo& info)
{
  int it,reached=0;
  VT& v =mem.v(0);
  VT& d =mem.v(1);

  for (it=0;!reached;it++)
  {
    double res  =  A.nl_residual(d,u);
    reached     =  info.check_residual(it,res);
    if (reached) break;
    v.reinit(u);
    A.nl_update(u,v,d);
  }
  if (reached<0) return 1;
  return 0;
}

//////////
template<class MT>
inline int nonlinear_light(MT& A, IterationInfo& info)
{
  int it,reached=0;

  info.reset();
  for (it=0;!reached;it++)
  {
    double res  =  A.nll_residual();
    reached     =  info.check_residual(it,res);
    if (reached) break;
    A.nll_update();
  }
  return reached;
}

#endif