算法导论 矩阵运算 LUP

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void swopInt(int *a,int *b)
{
int temp=(*a);
(*a)=(*b);
(*b)=temp;
}

void swopDouble(double *a,double *b)
{
double temp=(*a);
(*a)=(*b);
(*b)=temp;
}

void LUPDecomposition(double **A,int n,int *pi)
{
for(int i=0;i<n;i++)
{
pi[i]=i;
}
for(int k=0;k<n;k++)
{
double p=0.0;
int kk=0;
for(int i=k;i<n;i++)
{
if(abs(A[i][k])>p)
{
p=abs(A[i][k]);
kk=i;
}
}
if(p==0)
{
printf("singular matrix");
return;
}
swopInt(&pi[kk],&pi[k]);
for(int i=0;i<n;i++)
swopDouble(&A[k][i],&A[kk][i]);
for(int i=k+1;i<n;i++)
{
A[i][k]=A[i][k]/A[k][k];
for(int j=k+1;j<n;j++)
{
A[i][j]=A[i][j]-A[i][k]*A[k][j];
}
}
}
}

void printA(double **A,int n)
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
printf("%.1f ",A[i][j]);
}
printf("\n");
}
}

double* LUPSolve(double **L,double **U,int *pi,double *b,int n)
{
double *x=(double*)malloc(n*sizeof(double));
double *y=(double*)malloc(n*sizeof(double));
for(int i=0;i<n;i++)
{
double sum=0.0;
for(int j=0;j<=i-1;j++)
{
sum+=L[i][j]*y[j];
}
y[i]=b[pi[i]]-sum;
}
for(int i=n-1;i>=0;i--)
{
double sum=0.0;
for(int j=i+1;j<n;j++)
{
sum+=U[i][j]*x[j];
}
x[i]=(y[i]-sum)/U[i][i];
}
return x;
}

void generateLU(double **A,double **L,double **U,int n)
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
if(i==j)
L[i][j]=1;
if(i>j)
{
L[i][j]=A[i][j];
U[i][j]=0;
}
else
{
U[i][j]=A[i][j];
if(i!=j)
L[i][j]=0;
}
}
}
}

void printX(double *x,int n)
{
for(int i=0;i<n;i++)
{
printf("%.1f ",x[i]);
}
printf("\n");
}

void main()
{
int n=4;
double **A=(double**)malloc(n*sizeof(double*));
double **L=(double**)malloc(n*sizeof(double*));
double **U=(double**)malloc(n*sizeof(double*));
for(int i=0;i<n;i++)
{
A[i]=(double*)malloc(n*sizeof(double));
L[i]=(double*)malloc(n*sizeof(double));
U[i]=(double*)malloc(n*sizeof(double));
}
int *pi=(int*)malloc(n*sizeof(int));
A[0][0]=2;
A[0][1]=0;
A[0][2]=2;
A[0][3]=0.6;
A[1][0]=3;
A[1][1]=3;
A[1][2]=4;
A[1][3]=-2;
A[2][0]=5;
A[2][1]=5;
A[2][2]=4;
A[2][3]=2;
A[3][0]=-1;
A[3][1]=-2;
A[3][2]=3.4;
A[3][3]=-1;
double *b=(double*)malloc(n*sizeof(double));
b[0]=4.6;
b[1]=8;
b[2]=16;
b[3]=-0.6;
LUPDecomposition(A,n,pi);
generateLU(A,L,U,n);
double *x=LUPSolve(L,U,pi,b,n);
printA(U,n);
printX(x,n);
getchar();
}