matrix2.txt

几个算法程序

武汉白云黄鹤站∶精华区武汉白云黄鹤站∶精华区
发信人: sweeping (simon), 信区: Algorithm WWW-POST 
标  题: 一些矩阵源码2 
发信站: 武汉白云黄鹤站 (Thu Nov  8 17:40:26 2001) , 转信 
 
 
/////////////////////////////// 
// CMatrix.cpp 
// please refer to the book by Xu Shiliang 
 
#include "stdafx.h" 
#include <math.h> 
 
#include "Matrix.h" 
 
////////////////////// 
// methods for CMatrix 
 
CMatrix::CMatrix(int rows, int cols) 
{ 
int i, n; 
 
if ( rows<1 || cols<1 || rows>20000 || cols>20000 ) 
{ 
m_a=NULL; 
m_rows=0; 
m_cols=0; 
return; 
} 
m_rows=rows; 
m_cols=cols; 
n=rows*cols; 
m_a=new double[n]; 
if ( m_a==NULL ) 
{ 
m_rows=0; 
m_cols=0; 
return; 
} 
 
for ( i=0;  i<n;  i++ ) 
m_a[i]=0.0; 
} 
 
CMatrix::~CMatrix() 
{ 
if ( m_a!=NULL ) 
delete []m_a; 
} 
 
BOOL CMatrix::Resize(int rows, int cols) 
{ 
int i, n; 
 
if ( rows<1 || cols<1 || rows>20000 || cols>20000 ) 
return FALSE; 
if ( m_a!=NULL ) 
{ 
delete []m_a; 
m_a=NULL; 
} 
m_rows=rows; 
m_cols=cols; 
m_a=new double[rows*cols]; 
if ( m_a==NULL ) 
{ 
m_rows=0; 
m_cols=0; 
return FALSE; 
} 
 
n=rows*cols; 
for ( i=0;  i<n;  i++ ) 
m_a[i]=0.0; 
return TRUE; 
} 
 
void CMatrix::SetElement(int row, int col, double e) 
{ 
if ( m_rows<1 || m_cols<1 ) 
return; 
 
if ( row<0 || row>=m_rows || col<0 || col>=m_cols ) 
return; 
m_a[row*m_cols+col]=e; 
} 
 
void CMatrix::SetAllElements(double e) 
{ 
int i, n; 
 
if ( m_rows<1 || m_cols<1 ) 
return; 
 
n=m_rows*m_cols; 
for ( i=0;  i<n;  i++ ) 
m_a[i]=e; 
} 
 
double CMatrix::GetElement(int row, int col) 
{ 
if ( m_rows<1 || m_cols<1 ) 
return 0.0; 
 
if ( row<0 || row>=m_rows || col<0 || col>=m_cols ) 
return 0.0; 
return m_a[row*m_cols+col]; 
} 
 
int CMatrix::GetRows() 
{ 
return m_rows; 
} 
 
int CMatrix::GetColumns() 
{ 
return m_cols; 
} 
 
double *CMatrix::GetBuffer() 
{ 
return m_a; 
} 
 
CMatrix *CMatrix::Copy() 
{ 
CMatrix *pC; 
int             row, col; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
 
pC=new CMatrix(m_rows, m_cols); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_rows;  row++ ) 
for ( col=0;  col<m_cols;  col++ ) 
pC->SetElement(row, col, GetElement(row, col)); 
return pC; 
} 
 
CMatrix *CMatrix::AddMatrix(CMatrix *pB) 
{ 
CMatrix *pC; 
int             row, col; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
if ( pB==NULL ) 
return NULL; 
if ( m_rows!=pB->GetRows() && m_cols!=pB->GetColumns() ) 
return NULL; 
 
pC=new CMatrix(m_rows, m_cols); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_rows;  row++ ) 
for ( col=0;  col<m_cols;  col++ ) 
pC->SetElement(row, col, GetElement(row, col)+pB->GetElement(row, col)); 
 
return pC; 
} 
 
CMatrix *CMatrix::SubMatrix(CMatrix *pB) 
{ 
CMatrix *pC; 
int             row, col; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
if ( pB==NULL ) 
return NULL; 
if ( m_rows!=pB->GetRows() && m_cols!=pB->GetColumns() ) 
return NULL; 
 
pC=new CMatrix(m_rows, m_cols); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_rows;  row++ ) 
for ( col=0;  col<m_cols;  col++ ) 
pC->SetElement(row, col, GetElement(row, col)-pB->GetElement(row, col)); 
 
return pC; 
} 
 
CMatrix *CMatrix::Transpose() 
{ 
CMatrix *pC; 
int             row, col; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
 
pC=new CMatrix(m_cols, m_rows); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_rows;  row++ ) 
for ( col=0;  col<m_cols;  col++ ) 
pC->SetElement(row, col, GetElement(col, row)); 
 
return pC; 
} 
 
CMatrix *CMatrix::Multiply(CMatrix *pB) 
{ 
CMatrix *pC; 
int             row, col; 
int             rows, cols; 
int             i; 
double          sum; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
if ( pB==NULL ) 
return NULL; 
rows=pB->GetRows(); 
cols=pB->GetColumns(); 
if ( m_cols!=rows ) 
return NULL; 
 
pC=new CMatrix(m_rows, cols); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_rows;  row++ ) 
for ( col=0;  col<cols;  col++ ) 
{ 
sum=0.0; 
for ( i=0;  i<m_cols;  i++ ) 
sum += GetElement(row, i)*pB->GetElement(i, col); 
pC->SetElement(row, col, sum); 
} 
 
