mes_svd.h

在数学计算中经常需要用到矩阵运算

			}
			else{	
				/* use the smallest eigenvalue of the matrix */
				double ztmp=0.5*(real(t11)-real(t22));
				double diff=sqrt(ztmp*ztmp+real(t12*t21));
				double sum=0.5*(real(t11)+real(t22));
				double det = real(t11)*real(t22)-real(t12*t21);
				if(Abs(sum+diff) > Abs(sum-diff))
					shift=det/(sum+diff);
				else
					shift=det/(sum-diff);
			}

			/* perturb shift if convergence is slow */
			if(iter%10==0){
				shift += iter*0.02;
			}
			
			/* initial Givens' rotation */
			mes_givens(conj(dd[i_min])*dd[i_min]-shift, conj(dd[i_min])*ff[i_min], &c, &s);
			
			/* do initial Givens' rotations */
			dtmp        = c*dd[i_min] - s*ff[i_min];
			ff[i_min]   = c*ff[i_min] + conj(s)*dd[i_min];
			dd[i_min]   = dtmp;
			z           = -s*dd[i_min+1];
			dd[i_min+1] =  c*dd[i_min+1];
			if ( NumCols(V) )
				mes_rot_rows(V,i_min,i_min+1,c,conj(s),V);
			
			/* 2nd Givens' rotation */
			mes_givens(dd[i_min],z, &c, &s);
			dd[i_min]   = c*dd[i_min] - s*z;
			dtmp        = c*dd[i_min+1] + conj(s)*ff[i_min];
			ff[i_min]   = c*ff[i_min] - s*dd[i_min+1];
			dd[i_min+1] = dtmp;
			if ( i_min+1 < i_max ){
				z           = -s*ff[i_min+1];
				ff[i_min+1] =  c*ff[i_min+1];
			}
			if ( NumCols(U) )
				mes_rot_rows(U,i_min,i_min+1,c,s,U);
			
			for ( i = i_min+1; i < i_max; i++ ){
				/* get Givens' rotation for zeroing z */
				mes_givens(ff[i-1],z, &c, &s);
				ff[i-1] = c*ff[i-1] - s*z;
				dtmp    = c*dd[i] - s*ff[i];
				ff[i]   = c*ff[i] + conj(s)*dd[i];
				dd[i]   = dtmp;
				z       = -s*dd[i+1];
				dd[i+1] =  c*dd[i+1];
				if ( NumCols(V) )
					mes_rot_rows(V,i,i+1,c,conj(s),V);
				
				/* get 2nd Givens' rotation */
				mes_givens(dd[i],z, &c, &s);
				dd[i]   = c*dd[i] - s*z;
				dtmp    = c*dd[i+1] + conj(s)*ff[i];
				ff[i]   = c*ff[i] - s*dd[i+1];
				dd[i+1] = dtmp;
				if ( i+1 < i_max ){
					z       = -s*ff[i+1];
					ff[i+1] =  c*ff[i+1];
				}
				if ( NumCols(U) )
					mes_rot_rows(U,i,i+1,c,s,U);
			}
			
			/* should matrix be split? */
			for ( i = i_min; i < i_max; i++ ){
				if ( Abs(ff[i]) < (MACH_EPS*(Abs(dd[i])+Abs(dd[i+1]))) ){
					split = 1;
					ff[i] = 0.0;
				}
				else if ( Abs(dd[i]) < (MACH_EPS*size) ){
					split = 1;
					dd[i] = 0.0;
				}
			}
			
			if(iter==maxIter){
				cerr << "ERROR: No convergence in " << maxIter << " mes_bisvd iterations" << endl;
				isOK=0;
				assert(0);
				return isOK;
			}
		}
	}

	/* rotate dd[i] so it is real and rotate U & V appropriately */
	isOK=mes_rotsvd(dd,U,V);

	/* sort singular values and singular vectors in the descending order of singular values */
	isOK=mes_fixsvd(dd,U,V);
	
	return isOK;
}


template <typename T> int Matrix<T>::mes_rotsvd(Vector<T>& dd, Matrix<T>& U, Matrix<T>& V)
{
	int i,j;
	double size = mes_norm_inf(dd);
	for(i=1; i<=dd.Size();i++){
		if((Abs(dd[i])>MACH_EPS*size) && (Abs(imag(dd[i]))>MACH_EPS*size)){
			double theta=atan2( Abs(imag(dd[i])), Abs(real(dd[i])) );
			T dtmp1=T(mtl_sign(real(dd[i]))*cos(theta/2), -mtl_sign(imag(dd[i]))*sin(theta/2));
			T dtmp2=-dtmp1;
			dd[i]=dd[i]*T(dtmp1*dtmp2);
			for(j=1;j<=NumCols(U);j++){
				U[i][j]=U[i][j]*dtmp1;
			}
			for(j=1;j<=NumCols(V);j++){
				V[i][j]=V[i][j]*conj(dtmp2);
			}
		}
		if(Abs(imag(dd[i]))<MACH_EPS*size){
			dd[i]=T(real(dd[i]),0);
		}
		else{
			cerr << "ERROR: singular values must be real.\n";
			assert(0);
		}
	}
	return 1;
}


inline int Matrix<float>::mes_rotsvd(Vector<float>& dd, Matrix<float>& U, Matrix<float>& V)
{
	return 1;
}


