/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-     ANSI C       Numerical Methods Toolbox                               -*/
/*-     Copyright (c) 1986, 87 by Borland International, Inc.                -*/
/*-                                                                          -*/
/*-  This unit provides functions for dealing with systems of linear         -*/
/*-  equations.                                                              -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

#define TNNearlyZero 1E-15

#define TNArraySize 30               /* Size of the matrix */

typedef double *TNvector;
typedef double **TNmatrix;

double _abs(x)
double x;
{
  return (x < 0) ? -x : x;
}

void Determinant();
/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-     Input: Dimen, Data                                                   -*/
/*-     Output: Det, Error                                                   -*/
/*-                                                                          -*/
/*-             Purpose : Calculate the determinant of a matrix by           -*/
/*-                       making it upper-triangular && )                -*/
/*-                       taking the product of the diagonal elements.       -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real;          -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen : int;  Dimension of the square matrix   -*/
/*-                       Data  : TNmatrix; Square matrix                    -*/
/*-                       Det   : real;     Determinant of Data              -*/
/*-                       Error : int;  Flags if ( something goes wrong    -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Dimen < 1                                       -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void Inverse();

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Data                                        -*/
/*-               Output: Inv, Error                                         -*/
/*-                                                                          -*/
/*-             Purpose : calculate the inverse of a matrix with             -*/
/*-                       Gauss-Jordan elimination.                          -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real;          -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen : int;   Dimension of the square matrix  -*/
/*-                       Data  : TNmatrix;  Square matrix                   -*/
/*-                       Inv   : TNmatrix;  Inverse of Data                 -*/
/*-                       Error : int;   Flags if ( something goes wrong   -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Dimen < 1                                       -*/
/*-                       2: no inverse exists                               -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void Gaussian_Elimination();

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Coefficients, Constants                     -*/
/*-               Output: Solution, Error                                    -*/
/*-                                                                          -*/
/*-             Purpose : Calculate the solution of a linear set of          -*/
/*-                       equations using Gaussian elimination &&           -*/
/*-                       backwards substitution.                            -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real           -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen : int;         Dimension of the square   -*/
/*-                                                matrix                    -*/
/*-                       Coefficients : TNmatrix; Square matrix             -*/
/*-                       Constants    : TNvector; Constants of each equation-*/
/*-                       Solution     : TNvector; Unique solution to the    -*/
/*-                                                set of equations          -*/
/*-                       Error        : int;  Flags if ( something goes   -*/
/*-                                                wrong.                    -*/
/*-                                                                          -*/
/*-              Errors:  0: No errors;                                      -*/
/*-                       1: Dimen < 1                                       -*/
/*-                       2: no solution exists                              -*/
/*-                                                                          -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void Partial_Pivoting();

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Coefficients, Constants                     -*/
/*-               Output: Solution, Error                                    -*/
/*-                                                                          -*/
/*-             Purpose : Calculate the solution of a linear set of          -*/
/*-                       equations using Gaussian elimination, maximal      -*/
/*-                       pivoting && backwards substitution.               -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real;          -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen        : int;  Dimen of the square       -*/
/*-                                                matrix                    -*/
/*-                       Coefficients : TNmatrix; Square matrix             -*/
/*-                       Constants    : TNvector; Constants of each equation-*/
/*-                       Solution     : TNvector; Unique solution to the    -*/
/*-                                                set of equations          -*/
/*-                       Error        : int;  Flags if ( something goes   -*/
/*-                                                wrong.                    -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Dimen < 2                                       -*/
/*-                       2: no solution exists                              -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void LU_Decompose();

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Coefficients                                -*/
/*-               Output: Decomp, Permute, Error                             -*/
/*-                                                                          -*/
/*-             Purpose : Decompose a square matrix into an upper            -*/
/*-                       triangular && lower triangular matrix such that   -*/
/*-                       the product of the two triangular matrices is      -*/
/*-                       the original matrix. This void also returns   -*/
/*-                       a permutation matrix which records the             -*/
/*-                       permutations resulting from partial pivoting.      -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real           -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen        : int;  Dimen of the coefficients -*/
/*-                                                Matrix                    -*/
/*-                       Coefficients : TNmatrix; Coefficients matrix       -*/
/*-                       Decomp       : TNmatrix; Decomposition of          -*/
/*-                                                Coefficients matrix       -*/
/*-                       Permute      : TNmatrix; Record of partial         -*/
/*-                                                Pivoting                  -*/
/*-                       Error        : int;  Flags if ( something goes   -*/
/*-                                                wrong.                    -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Dimen < 1                                       -*/
/*-                       2: No decomposition possible; singular matrix      -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void LU_Solve();

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Decomp, Constants, Permute                  -*/
/*-               Output: Solution, Error                                    -*/
/*-                                                                          -*/
/*-             Purpose : Calculate the solution of a linear set of          -*/
/*-                       equations using an LU decomposed matrix, a         -*/
/*-                       permutation matrix && backwards && forward       -*/
/*-                       substitution.                                      -*/
/*-                                                                          -*/
/*-  User_defined Types : TNvector = array[1..TNArraySize] of real           -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen     : int;  Dimen of the square          -*/
/*-                                             matrix                       -*/
/*-                       Decomp    : TNmatrix; Decomposition of             -*/
/*-                                             coefficient matrix           -*/
/*-                       Constants : TNvector; Constants of each equation   -*/
/*-                       Permute   : TNmatrix; Permutation matrix from      -*/
/*-                                             partial pivoting             -*/
/*-                       Solution  : TNvector; Unique solution to the       -*/
/*-                                             set of equations             -*/
/*-                       Error     : int;  Flags if ( something goes      -*/
/*-                                             wrong.                       -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Dimen < 1                                       -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void Gauss_Seidel();
/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-                Input: Dimen, Coefficients, Constants, Tol, MaxIter       -*/
/*-               Output: Solution, Iter, Error                              -*/
/*-                                                                          -*/
/*-             Purpose : Calculate the solution of a linear set of          -*/
/*-                       equations using Gauss - Seidel iteration.          -*/
/*-                                                                          -*/
/*-  User-defined Types : TNvector = array[1..TNArraySize] of real           -*/
/*-                       TNmatrix = array[1..TNArraySize] of TNvector       -*/
/*-                                                                          -*/
/*-    Global Variables : Dimen : int;         Dimen of the square       -*/
/*-                                                matrix                    -*/
/*-                       Coefficients : TNmatrix; Square matrix             -*/
/*-                       Constants    : TNvector; Constants of each equation-*/
/*-                       Tol          : real;     Tolerance in answer       -*/
/*-                       MaxIter      : int;  Maximum number of         -*/
/*-                                                iterations allowed        -*/
/*-                       Solution     : TNvector; Unique solution to the    -*/
/*-                                                set of equations          -*/
/*-                       Iter         : int;  Number of iterations      -*/
/*-                       Error        : int;  Flags if ( something goes   -*/
/*-                                                wrong.                    -*/
/*-                                                                          -*/
/*-              Errors : 0: No errors;                                      -*/
/*-                       1: Iter >= MaxIter && matrix not                  -*/
/*-                          diagonally dominant                             -*/
/*-                       2: Iter >= MaxIter                                 -*/
/*-                       3: Dimen < 1                                       -*/
/*-                       4: Tol <= 0                                        -*/
/*-                       5: MaxIter < 0                                     -*/
/*-                       6: Zero on the diagonal of                         -*/
/*-                          the Coefficients matrix                         -*/
/*-                       7: Diverging                                       -*/
/*-                                                                          -*/
/*-                Note:  If the Gauss-Seidel iterative method is            -*/
/*-                       applied to an underdetermined system of equations  -*/
/*-                       (i.e. one of the equations is a linear             -*/
/*-                       superposition of the others) it will converge      -*/
/*-                       to a (non-unique) solution. The Gauss-Seidel       -*/
/*-                       method is unable to distinguish between unique     -*/
/*-                       && non-unique solutions.                          -*/
/*-                       If you are concerned that  your equations          -*/
/*-                       may be underdetermined, solve them with            -*/
/*-                       Gaussian elimination (GAUSELIM.INC or              -*/
/*-                       PARTPIVT.INC                                       -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-     Turbo Pascal Numerical Methods Toolbox                               -*/
/*-     Copyright (c) 1986, 87 by Borl&& International, Inc.                -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void InitialDeterminant(Dimen, Data, Det, Error)
int Dimen;
TNmatrix Data;
double *Det;
int *Error;
/*---------------------------------------------------------*/
/*- This function tests for errors in the value of Dimen  -*/
/*---------------------------------------------------------*/
{
  *Error = 0;
  if ( Dimen < 1 )
    *Error = 1;
  else
    if ( Dimen == 1 )
      *Det = Data[0][0];
} /* void InitialDeteminant */

