---

guma_newton.h

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#ifndef GUMA_NEWTON_H
#define GUMA_NEWTON_H

namespace guma {

using gul::Ptr;

/*
#define ALF  1.0e-4
#define TOLX 1.0e-7
#define MAXITS 200
#define TOLF 1.0e-4
#define TOLMIN 1.0e-6
#define STPMAX 100.0
*/


template< class T >
class newton_info : public gul::pool_object
{
public:
  int dim;
  T alf;    // average rate of decrease
  T tolx;   // convergence in domain
  
  gul::Ptr<T>             fvec;
  gul::Ptr< gul::Ptr<T> > fjac;

  newton_info( int in_dim, T in_alf, T in_tolx )
  {
    int i;
    
    dim = in_dim;
    alf = in_alf;
    tolx = in_tolx;

    fvec.reserve_pool(dim);
    fjac.reserve_pool(dim);
    for( i = 0; i < dim; i++ )
      fjac[i].reserve_pool(dim);
  }
  virtual ~newton_info() { }

  /*----------------------------------------------------------------------
    user provided, has to calculate the function value and Jacobian, and
    store it in fvec and fjac. It has to return, if convergence is
    achieved, too. 
  -----------------------------------------------------------------------*/
  virtual bool evaluate( const gul::Ptr<T>& dom ) = 0;

  /*------------------------------------
    calculates F, and returns 1/2*F*F
  ------------------------------------*/
  bool calcf( const gul::Ptr<T>& dom, T *f )
  {
    T sum;
    int i;
    bool found;

    found = evaluate(dom);

    for( sum = (T)0, i = 0; i < dim; i++ )
      sum += fvec[i] * fvec[i];
    
    *f = (T)0.5 * sum;

    return found;
  }
};

/*------------------------------------------------------------------------
  Searches a root of a function, with n-dimensional domain and 
  n-dimensional function values. (see "Numerical Recipes in C")
  in:
    domain:           [x1[0],x2[0]] * ... * [x1[n-1],x2[n-1]]
    start value:      (x[0],...,x[n-1])
  out:
    root:             (x[0],...,x[n-1])
-----------------------------------------------------------------------*/
template< class T >
bool Newton(
         Ptr<T>& x,
         const Ptr<T>& x1,
         const Ptr<T>& x2,
         newton_info<T>& I,
         int maxits );

}

#endif

---