h12.inl

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

} 

//////////////////////////////////////////////////////////////////////
// 内部函数，由SplitUV函数调用
//////////////////////////////////////////////////////////////////////
void CMatrix::ppp(double a[], double e[], double s[], double v[], int m, int n)
{ 
 int i,j,p,q;
    double d; 

    if (m>=n) 
  i=n;
    else 
  i=m; 

    for (j=1; j<=i-1; j++)
    { 
  a[(j-1)*n+j-1]=s[j-1];
        a[(j-1)*n+j]=e[j-1];
    }
    
 a[(i-1)*n+i-1]=s[i-1];
    if (m<n) 
  a[(i-1)*n+i]=e[i-1];
    
 for (i=1; i<=n-1; i++)
 {
  for (j=i+1; j<=n; j++)
  { 
   p=(i-1)*n+j-1; 
   q=(j-1)*n+i-1;
   d=v[p]; 
   v[p]=v[q]; 
   v[q]=d;
  }
 }
} 

//////////////////////////////////////////////////////////////////////
// 内部函数，由SplitUV函数调用
//////////////////////////////////////////////////////////////////////
void CMatrix::sss(double fg[2], double cs[2])
{ 
 double r,d;
    
 if ((fabs(fg[0])+fabs(fg[1]))==0.0)
    { 
  cs[0]=1.0; 
  cs[1]=0.0; 
  d=0.0;
 }
    else 
    { 
  d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
        if (fabs(fg[0])>fabs(fg[1]))
        { 
   d=fabs(d);
            if (fg[0]<0.0) 
    d=-d;
        }
        if (fabs(fg[1])>=fabs(fg[0]))
        { 
   d=fabs(d);
            if (fg[1]<0.0) 
    d=-d;
        }
        
  cs[0]=fg[0]/d; 
  cs[1]=fg[1]/d;
    }
    
 r=1.0;
    if (fabs(fg[0])>fabs(fg[1])) 
  r=cs[1];
    else if (cs[0]!=0.0) 
  r=1.0/cs[0]; 

    fg[0]=d; 
 fg[1]=r;
} 

//////////////////////////////////////////////////////////////////////
// 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值
//
// 参数：
// 1. CMatrix& mtxAP - 返回原矩阵的广义逆矩阵
// 2. CMatrix& mtxU - 返回分解后的U矩阵
// 3. CMatrix& mtxV - 返回分解后的V矩阵
// 4. double eps - 计算精度，默认值为0.000001
//
// 返回值：bool型，求解是否成功
//////////////////////////////////////////////////////////////////////
bool CMatrix::GInvertUV(CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/)
{ 
 int i,j,k,l,t,p,q,f; 

 // 调用奇异值分解
    if (! SplitUV(mtxU, mtxV, eps))
  return false; 

 int m = m_nNumRows;
 int n = m_nNumColumns; 

 // 初始化广义逆矩阵
 if (! mtxAP.Init(n, m))
  return false; 

 // 计算广义逆矩阵 

    j=n;
    if (m<n) 
  j=m;
    j=j-1;
    k=0;
    while ((k<=j)&&(m_pData[k*n+k]!=0.0)) 
  k=k+1; 

    k=k-1;
    for (i=0; i<=n-1; i++)
 {
  for (j=0; j<=m-1; j++)
  { 
   t=i*m+j; 
   mtxAP.m_pData[t]=0.0;
   for (l=0; l<=k; l++)
   { 
    f=l*n+i; 
    p=j*m+l; 
    q=l*n+l;
    mtxAP.m_pData[t]=mtxAP.m_pData[t]+mtxV.m_pData[f]*mtxU.m_pData[p]/m_pData[q];
   }
  }
 } 

    return true;
} 