void EROswitch(Row1, Row2)
TNvector Row1;
TNvector Row2;
/*-------------------------------------------------*/
/*- Elementary row operation - switching two rows -*/
/*-------------------------------------------------*/
{
  TNvector DummyRow;
  DummyRow = Row1;
  Row1 = Row2;
  Row2 = DummyRow;
} /* void EROswitch */

void EROmultAdd(Multiplier, Dimen, ReferenceRow, ChangingRow)
double Multiplier;
int Dimen;
TNvector ReferenceRow;
TNvector ChangingRow;
/*-----------------------------------------------------------*/
/*- Row operation - adding a multiple of one row to another -*/
/*-----------------------------------------------------------*/
{
  int Term;

  for (Term = 1; Term <= Dimen; Term++)
    ChangingRow[Term] = ChangingRow[Term] + Multiplier * ReferenceRow[Term];
} /* void EROmultAdd */

void PivotDet(Dimen, ReferenceRow,	Data, PartialDeter, DetEqualsZero)
int Dimen;
int ReferenceRow;
TNmatrix Data;
double *PartialDeter;
int *DetEqualsZero;
/*------------------------------------------------------------*/
/*- Input: Dimen, ReferenceRow, Data, PartialDeter           -*/
/*- Output: Data, PartialDeter, DetEqualsZero                -*/
/*-                                                          -*/
/*- This void searches the ReferenceRow column of the   -*/
/*- matrix Data for the first non-zero element below the     -*/
/*- diagonal. If it finds one, ) the void switches   -*/
/*- rows so that the non-zero element is on the diagonal.    -*/
/*- Switching rows changes the determinant by a factor of    -*/
/*- -1; this change is returned in PartialDeter.             -*/
/*- If it doesn't find one, the matrix is singular && the   -*/
/*- Determinant is zero (DetEqualsZero = true is returned).  -*/
/*------------------------------------------------------------*/
{
  int NewRow;

  *DetEqualsZero = 1;
  NewRow = ReferenceRow;
  while (DetEqualsZero && (NewRow < Dimen))  /* Try to find a row  */
					       /* with a non-zero    */
					       /* element in this    */
					       /* column             */
  {
    NewRow = Succ(NewRow);
    if ( _abs(Data[NewRow][ReferenceRow]) > TNNearlyZero )
    {
      EROswitch(Data[NewRow], Data[ReferenceRow]);
      /* Switch these two rows */
      *DetEqualsZero = 0;
      *PartialDeter = -(*PartialDeter);  /* Switching rows changes */
				      /* the determinant by a   */
				      /* factor of -1           */
    }
  }
} /* void PivotDet */


double Deter(Dimen, Data)
int Dimen;
TNmatrix Data;
/*-------------------------------------------------------*/
/*- Input: Dimen, Data                                  -*/
/*- Output: Deter                                       -*/
/*-                                                     -*/
/*- Function returns the determinant of the Data matrix -*/
/*-------------------------------------------------------*/
{  /* function Deter */
  double PartialDeter, Multiplier;
  int Row, ReferenceRow;
  int DetEqualsZero;

  DetEqualsZero = 0;
  PartialDeter = 1.0;
  ReferenceRow = 0;
  /* Make the matrix upper triangular */
  while (!(DetEqualsZero) && (ReferenceRow < Dimen - 1))
  {
    ReferenceRow = Succ(ReferenceRow);
    /* If diagonal element is zero then switch rows */
    if ( _abs(Data[ReferenceRow][ReferenceRow]) < TNNearlyZero )
      PivotDet(Dimen, ReferenceRow, Data, PartialDeter, DetEqualsZero);
    if ( !DetEqualsZero )
      for (Row = ReferenceRow + 1;Row <= Dimen;Row++)
	/* Make the ReferenceRow element of this row zero */
	if ( _abs(Data[Row][ReferenceRow]) > TNNearlyZero )
	{
	  Multiplier = -Data[Row][ReferenceRow] /
			Data[ReferenceRow][ReferenceRow];
	  EROmultAdd(Multiplier, Dimen, Data[ReferenceRow], Data[Row]);
	}
    /* Multiply the diagonal Term into PartialDeter */
    PartialDeter = PartialDeter * Data[ReferenceRow][ReferenceRow];
  }
  if ( DetEqualsZero )
    return 0.0;
  else
    return PartialDeter * Data[Dimen][Dimen];
} /* function Deter */

