# 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	<stdio.h>
-#include	<math.h>
-#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	<stdio.h>
-#include	<math.h>
-#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	<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;
-}
//GO.SYSIN DD arnoldi.c