算法导论LUP分解

#include<iostream>
#include<math.h>

using namespace std;

void Swap(double *a,double *b) {//交换函数 
	double temp=*a;
	*a=*b;
	*b=temp;
}
void Swap1(int *a,int *b) {//交换函数 
	int temp=*a;
	*a=*b;
	*b=temp;
}


double Round(double dVal, short iPlaces) {//四舍五入保留小数位数 
    double dRetval;
    double dMod = 0.0000001;
    if (dVal<0.0) dMod=-0.0000001;
    dRetval=dVal;
    dRetval+=(5.0/pow(10.0,iPlaces+1.0));
    dRetval*=pow(10.0,iPlaces);
    dRetval=floor(dRetval+dMod);
    dRetval/=pow(10.0,iPlaces);
    return(dRetval);
}


void Input(double **A,double *b,int n) {
	int i,j;
	cout<<"请输入A矩阵:"<<endl; 
	for(i=1;i<=n;i++) {
		for(j=1;j<=n;j++) {
			cin>>A[i][j];
		}
	}
	
	cout<<"请输入b数组:"<<endl;
	for(i=1;i<=n;i++) {
		cin>>b[i];
	} 
}

void LUP_DECOMPOSITION(double **A,double **L,double **U,int *P,int n) {
	int i,j,k,k1;
	double temp,temp1;;
	for(i=1;i<=n;i++) {
		P[i]=i;
	}
	
	for(k=1;k<=n;k++) {
		int p=0;
		for(i=k;i<=n;i++) {
			if(fabs(A[i][k])>p) {
				p=fabs(A[i][k]);
				k1=i;
			}
		}
		if(p==0) {
			cout<<"singular matrix"<<endl;
		} else {
			Swap1(&P[k],&P[k1]);
		}
		
		for(i=1;i<=n;i++) {
			Swap(&A[k][i],&A[k1][i]);
		} 
		
		for(i=k+1;i<=n;i++) {
			temp=A[i][k]/A[k][k];
			A[i][k]=Round(temp,2);
		for(j=k+1;j<=n;j++) {
			temp=A[i][k]*A[k][j];
			temp1=Round(temp,2);
			A[i][j]=A[i][j]-temp1;
		}
	}
}
        for(i=1;i<=n;i++) {
        	for(j=1;j<=n;j++) {
        		if(i==j) {
        			L[i][j]=1;
				} else {
					if(i<j) {
						L[i][j]=0;
					} else {
						L[i][j]=A[i][j];
					}
				}
			}
		}
		
		for(i=1;i<=n;i++) {
        	for(j=1;j<=n;j++) {
        		if(i<=j) {
        			U[i][j]=A[i][j];
			} else {
				U[i][j]=0;
			}
		}
	}	

}

void LUP_SOLVE(double **A,double **L,double **U,int *P,double *x,double *y,double *b,int n) {
	int i,j;
	double sum,temp;
	y[1]=b[P[1]];
	for(i=2;i<=n;i++) {
		sum=0.0;
		for(j=1;j<=i-1;j++) {
			temp=L[i][j]*y[j];
			sum+=Round(temp,2);
		}
		y[i]=b[P[i]]-sum;
	}
	
	temp=y[n]/U[n][n];
	x[n]=Round(temp,2);
	for(i=n-1;i>=1;i++) {
		sum=0.0;
		for(j=i+1;j<=n;j++) {
			temp=U[i][j]*x[j];
			sum+=Round(temp,2);
		}
		temp=(y[i]-sum)/U[i][i];
		x[i]=Round(temp,2);
	}
}

void Print(double *x,int n) {
	int i,j;
	for(i=1;i<=n;i++) {
		cout<<x[i]<<"\t";
	}
}

void PrintL(double **L,int n) {
	int i,j;
	for(i=1;i<=n;i++) {
		for(j=1;j<=n;j++) {
			cout<<L[i][j]<<"\t";
		}
		cout<<endl;
	}
}

void Delete(double **A,double **L,double **U,double *x,double *y,double *b,int *P,int n) {
	int i,j;
	for(i=1;i<=n;i++) {
		delete []L[i];
	}
	delete []L;
	
	for(i=1;i<=n;i++) {
		delete []U[i];
	}
	delete []U;
	
	for(i=1;i<=n;i++) {
		delete []A[i];
	}
	delete []A;
	
	delete []P;
	delete []x;
	delete []y;
	delete []b;
}
	
int main() {
	int n;
	double *x=NULL;
    double *y=NULL;
	double **L=NULL;
	double **U=NULL;
	int     *P=NULL;
	double **A=NULL;
	double *b=NULL;
	cin>>n;
	int i=0; 
	L=new double*[n+1];
	for(i=1;i<=n;i++) {
		L[i]=new double[n+1];
	}
	
	U=new double*[n+1];
	for(i=1;i<=n;i++) {
		U[i]=new double[n+1];
	}
	
	A=new double*[n+1];
	for(i=1;i<=n;i++) {
		A[i]=new double[n+1];
	}
	
	P=new int[n+1];
	x=new double[n+1];
	y=new double[n+1];
	b=new double[n+1];
	 
	Input(A,b,n);
	LUP_DECOMPOSITION(A,L,U,P,n);
	LUP_SOLVE(A,L,U,P,x,y,b,n);
	Print(x,n);
	PrintL(L,n);
	Delete(A,L,U,x,y,b,P,n);
	return 0;
}

/*
2 0 2 0.6
3 3 4 -2
5 5 4 2
-1 -2 3.4 -1

1 2 0
3 4 4
5 6 3
*/