/*
 * Matrix Inversion (Gauss-Jordan Elimination)
 *
 * Stephan Herhut (s.a.herhut@herts.ac.uk)
 * Frank Penczek  (f.penczek@herts.ac.uk)
 *
 * $Id$
 */


use Structures: all;
use StdIO: all;
use STAP_Mat_Invert: all;


complex[*] machtoc( int[*] data) 
{
  result = with {
             (. <= iv <= .) : toc( data[iv]);
           } : genarray( shape( data), toc(0));

  return( result);
}

void PrintImag( complex[*] B) 
{
  res = with {
          (. <= iv <= .) : imag( B[iv]);
        } : genarray( shape(B), 0.0);

  print(res);

  return();
}

void PrintReal( complex[*] B) 
{
  res = with {
          (. <= iv <= .) : real( B[iv]);
        } : genarray( shape(B), 0.0);

  print(res);

  return();
}

void PrintCplx( complex[*] B)
{
  PrintReal( B);
  PrintImag( B);
}

int main()
{

  A = machtoc( reshape( [2,2], [1,2,3,4]));

  printf( "A\n");
  PrintCplx(A);
  
  B = STAP_Mat_Invert( A);

  printf( "B\n");
  PrintCplx( B);

  return( 0);
}