matrix.cs

工程上可以直接使用的关于数值计算方法的原代码

					det=0.0; 
					return(det);
				}
	        
				if (nis!=k)
				{ 
					f=-f;
					for (j=k; j<=numColumns-1; j++)
					{ 
						u=k*numColumns+j; 
						v=nis*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,nis = 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; 
							nis=i; 
							js=j;
						}
					}
				}

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

				k=k+1;
				if (nis!=l)
				{ 
					for (j=l; j<=numColumns-1; j++)
					{ 
						u=l*numColumns+j; 
						v=nis*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 bool型，求解是否成功
		 */
		public bool ComputeDetCholesky(ref double realDetValue)
		{ 
			int i,j,k,u,l;
			double d;
	    
			// 不满足求解要求
			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];
				}
			}
	    
			// 行列式求值
			realDetValue=d*d;
		
			// 下三角矩阵
			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 bool型，求解是否成功
		 */
		public bool 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 bool型，求解是否成功
		 */
		public bool 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 bool型，求解是否成功
		 */
		public bool 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++)