操作矩阵的类 cmatrix.txt

一个操作矩阵的类CMatrix的算法

} 

////////////////////////////////////////////////////////////////////// 
// 重载运算符-，实现矩阵的减法 
// 
// 参数： 
// 1. const CMatrix& other - 与指定矩阵相减的矩阵 
// 
// 返回值：CMatrix型，指定矩阵与other相减之差 
////////////////////////////////////////////////////////////////////// 
CMatrix CMatrix::operator-(const CMatrix& other) const 
{ 
// 首先检查行列数是否相等 
ASSERT (m_nNumColumns == other.GetNumColumns() && m_nNumRows == 
other.GetNumRows()); 

// 构造目标矩阵 
CMatrix result(*this) ; // copy ourselves 
// 进行减法操作 
for (int i = 0 ; i < m_nNumRows ; ++i) 
{ 
for (int j = 0 ; j < m_nNumColumns; ++j) 
result.SetElement(i, j, result.GetElement(i, j) - 
other.GetElement(i, j)) ; 
} 

return result ; 
} 

////////////////////////////////////////////////////////////////////// 
// 重载运算符*，实现矩阵的数乘 
// 
// 参数： 
// 1. double value - 与指定矩阵相乘的实数 
// 
// 返回值：CMatrix型，指定矩阵与value相乘之积 
////////////////////////////////////////////////////////////////////// 
CMatrix CMatrix::operator*(double value) const 
{ 
// 构造目标矩阵 
CMatrix result(*this) ; // copy ourselves 
// 进行数乘 
for (int i = 0 ; i < m_nNumRows ; ++i) 
{ 
for (int j = 0 ; j < m_nNumColumns; ++j) 
result.SetElement(i, j, result.GetElement(i, j) * 
value) ; 
} 

return result ; 
} 

////////////////////////////////////////////////////////////////////// 
// 重载运算符*，实现矩阵的乘法 
// 
// 参数： 
// 1. const CMatrix& other - 与指定矩阵相乘的矩阵 
// 
// 返回值：CMatrix型，指定矩阵与other相乘之积 
////////////////////////////////////////////////////////////////////// 
CMatrix CMatrix::operator*(const CMatrix& other) const 
{ 
// 首先检查行列数是否符合要求 
ASSERT (m_nNumColumns == other.GetNumRows()); 

// construct the object we are going to return 
CMatrix result(m_nNumRows, other.GetNumColumns()) ; 

// 矩阵乘法，即 
// 
// [A][B][C] [G][H] [A*G + B*I + C*K][A*H + B*J + C*L] 
// [D][E][F] * [I][J] = [D*G + E*I + F*K][D*H + E*J + F*L] 
// [K][L] 
// 
double value ; 
for (int i = 0 ; i < result.GetNumRows() ; ++i) 
{ 
for (int j = 0 ; j < other.GetNumColumns() ; ++j) 
{ 
value = 0.0 ; 
for (int k = 0 ; k < m_nNumColumns ; ++k) 
{ 
value += GetElement(i, k) * 
other.GetElement(k, j) ; 
} 

result.SetElement(i, j, value) ; 
} 
} 

return result ; 
} 

////////////////////////////////////////////////////////////////////// 
// 复矩阵的乘法 
// 
// 参数： 
// 1. const CMatrix& AR - 左边复矩阵的实部矩阵 
// 2. const CMatrix& AI - 左边复矩阵的虚部矩阵 
// 3. const CMatrix& BR - 右边复矩阵的实部矩阵 
// 4. const CMatrix& BI - 右边复矩阵的虚部矩阵 
// 5. CMatrix& CR - 乘积复矩阵的实部矩阵 
// 6. CMatrix& CI - 乘积复矩阵的虚部矩阵 
// 
// 返回值：BOOL型，复矩阵乘法是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::CMul(const CMatrix& AR, const CMatrix& AI, const CMatrix& BR, 
const CMatrix& BI, CMatrix& CR, CMatrix& CI) const 
{ 
// 首先检查行列数是否符合要求 
if (AR.GetNumColumns() != AI.GetNumColumns() || 
AR.GetNumRows() != AI.GetNumRows() || 
BR.GetNumColumns() != BI.GetNumColumns() || 
BR.GetNumRows() != BI.GetNumRows() || 
AR.GetNumColumns() != BR.GetNumRows()) 
return FALSE; 

// 构造乘积矩阵实部矩阵和虚部矩阵 
CMatrix mtxCR(AR.GetNumRows(), BR.GetNumColumns()), 
mtxCI(AR.GetNumRows(), BR.GetNumColumns()); 
// 复矩阵相乘 
for (int i=0; i<AR.GetNumRows(); ++i) 
{ 
for (int j=0; j<BR.GetNumColumns(); ++j) 
{ 
double vr = 0; 
double vi = 0; 
for (int k =0; k<AR.GetNumColumns(); ++k) 
{ 
double p = AR.GetElement(i, k) * BR.GetElement(k, j); 
double q = AI.GetElement(i, k) * BI.GetElement(k, j); 
double s = (AR.GetElement(i, k) + AI.GetElement(i, k)) * 
(BR.GetElement(k, j) + BI.GetElement(k, j)); 
vr += p - q; 
vi += s - p - q; 
} 
mtxCR.SetElement(i, j, vr); 
mtxCI.SetElement(i, j, vi); 
} 
} 

CR = mtxCR; 
CI = mtxCI; 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 矩阵的转置 
// 
// 参数：无 
// 
// 返回值：CMatrix型，指定矩阵转置矩阵 
////////////////////////////////////////////////////////////////////// 
CMatrix CMatrix::Transpose() const 
{ 
// 构造目标矩阵 
CMatrix Trans(m_nNumColumns, m_nNumRows); 

// 转置各元素 
for (int i = 0 ; i < m_nNumRows ; ++i) 
{ 
for (int j = 0 ; j < m_nNumColumns ; ++j) 
Trans.SetElement(j, i, GetElement(i, j)) ; 
} 

return Trans; 
} 

////////////////////////////////////////////////////////////////////// 
// 实矩阵求逆的全选主元高斯－约当法 
// 
// 参数：无 
// 
// 返回值：BOOL型，求逆是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::InvertGaussJordan() 
{ 
int *pnRow, *pnCol,i,j,k,l,u,v; 
double d = 0, p = 0; 

// 分配内存 
pnRow = new int[m_nNumColumns]; 
pnCol = new int[m_nNumColumns]; 
if (pnRow == NULL || pnCol == NULL) 
return FALSE; 

// 消元 
for (k=0; k<=m_nNumColumns-1; k++) 
{ 
d=0.0; 
for (i=k; i<=m_nNumColumns-1; i++) 
{ 
for (j=k; j<=m_nNumColumns-1; j++) 
{ 
l=i*m_nNumColumns+j; p=fabs(m_pData[l]); 
if (p>d) 
{ 
d=p; 
pnRow[k]=i; 
pnCol[k]=j; 
} 
} 
} 

// 失败 
if (d == 0.0) 
{ 
delete[] pnRow; 
delete[] pnCol; 
return FALSE; 
} 

if (pnRow[k] != k) 
{ 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
u=k*m_nNumColumns+j; 
v=pnRow[k]*m_nNumColumns+j; 
p=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=p; 
} 
} 

if (pnCol[k] != k) 
{ 
for (i=0; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns+k; 
v=i*m_nNumColumns+pnCol[k]; 
p=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=p; 
} 
} 

l=k*m_nNumColumns+k; 
m_pData[l]=1.0/m_pData[l]; 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
if (j != k) 
{ 
u=k*m_nNumColumns+j; 
m_pData[u]=m_pData[u]*m_pData[l]; 
} 
} 

for (i=0; i<=m_nNumColumns-1; i++) 
{ 
if (i!=k) 
{ 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
if (j!=k) 
{ 
u=i*m_nNumColumns+j; 

m_pData[u]=m_pData[u]-m_pData[i*m_nNumColumns+k]*m_pData[k*m_nNumColumns+j]; 
} 
} 
} 
} 

for (i=0; i<=m_nNumColumns-1; i++) 
{ 
if (i!=k) 
{ 
u=i*m_nNumColumns+k; 
m_pData[u]=-m_pData[u]*m_pData[l]; 
} 
} 
} 

// 调整恢复行列次序 
for (k=m_nNumColumns-1; k>=0; k--) 
{ 
if (pnCol[k]!=k) 
{ 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
u=k*m_nNumColumns+j; 
v=pnCol[k]*m_nNumColumns+j; 
p=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=p; 
} 
} 

