vic/vic_run/src/newt_raph_func_fast.c
int
newt_raph(void (*vecfunc)(double x[], double fvec[], int n, int init, ...),
          double x[],
          int n)
{
    extern parameters_struct param;

    int                      k, i, Error;
    double                   errx, errf, fvec[MAX_NODES], p[MAX_NODES];
    double                   a[MAX_NODES], b[MAX_NODES], c[MAX_NODES];

    Error = 0;

    for (k = 0; k < param.NEWT_RAPH_MAXTRIAL; k++) {
        // calculate function value for all nodes, i.e. focus = -1
        (*vecfunc)(x, fvec, n, 0, -1);

        // stop if TOLF is satisfied
        errf = 0.0;
        for (i = 0; i < n; i++) {
            errf += fabs(fvec[i]);
        }
        if (errf <= param.NEWT_RAPH_TOLF) {
            return (Error);
        }

        // calculate the Jacobian
        fdjac3(x, fvec, a, b, c, vecfunc, n);

        for (i = 0; i < n; i++) {
            p[i] = -fvec[i];
        }

        // only the tri-diagnal part of the Jacobian matrix is significant
        // and most of the off-belt entries are encessially zeros
        tridiag(a, b, c, p, n);

        errx = 0.0;
        for (i = 0; i < n; i++) {
            errx += fabs(p[i]);

            if (k > 10 && k <= 20 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX1;
            }
            else if (k > 20 && k <= 60 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                     param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX2;
            }
            else if (k > 60 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                     param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX3;
            }
            else {
                x[i] += p[i];
            }
        }

        // stop if TOLX is satisfied
        if (errx <= param.NEWT_RAPH_TOLX) {
            return (Error);
        }
    }
    Error = 1;

    return (Error);
}
/******************************************************************************
 * @brief    Newton-Raphson method to solve non-linear system adapted from
 *           "Numerical Recipes"
 *****************************************************************************/
@brief    Newton-Raphson method to solve non-linear system adapted from
"Numerical Recipes"
@brief Newton - Raphson method to solve non - linear system adapted from " Numerical Recipes "
int
newt_raph(void (*vecfunc)(double x[], double fvec[], int n, int init, ...),
          double x[],
          int n)
{
    extern parameters_struct param;
    int                      k, i, Error;
    double                   errx, errf, fvec[MAX_NODES], p[MAX_NODES];
    double                   a[MAX_NODES], b[MAX_NODES], c[MAX_NODES];
    Error = 0;
    for (k = 0; k < param.NEWT_RAPH_MAXTRIAL; k++) {
        (*vecfunc)(x, fvec, n, 0, -1);
        errf = 0.0;
        for (i = 0; i < n; i++) {
            errf += fabs(fvec[i]);
        }
        if (errf <= param.NEWT_RAPH_TOLF) {
            return (Error);
        }
        fdjac3(x, fvec, a, b, c, vecfunc, n);
        for (i = 0; i < n; i++) {
            p[i] = -fvec[i];
        }
        tridiag(a, b, c, p, n);
        errx = 0.0;
        for (i = 0; i < n; i++) {
            errx += fabs(p[i]);
            if (k > 10 && k <= 20 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX1;
            }
            else if (k > 20 && k <= 60 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                     param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX2;
            }
            else if (k > 60 && x[i] < param.NEWT_RAPH_R_MAX && x[i] >
                     param.NEWT_RAPH_R_MIN) {
                x[i] += p[i] * param.NEWT_RAPH_RELAX3;
            }
            else {
                x[i] += p[i];
            }
        }
        if (errx <= param.NEWT_RAPH_TOLX) {
            return (Error);
        }
    }
    Error = 1;
    return (Error);
}
int newt_raph ( void ( * vecfunc ) ( double x [ ] , double fvec [ ] , int n , int init , ... ) , double x [ ] , int n ) { extern parameters_struct param ; int k , i , Error ; double errx , errf , fvec [ MAX_NODES ] , p [ MAX_NODES ] ; double a [ MAX_NODES ] , b [ MAX_NODES ] , c [ MAX_NODES ] ; Error = 0 ; for ( k = 0 ; k < param . NEWT_RAPH_MAXTRIAL ; k ++ ) { ( * vecfunc ) ( x , fvec , n , 0 , -1 ) ; errf = 0.0 ; for ( i = 0 ; i < n ; i ++ ) { errf += fabs ( fvec [ i ] ) ; } if ( errf <= param . NEWT_RAPH_TOLF ) { return ( Error ) ; } fdjac3 ( x , fvec , a , b , c , vecfunc , n ) ; for ( i = 0 ; i < n ; i ++ ) { p [ i ] = - fvec [ i ] ; } tridiag ( a , b , c , p , n ) ; errx = 0.0 ; for ( i = 0 ; i < n ; i ++ ) { errx += fabs ( p [ i ] ) ; if ( k > 10 && k <= 20 && x [ i ] < param . NEWT_RAPH_R_MAX && x [ i ] > param . NEWT_RAPH_R_MIN ) { x [ i ] += p [ i ] * param . NEWT_RAPH_RELAX1 ; } else if ( k > 20 && k <= 60 && x [ i ] < param . NEWT_RAPH_R_MAX && x [ i ] > param . NEWT_RAPH_R_MIN ) { x [ i ] += p [ i ] * param . NEWT_RAPH_RELAX2 ; } else if ( k > 60 && x [ i ] < param . NEWT_RAPH_R_MAX && x [ i ] > param . NEWT_RAPH_R_MIN ) { x [ i ] += p [ i ] * param . NEWT_RAPH_RELAX3 ; } else { x [ i ] += p [ i ] ; } } if ( errx <= param . NEWT_RAPH_TOLX ) { return ( Error ) ; } } Error = 1 ; return ( Error ) ; }
@brief    Newton-Raphson method to solve non-linear system adapted from
           "Numerical Recipes"
@brief Newton - Raphson method to solve non - linear system adapted from " Numerical Recipes "
// calculate function value for all nodes, i.e. focus = -1 // stop if TOLF is satisfied // calculate the Jacobian // only the tri-diagnal part of the Jacobian matrix is significant // and most of the off-belt entries are encessially zeros // stop if TOLX is satisfied