matrix.java

实现用于工程计算的的数学方法

	        
			for (i=1; i<=n; i++)
				mtxV.elements[(i-1)*n+kk-1]=0.0;
	        
			mtxV.elements[iz-n]=1.0;
	    }
	    
		for (i=1; i<=m; i++)
			for (j=1; j<=n; j++)
				elements[(i-1)*n+j-1]=0.0;
	    
		m1=mm; 
		it=60;
	    while (true)
	    { 
			if (mm==0)
	        { 
				ppp(elements,e,s,mtxV.elements,m,n);
	            return true;
	        }
	        if (it==0)
	        { 
				ppp(elements,e,s,mtxV.elements,m,n);
	            return false;
	        }
	        
			kk=mm-1;
			while ((kk!=0) && (Math.abs(e[kk-1])!=0.0))
	        { 
				d=Math.abs(s[kk-1])+Math.abs(s[kk]);
	            dd=Math.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++)
	                { 
						ix=(i-1)*n+kk-1; 
						mtxV.elements[ix]=-mtxV.elements[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.elements[ix]; 
								mtxV.elements[ix]=mtxV.elements[iy]; 
								mtxV.elements[iy]=d;
							}
						}

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

						kk=kk+1;
	            }
	            
				it=60;
	            mm=mm-1;
	        }
	        else
	        { 
				ks=mm;
	            while ((ks>kk) && (Math.abs(s[ks-1])!=0.0))
	            { 
					d=0.0;
	                if (ks!=mm) 
						d=d+Math.abs(e[ks-1]);
	                if (ks!=kk+1) 
						d=d+Math.abs(e[ks-2]);
	                
					dd=Math.abs(s[ks-1]);
	                if (dd>eps*d) 
						ks=ks-1;
	                else 
						s[ks-1]=0.0;
	            }
	            
				if (ks==kk)
	            { 
					kk=kk+1;
	                d=Math.abs(s[mm-1]);
	                t=Math.abs(s[mm-2]);
	                if (t>d) 
						d=t;
	                
					t=Math.abs(e[mm-2]);
	                if (t>d) 
						d=t;
	                
					t=Math.abs(s[kk-1]);
	                if (t>d) 
						d=t;
	                
					t=Math.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=Math.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.elements[ix]+cs[1]*mtxV.elements[iy];
								mtxV.elements[iy]=-cs[1]*mtxV.elements[ix]+cs[0]*mtxV.elements[iy];
								mtxV.elements[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.elements[ix]+cs[1]*mtxU.elements[iy];
									mtxU.elements[iy]=-cs[1]*mtxU.elements[ix]+cs[0]*mtxU.elements[iy];
									mtxU.elements[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.elements[ix]+cs[1]*mtxV.elements[iy];
									mtxV.elements[iy]=-cs[1]*mtxV.elements[ix]+cs[0]*mtxV.elements[iy];
									mtxV.elements[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.elements[ix]+cs[1]*mtxU.elements[iy];
									mtxU.elements[iy]=-cs[1]*mtxU.elements[ix]+cs[0]*mtxU.elements[iy];
									mtxU.elements[ix]=d;
	                            }
							}
	                    }
	                }
	            }
	        }
	    }
	}

	/**
	 * 内部函数，由SplitUV函数调用
	 */
	private void 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函数调用
	 */
	private void sss(double[] fg, double[] cs)
	{ 
		double r,d;
	    
		if ((Math.abs(fg[0])+Math.abs(fg[1]))==0.0)
	    { 
			cs[0]=1.0; 
			cs[1]=0.0; 
			d=0.0;
		}
	    else 
	    { 
			d=Math.sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
	        if (Math.abs(fg[0])>Math.abs(fg[1]))
	        { 
				d=Math.abs(d);
	            if (fg[0]<0.0) 
					d=-d;
	        }
	        if (Math.abs(fg[1])>=Math.abs(fg[0]))
	        { 
				d=Math.abs(d);
	            if (fg[1]<0.0) 
					d=-d;
	        }
	        
			cs[0]=fg[0]/d; 
			cs[1]=fg[1]/d;
	    }
	    
		r=1.0;
	    if (Math.abs(fg[0])>Math.abs(fg[1])) 
			r=cs[1];
	    else if (cs[0]!=0.0) 
			r=1.0/cs[0];

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

	/**
	 * 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值
	 * 
	 * @param mtxAP - 返回原矩阵的广义逆矩阵
	 * @param mtxU - 返回分解后的U矩阵
	 * @param mtxV - 返回分解后的V矩阵
	 * @param eps - 计算精度
	 * @return boolean型，求解是否成功
	 */
	public boolean invertUV(Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
	{ 
		int i,j,k,l,t,p,q,f;

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

		int m = numRows;
		int n = numColumns;

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

		// 计算广义逆矩阵

	    j=n;
	    if (m<n) 
			j=m;
	    j=j-1;
	    k=0;
	    while ((k<=j) && (elements[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.elements[t]=0.0;
				for (l=0; l<=k; l++)
				{ 
					f=l*n+i; 
					p=j*m+l; 
					q=l*n+l;
					mtxAP.elements[t]=mtxAP.elements[t]+mtxV.elements[f]*mtxU.elements[p]/elements[q];
				}
			}
		}

	    return true;
	}

	/**
	 * 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
	 * 
	 * @param mtxQ - 返回豪斯荷尔德变换的乘积矩阵Q
	 * @param mtxT - 返回求得的对称三对角阵
	 * @param dblB - 一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元素
	 * @param dblC - 一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的
	 *               次对角线元素
	 * @return boolean型，求解是否成功
	 */
	public boolean makeSymTri(Matrix mtxQ, Matrix mtxT, double[] dblB, double[] dblC)
	{ 
		int i,j,k,u;
	    double h,f,g,h2;
	    
		// 初始化矩阵Q和T
		if (! mtxQ.init(numColumns, numColumns) ||
			! mtxT.init(numColumns, numColumns))
			return false;

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

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

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

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

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

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

	                dblC[j]=g/h;
	                f=f+g*mtxQ.elements[j*numColumns+i];
	            }
	            
				h2=f/(h+h);
	            for (j=0; j<=i-1; j++)
	            { 
					f=mtxQ.elements[i*numColumns+j];
	                g=dblC[j]-h2*f;
	                dblC[j]=g;
	                for (k=0; k<=j; k++)
	                { 
						u=j*numColumns+k;
	                    mtxQ.elements[u]=mtxQ.elements[u]-f*dblC[k]-g*mtxQ.elements[i*numColumns+k];
	                }
	            }
	            
				dblB[i]=h;
	        }
	    }
	    
		for (i=0; i<=numColumns-2; i++) 
			dblC[i]=dblC[i+1];
	    
		dblC[numColumns-1]=0.0;
	    dblB[0]=0.0;
	    for (i=0; i<=numColumns-1; i++)
	    { 
			if ((dblB[i]!=(double)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.elements[i*numColumns+k]*mtxQ.elements[k*numColumns+j];

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

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

	    // 构造对称三对角矩阵
	    for (i=0; i<numColumns; ++i)
		{
		    for (j=0; j<numColumns; ++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;
	}

	/**
	 * 实对称三对角阵的全部特征值与特征向量的计算
	 * 
	 * @param dblB - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
	 *			     返回时存放全部特征值。
	 * @param dblC - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的
	 *               次对角线元素
	 * @param mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
	 *			     如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积
	 *               矩阵Q，则返回矩阵A的特征值向量矩阵。其中第i列为与数组dblB
	 *               中第j个特征值对应的特征向量。
	 * @param nMaxIt - 迭代次数
	 * @param eps - 计算精度
	 * @return boolean型，求解是否成功
	 */
	public boolean computeEvSymTri(double[] dblB, double[] dblC, Matrix mtxQ, int nMaxIt, double eps)
	{
		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;
	    
		// 迭代计算