return pC; 
} 
 
CSymmetricalMatrix *CMatrix::LeftTransposeMultiply() 
{ 
CSymmetricalMatrix *pC; 
int             row, col; 
int             i, j; 
double          sum; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
pC=new CSymmetricalMatrix(m_cols); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0;  row<m_cols;  row++ ) 
{ 
for ( col=0;  col<m_cols;  col++ ) 
{ 
sum=0.0; 
for ( i=0, j=0;  i<m_rows;  i++, j+=m_cols ) 
sum += m_a[j+row]*m_a[j+col]; 
pC->SetElement(row, col, sum); 
} 
} 
 
return pC; 
} 
 
CSymmetricalMatrix *CMatrix::RightTransposeMultiply() 
{ 
CSymmetricalMatrix *pC; 
int             row, col; 
int             i, j, k; 
double          sum; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
pC=new CSymmetricalMatrix(m_rows); 
if ( pC==NULL ) 
return NULL; 
 
for ( row=0, j=0;  row<m_rows;  row++, j+=m_cols ) 
{ 
for ( col=0, k=0;  col<m_rows;  col++, k+=m_cols ) 
{ 
sum=0.0; 
for ( i=0;  i<m_rows;  i++ ) 
sum += m_a[j+row]*m_a[k+col]; 
pC->SetElement(row, col, sum); 
} 
} 
 
return pC; 
} 
 
BOOL CMatrix::ReadFromFile(char *filename) 
{ 
int cols, rows, i, n; 
CFile RawFile; 
 
if ( !RawFile.Open(filename, CFile::modeRead) ) 
return FALSE; 
 
if ( RawFile.Read((void *)&rows, sizeof(int))!=sizeof(int) ) 
{ 
RawFile.Close(); 
return FALSE; 
} 
if ( RawFile.Read((void *)&cols, sizeof(int))!=sizeof(int) ) 
{ 
RawFile.Close(); 
return FALSE; 
} 
if ( rows<1 || cols<1 || rows>20000 || cols>20000 ) 
{ 
RawFile.Close(); 
return FALSE; 
} 
n=rows*cols; 
if ( n!=m_rows*m_cols ) 
{ 
if ( m_a!=NULL ) 
delete []m_a; 
m_a=new double[n]; 
if ( m_a==NULL ) 
{ 
m_rows=0; 
m_cols=0; 
RawFile.Close(); 
return FALSE; 
} 
} 
m_rows=rows; 
m_cols=cols; 
for ( i=0;  i<n;  i++ ) 
m_a[i]=0.0; 
for ( i=0;  i<n;  i++ ) 
{ 
if ( RawFile.Read((void *)&m_a[i], sizeof(double))!=sizeof(double) ) 
{ 
RawFile.Close(); 
return FALSE; 
} 
} 
RawFile.Close(); 
return TRUE; 
} 
 
BOOL CMatrix::WriteToFile(char *filename) 
{ 
int i, n; 
CFile RawFile(filename, CFile::modeCreate|CFile::modeWrite|CFile::modeRead); 
 
n=m_rows*m_cols; 
try 
{ 
RawFile.Write((void *)&m_rows, sizeof(int)); 
} 
catch(CFileException fe) 
{ 
RawFile.Close(); 
return FALSE; 
} 
try 
{ 
RawFile.Write((void *)&m_cols, sizeof(int)); 
} 
catch(CFileException fe) 
{ 
RawFile.Close(); 
return FALSE; 
} 
for ( i=0;  i<n;  i++ ) 
{ 
try 
{ 
RawFile.Write((void *)&m_a[i], sizeof(double)); 
} 
catch(CFileException fe) 
{ 
RawFile.Close(); 
return FALSE; 
} 
} 
return TRUE; 
} 
 
///////////////////// 
// methods for CSymmetricalMatrix 
 
// calculate the eigenvalue and eigenvictor of  
// nnxnn symmetric real matrix  
// the result eigenvictors are in returned matrix, 
// the result eigenvalues are in the diagonal of original matrix 
// original matrix is destroyed! 
// if original matrix is not symmetrical, 
// we will use the upper-right to fill out down-left to make it symmetrical 
CSymmetricalMatrix *CSymmetricalMatrix::Jacob() 
{ 
CSymmetricalMatrix *pC; 
int i, j; 
double a; 
 
if ( m_rows<1 || m_cols<1 ) 
return NULL; 
if ( m_rows!=m_cols ) 
return NULL; 
 
pC=new CSymmetricalMatrix(m_rows); 
if ( pC==NULL ) 
return NULL; 
 
for ( i=0;  i<m_rows-1;  i++ ) 
{ 
for ( j=i+1;  j<m_cols;  j++ ) 
{ 
a=GetElement(i, j); 
SetElement(j, i, a); 
} 
} 
 
jacobian(GetBuffer(), pC->GetBuffer(), m_rows); 
 
return pC; 
} 
 
// calculate the eigenvalue and eigenvector of  
// nnxnn symmetric real matrix aa 
// the result eigenvectors are ss, 
// the result eigenvalues are in the diagonal of aa 
void CSymmetricalMatrix::jacobian(double *aa, double *ss, int nn) 
{ 
int             i, j, k, l, n;