matrix.cpp

很好用的有关矩阵的各种操作

					
					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;
	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];