//////////////////////////////////////////////////////////////////////
// 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
//
// 参数：
// 1. CMatrix& mtxQ - 返回豪斯荷尔德变换的乘积矩阵Q
// 2. CMatrix& mtxT - 返回求得的对称三对角阵
// 3. double dblB[] - 一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元素
// 4. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的次对角线元素
//
// 返回值：bool型，求解是否成功
//////////////////////////////////////////////////////////////////////
bool CMatrix::MakeSymTri(CMatrix& mtxQ, CMatrix& mtxT, double dblB[], double dblC[])
{ 
 int i,j,k,u;
    double h,f,g,h2;
    
 // 初始化矩阵Q和T
 if (! mtxQ.Init(m_nNumColumns, m_nNumColumns) ||
  ! mtxT.Init(m_nNumColumns, m_nNumColumns))
  return false; 

 if (dblB == NULL || dblC == NULL)
  return false; 

 for (i=0; i<=m_nNumColumns-1; i++)
 {
  for (j=0; j<=m_nNumColumns-1; j++)
  { 
   u=i*m_nNumColumns+j; 
   mtxQ.m_pData[u]=m_pData[u];
  }
 } 

    for (i=m_nNumColumns-1; i>=1; i--)
    { 
  h=0.0;
        if (i>1)
  {
   for (k=0; k<=i-1; k++)
            { 
    u=i*m_nNumColumns+k; 
    h=h+mtxQ.m_pData[u]*mtxQ.m_pData[u];
   }
  } 

        if (h == 0.0)
        { 
   dblC[i]=0.0;
            if (i==1) 
    dblC[i]=mtxQ.m_pData[i*m_nNumColumns+i-1];
            dblB[i]=0.0;
        }
        else
        { 
   dblC[i]=sqrt(h);
            u=i*m_nNumColumns+i-1;
            if (mtxQ.m_pData[u]>0.0) 
    dblC[i]=-dblC[i]; 

            h=h-mtxQ.m_pData[u]*dblC[i];
            mtxQ.m_pData[u]=mtxQ.m_pData[u]-dblC[i];
            f=0.0;
            for (j=0; j<=i-1; j++)
            { 
    mtxQ.m_pData[j*m_nNumColumns+i]=mtxQ.m_pData[i*m_nNumColumns+j]/h;
                g=0.0;
                for (k=0; k<=j; k++)
     g=g+mtxQ.m_pData[j*m_nNumColumns+k]*mtxQ.m_pData[i*m_nNumColumns+k]; 

    if (j+1<=i-1)
     for (k=j+1; k<=i-1; k++)
      g=g+mtxQ.m_pData[k*m_nNumColumns+j]*mtxQ.m_pData[i*m_nNumColumns+k]; 

                dblC[j]=g/h;
                f=f+g*mtxQ.m_pData[j*m_nNumColumns+i];
            }
            
   h2=f/(h+h);
            for (j=0; j<=i-1; j++)
            { 
    f=mtxQ.m_pData[i*m_nNumColumns+j];
                g=dblC[j]-h2*f;
                dblC[j]=g;
                for (k=0; k<=j; k++)
                { 
     u=j*m_nNumColumns+k;
                    mtxQ.m_pData[u]=mtxQ.m_pData[u]-f*dblC[k]-g*mtxQ.m_pData[i*m_nNumColumns+k];
                }
            }
            
   dblB[i]=h;
        }
    }
    
 for (i=0; i<=m_nNumColumns-2; i++) 
  dblC[i]=dblC[i+1];
    
 dblC[m_nNumColumns-1]=0.0;
    dblB[0]=0.0;
    for (i=0; i<=m_nNumColumns-1; i++)
    { 
  if ((dblB[i]!=0.0)&&(i-1>=0))
  {
   for (j=0; j<=i-1; j++)
            { 
    g=0.0;
    for (k=0; k<=i-1; k++)
     g=g+mtxQ.m_pData[i*m_nNumColumns+k]*mtxQ.m_pData[k*m_nNumColumns+j]; 

    for (k=0; k<=i-1; k++)
                { 
     u=k*m_nNumColumns+j;
     mtxQ.m_pData[u]=mtxQ.m_pData[u]-g*mtxQ.m_pData[k*m_nNumColumns+i];
                }
            }
  } 

        u=i*m_nNumColumns+i;
        dblB[i]=mtxQ.m_pData[u]; mtxQ.m_pData[u]=1.0;
        if (i-1>=0)
  {
   for (j=0; j<=i-1; j++)
            { 
    mtxQ.m_pData[i*m_nNumColumns+j]=0.0; 
    mtxQ.m_pData[j*m_nNumColumns+i]=0.0;
   }
  }
    } 

    // 构造对称三对角矩阵
    for (i=0; i<m_nNumColumns; ++i)
 {
     for (j=0; j<m_nNumColumns; ++j)
  {
            mtxT.SetElement(i, j, 0);
            k = i - j;
            if (k == 0) 
             mtxT.SetElement(i, j, dblB[j]);
   else if (k == 1)
             mtxT.SetElement(i, j, dblC[j]);
   else if (k == -1)
             mtxT.SetElement(i, j, dblC[i]);
        }
    } 

 return true;
} 

