mes_factor.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_factor.h
// Revision:
//    1. Revised to support complex matrix factorization.
///////////////////////////////////////////////////////////////////////////////

#ifndef	_MES_FACTOR_H_
#define	_MES_FACTOR_H_	

/**************************************************************************
**
** 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.
**
***************************************************************************/

/* Most matrix factorisation routines are in-situ unless otherwise specified */
/* Matrix factorisation routines to work with the other matrix files.     */
/* Cholesky decomposition: A=R'*R                                         */
/* LU factorization: A=L*U                                                */
/* QR factorization: A=Q*R                                                */

///////////////////////////////////////////////////////////////////////////////
// CHfactor -- Cholesky L.L' factorisation of A in-situ
// Computes the Cholesky decompostion for a positive-definite symmetric matrix.
// A = L * L^T
// The L is returned in the lower triangle of A.
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_CHfactor(Matrix<T>& A){
	int	i, j, k, n;
	T  sum;
	if ( NumRows(A) != NumCols(A) )
		mes_error(E_SQUARE,"CHfactor");
	n = NumCols(A);
	for ( k=1; k<=n; k++ ){
		double norm=mes_norm_inf(A.Col(k));
		double eps=norm*Abs(mtl_numeric_limits<T>::epsilon());
		/* do diagonal element */
		for ( sum = A[k][k], j=1; j<k; j++ ) sum -= A[k][j]*conj(A[k][j]);
		if ( Abs(sum) < eps || Abs(imag(sum)) > eps ){
			mes_error(E_POSDEF,"CHfactor");
			return 0;
		}
		else{
			sum=(sum+conj(sum))/T(2); // remove imaginary part
		}
		A[k][k] = sqrt(sum);

		/* set values of column k */
		for ( i=k+1; i<=n; i++ ){
			for ( sum = A[i][k], j=1; j<k; j++ ) sum -= A[i][j]*conj(A[k][j]);
			A[i][j] = sum/A[k][k];
			A[j][i] = conj(A[i][j]);
		}
	}
	return 1;
}


/* CHsolve -- given a CHolesky factorization in A, solve A*x=b */
template <typename T> int	Matrix<T>::mes_CHsolve(const Matrix<T>& A,const Vector<T>& b,Vector<T>& x){
	if ( NumRows(A) != NumCols(A) || NumCols(A) != b.Size() )
		mes_error(E_SIZES,"CHsolve");
	x.Resize(b.Size());
	mes_Lsolve(A,b,x,0.0);
	mes_Usolve(A,x,x,0.0);

	return 1;
}


/* LDLfactor -- L.D.L' factorisation of A in-situ */
template <typename T> int	Matrix<T>::mes_LDLfactor(Matrix<T>& A){
	int	i, k, n, p;
	double d, sum;
	static Vector<T>	r;

	if ( NumRows(A) != NumCols(A) )
		mes_error(E_SQUARE,"LDLfactor");
	n = NumCols(A);
	r.Resize(n);
	for ( k = 1; k <= n; k++ ){
		sum = 0.0;
		for ( p = 1; p < k; p++ ){
		    r[p] = A[p][p]*A[k][p];
		    sum += r[p]*A[k][p];
		}
		d = A[k][k] -= sum;

		if ( d < mtl_numeric_limits<double>::epsilon() )
		    mes_error(E_SING,"LDLfactor");
		for ( i = k+1; i <= n; i++ ){
		    sum = 0.0;
		    for ( p = 1; p < k; p++ )
			sum += A[i][p]*r[p];
		    A[i][k] = (A[i][k] - sum)/d;
		}
	}

	return 1;
}


/* LDLsolve -- given an LDL factorisation in A, solve Ax=b */
template <typename T> int	Matrix<T>::mes_LDLsolve(const Matrix<T>& LDL,const Vector<T>& b,Vector<T>& x){
	if ( NumRows(LDL) != NumCols(LDL) )
		mes_error(E_SQUARE,"LDLsolve");
	if ( NumRows(LDL) != b.Size() )
		mes_error(E_SIZES,"LDLsolve");
	x.Resize(b.Size());
	mes_Lsolve(LDL,b,x,1.0);
	mes_Dsolve(LDL,x,x);
	mes_LTsolve(LDL,x,x,1.0);
	return 1;
}


