File : lin-ops.adb

with Messages;
use Messages;

package body Lin.Ops is

  function "-"(U: Ivec) return Ivec is
    W: IVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := -U(I);
    end loop;
    return W;
  end "-";

  function "+"(U,V: IVec) return IVec is
    W: IVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := U(I)+V(I);
    end loop;
    return W;
  end "+";

  function "-"(U,V: IVec) return IVec is
    W: IVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := U(I)-V(I);
    end loop;
    return W;
  end "-";

  function "-"(U: Rvec) return Rvec is
    W: RVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := -U(I);
    end loop;
    return W;
  end "-";

  function "+"(U,V: RVec) return RVec is
    W: RVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := U(I)+V(I);
    end loop;
    return W;
  end "+";

  function "-"(U,V: RVec) return RVec is
    W: RVec(U'Range);
  begin
    for I in U'Range loop
      W(I) := U(I)-V(I);
    end loop;
    return W;
  end "-";

  function "*"(U,V: RVec) return Rep is
    R: Rep := Zero;
  begin
    for I in U'Range loop
      R := R+U(I)*V(I);
    end loop;
    return R;
  end "*";

  function "-"(A: RMat) return RMat is
    B: RMat(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        B(I,J) := -A(I,J);
      end loop;
    end loop;
    return B;
  end "-";

  function "+"(A,B: RMat) return RMat is
    C: RMat(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        C(I,J) := A(I,J)+B(I,J);
      end loop;
    end loop;
    return C;
  end "+";

  function "-"(A,B: RMat) return RMat is
    C: RMat(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        C(I,J) := A(I,J)-B(I,J);
      end loop;
    end loop;
    return C;
  end "-";

  function "*"(A: RMat; U: RVec) return RVec is
    R: Rep;
    V: RVec(A'Range(1));
  begin
    for I in V'Range loop
      R := Zero;
      for J in U'Range loop
        R := R+A(I,J)*U(J);
      end loop;
      V(I) := R;
    end loop;
    return V;
  end "*";

  function "*"(A,B: RMat) return RMat is
    R: Rep;
    C: RMat(A'Range(1),B'Range(2));
  begin
    for I in C'Range(1) loop
      for J in C'Range(2) loop
        R := Zero;
        for K in B'Range(1) loop
          R := R+A(I,K)*B(K,J);
        end loop;
        C(I,J) := R;
      end loop;
    end loop;
    return C;
  end "*";

  function "-"(A: QMat) return Qmat is
    D: constant Positive := A.P'Last;
    B: QMat(D);
  begin
    for I in 1 .. D loop
      for J in 1 .. D loop
        B.P(I,J) := -A.P(I,J);
      end loop;
    end loop;
    B.Q := A.Q;
    return B;
  end "-";

  function "*"(A: Qmat; U: IVec) return IVec is
    D: constant Positive := A.P'Last;
    R: Integer;
    V: IVec(1 .. D);
  begin
    for I in 1 .. D loop
      R := 0;
      for J in 1 .. D loop
        R := R+A.P(I,J)*U(J);
      end loop;
      if (R mod A.Q) /= 0 then
        Message("Lin.Ops.*",Not_Implemented);
      end if;
      V(I) := R/A.Q;
    end loop;
    return V;
  end "*";

  function Gcd(P: Positive; Q: Natural) return Positive is
    Pt: Positive;
    Qt,R: Natural;
  begin
    if Q>P then Pt := Q; Qt := P;
    else Pt := P; Qt := Q; end if;
    while Qt>0 loop
      R := Pt mod Qt;
      Pt := Qt;
      Qt := R;
    end loop;
    return Pt;
  end Gcd;

  procedure Cleanup(A: in out Qmat) is
    D: constant Positive := A.P'Last;
    G: Positive := A.Q;
  begin
    for I in 1 .. D loop
      for J in 1 .. D loop
        G := Gcd(G,abs(A.P(I,J)));
      end loop;
    end loop;
    for I in 1 .. D loop
      for J in 1 .. D loop
        A.P(I,J) := A.P(I,J)/G;
      end loop;
    end loop;
    A.Q := A.Q/G;
  end Cleanup;

  function "*"(A,B: QMat) return QMat is
    D: constant Positive := A.P'Last;
    R: Integer;
    C: QMat(D);
  begin
    for I in 1 .. D loop
      for J in 1 .. D loop
        R := 0;
        for K in 1 .. D loop
          R := R+A.P(I,K)*B.P(K,J);
        end loop;
        C.P(I,J) := R;
      end loop;
    end loop;
    C.Q := A.Q*B.Q;
    if C.Q>1 then Cleanup(C); end if;
    return C;
  end "*";

  function Sqr(A: QMat) return QMat is
    D: constant Positive := A.P'Last;
    R: Integer;
    C: QMat(D);
  begin
    for I in 1 .. D loop
      for J in 1 .. D loop
        R := 0;
        for K in 1 .. D loop
          R := R+A.P(I,K)*A.P(K,J);
        end loop;
        C.P(I,J) := R;
      end loop;
    end loop;
    C.Q := A.Q*A.Q;
    if C.Q>1 then Cleanup(C); end if;
    return C;
  end Sqr;

  procedure Divide(A: in out QMat; N: in Positive) is
  begin
    A.Q := N*A.Q;
  end Divide;

   procedure dgefa(n: in Integer; a: in out RMat; ipvt: out IVec) is
-- translated from linpack routine dgefa, http://www.netlib.org/linpack/
      l: Integer;
      t: Rep;
   begin
      for k in 1 .. n-1 loop
         t := abs(a(n,k));
         l := n;
         for m in reverse k .. n-1 loop
            if (abs(a(m,k)) >= t) then
               t := abs(a(m,k));
               l := m;
            end if;
         end loop;
         ipvt(k) := l;
         if (a(l,k) /= Zero) then
            if (l /= k) then
               t := a(l,k);
               a(l,k) := a(k,k);
               a(k,k) := t;
            end if;
            t := -One/a(k,k);
            for m in k+1 .. n loop
               a(m,k) := t*a(m,k);
            end loop;
            for j in k+1 .. n loop
               t := a(l,j);
               if (l /= k) then
                  a(l,j) := a(k,j);
                  a(k,j) := t;
               end if;
               for m in k+1 .. n loop
                  a(m,j) := a(m,j) + t*a(m,k);
               end loop;
            end loop;
         end if;
      end loop;
      ipvt(n) := n;
   end dgefa;

  procedure Det(N: in Integer; A: in out RMat; R: out Rep) is
    Ipvt: IVec(1 .. N);
  begin
    Dgefa(N,A,Ipvt);
    R := One;
    for I in 1 .. N loop
      if Ipvt(I)=I then
        R := R*A(I,I);
      else
        R := -R*A(I,I);
      end if;
    end loop;
  end Det;

end Lin.Ops;