//////////////////////////////////////////////////////////////////////
// 实对称三对角阵的全部特征值与特征向量的计算
//
// 参数：
// 1. double dblB[] - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
//    返回时存放全部特征值。
// 2. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素
// 3. CMatrix& mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
//    如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积矩阵Q，则返回矩阵A的
//    特征值向量矩阵。其中第i列为与数组dblB中第j个特征值对应的特征向量。
// 4. int nMaxIt - 迭代次数，默认值为60
// 5. double eps - 计算精度，默认值为0.000001
//
// 返回值：bool型，求解是否成功
//////////////////////////////////////////////////////////////////////
bool CMatrix::SymTriEigenv(double dblB[], double dblC[], CMatrix& mtxQ, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
 int i,j,k,m,it,u,v;
    double d,f,h,g,p,r,e,s;
    
 // 初值
 int n = mtxQ.GetNumColumns();
 dblC[n-1]=0.0; 
 d=0.0; 
 f=0.0;
    
 // 迭代计算 

 for (j=0; j<=n-1; j++)
    { 
  it=0;
        h=eps*(fabs(dblB[j])+fabs(dblC[j]));
        if (h>d) 
   d=h;
        
  m=j;
        while ((m<=n-1)&&(fabs(dblC[m])>d)) 
   m=m+1;
        
  if (m!=j)
        { 
   do
            { 
    if (it==nMaxIt)
     return false; 

                it=it+1;
                g=dblB[j];
                p=(dblB[j+1]-g)/(2.0*dblC[j]);
                r=sqrt(p*p+1.0);
                if (p>=0.0) 
     dblB[j]=dblC[j]/(p+r);
                else 
     dblB[j]=dblC[j]/(p-r);
                
    h=g-dblB[j];
                for (i=j+1; i<=n-1; i++)
     dblB[i]=dblB[i]-h;
                
    f=f+h; 
    p=dblB[m]; 
    e=1.0; 
    s=0.0;
                for (i=m-1; i>=j; i--)
                { 
     g=e*dblC[i]; 
     h=e*p;
                    if (fabs(p)>=fabs(dblC[i]))
                    { 
      e=dblC[i]/p; 
      r=sqrt(e*e+1.0);
                        dblC[i+1]=s*p*r; 
      s=e/r; 
      e=1.0/r;
                    }
                    else
     { 
      e=p/dblC[i]; 
      r=sqrt(e*e+1.0);
                        dblC[i+1]=s*dblC[i]*r;
                        s=1.0/r; 
      e=e/r;
                    }
                    
     p=e*dblB[i]-s*g;
                    dblB[i+1]=h+s*(e*g+s*dblB[i]);
                    for (k=0; k<=n-1; k++)
                    { 
      u=k*n+i+1; 
      v=u-1;
                        h=mtxQ.m_pData[u]; 
      mtxQ.m_pData[u]=s*mtxQ.m_pData[v]+e*h;
                        mtxQ.m_pData[v]=e*mtxQ.m_pData[v]-s*h;
                    }
                }
                
    dblC[j]=s*p; 
    dblB[j]=e*p;
            
   } while (fabs(dblC[j])>d);
        }
        
  dblB[j]=dblB[j]+f;
    }
    
 for (i=0; i<=n-1; i++)
    { 
  k=i; 
  p=dblB[i];
        if (i+1<=n-1)
        { 
   j=i+1;
            while ((j<=n-1)&&(dblB[j]<=p))
            { 
    k=j; 
    p=dblB[j]; 
    j=j+1;
   }
        } 

        if (k!=i)
        { 
   dblB[k]=dblB[i]; 
   dblB[i]=p;
            for (j=0; j<=n-1; j++)
            { 
    u=j*n+i; 
    v=j*n+k;
                p=mtxQ.m_pData[u]; 
    mtxQ.m_pData[u]=mtxQ.m_pData[v]; 
    mtxQ.m_pData[v]=p;
            }
        }
    }
    
 return true;
} 

//////////////////////////////////////////////////////////////////////
// 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
//
// 参数：无
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CMatrix::MakeHberg()
{ 
 int i,j,k,u,v;
    double d,t; 

    for (k=1; k<=m_nNumColumns-2; k++)
    { 
  d=0.0;
        for (j=k; j<=m_nNumColumns-1; j++)
        { 
   u=j*m_nNumColumns+k-1; 
   t=m_pData[u];
            if (fabs(t)>fabs(d))
            { 
    d=t; 
    i=j;
   }
        }
        
  if (d != 0.0)
        { 
   if (i!=k)
            { 
    for (j=k-1; j<=m_nNumColumns-1; j++)
                { 
     u=i*m_nNumColumns+j; 
     v=k*m_nNumColumns+j;
                    t=m_pData[u]; 
     m_pData[u]=m_pData[v]; 
     m_pData[v]=t;
                }
                
    for (j=0; j<=m_nNumColumns-1; j++)
                { 
     u=j*m_nNumColumns+i; 
     v=j*m_nNumColumns+k;
                    t=m_pData[u]; 
     m_pData[u]=m_pData[v]; 
     m_pData[v]=t;
                }
            }
            
   for (i=k+1; i<=m_nNumColumns-1; i++)
            { 
    u=i*m_nNumColumns+k-1; 
    t=m_pData[u]/d; 
    m_pData[u]=0.0;
                for (j=k; j<=m_nNumColumns-1; j++)
                { 
     v=i*m_nNumColumns+j;
                    m_pData[v]=m_pData[v]-t*m_pData[k*m_nNumColumns+j];
                }
                
    for (j=0; j<=m_nNumColumns-1; j++)
                { 
     v=j*m_nNumColumns+k;
                    m_pData[v]=m_pData[v]+t*m_pData[j*m_nNumColumns+i];
                }
            }
        }
    }
}

//////////////////////////////////////////////////////////////////////
// 求赫申伯格矩阵全部特征值的QR方法
//
// 参数：
// 1. double dblU[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的实部
// 2. double dblV[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的虚部
// 3. int nMaxIt - 迭代次数，默认值为60
// 4. double eps - 计算精度，默认值为0.000001
//
// 返回值：bool型，求解是否成功
////////////////////////////