对角线方程组.cs

求解三对角线方程组的追赶法 return bool 型

/**
		 * 求解三对角线方程组的追赶法
		 * 
		 * @param mtxResult - Matrix对象，返回方程组解矩阵
		 * @return bool 型，方程组求解是否成功
		 */
		public bool GetRootsetTriDiagonal(Matrix mtxResult)
		{ 
			int k,j;
			double s;
	    
			// 将常数矩阵赋给解矩阵
			mtxResult.SetValue(mtxLEConst);
			double[] pDataConst = mtxResult.GetData();

			int n = GetNumUnknowns();
			if (mtxLECoef.GetNumRows() != n)
				return false;

			// 为系数矩阵对角线数组分配内存
			double[] pDiagData = new double[3*n-2];

			// 构造系数矩阵对角线元素数组
			k = j = 0;
			if (k == 0)
			{
				pDiagData[j++] = mtxLECoef.GetElement(k, k);
				pDiagData[j++] = mtxLECoef.GetElement(k, k+1);
			}
			for (k=1; k<n-1; ++k)
			{
				pDiagData[j++] = mtxLECoef.GetElement(k, k-1);
				pDiagData[j++] = mtxLECoef.GetElement(k, k);
				pDiagData[j++] = mtxLECoef.GetElement(k, k+1);
			}
			if (k == n-1)
			{
				pDiagData[j++] = mtxLECoef.GetElement(k, k-1);
				pDiagData[j++] = mtxLECoef.GetElement(k, k);
			}

			// 求解
			for (k=0;k<=n-2;k++)
			{ 
				j=3*k; 
				s=pDiagData[j];

				// 求解失败
				if (Math.Abs(s)+1.0==1.0)
				{
					return false;
				}

				pDiagData[j+1]=pDiagData[j+1]/s;
				pDataConst[k]=pDataConst[k]/s;
				pDiagData[j+3]=pDiagData[j+3]-pDiagData[j+2]*pDiagData[j+1];
				pDataConst[k+1]=pDataConst[k+1]-pDiagData[j+2]*pDataConst[k];
			}
	    
			s=pDiagData[3*n-3];
			if (s == 0.0)
			{
				return false;
			}
	    
			// 调整
			pDataConst[n-1]=pDataConst[n-1]/s;
			for (k=n-2;k>=0;k--)
				pDataConst[k]=pDataConst[k]-pDiagData[3*k+1]*pDataConst[k+1];
	    
			return true;
		}

		/**
		 * 一般带型方程组的求解
		 * 
		 * @param nBandWidth - 带宽
		 * @param mtxResult - Matrix对象，返回方程组解矩阵
		 * @return bool 型，方程组求解是否成功
		 */
		public bool GetRootsetBand(int nBandWidth, Matrix mtxResult)
		{ 
			int ls,k,i,j,nis=0,u,v;
			double p,t;
	    
			// 带宽必须为奇数
			if ((nBandWidth-1)%2 != 0)
				return false;

			// 将常数矩阵赋给解矩阵
			mtxResult.SetValue(mtxLEConst);
			double[] pDataConst = mtxResult.GetData();

			// 方程组属性
			int m = mtxLEConst.GetNumColumns();
			int n = GetNumUnknowns();
			if (mtxLECoef.GetNumRows() != n)
				return false;

			// 带宽数组：带型矩阵的有效数据
			double[] pBandData = new double[nBandWidth*n];

			// 半带宽
			ls = (nBandWidth-1)/2;

			// 构造带宽数组
			for (i=0; i<n; ++i)
			{
				j = 0;
				for (k=Math.Max(0, i-ls); k<Math.Max(0, i-ls)+nBandWidth; ++k)
				{
					if (k < n)
						pBandData[i*nBandWidth+j++] = mtxLECoef.GetElement(i, k);
					else
						pBandData[i*nBandWidth+j++] = 0;
				}
			}

			// 求解
			for (k=0;k<=n-2;k++)
			{ 
				p=0.0;
				for (i=k;i<=ls;i++)
				{ 
					t=Math.Abs(pBandData[i*nBandWidth]);
					if (t>p) 
					{
						p=t; 
						nis=i;
					}
				}
	        
				if (p == 0.0)
				{
					return false;
				}

				for (j=0;j<=m-1;j++)
				{ 
					u=k*m+j; 
					v=nis*m+j;
					t=pDataConst[u]; 
					pDataConst[u]=pDataConst[v]; 
					pDataConst[v]=t;
				}
	        
				for (j=0;j<=nBandWidth-1;j++)
				{ 
					u=k*nBandWidth+j; 
					v=nis*nBandWidth+j;
					t=pBandData[u]; 
					pBandData[u]=pBandData[v]; 
					pBandData[v]=t;
				}
	        
				for (j=0;j<=m-1;j++)
				{ 
					u=k*m+j; 
					pDataConst[u]=pDataConst[u]/pBandData[k*nBandWidth];
				}
	        
				for (j=1;j<=nBandWidth-1;j++)
				{ 
					u=k*nBandWidth+j; 
					pBandData[u]=pBandData[u]/pBandData[k*nBandWidth];
				}
	        
				for (i=k+1;i<=ls;i++)
				{ 
					t=pBandData[i*nBandWidth];
					for (j=0;j<=m-1;j++)
					{ 
						u=i*m+j; 
						v=k*m+j;
						pDataConst[u]=pDataConst[u]-t*pDataConst[v];
					}
	            
					for (j=1;j<=nBandWidth-1;j++)
					{ 
						u=i*nBandWidth+j; 
						v=k*nBandWidth+j;
						pBandData[u-1]=pBandData[u]-t*pBandData[v];
					}
	            
					u=i*nBandWidth+nBandWidth-1; pBandData[u]=0.0;
				}
	        
				if (ls!=(n-1)) 
					ls=ls+1;
			}
	    
			p=pBandData[(n-1)*nBandWidth];
			if (p == 0.0)
			{
				return false;
			}

			for (j=0;j<=m-1;j++)
			{ 
				u=(n-1)*m+j; 
				pDataConst[u]=pDataConst[u]/p;
			}
	    
			ls=1;
			for (i=n-2;i>=0;i--)
			{ 
				for (k=0;k<=m-1;k++)
				{ 
					u=i*m+k;
					for (j=1;j<=ls;j++)
					{ 
						v=i*nBandWidth+j; 
						nis=(i+j)*m+k;
						pDataConst[u]=pDataConst[u]-pBandData[v]*pDataConst[nis];
					}
				}
	        
				if (ls!=(nBandWidth-1)) 
					ls=ls+1;
			}
	    
			return true;
		}