---

guma_newton.cpp

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.
 */

#include "stdafx.h"

#include <iostream>

#include "gul_types.h"
#include "gul_matrix.h"
#include "guma_lineq.h"
#include "guma_newton.h"


namespace guma {

using gul::Ptr;
using gul::rtr;
using gul::Max;
using gul::MCopy;

/********************************************************************

  This file contains Newton-Iteration related functions,
  taken from "Numerical Recipes in C"

*********************************************************************/

/*-------------------------------------------------------------------
  Search along a direction 'P' for a point, where function values have
  decreased sufficiently (<ALF)
---------------------------------------------------------------------*/

template< class T >
bool LinearSearch( 
              const Ptr<T>& p,     /* Newton direction (domain) */
              const Ptr<T>& xold,  /* starting point (domain) */
              T fold,              /* function value at start point */
              const Ptr<T>& g,     /* gradient of function at 'xold' */
              const Ptr<T>& x1,    /* bounds for x */
              const Ptr<T>& x2,
              newton_info<T>& I,   /* defines function and tolerances */
 
              bool&    found,
              bool&    check,      /* returns if convergence must be checked */
                                   /* (i.e. |x_{i+1}-x_{i}| < TOLX)          */
              Ptr<T>&  x,          /* return: point in domain */
              T       *fx          /* return: function value there */
            )
{
  T sum, slope, test, temp, alamin, alam, alphamin, alpha, stpmax;
  T a, b, alam2, disc, f2, fold2, rhs1, rhs2, tmplam;
  int i, n;
  T f;
  
  found = false;
  n = I.dim;

  // maximale Schrittweite bestimmen

  alphamin = 1.0;
  
  for( i = 0; i < n; i++ )
  {
    if( xold[i] + p[i] < x1[i] )
      alpha = (x1[i] - xold[i]) / p[i];
    else if( xold[i] + p[i] > x2[i] )
      alpha = (x2[i] - xold[i]) / p[i];
    else
      alpha = 1.0;

    if( rtr<T>::fabs(alpha) == 0.0 )
    {
      p[i] = 0.0;
      alpha = 1.0;   
    }
  
    if( alpha < alphamin )
      alphamin = alpha;  
  }        
    
  for( sum = 0.0, i = 0; i < n; i++ )
    sum += p[i] * p[i];

  stpmax = rtr<T>::fabs(alphamin) * rtr<T>::sqrt(sum);

  /*-------------------------------------------------------------*/

  for( sum = 0.0, i = 0; i < n; i++ )
    sum += p[i] * p[i];
  sum = rtr<T>::sqrt( sum );  
  
  if( sum > stpmax )
    for( i = 0; i < n; i++ )
      p[i] *= stpmax / sum;      /* Skalieren, falls Schritt zu gross */
      
  for( slope = 0.0, i = 0; i < n; i++ )
    slope += g[i] * p[i];      
      
  test = 0.0;
  
  for( i = 0; i < n; i++ )
  {
    temp = rtr<T>::fabs(p[i]) / Max(rtr<T>::fabs(xold[i]),(T)1);     
    if( temp > test ) test = temp;
  }
  if( test == (T)0 ) test = rtr<T>::tiny();
  
  alamin = I.tolx / test;
  alam = (T)1;              /* Im ersten Durchgang volle Schrittweite */
      
  for( ;; )
  {
    for( i = 0; i < n; i++ )            
    {
      x[i] = xold[i] + alam * p[i];

      if( x[i] < x1[i] ) x[i] = x1[i];
      else if( x[i] > x2[i] ) x[i] = x2[i];
    } 

    found = I.calcf( x, &f );
    if( found ) { check = false; break; }
    
    if( alam < alamin )           /* Konvergenz von dx,                     */
    {                             /* aufrufendes Programm sollte Konvergenz */
      for( i = 0; i < n; i++ )    /* verifizieren                           */
        x[i] = xold[i];
      check = true;
      break;  
    }

    if( f <= fold + I.alf * alam * slope )  /* Funktionswert nimmt stark */
    { check = false;  break; }              /* genug ab */

    if( alam == (T)1 )                             /* falls 1.Durchgang */
      tmplam = -slope / ((T)2 * (f-fold-slope));
    else  
    {
      rhs1 = f - fold - alam * slope;
      rhs2 = f2 - fold2 - alam2 * slope;
      a = (rhs1/(alam*alam) - rhs2/(alam2*alam2)) / (alam-alam2);
      b = (-alam2*rhs1/(alam*alam) + alam*rhs2/(alam2*alam2)) / (alam-alam2);
       
      if( a == (T)0 )
        tmplam = -slope/((T)2*b);
      else
      {
        disc = b*b - (T)3*a*slope;

        if( disc < (T)0 )     /* roundoff problem */
          return false;    

        tmplam = (-b + rtr<T>::sqrt(disc)) / ((T)3*a);
      }
      
      if( tmplam > (T)0.5 * alam )
        tmplam = (T)0.5 * alam;        /* lambda <= 0.5 * lambda1  */            
    }        
    
    alam2 = alam;
    f2 = f;
    fold2 = fold;
    alam = Max( tmplam, (T)0.1*alam );   /* lambda >= 0.1 * lambda1  */ 
  }  
    
  *fx = f;
  return true;  
}    
// template instantiations
template bool LinearSearch( 
              const Ptr<float>& p,     /* Newton direction (domain) */
              const Ptr<float>& xold,  /* starting point (domain) */
              float fold,              /* function value at start point */
              const Ptr<float>& g,     /* gradient of function at 'xold' */
              const Ptr<float>& x1,    /* bounds for x */
              const Ptr<float>& x2,   
              newton_info<float>& I,   /* defines function and tolerances */

              bool& found,
              bool& check,         /* returns if convergence must be checked */
                                   /* (i.e. |x_{i+1}-x_{i}| < TOLX)          */
              Ptr<float>& x,           /* return: point in domain */
              float *fx                /* return: function value there */
            );



/*------------------------------------------------------------------------
  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 )
{
  int its, i, j, n;
  T f, fold, sum, temp, test;
  bool result, found, check;
  Ptr<int> indx;
  Ptr<T> g,p,xold;
  Ptr< Ptr<T> > A;
  int sign;

  n = I.dim;

  A.reserve_pool(n);
  for( i = 0; i < n; i++ )
    A[i].reserve_pool(n);

  indx.reserve_pool(n);
  g.reserve_pool(n);
  p.reserve_pool(n);
  xold.reserve_pool(n);
     
  /* ------- Hauptschleife: ---------------------------------------- */  

