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

    int                      j;
    double                   h, temp, f[MAX_NODES];

    for (j = 0; j < n; j++) {
        temp = x[j];
        h = param.NEWT_RAPH_EPS2 * fabs(temp);
        if (h == 0) {
            h = param.NEWT_RAPH_EPS2;
        }
        x[j] = temp + h;
        h = x[j] - temp;

        // only update column j-1, j and j+1, caused by change in x[j]
        (*vecfunc)(x, f, n, 0, j);

        x[j] = temp;

        b[j] = (f[j] - fvec[j]) / h;
        if (j != 0) {
            c[j - 1] = (f[j - 1] - fvec[j - 1]) / h;
        }
        if (j != n - 1) {
            a[j + 1] = (f[j + 1] - fvec[j + 1]) / h;
        }
    }
}
/******************************************************************************
 * @brief    Forward difference approx to Jacobian, adapted from "Numerical
 *           Recipes"
 *****************************************************************************/
@brief    Forward difference approx to Jacobian, adapted from "Numerical
Recipes"
@brief Forward difference approx to Jacobian adapted from " Numerical Recipes "
void
fdjac3(double x[],
       double fvec[],
       double a[],
       double b[],
       double c[],
       void (*vecfunc)(double x[], double fvec[], int n, int init, ...),
       int n)
{
    extern parameters_struct param;
    int                      j;
    double                   h, temp, f[MAX_NODES];
    for (j = 0; j < n; j++) {
        temp = x[j];
        h = param.NEWT_RAPH_EPS2 * fabs(temp);
        if (h == 0) {
            h = param.NEWT_RAPH_EPS2;
        }
        x[j] = temp + h;
        h = x[j] - temp;
        (*vecfunc)(x, f, n, 0, j);
        x[j] = temp;
        b[j] = (f[j] - fvec[j]) / h;
        if (j != 0) {
            c[j - 1] = (f[j - 1] - fvec[j - 1]) / h;
        }
        if (j != n - 1) {
            a[j + 1] = (f[j + 1] - fvec[j + 1]) / h;
        }
    }
}
void fdjac3 ( double x [ ] , double fvec [ ] , double a [ ] , double b [ ] , double c [ ] , void ( * vecfunc ) ( double x [ ] , double fvec [ ] , int n , int init , ... ) , int n ) { extern parameters_struct param ; int j ; double h , temp , f [ MAX_NODES ] ; for ( j = 0 ; j < n ; j ++ ) { temp = x [ j ] ; h = param . NEWT_RAPH_EPS2 * fabs ( temp ) ; if ( h == 0 ) { h = param . NEWT_RAPH_EPS2 ; } x [ j ] = temp + h ; h = x [ j ] - temp ; ( * vecfunc ) ( x , f , n , 0 , j ) ; x [ j ] = temp ; b [ j ] = ( f [ j ] - fvec [ j ] ) / h ; if ( j != 0 ) { c [ j - 1 ] = ( f [ j - 1 ] - fvec [ j - 1 ] ) / h ; } if ( j != n - 1 ) { a [ j + 1 ] = ( f [ j + 1 ] - fvec [ j + 1 ] ) / h ; } } }
@brief    Forward difference approx to Jacobian, adapted from "Numerical
           Recipes"
@brief Forward difference approx to Jacobian adapted from " Numerical Recipes "
// only update column j-1, j and j+1, caused by change in x[j]