/*************************    Program 4.2    ****************************/
/*                                                                      */
/************************************************************************/
/* 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 dtrm (a,n,d,indx)
/* Subroutine for evaluating the determinant of a matrix using 
   the partial pivoting Gaussian elimination scheme.
   Copyright (c) Tao Pang 1997. */
int n;
int indx[];
double *d;
double a[NMAX][NMAX];
{
int i,j,ir,msgn;
void elgs();

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

msgn = 1;
for (i = 0; i < n; ++i)
  {
  while (i /= indx[i])
    {
    msgn = -msgn;
    j = indx[i];
    indx[i] = indx[j];
    indx[j] = j;
    }
  }
*d = msgn*(*d);
}

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[i]][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];
      }
    }
  }
}