mes_eig.h

矩阵运算的模板类

// -*- c++ -*-
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2001 Oh-Wook Kwon, all rights reserved. [e-mail]
//
//                          Easy Matrix Template Library
// 
// This Easy Matrix Template 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
// Easy Matrix Template 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.
//
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// Filename: mes_eig.h
// Revision:
//    1. Revised to support complex matrices.
// Note:
//    1. I strongly recommend the routines in this file be called 
//       with double precision arguments.
///////////////////////////////////////////////////////////////////////////////

#ifndef	_MES_EIG_H_
#define	_MES_EIG_H_	

/* ----- Eigen Analysis for Symmetric Matrices and General Matrices ----- */
/**************************************************************************
**
** Copyright (C) 1993 David E. Stewart & 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.
**
***************************************************************************/


template <typename T> inline T Matrix<T>::mes_in_prod(const Vector<T>& a,const Vector<T>& b) {
	return mes_in_prod_offset(a,b,1);
}


/* in_prod_ -- inner product of two vectors from i0 downwards */
template <typename T> T Matrix<T>::mes_in_prod_offset(const Vector<T>& a,const Vector<T>& b,int i0){
	int	i,limit;
	limit = local_min(a.Size(),b.Size());
	assert ( i0 <= limit );
	T sum=T();
	for(i=i0;i<=limit;i++) sum += conj(a[i])*b[i];
	return sum;
}


/* mes_norm2_offset -- norm of vectors from i0 downwards */
template <typename T> double Matrix<T>::mes_norm2_offset(const Vector<T>& a,int i0){
	int	i,limit;
	limit = a.Size();
	assert ( i0 >= 1 && i0 <= limit );
	double sum=0, z=0;
	for(i=i0;i<=limit;i++) if(z<Abs(a[i])) z=Abs(a[i]);
	if(z==0) return 0;
	for(i=i0;i<=limit;i++) {
		double x=Abs(a[i]/T(z));
		sum += x*x;
	}
	return (z*sqrt(sum));
}


/* _v_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */
template <typename T> double	Matrix<T>::mes_norm_inf(const Vector<T>& x){
	return Norm(x,mtl_numeric_limits<int>::max());
}


/* v_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */
template <typename T> double	Matrix<T>::mes_norm2(const Vector<T>& x){
	return Norm(x,2);
}


/* get_row -- gets a specified row of a matrix and retruns it as a vector */
template <typename T> void	Matrix<T>::mes_get_row(const Matrix<T>& mat,int row,Vector<T>& vec){
   int	i;
   assert ( row <= NumRows(mat) );
   if ( vec.Size() != NumCols(mat) )
     vec.Resize(NumCols(mat));
   
   for ( i=1; i<=NumCols(mat); i++ )
     vec[i] = mat[row][i];
   
   return;
}


/* get_col -- gets a specified column of a matrix and retruns it as a vector */
template <typename T> void	Matrix<T>::mes_get_col(const Matrix<T>& mat,int col,Vector<T>& vec){
   int	i;
   assert ( col <= NumCols(mat) );
   if ( vec.Size() != NumRows(mat) )
     vec.Resize(NumRows(mat));
   
   for ( i=1; i<=NumRows(mat); i++ )
     vec[i] = mat[i][col];
   
   return;
}


/* set_row -- sets row of matrix to values given in vec (in situ) */
template <typename T> void	Matrix<T>::mes_set_row(Matrix<T>& mat,int row,Vector<T>& vec){
   int	i,lim;
   
   assert( row <= NumRows(mat) );
   lim = local_min(NumCols(mat),vec.Size());
   for ( i=1; i<=lim; i++ )
     mat[row][i] = vec[i];

}


/* set_col -- sets column of matrix to values given in vec (in situ) */
template <typename T> void	Matrix<T>::mes_set_col(Matrix<T>& mat,int col,Vector<T>& vec){
   int	i,lim;
   
   assert( col <= NumCols(mat) );
   lim = local_min(NumRows(mat),vec.Size());
   for ( i=1; i<=lim; i++ )
     mat[i][col] = vec[i];

}


/* v_copy_ -- copies vector into new area */
template <typename T> void	Matrix<T>::mes_copy_offset(const Vector<T>& in,Vector<T>& out,int i0){
	if ( &in[1]==&out[1] )
		return;
	if ( out.Size() != in.Size() )
		out.Resize(in.Size());
	for (int i=i0; i <= in.Size(); i++ )
		out[i] = in[i];
}