/* LUfactor -- gaussian elimination with scaled partial pivoting
		-- Note: returns LU matrix which is A */
template <typename T> int	Matrix<T>::mes_LUfactor(Matrix<T>& A,Vector<int>& pivot){
	int	i, j, k, k_max, m, n;
	int	i_max;
	double	max1, temp, tiny;
	static	Vector<double> scale;
	T atemp;

	pivot.Resize(NumRows(A));
	if ( pivot.Size() != NumRows(A) )
		mes_error(E_SIZES,"LUfactor");
	m = NumRows(A);	n = NumCols(A);
	scale.Resize(NumRows(A));

	tiny = 10.0/HUGE_VAL;

	/* initialise pivot with identity permutation */
	for ( i=1; i<=m; i++ )
		pivot[i] = i;

	/* set scale parameters */
	for ( i=1; i<=m; i++ ){
		max1 = 0.0;
		for ( j=1; j<=n; j++ ){
			temp = Abs(A[i][j]);
			max1 = local_max(max1,temp);
		}
		scale[i] = max1;
	}

	/* main loop */
	k_max = local_min(m,n);
	for ( k=1; k<=k_max; k++ ){
	    /* find best pivot row */
	    max1 = 0.0;	i_max = 0;
	    for ( i=k; i<=m; i++ )
		if ( Abs(scale[i]) >= tiny*Abs(A[i][k]) ){
		    temp = Abs(A[i][k])/scale[i];
		    if ( temp > max1 ){
				max1 = temp;	i_max = i;
			}
		}
	    
	    /* if no pivot then ignore column k... */
	    if ( i_max == 0 ){
		/* set pivot entry A[k][k] exactly to zero,
			rather than just "small" */
			A[k][k] = 0.0;
			continue;
		}
		
		/* do we pivot ? */
		if ( i_max != k ){	/* yes we do... */
			mes_px_transp(pivot,i_max,k);
			for ( j=1; j<=n; j++ ){
				atemp = A[i_max][j];
				A[i_max][j] = A[k][j];
				A[k][j] = atemp;
			}
		}
		
		/* row operations */
		for ( i=k+1; i<=m; i++ ){	/* for each row do... */
			/* Note: divide by zero should never happen */
			atemp = A[i][k] = A[i][k]/A[k][k];
			if ( k+1 <= n ){
				for ( j=k+1; j<=n; j++ )
					A[i][j] -= atemp*A[k][j];
			}
		}
	    
	}
	return 1;
}


/* LUsolve -- given an LU factorisation in A, solve A*x=b */
template <typename T> int	Matrix<T>::mes_LUsolve(Matrix<T>& A,Vector<int>& pivot,const Vector<T>& b,Vector<T>& x){
	if ( NumRows(A) != NumCols(A) || NumCols(A) != b.Size() )
		mes_error(E_SIZES,"LUsolve");
	Vector<T> btmp=b;
	x.Resize(b.Size());
	mes_px_vec(pivot,btmp,x);	/* x := P.b */
	mes_Lsolve(A,x,x,1.0);	/* implicit diagonal = 1 */
	mes_Usolve(A,x,x,0.0);	/* explicit diagonal */
	return 1;
}


/* LUTsolve -- given an LU factorisation in A, solve A^T*x=b */
template <typename T> int	Matrix<T>::mes_LUTsolve(Matrix<T>& LU,Vector<int>& pivot,const Vector<T>& b,Vector<T>& x)
{
	if ( NumRows(LU) != NumCols(LU) || NumCols(LU) != b.Size() )
		mes_error(E_SIZES,"LUTsolve");
	x = b;
	mes_UTsolve(LU,x,x,0.0);	/* explicit diagonal */
	mes_LTsolve(LU,x,x,1.0);	/* implicit diagonal = 1 */
	mes_pxinv_vec(pivot,x,x);	/* x := P^T.tmp */
	return 1;
}