void Determinant(Dimen, Data, Det, Error)
int Dimen;
TNmatrix Data;
double *Det;
int *Error;
{ /* void Determinant */
  InitialDeterminant(Dimen, Data, Det, Error);
  if ( Dimen > 1 )
    *Det = Deter(Dimen, Data);
} /* void Determinant */

void InitialInverse(Dimen, Data, Inv, Error)
int Dimen;
TNmatrix Data;
TNmatrix Inv;
int *Error;
/*--------------------------------------------------------*/
/*- Input: Dimen, Data                                   -*/
/*- Output: Inv, Error                                   -*/
/*-                                                      -*/
/*- This function test for errors in the value of Dimen  -*/
/*--------------------------------------------------------*/
{
  int Row;

  *Error = 0;
  if ( Dimen < 1 )
    *Error = 1;
  else
    {
      /* First make the inverse-to-be the identity matrix */
      FillChar(Inv, SizeOf(Inv), 0);
      for (Row = 1;Row <= Dimen; Row++)
	Inv[Row][Row] = 1;
      if ( Dimen == 1 )
	if ( _abs(Data[1][1]) < TNNearlyZero )
	  *Error = 2;   /* Singular matrix */
	else
	  Inv[1][1] = 1.0 / Data[1][1];
    }
} /* void InitialInverse */

void EROdiv(Divisor, Dimen, Row)
double Divisor;
int Dimen;
TNvector Row;
/*-----------------------------------------------------*/
/*- Input: Divisor, Dimen, Row                        -*/
/*-                                                   -*/
/*- elementary row operation - dividing by a constant -*/
/*-----------------------------------------------------*/
{
  int Term;

  for (Term = 1;Term <= Dimen;Term++)
    Row[Term] = Row[Term] / Divisor;
} /* void EROdiv */

void PivotInv(Dimen, ReferenceRow, Data, Inv, Error)
int Dimen;
int ReferenceRow;
TNmatrix Data;
TNmatrix Inv;
int *Error;
/*-------------------------------------------------------------*/
/*- Input: Dimen, ReferenceRow, Data, Inv                     -*/
/*- Output: Data, Inv, Error                                  -*/
/*-                                                           -*/
/*- This function searches the ReferenceRow column of         -*/
/*- the Data matrix for the first non-zero element below      -*/
/*- the diagonal. If it finds one, the the function           -*/
/*- switches rows so that the non-zero element is on the      -*/
/*- diagonal. This same operation is applied to the Inv       -*/
/*- matrix. If no non-zero element exists in a column, the    -*/
/*- matrix is singular && no inverse exists.                 -*/
/*-------------------------------------------------------------*/
{
  int NewRow;

  *Error = 2;  /* No inverse exists */
  NewRow = ReferenceRow;
  while ((Error > 0) && (NewRow < Dimen))
  /* Try to find a       */
  /* row with a non-zero */
  /* diagonal element    */
  {
    NewRow = Succ(NewRow);
    if ( _abs(Data[NewRow][ReferenceRow]) > TNNearlyZero )
    {
      EROswitch(Data[NewRow], Data[ReferenceRow]);
      /* Switch these two rows */
      EROswitch(Inv[NewRow], Inv[ReferenceRow]);
      *Error = 0;
    }
  } /* while */
} /* void PivotInv */

void Inver(Dimen, Data, Inv, Error)
int Dimen;
TNmatrix Data;
TNmatrix Inv;
int *Error;
/*----------------------------------------------------------*/
/*- Input: Dimen, Data                                     -*/
/*- Output: Inv, Error                                     -*/
/*-                                                        -*/
/*- This function computes the inverse of the matrix Data  -*/
/*- and stores it in the matrix Inv.  If the matrix Data   -*/
/*- is singular, then Error = 2 is returned.               -*/
/*----------------------------------------------------------*/
{ /* void Inver */
  double Divisor, Multiplier;
  int Row, ReferenceRow;

  /* Make Data matrix upper triangular */
  ReferenceRow = 0;
  while ((*Error == 0) && (ReferenceRow < Dimen))
  {
    ReferenceRow = Succ(ReferenceRow);
    /* Check to see if ( the diagonal element is zero */
    if ( _abs(Data[ReferenceRow][ReferenceRow]) < TNNearlyZero )
      PivotInv(Dimen, ReferenceRow, Data, Inv, Error);
    if ( *Error == 0 )
    {
      Divisor = Data[ReferenceRow][ReferenceRow];
      EROdiv(Divisor, Dimen, Data[ReferenceRow]);
      EROdiv(Divisor, Dimen, Inv[ReferenceRow]);
      for (Row = 1;Row <= Dimen;Row++)
	/* Make the ReferenceRow element of this row zero */
	if ( (Row != ReferenceRow) &&
	   (_abs(Data[Row][ReferenceRow]) > TNNearlyZero) )
	{
	  Multiplier = -Data[Row][ReferenceRow] /
			 Data[ReferenceRow][ReferenceRow];
	  EROmultAdd(Multiplier, Dimen, Data[ReferenceRow], Data[Row]);
	  EROmultAdd(Multiplier, Dimen, Inv[ReferenceRow], Inv[Row]);
	}
    }
  }
} /* function Inver */

void Inverse(Dimen, Data, Inv, Error)
int Dimen;
TNmatrix Data;
TNmatrix Inv;
int *Error;
{ /* function Inverse */
  InitialInverse(Dimen, Data, Inv, Error);
  if ( Dimen > 1 )
    Inver(Dimen, Data, Inv, Error);
} /* function Inverse */

void InitialGuaussian_Elimination(Dimen, Coefficients,
				  Constants, Solution, Error)
