cholesky.cpp

乔列斯基分解的C语言实现

#include<stdio.h>
#include<math.h>
#define ZERO 1.0E-10
#define true 1
#define false 0

void INPUT(int *, double A[][12], int *);
void OUTPUT(double L[][11], double A[][12], int N);
void BackX(double L[][11],double* X,int N,int FG);

main()
{
   double A[10][12], L[10][11],X[10];
   double sum;
   int i,j,k;
   int OK;
   int N;

   INPUT(&OK, A, &N);
   if(OK)
   {
	   for (i=1; i<=N; i++)
	   {
		   X[i-1] = 0.0;
		   for(j=1; j<=N+2; j++) L[i-1][j-1] = 0.0;//L、X阵初始化
	   }
   }
//将A矩阵分解为L、D矩阵
   for(i=1; i<=N; i++)
   {
	   for(j=1; j<=i; j++)
	   {
		   L[i-1][0] = A[i-1][0] / L[0][0];
		   L[0][0] = A[0][0];

		   if(i==j)//求解D矩阵
		   {
			   for(k=1,sum=0.0; k<=i-1; k++)
			   {
				   sum = sum + L[k-1][k-1]*L[i-1][k-1]*L[i-1][k-1]; 
				   L[i-1][i-1] = A[i-1][i-1] - sum;
			   }
		   }
		   else//求解L矩阵
		   {
			   for(k=1,sum=0.0; k<=j-1; k++)
			   {
				   sum = sum + L[k-1][k-1] * L[i-1][k-1] * L[j-1][k-1];
				   L[i-1][j-1] = (A[i-1][j-1] - sum) / L[j-1][j-1];
			   }
		   }
	   }
   }
//计算B矩阵第一列解向量
	for(i=1; i<=N; i++) L[i-1][N]=A[i-1][N];
	BackX(L, X, N, 1);
	for(i=1; i<=N; i++) L[i-1][N]=X[i-1];
	BackX(L, X, N, 2);
	for(i=1; i<=N; i++) L[i-1][N]=X[i-1];
	BackX(L, X, N, 3);
	for(i=1; i<=N; i++) A[i-1][N]=X[i-1];
//计算B矩阵第二列解向量
	for(i=1; i<=N; i++) L[i-1][N]=A[i-1][N+1];
	BackX(L, X, N, 1);
	for(i=1; i<=N; i++) L[i-1][N]=X[i-1];
	BackX(L, X, N, 2);
	for(i=1; i<=N; i++) L[i-1][N]=X[i-1];
	BackX(L, X, N, 3);
	for(i=1; i<=N; i++) A[i-1][N+1]=X[i-1];

	OUTPUT(L, A, N);

	return 0;
}

//回代求解函数
//L 三角矩阵，X-解向量矩阵，N-向量维数，FG-L、D、LT矩阵判断标志
void BackX(double L[][11],double* X,int N,int FG)
{
	int i,k;
	double sum = 0.0;

	if(FG==1)//解Lx=B形式
	{
		X[0] = L[0][N];
		for(i=1; i<=N; i++)
			for(k=1,sum=0.0; k<i; k++)
			{				
				sum = sum + L[i-1][k-1]*X[k-1];
				X[i-1] = L[i-1][N] - sum;
			}
	}
	else if(FG==2)//解Dx=B形式
	{
		for(i=1; i<=N; i++)
			X[i-1] = L[i-1][N] / L[i-1][i-1];
	}
	else//解LTx=B形式
	{
		X[N-1] = L[N-1][N];
		for(i=N; i>=1; i--)
			for(k=N,sum=0.0; k>i; k--)
			{
				sum = sum + L[k-1][i-1]*X[k-1];
				X[i-1] = L[i-1][N] - sum;
			}
	}
}


void INPUT(int *OK, double A[][12], int *N)
{
   int I, J, FLAG;
   char AA;
   char NAME[30];
   FILE *INP; 

   printf("乔累斯基分解法. \n");
   printf("系数阵请按以下格式写在txt文档中: \n");
   printf("a(1,1),a(1,2),...,b(1,1),b(1,2),a(2,1),a(2,2),...,b(2,1),b(2,2),...,a(N,1),a(N,2),...,b(N,1),b(N,2)\n");
   printf("所有元素都写在同一行内，元素之间用空格隔开.\n");
   printf("txt文档是否建好? - 键入 Y or N.\n");
   scanf("%c",&AA);
   if ((AA == 'Y') || (AA == 'y')) 
   {
      printf("输入文档名称，格式为 - drive:name.txt \n");
      printf("例如:   A:DATA.DTA\n");
      scanf("%s", NAME);
      INP = fopen(NAME, "r");
      *OK = false;
      while (!(*OK)) 
	  {
         printf("输入方程维数：.\n");
         scanf("%d", N);
         if (*N > 0) 
		 {
            for (I=1; I<=*N; I++) 
			{
               for (J=1; J<=*N+2; J++) fscanf(INP, "%lf", &A[I-1][J-1]);
               fscanf(INP, "\n");
            }
            *OK = true;
            fclose(INP);
         }
         else printf("维数必须为正整数.\n");
      }
   }
   else 
   {
      printf("没有找到该文档.\n");
      *OK = false;
   }
}


void OUTPUT(double L[][11], double A[][12], int N)
{
   int I, J;
   FILE *OUP;

   OUP = stdout;
   fprintf(OUP, "\n\n结果输出：\n");
   fprintf(OUP, "\n解向量1为:\n");
   for (I=1; I<=N; I++) 
   {
      fprintf(OUP, "  %12.8f", A[I-1][N]);
   }

   fprintf(OUP, "\n解向量2为:\n");
   for (I=1; I<=N; I++) 
   {
      fprintf(OUP, "  %12.8f", A[I-1][N+1]);
   }
   fprintf(OUP, "\n");

   fprintf(OUP, "\nLD矩阵合写为:\n");
   for (I=1; I<=N; I++) 
   {
	   for(J=1; J<=N; J++)
	   {
		   fprintf(OUP, "  %12.8f", L[I-1][J-1]);
	   }
	   fprintf(OUP,"\n");
   }

   fclose(OUP);
}