RISC DISC 1/RISC_DISC_1.iso/pd_share/code/meschach/!Meschach/c/arnoldi

Text File | 1994-01-13 | 4.7 KB | 188 lines

/**************************************************************************
**
** 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.
**
***************************************************************************/

/*
    Arnoldi method for finding eigenvalues of large non-symmetric
        matrices
*/
#include    <stdio.h>
#include    <math.h>
#include    "matrix.h"
#include    "matrix2.h"
#include    "sparse.h"

static char rcsid[] = "$Id: arnoldi.c,v 1.3 1994/01/13 05:45:40 des Exp $";

/* arnoldi -- an implementation of the Arnoldi method */
MAT    *arnoldi(A,A_param,x0,m,h_rem,Q,H)
VEC    *(*A)();
void    *A_param;
VEC    *x0;
int    m;
Real    *h_rem;
MAT    *Q, *H;
{
    static VEC    *v=VNULL, *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
    int    i;
    Real    h_val;

    if ( ! A || ! Q || ! x0 )
        error(E_NULL,"arnoldi");
    if ( m <= 0 )
        error(E_BOUNDS,"arnoldi");
    if ( Q->n != x0->dim ||    Q->m != m )
        error(E_SIZES,"arnoldi");

    m_zero(Q);
    H = m_resize(H,m,m);
    m_zero(H);
    u = v_resize(u,x0->dim);
    v = v_resize(v,x0->dim);
    r = v_resize(r,m);
    s = v_resize(s,m);
    tmp = v_resize(tmp,x0->dim);
    MEM_STAT_REG(u,TYPE_VEC);
    MEM_STAT_REG(v,TYPE_VEC);
    MEM_STAT_REG(r,TYPE_VEC);
    MEM_STAT_REG(s,TYPE_VEC);
    MEM_STAT_REG(tmp,TYPE_VEC);
    sv_mlt(1.0/v_norm2(x0),x0,v);
    for ( i = 0; i < m; i++ )
    {
        set_row(Q,i,v);
        u = (*A)(A_param,v,u);
        r = mv_mlt(Q,u,r);
        tmp = vm_mlt(Q,r,tmp);
        v_sub(u,tmp,u);
        h_val = v_norm2(u);
        /* if u == 0 then we have an exact subspace */
        if ( h_val == 0.0 )
        {
        *h_rem = h_val;
        return H;
        }
        /* iterative refinement -- ensures near orthogonality */
        do {
        s = mv_mlt(Q,u,s);
        tmp = vm_mlt(Q,s,tmp);
        v_sub(u,tmp,u);
        v_add(r,s,r);
        } while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
        /* now that u is nearly orthogonal to Q, update H */
        set_col(H,i,r);
        if ( i == m-1 )
        {
        *h_rem = h_val;
        continue;
        }
        /* H->me[i+1][i] = h_val; */
        m_set_val(H,i+1,i,h_val);
        sv_mlt(1.0/h_val,u,v);
    }

    return H;
}

/* sp_arnoldi -- uses arnoldi() with an explicit representation of A */
MAT    *sp_arnoldi(A,x0,m,h_rem,Q,H)
SPMAT    *A;
VEC    *x0;
int    m;
Real    *h_rem;
MAT    *Q, *H;
{    return arnoldi(sp_mv_mlt,A,x0,m,h_rem,Q,H);    }

/* gmres -- generalised minimum residual algorithm of Saad & Schultz
        SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
    -- y is overwritten with the solution */
VEC    *gmres(A,A_param,m,Q,R,b,tol,x)
VEC    *(*A)();
void    *A_param;
VEC    *b, *x;
int    m;
MAT    *Q, *R;
double    tol;
{
    static VEC    *v=VNULL, *u=VNULL, *r=VNULL, *tmp=VNULL, *rhs=VNULL;
    static VEC    *diag=VNULL, *beta=VNULL;
    int    i;
    Real    h_val, norm_b;

    if ( ! A || ! Q || ! b || ! R )
    error(E_NULL,"gmres");
    if ( m <= 0 )
    error(E_BOUNDS,"gmres");
    if ( Q->n != b->dim || Q->m != m )
    error(E_SIZES,"gmres");

    x = v_copy(b,x);
    m_zero(Q);
    R = m_resize(R,m+1,m);
    m_zero(R);
    u = v_resize(u,x->dim);
    v = v_resize(v,x->dim);
    tmp = v_resize(tmp,x->dim);
    rhs = v_resize(rhs,m+1);
    MEM_STAT_REG(u,TYPE_VEC);
    MEM_STAT_REG(v,TYPE_VEC);
    MEM_STAT_REG(r,TYPE_VEC);
    MEM_STAT_REG(tmp,TYPE_VEC);
    MEM_STAT_REG(rhs,TYPE_VEC);
    norm_b = v_norm2(x);
    if ( norm_b == 0.0 )
    error(E_RANGE,"gmres");
    sv_mlt(1.0/norm_b,x,v);

    for ( i = 0; i < m; i++ )
    {
    set_row(Q,i,v);
    tracecatch(u = (*A)(A_param,v,u),"gmres");
    r = mv_mlt(Q,u,r);
    tmp = vm_mlt(Q,r,tmp);
    v_sub(u,tmp,u);
    h_val = v_norm2(u);
    set_col(R,i,r);
    R->me[i+1][i] = h_val;
    sv_mlt(1.0/h_val,u,v);
    }

    /* use i x i submatrix of R */
    R = m_resize(R,i+1,i);
    rhs = v_resize(rhs,i+1);
    v_zero(rhs);
    rhs->ve[0] = norm_b;
    tmp = v_resize(tmp,i);
    diag = v_resize(diag,i+1);
    beta = v_resize(beta,i+1);
    MEM_STAT_REG(beta,TYPE_VEC);
    MEM_STAT_REG(diag,TYPE_VEC);
    QRfactor(R,diag /* ,beta */);
    tmp = QRsolve(R,diag, /* beta, */ rhs,tmp);
    v_resize(tmp,m);
    vm_mlt(Q,tmp,x);

    return x;
}