[分享]和矩阵运算相关的类代码
CMatrix::CMatrix(){
m_nNumColumns = 1;
m_nNumRows = 1;
m_pData = NULL;
BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
ASSERT(bSuccess);
}
//////////////////////////////////////////////////////////////////////
// 指定行列构造函数
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. int nCols - 指定的矩阵列数
//////////////////////////////////////////////////////////////////////
CMatrix::CMatrix(int nRows, int nCols)
{
m_nNumRows = nRows;
m_nNumColumns = nCols;
m_pData = NULL;
BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
ASSERT(bSuccess);
}
//////////////////////////////////////////////////////////////////////
// 指定值构造函数
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. int nCols - 指定的矩阵列数
// 3. double value[] - 一维数组，长度为nRows*nCols，存储矩阵各元素的值
//////////////////////////////////////////////////////////////////////
CMatrix::CMatrix(int nRows, int nCols, double value[])
{
m_nNumRows = nRows;
m_nNumColumns = nCols;
m_pData = NULL;
BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
ASSERT(bSuccess);
SetData(value);
}
//////////////////////////////////////////////////////////////////////
// 方阵构造函数
//
// 参数：
// 1. int nSize - 方阵行列数
//////////////////////////////////////////////////////////////////////
CMatrix::CMatrix(int nSize)
{
m_nNumRows = nSize;
m_nNumColumns = nSize;
m_pData = NULL;
BOOL bSuccess = Init(nSize, nSize);
ASSERT (bSuccess);
}
//////////////////////////////////////////////////////////////////////
// 方阵构造函数
//
// 参数：
// 1. int nSize - 方阵行列数
// 2. double value[] - 一维数组，长度为nRows*nRows，存储方阵各元素的值
//////////////////////////////////////////////////////////////////////
CMatrix::CMatrix(int nSize, double value[])
{
m_nNumRows = nSize;
m_nNumColumns = nSize;
m_pData = NULL;
BOOL bSuccess = Init(nSize, nSize);
ASSERT (bSuccess);
SetData(value);
}
//////////////////////////////////////////////////////////////////////
// 拷贝构造函数
//
// 参数：
// 1. const CMatrix& other - 源矩阵
//////////////////////////////////////////////////////////////////////
CMatrix::CMatrix(const CMatrix& other)
{
m_nNumColumns = other.GetNumColumns();
m_nNumRows = other.GetNumRows();
m_pData = NULL;
BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
ASSERT(bSuccess);
// copy the pointer
memcpy(m_pData, other.m_pData, sizeof(double)*m_nNumColumns*m_nNumRows);
}
//////////////////////////////////////////////////////////////////////
// 析构函数
//////////////////////////////////////////////////////////////////////
CMatrix::~CMatrix()
{
if (m_pData)
{
delete[] m_pData;
m_pData = NULL;
}
}
//////////////////////////////////////////////////////////////////////
// 初始化函数
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. int nCols - 指定的矩阵列数
//
// 返回值：BOOL 型，初始化是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::Init(int nRows, int nCols)
{
if (m_pData)
{
delete[] m_pData;
m_pData = NULL;
}
m_nNumRows = nRows;
m_nNumColumns = nCols;
int nSize = nCols*nRows;
if (nSize < 0)
return FALSE;
// 分配内存
m_pData = new double;
if (m_pData == NULL)
return FALSE; // 内存分配失败
if (IsBadReadPtr(m_pData, sizeof(double) * nSize))
return FALSE;
// 将各元素值置0
memset(m_pData, 0, sizeof(double) * nSize);
return TRUE;
}
//////////////////////////////////////////////////////////////////////
// 将方阵初始化为单位矩阵
//
// 参数：
// 1. int nSize - 方阵行列数
//
// 返回值：BOOL 型，初始化是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::MakeUnitMatrix(int nSize)
{
if (! Init(nSize, nSize))
return FALSE;
for (int i=0; i<nSize; ++i)
for (int j=0; j<nSize; ++j)
if (i == j)
SetElement(i, j, 1);
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:03 编辑 ] //////////////////////////////////////////////////////////////////////
// 将字符串转化为矩阵的值
//
// 参数：
// 1. CString s - 数字和分隔符构成的字符串
// 2. const CString& sDelim - 数字之间的分隔符，默认为空格
// 3. BOOL bLineBreak - 行与行之间是否有回车换行符，默认为真(有换行符)
// 当该参数为FALSE时，所有元素值都在一行中输入，字符串的第一个
// 数值应为矩阵的行数，第二个数值应为矩阵的列数
//
// 返回值：BOOL 型，转换是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::FromString(CString s, const CString& sDelim /*= " "*/, BOOL bLineBreak /*= TRUE*/)
{
if (s.IsEmpty())
return FALSE;
// 分行处理
if (bLineBreak)
{
CTokenizer tk(s, "\r\n");
CStringList ListRow;
CString sRow;
while (tk.Next(sRow))
{
sRow.TrimLeft();
sRow.TrimRight();
if (sRow.IsEmpty())
break;
ListRow.AddTail(sRow);
}
// 行数
m_nNumRows = ListRow.GetCount();
sRow = ListRow.GetHead();
CTokenizer tkRow(sRow, sDelim);
CString sElement;
// 列数
m_nNumColumns = 0;
while (tkRow.Next(sElement))
{
m_nNumColumns++;
}
// 初始化矩阵
if (! Init(m_nNumRows, m_nNumColumns))
return FALSE;
// 设置值
POSITION pos = ListRow.GetHeadPosition();
for (int i=0; i<m_nNumRows; i++)
{
sRow = ListRow.GetNext(pos);
int j = 0;
CTokenizer tkRow(sRow, sDelim);
while (tkRow.Next(sElement))
{
sElement.TrimLeft();
sElement.TrimRight();
double v = atof(sElement);
SetElement(i, j++, v);
}
}
return TRUE;
}
// 不分行（单行）处理
CTokenizer tk(s, sDelim);
CString sElement;
// 行数
tk.Next(sElement);
sElement.TrimLeft();
sElement.TrimRight();
m_nNumRows = atoi(sElement);
// 列数
tk.Next(sElement);
sElement.TrimLeft();
sElement.TrimRight();
m_nNumColumns = atoi(sElement);
// 初始化矩阵
if (! Init(m_nNumRows, m_nNumColumns))
return FALSE;
// 设置值
int i = 0, j = 0;
while (tk.Next(sElement))
{
sElement.TrimLeft();
sElement.TrimRight();
double v = atof(sElement);
SetElement(i, j++, v);
if (j == m_nNumColumns)
{
j = 0;
i++;
if (i == m_nNumRows)
break;
}
}
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:03 编辑 ] //////////////////////////////////////////////////////////////////////
// 将矩阵各元素的值转化为字符串
//
// 参数：
// 1. const CString& sDelim - 数字之间的分隔符，默认为空格
// 2 BOOL bLineBreak - 行与行之间是否有回车换行符，默认为真(有换行符)
//
// 返回值：CString 型，转换得到的字符串
//////////////////////////////////////////////////////////////////////
CString CMatrix::ToString(const CString& sDelim /*= " "*/, BOOL bLineBreak /*= TRUE*/) const
{
CString s="";
for (int i=0; i<m_nNumRows; ++i)
{
for (int j=0; j<m_nNumColumns; ++j)
{
CString ss;
ss.Format("%f", GetElement(i, j));
s += ss;
if (bLineBreak)
{
if (j != m_nNumColumns-1)
s += sDelim;
}
else
{
if (i != m_nNumRows-1 || j != m_nNumColumns-1)
s += sDelim;
}
}
if (bLineBreak)
if (i != m_nNumRows-1)
s += "\r\n";
}
return s;
}
//////////////////////////////////////////////////////////////////////
// 将矩阵指定行中各元素的值转化为字符串
//
// 参数：
// 1. int nRow - 指定的矩阵行，nRow = 0表示第一行
// 2. const CString& sDelim - 数字之间的分隔符，默认为空格
//
// 返回值：CString 型，转换得到的字符串
//////////////////////////////////////////////////////////////////////
CString CMatrix::RowToString(int nRow, const CString& sDelim /*= " "*/) const
{
CString s = "";
if (nRow >= m_nNumRows)
return s;
for (int j=0; j<m_nNumColumns; ++j)
{
CString ss;
ss.Format("%f", GetElement(nRow, j));
s += ss;
if (j != m_nNumColumns-1)
s += sDelim;
}
return s;
}
//////////////////////////////////////////////////////////////////////
// 将矩阵指定列中各元素的值转化为字符串
//
// 参数：
// 1. int nCol - 指定的矩阵行，nCol = 0表示第一列
// 2. const CString& sDelim - 数字之间的分隔符，默认为空格
//
// 返回值：CString 型，转换得到的字符串
//////////////////////////////////////////////////////////////////////
CString CMatrix::ColToString(int nCol, const CString& sDelim /*= " "*/) const
{
CString s = "";
if (nCol >= m_nNumColumns)
return s;
for (int i=0; i<m_nNumRows; ++i)
{
CString ss;
ss.Format("%f", GetElement(i, nCol));
s += ss;
if (i != m_nNumRows-1)
s += sDelim;
}
return s;
}
//////////////////////////////////////////////////////////////////////
// 设置矩阵各元素的值
//
// 参数：
// 1. double value[] - 一维数组，长度为m_nNumColumns*m_nNumRows，存储
// 矩阵各元素的值
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CMatrix::SetData(double value[])
{
// empty the memory
memset(m_pData, 0, sizeof(double) * m_nNumColumns*m_nNumRows);
// copy data
memcpy(m_pData, value, sizeof(double)*m_nNumColumns*m_nNumRows);
}
//////////////////////////////////////////////////////////////////////
// 设置指定元素的值
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. int nCols - 指定的矩阵列数
// 3. double value - 指定元素的值
//
// 返回值：BOOL 型，说明设置是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::SetElement(int nRow, int nCol, double value)
{
if (nCol < 0 || nCol >= m_nNumColumns || nRow < 0 || nRow >= m_nNumRows)
return FALSE; // array bounds error
if (m_pData == NULL)
return FALSE; // bad pointer error
m_pData = value;
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:04 编辑 ] //////////////////////////////////////////////////////////////////////
// 设置指定元素的值
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. int nCols - 指定的矩阵列数
//
// 返回值：double 型，指定元素的值
//////////////////////////////////////////////////////////////////////
double CMatrix::GetElement(int nRow, int nCol) const
{
ASSERT(nCol >= 0 && nCol < m_nNumColumns && nRow >= 0 && nRow < m_nNumRows); // array bounds error
ASSERT(m_pData); // bad pointer error
return m_pData ;
}
//////////////////////////////////////////////////////////////////////
// 获取矩阵的列数
//
// 参数：无
//
// 返回值：int 型，矩阵的列数
//////////////////////////////////////////////////////////////////////
int CMatrix::GetNumColumns() const
{
return m_nNumColumns;
}
//////////////////////////////////////////////////////////////////////
// 获取矩阵的行数
//
// 参数：无
//
// 返回值：int 型，矩阵的行数
//////////////////////////////////////////////////////////////////////
int CMatrix::GetNumRows() const
{
return m_nNumRows;
}
//////////////////////////////////////////////////////////////////////
// 获取矩阵的数据
//
// 参数：无
//
// 返回值：double型指针，指向矩阵各元素的数据缓冲区
//////////////////////////////////////////////////////////////////////
double* CMatrix::GetData() const
{
return m_pData;
}
//////////////////////////////////////////////////////////////////////
// 获取指定行的向量
//
// 参数：
// 1. int nRows - 指定的矩阵行数
// 2. double* pVector - 指向向量中各元素的缓冲区
//
// 返回值：int 型，向量中元素的个数，即矩阵的列数
//////////////////////////////////////////////////////////////////////
int CMatrix::GetRowVector(int nRow, double* pVector) const
{
if (pVector == NULL)
delete pVector;
pVector = new double;
ASSERT(pVector != NULL);
for (int j=0; j<m_nNumColumns; ++j)
pVector = GetElement(nRow, j);
return m_nNumColumns;
}
//////////////////////////////////////////////////////////////////////
// 获取指定列的向量
//
// 参数：
// 1. int nCols - 指定的矩阵列数
// 2. double* pVector - 指向向量中各元素的缓冲区
//
// 返回值：int 型，向量中元素的个数，即矩阵的行数
//////////////////////////////////////////////////////////////////////
int CMatrix::GetColVector(int nCol, double* pVector) const
{
if (pVector == NULL)
delete pVector;
pVector = new double;
ASSERT(pVector != NULL);
for (int i=0; i<m_nNumRows; ++i)
pVector = GetElement(i, nCol);
return m_nNumRows;
}
//////////////////////////////////////////////////////////////////////
// 重载运算符=，给矩阵赋值
//
// 参数：
// 1. const CMatrix& other - 用于给矩阵赋值的源矩阵
//
// 返回值：CMatrix型的引用，所引用的矩阵与other相等
//////////////////////////////////////////////////////////////////////
CMatrix& CMatrix::operator=(const CMatrix& other)
{
if (&other != this)
{
BOOL bSuccess = Init(other.GetNumRows(), other.GetNumColumns());
ASSERT(bSuccess);
// copy the pointer
memcpy(m_pData, other.m_pData, sizeof(double)*m_nNumColumns*m_nNumRows);
}
// finally return a reference to ourselves
return *this ;
}
//////////////////////////////////////////////////////////////////////
// 重载运算符==，判断矩阵是否相等
//
// 参数：
// 1. const CMatrix& other - 用于比较的矩阵
//
// 返回值：BOOL 型，两个矩阵相等则为TRUE，否则为FALSE
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::operator==(const CMatrix& other) const
{
// 首先检查行列数是否相等
if (m_nNumColumns != other.GetNumColumns() || m_nNumRows != other.GetNumRows())
return FALSE;
for (int i=0; i<m_nNumRows; ++i)
{
for (int j=0; j<m_nNumColumns; ++j)
{
if (GetElement(i, j) != other.GetElement(i, j))
return FALSE;
}
}
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:07 编辑 ]
回复：（风花雪月）[分享]和矩阵运算相关的类代码
//////////////////////////////////////////////////////////////////////// 重载运算符!=，判断矩阵是否不相等
//
// 参数：
// 1. const CMatrix& other - 用于比较的矩阵
//
// 返回值：BOOL 型，两个不矩阵相等则为TRUE，否则为FALSE
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::operator!=(const CMatrix& other) const
{
return !(*this == other);
}
//////////////////////////////////////////////////////////////////////
// 重载运算符+，实现矩阵的加法
//
// 参数：
// 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) ; // 拷贝构造
// 矩阵加法
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. 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()) ;
// 矩阵乘法，即
//
//
// * =
//
//
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 ;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:07 编辑 ]
回复：（风花雪月）[分享]和矩阵运算相关的类代码
//////////////////////////////////////////////////////////////////////// 复矩阵的乘法
//
// 参数：
// 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
{
for (int j=0; j
{
double vr = 0;
double vi = 0;
for (int k =0; 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;
pnCol = new int;
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);
if (p>d)
{
d=p;
pnRow=i;
pnCol=j;
}
}
}
// 失败
if (d == 0.0)
{
delete[] pnRow;
delete[] pnCol;
return FALSE;
}
if (pnRow != k)
{
for (j=0; j<=m_nNumColumns-1; j++)
{
u=k*m_nNumColumns+j;
v=pnRow*m_nNumColumns+j;
p=m_pData;
m_pData=m_pData;
m_pData=p;
}
}
if (pnCol != k)
{
for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+k;
v=i*m_nNumColumns+pnCol;
p=m_pData;
m_pData=m_pData;
m_pData=p;
}
}
l=k*m_nNumColumns+k;
m_pData=1.0/m_pData;
for (j=0; j<=m_nNumColumns-1; j++)
{
if (j != k)
{
u=k*m_nNumColumns+j;
m_pData=m_pData*m_pData;
}
}
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=m_pData-m_pData*m_pData;
}
}
}
}
for (i=0; i<=m_nNumColumns-1; i++)
{
if (i!=k)
{
u=i*m_nNumColumns+k;
m_pData=-m_pData*m_pData;
}
}
}
// 调整恢复行列次序
for (k=m_nNumColumns-1; k>=0; k--)
{
if (pnCol!=k)
{
for (j=0; j<=m_nNumColumns-1; j++)
{
u=k*m_nNumColumns+j;
v=pnCol*m_nNumColumns+j;
p=m_pData;
m_pData=m_pData;
m_pData=p;
}
}
if (pnRow!=k)
{
for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+k;
v=i*m_nNumColumns+pnRow;
p=m_pData;
m_pData=m_pData;
m_pData=p;
}
}
}
// 清理内存
delete[] pnRow;
delete[] pnCol;
// 成功返回
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:09 编辑 ]
回复：（风花雪月）[分享]和矩阵运算相关的类代码
//////////////////////////////////////////////////////////////////////// 复矩阵求逆的全选主元高斯－约当法
//
// 参数：
// 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;
pnCol = new int;
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*m_pData+mtxImag.m_pData*mtxImag.m_pData;
if (p>d)
{
d=p;
pnRow=i;
pnCol=j;
}
}
}
// 失败
if (d == 0.0)
{
delete[] pnRow;
delete[] pnCol;
return(0);
}
if (pnRow!=k)
{
for (j=0; j<=m_nNumColumns-1; j++)
{
u=k*m_nNumColumns+j;
v=pnRow*m_nNumColumns+j;
t=m_pData;
m_pData=m_pData;
m_pData=t;
t=mtxImag.m_pData;
mtxImag.m_pData=mtxImag.m_pData;
mtxImag.m_pData=t;
}
}
if (pnCol!=k)
{
for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+k;
v=i*m_nNumColumns+pnCol;
t=m_pData;
m_pData=m_pData;
m_pData=t;
t=mtxImag.m_pData;
mtxImag.m_pData=mtxImag.m_pData;
mtxImag.m_pData=t;
}
}
l=k*m_nNumColumns+k;
m_pData=m_pData/d; mtxImag.m_pData=-mtxImag.m_pData/d;
for (j=0; j<=m_nNumColumns-1; j++)
{
if (j!=k)
{
u=k*m_nNumColumns+j;
p=m_pData*m_pData;
q=mtxImag.m_pData*mtxImag.m_pData;
s=(m_pData+mtxImag.m_pData)*(m_pData+mtxImag.m_pData);
m_pData=p-q;
mtxImag.m_pData=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*m_pData;
q=mtxImag.m_pData*mtxImag.m_pData;
s=(m_pData+mtxImag.m_pData)*(m_pData+mtxImag.m_pData);
t=p-q;
b=s-p-q;
m_pData=m_pData-t;
mtxImag.m_pData=mtxImag.m_pData-b;
}
}
}
}
for (i=0; i<=m_nNumColumns-1; i++)
{
if (i!=k)
{
u=i*m_nNumColumns+k;
p=m_pData*m_pData;
q=mtxImag.m_pData*mtxImag.m_pData;
s=(m_pData+mtxImag.m_pData)*(m_pData+mtxImag.m_pData);
m_pData=q-p;
mtxImag.m_pData=p+q-s;
}
}
}
// 调整恢复行列次序
for (k=m_nNumColumns-1; k>=0; k--)
{
if (pnCol!=k)
{
for (j=0; j<=m_nNumColumns-1; j++)
{
u=k*m_nNumColumns+j;
v=pnCol*m_nNumColumns+j;
t=m_pData;
m_pData=m_pData;
m_pData=t;
t=mtxImag.m_pData;
mtxImag.m_pData=mtxImag.m_pData;
mtxImag.m_pData=t;
}
}
if (pnRow!=k)
{
for (i=0; i<=m_nNumColumns-1; i++)
{
u=i*m_nNumColumns+k;
v=i*m_nNumColumns+pnRow;
t=m_pData;
m_pData=m_pData;
m_pData=t;
t=mtxImag.m_pData;
mtxImag.m_pData=mtxImag.m_pData;
mtxImag.m_pData=t;
}
}
}
// 清理内存
delete[] pnRow;
delete[] pnCol;
// 成功返回
return TRUE;
}
[ 本帖最后由 风花雪月 于 2008-10-28 17:11 编辑 ]
回复：（风花雪月）[分享]和矩阵运算相关的类代码
Matrix.h部分代码://
class CMatrix
{
//
// 公有接口函数
//
public:
//
// 构造与析构
//
CMatrix(); // 基础构造函数
CMatrix(int nRows, int nCols); // 指定行列构造函数
CMatrix(int nRows, int nCols, double value[]); // 指定数据构造函数
CMatrix(int nSize); // 方阵构造函数
CMatrix(int nSize, double value[]); // 指定数据方阵构造函数
CMatrix(const CMatrix& other); // 拷贝构造函数
BOOL Init(int nRows, int nCols); // 初始化矩阵
BOOL MakeUnitMatrix(int nSize); // 将方阵初始化为单位矩阵
virtual ~CMatrix(); // 析构函数
//
// 输入与显示
//
// 将字符串转换为矩阵数据
BOOL FromString(CString s, const CString& sDelim = " ", BOOL bLineBreak = TRUE);
// 将矩阵转换为字符串
CString ToString(const CString& sDelim = " ", BOOL bLineBreak = TRUE) const;
// 将矩阵的指定行转换为字符串
CString RowToString(int nRow, const CString& sDelim = " ") const;
// 将矩阵的指定列转换为字符串
CString ColToString(int nCol, const CString& sDelim = " ") const;
//
// 元素与值操作
//
BOOL SetElement(int nRow, int nCol, double value); // 设置指定元素的值
double GetElement(int nRow, int nCol) const; // 获取指定元素的值
void SetData(double value[]); // 设置矩阵的值
int GetNumColumns() const; // 获取矩阵的列数
int GetNumRows() const; // 获取矩阵的行数
int GetRowVector(int nRow, double* pVector) const; // 获取矩阵的指定行矩阵
int GetColVector(int nCol, double* pVector) const; // 获取矩阵的指定列矩阵
double* GetData() const; // 获取矩阵的值
//
// 数学操作
//
CMatrix& operator=(const CMatrix& other);
BOOL operator==(const CMatrix& other) const;
BOOL operator!=(const CMatrix& other) const;
CMatrix operator+(const CMatrix& other) const;
CMatrix operator-(const CMatrix& other) const;
CMatrix operator*(double value) const;
CMatrix operator*(const CMatrix& other) const;
// 复矩阵乘法
BOOL CMul(const CMatrix& AR, const CMatrix& AI, const
CMatrix& BR, const CMatrix& BI, CMatrix& CR, CMatrix&
CI) const;
// 矩阵的转置
CMatrix Transpose() const;
//
// 算法
//
// 实矩阵求逆的全选主元高斯－约当法
BOOL InvertGaussJordan();
// 复矩阵求逆的全选主元高斯－约当法
BOOL InvertGaussJordan(CMatrix& mtxImag);
// 对称正定矩阵的求逆
BOOL InvertSsgj();
// 托伯利兹矩阵求逆的埃兰特方法
BOOL InvertTrench();
// 求行列式值的全选主元高斯消去法
double DetGauss();
// 求矩阵秩的全选主元高斯消去法
int RankGauss();
// 对称正定矩阵的乔里斯基分解与行列式的求值
BOOL DetCholesky(double* dblDet);
// 矩阵的三角分解
BOOL SplitLU(CMatrix& mtxL, CMatrix& mtxU);
// 一般实矩阵的QR分解
BOOL SplitQR(CMatrix& mtxQ);
// 一般实矩阵的奇异值分解
BOOL SplitUV(CMatrix& mtxU, CMatrix& mtxV, double eps = 0.000001);
// 求广义逆的奇异值分解法
BOOL GInvertUV(CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double eps = 0.000001);
// 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
BOOL MakeSymTri(CMatrix& mtxQ, CMatrix& mtxT, double dblB[], double dblC[]);
// 实对称三对角阵的全部特征值与特征向量的计算
BOOL SymTriEigenv(double dblB[], double dblC[], CMatrix& mtxQ, int nMaxIt = 60, double eps = 0.000001);
// 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
void MakeHberg();
// 求赫申伯格矩阵全部特征值的QR方法
BOOL HBergEigenv(double dblU[], double dblV[], int nMaxIt = 60, double eps = 0.000001);
// 求实对称矩阵特征值与特征向量的雅可比法
BOOL JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt = 60, double eps = 0.000001);
// 求实对称矩阵特征值与特征向量的雅可比过关法
BOOL JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, double eps = 0.000001);
//
// 保护性数据成员
//
protected:
int m_nNumColumns; // 矩阵列数
int m_nNumRows; // 矩阵行数
double* m_pData; // 矩阵数据缓冲区
//
// 内部函数
//
private:
void ppp(double a[], double e[], double s[], double v[], int m, int n);
void sss(double fg, double cs);
};
[ 本帖最后由 风花雪月 于 2008-10-28 17:11 编辑 ] 顶起来，谢谢，学习 我想求广义逆，不知好不好用，谢谢你。 谢谢了，我也学习了 谢谢，学习:lol :lol :lol
页:
[1]