C语言写矩阵库，适合做矩阵运算

/**************************************************************************
**
** 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 < m; 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