  found = I.calcf( x, &f );
  if( found ) { return true; }

  for( its = 1; its < maxits; its++ )
  { 
    for( i = 0; i < n; i++ )              /* Gradient von f=1/2*F*F */
    {                                     /* berechnen (=> df = F*J */
      for( sum = 0.0, j = 0; j < n; j++ )
        sum += I.fvec[j] * I.fjac[j][i];
      g[i] = sum;  
    }

    for( i = 0; i < n; i++ )
      xold[i] = x[i];
    fold = f;
    
    for( i = 0; i < n; i++ )
      p[i] = -I.fvec[i];              /* rechte Seite der Gleichung */

    MCopy(n,n,I.fjac,A);
    result = LUDecomposition( n, A, &sign, indx );
    if( !result ) return(false);
    LUBackSubstitution( n, indx, A, p );

    /* ------------------------------------------------------------------- */

    result = LinearSearch( p, xold, fold, g, x1, x2, I, found, check, x, &f );  
    // if( (!result) || check ) return false;
    if( !result ) return false;

    // convergence of dx may happen when there is no better solution then a
    // point on the border of the surface, so for the time beeing i am
    // simply ignoring the 'check' flag (i am not sure, if this is allright
    // in all cases)

    if( found ) { return true; }

    test = 0.0;
    for( i = 0; i < n; i++ )
    {
      temp = (rtr<T>::fabs(x[i]-xold[i])) / Max(rtr<T>::fabs(x[i]),(T)1);
      if( temp > test ) test = temp;
    }  
    if( test < I.tolx ) return true;
  }

  /* maximum iterations exceeded */
  return false;
}    
// template instantiations
template bool Newton(
         Ptr<float>& x,
         const Ptr<float>& x1,
         const Ptr<float>& x2,
         newton_info<float>& I,
         int maxits );
template bool Newton(
         Ptr<double>& x,
         const Ptr<double>& x1,
         const Ptr<double>& x2,
         newton_info<double>& I,
         int maxits );







#if 0

*******************************************************************
old version
*******************************************************************


/*------------------------------------------------------------------------
  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]] x ... x [x1[n-1],x2[n-1]]
    start value:      (x[0],...,x[n-1])
    function values:  fvec
    Jacobi Matrix:    fjac
  out:
    root:             (x[0],...,x[n-1])
    check:            if TRUE(!=0), try again with another start value
-----------------------------------------------------------------------*/

