操作矩阵的类 cmatrix.txt

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

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

m_pData[k]=m_pData[i*m_nNumColumns+j]-c[i]*m_pData[j+1]; 

m_pData[k]=m_pData[k]+c[m_nNumColumns-j-2]*m_pData[m_nNumColumns-i-1]; 
} 
} 

// 临时内存清理 
delete[] t; 
delete[] tt; 
delete[] c; 
delete[] r; 
delete[] p; 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 求行列式值的全选主元高斯消去法 
// 
// 参数：无 
// 
// 返回值：double型，行列式的值 
////////////////////////////////////////////////////////////////////// 
double CMatrix::DetGauss() 
{ 
int i,j,k,is,js,l,u,v; 
double f,det,q,d; 

// 初值 
f=1.0; 
det=1.0; 

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

if (q == 0.0) 
{ 
det=0.0; 
return(det); 
} 

if (is!=k) 
{ 
f=-f; 
for (j=k; j<=m_nNumColumns-1; j++) 
{ 
u=k*m_nNumColumns+j; 
v=is*m_nNumColumns+j; 
d=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=d; 
} 
} 

if (js!=k) 
{ 
f=-f; 
for (i=k; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns+js; 
v=i*m_nNumColumns+k; 
d=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=d; 
} 
} 

l=k*m_nNumColumns+k; 
det=det*m_pData[l]; 
for (i=k+1; i<=m_nNumColumns-1; i++) 
{ 
d=m_pData[i*m_nNumColumns+k]/m_pData[l]; 
for (j=k+1; j<=m_nNumColumns-1; j++) 
{ 
u=i*m_nNumColumns+j; 
m_pData[u]=m_pData[u]-d*m_pData[k*m_nNumColumns+j]; 
} 
} 
} 

// 求值 
det=f*det*m_pData[m_nNumColumns*m_nNumColumns-1]; 

return(det); 
} 

////////////////////////////////////////////////////////////////////// 
// 求矩阵秩的全选主元高斯消去法 
// 
// 参数：无 
// 
// 返回值：int型，矩阵的秩 
////////////////////////////////////////////////////////////////////// 
int CMatrix::RankGauss() 
{ 
int i,j,k,nn,is,js,l,ll,u,v; 
double q,d; 

// 秩小于等于行列数 
nn = m_nNumRows; 
if (m_nNumRows >= m_nNumColumns) 
nn = m_nNumColumns; 

k=0; 

// 消元求解 
for (l=0; l<=nn-1; l++) 
{ 
q=0.0; 
for (i=l; i<=m_nNumRows-1; i++) 
{ 
for (j=l; j<=m_nNumColumns-1; j++) 
{ 
ll=i*m_nNumColumns+j; 
d=fabs(m_pData[ll]); 
if (d>q) 
{ 
q=d; 
is=i; 
js=j; 
} 
} 
} 

if (q == 0.0) 
return(k); 

k=k+1; 
if (is!=l) 
{ 
for (j=l; j<=m_nNumColumns-1; j++) 
{ 
u=l*m_nNumColumns+j; 
v=is*m_nNumColumns+j; 
d=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=d; 
} 
} 
if (js!=l) 
{ 
for (i=l; i<=m_nNumRows-1; i++) 
{ 
u=i*m_nNumColumns+js; 
v=i*m_nNumColumns+l; 
d=m_pData[u]; 
m_pData[u]=m_pData[v]; 
m_pData[v]=d; 
} 
} 

ll=l*m_nNumColumns+l; 
for (i=l+1; i<=m_nNumColumns-1; i++) 
{ 
d=m_pData[i*m_nNumColumns+l]/m_pData[ll]; 
for (j=l+1; j<=m_nNumColumns-1; j++) 
{ 
u=i*m_nNumColumns+j; 
m_pData[u]=m_pData[u]-d*m_pData[l*m_nNumColumns+j]; 
} 
} 
} 

return(k); 
} 