if (pnRow[k]!=k) 
{ 
for (i=0; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns+k; 
v=i*m_nNumColumns+pnRow[k]; 
p=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=p; 
} 
} 
} 

// 清理内存 
delete[] pnRow; 
delete[] pnCol; 

// 成功返回 
return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 复矩阵求逆的全选主元高斯－约当法 
// 
// 参数： 
// 1. CMatrix& mtxImag - 复矩阵的虚部矩阵，当前矩阵为复矩阵的实部 
// 
// 返回值：BOOL型，求逆是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::InvertGaussJordan(CMatrix& mtxImag) 
{ 
int *pnRow,*pnCol,i,j,k,l,u,v,w; 
double p,q,s,t,d,b; 

// 分配内存 
pnRow = new int[m_nNumColumns]; 
pnCol = new int[m_nNumColumns]; 
if (pnRow == NULL || pnCol == NULL) 
return FALSE; 

// 消元 
for (k=0; k<=m_nNumColumns-1; k++) 
{ 
d=0.0; 
for (i=k; i<=m_nNumColumns-1; i++) 
{ 
for (j=k; j<=m_nNumColumns-1; j++) 
{ 
u=i*m_nNumColumns+j; 

p=m_pData[u]*m_pData[u]+mtxImag.m_pData[u]*mtxImag.m_pData[u]; 
if (p>d) 
{ 
d=p; 
pnRow[k]=i; 
pnCol[k]=j; 
} 
} 
} 

// 失败 
if (d == 0.0) 
{ 
delete[] pnRow; 
delete[] pnCol; 
return(0); 
} 

if (pnRow[k]!=k) 
{ 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
u=k*m_nNumColumns+j; 
v=pnRow[k]*m_nNumColumns+j; 
t=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=t; 
t=mtxImag.m_pData[u]; 
mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
mtxImag.m_pData[v]=t; 
} 
} 

if (pnCol[k]!=k) 
{ 
for (i=0; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns+k; 
v=i*m_nNumColumns+pnCol[k]; 
t=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=t; 
t=mtxImag.m_pData[u]; 
mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
mtxImag.m_pData[v]=t; 
} 
} 

l=k*m_nNumColumns+k; 
m_pData[l]=m_pData[l]/d; mtxImag.m_pData[l]=-mtxImag.m_pData[l]/d; 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
if (j!=k) 
{ 
u=k*m_nNumColumns+j; 
p=m_pData[u]*m_pData[l]; 
q=mtxImag.m_pData[u]*mtxImag.m_pData[l]; 

s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[l]+mtxImag.m_pData[l]); 
m_pData[u]=p-q; 
mtxImag.m_pData[u]=s-p-q; 
} 
} 

for (i=0; i<=m_nNumColumns-1; i++) 
{ 
if (i!=k) 
{ 
v=i*m_nNumColumns+k; 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
if (j!=k) 
{ 
u=k*m_nNumColumns+j; 
w=i*m_nNumColumns+j; 
p=m_pData[u]*m_pData[v]; 

q=mtxImag.m_pData[u]*mtxImag.m_pData[v]; 

s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[v]+mtxImag.m_pData[v]); 
t=p-q; 
b=s-p-q; 
m_pData[w]=m_pData[w]-t; 

mtxImag.m_pData[w]=mtxImag.m_pData[w]-b; 
} 
} 
} 
} 

for (i=0; i<=m_nNumColumns-1; i++) 
{ 
if (i!=k) 
{ 
u=i*m_nNumColumns+k; 
p=m_pData[u]*m_pData[l]; 
q=mtxImag.m_pData[u]*mtxImag.m_pData[l]; 

s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[l]+mtxImag.m_pData[l]); 
m_pData[u]=q-p; 
mtxImag.m_pData[u]=p+q-s; 
} 
} 
} 

// 调整恢复行列次序 
for (k=m_nNumColumns-1; k>=0; k--) 
{ 
if (pnCol[k]!=k) 
{ 
for (j=0; j<=m_nNumColumns-1; j++) 
{ 
u=k*m_nNumColumns+j; 
v=pnCol[k]*m_nNumColumns+j; 
t=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=t; 
t=mtxImag.m_pData[u]; 
mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
mtxImag.m_pData[v]=t; 
} 
} 

if (pnRow[k]!=k) 
{ 
for (i=0; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns+k; 
v=i*m_nNumColumns+pnRow[k]; 
t=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=t; 
t=mtxImag.m_pData[u]; 
mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
mtxImag.m_pData[v]=t; 
} 
} 
} 

