h12.inl

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

   m_pData[i*m_nNumColumns+j]=m_pData[j*m_nNumColumns+i];
 // 临时内存清理
 delete[] pTmp;
 return true;
}
//////////////////////////////////////////////////////////////////////
// 托伯利兹矩阵求逆的埃兰特方法
//
// 参数：无
//
// 返回值：bool型，求逆是否成功
//////////////////////////////////////////////////////////////////////
bool CMatrix::InvertTrench()
{ 
 int i,j,k;
    double a,s,*t,*tt,*c,*r,*p;
 // 上三角元素
 t = new double[m_nNumColumns];
 // 下三角元素
 tt = new double[m_nNumColumns];
 // 上、下三角元素赋值
 for (i=0; i<m_nNumColumns; ++i)
 {
  t[i] = GetElement(0, i);
     tt[i] = GetElement(i, 0);
 }
 // 临时缓冲区
 c = new double[m_nNumColumns];
 r = new double[m_nNumColumns];
 p = new double[m_nNumColumns];
 // 非Toeplitz矩阵，返回
    if (t[0] == 0.0)
    { 
  delete[] t;
  delete[] tt;
  delete[] c;
  delete[] r;
  delete[] p;
  return false;
    }
    a=t[0]; 
 c[0]=tt[1]/t[0]; 
 r[0]=t[1]/t[0];
    for (k=0; k<=m_nNumColumns-3; k++)
    { 
  s=0.0;
        for (j=1; j<=k+1; j++)
   s=s+c[k+1-j]*tt[j];
        s=(s-tt[k+2])/a;
  for (i=0; i<=k; i++)
   p[i]=c[i]+s*r[k-i];
        c[k+1]=-s;
        s=0.0;
        for (j=1; j<=k+1; j++)
          s=s+r[k+1-j]*t[j];
        
  s=(s-t[k+2])/a;
        for (i=0; i<=k; i++)
        { 
   r[i]=r[i]+s*c[k-i];
            c[k-i]=p[k-i];
        }
        r[k+1]=-s;
  a=0.0;
        for (j=1; j<=k+2; j++)
          a=a+t[j]*c[j-1];
        a=t[0]-a;
  // 求解失败
        if (a == 0.0)
  { 
   delete[] t;
   delete[] tt;
   delete[] c;
   delete[] r;
   delete[] p;
   return false;
  }
    }
    m_pData[0]=1.0/a;
    for (i=0; i<=m_nNumColumns-2; i++)
    { 
  k=i+1; 
  j=(i+1)*m_nNumColumns;
        m_pData[k]=-r[i]/a; 
  m_pData[j]=-c[i]/a;
    }
   for (i=0; i<=m_nNumColumns-2; i++)
 {
  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;