-- cxvec.adb  from 535 general eigenvalue, eigenvector Fortran
with Ada.Numerics.Long_Complex_Types;
use  Ada.Numerics.Long_Complex_Types;
with Ada.numerics.generic_elementary_functions;
with Complex_Arrays; -- mine
use  Complex_Arrays; -- mine
with Real_Arrays;
use  Real_Arrays;
with Ada.Text_IO.Complex_IO; -- generic
-- with Complex_Arrays_IO; -- generic
with Ada.Text_IO;
use  Ada.Text_IO;

procedure cxvec(N:Integer; A: Complex_Matrix; B:in out Complex_Matrix;
                Z:in out Complex_Matrix; Alf: Complex_Vector;
                Beta: Real_vector) is

--    THIS SUBROUTINE ACCEPTS A PAIR OF COMPLEX MATRICES IN UPPER
--    TRIANGULAR FORM, WHERE ONE OF THEM FURTHER MUST HAVE REAL DIAGONAL
--    ELEMENTS.  IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM
--    AND TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM.
--    IT IS USUALLY PRECEDED BY  CQZHES  AND  CQZVAL.
--
--    ON INPUT-
--       N IS THE ORDER OF THE MATRICES,
--       A=(AR,AI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX,
--       B=(BR,BI) CONTAINS A COMPLEX UPPER TRIANGULAR MATRIX WITH REAL
--         DIAGONAL ELEMENTS.  IN ADDITION, LOCATION BR(N,1) CONTAINS
--         THE TOLERANCE QUANTITY (EPSB) COMPUTED AND SAVED IN  CQZVAL,
--       ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE
--         RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED
--         EIGENVALUES.  THEY ARE USUALLY OBTAINED FROM  CQZVAL,
--       Z=(ZR,ZI) CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
--         REDUCTIONS BY  CQZHES  AND  CQZVAL, IF PERFORMED.
--         IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE
--         DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX.
--
--    ON OUTPUT-
--       A IS UNALTERED,
--       B HAS BEEN DESTROYED,
--       ALFR, ALFI, AND BETA ARE UNALTERED,
--       Z CONTAINS THE EIGENVECTORS.  EACH EIGENVECTOR IS NORMALIZED
--         SO THAT THE MODULUS OF ITS LARGEST COMPONENT IS 1.0 .

  M,NA : Integer;
  AR,AI,BR,BI,ZR,ZI : Real_Matrix(1..N,1..N);
  ALFR,ALFI: Real_Vector(1..N);
  R,T,RI,TI,XI,ALMI,ALMR,BETM,EPSB : Real;
  Z3 : Complex;

  package Elementary_Functions is new
          Ada.numerics.Generic_Elementary_Functions(Real);
  use Elementary_Functions;

begin
  -- setup reals
  for I in 1..N loop
    for J in 1..N loop
      AR(I,J) := Re(A(I,J));
      AI(I,J) := Im(A(I,J));
      BR(I,J) := Re(B(I,J));
      BI(I,J) := Im(B(I,J));
      ZR(I,J) := Re(Z(I,J));
      ZI(I,J) := Im(Z(I,J));
    end loop;
    ALFR(I) := Re(ALF(I));
    ALFI(I) := Im(ALF(I));
  end loop;

  if N <= 1 then return; end if;
  EPSB := BR(N,1);
  -- FOR EN:=N STEP -1 UNTIL 2 DO --
  -- for NN in 2..N loop -- 800
  --   EN := N + 2 - NN;
  for EN in reverse 2..N loop -- 800
    NA := EN - 1;
    ALMR := ALFR(EN);
    ALMI := ALFI(EN);
    BETM := BETA(EN);
    -- FOR I = EN-1 STEP -1 UNTIL 1 DO --
    -- for II in 1..NA loop -- 700
    --   I := EN - II;
    for I in reverse 1..EN-1 loop -- 700
      R := 0.0;
      RI := 0.0;
      M := I + 1;
      for J in M..EN loop -- 610
        T := BETM * AR(I,J) - ALMR * BR(I,J) + ALMI * BI(I,J);
        TI := BETM * AI(I,J) - ALMR * BI(I,J) - ALMI * BR(I,J);
        if J /= EN then -- 605
          XI := T * BI(J,EN) + TI * BR(J,EN);
          T := T * BR(J,EN) - TI * BI(J,EN);
          TI := XI;
        end if; -- 605
        R := R + T;
        RI := RI + TI;
      end loop; -- 610

      T := ALMR * BETA(I) - BETM * ALFR(I);
      TI := ALMI * BETA(I) - BETM * ALFI(I);
      if T = 0.0 and TI = 0.0 then T := EPSB; end if;
      Z3 := (R,RI) / (T,TI);
      BR(I,EN) := Re(Z3);
      BI(I,EN) := Im(Z3);
    end loop; -- 700
  end loop; -- 800
  --  END BACK SUBSTITUTION.
  --  TRANSFORM TO ORIGINAL COORDINATE SYSTEM.
  -- FOR J:=N STEP -1 UNTIL 2 DO --
  -- for JJ in 2..N loop -- 880
  --   J := N + 2 - JJ
  for J in reverse 2..N loop -- 880
    M := J - 1;
    for I in 1..N loop -- 880
      for K in 1..M loop -- 860
        ZR(I,J) := ZR(I,J) + ZR(I,K) * BR(K,J) - ZI(I,K) * BI(K,J);
        ZI(I,J) := ZI(I,J) + ZR(I,K) * BI(K,J) + ZI(I,K) * BR(K,J);
      end loop; -- 860
    end loop; -- 880
  end loop; -- 880
  --  NORMALIZE SO THAT MODULUS OF LARGEST
  --  COMPONENT OF EACH VECTOR IS 1
  for J in 1..N loop -- 950
    T := 0.0;
    for I in 1..N loop -- 930
      R := sqrt(ZR(I,J)**2+ZI(I,J)**2);
      if R > T then T := R; end if;
    end loop; -- 930
    for I in 1..N loop -- 940
      ZR(I,J) := ZR(I,J) / T;
      ZI(I,J) := ZI(I,J) / T;
    end loop; -- 940
  end loop; -- 950

  for I in 1..N loop
    for J in 1..N loop
      Z(I,J) := (ZR(I,J),ZI(I,J));
    end loop;
  end loop;
end Cxvec;