技术文章 > 数值分析-Jacobi旋转矩阵

数值分析-Jacobi旋转矩阵

2018-08-22 01:40

文档管理软件,文档管理系统,知识管理系统,档案管理系统的技术资料:
/*****************************************************
* 文件: JacobiRotate.h
* 描述: Jacobi旋转矩阵
* 创建: ZBM 2003.12.17
*****************************************************/
//---------------------------------------------------------------------------
#ifndef JacobiRotateH
#define JacobiRotateH
//---------------------------------------------------------------------------
class CJacobiRotateMatrix
{
private:
unsigned int m_iDimension; // 矩阵维数

double **m_ppdOrignalMatrix; // 原始矩阵
double **m_ppdFeatureVectors; // 特征向量
double **m_ppdAn; // 中间结果阵
double **m_ppdJacobiMatrixN; // Jacobi阵
unsigned int m_iRowNo;
unsigned int m_iColNo;
unsigned int m_iMaxCount;
double m_dEposina;
private:
void GetTargetElemCoord( double **ppdMA );
double GetTheta( double **ppdM );
void InitDblMatrix( double **ppdM );
double **NewDblMatrix();
void DeleteDblMatrix( double **ppdM );
void CopyMatrix( double **pdDest, const double **ppdSrc );
void MatrixMultiply( const double **ppdSrcM1, const double **ppdSrcM2, double **ppdResult );
void MatrixTranslation( double **ppdM );
bool IsCompleted( double **ppdM );
void PrintMatrix( char *sName, double **ppdM );

public:
CJacobiRotateMatrix( unsigned int iDimension, double dEposina = 1e-6, unsigned int iMaxCount = 100 );
~CJacobiRotateMatrix();
void LoadMatrix( const double *pdMatrix ); // 加载矩阵
void UnloadResult( double **ppdMatrix ); // 导出结果矩阵
bool CheckValid();
void Calculate( unsigned int *piIterateCount );
};
#endif