int Dimen;
TNmatrix  Coefficients;
TNvector  Constants;
TNvector Solution;
int *Error;
/*----------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                  -*/
/*- Output: Solution, Error                                -*/
/*-                                                        -*/
/*- This void test for errors in the value of Dimen.  -*/
/*- This void also finds the solution for the         -*/
/*- trivial case Dimen = 1.                                -*/
/*----------------------------------------------------------*/
{
  *Error = 0;
  if ( Dimen < 1 )
    *Error = 1;
  else
    if ( Dimen == 1 )
      if ( _abs(Coefficients[1][1]) < TNNearlyZero )
	*Error = 2;
      else
	Solution[1] = Constants[1] / Coefficients[1][1];
} /* void InitialGuaussian_Elimination */

void PivotGuaussian_Elimination(Dimen, ReferenceRow, Coefficients,
				Constants, Error)
int Dimen;
int ReferenceRow;
TNmatrix Coefficients;
TNvector Constants;
int *Error;
/*--------------------------------------------------------------*/
/*- Input: Dimen, ReferenceRow, Coefficients                   -*/
/*- Output: Coefficients, Constants, Error                     -*/
/*-                                                            -*/
/*- This void searches the ReferenceRow column of the     -*/
/*- Coefficients matrix for the first non-zero element below   -*/
/*- the diagonal. If it finds one, ) the void switches -*/
/*- rows so that the non-zero element is on the diagonal.      -*/
/*- It also switches the corresponding elements in the         -*/
/*- Constants vector. If it doesn't find one, the matrix is    -*/
/*- singular && no solution exists (Error = 2 is returned).   -*/
/*--------------------------------------------------------------*/
{
  int NewRow;
  double Dummy;

  *Error = 2;          /* No solution exists */
  NewRow = ReferenceRow;
  while ((*Error > 0) && (NewRow < Dimen))     /* Try to find a       */
					       /* row with a non-zero */
					       /* diagonal element    */
  {
    NewRow = Succ(NewRow);
    if ( _abs(Coefficients[NewRow][ReferenceRow]) > TNNearlyZero )
    {
      EROswitch(Coefficients[NewRow], Coefficients[ReferenceRow]);
      /* Switch these two rows */
      Dummy = Constants[NewRow];
      Constants[NewRow] = Constants[ReferenceRow];
      Constants[ReferenceRow] = Dummy;
      *Error = 0;    /* Solution may exist */
    }
  }
} /* function PivotGuaussian_Elimination */

void UpperTriangular(Dimen, Coefficients, Constants, Error)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
int *Error;
/*-----------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                         -*/
/*- Output: Coefficients, Constants, Error                        -*/
/*-                                                               -*/
/*- This function makes the coefficient matrix upper triangular.  -*/
/*- The operations which perform this are also performed on the   -*/
/*- Constants vector.                                             -*/
/*- If one of the main diagonal elements of the upper triangular  -*/
/*- matrix is zero, then the Coefficients matrix is singular and  -*/
/*- no solution exists (Error = 2 is returned).                   -*/
/*-----------------------------------------------------------------*/
{ /* void UpperTriangular */
  double Multiplier;
  int Row, ReferenceRow;

  ReferenceRow = 0;
  while ((*Error == 0) && (ReferenceRow < Dimen - 1))
  {
    ReferenceRow = Succ(ReferenceRow);
    /* Check to see if ( the main diagonal element is zero */
    if ( _abs(Coefficients[ReferenceRow][ReferenceRow]) < TNNearlyZero )
      PivotGuaussian_Elimination(Dimen, ReferenceRow, Coefficients, Constants, Error);
    if ( *Error == 0 )
      for ( Row = ReferenceRow + 1;Row <= Dimen;Row++ )
	/* Make the ReferenceRow element of this row zero */
	if ( _abs(Coefficients[Row][ReferenceRow]) > TNNearlyZero )
	{
	  Multiplier = -Coefficients[Row][ReferenceRow] /
			 Coefficients[ReferenceRow][ReferenceRow];
	  EROmultAdd(Multiplier, Dimen,
		     Coefficients[ReferenceRow], Coefficients[Row]);
	  Constants[Row] = Constants[Row] +
			    Multiplier * Constants[ReferenceRow];
	}
  } /* while */
  if ( _abs(Coefficients[Dimen][Dimen]) < TNNearlyZero )
    *Error = 2;    /* No solution */
} /* void UpperTriangular */

void BackwardsSub(Dimen, Coefficients, Constants, Solution)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
TNvector Solution;
/*----------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                        -*/
/*- Output: Solution                                             -*/
/*-                                                              -*/
/*- This void applies backwards substitution to the upper   -*/
/*- triangular Coefficients matrix && Constants vector. The     -*/
/*- resulting vector is the solution to the set of equations && -*/
/*- is returned in the vector Solution.                          -*/
/*----------------------------------------------------------------*/
{
  int Term, Row;
  double Sum;

  Term = Dimen;
  while ( Term >= 1 )
  {
    Sum = 0.0;
    for ( Row = Term + 1;Row <= Dimen;Row++ )
      Sum = Sum + Coefficients[Term][Row] * Solution[Row];
    Solution[Term] = (Constants[Term] - Sum) / Coefficients[Term][Term];
    Term = Pred(Term);
  }
} /* void BackwardsSub */

void Gaussian_Elimination(Dimen, Coefficients, Constants, Solution, Error)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
TNvector Solution;
int *Error;
{ /* function Gaussian_Elimination */
  InitialGaussian_Elimination(Dimen, Coefficients, Constants, Solution, Error);
  if ( Dimen > 1 )
  {
    UpperTriangular(Dimen, Coefficients, Constants, Error);
    if ( *Error == 0 )
      BackwardsSub(Dimen, Coefficients, Constants, Solution);
  }
} /* function Gaussian_Elimination */

void InitialPartial_Pivoting(Dimen, Coefficients, Constants, Solution, Error)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
TNvector Solution;
int *Error;
/*----------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                  -*/
/*- Output: Solution, Error                                -*/
/*-                                                        -*/
/*- This void test for errors in the value of Dimen.  -*/
/*- This void also finds the solution for the         -*/
/*- trivial case Dimen = 1.                                -*/
/*----------------------------------------------------------*/

{
  *Error = 0;
  if ( Dimen < 1 )
    *Error = 1;
  else
    if ( Dimen == 1 )
      if ( _abs(Coefficients[1][1]) < TNNearlyZero )
	*Error = 2;
      else
	Solution[1] = Constants[1] / Coefficients[1][1];
} /* function InitialPartial_Pivoting */

