//改進的平方根法解線性方程組
//Ax=b;  把A進行A=LDL’分解；L是一個單位下三角陣，D是一個對角陣；
//把Ax=b改寫為LDL'x=b,他等價于：
//    Ly=b;   DL'x=y;


#include<iostream.h>
#include<stdlib.h>
#define N 100
void main()
{
	cout<<"用改進的平方根法解線性方程組!\n"<<"Ax=b\n"<<"其中A是一個n*n的方陣，x和b是n維的列向量\n";
	int i,j,k;
	int n;   //n是A的維數
	double a[N][N],b[N],l[N][N],temp=0;
	
	cout<<"請輸入方陣A的階數n：";
	cin>>n;
	
	cout<<"請輸入方陣A：\n";
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			cin>>a[i][j];
		}
	
	cout<<"請輸入列向量b：\n";
    for(i=0;i<n;i++)
		cin>>b[i];
	if(a[0][0]==0)
	{
		cout<<"A的各介順序主子不全為非零；\n"<<"Ax=b無解";
		exit(0);
	}
	
	for(i=1;i<n;i++)
	{   
		for(j=0;j<i;j++)
		{   
			temp=0;
			for(k=0;k<j;k++)
				temp+=a[i][k]*l[j][k];
			a[i][j]=a[i][j]-temp;
			
			l[i][j]=a[i][j]/a[j][j];
		}
		temp=0;
		for(k=0;k<i;k++)
			temp+=a[i][k]*l[i][k];
		a[i][i]-=temp;
		if(a[i][i]==0)
		{
			cout<<"A的各介順序主子不全為非零；\n"<<"Ax=b無解";
		    exit(0);
		}
		
	}


	for(i=0;i<n;i++)
	{
		temp=0;
		for(k=0;k<i;k++)
			temp+=l[i][k]*b[k];
		b[i]=b[i]-temp;
	}
	for(i=n-1;i>=0;i--)
	{
		temp=0;
		for(k=i+1;k<n;k++)
			temp+=l[k][i]*b[k];
		b[i]=b[i]/a[i][i]-temp;
	}

	cout<<"Ax=b的解為：\n";
	for(i=0;i<n;i++)
		cout<<"x"<<i+1<<"="<<b[i]<<"\n";

}