with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Natural_Vectors;
with Standard_Floating_Vectors;
with Standard_Complex_Vectors;
with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;
with Standard_Complex_QR_Decomposition;  use Standard_Complex_QR_Decomposition;

package body Standard_Random_Matrices is

  function Random_Matrix ( n,m : natural; low,upp : integer )
                         return Standard_Integer_Matrices.Matrix is

    res : Standard_Integer_Matrices.Matrix(1..n,1..m);

  begin
    for i in 1..n loop
      for j in 1..m loop
        res(i,j) := Random(low,upp);
      end loop;
    end loop;
    return res;
  end Random_Matrix;

  function Random_Matrix ( n,m : natural )
                         return Standard_Floating_Matrices.Matrix is

    res : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    for i in 1..n loop
      for j in 1..m loop
        res(i,j) := Random;
      end loop;
    end loop;
    return res;
  end Random_Matrix;

  function Orthogonalize ( mat : Standard_Floating_Matrices.Matrix )
                         return Standard_Floating_Matrices.Matrix is

    n : constant natural := mat'length(1);
    m : constant natural := mat'length(2);
    res : Standard_Floating_Matrices.Matrix(1..n,1..m);
    wrk : Standard_Floating_Matrices.Matrix(1..n,1..m);
    bas : Standard_Floating_Matrices.Matrix(1..n,1..n);
    qra : Standard_Floating_Vectors.Vector(1..n) := (1..n => 0.0);
    pvt : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0); 

  begin
    wrk := mat;
    QRD(wrk,qra,pvt,false);
    for i in wrk'range(1) loop
      for j in wrk'range(2) loop
        bas(i,j) := wrk(i,j);
      end loop;
      for j in m+1..n loop
        bas(i,j) := 0.0;
	  end loop;
    end loop;
    Basis(bas,mat); 
    for i in res'range(1) loop
      for j in res'range(2) loop
        res(i,j) := bas(i,j);
      end loop;
    end loop;
    return res;
  end Orthogonalize;

  function Random_Orthogonal_Matrix
             ( n,m : natural ) return Standard_Floating_Matrices.Matrix is

    res : Standard_Floating_Matrices.Matrix(1..n,1..m)
        := Orthogonalize(Random_Matrix(n,m));

  begin
    return res;
  end Random_Orthogonal_Matrix;

  function Random_Matrix ( n,m : natural )
                         return Standard_Complex_Matrices.Matrix is

    res : Standard_Complex_Matrices.Matrix(1..n,1..m);

  begin
    for i in 1..n loop
      for j in 1..m loop
        res(i,j) := Random;
      end loop;
    end loop;
    return res;
  end Random_Matrix;

  function Orthogonalize ( mat : Standard_Complex_Matrices.Matrix )
                         return Standard_Complex_Matrices.Matrix is

    n : constant natural := mat'length(1);
    m : constant natural := mat'length(2);
    res : Standard_Complex_Matrices.Matrix(1..n,1..m);
    wrk : Standard_Complex_Matrices.Matrix(1..n,1..m);
    bas : Standard_Complex_Matrices.Matrix(1..n,1..n);
    qra : Standard_Complex_Vectors.Vector(1..n) := (1..n => Create(0.0));
    pvt : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0); 

  begin
    wrk := mat;
    QRD(wrk,qra,pvt,false);
    for i in wrk'range(1) loop
      for j in wrk'range(2) loop
        bas(i,j) := wrk(i,j);
      end loop;
      for j in m+1..n loop
        bas(i,j) := Create(0.0);
	  end loop;
    end loop;
    Basis(bas,mat); 
    for i in res'range(1) loop
      for j in res'range(2) loop
        res(i,j) := bas(i,j);
      end loop;
    end loop;
    return res;
  end Orthogonalize;

  function Random_Orthogonal_Matrix
                ( n,m : natural ) return Standard_Complex_Matrices.Matrix is

    res : Standard_Complex_Matrices.Matrix(1..n,1..m)
        := Orthogonalize(Random_Matrix(n,m));

  begin
    return res;
  end Random_Orthogonal_Matrix;

end Standard_Random_Matrices;