Code C/C++: Giải hệ phương trình tuyến tính dựa vào phân rã LU

Ý tưởng thuật toán: cho hệ phương trình tuyến tính tổng quát A.X=B. Ta tiến hành phân rã A=L.U. Trong đó, L là ma trận tam giác dưới và U là ma trận tam giác trên.
Khi đó,

Cài đặt chương trình:
#include <conio.h>
#include <iostream>
#define max 100
using namespace std;
/* Nhập ma trận hệ số */
void Nhap(float A[max][max],int n) {
for(int i = 0; i<n; i++)
for(int j = 0; j<n; j++){
cout<<"a["<<i<<"]["<<j<<"] = ";
cin>>A[i][j];
}
}
/* Nhập ma trận hệ số tự do */
void Nhap(float B[max],int n) {
for(int i = 0; i<n; i++) {
cout<<"b["<<i<<"] = ";
cin>>B[i];
}
}
/* Xuất ma trận hệ số tự do */
void Xuat(float B[max],int n) {
cout<<"(";
for(int i = 0; i<n-1; i++)
cout<<B[i]<<",";
cout<<B[n-1]<<")";
}
/* Xuất ma trận */
void Xuat(float A[max][max], int n) {
cout<<"\n";
for(int i=0 ; i<n; i++){
cout<<endl;
for(int j=0 ; j<n; j++)
cout<<A[i][j]<<"\t";
}
}
/* Xuất nghiệm */
void XuatNghiem(float X[], int n, char * s)  {
cout<<"\nNghiem cua he PTTT";
for(int i=0; i<n; i++)
cout<<s<<i+1<<"="<<X[i];
}
char HeTamGiacDuoi (float A[max][max], float X[max], float B[max], int n ) {
for(int i = 0; i<n; i++) {
if (A[i][i]!=0) {
if (i==0)
X[i] = B[i]/A[i][i];
else {
X[i] = B[i];
for(int j=0; j<i; j++)
X[i]=X[i]-A[i][j]*X[j];
X[i] = X[i]/A[i][i];
}
} else
return 0;
}
return 1;
}
char HeTamGiacTren (float A[max][max], float X[max], float B[max], int n ) {
for(int i = n-1; i>=0; i--) {
if (A[i][i]!=0) {
if (i==n-1)
X[i] = B[i]/A[i][i];
else {
X[i] = B[i];
for(int j=i+1; j<n; j++)
X[i]=X[i]-A[i][j]*X[j];
X[i] = X[i]/A[i][i];
}
} else
return 0;
}
return 1;
}
void PhanRaLU(float A[max][max], float L[max][max], float U[max][max], int n){
for(int k =0; k<n; k++) {
U[k][k] = A[k][k];
L[k][k] = 1;
for(int i=k+1; i<n; i++) { 
L[i][k] = A[i][k]/U[k][k];
U[k][i] = A[k][i];
U[i][k] = 0;
L[k][i] = 0;
}
for(int i = k+1; i<n; i++)
for(int j = k+1; j<n; j++)
A[i][j] = A[i][j]-L[i][k]*U[k][j];
}
}
/* Giải hệ phương trình tổng quát LUX=B*/
void GiaiHePTTT(float A[max][max], float X[max], float B[max], int n) {
float L[max][max],U[max][max], Y[max];
cout<<"Phan ra A = L.U\n";
PhanRaLU(A,L,U,n);
cout<<"Ma tran L";
Xuat(L,n);
cout<<"\nMa tran U";
Xuat(U,n);
cout<<"\nGiai LY = B. Nghiem Y";
if(HeTamGiacDuoi(L,Y,B,n)) {
XuatNghiem(Y,n,"\ny");
cout<<"\nGiai UX = Y. Nghiem X";
if(HeTamGiacTren(U,X,Y,n))
XuatNghiem(X,n,"\nx");
else
cout<<"\nHe pttt k co nghiem duy nhat";
} else
cout<<"\nHe pttt k co nghiem duy nhat";
}
int main() {
int n;
float A[max][max],B[max], X[max];
cout<<"Nhap so he phuong trinh n = "; cin>>n;
cout<<"Nhap ma tran he so A\n";
Nhap(A,n);
cout<<"Nhap ma tran he so tu do B\n";
Nhap(B,n);
GiaiHePTTT(A,X,B,n);
getch();
return 0;

} 
Kết quả chạy chương trình: