// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "Matrix.h"
#include "math.h"
#include "string.h"
#include "fstream.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
// 无参数构造函数
CMatrix::CMatrix()
{
m_nRows = 1;
m_nCols = 1;
// if(m_pData != NULL)
// delete [] m_pData;
m_pData = NULL;
if((m_pData = new double[m_nRows * m_nCols]) == NULL)
{
AfxMessageBox("没有空间分配内存");
exit(1);
}
}
// 给定矩阵大小后的构造函数
CMatrix::CMatrix(int nRows, int nCols)
{
m_nRows = nRows;
m_nCols = nCols;
// if(m_pData != NULL)
// delete [] m_pData;
m_pData = NULL;
if((m_pData = new double[m_nRows * m_nCols]) == NULL)
{
AfxMessageBox("没有空间分配内存");
exit(1);
}
memset(m_pData, 0, sizeof(double) * m_nRows * m_nCols);
}
// 给定行列数及矩阵元素值后的构造函数
CMatrix::CMatrix(int nRows, int nCols, double dvalue[])
{
m_nRows = nRows;
m_nCols = nCols;
// if(m_pData != NULL)
// delete [] m_pData;
m_pData = NULL;
if((m_pData = new double[m_nRows * m_nCols]) == NULL)
{
AfxMessageBox("没有空间分配内存");
exit(1);
}
memset(m_pData, 0, sizeof(double) * m_nRows * m_nCols);
memcpy(m_pData, dvalue, sizeof(double) * m_nRows * m_nCols);
}
CMatrix::CMatrix(CMatrix *mat)
{
m_nRows = mat -> m_nRows;
m_nCols = mat -> m_nCols;
// if(m_pData != NULL)
// delete [] m_pData;
m_pData = NULL;
if((m_pData = new double[m_nRows * m_nCols]) == NULL)
{
AfxMessageBox("没有空间分配内存");
exit(1);
}
memset(m_pData, 0, sizeof(double) * m_nRows * m_nCols);
memcpy(m_pData, mat -> m_pData, sizeof(double) * m_nRows * m_nCols);
}
//虚构函数
CMatrix::~CMatrix()
{
if(m_pData != NULL)
delete [] m_pData;
}
bool CMatrix::Init( int nRows, int nCols )
{
if (m_pData)
{
delete [] m_pData;
m_pData = NULL;
}
m_nRows = nRows;
m_nCols = nCols;
int nSize = m_nRows*m_nCols;
if(nSize<0)
{
return false;
}
if((m_pData = new double[nSize]) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
if(m_pData == NULL)
{
return false;
}
memset(m_pData, 0, sizeof(double) * m_nRows * m_nCols);
return true;
}
int CMatrix::GetColNum()
{
return m_nCols;
}
int CMatrix::GetRowNum()
{
return m_nRows;
}
double CMatrix::GetElem(int nx, int ny)
{
double Num = m_pData[nx * m_nCols + ny];
return Num;
}
bool CMatrix::SetElem(int nx, int ny, double dNum)
{
m_pData[nx * m_nCols + ny] = dNum;
return true;
}
bool CMatrix::MatrixMul(int nARows,int nACols,int nBRows,int nBCols,double dAVal[],double dBVal[])
{
int i,j,l,u;
int nAColsandBRows = nBRows;
if (nACols != nBRows)
{
AfxMessageBox("乘法错误");
return false;
}
Init(nARows,nBCols);
for (i = 0; i < nARows; i++)
{ for (j = 0; j < nBCols; j++)
{
u = i * nBCols + j;
m_pData[u] = 0;
for (l = 0; l < nAColsandBRows; l++)
m_pData[u] = m_pData[u] + dAVal[i * nAColsandBRows + l] * dBVal[l * nBCols + j];
}
}
return true;
}
bool CMatrix::MatrixInverse()
{
if (m_nRows != m_nCols)
{
AfxMessageBox("不是方阵");
return false;
}
int n = m_nRows;
int *is,*js;
int i,j,k,l,u,v;
k = 0;
double d,p;
if((is = new int(n*sizeof(int))) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
if((js = new int(n*sizeof(int))) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
for (k=0; k<=n-1; k++)
{
d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs(m_pData[l]);
if (p>d)
{ d=p;
is[k]=i;
js[k]=j;}
}
if (d+1.0==1.0)
{
delete [] is;
delete [] js; //释放总出错，没解决
AfxMessageBox("奇异阵");
return false;
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is[k]*n+j;
p=m_pData[u]; m_pData[u]=m_pData[v]; m_pData[v]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
p=m_pData[u];
m_pData[u]=m_pData[v];
m_pData[v]=p;
}
l=k*n+k;
m_pData[l]=1.0/m_pData[l];
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j;
m_pData[u]=m_pData[u]*m_pData[l];
}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
m_pData[u]=m_pData[u]-m_pData[i*n+k]*m_pData[k*n+j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; m_pData[u]=-m_pData[u]*m_pData[l];}
}
for (k=n-1; k>=0; k--)
{ if (js[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js[k]*n+j;
p=m_pData[u]; m_pData[u]=m_pData[v]; m_pData[v]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is[k];
p=m_pData[u]; m_pData[u]=m_pData[v]; m_pData[v]=p;
}
}
delete is;
delete js;//释放总出错，没解决
return true;
}
bool CMatrix::MatrixTrs(int nRows,int nCols,double dval[])
{
Init(nCols,nRows);
for(int i = 0;i < nRows;i++)
{
for(int j = 0;j < nCols;j++)
{ m_pData[j * m_nCols + i] = dval[i * nCols + j];}
}
return true;
}
double CMatrix::VectorColMax()
{
if (m_nCols != 1)
{
AfxMessageBox("不是列向量");
return false;
}
double temp = 0;
for (int i = 0;i < m_nRows;i++)
{
temp = max(temp,fabs(m_pData[i]));
}
return temp;
}
int CMatrix::VectorColMaxPosition()
{
if (m_nCols != 1)
{
AfxMessageBox("不是列向量");
return false;
}
double temp = 0;
int times = 0;
for (int i = 0;i < m_nRows;i++)
{
if (fabs(m_pData[i]) > temp)
{
temp = fabs(m_pData[i]);
times = i;
}
}
return times;
}
bool CMatrix::MatrixSVD( CMatrix& mtxU, CMatrix& mtxV, double eps )
{
int i,j;
int k; // k=min(n,m-1)
int l; // l=min(m,n-2)
int it; // 迭代次数
int 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_nRows;
int n = m_nCols;
// 初始化U, V矩阵
if (! mtxU.Init(m, m) || ! mtxV.Init(n, n))
{
return false;
}
// 临时缓冲区
int ka = max(m, n) + 1;
if((s = new double[ka]) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
if((e = new double[ka]) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
if((w = new double[ka]) == NULL)
{
AfxMessageBox("没有空间分配内存");
return false;
}
// 指定迭代次数为60
it=60;
// 确定变量k,l的值
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>
