AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

frprmn.h

Go to the documentation of this file.

#ifndef NR_FRPRMN
#define NR_FRPRMN

#include <iostream>
#include "utils.h"
#include "linmin.h"
#include "dlinmin.h"

namespace nr {

template < class TraitsT = Old_Traits<> >
struct frprmn
{

    typedef TraitsT Traits;
    typedef typename Traits::Real Real;
    typedef typename Traits::Vector Vector;

    
    template < class Func, class DFunc >
    frprmn (
        Vector& p, 
        unsigned int n,
        Real ftol, 
        unsigned int *iter, 
        Real *fret,
        Func *funct,
        DFunc *gradient,
        unsigned int ITMAX = 200,
        Real TOL = 2.0e-4,
        Real EPS = 1.0e-10
    )
    {
        std::size_t j,its;
        Real gg,gam,fp,dgg;
        Vector g, h, xi;

        Traits::alloc_vector( g, n );
        Traits::alloc_vector( h, n );
        Traits::alloc_vector( xi, n );

        fp=(*funct)(p);
        (*gradient)(p,xi);
        for (j=0;j<n;j++) {
            g[j] = -xi[j];
            xi[j]=h[j]=g[j];
        }
        for (its=1;its<=ITMAX;its++) {
            *iter=its;
            dlinmin<Traits>(p,xi,n,fret,funct,gradient,TOL);
            //linmin<Traits>(p,xi,n,fret,funct,TOL);
            //std::cout<<std::endl<<" ===============frprmn: fp = "<< fp <<", fret = "<< *fret << ", p = " << p << std::endl;
            if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
                //std::cerr << "frprmn: precision reached " << endl;
                Traits::free_vector( g );
                Traits::free_vector( h );
                Traits::free_vector( xi );
                //cerr << "frprmn:: vectors deallocated" << endl;
                return;
            }
            //cerr <<"continue frprmn " << endl;
            fp=(*funct)(p);
            (*gradient)(p,xi);
            dgg=gg=0.0;
            //cerr <<"continue frprmn after next computation" << endl;
            for (j=0;j<n;j++) {
                gg += g[j]*g[j];
                dgg += (xi[j]+g[j])*xi[j];
            }
            //cerr << "frprmn:: linear combination done " << endl;
            if ( gg == 0.0 ) {
                //cerr << "frprmn: gg==0.0 " << endl;
                Traits::free_vector( g );
                Traits::free_vector( h );
                Traits::free_vector( xi );
                //cerr << "frprmn:: vectors deallocated" << endl;
                return;
            }
            gam=dgg/gg;
            for (j=0;j<n;j++) {
                g[j] = -xi[j];
                xi[j]=h[j]=g[j]+gam*h[j];
            }
        //cerr << "frprmn:: end loop " << endl;
        }
//      nrerror("Too many iterations in frprmn");
        //throw TooManyIterations( ITMAX, *fret );
    }
    
    

    
    // Exception: too many iterations
    struct TooManyIterations 
    {
        unsigned int nbIterations;  
        Real functionValue;         

        TooManyIterations ( unsigned int n, Real f )
            : nbIterations( n )
            , functionValue( f )
        {}
        

        void print ( std::ostream& str ) const
        {
            str << "frprmn gives up after " << nbIterations
                << " iterations, function value = " << functionValue
                << std::endl;
        }
        
        
    };

};


} // end namespace nr


#endif

---