/*****************************************************
* 文件: JacobiRotate.cpp
* 描述: Jacobi旋转矩阵
* 创建: ZBM 2003.12.17
*****************************************************/
//---------------------------------------------------------------------------
#pragma hdrstop
#include
#include
#include
#include
#include "JacobiRotate.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
// 定义浮点的0
#define DBL_ZERO (1e-5)
#define PI (3.14159265)
//---------------------------------------------------------------------------
// 构造
CJacobiRotateMatrix::CJacobiRotateMatrix( unsigned int iDimension, double dEposina, unsigned int iMaxCount )
{
m_iMaxCount = iMaxCount;
m_dEposina = dEposina;
m_iDimension = iDimension;
m_ppdOrignalMatrix = ( m_iDimension > 0 ) ? NewDblMatrix() : NULL;
m_ppdJacobiMatrixN = ( m_iDimension > 0 ) ? NewDblMatrix() : NULL;
m_ppdAn = ( m_iDimension > 0 ) ? NewDblMatrix() : NULL;
m_ppdFeatureVectors = ( m_iDimension > 0 ) ? NewDblMatrix() : NULL;
m_iRowNo = 0;
m_iColNo = 0;
}
//---------------------------------------------------------------------------
// 析构
CJacobiRotateMatrix::~CJacobiRotateMatrix()
{
if( m_iDimension > 0 )
{
DeleteDblMatrix( m_ppdOrignalMatrix );
DeleteDblMatrix( m_ppdJacobiMatrixN );
DeleteDblMatrix( m_ppdAn );
DeleteDblMatrix( m_ppdFeatureVectors );
}
getch();
}
//---------------------------------------------------------------------------
// 初始化矩阵元素(置0)
void CJacobiRotateMatrix::InitDblMatrix( double **ppdM )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
ppdM[ m ][ n ] = 0.0;
}
}
}
//---------------------------------------------------------------------------
// 分配新矩阵空间
double **CJacobiRotateMatrix::NewDblMatrix()
{
double **ppdMNew = NULL;
if( m_iDimension > 0 )
{
ppdMNew = new double*[ m_iDimension ];
for( unsigned int i = 0; i < m_iDimension; i++ )
{
ppdMNew[ i ] = new double[ m_iDimension ];
}
}
InitDblMatrix( ppdMNew ); // 初始化
return ppdMNew;
}
//---------------------------------------------------------------------------
// 释放矩阵空间
void CJacobiRotateMatrix::DeleteDblMatrix( double **ppdM )
{
if( m_iDimension > 0 )
{
for( unsigned int i = 0; i < m_iDimension; i++ )
{
delete[] ppdM[ i ];
}
delete[] ppdM;
}
}
//---------------------------------------------------------------------------
// 装载矩阵
void CJacobiRotateMatrix::LoadMatrix( const double *pdMatrix )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
m_ppdOrignalMatrix[ m ][ n ] = *(pdMatrix+ m*m_iDimension + n);
}
}
}
//---------------------------------------------------------------------------
// 导出结果矩阵
void CJacobiRotateMatrix::UnloadResult( double **ppdMatrix )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
ppdMatrix[ m ][ n ] = m_ppdFeatureVectors[ m ][ n ];
}
}
}
//---------------------------------------------------------------------------
// 检查矩阵是否适用(是否是对称矩阵)
bool CJacobiRotateMatrix::CheckValid()
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
if( fabs( m_ppdOrignalMatrix[ m ][ m ] ) < DBL_ZERO )
{
return false;
}
for( unsigned int n = m + 1; n < m_iDimension; n++ )
{
if( fabs( m_ppdOrignalMatrix[ m ][ n ] - m_ppdOrignalMatrix[ n ][ m ] ) > DBL_ZERO )
{
return false;
}
}
}
return true;
}
//---------------------------------------------------------------------------
// 取得目标元位置和值
void CJacobiRotateMatrix::GetTargetElemCoord( double **ppdMA )
{
double dMaxVal = 0.0;
m_iRowNo = 0;
m_iRowNo = 0;

for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = m + 1; n < m_iDimension; n++ )
{
if( fabs( ppdMA[ m ][ n ] ) > dMaxVal &&
fabs( fabs( ppdMA[ m ][ n ] ) - dMaxVal ) > DBL_ZERO )
{
dMaxVal = fabs( ppdMA[ m ][ n ] );
m_iRowNo = m;
m_iColNo = n;
}
}
}
}
//---------------------------------------------------------------------------
// 计算Theta角度
double CJacobiRotateMatrix::GetTheta( double **ppdM )
{
double dTmp = ppdM[ m_iRowNo ][ m_iRowNo ] - ppdM[ m_iColNo ][ m_iColNo ];
if( fabs( dTmp ) > DBL_ZERO )
{
dTmp = 2.0*ppdM[ m_iRowNo ][ m_iColNo ] / dTmp;
dTmp = atan( dTmp )/2.0;
}
else
{
dTmp = PI/4.0;
}
return dTmp;
}
//---------------------------------------------------------------------------
// 复制矩阵
void CJacobiRotateMatrix::CopyMatrix( double **pdDest, const double **ppdSrc )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
pdDest[ m ][ n ] = ppdSrc[ m ][ n ];
}
}
}
//---------------------------------------------------------------------------
// 矩阵相乘
void CJacobiRotateMatrix::MatrixMultiply( const double **ppdSrcM1, const double **ppdSrcM2, double **ppdResult )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
ppdResult[ m ][ n ] = 0.0;
for( unsigned int i = 0; i < m_iDimension; i++ )
{
ppdResult[ m ][ n ] += ppdSrcM1[ m ][ i ]*ppdSrcM2[ i ][ n ];
}
if( fabs( ppdResult[ m ][ n ] ) < DBL_ZERO )
{
ppdResult[ m ][ n ] = 0.0;
}
}
}
}
//---------------------------------------------------------------------------
// 矩阵转置
void CJacobiRotateMatrix::MatrixTranslation( double **ppdM )
{
double dTmp = 0.0;

for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = m + 1; n < m_iDimension; n++ )
{
dTmp = ppdM[ m ][ n ];
ppdM[ m ][ n ] = ppdM[ n ][ m ];
ppdM[ n ][ m ] = dTmp;
}
}
}
//---------------------------------------------------------------------------
// 判断是否计算完成
bool CJacobiRotateMatrix::IsCompleted( double **ppdM )
{
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
if( m != n )
{
if( fabs( ppdM[ m ][ n ] ) > m_dEposina )
{
return false;
}
}
}
}
return true;
}
//---------------------------------------------------------------------------
// 在屏幕上打印出矩阵
void CJacobiRotateMatrix::PrintMatrix( char *sName, double **ppdM )
{
printf( "%s:\n", sName );
for( unsigned int m = 0; m < m_iDimension; m++ )
{
for( unsigned int n = 0; n < m_iDimension; n++ )
{
printf( "%4.7lf", ppdM[ m ][ n ] );
if( n < m_iDimension - 1 )
{
printf(", ");
}
}
printf( "\n" );
}
}
//---------------------------------------------------------------------------
// 计算
void CJacobiRotateMatrix::Calculate( unsigned int *piIterateCount )
{
double **ppdMAnTmp = NewDblMatrix();
double **ppdMJnTmp = NewDblMatrix();
CopyMatrix( m_ppdAn, (const double **)m_ppdOrignalMatrix );
for( ( *piIterateCount ) = 0; ( *piIterateCount ) < m_iMaxCount; ( *piIterateCount )++ )
{
printf( "------------------------ %d --------------------------\n", ( *piIterateCount )+1 );
// 计算Jn
GetTargetElemCoord( m_ppdAn );
double dTheta = GetTheta( m_ppdAn );
printf( "i= %d, j= %d, theta= %3.2lf degree, arc = %lf \n", m_iRowNo, m_iColNo, dTheta*180.0/PI, dTheta );

InitDblMatrix( m_ppdJacobiMatrixN );
for( unsigned int i = 0; i < m_iDimension; i++ )
{
if( i != m_iRowNo && i != m_iColNo )
{
m_ppdJacobiMatrixN[ i ][ i ] = 1.0;
}
}
m_ppdJacobiMatrixN[ m_iRowNo ][ m_iRowNo ] = cos( dTheta );
m_ppdJacobiMatrixN[ m_iColNo ][ m_iColNo ] = m_ppdJacobiMatrixN[ m_iRowNo ][ m_iRowNo ];
m_ppdJacobiMatrixN[ m_iRowNo ][ m_iColNo ] = sin( dTheta );
m_ppdJacobiMatrixN[ m_iColNo ][ m_iRowNo ] = -m_ppdJacobiMatrixN[ m_iRowNo ][ m_iColNo ];
PrintMatrix( "Jn", m_ppdJacobiMatrixN );

// 计算An
// AnTmp = Jn*An;
MatrixMultiply( (const double **)m_ppdJacobiMatrixN, (const double **)m_ppdAn, ppdMAnTmp );
//PrintMatrix( "An_Tmp", ppdMAnTmp );
// Jn = (Jn)`
MatrixTranslation( m_ppdJacobiMatrixN );
//PrintMatrix( "Jn`", m_ppdJacobiMatrixN );

// m_ppdAn = AnTmp*Jn
MatrixMultiply( (const double **)ppdMAnTmp, (const double **)m_ppdJacobiMatrixN, m_ppdAn );
// Qn
if( ( *piIterateCount ) > 0 )
{
MatrixMultiply( (const double **)m_ppdFeatureVectors, (const double **)m_ppdJacobiMatrixN, ppdMJnTmp );
CopyMatrix( m_ppdFeatureVectors, (const double **)ppdMJnTmp );
}
else
{
CopyMatrix( m_ppdFeatureVectors, (const double **)m_ppdJacobiMatrixN );
}
PrintMatrix( "An", m_ppdAn );

if( IsCompleted( m_ppdAn ) )
{
printf( "========== completed ! ==============\n" );
PrintMatrix( "Q", m_ppdFeatureVectors );
getch();
break;
}
getch();
}
DeleteDblMatrix( ppdMAnTmp );
DeleteDblMatrix( ppdMJnTmp );
}
//---------------------------------------------------------------------------

