#  include "stdio.h"
#  include "math.h"
#  include <iostream.h>
#  include "fstream.h"
#  include "stdlib.h"
#  define MAX 200
#  define N   200
double U[N][N],L[N][N];
void InputData(int *p,double a[MAX][MAX],double b[MAX])
{
	int i,j;
	ifstream infile("data.txt",ios::in|ios::nocreate);
	if(!infile)
	{
		cout<<"ERROR:不能打開輸入文件!";
		exit(1);
	}
	infile>>*p;
    printf("從Data.txt中讀入的方程的階數(shù)為%d:\n",*p);
	printf("方程的系數(shù):\n");

    for(i=0;i<*p;i++)
	{	
		for(j=0;j<*p;j++)
		{
           infile>>a[i][j];
		   printf("%lf  ",a[i][j]);
        }
		printf("\n");
	}
    printf("向量b:\n");
    for(i=0;i<*p;i++)
	{
		infile>>b[i];
		printf("%lf  ",b[i]);
    }
    printf("\n");
}

void lut(int n,double a[MAX][MAX],double b[MAX])
{   
    double y[N];
    double x[N];
    int i,j,k;
    double sum;

	for(i=0;i<N;i++)
	for(j=0;j<N;j++)
	{
			L[i][j]=0;
			U[i][j]=0;
	}
    if (a[0][0]==0) { printf("No answer"); getchar();}
    for(j=0;j<n;j++)
        U[0][j]=a[0][j];

    for(j=1;j<n;j++)
        L[j][0]=a[j][0]/U[0][0];
 for(i=1;i<n;i++)
 {
   sum=0;
   for(k=0;k<i;k++) sum=sum+L[i][k]*U[k][i];
   U[i][i]=a[i][i]-sum;
   if(U[i][i]==0)  { printf("No answer"); getchar();}
        for(j=i+1;j<n;j++)
        {
             sum=0;
             for(k=0;k<i;k++) sum=sum+L[i][k]*U[k][j];
             U[i][j]=a[i][j]-sum;
             sum=0;
             for(k=0;k<i;k++) sum=sum+L[j][k]*U[k][i];
             L[j][i]=(a[j][i]-sum)/U[i][i];
        }

 }
   for (i=0;i<n;i++)
   {         
	         L[i][i]=1;
             y[0]=b[0]/L[0][0];
             sum=0  ;
             for (k=0;k<i;k++)   sum=sum+L[i][k]*y[k];
             y[i]=(b[i]-sum)/L[i][i];
    }

   x[n-1]=y[n-1]/U[n-1][n-1];

      for (i=1;i<n;i++)
      {  
	     k=n-1-i;
         sum=0;
         for(j=k+1;j<n+1;j++)
             sum=sum+U[k][j]*x[j];
         x[k]=(y[k]-sum)/U[k][k];
      }
	  //printf("\nx[N]:\n");
	  for (i=0;i<n;i++)
      {		  b[i]=x[i];
              //printf("%lf ",x[i]);
	  }
}



