with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;

package body Standard_Complex_Least_Squares is

-- AUXILIARIES :

  function min0 ( a,b : integer ) return integer is

  -- DESCRIPTION : returns the minimum of a and b.

  begin
    if a <= b
     then return a;
     else return b;
	end if;
  end min0;                         

  procedure zcopy ( n,start : in natural;
                    x : in Standard_Complex_Vectors.Vector;
                    y : out Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Copies the n entries of x to y, beginning at start.

  begin
    for i in start..start+n-1 loop
      y(i) := x(i);
    end loop;
  end zcopy;

  function zdotc ( row : natural;
                   x : Standard_Complex_Matrices.Matrix;
                   y : Standard_Complex_Vectors.Vector )
                 return Complex_Number is

  -- DESCRIPTION :
  --   Dot product of two vectors : x(row..x'last(1),row)*y(row..y'last).

    res : Complex_Number := Create(0.0);

  begin
    for i in row..y'last loop
      res := res + Conjugate(x(i,row))*y(i);
    end loop;
    return res;
  end zdotc;

  procedure zaxpy ( n,row,col : in natural; f : in Complex_Number;
                    x : in Standard_Complex_Matrices.Matrix;
                    y : in out Standard_Complex_Vectors.Vector ) is

  -- DESCRIPTION :
  --   y(i) := y(i) + f*x(i,col) for i in row..row+n-1.

  begin
    for i in row..row+n-1 loop
      y(i) := y(i) + f*x(i,col);
    end loop;
  end zaxpy;

  procedure QRLS ( x : in out Standard_Complex_Matrices.Matrix;
                   ldx,n,k : in natural;
                   qraux,y : in Standard_Complex_Vectors.Vector;
                   qy,qty,b,rsd,xb : out Standard_Complex_Vectors.Vector;
                   job : in natural; info : out natural ) is

    cb,cqy,cqty,cr,cxb : boolean;
    jj,ju,kp1 : integer;
    t,temp : Complex_Number;

  begin
    info := 0;                                               -- set info flag
    cqy := (job/10000 /= 0);                     -- determine what to compute
    cqty := (job mod 10000 /= 0);
    cb := ((job mod 1000)/100 /= 0);
    cr := ((job mod 100)/10 /= 0);
    cxb := ((job mod 10) /= 0);
    ju := min0(k,n-1);
    if ju = 0                                    -- special action when n = 1
     then if cqy then qy(1) := y(1); end if;
          if cqty then qty(1) := y(1); end if;
          if cxb then xb(1) := y(1); end if;
          if cb
           then if AbsVal(x(1,1)) = 0.0
                 then info := 1;
                 else b(1) := y(1)/x(1,1);
                end if;
          end if;
          if cr then rsd(1) := Create(0.0); end if;
          return;
    end if;
    if cqy                                     -- set up to compute qy or qty
     then zcopy(n,y'first,y,qy);
    end if;
    if cqty
     then zcopy(n,y'first,y,qty);
    end if;
    if cqy                                                      -- compute qy
     then for j in 1..ju loop
            jj := ju - j + 1;
            if AbsVal(qraux(jj)) /= 0.0
             then temp := x(jj,jj);
                  x(jj,jj) := qraux(jj);
                  t := -zdotc(jj,x,qy)/x(jj,jj);
                  zaxpy(n-jj+1,jj,jj,t,x,qy);
                  x(jj,jj) := temp;
            end if;
          end loop;
    end if;
    if cqty                                             -- compute trans(q)*y
     then for j in 1..ju loop
            if AbsVal(qraux(j)) /= 0.0
             then temp := x(j,j);
                  x(j,j) := qraux(j);
                  t := -zdotc(j,x,qty)/x(j,j);
                  zaxpy(n-j+1,j,j,t,x,qty);
                  x(j,j) := temp;
             end if;
           end loop;
    end if;
    if cb                                   -- set up to compute b,rsd, or xb
     then zcopy(k,qty'first,qty,b);
    end if;
    kp1 := k + 1;
    if cxb then zcopy(k,qty'first,qty,xb); end if;
    if (cr and (k < n))
     then zcopy(n-k,kp1,qty,rsd);
    end if;
    if (cxb and (kp1 <= n))
     then for i in kp1..n loop
            xb(i) := Create(0.0);
          end loop;
    end if;
    if cr
     then for i in 1..k loop
            rsd(i) := Create(0.0);
          end loop;
    end if;
    if cb                                                        -- compute b
     then for j in 1..k loop
            jj := k - j + 1;
            if AbsVal(x(jj,jj)) = 0.0
             then info := jj; exit;
            end if;
            b(jj) := b(jj)/x(jj,jj);
            if jj /= 1
             then t := -b(jj);
                  zaxpy(jj-1,1,jj,t,x,b);
            end if;
		  end loop;
    end if;
    if cr or cxb                             -- compute rsd or xb as requested
	 then for j in 1..ju loop
            jj := ju - j + 1;
            if AbsVal(qraux(jj)) /= 0.0
             then temp := x(jj,jj);
                  x(jj,jj) := qraux(jj);
                  if cr
                   then t := -zdotc(jj,x,rsd)/x(jj,jj);
                        zaxpy(n-jj+1,jj,jj,t,x,rsd);
                  end if;
                  if cxb
                   then t := -zdotc(jj,x,xb)/x(jj,jj);
                        zaxpy(n-jj+1,jj,jj,t,x,xb);
                  end if;
                  x(jj,jj) := temp;
            end if;
          end loop;
    end if;
  end QRLS;

end Standard_Complex_Least_Squares;