# to unbundle, sh this file (in an empty directory) echo conjgrad.c 1>&2 sed >conjgrad.c <<'//GO.SYSIN DD conjgrad.c' 's/^-//' - -/************************************************************************** -** -** 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 -#include -#include "matrix.h" -#include "sparse.h" -static char rcsid[] = "$Id: conjgrad.c,v 1.4 1994/01/13 05:36:45 des 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); } - //GO.SYSIN DD conjgrad.c echo lanczos.c 1>&2 sed >lanczos.c <<'//GO.SYSIN DD lanczos.c' 's/^-//' - -/************************************************************************** -** -** 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. -** -***************************************************************************/ - - -/* - File containing Lanczos type routines for finding eigenvalues - of large, sparse, symmetic matrices -*/ - -#include -#include -#include "matrix.h" -#include "sparse.h" - -static char rcsid[] = "$Id: lanczos.c,v 1.4 1994/01/13 05:28:24 des Exp $"; - -#ifdef ANSI_C -extern VEC *trieig(VEC *,VEC *,MAT *); -#else -extern VEC *trieig(); -#endif - -/* lanczos -- raw lanczos algorithm -- no re-orthogonalisation - -- creates T matrix of size == m, - but no larger than before beta_k == 0 - -- uses passed routine to do matrix-vector multiplies */ -void lanczos(A_fn,A_params,m,x0,a,b,beta2,Q) -VEC *(*A_fn)(); /* VEC *(*A_fn)(void *A_params,VEC *in, VEC *out) */ -void *A_params; -int m; -VEC *x0, *a, *b; -Real *beta2; -MAT *Q; -{ - int j; - VEC *v, *w, *tmp; - Real alpha, beta; - - if ( ! A_fn || ! x0 || ! a || ! b ) - error(E_NULL,"lanczos"); - if ( m <= 0 ) - error(E_BOUNDS,"lanczos"); - if ( Q && ( Q->m < x0->dim || Q->n < m ) ) - error(E_SIZES,"lanczos"); - - a = v_resize(a,(u_int)m); b = v_resize(b,(u_int)(m-1)); - v = v_get(x0->dim); - w = v_get(x0->dim); - tmp = v_get(x0->dim); - - beta = 1.0; - /* normalise x0 as w */ - sv_mlt(1.0/v_norm2(x0),x0,w); - - (*A_fn)(A_params,w,v); - - for ( j = 0; j < m; j++ ) - { - /* store w in Q if Q not NULL */ - if ( Q ) - set_col(Q,j,w); - - alpha = in_prod(w,v); - a->ve[j] = alpha; - v_mltadd(v,w,-alpha,v); - beta = v_norm2(v); - if ( beta == 0.0 ) - { - v_resize(a,(u_int)j+1); - v_resize(b,(u_int)j); - *beta2 = 0.0; - if ( Q ) - Q = m_resize(Q,Q->m,j+1); - return; - } - if ( j < m-1 ) - b->ve[j] = beta; - v_copy(w,tmp); - sv_mlt(1/beta,v,w); - sv_mlt(-beta,tmp,v); - (*A_fn)(A_params,w,tmp); - v_add(v,tmp,v); - } - *beta2 = beta; - - - V_FREE(v); V_FREE(w); V_FREE(tmp); -} - -extern double frexp(), ldexp(); - -/* product -- returns the product of a long list of numbers - -- answer stored in mant (mantissa) and expt (exponent) */ -static double product(a,offset,expt) -VEC *a; -double offset; -int *expt; -{ - Real mant, tmp_fctr; - int i, tmp_expt; - - if ( ! a ) - error(E_NULL,"product"); - - mant = 1.0; - *expt = 0; - if ( offset == 0.0 ) - for ( i = 0; i < a->dim; i++ ) - { - mant *= frexp(a->ve[i],&tmp_expt); - *expt += tmp_expt; - if ( ! (i % 10) ) - { - mant = frexp(mant,&tmp_expt); - *expt += tmp_expt; - } - } - else - for ( i = 0; i < a->dim; i++ ) - { - tmp_fctr = a->ve[i] - offset; - tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset : - MACHEPS*offset; - mant *= frexp(tmp_fctr,&tmp_expt); - *expt += tmp_expt; - if ( ! (i % 10) ) - { - mant = frexp(mant,&tmp_expt); - *expt += tmp_expt; - } - } - - mant = frexp(mant,&tmp_expt); - *expt += tmp_expt; - - return mant; -} - -/* product2 -- returns the product of a long list of numbers - -- answer stored in mant (mantissa) and expt (exponent) */ -static double product2(a,k,expt) -VEC *a; -int k; /* entry of a to leave out */ -int *expt; -{ - Real mant, mu, tmp_fctr; - int i, tmp_expt; - - if ( ! a ) - error(E_NULL,"product2"); - if ( k < 0 || k >= a->dim ) - error(E_BOUNDS,"product2"); - - mant = 1.0; - *expt = 0; - mu = a->ve[k]; - for ( i = 0; i < a->dim; i++ ) - { - if ( i == k ) - continue; - tmp_fctr = a->ve[i] - mu; - tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu; - mant *= frexp(tmp_fctr,&tmp_expt); - *expt += tmp_expt; - if ( ! (i % 10) ) - { - mant = frexp(mant,&tmp_expt); - *expt += tmp_expt; - } - } - mant = frexp(mant,&tmp_expt); - *expt += tmp_expt; - - return mant; -} - -/* dbl_cmp -- comparison function to pass to qsort() */ -static int dbl_cmp(x,y) -Real *x, *y; -{ - Real tmp; - - tmp = *x - *y; - return (tmp > 0 ? 1 : tmp < 0 ? -1: 0); -} - -/* lanczos2 -- lanczos + error estimate for every e-val - -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978 - -- returns multiple e-vals where multiple e-vals may not exist - -- returns evals vector */ -VEC *lanczos2(A_fn,A_params,m,x0,evals,err_est) -VEC *(*A_fn)(); -void *A_params; -int m; -VEC *x0; /* initial vector */ -VEC *evals; /* eigenvalue vector */ -VEC *err_est; /* error estimates of eigenvalues */ -{ - VEC *a; - static VEC *b=VNULL, *a2=VNULL, *b2=VNULL; - Real beta, pb_mant, det_mant, det_mant1, det_mant2; - int i, pb_expt, det_expt, det_expt1, det_expt2; - - if ( ! A_fn || ! x0 ) - error(E_NULL,"lanczos2"); - if ( m <= 0 ) - error(E_RANGE,"lanczos2"); - - a = evals; - a = v_resize(a,(u_int)m); - b = v_resize(b,(u_int)(m-1)); - MEM_STAT_REG(b,TYPE_VEC); - - lanczos(A_fn,A_params,m,x0,a,b,&beta,MNULL); - - /* printf("# beta =%g\n",beta); */ - pb_mant = 0.0; - if ( err_est ) - { - pb_mant = product(b,(double)0.0,&pb_expt); - /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */ - } - - /* printf("# diags =\n"); out_vec(a); */ - /* printf("# off diags =\n"); out_vec(b); */ - a2 = v_resize(a2,a->dim - 1); - b2 = v_resize(b2,b->dim - 1); - MEM_STAT_REG(a2,TYPE_VEC); - MEM_STAT_REG(b2,TYPE_VEC); - for ( i = 0; i < a2->dim - 1; i++ ) - { - a2->ve[i] = a->ve[i+1]; - b2->ve[i] = b->ve[i+1]; - } - a2->ve[a2->dim-1] = a->ve[a2->dim]; - - trieig(a,b,MNULL); - - /* sort evals as a courtesy */ - qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp); - - /* error estimates */ - if ( err_est ) - { - err_est = v_resize(err_est,(u_int)m); - - trieig(a2,b2,MNULL); - /* printf("# a =\n"); out_vec(a); */ - /* printf("# a2 =\n"); out_vec(a2); */ - - for ( i = 0; i < a->dim; i++ ) - { - det_mant1 = product2(a,i,&det_expt1); - det_mant2 = product(a2,(double)a->ve[i],&det_expt2); - /* printf("# det_mant1=%g, det_expt1=%d\n", - det_mant1,det_expt1); */ - /* printf("# det_mant2=%g, det_expt2=%d\n", - det_mant2,det_expt2); */ - if ( det_mant1 == 0.0 ) - { /* multiple e-val of T */ - err_est->ve[i] = 0.0; - continue; - } - else if ( det_mant2 == 0.0 ) - { - err_est->ve[i] = HUGE; - continue; - } - if ( (det_expt1 + det_expt2) % 2 ) - /* if odd... */ - det_mant = sqrt(2.0*fabs(det_mant1*det_mant2)); - else /* if even... */ - det_mant = sqrt(fabs(det_mant1*det_mant2)); - det_expt = (det_expt1+det_expt2)/2; - err_est->ve[i] = fabs(beta* - ldexp(pb_mant/det_mant,pb_expt-det_expt)); - } - } - - return a; -} - -/* sp_lanczos -- version that uses sparse matrix data structure */ -void sp_lanczos(A,m,x0,a,b,beta2,Q) -SPMAT *A; -int m; -VEC *x0, *a, *b; -Real *beta2; -MAT *Q; -{ lanczos(sp_mv_mlt,A,m,x0,a,b,beta2,Q); } - -/* sp_lanczos2 -- version of lanczos2() that uses sparse matrix data - structure */ -VEC *sp_lanczos2(A,m,x0,evals,err_est) -SPMAT *A; -int m; -VEC *x0; /* initial vector */ -VEC *evals; /* eigenvalue vector */ -VEC *err_est; /* error estimates of eigenvalues */ -{ return lanczos2(sp_mv_mlt,A,m,x0,evals,err_est); } - //GO.SYSIN DD lanczos.c echo arnoldi.c 1>&2 sed >arnoldi.c <<'//GO.SYSIN DD arnoldi.c' 's/^-//' - -/************************************************************************** -** -** 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 -#include -#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; -} //GO.SYSIN DD arnoldi.c .