matrix.java

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

						js=j;
					}
				}
			}

	        if (q == 0.0) 
				return(k);

	        k=k+1;
	        if (is!=l)
	        { 
				for (j=l; j<=numColumns-1; j++)
	            { 
					u=l*numColumns+j; 
					v=is*numColumns+j;
	                d=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=d;
	            }
	        }
	        if (js!=l)
	        { 
				for (i=l; i<=numRows-1; i++)
	            { 
					u=i*numColumns+js; 
					v=i*numColumns+l;
	                d=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=d;
	            }
	        }
	        
			ll=l*numColumns+l;
	        for (i=l+1; i<=numColumns-1; i++)
	        { 
				d=elements[i*numColumns+l]/elements[ll];
	            for (j=l+1; j<=numColumns-1; j++)
	            { 
					u=i*numColumns+j;
	                elements[u]=elements[u]-d*elements[l*numColumns+j];
	            }
	        }
	    }
	    
		return(k);
	}

	/**
	 * 对称正定矩阵的乔里斯基分解与行列式的求值
	 * 
	 * @param realDetValue - 返回行列式的值
	 * @return boolean型，求解是否成功
	 */
	public boolean computeDetCholesky(Real realDetValue)
	{ 
		int i,j,k,u,l;
	    double d, dblDet;
	    
		// 不满足求解要求
		if (elements[0] <= 0.0)
			return false;

		// 乔里斯基分解

	    elements[0]=Math.sqrt(elements[0]);
	    d=elements[0];

	    for (i=1; i<=numColumns-1; i++)
	    { 
			u=i*numColumns; 
			elements[u]=elements[u]/elements[0];
		}
	    
		for (j=1; j<=numColumns-1; j++)
	    { 
			l=j*numColumns+j;
	        for (k=0; k<=j-1; k++)
	        { 
				u=j*numColumns+k; 
				elements[l]=elements[l]-elements[u]*elements[u];
			}
	        
			if (elements[l] <= 0.0)
				return false;

	        elements[l]=Math.sqrt(elements[l]);
	        d=d*elements[l];
	        
			for (i=j+1; i<=numColumns-1; i++)
	        { 
				u=i*numColumns+j;
	            for (k=0; k<=j-1; k++)
					elements[u]=elements[u]-elements[i*numColumns+k]*elements[j*numColumns+k];
	            
				elements[u]=elements[u]/elements[l];
	        }
	    }
	    
		// 行列式求值
		dblDet=d*d;
		realDetValue.setValue(dblDet);
		
		// 下三角矩阵
	    for (i=0; i<=numColumns-2; i++)
			for (j=i+1; j<=numColumns-1; j++)
				elements[i*numColumns+j]=0.0;

		return true;
	}

	/**
	 * 矩阵的三角分解，分解成功后，原矩阵将成为Q矩阵
	 * 
	 * @param mtxL - 返回分解后的L矩阵
	 * @param mtxU - 返回分解后的U矩阵
	 * @return boolean型，求解是否成功
	 */
	public boolean splitLU(Matrix mtxL, Matrix mtxU)
	{ 
		int i,j,k,w,v,ll;
	    
		// 初始化结果矩阵
		if (! mtxL.init(numColumns, numColumns) ||
			! mtxU.init(numColumns, numColumns))
			return false;

		for (k=0; k<=numColumns-2; k++)
	    { 
			ll=k*numColumns+k;
			if (elements[ll] == 0.0)
				return false;

	        for (i=k+1; i<=numColumns-1; i++)
			{ 
				w=i*numColumns+k; 
				elements[w]=elements[w]/elements[ll];
			}

	        for (i=k+1; i<=numColumns-1; i++)
	        { 
				w=i*numColumns+k;
	            for (j=k+1; j<=numColumns-1; j++)
	            { 
					v=i*numColumns+j;
	                elements[v]=elements[v]-elements[w]*elements[k*numColumns+j];
	            }
	        }
	    }
	    
		for (i=0; i<=numColumns-1; i++)
	    {
			for (j=0; j<i; j++)
	        { 
				w=i*numColumns+j; 
				mtxL.elements[w]=elements[w]; 
				mtxU.elements[w]=0.0;
			}

	        w=i*numColumns+i;
	        mtxL.elements[w]=1.0; 
			mtxU.elements[w]=elements[w];
	        
			for (j=i+1; j<=numColumns-1; j++)
	        { 
				w=i*numColumns+j; 
				mtxL.elements[w]=0.0; 
				mtxU.elements[w]=elements[w];
			}
	    }

		return true;
	}

	/**
	 * 一般实矩阵的QR分解，分解成功后，原矩阵将成为R矩阵
	 * 
	 * @param mtxQ - 返回分解后的Q矩阵
	 * @return boolean型，求解是否成功
	 */
	public boolean splitQR(Matrix mtxQ)
	{ 
		int i,j,k,l,nn,p,jj;
	    double u,alpha,w,t;
	    
		if (numRows < numColumns)
			return false;

		// 初始化Q矩阵
		if (! mtxQ.init(numRows, numRows))
			return false;

		// 对角线元素单位化
	    for (i=0; i<=numRows-1; i++)
		{
			for (j=0; j<=numRows-1; j++)
			{ 
				l=i*numRows+j; 
				mtxQ.elements[l]=0.0;
				if (i==j) 
					mtxQ.elements[l]=1.0;
			}
		}

		// 开始分解

	    nn=numColumns;
	    if (numRows == numColumns) 
			nn=numRows-1;

	    for (k=0; k<=nn-1; k++)
	    { 
			u=0.0; 
			l=k*numColumns+k;
	        for (i=k; i<=numRows-1; i++)
	        { 
				w=Math.abs(elements[i*numColumns+k]);
	            if (w>u) 
					u=w;
	        }
	        
			alpha=0.0;
	        for (i=k; i<=numRows-1; i++)
	        { 
				t=elements[i*numColumns+k]/u; 
				alpha=alpha+t*t;
			}

	        if (elements[l]>0.0) 
				u=-u;

	        alpha=u*Math.sqrt(alpha);
	        if (alpha == 0.0)
				return false;

	        u=Math.sqrt(2.0*alpha*(alpha-elements[l]));
	        if ((u+1.0)!=1.0)
	        { 
				elements[l]=(elements[l]-alpha)/u;
	            for (i=k+1; i<=numRows-1; i++)
	            { 
					p=i*numColumns+k; 
					elements[p]=elements[p]/u;
				}
	            
				for (j=0; j<=numRows-1; j++)
	            { 
					t=0.0;
	                for (jj=k; jj<=numRows-1; jj++)
						t=t+elements[jj*numColumns+k]*mtxQ.elements[jj*numRows+j];

	                for (i=k; i<=numRows-1; i++)
	                { 
						p=i*numRows+j; 
						mtxQ.elements[p]=mtxQ.elements[p]-2.0*t*elements[i*numColumns+k];
					}
	            }
	            
				for (j=k+1; j<=numColumns-1; j++)
	            { 
					t=0.0;
	                
					for (jj=k; jj<=numRows-1; jj++)
						t=t+elements[jj*numColumns+k]*elements[jj*numColumns+j];
	                
					for (i=k; i<=numRows-1; i++)
	                { 
						p=i*numColumns+j; 
						elements[p]=elements[p]-2.0*t*elements[i*numColumns+k];
					}
	            }
	            
				elements[l]=alpha;
	            for (i=k+1; i<=numRows-1; i++)
					elements[i*numColumns+k]=0.0;
	        }
	    }
	    
		// 调整元素
		for (i=0; i<=numRows-2; i++)
		{
			for (j=i+1; j<=numRows-1;j++)
			{ 
				p=i*numRows+j; 
				l=j*numRows+i;
				t=mtxQ.elements[p]; 
				mtxQ.elements[p]=mtxQ.elements[l]; 
				mtxQ.elements[l]=t;
			}
		}

		return true;
	}

	/**
	 * 一般实矩阵的奇异值分解，分解成功后，原矩阵对角线元素就是矩阵的奇异值
	 * 
	 * @param mtxU - 返回分解后的U矩阵
	 * @param mtxV - 返回分解后的V矩阵
	 * @param eps - 计算精度
	 * @return boolean型，求解是否成功
	 */
	public boolean splitUV(Matrix mtxU, Matrix mtxV, double eps)
	{ 
		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;
	    double[] fg = new double[2];
	    double[] cs = new double[2];

		int m = numRows;
		int n = numColumns;

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

		// 临时缓冲区
		int ka = Math.max(m, n) + 1;
		double[] s = new double[ka];
		double[] e = new double[ka];
		double[] 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+elements[ix]*elements[ix];
					}

	                s[kk-1]=Math.sqrt(d);
	                if (s[kk-1]!=0.0)
	                { 
						ix=(kk-1)*n+kk-1;
	                    if (elements[ix]!=0.0)
	                    { 
							s[kk-1]=Math.abs(s[kk-1]);
	                        if (elements[ix]<0.0) 
								s[kk-1]=-s[kk-1];
	                    }
	                    
						for (i=kk; i<=m; i++)
	                    { 
							iy=(i-1)*n+kk-1;
	                        elements[iy]=elements[iy]/s[kk-1];
	                    }
	                    
						elements[ix]=1.0+elements[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+elements[ix]*elements[iy];
	                        }
	                        
							d=-d/elements[(kk-1)*n+kk-1];
	                        for (i=kk; i<=m; i++)
	                        { 
								ix=(i-1)*n+j-1;
	                            iy=(i-1)*n+kk-1;
	                            elements[ix]=elements[ix]+d*elements[iy];
	                        }
	                    }
	                    
						e[j-1]=elements[(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.elements[ix]=elements[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]=Math.sqrt(d);
	                if (e[kk-1]!=0.0)
	                { 
						if (e[kk]!=0.0)
	                    { 
							e[kk-1]=Math.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)&& (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]*elements[(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;
								elements[ix]=elements[ix]-w[i-1]*e[j-1]/e[kk];
	                        }
						}
	                }
	                
					for (i=kk+1; i<=n; i++)
	                  mtxV.elements[(i-1)*n+kk-1]=e[i-1];
	            }
	        }
	    }
	    
		mm=n;
	    if (m+1<n) 
			mm=m+1;
	    if (k<n) 
			s[k]=elements[k*n+k];
	    if (m<mm) 
			s[mm-1]=0.0;
	    if (l+1<mm) 
			e[l]=elements[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.elements[(i-1)*m+j-1]=0.0;
	            mtxU.elements[(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.elements[ix]*mtxU.elements[iy]/mtxU.elements[iz];
							}

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

					mtxU.elements[iz]=1.0+mtxU.elements[iz];
					if (kk-1>=1)
					{
						for (i=1; i<=kk-1; i++)
							mtxU.elements[(i-1)*m+kk-1]=0.0;
					}
				}
	            else
	            { 
					for (i=1; i<=m; i++)
						mtxU.elements[(i-1)*m+kk-1]=0.0;
	                mtxU.elements[(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.elements[ix]*mtxV.elements[iy]/mtxV.elements[iz];
	                }
	                
					d=-d;
	                for (i=kk+1; i<=n; i++)
	                { 
						ix=(i-1)*n+j-1; 
						iy=(i-1)*n+kk-1;
	                    mtxV.elements[ix]=mtxV.elements[ix]+d*mtxV.elements[iy];
	                }
	            }
	        }