conjgrad.c

Go to the documentation of this file.

#include <../../nrnconf.h>

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**               Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
    Conjugate gradient routines file
    Uses sparse matrix input & sparse Cholesky factorisation in pccg().

    All the following routines use routines to define a matrix
        rather than use any explicit representation
        (with the exeception of the pccg() pre-conditioner)
    The matrix A is defined by

        VEC *(*A)(void *params, VEC *x, VEC *y)

    where y = A.x on exit, and y is returned. The params argument is
    intended to make it easier to re-use & modify such routines.

    If we have a sparse matrix data structure
        SPMAT   *A_mat;
    then these can be used by passing sp_mv_mlt as the function, and
    A_mat as the param.
*/

#include    <stdio.h>
#include    <math.h>
#include    "matrix.h"
#include    "sparse.h"
static char rcsid[] = "conjgrad.c,v 1.1 1997/12/04 17:55:16 hines Exp";


/* #define  MAX_ITER    10000 */
static  int max_iter = 10000;
int cg_num_iters;

/* matrix-as-routine type definition */
/* #ifdef ANSI_C */
/* typedef VEC  *(*MTX_FN)(void *params, VEC *x, VEC *out); */
/* #else */
typedef VEC *(*MTX_FN)();
/* #endif */
#ifdef ANSI_C
VEC *spCHsolve(SPMAT *,VEC *,VEC *);
#else
VEC *spCHsolve();
#endif

/* cg_set_maxiter -- sets maximum number of iterations if numiter > 1
    -- just returns current max_iter otherwise
    -- returns old maximum */
int cg_set_maxiter(numiter)
int numiter;
{
    int temp;

    if ( numiter < 2 )
        return max_iter;
    temp = max_iter;
    max_iter = numiter;
    return temp;
}


/* pccg -- solves A.x = b using pre-conditioner M
            (assumed factored a la spCHfctr())
    -- results are stored in x (if x != NULL), which is returned */
VEC *pccg(A,A_params,M_inv,M_params,b,eps,x)
MTX_FN  A, M_inv;
VEC *b, *x;
double  eps;
void    *A_params, *M_params;
{
    VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
    int k;
    Real    alpha, beta, ip, old_ip, norm_b;

    if ( ! A || ! b )
        error(E_NULL,"pccg");
    if ( x == b )
        error(E_INSITU,"pccg");
    x = v_resize(x,b->dim);
    if ( eps <= 0.0 )
        eps = MACHEPS;

    r = v_get(b->dim);
    p = v_get(b->dim);
    q = v_get(b->dim);
    z = v_get(b->dim);

    norm_b = v_norm2(b);

    v_zero(x);
    r = v_copy(b,r);
    old_ip = 0.0;
    for ( k = 0; ; k++ )
    {
        if ( v_norm2(r) < eps*norm_b )
            break;
        if ( k > max_iter )
            error(E_ITER,"pccg");
        if ( M_inv )
            (*M_inv)(M_params,r,z);
        else
            v_copy(r,z);    /* M == identity */
        ip = in_prod(z,r);
        if ( k )    /* if ( k > 0 ) ... */
        {
            beta = ip/old_ip;
            p = v_mltadd(z,p,beta,p);
        }
        else        /* if ( k == 0 ) ... */
        {
            beta = 0.0;
            p = v_copy(z,p);
            old_ip = 0.0;
        }
        q = (*A)(A_params,p,q);
        alpha = ip/in_prod(p,q);
        x = v_mltadd(x,p,alpha,x);
        r = v_mltadd(r,q,-alpha,r);
        old_ip = ip;
    }
    cg_num_iters = k;

    V_FREE(p);
    V_FREE(q);
    V_FREE(r);
    V_FREE(z);

    return x;
}

/* sp_pccg -- a simple interface to pccg() which uses sparse matrix
        data structures
    -- assumes that LLT contains the Cholesky factorisation of the
        actual pre-conditioner */
VEC *sp_pccg(A,LLT,b,eps,x)
SPMAT   *A, *LLT;
VEC *b, *x;
double  eps;
{   return pccg(sp_mv_mlt,A,spCHsolve,LLT,b,eps,x);     }


/*
    Routines for performing the CGS (Conjugate Gradient Squared)
    algorithm of P. Sonneveld:
        "CGS, a fast Lanczos-type solver for nonsymmetric linear
        systems", SIAM J. Sci. & Stat. Comp. v. 10, pp. 36--52
*/

/* cgs -- uses CGS to compute a solution x to A.x=b
    -- the matrix A is not passed explicitly, rather a routine
        A is passed where A(x,Ax,params) computes
        Ax = A.x
    -- the computed solution is passed */
VEC *cgs(A,A_params,b,r0,tol,x)
MTX_FN  A;
VEC *x, *b;
VEC *r0;        /* tilde r0 parameter -- should be random??? */
double  tol;        /* error tolerance used */
void    *A_params;
{
    VEC *p, *q, *r, *u, *v, *tmp1, *tmp2;
    Real    alpha, beta, norm_b, rho, old_rho, sigma;
    int iter;

    if ( ! A || ! x || ! b || ! r0 )
        error(E_NULL,"cgs");
    if ( x->dim != b->dim || r0->dim != x->dim )
        error(E_SIZES,"cgs");
    if ( tol <= 0.0 )
        tol = MACHEPS;

    p = v_get(x->dim);
    q = v_get(x->dim);
    r = v_get(x->dim);
    u = v_get(x->dim);
    v = v_get(x->dim);
    tmp1 = v_get(x->dim);
    tmp2 = v_get(x->dim);

    norm_b = v_norm2(b);
    (*A)(A_params,x,tmp1);
    v_sub(b,tmp1,r);
    v_zero(p);  v_zero(q);
    old_rho = 1.0;

    iter = 0;
    while ( v_norm2(r) > tol*norm_b )
    {
        if ( ++iter > max_iter ) break;
        /*    error(E_ITER,"cgs");  */
        rho = in_prod(r0,r);
        if ( old_rho == 0.0 )
            error(E_SING,"cgs");
        beta = rho/old_rho;
        v_mltadd(r,q,beta,u);
        v_mltadd(q,p,beta,tmp1);
        v_mltadd(u,tmp1,beta,p);

        (*A)(A_params,p,v);

        sigma = in_prod(r0,v);
        if ( sigma == 0.0 )
            error(E_SING,"cgs");
        alpha = rho/sigma;
        v_mltadd(u,v,-alpha,q);
        v_add(u,q,tmp1);

        (*A)(A_params,tmp1,tmp2);

        v_mltadd(r,tmp2,-alpha,r);
        v_mltadd(x,tmp1,alpha,x);

        old_rho = rho;
    }
    cg_num_iters = iter;

    V_FREE(p);  V_FREE(q);  V_FREE(r);
    V_FREE(u);  V_FREE(v);
    V_FREE(tmp1);   V_FREE(tmp2);

    return x;
}

/* sp_cgs -- simple interface for SPMAT data structures */
VEC *sp_cgs(A,b,r0,tol,x)
SPMAT   *A;
VEC *b, *r0, *x;
double  tol;
{   return cgs(sp_mv_mlt,A,b,r0,tol,x); }

/*
    Routine for performing LSQR -- the least squares QR algorithm
    of Paige and Saunders:
        "LSQR: an algorithm for sparse linear equations and
        sparse least squares", ACM Trans. Math. Soft., v. 8
        pp. 43--71 (1982)
*/
/* lsqr -- sparse CG-like least squares routine:
    -- finds min_x ||A.x-b||_2 using A defined through A & AT
    -- returns x (if x != NULL) */
VEC *lsqr(A,AT,A_params,b,tol,x)
MTX_FN  A, AT;  /* AT is A transposed */
VEC *x, *b;
double  tol;        /* error tolerance used */
void    *A_params;
{
    VEC *u, *v, *w, *tmp;
    Real    alpha, beta, norm_b, phi, phi_bar,
                rho, rho_bar, rho_max, theta;
    Real    s, c;   /* for Givens' rotations */
    int iter, m, n;

    if ( ! b || ! x )
        error(E_NULL,"lsqr");
    if ( tol <= 0.0 )
        tol = MACHEPS;

    m = b->dim; n = x->dim;
    u = v_get((u_int)m);
    v = v_get((u_int)n);
    w = v_get((u_int)n);
    tmp = v_get((u_int)n);
    norm_b = v_norm2(b);

    v_zero(x);
    beta = v_norm2(b);
    if ( beta == 0.0 )
        return x;
    sv_mlt(1.0/beta,b,u);
    tracecatch((*AT)(A_params,u,v),"lsqr");
    alpha = v_norm2(v);
    if ( alpha == 0.0 )
        return x;
    sv_mlt(1.0/alpha,v,v);
    v_copy(v,w);
    phi_bar = beta;     rho_bar = alpha;

    rho_max = 1.0;
    iter = 0;
    do {
        if ( ++iter > max_iter )
            error(E_ITER,"lsqr");

        tmp = v_resize(tmp,m);
        tracecatch((*A) (A_params,v,tmp),"lsqr");

        v_mltadd(tmp,u,-alpha,u);
        beta = v_norm2(u);  sv_mlt(1.0/beta,u,u);

        tmp = v_resize(tmp,n);
        tracecatch((*AT)(A_params,u,tmp),"lsqr");
        v_mltadd(tmp,v,-beta,v);
        alpha = v_norm2(v); sv_mlt(1.0/alpha,v,v);

        rho = sqrt(rho_bar*rho_bar+beta*beta);
        if ( rho > rho_max )
            rho_max = rho;
        c   = rho_bar/rho;
        s   = beta/rho;
        theta   =  s*alpha;
        rho_bar = -c*alpha;
        phi     =  c*phi_bar;
        phi_bar =  s*phi_bar;

        /* update x & w */
        if ( rho == 0.0 )
            error(E_SING,"lsqr");
        v_mltadd(x,w,phi/rho,x);
        v_mltadd(v,w,-theta/rho,w);
    } while ( fabs(phi_bar*alpha*c) > tol*norm_b/rho_max );

    cg_num_iters = iter;

    V_FREE(tmp);    V_FREE(u);  V_FREE(v);  V_FREE(w);

    return x;
}

/* sp_lsqr -- simple interface for SPMAT data structures */
VEC *sp_lsqr(A,b,tol,x)
SPMAT   *A;
VEC *b, *x;
double  tol;
{   return lsqr(sp_mv_mlt,sp_vm_mlt,A,b,tol,x); }