/* m_inverse -- returns inverse of A, provided A is not too rank deficient
	-- uses LU factorisation */
template <typename T> int	Matrix<T>::mes_inverse(const Matrix<T>& A,Matrix<T>& out)
{
	int	i;
	static Vector<T> tmp, tmp2;
	static Matrix<T> A_cp;
	static Vector<int> pivot;

	if ( NumRows(A) != NumCols(A) )
	    mes_error(E_SQUARE,"m_inverse");
	if ( NumRows(out) != NumRows(A) || NumCols(out) != NumCols(A) )
	    out.Resize(NumRows(A),NumCols(A));

	A_cp = A;
	tmp.Resize(NumRows(A));
	tmp2.Resize(NumRows(A));
	pivot.Resize(NumRows(A));
	mes_LUfactor(A_cp,pivot);
	for ( i = 1; i <= NumCols(A); i++ ){
	    tmp = T();
	    tmp[i] = 1.0;
	    mes_LUsolve(A_cp,pivot,tmp,tmp2);
	    mes_set_col(out,i,tmp2);
	}
	return 1;
}


/* LUcondest -- returns an estimate of the condition number of LU given the
	LU factorisation in compact form */
template <typename T> double	Matrix<T>::mes_LUcondest(Matrix<T>& LU,Vector<int>& pivot)
{
    static	Vector<T> y, z;
    double	cond_est, L_norm, U_norm, norm;
	T sum;
    int		i, j, n;

    if ( NumRows(LU) != NumCols(LU) )
	mes_error(E_SQUARE,"LUcondest");
    if ( NumCols(LU) != pivot.Size() )
	mes_error(E_SIZES,"LUcondest");

    n = NumCols(LU);
    y.Resize(n);
    z.Resize(n);
    for ( i = 1; i <= n; i++ ){
		sum = T();
		for ( j = 1; j < i; j++ )
			sum -= LU[j][i]*y[j];
		sum += mtl_sign(sum); // To support complex
		if ( mtl_iszero(LU[i][i]) )
			return HUGE_VAL;
		y[i] = sum / LU[i][i];
    }
	
    if(mes_LTsolve(LU,y,y,1.0) == 0)
		return HUGE_VAL;
	if(mes_LUsolve(LU,pivot,y,z) == 0)
		return HUGE_VAL;

    /* now estimate norm of A (even though it is not directly available) */
    /* actually computes ||L||_inf.||U||_inf */
    U_norm = 0.0;
    for ( i = 1; i <= n; i++ ){
		norm = 0.0;
		for ( j = i; j <= n; j++ )
			norm += Abs(LU[i][j]);
		if ( norm > U_norm )
			U_norm = norm;
    }
    L_norm = 0.0;
    for ( i = 1; i <= n; i++ ){
		norm = 1.0;
		for ( j = 1; j < i; j++ )
			norm += Abs(LU[i][j]);
		if ( norm > L_norm )
			L_norm = norm;
    }
	
    cond_est = U_norm*L_norm*mes_norm_inf(z)/mes_norm_inf(y);	
    return cond_est;
}


/* Note: The usual representation of a Householder transformation is taken
   to be:
   P = I - beta.u.uT
   where beta = 2/(uT.u) and u is called the Householder vector
*/

/* QRfactor -- forms the QR factorisation of A -- factorisation stored in
   compact form as described above ( not quite standard format ) */
template <typename T> int	Matrix<T>::mes_QRfactor(Matrix<T>& A,Vector<T>& diag)
{
    int	k,limit;
    double	beta;
    static	Vector<T> tmp1;
    
    limit = local_min(NumRows(A),NumCols(A));
	diag.Resize(limit);
    if ( diag.Size() < limit )
	mes_error(E_SIZES,"QRfactor");
    
    tmp1.Resize(NumRows(A));
    
    for ( k=1; k<=limit; k++ ){
		/* get H/holder vector for the k-th column */
		mes_get_col(A,k,tmp1);
		mes_hhvec(tmp1,k,&beta,tmp1,&A[k][k]);
		diag[k] = tmp1[k];
		
		/* apply H/holder vector to remaining columns */
		mes_hhtrcols(A,k,k+1,tmp1,beta);
    }

    return 1;
}


