vic/vic_run/src/newt_raph_func_fast.c
void
tridiag(double   a[],
        double   b[],
        double   c[],
        double   r[],
        unsigned n)
{
    unsigned i;
    int      j;
    double   factor;

    /* forward substitution */
    factor = b[0];
    b[0] = 1.0;
    c[0] = c[0] / factor;
    r[0] = r[0] / factor;

    for (i = 1; i < n; i++) {
        factor = a[i];
        a[i] = a[i] - b[i - 1] * factor;
        b[i] = b[i] - c[i - 1] * factor;
        r[i] = r[i] - r[i - 1] * factor;

        factor = b[i];
        b[i] = 1.0;
        c[i] = c[i] / factor;
        r[i] = r[i] / factor;
    }

    /* backward substitution */
    for (j = n - 2; j >= 0; j--) {
        factor = c[j];
        c[j] = c[j] - b[j + 1] * factor;
        r[j] = r[j] - r[j + 1] * factor;

        factor = b[j];
        r[j] = r[j] / factor;
    }
}
/******************************************************************************
 * @brief    function to solve tridiagonal linear system adapted from
 *           "Numerical Recipes"
 *****************************************************************************/
@brief    function to solve tridiagonal linear system adapted from
"Numerical Recipes"
@brief function to solve tridiagonal linear system adapted from " Numerical Recipes "
void
tridiag(double   a[],
        double   b[],
        double   c[],
        double   r[],
        unsigned n)
{
    unsigned i;
    int      j;
    double   factor;
    factor = b[0];
    b[0] = 1.0;
    c[0] = c[0] / factor;
    r[0] = r[0] / factor;
    for (i = 1; i < n; i++) {
        factor = a[i];
        a[i] = a[i] - b[i - 1] * factor;
        b[i] = b[i] - c[i - 1] * factor;
        r[i] = r[i] - r[i - 1] * factor;
        factor = b[i];
        b[i] = 1.0;
        c[i] = c[i] / factor;
        r[i] = r[i] / factor;
    }
    for (j = n - 2; j >= 0; j--) {
        factor = c[j];
        c[j] = c[j] - b[j + 1] * factor;
        r[j] = r[j] - r[j + 1] * factor;
        factor = b[j];
        r[j] = r[j] / factor;
    }
}
void tridiag ( double a [ ] , double b [ ] , double c [ ] , double r [ ] , unsigned n ) { unsigned i ; int j ; double factor ; factor = b [ 0 ] ; b [ 0 ] = 1.0 ; c [ 0 ] = c [ 0 ] / factor ; r [ 0 ] = r [ 0 ] / factor ; for ( i = 1 ; i < n ; i ++ ) { factor = a [ i ] ; a [ i ] = a [ i ] - b [ i - 1 ] * factor ; b [ i ] = b [ i ] - c [ i - 1 ] * factor ; r [ i ] = r [ i ] - r [ i - 1 ] * factor ; factor = b [ i ] ; b [ i ] = 1.0 ; c [ i ] = c [ i ] / factor ; r [ i ] = r [ i ] / factor ; } for ( j = n - 2 ; j >= 0 ; j -- ) { factor = c [ j ] ; c [ j ] = c [ j ] - b [ j + 1 ] * factor ; r [ j ] = r [ j ] - r [ j + 1 ] * factor ; factor = b [ j ] ; r [ j ] = r [ j ] / factor ; } }
@brief    function to solve tridiagonal linear system adapted from
           "Numerical Recipes"
@brief function to solve tridiagonal linear system adapted from " Numerical Recipes "
/* forward substitution */ /* backward substitution */