// 清理内存 
delete[] pnRow; 
delete[] pnCol; 

// 成功返回 
return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 对称正定矩阵的求逆 
// 
// 参数：无 
// 
// 返回值：BOOL型，求逆是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::InvertSsgj() 
{ 
int i, j ,k, m; 
double w, g, *pTmp; 

// 临时内存 
pTmp = new double[m_nNumColumns]; 

// 逐列处理 
for (k=0; k<=m_nNumColumns-1; k++) 
{ 
w=m_pData[0]; 
if (w == 0.0) 
{ 
delete[] pTmp; 
return FALSE; 
} 

m=m_nNumColumns-k-1; 
for (i=1; i<=m_nNumColumns-1; i++) 
{ 
g=m_pData[i*m_nNumColumns]; 
pTmp[i]=g/w; 
if (i<=m) 
pTmp[i]=-pTmp[i]; 
for (j=1; j<=i; j++) 

m_pData[(i-1)*m_nNumColumns+j-1]=m_pData[i*m_nNumColumns+j]+g*pTmp[j]; 
} 

m_pData[m_nNumColumns*m_nNumColumns-1]=1.0/w; 
for (i=1; i<=m_nNumColumns-1; i++) 
m_pData[(m_nNumColumns-1)*m_nNumColumns+i-1]=pTmp[i]; 
} 

// 行列调整 
for (i=0; i<=m_nNumColumns-2; i++) 
for (j=i+1; j<=m_nNumColumns-1; j++) 
m_pData[i*m_nNumColumns+j]=m_pData[j*m_nNumColumns+i]; 

// 临时内存清理 
delete[] pTmp; 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 托伯利兹矩阵求逆的埃兰特方法 
// 
// 参数：无 
// 
// 返回值：BOOL型，求逆是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::InvertTrench() 
{ 
int i,j,k; 
double a,s,*t,*tt,*c,*r,*p; 

// 上三角元素 
t = new double[m_nNumColumns]; 
// 下三角元素 
tt = new double[m_nNumColumns]; 

// 上、下三角元素赋值 
for (i=0; i<m_nNumColumns; ++i) 
{ 
t[i] = GetElement(0, i); 
tt[i] = GetElement(i, 0); 
} 

// 临时缓冲区 
c = new double[m_nNumColumns]; 
r = new double[m_nNumColumns]; 
p = new double[m_nNumColumns]; 

// 非Toeplitz矩阵，返回 
if (t[0] == 0.0) 
{ 
delete[] t; 
delete[] tt; 
delete[] c; 
delete[] r; 
delete[] p; 
return FALSE; 
} 

a=t[0]; 
c[0]=tt[1]/t[0]; 
r[0]=t[1]/t[0]; 

for (k=0; k<=m_nNumColumns-3; k++) 
{ 
s=0.0; 
for (j=1; j<=k+1; j++) 
s=s+c[k+1-j]*tt[j]; 

s=(s-tt[k+2])/a; 
for (i=0; i<=k; i++) 
p[i]=c[i]+s*r[k-i]; 

c[k+1]=-s; 
s=0.0; 
for (j=1; j<=k+1; j++) 
s=s+r[k+1-j]*t[j]; 

s=(s-t[k+2])/a; 
for (i=0; i<=k; i++) 
{ 
r[i]=r[i]+s*c[k-i]; 
c[k-i]=p[k-i]; 
} 

r[k+1]=-s; 
a=0.0; 
for (j=1; j<=k+2; j++) 
a=a+t[j]*c[j-1]; 

a=t[0]-a; 

// 求解失败 
if (a == 0.0) 
{ 
delete[] t; 
delete[] tt; 
delete[] c; 
delete[] r; 
delete[] p; 
return FALSE; 
} 
} 

m_pData[0]=1.0/a; 
for (i=0; i<=m_nNumColumns-2; i++) 
{ 
k=i+1; 
j=(i+1)*m_nNumColumns; 
m_pData[k]=-r[i]/a; 
m_pData[j]=-c[i]/a; 
} 

for (i=0; i<=m_nNumColumns-2; i++) 
{