/* QRCPfactor -- forms the QR factorisation of A with column pivoting
   -- factorisation stored in compact form as described above
   ( not quite standard format )				*/
template <typename T> int	Matrix<T>::mes_QRCPfactor(Matrix<T>& A,Vector<T>& diag,Vector<int>& px)
{
    int	i, i_max, j, k, limit;
    static	Vector<T> gamma, tmp1, tmp2;
    double	beta, maxgamma, sum, tmp;
    
    limit = local_min(NumRows(A),NumCols(A));
	diag.Resize(limit);
	px.Resize(NumCols(A));
    if ( diag.Size() < limit || px.Size() != NumCols(A) )
	mes_error(E_SIZES,"QRCPfactor");
    
    tmp1.Resize(NumRows(A));
    tmp2.Resize(NumRows(A));
    gamma.Resize(NumCols(A));
    
    /* initialise gamma and px */
    for ( j=1; j<=NumCols(A); j++ ){
		px[j] = j;
		sum = 0.0;
		for ( i=1; i<=NumRows(A); i++ )
			sum += (A[i][j]*A[i][j]);
		gamma[j] = sum;
    }
    
    for ( k=1; k<=limit; k++ ){
		/* find "best" column to use */
		i_max = k;	maxgamma = gamma[k];
		for ( i=k+1; i<=NumCols(A); i++ ){
		/* Loop invariant:maxgamma=gamma[i_max]
	       >=gamma[l];l=k,...,i-1 */
		   if ( gamma[i] > maxgamma ){
			   maxgamma = gamma[i]; i_max = i;
		   }
		}
		   
		/* swap columns if necessary */
		if ( i_max != k ){
			/* swap gamma values */
			tmp = gamma[k];
			gamma[k] = gamma[i_max];
			gamma[i_max] = tmp;
			
			/* update column permutation */
			mes_px_transp(px,k,i_max);
			
			/* swap columns of A */
			for ( i=1; i<=NumRows(A); i++ ){
				tmp = A[i][k];
				A[i][k] = A[i][i_max];
				A[i][i_max] = tmp;
			}
		}
		
		/* get H/holder vector for the k-th column */
		mes_get_col(A,k,tmp1);
		mes_hhvec(tmp1,k,&beta,tmp1,&A[k][k]);
		diag[k] = tmp1[k];
		
		/* apply H/holder vector to remaining columns */
		mes_hhtrcols(A,k,k+1,tmp1,beta);
		
		/* update gamma values */
		for ( j=k+1; j<=NumCols(A); j++ )
			gamma[j] -= (A[k][j]*A[k][j]);
    }	
    return 1;
}


/* Qsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact
   form a la QRfactor() -- may be in-situ */
template <typename T> int	Matrix<T>::mes_Qsolve_(Matrix<T>& QR,Vector<T>& diag,const Vector<T>& b,Vector<T>& x,Vector<T>& tmp)
{
    int		k, limit;
    double	beta, r_ii, tmp_val;
    
    limit = local_min(NumRows(QR),NumCols(QR));
    if ( diag.Size() < limit || b.Size() != NumRows(QR) )
		mes_error(E_SIZES,"Qsolve_");
    x.Resize(NumRows(QR));
    tmp.Resize(NumRows(QR));
    
    /* apply H/holder transforms in normal order */
    x = b;
    for ( k = 1 ; k <= limit ; k++ ){
		mes_get_col(QR,k,tmp);
		r_ii = Abs(tmp[k]);
		tmp[k] = diag[k];
		tmp_val = (r_ii*Abs(diag[k]));
		beta = ( tmp_val == 0.0 ) ? 0.0 : 1.0/tmp_val;
		mes_hhtrvec(tmp,beta,k,x,x);