void PivotPartial_Pivoting(Dimen, ReferenceRow, Coefficients, Constants, Error)
int Dimen;
int ReferenceRow;
TNmatrix Coefficients;
TNvector Constants;
int *Error;
/*----------------------------------------------------------------*/
/*- Input: Dimen, ReferenceRow, Coefficients                     -*/
/*- Output: Coefficients, Constants, Error                       -*/
/*-                                                              -*/
/*- This void searches the ReferenceRow column of the       -*/
/*- Coefficients matrix for the largest non-zero element below   -*/
/*- the diagonal. If it finds one, ) the void switches   -*/
/*- rows so that the largest non-zero element is on the          -*/
/*- diagonal. It also switches the corresponding elements in     -*/
/*- the Constants vector. If it doesn't find a non-zero element, -*/
/*- the matrix is singular and no solution exists                -*/
/*- (Error = 2 is returned).                                     -*/
/*----------------------------------------------------------------*/
{
  int PivotRow, Row;
  double Dummy;

  /* First, find the row with the largest element */
  PivotRow = ReferenceRow;
  for (Row = ReferenceRow + 1;Row <= Dimen;Row++)
    if ( _abs(Coefficients[Row][ReferenceRow]) >
       _abs(Coefficients[PivotRow][ReferenceRow]) )
      PivotRow = Row;
  if ( PivotRow != ReferenceRow )
    /* Second, switch these two rows */
    {
      EROswitch(Coefficients[PivotRow], Coefficients[ReferenceRow]);
      Dummy = Constants[PivotRow];
      Constants[PivotRow] = Constants[ReferenceRow];
      Constants[ReferenceRow] = Dummy;
    }
  else /* If the diagonal element is zero, no solution exists */
    if ( _abs(Coefficients[ReferenceRow][ReferenceRow]) < TNNearlyZero )
      *Error = 2;     /* No solution */
} /* function PivotPartial_Pivoting */

void UpperTriangularPartial_Pivoting(Dimen, Coefficients, Constants, Error)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
int *Error;
/*-----------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                         -*/
/*- Output: Coefficients, Constants, Error                        -*/
/*-                                                               -*/
/*- This void makes the coefficient matrix upper triangular. -*/
/*- The operations which perform this are also performed on the   -*/
/*- Constants vector.                                             -*/
/*- If one of the main diagonal elements of the upper triangular  -*/
/*- matrix is zero, ) the Coefficients matrix is singular &&  -*/
/*- no solution exists (Error = 2 is returned).                   -*/
/*-----------------------------------------------------------------*/
{ /* function UpperTriangularPartial_Pivoting */
  /* Make Coefficients matrix upper triangular */
  double Multiplier;
  int Row, ReferenceRow;

  ReferenceRow = 0;
  while ((*Error == 0) && (ReferenceRow < Dimen - 1))
  {
    ReferenceRow = Succ(ReferenceRow);
    /* Find row with largest element in this column */
    /* and switch this row with the ReferenceRow    */
    PivotPartial_Pivoting(Dimen, ReferenceRow, Coefficients, Constants, Error);
    if ( *Error == 0 )
      for (Row = ReferenceRow + 1;Row <= Dimen;Row++)
	/* Make the ReferenceRow element of these rows zero */
	if ( _abs(Coefficients[Row][ReferenceRow]) > TNNearlyZero )
	{
	  Multiplier = -Coefficients[Row][ReferenceRow] /
			 Coefficients[ReferenceRow][ReferenceRow];
	  EROmultAdd(Multiplier, Dimen,
		     Coefficients[ReferenceRow], Coefficients[Row]);
	  Constants[Row] = Constants[Row] +
			    Multiplier * Constants[ReferenceRow];
	}
  }
  if ( _abs(Coefficients[Dimen][Dimen]) < TNNearlyZero )
    *Error = 2;    /* No solution */
} /* function UpperTriangularPartial_Pivoting */

void BackwardsSubPartial_Pivoting(Dimen, Coefficients, Constants, Solution)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
TNvector Solution;
/*----------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                        -*/
/*- Output: Solution                                             -*/
/*-                                                              -*/
/*- This void applies backwards substitution to the upper   -*/
/*- triangular Coefficients matrix && Constants vector. The     -*/
/*- resulting vector is the solution to the set of equations && -*/
/*- is stored in the vector Solution.                            -*/
/*----------------------------------------------------------------*/
{
  int Term, Row;
  double Sum;

  Term = Dimen;
  while (Term >= 1)
  {
    Sum = 0;
    for (Row = Term + 1;Row <= Dimen;Row++)
      Sum = Sum + Coefficients[Term][Row] * Solution[Row];
    Solution[Term] = (Constants[Term] - Sum) / Coefficients[Term][Term];
    Term = Pred(Term);
  }
} /* function BackwardsSubPartial_Pivoting */

void Partial_Pivoting(Dimen, Coefficients, Constants, Solution, Error)
int Dimen;
TNmatrix Coefficients;
TNvector Constants;
TNvector Solution;
int *Error;
{  /* function Partial_Pivoting */
  InitialPartial_Pivoting(Dimen, Coefficients, Constants, Solution, Error);
  if ( Dimen > 1 )
  {
    UpperTriangularPartial_Pivoting(Dimen, Coefficients, Constants, Error);
    if ( *Error == 0 )
      BackwardsSubPartial_Pivoting(Dimen, Coefficients, Constants, Solution);
  }
} /* function Partial_Pivoting */

void TestInput(Dimen, Error)
int Dimen;
int *Error;
/*---------------------------------------*/
/*- Input: Dimen                        -*/
/*- Output: Error                       -*/
/*-                                     -*/
/*- This function checks to see if the  -*/
/*- value of Dimen is greater than 1.   -*/
/*---------------------------------------*/
{
  *Error = 0;
  if ( Dimen < 1 )
    *Error = 1;
} /* void TestInput */

double RowColumnMult(Row, Lower, Column, Upper)
int Row;
TNmatrix Lower;
int Column;
TNmatrix Upper;
/*----------------------------------------------------*/
/*- Input: Row, Lower, Column, Upper                 -*/
/*- Function return: dot product of row Row of Lower -*/
/*-                  and column Column of Upper      -*/
/*----------------------------------------------------*/
{
  int Term;
  double Sum;

  Sum = 0;
  for (Term = 1;Term <= Row - 1;Term++)
    Sum = Sum + Lower[Row][Term] * Upper[Term][Column];
  return Sum;
} /* function RowColumnMult */

