mes_eig.h

矩阵运算的模板类

					magdet = det_re*det_re+det_im*det_im;
				}
				invdet_re =   det_re/magdet;
				invdet_im = - det_im/magdet;
				tmp_val1_re = t22_re*val1_re-t22_im*val1_im-t12*val2_re;
				tmp_val1_im = t22_im*val1_re+t22_re*val1_im-t12*val2_im;
				tmp_val2_re = t11_re*val2_re-t11_im*val2_im-t21*val1_re;
				tmp_val2_im = t11_im*val2_re+t11_re*val2_im-t21*val1_im;
				tmp1_re[j-1] = invdet_re*tmp_val1_re - invdet_im*tmp_val1_im;
				tmp1_im[j-1] = invdet_im*tmp_val1_re + invdet_re*tmp_val1_im;
				tmp1_re[j]   = invdet_re*tmp_val2_re - invdet_im*tmp_val2_im;
				tmp1_im[j]   = invdet_im*tmp_val2_re + invdet_re*tmp_val2_im;
				j -= 2;
			}
			else{
				t11_re = TT[j][j] - l_re;
				t11_im = - l_im;
				magdet = t11_re*t11_re + t11_im*t11_im;
				scale = Abs(TT[j][j]) + Abs(l_re);
				if ( sqrt(magdet) < MACH_EPS*scale ){
					t11_re = MACH_EPS*scale;
					magdet = t11_re*t11_re + t11_im*t11_im;
				}
				invdet_re =   t11_re/magdet;
				invdet_im = - t11_im/magdet;
				for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_re[k]*TT[j][k]; 
				val1_re = tmp1_re[j] - sum;
				for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_im[k]*TT[j][k]; 
				val1_im = tmp1_im[j] - sum;
				tmp1_re[j] = invdet_re*val1_re - invdet_im*val1_im;
				tmp1_im[j] = invdet_im*val1_re + invdet_re*val1_im;
				j -= 1;
			}
		}
		
		norm = mes_norm_inf(tmp1_re) + mes_norm_inf(tmp1_im);
		tmp1_re *= (1/norm);
		if ( l_im != 0.0 )
			tmp1_im *= (1/norm);
		tmp2_re=QQ*tmp1_re;
		if ( l_im != 0.0 )
			tmp2_im=QQ*tmp1_im;
		if ( l_im != 0.0 )
			norm = pythagoras(mes_norm2(tmp2_re),mes_norm2(tmp2_im));
		else
			norm = mes_norm2(tmp2_re);
		tmp2_re *= (1/norm);
		if ( l_im != 0.0 )
			tmp2_im *= (1/norm);
		
		if ( l_im != 0.0 ){
			mes_set_col(XX_re,i,tmp2_re);
			mes_set_col(XX_im,i,tmp2_im);
			tmp2_im *= -1.0;
			mes_set_col(XX_re,i+1,tmp2_re);
			mes_set_col(XX_im,i+1,tmp2_im);
			i += 2;
		}
		else{
			mes_set_col(XX_re,i,tmp2_re);
			mes_set_col(XX_im,i,tmp1_im);	/* zero vector */
			i += 1;
		}
	}
		
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// Hessenberg factorisations
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// Hfactor -- compute Hessenberg factorisation in compact form.
//	-- factorisation performed in situ
//	-- for details of the compact form see QRfactor.c and matrix2.doc
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_Hfactor(Matrix<T>& A, Vector<T>& diag, Vector<T>& beta){
	static	Vector<double>	tmp1;
	int	k, limit;
	double beta1;

	if(A.IsDouble()==false)
		mes_warn(WARN_SINGLE_PRECISION,"mes_Hfactor");

	if ( diag.Size() < NumRows(A) - 1 || beta.Size() < NumRows(A) - 1 )
		mes_error(E_SIZES,"Hfactor");
	if ( NumRows(A) != NumCols(A) )
		mes_error(E_SQUARE,"Hfactor");
	limit = NumRows(A);
	tmp1.Resize(NumRows(A));
	for ( k = 1; k < limit; k++ ){
		mes_get_col(A,(int)k,tmp1);
		beta1 = beta[k];
		mes_hhvec(tmp1,k+1,&beta1,tmp1,&A[k+1][k]);
		beta[k] = beta1;
		diag[k] = tmp1[k+1];
		mes_hhtrcols(A,k+1,k+1,tmp1,beta1);
		mes_hhtrrows(A,1  ,k+1,tmp1,beta1);
	}
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// makeHQ -- construct the Hessenberg orthogonalising matrix Q;
//	-- i.e. Hess M = Q * M *Q^T
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_makeHQ(const Matrix<T>& H, Vector<T>& diag, Vector<T>& beta, Matrix<T>& Qout){
	int	i, j, limit;
	static	Vector<T>	tmp1, tmp2;

	limit = NumRows(H);
	if ( diag.Size() < limit || beta.Size() < limit )
		mes_error(E_SIZES,"makeHQ");
	if ( NumRows(H) != NumCols(H) )
		mes_error(E_SQUARE,"makeHQ");
	Qout.Resize(NumRows(H),NumRows(H));

	tmp1.Resize(NumRows(H));
	tmp2.Resize(NumRows(H));

	for ( i = 1; i <= NumRows(H); i++ ){
		/* tmp1 = i'th basis vector */
		for ( j = 1; j <= NumRows(H); j++ )
			tmp1[j] = 0.0;
		tmp1[i] = 1.0;

		/* apply H/h transforms in reverse order */
		for ( j = limit-1; j > 0; j-- ){
			mes_get_col(H,(int)j,tmp2);
			tmp2[j+1] = diag[j];
			mes_hhtrvec(tmp2,beta[j],j+1,tmp1,tmp1);
		}

		/* insert into Qout */
		mes_set_col(Qout,(int)i,tmp1);
	}

	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// makeH -- construct actual Hessenberg matrix
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_makeH(const Matrix<T>& H,Matrix<T>& Hout){
	int	i, j, limit;

	if ( NumRows(H) != NumCols(H) )
		mes_error(E_SQUARE,"makeH");
	Hout.Resize(NumRows(H),NumRows(H));
	Hout=H;

	limit = NumRows(H);
	for ( i = 2; i <= limit; i++ )
		for ( j = 1; j < i-1; j++ )
			Hout[i][j] = 0.0;

	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// Householder transformation
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// hhvec -- calulates Householder vector to eliminate all entries after the
//	i0 entry of the vector vec. It is returned as out. May be in-situ */
// Note:
//   Do not change beta to float type.
//   The stability of the Householder transformation is very sensitive
//   to the precision of beta in hhvec() and hhtrrows().
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_hhvec(Vector<T>& vec,int i0,double* beta,Vector<T>& out,T* newval){
	double	norm;

	mes_copy_offset(vec,out,i0);
	norm = mes_norm2_offset(out,i0);
	if ( norm == 0.0 ){
		*beta = 0.0;
		*newval = out[i0];
		return 1;
	}
	double abs_val=Abs(out[i0]);
	*beta = 1.0/(norm * (norm+abs_val));
	if(abs_val == 0.0){
		*newval = norm;
	}
	else{
		*newval = -mtl_sign(out[i0])*T(norm);
	}
	out[i0] -= *newval;

	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// hhtrvec -- apply Householder transformation to vector -- may be in-situ
// hh = Householder vector
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_hhtrvec(Vector<T>& hh,double beta,int i0,Vector<T>& in,Vector<T>& out){
	T	scale;
	int	i;
	if ( in.Size() != hh.Size() )
		mes_error(E_SIZES,"hhtrvec");
	if ( i0 > in.Size() )
		mes_error(E_BOUNDS,"hhtrvec");
	scale = -T(beta)*mes_in_prod_offset(hh,in,i0);
	out = in;
	for ( i=i0; i<=out.Size(); i++ )
		out[i] += scale*hh[i];
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// hhtrrows -- transform a matrix by a Householder vector by rows
//	starting at row i0 from column j0 -- in-situ
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_hhtrrows(Matrix<T>& M,int i0,int j0,Vector<T>& hh,double beta){
	T	ip, scale;
	int	i;

	if ( NumCols(M) != hh.Size() )
		mes_error(E_RANGE,"hhtrrows");
	if ( i0 > NumRows(M)+1 || j0 > NumCols(M)+1 )
		mes_error(E_BOUNDS,"hhtrrows");

	if ( beta == 0.0 )	return 1;

	/* for each row ... */
	for ( i = i0; i <= NumRows(M); i++ ){	/* compute inner product */
		int j;
		ip = T();
		for ( j = j0; j <= NumCols(M); j++ )
			ip += M[i][j]*hh[j];
		scale = -beta*ip;
		if ( mtl_iszero(scale) )
		    continue;

		/* do operation */
		for ( j = j0; j <= NumCols(M); j++ )
			M[i][j] += scale*conj(hh[j]);
	}

	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// hhtrcols -- transform a matrix by a Householder vector by columns
//	starting at row i0 from column j0 -- in-situ
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_hhtrcols(Matrix<T>& M,int i0,int j0,Vector<T>& hh,double beta){
	int	i,k;
	T scale;
	static	Vector<T>	w;

	if ( NumRows(M) != hh.Size() )
		mes_error(E_SIZES,"hhtrcols");
	if ( i0 > NumRows(M)+1 || j0 > NumCols(M)+1 )
		mes_error(E_BOUNDS,"hhtrcols");

	if ( beta == 0.0 )	return 1;

	w.Resize(NumCols(M));
	w=T();

	for ( i = i0; i <= NumRows(M); i++ ){
		if ( !mtl_iszero(hh[i]) ){
			for(k=j0;k<=NumCols(M);k++){
				w[k]+=conj(M[i][k])*hh[i];
			}
		}		
	}
	for ( i = i0; i <= NumRows(M); i++ ){
		if ( !mtl_iszero(hh[i]) ){
			for(k=j0;k<=NumCols(M);k++){
				scale = -T(beta)*hh[i];
				M[i][k]+=conj(w[k])*scale;
			}
		}
	}
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// Givens rotations
//
// Complex Givens rotation matrix:
//   [ c   -s ]
//   [ s*   c ]
// Note that c is real and s is complex
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////
// givens -- returns c,s parameters for Givens rotation to
//		eliminate y in the vector [ x y ]'
///////////////////////////////////////////////////////////////////////////////
template <typename T> void	Matrix<T>::mes_givens(T x,T y,double* c,T* s){
	double norm = pythagoras(x,y);
	double inv_norm;
	if ( norm == 0.0 ){ /* identity */
		*c = 1.0;
		*s = 0.0;
	}	
	else{ // To support complex
		double absx=Abs(x);
		T xnorm;
		if(absx==0){
			xnorm=1;
		}
		else{
			xnorm = x/absx; // normalize x = x/|x|
		}
		inv_norm = 1/norm; // inv_norm = 1/||[x,y]||
		*c = inv_norm*absx;
		*s = -inv_norm*((xnorm)*conj(y)); // For real signals: c=x/norm, s=-y/norm
	}
}


///////////////////////////////////////////////////////////////////////////////
// rot_vec -- apply Givens rotation to x's i & k components
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_rot_vec(Vector<T>& x,int i,int k,double c,T s,Vector<T>& out){
	if ( i > x.Size() || k > x.Size() )
		mes_error(E_RANGE,"rot_vec");
	if(&x[1] != &out[1])
		out = x;
	T temp=c*out[i] - s*out[k]; // To support complex.
	out[k] = c*out[k] + conj(s)*out[i];
	out[i]=temp;
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// rot_rows -- premultiply mat by givens rotation described by c,s
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_rot_rows(Matrix<T>& mat,int i,int k,double c,T s,Matrix<T>& out){
	int	j;
	if ( i > NumRows(mat) || k > NumRows(mat) )
		mes_error(E_RANGE,"rot_rows");
	out = mat;
	for ( j=1; j<=NumCols(mat); j++ ){
		T temp=c*out[i][j] - s*out[k][j]; // To support complex.
		out[k][j] = c*out[k][j] + conj(s)*out[i][j];
		out[i][j] = temp;
	}

	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// rot_cols -- postmultiply mat by givens rotation described by c,s
///////////////////////////////////////////////////////////////////////////////
template <typename T> int	Matrix<T>::mes_rot_cols(Matrix<T>& mat,int i,int k,double c,T s,Matrix<T>& out){
	int	j;
	if ( i > NumCols(mat) || k > NumCols(mat) )
		mes_error(E_RANGE,"rot_cols");
	out = mat;
	for ( j=1; j<=NumRows(mat); j++ ){
		T temp=c*out[j][i] - conj(s)*out[j][k]; // To support complex.
		out[j][k] = c*out[j][k] + s*out[j][i];
		out[j][i] = temp;
	}
	return 1;
}


#endif