全选主元高斯消去算法

/*
全选主元高斯消去算法,求解线性方程组 A * x = b
By Kingwei 2005.2.27
*/

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

#define EPS 1e-11

#define MAX_DIM 20

int GaussLEquation(double x[], int dim, double A[][MAX_DIM], double b[])
{
int maxColPos[MAX_DIM], maxRowPos;
double maxElem, absElem, temp;
int k, i, j, status = 1;

for (k=0; k<dim; k++)
{
maxElem = 0.0;
for (i=k; i<dim; i++)
for (j=k; j<dim; j++)
{
absElem = fabs(A[i][j]);
if (maxElem < absElem)
{
maxElem = absElem;
maxRowPos = i;
maxColPos[k] = j;
}
}

if (maxElem <= EPS)
{
for (i=k; i<dim; i++)
if (fabs(b[i]) > EPS)
return 0;
for (i=k; i<dim; i++)
maxColPos[i] = i;
status = 2;
break;
}

if (maxRowPos != k)
{
for (j=k; j<dim; j++)
{
temp = A[k][j];
A[k][j] = A[maxRowPos][j];
A[maxRowPos][j] = temp;
}
temp = b[k];
b[k] = b[maxRowPos];
b[maxRowPos] = temp;
}

if (maxColPos[k] != k)
{
for (i=0; i<dim; i++)
{
temp = A[i][k];
A[i][k] = A[i][maxColPos[k]];
A[i][maxColPos[k]] = temp;
}
}

temp = A[k][k];
for (j=k; j<dim; j++)
A[k][j] /= temp;
b[k] /= temp;

for (i=k+1; i<dim; i++)
{
temp = A[i][k];
for (j=k; j<dim; j++)
A[i][j] -= A[k][j]*temp;
b[i] -= b[k]*temp;
}
}

x[dim-1] = b[dim-1];
for (i=dim-1; i>=0; i--)
{
temp = 0.0;
for (j=dim-1; j>i; j--)
temp += x[j]*A[i][j];
x[i] = b[i]-temp;
}

for (i=dim-1; i>=0; i--)
if (maxColPos[i] != i)
{
temp = x[i];
x[i] = x[maxColPos[i]];
x[maxColPos[i]] = temp;
}

return status;
}

int ReadLEquation(int *dim, double A[][MAX_DIM], double b[])
{
int i, j;
if (scanf("%d", dim) == EOF)
return 0;
for (i=0; i<*dim; i++)
{
for (j=0; j<*dim; j++)
scanf("%lf", &A[i][j]);
scanf("%lf", &b[i]);
}
return 1;
}

void DisplayRoot(double x[], int dim)
{
int i;
for (i=0; i<dim; i++)
printf("%lf ", x[i]);
printf("/n/n");
}

int main()
{
int dim;
double A[MAX_DIM][MAX_DIM], b[MAX_DIM], x[MAX_DIM];

while (ReadLEquation(&dim, A, b))
{
switch (GaussLEquation(x, dim, A, b))
{
case 0: printf("No Solution!/n/n");
break;
case 1: printf("Succeseful Solved!/n");
DisplayRoot(x, dim);
break;
case 2: printf("Unnumberable Solutions!/n");
DisplayRoot(x, dim);
break;
default:break;
}
}

return 0;
}