matrix.java

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

	public Matrix transpose() 
	{
		// 构造目标矩阵
		Matrix	Trans = new Matrix(numColumns, numRows);

		// 转置各元素
		for (int i = 0 ; i < numRows ; ++i)
		{
			for (int j = 0 ; j < numColumns ; ++j)
				Trans.setElement(j, i, getElement(i, j)) ;
		}

		return Trans;
	}

	/**
	 * 实矩阵求逆的全选主元高斯－约当法
	 * 
	 * @return boolean型，求逆是否成功
	 */
	public boolean invertGaussJordan()
	{
		int i,j,k,l,u,v;
	    double d = 0, p = 0;

		// 分配内存
	    int[] pnRow = new int[numColumns];
	    int[] pnCol = new int[numColumns];

		// 消元
	    for (k=0; k<=numColumns-1; k++)
	    { 
			d=0.0;
	        for (i=k; i<=numColumns-1; i++)
			{
				for (j=k; j<=numColumns-1; j++)
				{ 
					l=i*numColumns+j; p=Math.abs(elements[l]);
					if (p>d) 
					{ 
						d=p; 
						pnRow[k]=i; 
						pnCol[k]=j;
					}
				}
			}
	        
			// 失败
			if (d == 0.0)
			{
				return false;
			}

	        if (pnRow[k] != k)
			{
				for (j=0; j<=numColumns-1; j++)
				{ 
					u=k*numColumns+j; 
					v=pnRow[k]*numColumns+j;
					p=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=p;
				}
			}
	        
			if (pnCol[k] != k)
			{
				for (i=0; i<=numColumns-1; i++)
	            { 
					u=i*numColumns+k; 
					v=i*numColumns+pnCol[k];
					p=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=p;
	            }
			}

	        l=k*numColumns+k;
	        elements[l]=1.0/elements[l];
	        for (j=0; j<=numColumns-1; j++)
			{
				if (j != k)
	            { 
					u=k*numColumns+j; 
					elements[u]=elements[u]*elements[l];
				}
			}

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

	        for (i=0; i<=numColumns-1; i++)
			{
				if (i!=k)
	            { 
					u=i*numColumns+k; 
					elements[u]=-elements[u]*elements[l];
				}
			}
	    }

	    // 调整恢复行列次序
	    for (k=numColumns-1; k>=0; k--)
	    { 
			if (pnCol[k]!=k)
			{
				for (j=0; j<=numColumns-1; j++)
	            { 
					u=k*numColumns+j; 
					v=pnCol[k]*numColumns+j;
					p=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=p;
	            }
			}

	        if (pnRow[k]!=k)
			{
				for (i=0; i<=numColumns-1; i++)
	            { 
					u=i*numColumns+k; 
					v=i*numColumns+pnRow[k];
					p=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=p;
	            }
			}
	    }

		// 成功返回
		return true;
	}

	/**
	 * 复矩阵求逆的全选主元高斯－约当法
	 * 
	 * @param mtxImag - 复矩阵的虚部矩阵，当前矩阵为复矩阵的实部
	 * @return boolean型，求逆是否成功
	 */
	public boolean invertGaussJordan(Matrix mtxImag)
	{
		int i,j,k,l,u,v,w;
	    double p,q,s,t,d,b;

		// 分配内存
	    int[] pnRow = new int[numColumns];
	    int[] pnCol = new int[numColumns];

		// 消元
	    for (k=0; k<=numColumns-1; k++)
	    { 
			d=0.0;
	        for (i=k; i<=numColumns-1; i++)
			{
				for (j=k; j<=numColumns-1; j++)
				{ 
					u=i*numColumns+j;
					p=elements[u]*elements[u]+mtxImag.elements[u]*mtxImag.elements[u];
					if (p>d) 
					{ 
						d=p; 
						pnRow[k]=i; 
						pnCol[k]=j;
					}
				}
			}

			// 失败
	        if (d == 0.0)
	        { 
	            return false;
	        }

	        if (pnRow[k]!=k)
			{
				for (j=0; j<=numColumns-1; j++)
	            { 
					u=k*numColumns+j; 
					v=pnRow[k]*numColumns+j;
					t=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=t;
					t=mtxImag.elements[u]; 
					mtxImag.elements[u]=mtxImag.elements[v]; 
					mtxImag.elements[v]=t;
	            }
			}

	        if (pnCol[k]!=k)
			{
				for (i=0; i<=numColumns-1; i++)
	            { 
					u=i*numColumns+k; 
					v=i*numColumns+pnCol[k];
					t=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=t;
					t=mtxImag.elements[u]; 
					mtxImag.elements[u]=mtxImag.elements[v]; 
					mtxImag.elements[v]=t;
	            }
			}

	        l=k*numColumns+k;
	        elements[l]=elements[l]/d; mtxImag.elements[l]=-mtxImag.elements[l]/d;
	        for (j=0; j<=numColumns-1; j++)
			{
				if (j!=k)
	            { 
					u=k*numColumns+j;
					p=elements[u]*elements[l]; 
					q=mtxImag.elements[u]*mtxImag.elements[l];
					s=(elements[u]+mtxImag.elements[u])*(elements[l]+mtxImag.elements[l]);
					elements[u]=p-q; 
					mtxImag.elements[u]=s-p-q;
	            }
			}

	        for (i=0; i<=numColumns-1; i++)
			{
				if (i!=k)
	            { 
					v=i*numColumns+k;
					for (j=0; j<=numColumns-1; j++)
					{
						if (j!=k)
						{ 
							u=k*numColumns+j;  
							w=i*numColumns+j;
							p=elements[u]*elements[v]; 
							q=mtxImag.elements[u]*mtxImag.elements[v];
							s=(elements[u]+mtxImag.elements[u])*(elements[v]+mtxImag.elements[v]);
							t=p-q; 
							b=s-p-q;
							elements[w]=elements[w]-t;
							mtxImag.elements[w]=mtxImag.elements[w]-b;
						}
					}
	            }
			}

	        for (i=0; i<=numColumns-1; i++)
			{
				if (i!=k)
	            { 
					u=i*numColumns+k;
					p=elements[u]*elements[l]; 
					q=mtxImag.elements[u]*mtxImag.elements[l];
					s=(elements[u]+mtxImag.elements[u])*(elements[l]+mtxImag.elements[l]);
					elements[u]=q-p; 
					mtxImag.elements[u]=p+q-s;
	            }
			}
	    }

	    // 调整恢复行列次序
	    for (k=numColumns-1; k>=0; k--)
	    { 
			if (pnCol[k]!=k)
			{
				for (j=0; j<=numColumns-1; j++)
	            { 
					u=k*numColumns+j; 
					v=pnCol[k]*numColumns+j;
					t=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=t;
					t=mtxImag.elements[u]; 
					mtxImag.elements[u]=mtxImag.elements[v]; 
					mtxImag.elements[v]=t;
	            }
			}

	        if (pnRow[k]!=k)
			{
				for (i=0; i<=numColumns-1; i++)
	            { 
					u=i*numColumns+k; 
					v=i*numColumns+pnRow[k];
					t=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=t;
					t=mtxImag.elements[u]; 
					mtxImag.elements[u]=mtxImag.elements[v]; 
					mtxImag.elements[v]=t;
	            }
			}
	    }

		// 成功返回
		return true;
	}

	/**
	 * 对称正定矩阵的求逆
	 * 
	 * @return boolean型，求逆是否成功
	 */
	public boolean invertSsgj()
	{ 
		int i, j ,k, m;
	    double w, g;

		// 临时内存
	    double[] pTmp = new double[numColumns];

		// 逐列处理
	    for (k=0; k<=numColumns-1; k++)
	    { 
			w=elements[0];
	        if (w == 0.0)
	        { 
				return false;
			}

	        m=numColumns-k-1;
	        for (i=1; i<=numColumns-1; i++)
	        { 
				g=elements[i*numColumns]; 
				pTmp[i]=g/w;
	            if (i<=m) 
					pTmp[i]=-pTmp[i];
	            for (j=1; j<=i; j++)
	              elements[(i-1)*numColumns+j-1]=elements[i*numColumns+j]+g*pTmp[j];
	        }

	        elements[numColumns*numColumns-1]=1.0/w;
	        for (i=1; i<=numColumns-1; i++)
				elements[(numColumns-1)*numColumns+i-1]=pTmp[i];
	    }

		// 行列调整
	    for (i=0; i<=numColumns-2; i++)
			for (j=i+1; j<=numColumns-1; j++)
				elements[i*numColumns+j]=elements[j*numColumns+i];

		return true;
	}

	/**
	 * 托伯利兹矩阵求逆的埃兰特方法
	 * 
	 * @return boolean型，求逆是否成功
	 */
	public boolean invertTrench()
	{ 
		int i,j,k;
	    double a,s;

		// 上三角元素
		double[] t = new double[numColumns];
		// 下三角元素
		double[] tt = new double[numColumns];

		// 上、下三角元素赋值
		for (i=0; i<numColumns; ++i)
		{
			t[i] = getElement(0, i);
		    tt[i] = getElement(i, 0);
		}

		// 临时缓冲区
		double[] c = new double[numColumns];
		double[] r = new double[numColumns];
		double[] p = new double[numColumns];

		// 非Toeplitz矩阵，返回
	    if (t[0] == 0.0)
	    { 
			return false;
	    }

	    a=t[0]; 
		c[0]=tt[1]/t[0]; 
		r[0]=t[1]/t[0];

	    for (k=0; k<=numColumns-3; k++)
	    { 
			s=0.0;
	        for (j=1; j<=k+1; j++)
				s=s+c[k+1-j]*tt[j];

	        s=(s-tt[k+2])/a;
			for (i=0; i<=k; i++)
				p[i]=c[i]+s*r[k-i];

	        c[k+1]=-s;
	        s=0.0;
	        for (j=1; j<=k+1; j++)
	          s=s+r[k+1-j]*t[j];
	        
			s=(s-t[k+2])/a;
	        for (i=0; i<=k; i++)
	        { 
				r[i]=r[i]+s*c[k-i];
	            c[k-i]=p[k-i];
	        }

	        r[k+1]=-s;
			a=0.0;
	        for (j=1; j<=k+2; j++)
	          a=a+t[j]*c[j-1];

	        a=t[0]-a;

			// 求解失败
	        if (a == 0.0)
			{ 
				return false;
			}
	    }

	    elements[0]=1.0/a;
	    for (i=0; i<=numColumns-2; i++)
	    { 
			k=i+1; 
			j=(i+1)*numColumns;
	        elements[k]=-r[i]/a; 
			elements[j]=-c[i]/a;
	    }

	   for (i=0; i<=numColumns-2; i++)
		{
			for (j=0; j<=numColumns-2; j++)
			{ 
				k=(i+1)*numColumns+j+1;
				elements[k]=elements[i*numColumns+j]-c[i]*elements[j+1];
				elements[k]=elements[k]+c[numColumns-j-2]*elements[numColumns-i-1];
			}
		}

		return true;
	}

	/**
	 * 求行列式值的全选主元高斯消去法
	 * 
	 * @return double型，行列式的值
	 */
	public double computeDetGauss()
	{ 
		int i,j,k,is = 0,js = 0,l,u,v;
	    double f,det,q,d;
	    
		// 初值
		f=1.0; 
		det=1.0;
	    
		// 消元
		for (k=0; k<=numColumns-2; k++)
	    { 
			q=0.0;
	        for (i=k; i<=numColumns-1; i++)
			{
				for (j=k; j<=numColumns-1; j++)
				{ 
					l=i*numColumns+j; 
					d=Math.abs(elements[l]);
					if (d>q) 
					{ 
						q=d; 
						is=i; 
						js=j;
					}
				}
			}

	        if (q == 0.0)
	        { 
				det=0.0; 
				return(det);
			}
	        
			if (is!=k)
	        { 
				f=-f;
	            for (j=k; j<=numColumns-1; j++)
	            { 
					u=k*numColumns+j; 
					v=is*numColumns+j;
	                d=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=d;
	            }
	        }
	        
			if (js!=k)
	        { 
				f=-f;
	            for (i=k; i<=numColumns-1; i++)
	            {
					u=i*numColumns+js; 
					v=i*numColumns+k;
	                d=elements[u]; 
					elements[u]=elements[v]; 
					elements[v]=d;
	            }
	        }

	        l=k*numColumns+k;
	        det=det*elements[l];
	        for (i=k+1; i<=numColumns-1; i++)
	        { 
				d=elements[i*numColumns+k]/elements[l];
	            for (j=k+1; j<=numColumns-1; j++)
	            { 
					u=i*numColumns+j;
	                elements[u]=elements[u]-d*elements[k*numColumns+j];
	            }
	        }
	    }
	    
		// 求值
		det=f*det*elements[numColumns*numColumns-1];

	    return(det);
	}

	/**
	 * 求矩阵秩的全选主元高斯消去法
	 * 
	 * @return int型，矩阵的秩
	 */
	public int computeRankGauss()
	{ 
		int i,j,k,nn,is = 0,js = 0,l,ll,u,v;
	    double q,d;
	    
		// 秩小于等于行列数
		nn = numRows;
	    if (numRows >= numColumns) 
			nn = numColumns;

	    k=0;

		// 消元求解
	    for (l=0; l<=nn-1; l++)
	    { 
			q=0.0;
	        for (i=l; i<=numRows-1; i++)
			{
				for (j=l; j<=numColumns-1; j++)
				{ 
					ll=i*numColumns+j; 
					d=Math.abs(elements[ll]);
					if (d>q) 
					{ 
						q=d; 
						is=i;