* @param eps - 計算精度
		 * @return bool型，求解是否成功
		 */
		public bool 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;
	    
			// 迭代計算

			for (j=0; j<=n-1; j++)
			{ 
				it=0;
				h=eps*(Math.Abs(dblB[j])+Math.Abs(dblC[j]));
				if (h>d) 
					d=h;
	        
				m=j;
				while ((m<=n-1) && (Math.Abs(dblC[m])>d)) 
					m=m+1;
	        
				if (m!=j)
				{ 
					do
					{ 
						if (it==nMaxIt)
							return false;

						it=it+1;
						g=dblB[j];
						p=(dblB[j+1]-g)/(2.0*dblC[j]);
						r=Math.Sqrt(p*p+1.0);
						if (p>=0.0) 
							dblB[j]=dblC[j]/(p+r);
						else 
							dblB[j]=dblC[j]/(p-r);
	                
						h=g-dblB[j];
						for (i=j+1; i<=n-1; i++)
							dblB[i]=dblB[i]-h;
	                
						f=f+h; 
						p=dblB[m]; 
						e=1.0; 
						s=0.0;
						for (i=m-1; i>=j; i--)
						{ 
							g=e*dblC[i]; 
							h=e*p;
							if (Math.Abs(p)>=Math.Abs(dblC[i]))
							{ 
								e=dblC[i]/p; 
								r=Math.Sqrt(e*e+1.0);
								dblC[i+1]=s*p*r; 
								s=e/r; 
								e=1.0/r;
							}
							else
							{ 
								e=p/dblC[i]; 
								r=Math.Sqrt(e*e+1.0);
								dblC[i+1]=s*dblC[i]*r;
								s=1.0/r; 
								e=e/r;
							}
	                    
							p=e*dblB[i]-s*g;
							dblB[i+1]=h+s*(e*g+s*dblB[i]);
							for (k=0; k<=n-1; k++)
							{ 
								u=k*n+i+1; 
								v=u-1;
								h=mtxQ.elements[u]; 
								mtxQ.elements[u]=s*mtxQ.elements[v]+e*h;
								mtxQ.elements[v]=e*mtxQ.elements[v]-s*h;
							}
						}
	                
						dblC[j]=s*p; 
						dblB[j]=e*p;
	            
					} while (Math.Abs(dblC[j])>d);
				}
	        
				dblB[j]=dblB[j]+f;
			}
	    
			for (i=0; i<=n-1; i++)
			{ 
				k=i; 
				p=dblB[i];
				if (i+1<=n-1)
				{ 
					j=i+1;
					while ((j<=n-1) && (dblB[j]<=p))
					{ 
						k=j; 
						p=dblB[j]; 
						j=j+1;
					}
				}

				if (k!=i)
				{ 
					dblB[k]=dblB[i]; 
					dblB[i]=p;
					for (j=0; j<=n-1; j++)
					{ 
						u=j*n+i; 
						v=j*n+k;
						p=mtxQ.elements[u]; 
						mtxQ.elements[u]=mtxQ.elements[v]; 
						mtxQ.elements[v]=p;
					}
				}
			}
	    
			return true;
		}

		/**
		 * 約化一般實矩陣為赫申伯格矩陣的初等相似變換法
		 */
		public void MakeHberg()
		{ 
			int i = 0,j,k,u,v;
			double d,t;

			for (k=1; k<=numColumns-2; k++)
			{ 
				d=0.0;
				for (j=k; j<=numColumns-1; j++)
				{ 
					u=j*numColumns+k-1; 
					t=elements[u];
					if (Math.Abs(t)>Math.Abs(d))
					{ 
						d=t; 
						i=j;
					}
				}
	        
				if (d != 0.0)
				{ 
					if (i!=k)
					{ 
						for (j=k-1; j<=numColumns-1; j++)
						{ 
							u=i*numColumns+j; 
							v=k*numColumns+j;
							t=elements[u]; 
							elements[u]=elements[v]; 
							elements[v]=t;
						}
	                
						for (j=0; j<=numColumns-1; j++)
						{ 
							u=j*numColumns+i; 
							v=j*numColumns+k;
							t=elements[u]; 
							elements[u]=elements[v]; 
							elements[v]=t;
						}
					}
	            
					for (i=k+1; i<=numColumns-1; i++)
					{ 
						u=i*numColumns+k-1; 
						t=elements[u]/d; 
						elements[u]=0.0;
						for (j=k; j<=numColumns-1; j++)
						{ 
							v=i*numColumns+j;
							elements[v]=elements[v]-t*elements[k*numColumns+j];
						}
	                
						for (j=0; j<=numColumns-1; j++)
						{ 
							v=j*numColumns+k;
							elements[v]=elements[v]+t*elements[j*numColumns+i];
						}
					}
				}
			}
		}

		/**
		 * 求赫申伯格矩陣全部特征值的QR方法
		 * 
		 * @param dblU - 一維數組，長度為矩陣的階數，返回時存放特征值的實部
		 * @param dblV - 一維數組，長度為矩陣的階數，返回時存放特征值的虛部
		 * @param nMaxIt - 迭代次數
		 * @param eps - 計算精度
		 * @return bool型，求解是否成功
		 */
		public bool ComputeEvHBerg(double[] dblU, double[] dblV, int nMaxIt, double eps)
		{ 
			int m,it,i,j,k,l,ii,jj,kk,ll;
			double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
	    
			int n = numColumns;

			it=0; 
			m=n;
			while (m!=0)
			{ 
				l=m-1;
				while ((l>0) && (Math.Abs(elements[l*n+l-1]) > 
					eps*(Math.Abs(elements[(l-1)*n+l-1])+Math.Abs(elements[l*n+l])))) 
					l=l-1;

				ii=(m-1)*n+m-1; 
				jj=(m-1)*n+m-2;
				kk=(m-2)*n+m-1; 
				ll=(m-2)*n+m-2;
				if (l==m-1)
				{ 
					dblU[m-1]=elements[(m-1)*n+m-1]; 
					dblV[m-1]=0.0;
					m=m-1; 
					it=0;
				}
				else if (l==m-2)
				{ 
					b=-(elements[ii]+elements[ll]);
					c=elements[ii]*elements[ll]-elements[jj]*elements[kk];
					w=b*b-4.0*c;
					y=Math.Sqrt(Math.Abs(w));
					if (w>0.0)
					{ 
						xy=1.0;
						if (b<0.0) 
							xy=-1.0;
						dblU[m-1]=(-b-xy*y)/2.0;
						dblU[m-2]=c/dblU[m-1];
						dblV[m-1]=0.0; dblV[m-2]=0.0;
					}
					else
					{ 
						dblU[m-1]=-b/2.0; 
						dblU[m-2]=dblU[m-1];
						dblV[m-1]=y/2.0; 
						dblV[m-2]=-dblV[m-1];
					}
	            
					m=m-2; 
					it=0;
				}
				else
				{ 
					if (it>=nMaxIt)
						return false;

					it=it+1;
					for (j=l+2; j<=m-1; j++)
						elements[j*n+j-2]=0.0;
					for (j=l+3; j<=m-1; j++)
						elements[j*n+j-3]=0.0;
					for (k=l; k<=m-2; k++)
					{ 
						if (k!=l)
						{ 
							p=elements[k*n+k-1]; 
							q=elements[(k+1)*n+k-1];
							r=0.0;
							if (k!=m-2) 
								r=elements[(k+2)*n+k-1];
						}
						else
						{ 
							x=elements[ii]+elements[ll];
							y=elements[ll]*elements[ii]-elements[kk]*elements[jj];
							ii=l*n+l; 
							jj=l*n+l+1;
							kk=(l+1)*n+l; 
							ll=(l+1)*n+l+1;
							p=elements[ii]*(elements[ii]-x)+elements[jj]*elements[kk]+y;
							q=elements[kk]*(elements[ii]+elements[ll]-x);
							r=elements[kk]*elements[(l+2)*n+l+1];
						}
	                
						if ((Math.Abs(p)+Math.Abs(q)+Math.Abs(r))!=0.0)
						{ 
							xy=1.0;
							if (p<0.0) 
								xy=-1.0;
							s=xy*Math.Sqrt(p*p+q*q+r*r);
							if (k!=l) 
								elements[k*n+k-1]=-s;
							e=-q/s; 
							f=-r/s; 
							x=-p/s;
							y=-x-f*r/(p+s);
							g=e*r/(p+s);
							z=-x-e*q/(p+s);
							for (j=k; j<=m-1; j++)
							{ 
								ii=k*n+j; 
								jj=(k+1)*n+j;
								p=x*elements[ii]+e*elements[jj];
								q=e*elements[ii]+y*elements[jj];
								r=f*elements[ii]+g*elements[jj];
								if (k!=m-2)
								{ 
									kk=(k+2)*n+j;
									p=p+f*elements[kk];
									q=q+g*elements[kk];
									r=r+z*elements[kk]; 
									elements[kk]=r;
								}
	                        
								elements[jj]=q; elements[ii]=p;
							}
	                    
							j=k+3;
							if (j>=m-1) 
								j=m-1;
	                    
							for (i=l; i<=j; i++)
							{ 
								ii=i*n+k; 
								jj=i*n+k+1;
								p=x*elements[ii]+e*elements[jj];
								q=e*elements[ii]+y*elements[jj];
								r=f*elements[ii]+g*elements[jj];
								if (k!=m-2)
								{ 
									kk=i*n+k+2;
									p=p+f*elements[kk];
									q=q+g*elements[kk];
									r=r+z*elements[kk]; 
									elements[kk]=r;
								}
	                        
								elements[jj]=q; 
								elements[ii]=p;
							}
						}
					}
				}
			}
	    
			return true;
		}

		/**
		 * 求實對稱矩陣特征值與特征向量的雅可比法
		 * 
		 * @param dblEigenValue - 一維數組，長度為矩陣的階數，返回時存放特征值
		 * @param mtxEigenVector - 返回時存放特征向量矩陣，其中第i列為與數組
		 *                         dblEigenValue中第j個特征值對應的特征向量
		 * @param nMaxIt - 迭代次數
		 * @param eps - 計算精度
		 * @return bool型，求解是否成功
		 */
		public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, int nMaxIt, double eps)
		{ 
			int i,j,p = 0,q = 0,u,w,t,s,l;
			double fm,cn,sn,omega,x,y,d;
	    
			if (! mtxEigenVector.Init(numColumns, numColumns))
				return false;

			l=1;
			for (i=0; i<=numColumns-1; i++)
			{ 
				mtxEigenVector.elements[i*numColumns+i]=1.0;
				for (j=0; j<=numColumns-1; j++)
					if (i!=j) 
						mtxEigenVector.elements[i*numColumns+j]=0.0;
			}
	    
			while (true)
			{ 
				fm=0.0;
				for (i=1; i<=numColumns-1; i++)
				{
					for (j=0; j<=i-1; j++)
					{ 
						d=Math.Abs(elements[i*numColumns+j]);
						if ((i!=j) && (d>fm))
						{ 
							fm=d; 
							p=i; 
							q=j;
						}
					}
				}

				if (fm<eps)
				{
					for (i=0; i<numColumns; ++i)
						dblEigenValue[i] = GetElement(i,i);
					return true;
				}

				if (l>nMaxIt)  
					return false;
	        
				l=l+1;
				u=p*numColumns+q; 
				w=p*numColumns+p; 
				t=q*numColumns+p; 
				s=q*numColumns+q;
				x=-elements[u]; 
				y=(elements[s]-elements[w])/2.0;
				omega=x/Math.Sqrt(x*x+y*y);

				if (y<0.0) 
					omega=-omega;

				sn=1.0+Math.Sqrt(1.0-omega*omega);
				sn=omega/Math.Sqrt(2.0*sn);
				cn=Math.Sqrt(1.0-sn*sn);
				fm=elements[w];
				elements[w]=fm*cn*cn+elements[s]*sn*sn+elements[u]*omega;
				elements[s]=fm*sn*sn+elements[s]*cn*cn-elements[u]*omega;
				elements[u]=0.0; 
				elements[t]=0.0;
				for (j=0; j<=numColumns-1; j++)
				{
					if ((j!=p) && (j!=q))
					{ 
						u=p*numColumns+j; w=q*numColumns+j;
						fm=elements[u];
						elements[u]=fm*cn+elements[w]*sn;
						elements[w]=-fm*sn+elements[w]*cn;
					}
				}

				for (i=0; i<=numColumns-1; i++)
				{
					if ((i!=p) && (i!=q))
					{ 
						u=i*numColumns+p; 
						w=i*numColumns+q;
						fm=elements[u];
						elements[u]=fm*cn+elements[w]*sn;
						elements[w]=-fm*sn+elements[w]*cn;
					}
				}

				for (i=0; i<=numColumns-1; i++)
				{ 
					u=i*numColumns+p; 
					w=i*numColumns+q;
					fm=mtxEigenVector.elements[u];
					mtxEigenVector.elements[u]=fm*cn+mtxEigenVector.elements[w]*sn;
					mtxEigenVector.elements[w]=-fm*sn+mtxEigenVector.elements[w]*cn;
				}
			}
		}

		/**
		 * 求實對稱矩陣特征值與特征向量的雅可比過關法
		 * 
		 * @param dblEigenValue - 一維數組，長度為矩陣的階數，返回時存放特征值
		 * @param mtxEigenVector - 返回時存放特征向量矩陣，其中第i列為與數組
		 *                         dblEigenValue中第j個特征值對應的特征向量
		 * @param eps - 計算精度
		 * @return bool型，求解是否成功
		 */
		public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, double eps)
		{ 
			int i,j,p,q,u,w,t,s;
			double ff,fm,cn,sn,omega,x,y,d;
	    
			if (! mtxEigenVector.Init(numColumns, numColumns))
				return false;

			for (i=0; i<=numColumns-1; i++)
			{ 
				mtxEigenVector.elements[i*numColumns+i]=1.0;
				for (j=0; j<=numColumns-1; j++)
					if (i!=j) 
						mtxEigenVector.elements[i*numColumns+j]=0.0;
			}
	    
			ff=0.0;
			for (i=1; i<=numColumns-1; i++)
			{
				for (j=0; j<=i-1; j++)
				{ 
					d=elements[i*numColumns+j]; 
					ff=ff+d*d; 
				}
			}

			ff=Math.Sqrt(2.0*ff);
			ff=ff/(1.0*numColumns);

			bool nextLoop = false;
			while (true)
			{
				for (i=1; i<=numColumns-1; i++)
				{
					for (j=0; j<=i-1; j++)
					{ 
						d=Math.Abs(elements[i*numColumns+j]);
						if (d>ff)
						{ 
							p=i; 
							q=j;

							u=p*numColumns+q; 
							w=p*numColumns+p; 
							t=q*numColumns+p; 
							s=q*numColumns+q;
							x=-elements[u]; 
							y=(elements[s]-elements[w])/2.0;
							omega=x/Math.Sqrt(x*x+y*y);
							if (y<0.0) 
								omega=-omega;
					    
							sn=1.0+Math.Sqrt(1.0-omega*omega);
							sn=omega/Math.Sqrt(2.0*sn);
							cn=Math.Sqrt(1.0-sn*sn);
							fm=elements[w];
							elements[w]=fm*cn*cn+elements[s]*sn*sn+elements[u]*omega;
							elements[s]=fm*sn*sn+elements[s]*cn*cn-elements[u]*omega;
							elements[u]=0.0; elements[t]=0.0;
					    
							for (j=0; j<=numColumns-1; j++)
							{
								if ((j!=p)&&(j!=q))
								{ 
									u=p*numColumns+j; 
									w=q*numColumns+j;
									fm=elements[u];
									elements[u]=fm*cn+elements[w]*sn;
									elements[w]=-fm*sn+elements[w]*cn;
								}
							}

							for (i=0; i<=numColumns-1; i++)
							{
								if ((i!=p)&&(i!=q))
								{ 
									u=i*numColumns+p; 
									w=i*numColumns+q;
									fm=elements[u];
									elements[u]=fm*cn+elements[w]*sn;
									elements[w]=-fm*sn+elements[w]*cn;
								}
							}
					    
							for (i=0; i<=numColumns-1; i++)
							{ 
								u=i*numColumns+p; 
								w=i*numColumns+q;
								fm=mtxEigenVector.elements[u];
								mtxEigenVector.elements[u]=fm*cn+mtxEigenVector.elements[w]*sn;
								mtxEigenVector.elements[w]=-fm*sn+mtxEigenVector.elements[w]*cn;
							}

							nextLoop = true;
							break;
						}
					}

					if (nextLoop)
						break;
				}
		        
				if (nextLoop)
				{
					nextLoop = false;
					continue;
				}

				nextLoop = false;

				// 如果達到精度要求，退出循環，返回結果
				if (ff<eps) 
				{
					for (i=0; i<numColumns; ++i)
						dblEigenValue[i] = GetElement(i,i);
					return true;
				}
		    
				ff=ff/(1.0*numColumns);
			}
		}
	}
}