/**************************************************************************
**
**?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.cv?1.3?1994/01/13?05:45:40?des?Exp?$“;


/*?arnoldi?--?an?implementation?of?the?Arnoldi?method?*/
MAT	*arnoldi（AA_paramx0mh_remQH）
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（Hmm）;
	m_zero（H）;
	u?=?v_resize（ux0->dim）;
	v?=?v_resize（vx0->dim）;
	r?=?v_resize（rm）;
	s?=?v_resize（sm）;
	tmp?=?v_resize（tmpx0->dim）;
	MEM_STAT_REG（uTYPE_VEC）;
	MEM_STAT_REG（vTYPE_VEC）;
	MEM_STAT_REG（rTYPE_VEC）;
	MEM_STAT_REG（sTYPE_VEC）;
	MEM_STAT_REG（tmpTYPE_VEC）;
	sv_mlt（1.0/v_norm2（x0）x0v）;
	for?（?i?=?0;?i?	{
	????set_row（Qiv）;
	????u?=?（*A）（A_paramvu）;
	????r?=?mv_mlt（Qur）;
	????tmp?=?vm_mlt（Qrtmp）;
	????v_sub（utmpu）;
	????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（Qus）;
		tmp?=?vm_mlt（Qstmp）;
		v_sub（utmpu）;
		v_add（rsr）;
	????}?while?（?v_norm2（s）?>?0.1*（h_val?=?v_norm2（u））?）;
	????/*?now?that?u?is?nearly?orthogonal?to?Q?update?H?*/
	????set_col（Hir）;
	????if?（?i?==?m-1?）
	????{
		*h_rem?=?h_val;
		continue;
	????}
	????/*?H->me[i+1][i]?=?h_val;?*/
	????m_set_val（Hi+1ih_val）;
	????sv_mlt（1.0/h_valuv）;
	}

#ifde