////////////////////////////////////////////////////////////////////// 
// 对称正定矩阵的乔里斯基分解与行列式的求值 
// 
// 参数： 
// 1. double* dblDet - 返回行列式的值 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::DetCholesky(double* dblDet) 
{ 
int i,j,k,u,l; 
double d; 

// 不满足求解要求 
if (m_pData[0] <= 0.0) 
return FALSE; 

// 乔里斯基分解 

m_pData[0]=sqrt(m_pData[0]); 
d=m_pData[0]; 

for (i=1; i<=m_nNumColumns-1; i++) 
{ 
u=i*m_nNumColumns; 
m_pData[u]=m_pData[u]/m_pData[0]; 
} 

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

if (m_pData[l] <= 0.0) 
return FALSE; 

m_pData[l]=sqrt(m_pData[l]); 
d=d*m_pData[l]; 

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

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

m_pData[u]=m_pData[u]/m_pData[l]; 
} 
} 

// 行列式求值 
*dblDet=d*d; 

// 下三角矩阵 
for (i=0; i<=m_nNumColumns-2; i++) 
for (j=i+1; j<=m_nNumColumns-1; j++) 
m_pData[i*m_nNumColumns+j]=0.0; 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 矩阵的三角分解，分解成功后，原矩阵将成为Q矩阵 
// 
// 参数： 
// 1. CMatrix& mtxL - 返回分解后的L矩阵 
// 2. CMatrix& mtxU - 返回分解后的U矩阵 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::SplitLU(CMatrix& mtxL, CMatrix& mtxU) 
{ 
int i,j,k,w,v,ll; 

// 初始化结果矩阵 
if (! mtxL.Init(m_nNumColumns, m_nNumColumns) || 
! mtxU.Init(m_nNumColumns, m_nNumColumns)) 
return FALSE; 

for (k=0; k<=m_nNumColumns-2; k++) 
{ 
ll=k*m_nNumColumns+k; 
if (m_pData[ll] == 0.0) 
return FALSE; 

for (i=k+1; i<=m_nNumColumns-1; i++) 
{ 
w=i*m_nNumColumns+k; 
m_pData[w]=m_pData[w]/m_pData[ll]; 
} 

for (i=k+1; i<=m_nNumColumns-1; i++) 
{ 
w=i*m_nNumColumns+k; 
for (j=k+1; j<=m_nNumColumns-1; j++) 
{ 
v=i*m_nNumColumns+j; 
m_pData[v]=m_pData[v]-m_pData[w]*m_pData[k*m_nNumColumns+j]; 
} 
} 
} 

for (i=0; i<=m_nNumColumns-1; i++) 
{ 
for (j=0; j<i; j++) 
{ 
w=i*m_nNumColumns+j; 
mtxL.m_pData[w]=m_pData[w]; 
mtxU.m_pData[w]=0.0; 
} 

w=i*m_nNumColumns+i; 
mtxL.m_pData[w]=1.0; 
mtxU.m_pData[w]=m_pData[w]; 

for (j=i+1; j<=m_nNumColumns-1; j++) 
{ 
w=i*m_nNumColumns+j; 
mtxL.m_pData[w]=0.0; 
mtxU.m_pData[w]=m_pData[w]; 
} 
} 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 一般实矩阵的QR分解，分解成功后，原矩阵将成为R矩阵 
// 
// 参数： 
// 1. CMatrix& mtxQ - 返回分解后的Q矩阵 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::SplitQR(CMatrix& mtxQ) 
{ 
int i,j,k,l,nn,p,jj; 
double u,alpha,w,t; 

if (m_nNumRows < m_nNumColumns) 
return FALSE; 

// 初始化Q矩阵 
if (! mtxQ.Init(m_nNumRows, m_nNumRows)) 
return FALSE; 

// 对角线元素单位化 
for (i=0; i<=m_nNumRows-1; i++) 
{ 
for (j=0; j<=m_nNumRows-1; j++) 
{ 
l=i*m_nNumRows+j; 
mtxQ.m_pData[l]=0.0; 
if (i==j) 
mtxQ.m_pData[l]=1.0; 
} 
} 

// 开始分解 

nn=m_nNumColumns; 
if (m_nNumRows == m_nNumColumns) 
nn=m_nNumRows-1; 

for (k=0; k<=nn-1; k++) 
{ 
u=0.0; 
l=k*m_nNumColumns+k; 
for (i=k; i<=m_nNumRows-1; i++) 
{ 
w=fabs(m_pData[i*m_nNumColumns+k]); 
if (w>u) 
u=w; 
} 

alpha=0.0; 
for (i=k; i<=m_nNumRows-1; i++) 
{ 
t=m_pData[i*m_nNumColumns+k]/u; 
alpha=alpha+t*t; 
} 

if (m_pData[l]>0.0) 
u=-u; 

alpha=u*sqrt(alpha); 
if (alpha == 0.0) 
return FALSE; 

u=sqrt(2.0*alpha*(alpha-m_pData[l])); 
if ((u+1.0)!=1.0) 
{ 
m_pData[l]=(m_pData[l]-alpha)/u; 
for (i=k+1; i<=m_nNumRows-1; i++) 
{ 
p=i*m_nNumColumns+k; 
m_pData[p]=m_pData[p]/u; 
} 

for (j=0; j<=m_nNumRows-1; j++) 
{ 
t=0.0; 
for (jj=k; jj<=m_nNumRows-1; jj++) 

t=t+m_pData[jj*m_nNumColumns+k]*mtxQ.m_pData[jj*m_nNumRows+j]; 

for (i=k; i<=m_nNumRows-1; i++) 
{ 
p=i*m_nNumRows+j; 

mtxQ.m_pData[p]=mtxQ.m_pData[p]-2.0*t*m_pData[i*m_nNumColumns+k]; 
} 
} 

for (j=k+1; j<=m_nNumColumns-1; j++) 
{ 
t=0.0; 

for (jj=k; jj<=m_nNumRows-1; jj++) 

t=t+m_pData[jj*m_nNumColumns+k]*m_pData[jj*m_nNumColumns+j]; 

for (i=k; i<=m_nNumRows-1; i++) 
{ 
p=i*m_nNumColumns+j; 

m_pData[p]=m_pData[p]-2.0*t*m_pData[i*m_nNumColumns+k]; 
} 
} 

m_pData[l]=alpha; 
for (i=k+1; i<=m_nNumRows-1; i++) 
m_pData[i*m_nNumColumns+k]=0.0; 
} 
} 

// 调整元素 
for (i=0; i<=m_nNumRows-2; i++) 
{ 
for (j=i+1; j<=m_nNumRows-1;j++) 
{ 
p=i*m_nNumRows+j; 
l=j*m_nNumRows+i; 
t=mtxQ.m_pData[p]; 
mtxQ.m_pData[p]=mtxQ.m_pData[l]; 
mtxQ.m_pData[l]=t; 
} 
} 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 一般实矩阵的奇异值分解，分解成功后，原矩阵对角线元素就是矩阵的奇异值 
// 
// 参数： 
// 1. CMatrix& mtxU - 返回分解后的U矩阵 
// 2. CMatrix& mtxV - 返回分解后的V矩阵 
// 3. double eps - 计算精度，默认值为0.000001 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::SplitUV(CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/) 
{ 
int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks; 
double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2]; 
double *s,*e,*w; 

int m = m_nNumRows; 
int n = m_nNumColumns; 

// 初始化U, V矩阵 
if (! mtxU.Init(m, m) || ! mtxV.Init(n, n)) 
return FALSE; 

// 临时缓冲区 
int ka = max(m, n) + 1; 
s = new double[ka]; 
e = new double[ka]; 
w = new double[ka]; 

// 指定迭代次数为60 
it=60; 
k=n; 

if (m-1<n) 
k=m-1; 

l=m; 
if (n-2<m) 
l=n-2; 
if (l<0) 
l=0; 

// 循环迭代计算 
ll=k; 
if (l>k) 
ll=l; 
if (ll>=1) 
{ 
for (kk=1; kk<=ll; kk++) 
{ 
if (kk<=k) 
{ 
d=0.0; 
for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*n+kk-1; 
d=d+m_pData[ix]*m_pData[ix]; 
} 

s[kk-1]=sqrt(d); 
if (s[kk-1]!=0.0) 
{ 
ix=(kk-1)*n+kk-1; 
if (m_pData[ix]!=0.0) 
{ 
s[kk-1]=fabs(s[kk-1]); 
if (m_pData[ix]<0.0) 
s[kk-1]=-s[kk-1]; 
} 

for (i=kk; i<=m; i++) 
{ 
iy=(i-1)*n+kk-1; 
m_pData[iy]=m_pData[iy]/s[kk-1]; 
} 

m_pData[ix]=1.0+m_pData[ix]; 
} 

s[kk-1]=-s[kk-1]; 
} 

if (n>=kk+1) 
{ 
for (j=kk+1; j<=n; j++) 
{ 
if ((kk<=k)&&(s[kk-1]!=0.0)) 
{ 
d=0.0; 
for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*n+kk-1; 
iy=(i-1)*n+j-1; 
d=d+m_pData[ix]*m_pData[iy]; 
} 

d=-d/m_pData[(kk-1)*n+kk-1]; 
for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*n+j-1; 
iy=(i-1)*n+kk-1; 
m_pData[ix]=m_pData[ix]+d*m_pData[iy]; 
} 
} 

e[j-1]=m_pData[(kk-1)*n+j-1]; 
} 
} 

if (kk<=k) 
{ 
for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*m+kk-1; 
iy=(i-1)*n+kk-1; 
mtxU.m_pData[ix]=m_pData[iy]; 
} 
} 