///////////////////////////////////////////////////////////////////////////////
// trieig -- finds eigenvalues of symmetric tridiagonal matrices
//	-- matrix represented by a pair of vectors a (diag entries)
//		and b (sub- & super-diag entries)
//	-- eigenvalues in a on return	
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_trieig(Vector<T>& a,Vector<T>& b,Matrix<T>& Q,int maxIter){
	int	i, i_min, i_max, n, split, iter, isOK=1;
	double	bk, ak1, bk1, ak2, bk2, z;
	double	c, c2, cs, s, s2, d, mu;
	static double  local_SQRT2=1.4142135623730949;

	if ( a.Size() != b.Size() + 1 || ( NumRows(Q) != a.Size() ) )
		mes_error(E_SIZES,"trieig");
	if ( NumRows(Q) != NumCols(Q) )
		mes_error(E_SQUARE,"trieig");

	n = a.Size();

	i_min = 1;
	while ( i_min <= n ){		/* outer while loop */
		/* find i_max to suit;
			submatrix i_min..i_max should be irreducible */
		i_max = n;
		for ( i = i_min; i < n; i++ )
			if ( b[i] == 0.0 ){
				i_max = i;	break;
			}
		if ( i_max <= i_min ){
		    i_min = i_max + 1;
		    continue;	/* outer while loop */
		}

		/* repeatedly perform QR method until matrix splits */
		split = false;
		iter=0;
		while ( ! split ){		/* inner while loop */
			int sign_d;
			iter++;

			/* find Wilkinson shift */
			d = (a[i_max-1] - a[i_max])/2;
			sign_d = (d >= 0 ? 1 : -1 );
			mu = a[i_max] - (b[i_max-1]*b[i_max-1])/(d + sign_d*pythagoras(d,b[i_max-1]));
			
			/* initial Givens' rotation */
			mes_givens(a[i_min]-mu,b[i_min],&c,&s);

			if ( Abs(c) < local_SQRT2 ){
				c2 = c*c;	s2 = 1-c2;
			}
			else{	
				s2 = s*s;	c2 = 1-s2;
			}
			cs = c*s;
			ak1 = c2*a[i_min]+s2*a[i_min+1]-2*cs*b[i_min];
			bk1 = cs*(a[i_min]-a[i_min+1]) + (c2-s2)*b[i_min];
			ak2 = s2*a[i_min]+c2*a[i_min+1]+2*cs*b[i_min];
			bk2 = ( i_min < i_max-1 ) ? c*b[i_min+1] : 0.0;
			z  = ( i_min < i_max-1 ) ? -s*b[i_min+1] : 0.0;
			a[i_min] = ak1;
			a[i_min+1] = ak2;
			b[i_min] = bk1;
			if ( i_min < i_max-1 )
				b[i_min+1] = bk2;
			if ( NumCols(Q) )
				mes_rot_cols(Q,i_min,i_min+1,c,s,Q);
			
			for ( i = i_min+1; i < i_max; i++ ){
				/* get Givens' rotation for sub-block -- k == i-1 */
				mes_givens(b[i-1],z,&c,&s);
				
				/* perform Givens' rotation on sub-block */
				if ( Abs(c) < local_SQRT2 ){
					c2 = c*c;	s2 = 1-c2;
				}
				else{
					s2 = s*s;	c2 = 1-s2;
				}
				cs = c*s;
				bk  = c*b[i-1] - s*z;
				ak1 = c2*a[i]+s2*a[i+1]-2*cs*b[i];
				bk1 = cs*(a[i]-a[i+1]) +
					(c2-s2)*b[i];
				ak2 = s2*a[i]+c2*a[i+1]+2*cs*b[i];
				bk2 = ( i+1 < i_max ) ? c*b[i+1] : 0.0;
				z  = ( i+1 < i_max ) ? -s*b[i+1] : 0.0;
				a[i] = ak1;	a[i+1] = ak2;
				b[i] = bk1;
				if ( i < i_max-1 )
					b[i+1] = bk2;
				if ( i > i_min )
					b[i-1] = bk;
				if ( NumCols(Q) )
					mes_rot_cols(Q,i,i+1,c,s,Q);
			}
			
			/* test to see if matrix should be split */
			for ( i = i_min; i < i_max; i++ ){
				if ( (Abs(b[i])) < (MACH_EPS*(Abs(a[i])+Abs(a[i+1]))) ){
					b[i] = 0.0;
					split = true;
				}
			}
	
			if(iter==maxIter){
				cerr << "ERROR: Too many iterations in mes_trieig().\n";
				isOK=0;
				assert(0);
				return isOK;
			}
		}
	}
	
	return isOK;
}


