/*
**  GE1
**
**  Purpose:
**    Gauss elimination, find upper triangular matrix
**
**  Date:
**    20/8/2012
**
**  Version:
**    0     Start program, prepare the extended matrix
**    1     Start eliminating row 1
**
**  Author:
**    Charles Bos
*/
#include <oxstd.h>

/*
**  EliminateRow(const amAB, const j, const k)
**
**  Purpose:
**    Eliminate one row from the extended matrix mAB
**
**  Inputs:
**    amAB      address of the iK x iK+1 extended matrix
**    j         integer, row to eliminate
**    k         integer, pivot number, base row
**
**  Output:
**    amAB      address of iK x iK+1 extended matrix, with element [j][k] put to
**              zero
**
**  Return value:
**    none       
*/
EliminateRow(const amAB, const j, const k)
{
  decl i, iK, dMult;
  
  iK= rows(amAB[0]);
  dMult= amAB[0][j][k]/amAB[0][k][k];
  for (i= k; i < iK+1; ++i)
    amAB[0][j][i]= amAB[0][j][i] - dMult * amAB[0][k][i];
}

/*
**  EliminateColumn(const amAB, const k)
**
**  Purpose:
**    Eliminate column k of the extended matrix
**
**  Inputs:
**    amAB      address of the iK x iK+1 extended matrix
**    k         integer, pivot number, base row
**
**  Output:
**    amAB      address of iK x iK+1 extended matrix, with elements
**              below [k][k] all put to zero
**
**  Return value:
**    none
*/
EliminateColumn(const amAB, const k)
{
  decl iK, j;
  
  iK= rows(amAB[0]);
  for (j= k+1; j < iK; ++j)
    EliminateRow(amAB, j, k);
}

/*
**  EliminateMatrix(const amAB)
**
**  Purpose:
**    Eliminate extended matrix
**
**  Inputs:
**    amAB      address of the iK x iK+1 extended matrix
**
**  Output:
**    amAB      address of iK x iK+1 extended matrix, eliminated
**
**  Return value:
**    none
*/
EliminateMatrix(const amAB)
{
  decl iK, k;
  
  iK= rows(amAB[0]);
  for (k= 0; k < iK; ++k)
    EliminateColumn(amAB, k);
}




main()
{
  decl mA, vB, mAB, j, k;
  
  // 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;
  
  k= 0;
  j= 1;
  // Check
  print ("Extended matrix mAB is, before: ", mAB);
  EliminateMatrix(&mAB);
  
  // Check
  print ("Extended matrix mAB is, after: ", mAB);
}