h12.inl

C环境下矩阵运算类,包含矩阵的基本操作,另有其范数,均值等特殊函数

  }
 }
 return true;
}
//////////////////////////////////////////////////////////////////////
// 重载运算符!=，判断矩阵是否不相等
//
// 参数：
// 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()) ;
 // 矩阵乘法，即
 //
 // [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++)