/*
**  GE5
**
**  Purpose:
**    Gauss elimination, find upper triangular matrix
**
**  Date:
**    6/8/09
**
**  Version:
**    0     Start program, prepare the extended matrix
**    1     Eliminate one row
**    2     Introduce a loop, over rows
**    3     Introduce a function
**    4     Do multiple columns
**    5     Use function for columns
**
**  Author:
**    Charles Bos
*/
#include <oxstd.h>

/*
**  EliminateRow(const amAB, const j, const k)
**
**  Purpose:
**    Eliminate element j of column k
**
**  Inputs:
**    amAB      address of iK x iK+1 extended matrix
**    j         integer, row number
**    k         integer, pivot/column number
**
**  Outputs:
**    amAB      address of iK x iK+1 extended matrix, with zero
**              at element (j,k)
**
**  Return value:
**    TRUE if pivot not zero
*/
EliminateRow(const amAB, const j, const k)
{
  decl dF, ir;
  
  ir= amAB[0][k][k] == 0;
  
  dF= amAB[0][j][k] / amAB[0][k][k];

  println ("Multiplicity is ", dF);

  amAB[0][j][]= amAB[0][j][] - dF * amAB[0][k][]; 
  
  return !ir;
}

/*
**  EliminateColumn(const amAB, const k)
**
**  Purpose:
**    Eliminate all elements in column k below pivot
**
**  Inputs:
**    amAB      address of iK x iK+1 extended matrix
**    k         integer, pivot/column number
**
**  Outputs:
**    amAB      address of iK x iK+1 extended matrix, with zeros
**              at elements (k+1,k), ..., (iK-1,k)
**
**  Return value:
**    TRUE if ???
*/
EliminateColumn(const amAB, const k)
{
  decl j, iK, ir;
  
  iK= rows(amAB[0]);
  for (j= k+1; j < iK; ++j)
    {
      // Eliminate element j,k
      ir= EliminateRow(amAB, j, k);
    }
  print ("Extended matrix mAB after elimination step ", k, 
         " is: ", amAB[0]);
         
  return ir;         
}

/*
**  EliminateMatrix(const amAB)
**
**  Purpose:
**    Eliminate a matirx
**
**  Inputs:
**    amAB      address of iK x iK+iB matrix
**
**  Outputs:
**    amAB      address of matrix, triangularized
**
**  Return value:
**    None
*/
EliminateMatrix(const amAB)
{
  decl iK, k;
  
  // Get the number of rows
  iK= rows(amAB[0]);

  // Indicate the pivot number
  for (k= 0; k < iK-1; ++k)
    EliminateColumn(amAB, k);
}


main()
{
  decl mA, vB, mAB;
  
  // Magic numbers
  mA= < 6,  -2, 2,   4; 
       12,  -8, 6,  10; 
        3, -13, 9,   3; 
       -6,   4, 1, -18>; 
  vB= <16; 26; -19; -34>;
    
  // Initialisation
  mAB= mA~vB;
  
  print ("Extended matrix mAB is: ", mAB);
  
  EliminateMatrix(&mAB);
}