inline int Matrix<double>::mes_rotsvd(Vector<double>& dd, Matrix<double>& U, Matrix<double>& V)
{
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// bifactor -- perform preliminary factorisation for bisvd
//	-- updates U and/or V, which ever is defined
// Note:
//   Do not change beta to single precision (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_bifactor(Matrix<T>& A, Matrix<T>& U, Matrix<T>& V)
{
	int compU=1, compV=1;
	int	i, k, m, n;
	static Vector<T> tmp1, tmp2;
	double	beta;

	m=NumRows(A);
	n=NumCols(A);
	if ( ( NumRows(U) != NumCols(U) ) || ( NumRows(V) != NumCols(V) ) )
		mes_error(E_SQUARE,"bifactor");
	if ( ( NumRows(U) != m ) || ( NumRows(V) != n ) )
		mes_error(E_SIZES,"bifactor");
	tmp1.Resize(m);
	tmp2.Resize(n);
	
	if(m>=n){
		for ( k = 1; k <= n; k++ ){
			mes_get_col(A,k,tmp1);
			mes_hhvec(tmp1,k,&beta,tmp1,&A[k][k]);
			mes_hhtrcols(A,k,k+1,tmp1,beta);
			if ( NumCols(U) )
				mes_hhtrcols(U,k,1,tmp1,beta);
			if ( k >= n-1 )
				continue;
			mes_get_row(A,k,tmp2);
			mes_hhvec(tmp2,k+1,&beta,tmp2,&A[k][k+1]);
			for(i=1;i<=tmp2.Size();i++) // conujugate householder vector for V
				tmp2[i]=conj(tmp2[i]);
			mes_hhtrrows(A,k+1,k+1,tmp2,beta);
			if ( NumCols(V) )
				mes_hhtrcols(V,k+1,1,tmp2,beta);
		}
	}
	else{
		for ( k = 1; k <= m; k++ ){
			mes_get_row(A,k,tmp1);
			for(i=1;i<=tmp1.Size();i++) // conujugate column vector of A
				tmp1[i]=conj(tmp1[i]);
			mes_hhvec(tmp1,k,&beta,tmp1,&A[k][k]);
			mes_hhtrrows(A,k+1,k,tmp1,beta);
			if ( NumCols(V) )
				mes_hhtrcols(V,k,1,tmp1,beta);
			if ( k >= m-1 )
				continue;
			mes_get_col(A,k,tmp2);
			mes_hhvec(tmp2,k+1,&beta,tmp2,&A[k+1][k]);
			mes_hhtrcols(A,k+1,k+1,tmp2,beta);
			if ( NumCols(U) )
				mes_hhtrcols(U,k+1,1,tmp2,beta);
		}
	}
	return 1;
}


///////////////////////////////////////////////////////////////////////////////
// svd -- returns vector of singular values in d
//	-- also updates U and/or V, if one or the other is defined
//	-- destroys A
// Note:
//    Singular values for MxN (M<N) can be computed from SVD of A'
//    by conjugating A and interchanging the roles of U and V.
//      A = U' * S * V
//      A' = V' * S * U
///////////////////////////////////////////////////////////////////////////////
template <typename T> int Matrix<T>::mes_svd(const Matrix<T>& A, Matrix<T>& U, Matrix<T>& V, Vector<T>& dd)
{
	static Vector<T>	ff;
	int	i, limit, tryX, isOK;
	
#if 0
	if(NumRows(A)<NumCols(A)){ // This implementation wastes memory.
		return mes_svd(Transpose(A),V,U,dd);
	}
#endif
	if(A.IsDouble()==false)
		mes_warn(WARN_SINGLE_PRECISION,"mes_svd");
	int m = NumRows(A);
	int n = NumCols(A);
	U.Resize(m,m);
	V.Resize(n,n);
	if ( (NumRows(U) != NumCols(U)) || 
		   (NumRows(V) != NumCols(V)) )
	   mes_error(E_SQUARE,"svd");
	if ( (NumRows(U) != NumRows(A) ) || ( NumRows(V) != NumCols(A) ) )
		mes_error(E_SIZES,"svd");

	int maxIter=100;
	Matrix<T> A_tmp;
	limit = local_min(m,n);
	dd.Resize(limit);
	ff.Resize(limit-1);
	for(tryX=1;tryX<=2;tryX++){
		A_tmp=A;
		U.SetIdentity();
		V.SetIdentity();		
		isOK=mes_bifactor(A_tmp,U,V);
		if(isOK==0)
			break;
		if(m>=n){
			for ( i = 1; i <= limit; i++ ){
				dd[i] = (A_tmp[i][i]);
				if ( i+1 <= limit )
					ff[i] = (A_tmp[i][i+1]);
			}
		}
		else{
			for ( i = 1; i <= limit; i++ ){
				dd[i] = (A_tmp[i][i]);
				if ( i+1 <= limit )
					ff[i] = conj(A_tmp[i+1][i]);
			}
		}
		
		if(tryX==2) maxIter=5*maxIter;
		if(m>=n)
			isOK=mes_bisvd(dd,ff,U,V,maxIter);
		else
			isOK=mes_bisvd(dd,ff,V,U,maxIter);

		if(isOK != 0) // Succeeded!
			break;
	}
	if(isOK == 0){
		cerr << "ERROR: No convergence in " << maxIter << " bisvd iterations" << endl;
		assert(0);
	}
	return isOK;
}


#endif