/*************************    Program 4.3    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 100

void legs (a,n,b,x,indx)
/* Subroutine for solving the equation A*X = B with the partial
   pivoting Gaussian elimination.  Copyright (c) Tao Pang 1997. */
int n;
int indx[];
double b[],x[];
double a[NMAX][NMAX];
{
int i,j,ir;
void elgs();

ir = 0;
elgs (a,n,ir,indx);
for(i = 0; i < n-1; ++i)
  {
  for(j = i+1; j < n; ++j)
    {
    b[indx[j]] = b[indx[j]]-a[indx[j]][i]*b[indx[i]];
    }
  }

x[n-1] = b[indx[n-1]]/a[indx[n-1]][n-1];
for (i = n-2; i>=0; i = i-1)
  {
  x[i] = b[indx[i]];
  for (j = i+1; j < n; ++j)
    {
    x[i] = x[i]-a[indx[i]][j]*x[j];
    }
  x[i] = x[i]/a[indx[i]][i];
  }
}

void elgs (a,n,ir,indx)
/* Subroutine for the partial pivoting Gaussian elimination.  A is
   the coefficient matrix in the input and the transformed matrix
   in the output.  INDX records the pivoting order.  IR is the
   indicator for the rescaling.  Copyright (c) Tao Pang 1997. */
int n,ir;
int indx[];
double a[NMAX][NMAX];
{
int i,j,k,itmp;
double c1,pi,pi1,pj;
double c[NMAX];

if (n > NMAX)
  {
  printf("The matrix dimension is too large.\n");
  exit(1);
  }

/* Initialize the index */

for (i = 0; i < n; ++i)
  {
  indx[i] = i;
  }

/* Rescale the coefficients */
 
for (i = 0; i < n; ++i)
  {
  c1 = 0;
  for (j = 0; j < n; ++j)
    {
    if (fabs(a[i][j]) > c1) c1 = fabs(a[i][j]);
    }
  c[i] = c1;
  }
for (i = 0; i < n; ++i)
  {
  for (j = 0; j < n; ++j)
    {
    a[i][j] = a[i][j]/c[i];
    }
  }

/* Search the pivoting elements */

for (j = 0; j < n-1; ++j)
  {
  pi1 = 0;
  for (i = j; i < n; ++i)
    {
    pi = fabs(a[indx[i]][j]);
    if (pi > pi1)
      {
      pi1 = pi;
      k = i;
      }
    }

/* Eliminate the elements below the diagonal */

  itmp = indx[j];
  indx[j] = indx[k];
  indx[k] = itmp;
  for (i = j+1; i < n; ++i)
    {
    pj = a[indx[i]][j]/a[indx[j]][j];
    a[indx[i]][j] = pj;
    for (k = j+1; k < n; ++k)
      {
      a[indx[i]][k] = a[indx[i]][k]-pj*a[indx[j]][k];
      }
    }
  }

/* Rescale the coefficients back to the original magnitudes */

if (ir == 1)
  {
  for (i = 0; i < n; ++i)
    {
    for (j = 0; j < n; ++j)
      {
      a[i][j] = c[i]*a[i][j];
      }
    }
  }
}