/*********************************************************
 Note: This example runs only on Linux and Solaris platforms. 
 This example shows how to invert a non sparse matrix of 500x500
  elements of type complex with cartesian co-oprdinates with fields of 
  type real
 The matrix is initialized with the values of the Fourier transformation 
 and its inverse is compared to the known result. The distance between the inverse and the known result is displayed.
***********************************************************/
module top;
`define mod(val, n)  (val - (val/n)*n)

task ComputeDistance;
  begin
    distance = 0;
    for (i = 0; i < SIZE; i = i + 1)
    for (j = 0; j < SIZE; j = j + 1)
      begin
        WFC[i][j].Re = WFC[i][j].Re * SIZE;
        WFC[i][j].Im = WFC[i][j].Im * SIZE;
        polar = WFC[i][j]- WFCI[i][j];
        distance = (polar.Mag > distance) ? polar.Mag : distance;
      end
  end
endtask

parameter SIZE = 500;

VpFComplex WFC[SIZE-1:0][SIZE-1:0];
VpFComplex WFCI[SIZE-1:0][SIZE-1:0];
VpFPolar polar;

real phi0, distance;
integer i, j, mone;

initial begin
  phi0 = 2*$VpGetPi()/SIZE;
  $display("phi0 = %e\n", phi0);
  $InitM(WFC, $VpCos(phi0*`mod($I1*$I2, SIZE)),
              $VpSin(phi0*`mod($I1*$I2, SIZE)));
/*  $PrintM(WFC, "%e");*/
  $InitM(WFCI, $VpCos(phi0*`mod($I1*$I2, SIZE)),
              -$VpSin(phi0*`mod($I1*$I2, SIZE)));
/*  $PrintM(WFCI, "%e");*/

  mone = -1.0;
  WFC = WFC**(mone);
/*  $PrintM(WFC, "%e");*/
  ComputeDistance;
  $display("distance = %e\n", distance);
end
endmodule