File : scalvectors-ops.adb

package body ScalVectors.Ops is

  procedure Dgefa(A: in out Matrix; Pvt: out Pivot) is
--- adapted from linpack routine dgefa, http://www.netlib.org/linpack/
--- factors a matrix by Gaussian elimination, see Linear.Ops.dgefa for details
    N: constant Integer := Pvt'Last;
    L: Integer;
    X,Y: Radius;
    O,T: Scalar;
  begin
    SetZero(O);
    for K in 1 .. N-1 loop
      Y := Inf(abs(A(N,K)));
      L := N;
      for M in reverse K .. N-1 loop
        X := Inf(abs(A(M,K)));
        if X >= Y then Y := X; L := M; end if;
      end loop;
      Pvt(K) := L;
      if A(L,K) /= O then
        if L /= K then
          T := A(L,K);
          A(L,K) := A(K,K);
          A(K,K) := T;
        end if;
        T := -Inv(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;
    Pvt(N) := N;
  end Dgefa;

  procedure Dgesl(Pvt: in Pivot; A: in Matrix; B: in out Vector) is
--- adapted from linpack routine dgesl, http://www.netlib.org/linpack/
--- solves the system AX=B, see Linear.Ops.dgesl for details
    N: constant Integer := Pvt'Last;
    L: Integer;
    T: Scalar;
  begin
    for K in 1 .. N-1 loop
      L := Pvt(K);
      T := B(L);
      if L /= K then
        B(L) := B(K);
        B(K) := T;
      end if;
      for M in K+1 .. N loop
        B(M) := B(M)+T*A(M,K);
      end loop;
    end loop;
    for K in reverse 1 .. N loop
      B(K) := B(K)/A(K,K);
      T := -B(K);
      for M in 1 .. K-1 loop
        B(M) := B(M)+T*A(M,K);
      end loop;
    end loop;
  end Dgesl;

  procedure Dgedi(Pvt: in Pivot; A: in out Matrix; B: in out Vector) is
--- adapted from linpack routine dgedi, http://www.netlib.org/linpack/
--- computes the inverse of a matrix, see Linear.Ops.dgedi for details
    N: constant Integer := Pvt'Last;
    L: integer;
    O,T: Scalar;
  begin
    SetZero(O);
    for K in 1 .. N loop
      A(K,K) := Inv(A(K,K));
      T := -A(K,K);
      for M in 1 .. K-1 loop
        A(M,K) := T*A(M,K);
      end loop;
      for J in K+1 .. N loop
        T := A(K,J);
        A(K,J) := O;
        for M in 1 .. K loop
          A(M,J) := A(M,J)+T*A(M,K);
        end loop;
      end loop;
    end loop;
    for K in reverse 1 .. N-1 loop
      for I in K+1 .. N loop
        B(I) := A(I,K);
        A(I,K) := O;
      end loop;
      for J in K+1 .. N loop
        T := B(J);
        for M in 1 .. N loop
          A(M,K) := A(M,K)+T*A(M,J);
        end loop;
      end loop;
      L := Pvt(K);
      if L /= K  then
        for M in 1 .. N loop
          T := A(M,K);
          A(M,K) := A(M,L);
          A(M,L) := T;
        end loop;
      end if;
    end loop;
  end Dgedi;

  procedure Solve(A: in Matrix; B: in Vector; X: out Vector) is
    Pvt: Pivot(A'Range(1));
    Ar: Matrix(A'Range(1),A'Range(2));
  begin
    Ar := A;
    X := B;
    Dgefa(Ar,Pvt);
    Dgesl(Pvt,Ar,X);
  end Solve;

  procedure Invert(A: in out Matrix) is
    Pvt: Pivot(A'Range(2));
    B: Vector(A'Range(2));
  begin
    Dgefa(A,Pvt);
    Dgedi(Pvt,A,B);
  end Invert;

end ScalVectors.Ops;