Return to standard_floating_least_squares.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Standard_Floating_Numbers;          use Standard_Floating_Numbers;

package body Standard_Floating_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 dcopy ( n,start : in natural;
                    x : in Standard_Floating_Vectors.Vector;
                    y : out Standard_Floating_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 dcopy;

  function ddot ( row : natural;
                  x : Standard_Floating_Matrices.Matrix;
                  y : Standard_Floating_Vectors.Vector )
                return double_float is

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

    res : double_float := 0.0;

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

  procedure daxpy ( n,row,col : in natural; f : in double_float;
                    x : in Standard_Floating_Matrices.Matrix;
                    y : in out Standard_Floating_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 daxpy;

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

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

  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 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) := 0.0; end if;
          return;
    end if;
    if cqy                                     -- set up to compute qy or qty
     then dcopy(n,y'first,y,qy);
    end if;
    if cqty
     then dcopy(n,y'first,y,qty);
    end if;
    if cqy                                                      -- compute qy
     then for j in 1..ju loop
            jj := ju - j + 1;
            if qraux(jj) /= 0.0
             then temp := x(jj,jj);
                  x(jj,jj) := qraux(jj);
                  t := -ddot(jj,x,qy)/x(jj,jj);
                  daxpy(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 qraux(j) /= 0.0
             then temp := x(j,j);
                  x(j,j) := qraux(j);
                  t := -ddot(j,x,qty)/x(j,j);
                  daxpy(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 dcopy(k,qty'first,qty,b);
    end if;
    kp1 := k + 1;
    if cxb then dcopy(k,qty'first,qty,xb); end if;
    if (cr and (k < n))
     then dcopy(n-k,kp1,qty,rsd);
    end if;
    if (cxb and (kp1 <= n))
     then for i in kp1..n loop
            xb(i) := 0.0;
          end loop;
    end if;
    if cr
     then for i in 1..k loop
            rsd(i) := 0.0;
          end loop;
    end if;
    if cb                                                        -- compute b
     then for j in 1..k loop
            jj := k - j + 1;
            if 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);
                  daxpy(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 qraux(jj) /= 0.0
             then temp := x(jj,jj);
                  x(jj,jj) := qraux(jj);
                  if cr
                   then t := -ddot(jj,x,rsd)/x(jj,jj);
                        daxpy(n-jj+1,jj,jj,t,x,rsd);
                  end if;
                  if cxb
                   then t := -ddot(jj,x,xb)/x(jj,jj);
                        daxpy(n-jj+1,jj,jj,t,x,xb);
                  end if;
                  x(jj,jj) := temp;
            end if;
          end loop;
    end if;
  end QRLS;

end Standard_Floating_Least_Squares;