?? 操作矩陣的類 cmatrix.txt
字號:
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;
}
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -