leqs.cs

对称正定矩阵LDLT分解算法（带应用示例）。

using System;

namespace 对称正定矩阵LDLT分解
{
	class LEQS
	{
		#region LDLT改进平方根法分解对称正定矩阵
		// N----矩阵的阶数
		// A----待分解和分解后的对称正定矩阵（下三解一维存贮）
		public static void LDLT(int N,double[] A)
		{
			double[] Tmp = new double[N];			//中间临时数组
			for(int i=1;i<N;i++)								//对每行进行计算
			{
				int j=0;							//列号
				int ii=i*(i+1)/2;				//元素在数组中的序号
				double G=0;
				while(j<i)						//对该行的每列元素计算（不含对角元素）
				{
					Tmp[j] = A[ii];
					int k=j*(j+3)/2;			//所在列的对角元素序号
					double E=0;
					for(int jj=1; jj<=j; jj++)	E+=Tmp[j-jj]*A[k-jj];
					Tmp[j] -= E;
					A[ii] = Tmp[j]/A[k];
					G += Tmp[j]*A[ii];
					j++; ii++;
				}
				A[ii] -= G;		//对角元素计算
			}
		}
		#endregion LDLT改进平方根法分解对称正定矩阵

		#region 利用LDLT求解方程组
		// N----系数矩阵的阶数
		// A----分解后的系数矩阵（下三解一维存贮）
		// B----左端向量及求解结果
		public static void LEQS_LDLT(int N,double[] A,double[] B)
		{
			//下三角求解
			for(int i=0;i<N;i++)
			{
				int ii=i*(i+1)/2;
				for(int j=0;j<i;j++) B[i]-=A[ii+j]*B[j];
			}
			//主对角求解
			for(int i=0;i<N;i++) B[i]/=A[(i+1)*(i+2)/2-1];
			//上三角求解
			for(int i=N-1;i>=0;i--)
			{
				int ii=i*(i+1)/2;
				for(int j=0;j<i;j++) B[j]-=A[ii+j]*B[i];
			}
		}
		#endregion 利用LDLT求解方程组

		#region 矩阵乘法
		public static void MMult(int N,double[,] A,double[,] B,double[,] C)
		{
			Array.Clear(C,0,N*N);
			for(int i=0;i<N;i++) for(int j=0;j<N;j++) for(int k=0;k<N;k++) C[i,j]+=A[i,k]*B[k,j];
		}
		#endregion 矩阵乘法

		#region 测试用例
		static void Main()
		{
			/*例1
			int N=3;
			double[] A = new double[]
			{
				 4,	
				 2,  2,	
				-2, -3, 14
			};
			double[] B=new double[]{2, -3, 34};
			/*例1*/

			//*例2
			int N=4;
			double[] A = new double[]
			{
				5, 
				-4,  6, 
				1, -4,  6, 
				0,  1, -4,  5 
			};
			double[] B=new double[]{0, 0, -5, 10};
			/*例2*/

			//LDLT分解
			LDLT(N,A);
			for(int i=0;i<N;i++)
			{
				for(int j=i*(i+1)/2;j<(i+1)*(i+2)/2;j++) Console.Write("{0:n2}\t",A[j]);
				Console.WriteLine();
			}
			Console.WriteLine();
			//还原
			double[,] L=new double[N,N];
			double[,] D=new double[N,N];
			double[,] LT=new double[N,N];
			for(int i=0;i<N;i++)
			{
				D[i,i]=A[(i+1)*(i+2)/2-1];
				LT[i,i]=L[i,i]=1.0;
				for(int j=0;j<i;j++) LT[j,i]=L[i,j]=A[i*(i+1)/2+j];
			}
			double[,] LD=new double[N,N];	MMult(N,L,D,LD);
			double[,] LDLt=new double[N,N];	MMult(N,LD,LT,LDLt);
			for(int i=0;i<N;i++)
			{
				for(int j=0;j<N;j++) Console.Write("{0:n2}\t",LDLt[i,j]);
				Console.WriteLine();
			}
			Console.WriteLine();
			LEQS_LDLT(N,A,B);
			for(int i=0;i<N;i++) Console.Write("{0:n2}\t",B[i]);
			Console.WriteLine();
		}
		#endregion 测试用例
	}
}