void main()
{
int xuanze;
int i,j,k,n;
double c;
double a[MAX][MAX];
double b[MAX];
int state=0;
double lanmuda,d,e,Nmax;
double r[MAX],x[MAX];
memu:
printf("                        ╭═══════════════╮\n");
printf("                        ║       LU分解求解方程         ║\n");
printf("    ╭═════════┤                              ├═════════╮\n");
printf("    ║                  ║     沈紅麗,高蓮,胡玉貴       ║                  ║\n ");
printf("   ║                  ╰═══════════════╯                  ║\n ");
printf("   ║                                                                      ║\n");
printf("    ║                【 1 】.利用LU分解法求方程的解(Data.txt)              ║\n");
printf("    ║                【 2 】.計算希爾伯特矩陣的解                          ║\n");
printf("    ║                【 3 】.退出程序                                      ║\n"); 
printf("    ║                  (輸入數(shù)字選擇相應程序)                              ║\n");
printf("    ╰═══════════════════════════════════╯\n");
printf("請輸入選擇:");
scanf("%d",&xuanze);
if(xuanze==1)
{

    //printf("輸入系數(shù)矩陣階數(shù)n:");
    //scanf("%d",&n);

    //if(n>MAX)                        ｜
     //{printf("超出范圍\n");while(1);}
    //else
     
    /* printf("輸入系數(shù):\n");
     for(i=0;i<n;i++)
       for(j=0;j<n;j++)
       scanf("%lf",&a[i][j]);

      printf("輸入向量b:\n");
      for(i=0;i<n;i++)
      scanf("%lf",&b[i]);



     printf("輸入的方程組為:\n");
     for(i=0;i<n;i++)
     {
       for(j=0;j<n;j++)
         printf(" %lf  ",a[i][j]);
         printf("    %lf \n",b[i]);
     }*/


	 InputData(&n,a,b);
     lut(n,a,b);
     printf("\n分解后的L矩陣:\n");
     for(i=0;i<n;i++)
     {
       for(j=0;j<n;j++)
           printf(" %lf  ",L[i][j]);
       printf("\n");
     }

	 printf("\n分解后的U矩陣:\n");
     for(i=0;i<n;i++)
     {
       for(j=0;j<n;j++)
           printf(" %lf  ",U[i][j]);
       printf("\n");
     }
     printf("\n輸入方程組的解為:\n");
	 for(i=0;i<n;i++)
	 {
		 printf("x[%d]=%lf\n",i,b[i]);
	 }
	 goto memu;

   

}




else if(xuanze==2)
{


   printf("輸入hilbert矩陣的階數(shù)n: ");
   scanf("%d",&n);
   d=1/(double)n/2;

   printf("希爾伯特矩陣系數(shù)a:\n");
   for(i=0;i<n;i++)
   {  b[i]=0;
      for(j=0;j<n;j++)
      {a[i][j]=1/(double)(i+j+2);
       b[i]=b[i]+1/(double)(i+j+2);
       printf("  %lf",a[i][j]);
	   if((i*n+j+1)%7==0)printf("\n");
      }
      
   }
   printf("\n希爾伯特矩陣系數(shù)b:\n");
   for(i=0;i<n;i++)
   {   
	  printf("  %lf",b[i]);
      if((i+1)%7==0)printf("\n");
   }
   //printf("shu ru lanmuda( 0<lanmuda<=%lf )\n",d*d/9);
   //scanf("%lf",&lanmuda);
   lanmuda=d*d/9;
   printf("\n輸入最大迭代次數(shù)Nmax: ");
   scanf("%lf",&Nmax);
   printf("輸入計算精度e: ");
   scanf("%lf",&e);

   for(i=0;i<n;i++)
        a[i][i]=a[i][i]+lanmuda ;
   //liezhuyuan(n,a,b);
   lut(n,a,b);
   for(i=0;i<n;i++)
      x[i]=b[i];

   for(k=1;k<=Nmax;k++)
   {


   for(i=0;i<n;i++)
   {  b[i]=0;
      for(j=0;j<n;j++)
      {
	   a[i][j]=1/(double)(i+j+2);
       b[i]=b[i]+1/(double)(i+j+2);
      }
   }

      for(i=0;i<n;i++)
      {
         r[i]=b[i];
         for(j=0;j<n;j++)
         r[i]=r[i]-a[i][j]*x[j];
       /*  printf("  %lf  \n",r[i]);  */

      }
      for(i=0;i<n;i++)
        a[i][i]=a[i][i]+lanmuda;
      lut(n,a,b);

      c=0;
      for(i=0;i<n;i++)
      {
      x[i]=x[i]+r[i];
      if(c<fabs(r[i]))
      c=fabs(r[i]);
      }
      if(c<e)
        {
		  printf("\n計算結(jié)果:\n");
		  printf("======================================================================\n");
          for(i=0;i<n;i++)
		  {
               
			  printf("x[%3d]=%lf   ",i+1,x[i]);
			  if((i+1)%4==0)printf("\n");
          }

		  printf("\n=====================================================================\n");
		  //printf("\n迭代的次數(shù)為%d\n",k);
        //while(1);
        break;
        }
        
   }
   //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
   goto memu; 
   getchar();
}
else 
{
	exit(1);//退出;
}
 
}