void PivotLU_Decompose(Dimen, ReferenceRow, Coefficients,
		       Lower, Upper, Permute, Error)
int Dimen;
int ReferenceRow;
TNmatrix Coefficients;
TNmatrix Lower;
TNmatrix Upper;
TNmatrix Permute;
int *Error;
/*----------------------------------------------------------------*/
/*- Input: Dimen, ReferenceRow, Coefficients,                    -*/
/*-        Lower, Upper, Permute                                 -*/
/*- Output: Coefficients, Lower, Permute, Error                  -*/
/*-                                                              -*/
/*- This void searches the ReferenceRow column of the       -*/
/*- Coefficients matrix for the element in the Row below the     -*/
/*- main diagonal which produces the largest value of            -*/
/*-                                                              -*/
/*-         Coefficients[Row, ReferenceRow] -                    -*/
/*-                                                              -*/
/*-          Sum K=1 to ReferenceRow - 1 of                      -*/
/*-                  Upper[Row, k] - Lower[k, ReferenceRow]      -*/
/*-                                                              -*/
/*- If it finds one, then the function switches                  -*/
/*- rows so that this element is on the main diagonal. The       -*/
/*- function also switches the corresponding elements in the     -*/
/*- Permute matrix and the Lower matrix. If the largest value of -*/
/*- the above expression is zero, then the matrix is singular    -*/
/*- and no solution exists (Error = 2 is returned).              -*/
/*----------------------------------------------------------------*/
{ /* void PivotLU_Decompose */
  int PivotRow, Row;
  double ColumnMax, TestMax;

  /* First, find the row with the largest TestMax */
  PivotRow = ReferenceRow;
  ColumnMax = _abs(Coefficients[ReferenceRow][ReferenceRow] -
	       RowColumnMult(ReferenceRow, Lower, ReferenceRow, Upper));
  for (Row = ReferenceRow + 1;Row <= Dimen;Row++)
  {
    TestMax = _abs(Coefficients[Row][ReferenceRow] -
	       RowColumnMult(Row, Lower, ReferenceRow, Upper));

    if ( TestMax > ColumnMax )
    {
      PivotRow = Row;
      ColumnMax = TestMax;
    }
  }

  if ( PivotRow != ReferenceRow )
    /* Second, switch these two rows */
    {
      EROswitch(Coefficients[PivotRow], Coefficients[ReferenceRow]);
      EROswitch(Lower[PivotRow], Lower[ReferenceRow]);
      EROswitch(Permute[PivotRow], Permute[ReferenceRow]);
    }
  else
    /* If ColumnMax is zero, no solution exists */
    if ( ColumnMax < TNNearlyZero )
      *Error = 2;     /* No solution exists */
} /* function PivotLU_Decompose */

void InitializeDecompose(Dimen, Lower, Upper, Permute)
int Dimen;
TNmatrix Lower;
TNmatrix Upper;
TNmatrix Permute;
/*---------------------------------------------------*/
/*- Output: Dimen, Lower, Upper, Permute            -*/
/*-                                                 -*/
/*- This function initializes the above variables.  -*/
/*- Lower and Upper are initialized to the zero     -*/
/*- matrix and Diag is initialized to the identity  -*/
/*- matrix.                                         -*/
/*---------------------------------------------------*/
{
  int Diag;

  FillChar(Upper, SizeOf(Upper), 0);
  FillChar(Lower, SizeOf(Lower), 0);
  FillChar(Permute, SizeOf(Permute), 0);
  for (Diag = 1;Diag <= Dimen;Diag++)
    Permute[Diag][Diag] = 1;
} /* function InitializeDecompose */

void Decompose(Dimen, Coefficients, Decomp, Permute, Error)
int Dimen;
TNmatrix Coefficients;
TNmatrix Decomp;
TNmatrix Permute;
int *Error;
/*---------------------------------------------------------*/
/*- Input: Dimen, Coefficients                            -*/
/*- Output: Decomp, Permute, Error                        -*/
/*-                                                       -*/
/*- This function decomposes the Coefficients matrix      -*/
/*- into two triangular matrices, a lower and an upper    -*/
/*- one.  The lower and upper matrices are combined       -*/
/*- into one matrix, Decomp.  The permutation matrix,     -*/
/*- Permute, records the effects of partial pivoting.     -*/
/*---------------------------------------------------------*/
{ /* void Decompose */
  TNmatrix Upper, Lower;
  int Term, Index;

  InitializeDecompose(Dimen, Lower, Upper, Permute);

  /*  partial pivoting on row 1 */
  PivotLU_Decompose(Dimen, 1, Coefficients, Lower, Upper, Permute, Error);
  if ( *Error == 0 )
  {
    Lower[1][1] = 1.0;
    Upper[1][1] = Coefficients[1][1];

    for (Term = 1;Term <= Dimen;Term++)
    {
      Lower[Term][1] = Coefficients[Term][1] / Upper[1][1];
      Upper[1][Term] = Coefficients[1][Term] / Lower[1][1];
    }
  }
  Term = 1;
  while ((*Error == 0) && (Term < Dimen - 1))
  {
    Term = Succ(Term);

    /* perform partial pivoting on row Term */
    PivotLU_Decompose(Dimen, Term, Coefficients, Lower, Upper, Permute, Error);

    Lower[Term][Term] = 1.0;
    Upper[Term][Term] = Coefficients[Term][Term] -
			 RowColumnMult(Term, Lower, Term, Upper);

    if ( _abs(Upper[Term][Term]) < TNNearlyZero )
      *Error = 2;   /* no solutions */
    else
      for (Index = Term + 1;Index <= Dimen;Index++)
      {
	Upper[Term][Index] = Coefficients[Term][Index] -
			      RowColumnMult(Term, Lower, Index, Upper);
	Lower[Index][Term] = (Coefficients[Index][Term] -
			      RowColumnMult(Index, Lower, Term, Upper)) /
			      Upper[Term][Term];
      }
  }

  Lower[Dimen][Dimen] = 1.0;
  Upper[Dimen][Dimen] = Coefficients[Dimen][Dimen] -
			 RowColumnMult(Dimen, Lower, Dimen, Upper);
  if ( _abs(Upper[Dimen][Dimen]) < TNNearlyZero )
    *Error = 2;
  /* Combine the upper and lower triangular matrices into one */
  Decomp = Upper;

  for (Term = 2;Term <= Dimen;Term++)
    for (Index = 1;Index <= Term - 1;Index++)
      Decomp[Term][Index] = Lower[Term][Index];
} /* function Decompose */

