h12.inl

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

 int n = m_nNumColumns; 

 // 初始化U, V矩阵
 if (! mtxU.Init(m, m) || ! mtxV.Init(n, n))
  return false; 

 // 临时缓冲区
 int ka = __max(m, n) + 1;
    s = new double[ka];
    e = new double[ka];
    w = new double[ka]; 

 // 指定迭代次数为60
    it=60; 
 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>=kk+1)
            { 
    for (j=kk+1; j<=n; j++)
                { 
     if ((kk<=k)&&(s[kk-1]!=0.0))
                    { 
      d=0.0;
                        for (i=kk; i<=m; i++)
                        { 
       ix=(i-1)*n+kk-1;
                            iy=(i-1)*n+j-1;
                            d=d+m_pData[ix]*m_pData[iy];
                        }
                        
      d=-d/m_pData[(kk-1)*n+kk-1];
                        for (i=kk; i<=m; i++)
                        { 
       ix=(i-1)*n+j-1;
                            iy=(i-1)*n+kk-1;
                            m_pData[ix]=m_pData[ix]+d*m_pData[iy];
                        }
                    }
                    
     e[j-1]=m_pData[(kk-1)*n+j-1];
                }
            }
            
   if (kk<=k)
            { 
    for (i=kk; i<=m; i++)
                { 
     ix=(i-1)*m+kk-1; 
     iy=(i-1)*n+kk-1;
                    mtxU.m_pData[ix]=m_pData[iy];
                }
            }
            
   if (kk<=l)
            { 
    d=0.0;
                for (i=kk+1; i<=n; i++)
     d=d+e[i-1]*e[i-1];
                
    e[kk-1]=sqrt(d);
                if (e[kk-1]!=0.0)
                { 
     if (e[kk]!=0.0)
                    { 
      e[kk-1]=fabs(e[kk-1]);
                        if (e[kk]<0.0) 
       e[kk-1]=-e[kk-1];
                    } 

                    for (i=kk+1; i<=n; i++)
                      e[i-1]=e[i-1]/e[kk-1];
                    
     e[kk]=1.0+e[kk];
                }
                
    e[kk-1]=-e[kk-1];
                if ((kk+1<=m)&&(e[kk-1]!=0.0))
                { 
     for (i=kk+1; i<=m; i++) 
      w[i-1]=0.0;
                    
     for (j=kk+1; j<=n; j++)
      for (i=kk+1; i<=m; i++)
       w[i-1]=w[i-1]+e[j-1]*m_pData[(i-1)*n+j-1];
                    
     for (j=kk+1; j<=n; j++)
     {
      for (i=kk+1; i<=m; i++)
                        { 
       ix=(i-1)*n+j-1;
       m_pData[ix]=m_pData[ix]-w[i-1]*e[j-1]/e[kk];
                        }
     }
                }
                
    for (i=kk+1; i<=n; i++)
                  mtxV.m_pData[(i-1)*n+kk-1]=e[i-1];
            }
        }
    }
    
 mm=n;
    if (m+1<n) 
  mm=m+1;
    if (k<n) 
  s[k]=m_pData[k*n+k];
    if (m<mm) 
  s[mm-1]=0.0;
    if (l+1<mm) 
  e[l]=m_pData[l*n+mm-1]; 

    e[mm-1]=0.0;
    nn=m;
    if (m>n) 
  nn=n;
    if (nn>=k+1)
    { 
  for (j=k+1; j<=nn; j++)
        { 
   for (i=1; i<=m; i++)
    mtxU.m_pData[(i-1)*m+j-1]=0.0;
            mtxU.m_pData[(j-1)*m+j-1]=1.0;
        }
    }
    
 if (k>=1)
    { 
  for (ll=1; ll<=k; ll++)
        { 
   kk=k-ll+1; 
   iz=(kk-1)*m+kk-1;
            if (s[kk-1]!=0.0)
            { 
    if (nn>=kk+1)
    {
     for (j=kk+1; j<=nn; j++)
     { 
      d=0.0;
      for (i=kk; i<=m; i++)
      { 
       ix=(i-1)*m+kk-1;
       iy=(i-1)*m+j-1;
       d=d+mtxU.m_pData[ix]*mtxU.m_pData[iy]/mtxU.m_pData[iz];
      } 

      d=-d;
      for (i=kk; i<=m; i++)
      { 
       ix=(i-1)*m+j-1;
       iy=(i-1)*m+kk-1;
       mtxU.m_pData[ix]=mtxU.m_pData[ix]+d*mtxU.m_pData[iy];
      }
     }
    }
                  
    for (i=kk; i<=m; i++)
    { 
     ix=(i-1)*m+kk-1; 
     mtxU.m_pData[ix]=-mtxU.m_pData[ix];
    } 

    mtxU.m_pData[iz]=1.0+mtxU.m_pData[iz];
    if (kk-1>=1)
    {
     for (i=1; i<=kk-1; i++)
      mtxU.m_pData[(i-1)*m+kk-1]=0.0;
    }
   }
            else
            { 
    for (i=1; i<=m; i++)
     mtxU.m_pData[(i-1)*m+kk-1]=0.0;
                mtxU.m_pData[(kk-1)*m+kk-1]=1.0;
            }
  }
    } 

    for (ll=1; ll<=n; ll++)
    { 
  kk=n-ll+1; 
  iz=kk*n+kk-1;
        
  if ((kk<=l)&&(e[kk-1]!=0.0))
        { 
   for (j=kk+1; j<=n; j++)
            { 
    d=0.0;
                for (i=kk+1; i<=n; i++)
                { 
     ix=(i-1)*n+kk-1; 
     iy=(i-1)*n+j-1;
                    d=d+mtxV.m_pData[ix]*mtxV.m_pData[iy]/mtxV.m_pData[iz];
                }
                
    d=-d;
                for (i=kk+1; i<=n; i++)
                { 
     ix=(i-1)*n+j-1; 
     iy=(i-1)*n+kk-1;
                    mtxV.m_pData[ix]=mtxV.m_pData[ix]+d*mtxV.m_pData[iy];
                }
            }
        }
        
  for (i=1; i<=n; i++)
   mtxV.m_pData[(i-1)*n+kk-1]=0.0;
        
  mtxV.m_pData[iz-n]=1.0;
    }
    
 for (i=1; i<=m; i++)
  for (j=1; j<=n; j++)
   m_pData[(i-1)*n+j-1]=0.0;
    
 m1=mm; 
 it=60;
    while (true)
    { 
  if (mm==0)
        { 
   ppp(m_pData,e,s,mtxV.m_pData,m,n);
            return true;
        }
        if (it==0)
        { 
   ppp(m_pData,e,s,mtxV.m_pData,m,n);
            return false;
        }
        
  kk=mm-1;
  while ((kk!=0)&&(fabs(e[kk-1])!=0.0))
        { 
   d=fabs(s[kk-1])+fabs(s[kk]);
            dd=fabs(e[kk-1]);
            if (dd>eps*d) 
    kk=kk-1;
            else 
    e[kk-1]=0.0;
        }
        
  if (kk==mm-1)
        { 
   kk=kk+1;
            if (s[kk-1]<0.0)
            { 
    s[kk-1]=-s[kk-1];
                for (i=1; i<=n; i++)
                { 
     ix=(i-1)*n+kk-1; 
     mtxV.m_pData[ix]=-mtxV.m_pData[ix];}
    }
    
    while ((kk!=m1)&&(s[kk-1]<s[kk]))
    { 
     d=s[kk-1]; 
     s[kk-1]=s[kk]; 
     s[kk]=d;
     if (kk<n)
     {
      for (i=1; i<=n; i++)
      { 
       ix=(i-1)*n+kk-1; 
       iy=(i-1)*n+kk;
       d=mtxV.m_pData[ix]; 
       mtxV.m_pData[ix]=mtxV.m_pData[iy]; 
       mtxV.m_pData[iy]=d;
      }
     } 

     if (kk<m)
     {
      for (i=1; i<=m; i++)
      { 
       ix=(i-1)*m+kk-1; 
       iy=(i-1)*m+kk;
       d=mtxU.m_pData[ix]; 
       mtxU.m_pData[ix]=mtxU.m_pData[iy]; 
       mtxU.m_pData[iy]=d;
      }
     } 

     kk=kk+1;
            }
            
   it=60;
            mm=mm-1;
        }
        else
        { 
   ks=mm;
            while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
            { 
    d=0.0;
                if (ks!=mm) 
     d=d+fabs(e[ks-1]);
                if (ks!=kk+1) 
     d=d+fabs(e[ks-2]);
                
    dd=fabs(s[ks-1]);
                if (dd>eps*d) 
     ks=ks-1;
                else 
     s[ks-1]=0.0;
            }
            
   if (ks==kk)
            { 
    kk=kk+1;
                d=fabs(s[mm-1]);
                t=fabs(s[mm-2]);
                if (t>d) 
     d=t;
                
    t=fabs(e[mm-2]);
                if (t>d) 
     d=t;
                
    t=fabs(s[kk-1]);
                if (t>d) 
     d=t;
                
    t=fabs(e[kk-1]);
                if (t>d) 
     d=t;
                
    sm=s[mm-1]/d; 
    sm1=s[mm-2]/d;
                em1=e[mm-2]/d;
                sk=s[kk-1]/d; 
    ek=e[kk-1]/d;
                b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
                c=sm*em1; 
    c=c*c; 
    shh=0.0; 

                if ((b!=0.0)||(c!=0.0))
                { 
     shh=sqrt(b*b+c);
                    if (b<0.0) 
      shh=-shh; 

                    shh=c/(b+shh);
                }
                
    fg[0]=(sk+sm)*(sk-sm)-shh;
                fg[1]=sk*ek;
                for (i=kk; i<=mm-1; i++)
                { 
     sss(fg,cs);
                    if (i!=kk) 
      e[i-2]=fg[0]; 

                    fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
                    e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
                    fg[1]=cs[1]*s[i];
                    s[i]=cs[0]*s[i]; 

                    if ((cs[0]!=1.0)||(cs[1]!=0.0))
     {
      for (j=1; j<=n; j++)
                        { 
       ix=(j-1)*n+i-1;
       iy=(j-1)*n+i;
       d=cs[0]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy];
       mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy];
       mtxV.m_pData[ix]=d;
                        }
     } 

                    sss(fg,cs);
                    s[i-1]=fg[0];
                    fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
                    s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
                    fg[1]=cs[1]*e[i];
                    e[i]=cs[0]*e[i]; 

                    if (i<m)
     {
      if ((cs[0]!=1.0)||(cs[1]!=0.0))
      {
       for (j=1; j<=m; j++)
       { 
        ix=(j-1)*m+i-1;
        iy=(j-1)*m+i;
        d=cs[0]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy];
        mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy];
        mtxU.m_pData[ix]=d;
       }
      }
     }
                }
                
    e[mm-2]=fg[0];
                it=it-1;
            }
            else
            { 
    if (ks==mm)
                { 
     kk=kk+1;
                    fg[1]=e[mm-2]; 
     e[mm-2]=0.0;
                    for (ll=kk; ll<=mm-1; ll++)
                    { 
      i=mm+kk-ll-1;
                        fg[0]=s[i-1];
                        sss(fg,cs);
                        s[i-1]=fg[0];
                        if (i!=kk)
                        { 
       fg[1]=-cs[1]*e[i-2];
                            e[i-2]=cs[0]*e[i-2];
                        }
                        
      if ((cs[0]!=1.0)||(cs[1]!=0.0))
      {
       for (j=1; j<=n; j++)
                            { 
        ix=(j-1)*n+i-1;
        iy=(j-1)*n+mm-1;
        d=cs[0]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy];
        mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy];
        mtxV.m_pData[ix]=d;
                            }
      }
                    }
                }
                else
                { 
     kk=ks+1;
                    fg[1]=e[kk-2];
                    e[kk-2]=0.0;
                    for (i=kk; i<=mm; i++)
                    { 
      fg[0]=s[i-1];
                        sss(fg,cs);
                        s[i-1]=fg[0];
                        fg[1]=-cs[1]*e[i-1];
                        e[i-1]=cs[0]*e[i-1];
                        if ((cs[0]!=1.0)||(cs[1]!=0.0))
      {
       for (j=1; j<=m; j++)
                            { 
        ix=(j-1)*m+i-1;
        iy=(j-1)*m+kk-2;
        d=cs[0]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy];
        mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy];
        mtxU.m_pData[ix]=d;
                            }
      }
                    }
                }
            }
        }
    }
    
 return true;