/* Updated 10/24/2001. */

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

  elgs(a,n,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,indx)

/* Function to perform the partial-pivoting Gaussian elimination.
   a[][] is the original matrix in the input and transformed
   matrix plus the pivoting element ratios below the diagonal
   in the output.  indx[] records the pivoting order.
   Copyright (c) Tao Pang 2001. */

int n;
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;
  }

/* Find the rescaling factors, one from each row */
 
  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;
  }

/* Search the pivoting (largest) element from each column */ 

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

/* Interchange the rows via indx[] to record pivoting order */

    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];

/* Record pivoting ratios below the diagonal */

      a[indx[i]][j] = pj;

/* Modify other elements accordingly */

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