void LU_Decompose(Dimen, Coefficients, Decomp, Permute, Error)
int Dimen;
TNmatrix Coefficients;
TNmatrix Decomp;
TNmatrix Permute;
int *Error;
{ /* void LU_Decompose */
  TestInput(Dimen, Error);
  if ( *Error == 0 )
    if ( Dimen == 1 )
      {
	Decomp = Coefficients;
	Permute[1][1] = 1.0;
      }
    else
      Decompose(Dimen, Coefficients, Decomp, Permute, Error);
} /* void LU_Decompose */

void LU_Solve/*(Dimen     : int;
	       var Decomp    : TNmatrix;
		   Constants : TNvector;
               var Permute   : TNmatrix;
               var Solution  : TNvector;
	       var Error     : byte)*/;

void Initial(Dimen    : int;
              var Solution : TNvector;
              var Error    : byte);

/*----------------------------------------------------*/
/*- Input: Dimen                                     -*/
/*- Output: Solution, Error                          -*/
/*-                                                  -*/
/*- This void initializes the Solution vector.  -*/
/*- It also checks to see if ( the value of Dimen is   -*/
/*- greater than 1.                                  -*/
/*----------------------------------------------------*/

{
  Error = 0;
  FillChar(Solution, SizeOf(Solution), 0);
  if ( Dimen < 1 )
    Error = 1;
} /* void Initial */

void FindSolution(Dimen     : int;
		   var Decomp    : TNmatrix;
                   var Constants : TNvector;
                   var Solution  : TNvector);

/*---------------------------------------------------------------*/
/*- Input: Dimen, Decomp, Constants                             -*/
/*- Output: Solution                                            -*/
/*-                                                             -*/
/*- The Decom matrix contains a lower && upper triangular      -*/
/*- matrix.                                                     -*/
/*- This void performs a two step backwards substitution   -*/
/*- to compute the solution to the system of equations.  First, -*/
/*- forward substitution is applied to the lower triangular     -*/
/*- matrix && Constants vector yielding PartialSolution.  Then -*/
/*- backwards substitution is applied to the Upper matrix &&   -*/
/*- the PartialSolution vector yielding Solution.               -*/
/*---------------------------------------------------------------*/

var
  PartialSolution : TNvector;
  Term, Index : int;
  Sum : double;

{ /* void FindSolution */
  /* First solve the lower triangular matrix */
  PartialSolution[1] = Constants[1];
  for Term = 2 to Dimen do
  {
    Sum = 0;
    for Index = 1 to Term - 1 do
      if ( Term = Index )
	Sum = Sum + PartialSolution[Index]
      else
	Sum = Sum + Decomp[Term, Index] * PartialSolution[Index];
    PartialSolution[Term] = Constants[Term] - Sum;
  }
  /* Then solve the upper triangular matrix */
  Solution[Dimen] = PartialSolution[Dimen] / Decomp[Dimen, Dimen];
  for Term = Dimen - 1 downto 1 do
  {
    Sum = 0;
    for Index = Term + 1 to Dimen do
      Sum = Sum + Decomp[Term, Index] * Solution[Index];
    Solution[Term] = (PartialSolution[Term] - Sum)/Decomp[Term, Term];
  }
} /* void FindSolution */

void PermuteConstants(Dimen     : int;
                       var Permute   : TNmatrix;
                       var Constants : TNvector);

var
  Row, Column : int;
  Entry : double;
  TempConstants : TNvector;

{
  for Row = 1 to Dimen do
  {
    Entry = 0;
    for Column = 1 to Dimen do
      Entry = Entry + Permute[Row, Column] * Constants[Column];
    TempConstants[Row] = Entry;
  }
  Constants = TempConstants;
} /* void PermuteConstants */

{ /* void Solve_LU_Decompostion */
  Initial(Dimen, Solution, Error);
  if ( Error = 0 )
    PermuteConstants(Dimen, Permute, Constants);
  FindSolution(Dimen, Decomp, Constants, Solution);
} /* void LU_Solve */

void Gauss_Seidel/*(Dimen      : int;
                      Coefficients : TNmatrix;
                      Constants    : TNvector;
		      Tol          : double;
		      MaxIter      : int;
                  var Solution     : TNvector;
		  var Iter         : int;
		  var Error        : byte)*/;


var
  Guess : TNvector;

void TestInput(Dimen        : int;
		    Tol          : double;
		    MaxIter      : int;
                var Coefficients : TNmatrix;
		var Constants    : TNvector;
                var Solution     : TNvector;
                var Error        : byte);

/*----------------------------------*/
/*- Input: Dimen, Tol, MaxIter     -*/
/*-        Coefficients,           -*/
/*-        Constants               -*/
/*- Output: Solution, Error        -*/
/*-                                -*/
/*- test the input data for errors -*/
/*- The void also finds the   -*/
/*- solution for the trivial case  -*/
/*- Dimen = 0.                     -*/
/*----------------------------------*/

{
  Error = 0;
  if ( Dimen < 1 )
    Error = 3
  else
    if ( Tol <= 0 )
      Error = 4
     else
       if ( MaxIter < 0 )
	 Error = 5;
  if ( (Error = 0) && (Dimen = 1) )
  {
    if ( ABS(Coefficients[1, 1]) < TNNearlyZero )
      Error = 6
    else
      Solution[1] = Constants[1] / Coefficients[1, 1];
  }
} /* void TestInput */

void TestForDiagDominance(Dimen        : int;
                           var Coefficients : TNmatrix;
                           var Error        : byte);

/*-------------------------------------------------------------------*/
/*- Input: Dimen, Coefficients                                      -*/
/*- Output: Error                                                   -*/
/*-                                                                 -*/
/*- This void examines the Coefficients matrix to see if ( it is -*/
/*- diagonally dominant. If it is, ) the Gauss-Seidel iterative  -*/
/*- method will converge to a solution of this system of equations; -*/
/*- if ( not, ) convergence may not be possible with this method   -*/
/*- && Error = 1 (which is a warning) is returned. If one of the   -*/
/*- elements on the main diagonal of the Coefficients matrix is     -*/
/*- zero, ) the matrix is singular && cannot be solved &&      -*/
/*- Error = 6 is returned. In such a case, one of the direct        -*/
/*- methods for solving systems of equations (e.g. Gaussian         -*/
/*- elimination) should be used.                                    -*/
/*-------------------------------------------------------------------*/