template< class T >
bool Newton(
         Ptr<T>& x,
         const Ptr<T>& x1,
         const Ptr<T>& x2,
         newton_info<T>& I,
         int maxits,
         bool& check
       )
{
  T a;
  int its, i, j, n;
  T den, f, fold, stpmax, sum, temp, test;
  T alpha, alphamin;
  bool result;
  Ptr<int> indx;
  Ptr<T> g,p,xold;
  Ptr< Ptr<T> > A;
  int sign;

  n = I.dim;

  A.reserve_pool(n);
  for( i = 0; i < n; i++ )
    A[i].reserve_pool(n);

  indx.reserve_pool(n);
  g.reserve_pool(n);
  p.reserve_pool(n);
  xold.reserve_pool(n);

  f = I.calcf( x );

  for( test = (T)0, i = 0; i < n; i++ ) /* Pruefen ob Startwert Nullstelle */  
  {
    a = rtr<T>::fabs(I.fvec[i]);
    if( a > test ) test = a;
  }
  if( test < (T)0.01 * I.tolf )             /* falls ja */
  {
    check = false;
    return true;    
  }      

  /*
  // maximale Schrittweite fuer LinearSearch bestimmen:
  for( sum = 0.0, i = 0; i < n; i++ )
  sum += x[i] * x[i];
    
  stpmax = STPMAX * NUAR_MAX( sqrt(sum), (NURational)n );
  */
     
  /* ------- Hauptschleife: ---------------------------------------- */  

  for( its = 1; its < maxits; its++ )
  {
    I.Fjac( x );     /* Jakobi Matrix berechnen */  

    for( i = 0; i < n; i++ )              /* Gradient von f=1/2*F*F */
    {                                     /* berechnen (=> df = F*J */
      for( sum = 0.0, j = 0; j < n; j++ )
        sum += I.fvec[j] * I.fjac[j][i];
      g[i] = sum;  
    }

    for( i = 0; i < n; i++ )
      xold[i] = x[i];
    fold = f;
    
    for( i = 0; i < n; i++ )
      p[i] = -I.fvec[i];              /* rechte Seite der Gleichung */

    MCopy(n,n,I.fjac,A);
    result = LUDecomposition( n, A, &sign, indx );
    if( !result ) return(false);
    LUBackSubstitution( n, indx, A, p );

    /* --- maximale Schrittweite fuer LinearSearch bestimmen: ----------- */ 

    alphamin = 1.0;
  
    for( i = 0; i < n; i++ )
    {
      if( xold[i] + p[i] < x1[i] )
        alpha = (x1[i] - xold[i]) / p[i];
      else if( xold[i] + p[i] > x2[i] )
        alpha = (x2[i] - xold[i]) / p[i];
      else
        alpha = 1.0;

      if( rtr<T>::fabs(alpha) == 0.0 )
      {
        p[i] = 0.0;
        alpha = 1.0;   
      }
  
      if( alpha < alphamin )
        alphamin = alpha;  
    }        
    
    for( sum = 0.0, i = 0; i < n; i++ )
      sum += p[i] * p[i];

    stpmax = rtr<T>::fabs(alphamin) * rtr<T>::sqrt(sum);

    /* ------------------------------------------------------------------- */

    result = LinearSearch( p, xold, fold, g, stpmax, I, check, x, &f );  
    if( !result ) return false;
    
    for( i = 0; i < n; i++ )
    {
      if( x[i] < x1[i] )
        x[i] = x1[i];
      else if( x[i] > x2[i] )
        x[i] = x2[i];  
    }

    for( test = 0.0, i = 0; i < n; i++ ) /* Pruefen: Funktionswerte konvergent*/
      if( fabs(I.fvec[i]) > test )
        test = rtr<T>::fabs(I.fvec[i]);
    if( test < 0.01 * I.tolf )  
    {
      check = false;              /* falls ja, Erfolg, Ende */
      return true;    
    }      

    if( check )                /* Pruefen ob Gradient von f = 1/2*F*F = 0 ist,*/
    {                          /* falls ja, besteht die Moeglichkeit, dass das*/
      test = 0.0;              /* Minimum von f nicht die gesuchte Nullstelle */
      den = Max( f, (T)0.5*(T)n );   /* von F ist !!!!!!                     */
      
      for( i = 0; i < n; i++ )
      {
        temp = rtr<T>::fabs(g[i])*Max(rtr<T>::fabs(x[i]), (T)1) / den;
        if( temp > test ) test = temp;
      }
      check = (test < I.tolmin ? true : false);
      return true;   
    }

    test = 0.0;
    for( i = 0; i < n; i++ )
    {
      temp = (rtr<T>::fabs(x[i]-xold[i])) / Max(rtr<T>::fabs(x[i]),(T)1);
      if( temp > test ) test = temp;
    }  
    if( test < I.tolx ) return true;
  }

  /* maximum iterations exceeded */
  return false;
}    
// template instantiations
template bool Newton(
         Ptr<float>& x,
         const Ptr<float>& x1,
         const Ptr<float>& x2,
         newton_info<float>& I,
         int maxits );
#endif



}

---