if (kk<=l) 
{ 
d=0.0; 
for (i=kk+1; i<=n; i++) 
d=d+e[i-1]*e[i-1]; 

e[kk-1]=sqrt(d); 
if (e[kk-1]!=0.0) 
{ 
if (e[kk]!=0.0) 
{ 
e[kk-1]=fabs(e[kk-1]); 
if (e[kk]<0.0) 
e[kk-1]=-e[kk-1]; 
} 

for (i=kk+1; i<=n; i++) 
e[i-1]=e[i-1]/e[kk-1]; 

e[kk]=1.0+e[kk]; 
} 

e[kk-1]=-e[kk-1]; 
if ((kk+1<=m)&&(e[kk-1]!=0.0)) 
{ 
for (i=kk+1; i<=m; i++) 
w[i-1]=0.0; 

for (j=kk+1; j<=n; j++) 
for (i=kk+1; i<=m; i++) 

w[i-1]=w[i-1]+e[j-1]*m_pData[(i-1)*n+j-1]; 

for (j=kk+1; j<=n; j++) 
{ 
for (i=kk+1; i<=m; i++) 
{ 
ix=(i-1)*n+j-1; 

m_pData[ix]=m_pData[ix]-w[i-1]*e[j-1]/e[kk]; 
} 
} 
} 

for (i=kk+1; i<=n; i++) 
mtxV.m_pData[(i-1)*n+kk-1]=e[i-1]; 
} 
} 
} 

mm=n; 
if (m+1<n) 
mm=m+1; 
if (k<n) 
s[k]=m_pData[k*n+k]; 
if (m<mm) 
s[mm-1]=0.0; 
if (l+1<mm) 
e[l]=m_pData[l*n+mm-1]; 

e[mm-1]=0.0; 
nn=m; 
if (m>n) 
nn=n; 
if (nn>=k+1) 
{ 
for (j=k+1; j<=nn; j++) 
{ 
for (i=1; i<=m; i++) 
mtxU.m_pData[(i-1)*m+j-1]=0.0; 
mtxU.m_pData[(j-1)*m+j-1]=1.0; 
} 
}