matrix.inl

vc++实现矩阵运算 1高斯-约当法求逆矩阵 2对称正定矩阵的逆矩阵 3托伯利兹矩阵的逆矩阵 4实矩阵的三角(LU)分解 5一般实矩阵的QR分解 6对称正定矩阵的乔里斯基分

        {
			_Ty t=rhs(i,k)/u;
			alpha=alpha+t*t;
		}
        
		if(rhs(k,k)>0.0) u=-u;
        
		alpha=u*sqrt(alpha);
        
		if(FloatEqual(alpha,0.0)) return(0);

        u=sqrt(2.0*alpha*(alpha-rhs(k,k)));

        if(FloatNotEqual(u,0.0))
        { 
			rhs(k,k)=(rhs(k,k)-alpha)/u;
            
			for(i=k+1; i<stRow; i++)
            	rhs(i,k) /= u;

			for(size_t j=0; j<stRow; j++)
            {
				_Ty t=0.0;

                for(size_t jj=k; jj<stRow; jj++)
					t=t+rhs(jj,k)*rhq(jj,j);

                for(i=k; i<stRow; i++)
					rhq(i,j)=rhq(i,j)-2.0*t*rhs(i,k);
            }
            
			for(j=k+1; j<stCol; j++)
            { 
				_Ty t=0.0;
            
				for(size_t jj=k; jj<stRow; jj++)
					t=t+rhs(jj,k)*rhs(jj,j);

                for(i=k; i<stRow; i++)
            		rhs(i,j)=rhs(i,j)-2.0*t*rhs(i,k);
			}

            rhs(k,k)=alpha;

            for(i=k+1; i<stRow; i++)
              rhs(i,k)=0.0;
          }
    }

	for(i=0; i<stRow-1; i++)
		for(size_t j=i+1; j<stRow;j++)
				swap(rhq(i,j), rhq(j,i));
	
    return (1);

}

//对称正定阵的乔里斯基(Cholesky)分解及求其行列式值
//矩阵与返回值类型必须是浮点型
template <class _Ty>
long double MatrixSymmetryRegularCholesky(matrix<_Ty>& rhs)
{
	int iReValue= MatrixSymmetryRegular(rhs, 1);	//判别对称正定
	if(iReValue < 2)
		return long double(0);		//rhs不是对称正定阵

	size_t stRank = rhs.GetColNum();	// 矩阵阶数

	matrix<_Ty> m(rhs);				//生成一matrix对象，用rhs初始化

	long double Det = m(0,0);					//计算行列式值
	m(0,0) = sqrt(m(0,0)); 
	
	for(size_t i=1; i<stRank; i++)
		m(i, 0) /= m(0, 0);
	for(size_t j=1; j<stRank; j++)
	{
		for(size_t k=0; k<j; k++)
			m(j,j) = m(j,j) - m(j,k) * m(j,k);
		Det *= m(j,j);
		m(j,j) = sqrt(m(j,j));
		for(i=j+1; i<stRank; i++)
		{
			for(k=0; k<j; k++)
				m(i,j) = m(i,j) -m(i,k) * m(j,k);
			m(i,j) /= m(j,j);
		}
	}
	for(i=0; i<stRank-1; i++)
	    for(j=i+1; j<stRank; j++)
			m(i,j)=0;

	rhs = m;		//返回Cholesky阵，原矩阵将被复盖
	return Det;		//返回行列式值
}

