File : vectors.adb

with Ints, Ints.IO;
use Ints;

pragma Elaborate (Ints);

package body Vectors is

  use Ints.IO;

  function IsZero(V: Vector) return Boolean is
  begin
    for N in V'Range loop
      if not IsZero(V(N)) then return False; end if;
    end loop;
    return True;
  end IsZero;

  procedure SetZero(V: out Vector) is
  begin
    for N in V'Range loop
      SetZero(V(N));
    end loop;
  end SetZero;

  function EffLast(V: Vector) return Integer is
  begin
    for N in reverse V'First+1 .. V'Last loop
      if not IsZero(V(N)) then return N; end if;
    end loop;
    return V'First;
  end EffLast;

  procedure Neg(V: in out Vector) is
  begin
    for N in V'First .. EffLast(V) loop
      Neg(V(N));
    end loop;
  end Neg;

  procedure Neg(V1: in Vector; V2: out Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
  begin
    for N in NFirst .. NLast loop
      Neg(V1(N),V2(N));
    end loop;
  end Neg;

  function "-"(V: Vector) return Vector is
    W: Vector(V'Range);
  begin
    for N in V'First .. EffLast(V) loop
      Neg(V(N),W(N));
    end loop;
    return W;
  end "-";

  procedure Add(V1: in Vector; V2: in out Vector) is
  begin
    for N in V1'First .. EffLast(V1) loop
      Add(V1(N),V2(N));
    end loop;
  end Add;

  procedure Sum(V1,V2: in Vector; V3: out Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First,V3'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last,V3'Last);
  begin
    for N in NFirst .. NLast loop
      Sum(V1(N),V2(N),V3(N));
    end loop;
  end Sum;

  function "+"(V1,V2: Vector) return Vector is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
    V3: Vector(V1'Range);
  begin
    for N in NFirst .. NLast loop
      Sum(V1(N),V2(N),V3(N));
    end loop;
    return V3;
  end "+";

  procedure Sub(V1: in Vector; V2: in out Vector) is
  begin
    for N in V1'First .. EffLast(V1) loop
      Sub(V1(N),V2(N));
    end loop;
  end Sub;

  procedure Diff(V1,V2: in Vector; V3: out Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First,V3'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last,V3'Last);
  begin
    for N in NFirst .. NLast loop
      Diff(V1(N),V2(N),V3(N));
    end loop;
  end Diff;

  function "-"(V1,V2: Vector) return Vector is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
    V3: Vector(V1'Range);
  begin
    for N in NFirst .. NLast loop
      Diff(V1(N),V2(N),V3(N));
    end loop;
    return V3;
  end "-";

  procedure Mult(K: in Integer; V: in out Vector) is
  begin
    for N in V'First .. EffLast(V) loop
      Mult(K,V(N));
    end loop;
  end Mult;

  procedure Prod(K: in Integer; V1: in Vector; V2: out Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
  begin
    for N in NFirst .. NLast loop
      Prod(K,V1(N),V2(N));
    end loop;
  end Prod;

  function "*"(K: Integer; V: Vector) return Vector is
    W: Vector(V'Range);
  begin
    for N in V'First .. EffLast(V) loop
      Prod(K,V(N),W(N));
    end loop;
    return W;
  end "*";

----------------

  procedure Mult(S: in Scalar; V: in out Vector) is
  begin
    for N in V'First .. EffLast(V) loop
      Mult(S,V(N));
    end loop;
  end Mult;

  procedure Prod(S: in Scalar; V1: in Vector; V2: out Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
  begin
    for N in NFirst .. NLast loop
      Prod(S,V1(N),V2(N));
    end loop;
  end Prod;

  function "*"(S: Scalar; V: Vector) return Vector is
    W: Vector(V'Range);
  begin
    for N in V'First .. EffLast(V) loop
      Prod(S,V(N),W(N));
    end loop;
    return W;
  end "*";

  procedure MultAdd(S: Scalar; V1: in Vector; V2: in out Vector) is
  begin
    for N in V1'First .. EffLast(V1) loop
      MultAdd(S,V1(N),V2(N));
    end loop;
  end MultAdd;

  procedure Show1(N: in String; V: in Vector) is
    W: constant Positive := 5;
  begin
    for K in V'Range loop
      if not IsZero(V(K)) then
        Show1(N & Strng(K,W),V(K));
      end if;
    end loop;
    New_Line;
  end Show1;

  procedure Show2(N: in String; V1,V2: in Vector) is
    NFirst: constant Integer := IMin(V1'First,V2'First);
    NLast:  constant Integer := IMax(V1'Last,V2'Last);
    W: constant Positive := 5;
  begin
    for K in NFirst .. NLast loop
      if not (IsZero(V1(K)) and IsZero(V2(K))) then
        Show2(N & Strng(K,W),V1(K),V2(K));
      end if;
    end loop;
    New_Line;
  end Show2;

  procedure Get(F: in File_Type; V: out Vector; Decimal: in Boolean := False) is
    Trunc: Boolean := False;
    I,L: Integer;
    C: Component;
  begin
    SetZero(V);
    Get(F,L);
    I := V'First;
    for M in 1 .. L loop
      Get(F,C,Decimal);
      if I <= V'Last then V(I) := C; else Trunc := True; end if;
      I := I+1;
    end loop;
    if Trunc then
      Put_Line("Vectors.Get Warning: input vector truncated");
    end if;
  end Get;

  procedure Put(F: in File_Type; V: in Vector; Decimal: in Boolean := False) is
    VLast: constant Integer := EffLast(V);
    VLength: constant Integer := VLast-V'First+1;
  begin
    Put(F,VLength);
    for I in V'First .. VLast loop
      Put(F,V(I),Decimal);
    end loop;
  end Put;

  procedure Read(Name: in String; V: out Vector; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Reading " & Name);
    Open(F,In_File,Name);
    Get(F,V,Decimal);
    Close(F);
  end Read;

  procedure Write(Name: in String; V: in Vector; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Writing " & Name);
    Create(F,Out_File,Name);
    Put(F,V,Decimal);
    Close(F);
  end Write;

----------------------------------------------

  function IsZero(A: Matrix) return Boolean is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        if not IsZero(A(I,J)) then return False; end if;
      end loop;
    end loop;
    return True;
  end IsZero;

  procedure SetZero(A: out Matrix) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        SetZero(A(I,J));
      end loop;
    end loop;
  end SetZero;

  procedure ZeroRow(I: in Integer; A: in out Matrix) is
  begin
    for J in A'Range(2) loop
      SetZero(A(I,J));
    end loop;
  end ZeroRow;

  procedure ZeroColumn(J: in Integer; A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      SetZero(A(I,J));
    end loop;
  end ZeroColumn;

  procedure AddConst(R: Rep; A: in out Matrix) is
    AFirst: constant Integer := IMin(A'First(1),A'First(2));
    ALast: constant Integer := IMax(A'Last(1),A'Last(2));
  begin
    for K in AFirst .. ALast loop
      Add(Scal(R),A(K,K));
    end loop;
  end AddConst;

  function EffLast1(A: Matrix) return Integer is
  begin
    for I in reverse A'First(1)+1 .. A'Last(1) loop
      for J in A'Range(2) loop
        if not IsZero(A(I,J)) then return I; end if;
      end loop;
    end loop;
    return A'First(1);
  end EffLast1;

  function EffLast2(A: Matrix) return Integer is
  begin
    for J in reverse A'First(2)+1 .. A'Last(2) loop
      for I in A'Range(1) loop
        if not IsZero(A(I,J)) then return J; end if;
      end loop;
    end loop;
    return A'First(2);
  end EffLast2;

  procedure Transpose(A: in out Matrix) is
    AFirst: constant Integer := IMin(A'First(1),A'First(2));
    ALast: constant Integer := IMax(A'Last(1),A'Last(2));
    S: Scalar;
  begin
    for I in AFirst .. ALast loop
      for J in AFirst .. I-1 loop
        S := A(I,J);
        A(I,J) := A(J,I);
        A(J,I) := S;
      end loop;
    end loop;
  end Transpose;

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

  procedure Neg(A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Neg(A(I,J));
      end loop;
    end loop;
  end Neg;

  procedure Neg(A1: in Matrix; A2: out Matrix) is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Neg(A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Neg;

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

  procedure Add(A1: in Matrix; A2: in out Matrix) is
  begin
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        Add(A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Add;

  procedure Sum(A1,A2: in Matrix; A3: out Matrix) is
    IFirst: constant Integer := IMin(A1'First(1),A2'First(1),A3'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1),A3'Last(1));
    JFirst: constant Integer := IMin(A1'First(2),A2'First(2),A3'First(2));
    JLast:  constant Integer := IMax(A1'Last(2),A2'Last(2),A3'Last(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Sum(A1(I,J),A2(I,J),A3(I,J));
      end loop;
    end loop;
  end Sum;

  function "+"(A1,A2: Matrix) return Matrix is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
    A3: Matrix(A1'Range(1),A1'Range(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Sum(A1(I,J),A2(I,J),A3(I,J));
      end loop;
    end loop;
    return A3;
  end "+";

  procedure Sub(A1: in Matrix; A2: in out Matrix) is
  begin
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        Sub(A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Sub;

  procedure Diff(A1,A2: in Matrix; A3: out Matrix) is
    IFirst: constant Integer := IMin(A1'First(1),A2'First(1),A3'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1),A3'Last(1));
    JFirst: constant Integer := IMin(A1'First(2),A2'First(2),A3'First(2));
    JLast:  constant Integer := IMax(A1'Last(2),A2'Last(2),A3'Last(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Diff(A1(I,J),A2(I,J),A3(I,J));
      end loop;
    end loop;
  end Diff;

  function "-"(A1,A2: Matrix) return Matrix is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
    A3: Matrix(A1'Range(1),A1'Range(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Diff(A1(I,J),A2(I,J),A3(I,J));
      end loop;
    end loop;
    return A3;
  end "-";

  procedure Mult(K: in Integer; A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Mult(K,A(I,J));
      end loop;
    end loop;
  end Mult;

  procedure Prod(K: in Integer; A1: in Matrix; A2: out Matrix) is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Prod(K,A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Prod;

  function "*"(K: Integer; A: Matrix) return Matrix is
    B: Matrix(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Prod(K,A(I,J),B(I,J));
      end loop;
    end loop;
    return B;
  end "*";

---------------

  procedure Mult(S: in Scalar; A: in out Matrix) is
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Mult(S,A(I,J));
      end loop;
    end loop;
  end Mult;

  procedure Prod(S: in Scalar; A1: in Matrix; A2: out Matrix) is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        Prod(S,A1(I,J),A2(I,J));
      end loop;
    end loop;
  end Prod;

  function "*"(S: Scalar; A: Matrix) return Matrix is
    B: Matrix(A'Range(1),A'Range(2));
  begin
    for I in A'Range(1) loop
      for J in A'Range(2) loop
        Prod(S,A(I,J),B(I,J));
      end loop;
    end loop;
    return B;
  end "*";

  procedure MultAdd(S: in Scalar; A1: in Matrix; A2: in out Matrix) is
  begin
    for I in A1'Range(1) loop
      for J in A1'Range(2) loop
        MultAdd(S,A1(I,J),A2(I,J));
      end loop;
    end loop;
  end MultAdd;

------------

  procedure Mult(A: in Matrix; V: in out Vector) is
    W: constant Vector := V;
  begin
    Prod(A,W,V);
  end Mult;

  procedure Prod(A: in Matrix; V1: in Vector; V2: out Vector) is
    IFirst: constant Integer := IMin(A'First(1),V2'First);
    ILast:  constant Integer := IMax(A'Last(1),V2'Last);
    JFirst: constant Integer := IMin(A'First(2),V1'First);
    JLast:  constant Integer := IMax(A'Last(2),V1'Last);
    S: Component;
  begin
    for I in IFirst .. ILast loop
      SetZero(S);
      for J in JFirst .. JLast loop
        MultAdd(A(I,J),V1(J),S);
      end loop;
      V2(I) := S;
    end loop;
  end Prod;

  function "*"(A: Matrix; V: Vector) return Vector is
    JFirst: constant Integer := IMin(A'First(2),V'First);
    JLast:  constant Integer := IMax(A'Last(2),V'Last);
    S: Component;
    W: Vector(A'Range(1));
  begin
    for I in A'Range(1) loop
      SetZero(S);
      for J in JFirst .. JLast loop
        MultAdd(A(I,J),V(J),S);
      end loop;
      W(I) := S;
    end loop;
    return W;
  end "*";

  procedure MultAdd(A: in Matrix; V1: in Vector; V2: in out Vector) is
    JFirst: constant Integer := IMin(A'First(2),V1'First);
    JLast:  constant Integer := IMax(A'Last(2),V1'Last);
    S: Component;
  begin
    for I in A'Range(1) loop
      S := V2(I);
      for J in JFirst .. JLast loop
        MultAdd(A(I,J),V1(J),S);
      end loop;
      V2(I) := S;
    end loop;
  end MultAdd;

  procedure Mult(A1: in Matrix; A2: in out Matrix) is
    B: constant Matrix := A2;
  begin
    Prod(A1,B,A2);
  end Mult;

  procedure Prod(A1,A2: in Matrix; A3: out Matrix) is
    IFirst: constant Integer := IMin(A1'First(1),A3'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A3'Last(1));
    JFirst: constant Integer := IMin(A2'First(2),A3'First(2));
    JLast:  constant Integer := IMax(A2'Last(2),A3'Last(2));
    KFirst: constant Integer := IMin(A1'First(2),A2'First(1));
    KLast:  constant Integer := IMax(A1'Last(2),A2'Last(1));
    S: Scalar;
  begin
    for I in IFirst .. ILast loop
      for J in JFirst .. JLast loop
        SetZero(S);
        for K in KFirst .. KLast loop
          MultAdd(A1(I,K),A2(K,J),S);
        end loop;
        A3(I,J) := S;
      end loop;
    end loop;
  end Prod;

  function "*"(A1,A2: Matrix) return Matrix is
    KFirst: constant Integer := IMin(A1'First(2),A2'First(1));
    KLast:  constant Integer := IMax(A1'Last(2),A2'Last(1));
    S: Scalar;
    A3: Matrix(A1'Range(1),A2'Range(2));
  begin
    for I in A1'Range(1) loop
      for J in A2'Range(2) loop
        SetZero(S);
        for K in KFirst .. KLast loop
          MultAdd(A1(I,K),A2(K,J),S);
        end loop;
        A3(I,J) := S;
      end loop;
    end loop;
    return A3;
  end "*";

  procedure MultAdd(A1,A2: in Matrix; A3: in out Matrix) is
    KFirst: constant Integer := IMin(A1'First(2),A2'First(1));
    KLast:  constant Integer := IMax(A1'Last(2),A2'Last(1));
    S: Scalar;
  begin
    for I in A1'Range(1) loop
      for J in A2'Range(2) loop
        S := A3(I,J);
        for K in KFirst .. KLast loop
          MultAdd(A1(I,K),A2(K,J),S);
        end loop;
        A3(I,J) := S;
      end loop;
    end loop;
  end MultAdd;

-------

  procedure Show1(N: in String; A: in Matrix) is
    W: constant Positive := 3;
    Z: Boolean;
  begin
    for I in A'Range(1) loop
      Z := True;
      for J in A'Range(2) loop
        if not IsZero(A(I,J)) then
          Z := False;
          Show1(N & Strng(I,W) & Strng(J,W),A(I,J));
        end if;
      end loop;
      if not Z then New_Line; end if;
    end loop;
  end Show1;

  procedure Show2(N: in String; A1,A2: in Matrix) is
    IFirst: constant Integer := Imin(A1'First(1),A2'First(1));
    ILast:  constant Integer := IMax(A1'Last(1),A2'Last(1));
    JFirst: constant Integer := Imin(A1'First(2),A2'First(2));
    JLast:  constant Integer := Imax(A1'Last(2),A2'Last(2));
    W: constant Positive := 5;
    Z: Boolean;
  begin
    for I in IFirst .. ILast loop
      Z := True;
      for J in JFirst .. JLast loop
        if not (IsZero(A1(I,J)) and IsZero(A2(I,J))) then
          Z := False;
          Show2(N & Strng(I,W) & Strng(J,W),A1(I,J),A2(I,J));
        end if;
      end loop;
      if not Z then New_Line; end if;
    end loop;
  end Show2;

  procedure Get(F: in File_Type; A: out Matrix; Decimal: in Boolean := True) is
    Trunc: Boolean := False;
    I,J,LM,LN: Integer;
    S: Scalar;
  begin
    SetZero(A);
    Get(F,LM);
    Get(F,LN);
    I := A'First(1);
    for M in 1 .. LM loop
      J := A'First(2);
      for N in 1 .. LN loop
        Get(F,S,Decimal);
        if I <= A'Last(1) and then J <= A'Last(2) then A(I,J) := S; else Trunc := True; end if;
        J := J+1;
      end loop;
      I := I+1;
    end loop;
    if Trunc then
      Put_Line("Vectors.Get Warning: input matrix truncated");
    end if;
  end Get;

  procedure Put(F: in File_Type; A: in Matrix; Decimal: in Boolean := True) is
    ALast1: constant Integer := EffLast1(A);
    ALast2: constant Integer := EffLast2(A);
    ALength1: constant Integer := ALast1-A'First(1)+1;
    ALength2: constant Integer := ALast2-A'First(2)+1;
  begin
    Put(F,ALength1);
    Put(F,ALength2);
    for I in A'First(1) .. ALast1 loop
      for J in A'First(2) .. ALast2 loop
        Put(F,A(I,J),Decimal);
      end loop;
    end loop;
  end Put;

  procedure Read(Name: in String; A: out Matrix; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Reading " & Name);
    Open(F,In_File,Name);
    Get(F,A,Decimal);
    Close(F);
  end Read;

   procedure Write(Name: in String; A: in Matrix; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Put_Line("Writing " & Name);
    Create(F,Out_File,Name);
    Put(F,A,Decimal);
    Close(F);
  end Write;

end Vectors;