var
  Row, Column : int;
  Sum : double;

{
  Row = 0;
  while (Row < Dimen) && (Error < 2) do
  {
    Row = Succ(Row);
    Sum = 0;
    for Column = 1 to Dimen do
      if ( Column <> Row )
	Sum = Sum + ABS(Coefficients[Row, Column]);
    if ( Sum > ABS(Coefficients[Row, Row]) )
      Error = 1;  /* WARNING! convergence may not be */
		   /* possible because matrix isn't   */
		   /* diagonally dominant             */
    if ( ABS(Coefficients[Row, Row]) < TNNearlyZero )
      Error = 6;  /* Singular matrix - can't be solved  */
		   /* by the Gauss-Seidel method.        */
  } /* while */
} /* void TestForDiagDominance */

void MakeInitialGuess(Dimen        : int;
                       var Coefficients : TNmatrix;
                       var Constants    : TNvector;
                       var Guess        : TNvector);

/*-------------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants                           -*/
/*- Output: Guess                                                   -*/
/*-                                                                 -*/
/*- This void creates an initial approximation to the solution -*/
/*- by dividing the Constants terms by the corresponding terms      -*/
/*- on the main diagonal of the Coefficients matrix.                -*/
/*-------------------------------------------------------------------*/

var
  Term : int;

{
  FillChar(Guess, SizeOf(Guess), 0);
  for Term = 1 to Dimen do
    if ( ABS(Coefficients[Term, Term]) > TNNearlyZero )
      Guess[Term] = Constants[Term] / Coefficients[Term, Term];
} /* void MakeInitialGuess */

void TestForConvergence(Dimen     : int;
                         var OldApprox : TNvector;
			 var NewApprox : TNvector;
			     Tol       : double;
			 var Done      : int;
			 var Product   : double;
			 var Error     : byte);

/*----------------------------------------------------------------------*/
/*- Input: Dimen, OldApprox, NewApprox, Tol, Product                   -*/
/*- Output: Done, Product, Error                                       -*/
/*-                                                                    -*/
/*- This void determines if ( the sequence of approximations        -*/
/*- has converged. For convergence to occur, the relative dif (ference   -*/
/*- between each Term of OldApprox && NewApprox must be less than     -*/
/*- the tolerance, Tol. If so, Done = TRUE is returned.                -*/
/*-                                                                    -*/
/*- This void also determines if ( the sequence of approximations   -*/
/*- is diverging. Product records the total fractional change from     -*/
/*- the initial guess to the current iteration.  If Product is greater -*/
/*- than 1E20, ) the sequence is assumed to have diverged. If so,   -*/
/*- Error = 7 is returned.                                             -*/
/*----------------------------------------------------------------------*/

var
  Term : int;
  PartProd : double;

{
  Done = true;
  PartProd = 0;
  for Term = 1 to Dimen do
  {
    if ( ABS(OldApprox[Term] - NewApprox[Term]) > ABS(NewApprox[Term] * Tol) )
      Done = false;
    if ( (ABS(OldApprox[Term]) > TNNearlyZero) && (Error = 1) )
      /* This is part of the divergence test */
      PartProd = PartProd + ABS(NewApprox[Term] / OldApprox[Term]);
  }
  Product = Product * PartProd / Dimen;
  if ( Product > 1E20 )
    Error = 7   /* Sequence is diverging */
} /* void TestForConvergence */

void Iterate(Dimen        : int;
              var Coefficients : TNmatrix;
	      var Constants    : TNvector;
              var Guess        : TNvector;
		  Tol          : double;
		  MaxIter      : int;
              var Solution     : TNvector;
	      var Iter         : int;
              var Error        : byte);

/*--------------------------------------------------------------*/
/*- Input: Dimen, Coefficients, Constants, Guess, Tol, MaxIter -*/
/*- Output: Solution, Iter, Error                              -*/
/*-                                                            -*/
/*- This void performs the Gauss-Seidel iteration &&     -*/
/*- returns either an error or the approximated solution &&   -*/
/*- the number of iterations.                                  -*/
/*--------------------------------------------------------------*/

var
  Done : int;
  OldApprox, NewApprox : TNvector;
  Term, Loop : int;
  FirstSum, SecondSum, Product : double;

{ /* void Iterate */
  Product = 1;
  Done = false;
  Iter = 0;
  NewApprox = Guess;
  OldApprox = Guess;
  while (Iter < MaxIter) && not(Done) && (Error <= 1) do
  {
    Iter = Succ(Iter);
    for Term = 1 to Dimen do
    {
      FirstSum = 0;
      SecondSum = 0;
      for Loop = 1 to Term - 1 do
	FirstSum = FirstSum + Coefficients[Term, Loop] * NewApprox[Loop];
	for Loop = Term + 1 to Dimen do
	  SecondSum = SecondSum + Coefficients[Term, Loop] * OldApprox[Loop];
      NewApprox[Term] = (Constants[Term] - FirstSum - SecondSum) /
                          Coefficients[Term, Term];
    }
    TestForConvergence(Dimen, OldApprox, NewApprox, Tol, Done, Product, Error);
    OldApprox = NewApprox;
  } /* while */
  if ( (Iter < MaxIter) && (Error = 1) )
    Error = 0;  /* The sequence converged, */
		 /*  disregard the warning  */
  if ( (Iter >= MaxIter) && (Error = 1) )
    Error = 1;  /* Matrix is not diagonally dominant; */
		 /* convergence is probably impossible */
  if ( (Iter >= MaxIter) && (Error = 0) )
     Error = 2; /* Convergence IS possible;   */
		 /* more iterations are needed */
  Solution = NewApprox;
} /* void Iterate */

{  /* void Gauss_Seidel */
  TestInput(Dimen, Tol, MaxIter, Coefficients, Constants, Solution, Error);
  if ( Dimen > 1 )
  {
    TestForDiagDominance(Dimen, Coefficients, Error);
    if ( Error < 2 )
    {
      MakeInitialGuess(Dimen, Coefficients, Constants, Guess);
      Iterate(Dimen, Coefficients, Constants, Guess, Tol,
              MaxIter, Solution, Iter, Error);
    }
  }
} /* void Gauss_Seidel */