//一般实矩阵的奇异值分解
template <class _Ty>
int MatrixSingularValue(matrix<_Ty>& a, matrix<_Ty>& u, 
											matrix<_Ty>& v, _Ty eps)
{
	int i, it(60), kk, mm, nn, m1, ks, ka;
    _Ty d,dd,t,sm,sm1,em1,sk,ek,b,c,shh;

	int m = a.GetRowNum();
	int n = a.GetColNum();

	for(int j=0; j<m; j++) u(j, m-1) = _Ty(0);

	if(m > n) ka = m + 1;
	else ka = n + 1;
	
	valarray<_Ty> s(ka), e(ka), w(ka), fg(2), cs(2);
    
	int k = n;
    if(m-1<n) k=m-1;
    int l = m;
    if(n-2<m) l=n-2;
    if(l<0) l=0;
    int 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++)
					d=d+a(i-1,kk-1)*a(i-1,kk-1);
                s[kk-1]=sqrt(d);
                if(FloatNotEqual(s[kk-1],0.0))
                {
                    if(FloatNotEqual(a(kk-1,kk-1),0.0))
                    {
						s[kk-1]=Abs(s[kk-1]);
                        if(a(kk-1,kk-1)<0.0) s[kk-1]=-s[kk-1];
                    }
                    for(i=kk; i<=m; i++)	a(i-1,kk-1)=a(i-1,kk-1)/s[kk-1];
                    a(kk-1,kk-1)=1.0+a(kk-1,kk-1);
                }
                s[kk-1]=-s[kk-1];
            }
            if(n>=kk+1)
            { 
				for(j=kk+1; j<=n; j++)
                {
					if((kk<=k)&&(FloatNotEqual(s[kk-1],0.0)))
                    {
						d=0.0;
                        for(i=kk; i<=m; i++)	d=d+a(i-1,kk-1)*a(i-1,j-1);
                        d=-d/a(kk-1,kk-1);
                        for(i=kk; i<=m; i++)	a(i-1,j-1)=a(i-1,j-1)+d*a(i-1,kk-1);
                    }
                    e[j-1]=a(kk-1,j-1);
                }
            }
            if(kk<=k)
				for(i=kk; i<=m; i++)	u(i-1,kk-1)=a(i-1,kk-1);
            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(FloatNotEqual(e[kk-1],0.0))
                {
					if(FloatNotEqual(e[kk],0.0))
                    {
						e[kk-1]=Abs(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)&&(FloatNotEqual(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]*a(i-1,j-1);
                    for(j=kk+1; j<=n; j++)
                      for(i=kk+1; i<=m; i++)
                          a(i-1,j-1)=a(i-1,j-1)-w[i-1]*e[j-1]/e[kk];
                }
                for(i=kk+1; i<=n; i++)
                  v(i-1,kk-1)=e[i-1];
            }
        }
    }
    mm=n;
    if(m+1<n) mm=m+1;
    if(k<n) s[k]=a(k,k);
    if(m<mm) s[mm-1]=0.0;
    if(l+1<mm) e[l]=a(l,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++)	u(i-1,j-1)=0.0;
            u(j-1,j-1)=1.0;
        }
    }
    if(k>=1)
    { 
		for(ll=1; ll<=k; ll++)
        {
			kk=k-ll+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++)
						d=d+u(i-1,kk-1)*u(i-1,j-1)/u(kk-1,kk-1);
                      d=-d;
                      for(i=kk; i<=m; i++)
                          u(i-1,j-1)=u(i-1,j-1)+d*u(i-1,kk-1);
                  }
                  for(i=kk; i<=m; i++)
					  u(i-1,kk-1)=-u(i-1,kk-1);
                  u(kk-1,kk-1)=1.0+u(kk-1,kk-1);
                  if(kk-1>=1)
                    for(i=1; i<kk; i++)	u(i-1,kk-1)=0.0;
            }
            else
            { 
				for(i=1; i<=m; i++)	u(i-1,kk-1)=0.0;
                u(kk-1,kk-1)=1.0;
            }
        }
    }
    for(ll=1; ll<=n; ll++)
    { 
		kk=n-ll+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++)
                    d=d+v(i-1,kk-1)*v(i-1,j-1)/v(kk,kk-1);
                d=-d;
                for(i=kk+1; i<=n; i++)
                    v(i-1,j-1)=v(i-1,j-1)+d*v(i-1,kk-1);
            }
        }
        for(i=1; i<=n; i++)	v(i-1,kk-1)=0.0;
        v(kk-1,kk-1)=1.0;
    }
    for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)	a(i-1,j-1)=0.0;
    m1=mm; 
	it=60;
    while(1)
    { 
		if(mm==0)
        {
			_MSV_1(a,e,s,v,m,n);
			return(1);
        }
        if(it==0)
        { 
			_MSV_1(a,e,s,v,m,n);
			return(0);
        }
        kk=mm-1;
	while((kk!=0)&&(FloatNotEqual(e[kk-1],0.0)))
    { 
		d=Abs(s[kk-1])+Abs(s[kk]);
        dd=Abs(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++)		v(i-1,kk-1)=-v(i-1,kk-1);
        }
        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++)
                {
                  d=v(i-1,kk-1);
				  v(i-1,kk-1)=v(i-1,kk); 
				  v(i-1,kk)=d;
                }
             if(kk<m)
                  for(i=1; i<=m; i++)
                  { 
                      d=u(i-1,kk-1);
					  u(i-1,kk-1)=u(i-1,kk);
					  u(i-1,kk)=d;
                  }
                kk=kk+1;
            }
            it=60;
            mm=mm-1;
        }
        else
        { 
			ks=mm;
            while((ks>kk)&&(Abs(s[ks-1])!=0.0))
            {
				d=0.0;
                if(ks!=mm) d=d+Abs(e[ks-1]);
                if(ks!=kk+1) d=d+Abs(e[ks-2]);
                dd=Abs(s[ks-1]);
                if(dd>eps*d) ks=ks-1;
                else s[ks-1]=0.0;
            }
            if(ks==kk)
            { 
				kk=kk+1;
                d=Abs(s[mm-1]);
                t=Abs(s[mm-2]);
                if(t>d) d=t;
                t=Abs(e[mm-2]);
                if(t>d) d=t;
                t=Abs(s[kk-1]);
                if(t>d) d=t;
                t=Abs(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++)
                { 
					_MSV_2(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++)
                      {
                          d=cs[0]*v(j-1,i-1)+cs[1]*v(j-1,i);
                          v(j-1,i)=-cs[1]*v(j-1,i-1)+cs[0]*v(j-1,i);
                          v(j-1,i-1)=d;
                      }
                    _MSV_2(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++)
                        { 
                            d=cs[0]*u(j-1,i-1)+cs[1]*u(j-1,i);
                            u(j-1,i)=-cs[1]*u(j-1,i-1)+cs[0]*u(j-1,i);
                            u(j-1,i-1)=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];
                        _MSV_2(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++)
                          { 
                              d=cs[0]*v(j-1,i-1)+cs[1]*v(j-1,mm-1);
                              v(j-1,mm-1)=-cs[1]*v(j-1,i-1)+cs[0]*v(j-1,mm-1);
                              v(j-1,i-1)=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];
                        _MSV_2(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++)
                          {
                              d=cs[0]*u(j-1,i-1)+cs[1]*u(j-1,kk-2);
                              u(j-1,kk-2)=-cs[1]*u(j-1,i-1)+cs[0]*u(j-1,kk-2);
                              u(j-1,i-1)=d;
                          }
                    }
                }
            }
        }
    }
    return(1);
}