#include
#include "JacobiRotate.h"
int main(int argc, char* argv[])
{const unsigned int kiDimension = 3;
double dM[ kiDimension ][ kiDimension ];
/*
dM[ 0 ][ 0 ] = 2.0;
dM[ 0 ][ 1 ] = -1.0;
dM[ 0 ][ 2 ] = 0.0;
dM[ 1 ][ 0 ] = -1.0;
dM[ 1 ][ 1 ] = 2.0;
dM[ 1 ][ 2 ] = -1.0;
dM[ 2 ][ 0 ] = 0.0;
dM[ 2 ][ 1 ] = -1.0;
dM[ 2 ][ 2 ] = 2.0;
*/
dM[ 0 ][ 0 ] = 2.0;
dM[ 0 ][ 1 ] = 0.0;
dM[ 0 ][ 2 ] = 0.0;
dM[ 1 ][ 0 ] = 0.0;
dM[ 1 ][ 1 ] = 3.0;
dM[ 1 ][ 2 ] = 1.0;
dM[ 2 ][ 0 ] = 0.0;
dM[ 2 ][ 1 ] = 1.0;
dM[ 2 ][ 2 ] = 4.0;

unsigned int iCount = 0;

CJacobiRotateMatrix JacobiM( kiDimension );
JacobiM.LoadMatrix( (const double *)dM );
JacobiM.CheckValid();
JacobiM.Calculate( &iCount );
return 0;

}