///////////////////////////////////////////////////////////////////////////////
// symmeig -- computes eigenvalues of a dense symmetric matrix
//	-- A **must** be symmetric on entry
//	-- eigenvalues stored in out
//	-- Q contains orthogonal matrix of eigenvectors
//	-- returns vector of eigenvalues
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_symmeig(const Matrix<T>& A,Matrix<T>& QQ,Vector<T>& oo){
	int	i, isOK;
	static Vector<double> bb, diag, beta;

	if(A.IsComplex())
		mes_error(E_NEED_REAL,"mes_symmeig");
	if(A.IsDouble()==false)
		mes_warn(WARN_SINGLE_PRECISION,"mes_symmeig");

	Matrix<double> AA=A;
	if ( NumRows(AA) != NumCols(AA) )
		mes_error(E_SQUARE,"symmeig");
	if ( oo.Size() != NumRows(AA) )
		oo.Resize(NumRows(AA));

	bb  .Resize(NumRows(AA) - 1);
	diag.Resize((int)NumRows(AA));
	beta.Resize((int)NumRows(AA));

	isOK=Matrix<double>::mes_Hfactor(AA,diag,beta);
	if(isOK==0) return isOK;
	if ( NumCols(QQ) ) {
		Matrix<double>::mes_makeHQ(AA,diag,beta,QQ);
		if(isOK==0) return isOK;
	}

	for ( i = 1; i <= NumRows(AA) - 1; i++ ){
		oo[i] = AA[i][i];
		bb[i] = AA[i][i+1];
	}
	oo[i] = AA[i][i];
	isOK=Matrix<double>::mes_trieig(oo,bb,QQ);

	return isOK;
}


///////////////////////////////////////////////////////////////////////////////
// Schur decomposition of a real non-symmetric matrix
///////////////////////////////////////////////////////////////////////////////

template <typename T> inline void	Matrix<T>::mes_hhldr3(double x, double y, double z, double *nu1, double *beta, double *newval){
	double	alpha;

	if ( x >= 0.0 )
		alpha = pythagoras3(x,y,z);
	else
		alpha = -pythagoras3(x,y,z);
	*nu1 = x + alpha;
	*beta = 1.0/(alpha*(*nu1));
	*newval = alpha;
}


template <typename T> 	int	Matrix<T>::mes_hhldr3cols(Matrix<T>& A, int k, int j0, double beta, double nu1, double nu2, double nu3){
	double	ip, prod;
	int	j, n;

	if ( k <= 0 || k+3 > NumRows(A)+1 || j0 <= 0 )
		mes_error(E_BOUNDS,"hhldr3cols");
	n = NumCols(A);

	for ( j = j0; j <= n; j++ ){
	    ip = nu1*A(k,j)+nu2*A(k+1,j)+nu3*A(k+2,j);
	    prod = ip*beta;	    
		A[k][j] += (- prod*nu1);
	    A[k+1][j] += (- prod*nu2);
	    A[k+2][j] += (- prod*nu3);
	}
	return 1;
}


template <typename T> 	int	Matrix<T>::mes_hhldr3rows(Matrix<T>& A, int k, int i0, double beta, double nu1, double nu2, double nu3){
	double	ip, prod;
	int	i, m;

	if ( k <= 0 || k+3 > NumCols(A)+1 )
		mes_error(E_BOUNDS,"hhldr3rows");
	m = NumRows(A);
	i0 = local_min(i0,m);

	for ( i = 1; i <= i0; i++ ){
	    ip = nu1*A(i,k)+nu2*A(i,k+1)+nu3*A(i,k+2);
	    prod = ip*beta;
		A[i][k] += (- prod*nu1);
	    A[i][k+1] += (- prod*nu2);
	    A[i][k+2] += (- prod*nu3);
	}
	return 1;