//一般实矩阵的奇异值分解辅助函数
template <class _Ty>
void _MSV_1(matrix<_Ty>& a, valarray<_Ty>& e, valarray<_Ty>& s, 
											matrix<_Ty>& v, int m, int n)
{
	int i,j;
    _Ty d;
    if(m>=n) i=n;
    else i=m;
    for(j=1; j<i; j++)
    {
		a(j-1,j-1)=s[j-1];
        a(j-1,j)=e[j-1];
    }
    a(i-1,i-1)=s[i-1];
    if(m<n) a(i-1,i)=e[i-1];
    for(i=1; i<n; i++)
    for(j=i+1; j<=n; j++)
    { 
        d=v(i-1,j-1);
		v(i-1,j-1)=v(j-1,i-1); 
		v(j-1,i-1)=d;
    }
}

//一般实矩阵的奇异值分解辅助函数
template <class _Ty>
void _MSV_2(valarray<_Ty>& fg, valarray<_Ty>& cs)
{
	_Ty r,d;
	r = Abs(fg[0])+Abs(fg[1]);
    if(FloatEqual(r,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(Abs(fg[0])>Abs(fg[1]))
        {
			d=Abs(d);
            if(fg[0]<0.0) d=-d;
        }
        if(Abs(fg[1])>Abs(fg[0])||FloatEqual(Abs(fg[1]),Abs(fg[0])))
        { 
			d=Abs(d);
            if(fg[1]<0.0) d=-d;
        }
        cs[0]=fg[0]/d; 
		cs[1]=fg[1]/d;
    }
    r=1.0;
    if(Abs(fg[0])>Abs(fg[1])) r=cs[1];
    else
      if(FloatNotEqual(cs[0],0.0)) r=1.0/cs[0];
    fg[0]=d;
	fg[1]=r;
}

//广义逆的奇异值分解
template <class _Ty>
int GeneralizedInversionSingularValue(matrix<_Ty>& a, matrix<_Ty>& aa, 
							_Ty eps, matrix<_Ty>& u, matrix<_Ty>& v)
{
	int k(0), l, ka;
	int m = a.GetRowNum();
	int n = a.GetColNum();

	if(m > n) ka = m + 1;
	else ka = n + 1;
  
	int i=MatrixSingularValue(a,u,v,eps);
    if (i<0) return(0);
    int j = n;
    if(m < n) j = m;
    j = j - 1;
    while((k<=j)&&(FloatNotEqual(a(k,k),0.0))) k = k + 1;
    k = k - 1;
    for(i=0; i<n; i++)
    for(j=0; j<m; j++)
    {
		aa(i,j)=0.0;
        for(l=0; l<=k; l++)
        {
            aa(i,j)=aa(i,j)+v(l,i)*u(j,l)/a(l,l);
        }
    }
    return(1);
}

#endif			// _MATRIX_INL