-- cxcheck.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;

procedure cxcheck(N:Integer; A: Complex_Matrix; B:in out Complex_Matrix;
                Z:in out Complex_Matrix; Alf: Complex_Vector;
                Beta: Real_Vector; Norm: out Real_Vector; Resdul: out Real) is

--   THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
--   A*Z-B*Z*DIAG(W) WHERE A AND B ARE COMPLEX GENERAL MATRICES, Z IS
--   A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM
--   A*Z - B*Z*DIAG(W), AND  W  STANDS FOR A VECTOR OF CORRESPONDING
--   EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS  ALFR,
--   ALFI,  AND  BETA  BY THE CORRESPONDENCES
--      W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) .
--   ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
--   INPUT-
--      N IS THE ORDER OF THE MATRICES  A  AND  B,
--      A(NM,N),AI(NM,N),B(NM,N), AND BI(NM,N)  ARE ARRAYS WHICH
--        CONTAIN THE MATRICES OF THE SYSTEM,
--      Z(NM,N) AND ZI(NM,N)  ARE ARRAYS WHICH CONTAIN THE
--        EIGENVECTORS OF THE SYSTEM,
--      ALFR(N), ALFI(N), AND BETA(N)  ARE ARRAYS CONTAINING THE
--        COMPONENTS OF THE EIGENVALUES OF THE SYSTEM.
--   OUTPUT-
--      Z(NM,N) AND ZI(NM,N)  ARE ARRAYS WHICH CONTAIN THE NORMALIZED
--        APPROXIMATE EIGENVECTORS OF THE SYSTEM.  THE EIGENVECTORS
--        ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST
--        ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
--        EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE,
--      NORM(N) IS AN ARRAY SUCH THAT FOR EACH K,
--                   ..BETA(K)*A*Z(K)-ALFA*B*Z(K)..
--        NORM(K) = -----------------------------------------------
--                   ..Z(K)..*(BETA(K)*..A..+.ALFA(K).*..B..)
--
--        WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI),
--      RESDUL IS THE REAL NUMBER
--        ..BETA*A*Z-ALFA*B*Z../(..Z..*(BETA*..A..+.ALFA.*..B..))

  AR,AI,BR,BI,ZR,ZI : Real_Matrix(1..N,1..N);
  ALFR,ALFI: Real_Vector(1..N);
  XR, XI, S, SUMZ, SUMR, SUMI, NORMAB,SUMA,SUMB,NORMA,NORMB : Long_Float;
  SUMR2,SUMI2,SUMR3,SUMI3 : Long_Float;
  C, C1 : Complex;

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

  function Max(X,Y: Long_Float) return Long_Float is
  begin
    if X>Y then
      return X;
    end if;
    return Y;
  end Max;

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;

  NORMB := 0.0;
  NORMA := 0.0;
  RESDUL := 0.0;
  for I in 1..N loop -- 20
    SUMA := 0.0;
    SUMB := 0.0;
    for L in 1..N loop -- 10
      SUMA := SUMA + abs(AR(L,I)) + abs(AI(L,I));
      SUMB := SUMB + abs(BR(L,I)) + abs(BI(L,I));
    end loop; -- 10
    NORMA := max(NORMA,SUMA);
    NORMB := max(NORMB,SUMB);
  end loop; -- 20

  for I in 1..N loop -- 160
    S := 0.0;
    SUMZ := 0.0;
    for L in 1..N loop -- 110
      SUMR := 0.0;
      SUMI := 0.0;
      SUMR2 := 0.0;
      SUMI2 := 0.0;
      SUMZ := SUMZ + sqrt(ZR(L,I)**2+ZI(L,I)**2);
      for K in 1..N loop -- 100
        SUMR := SUMR + BR(L,K)*ZR(K,I) - BI(L,K)*ZI(K,I);
        SUMI := SUMI + BR(L,K)*ZI(K,I) + BI(L,K)*ZR(K,I);
        SUMR2 := SUMR2 + AR(L,K)*ZR(K,I) - AI(L,K)*ZI(K,I);
        SUMI2 := SUMI2 + AR(L,K)*ZI(K,I) + AI(L,K)*ZR(K,I);
      end loop; -- 100
      SUMR3 := -ALFR(I)*SUMR + ALFI(I)*SUMI;
      SUMI3 := -ALFI(I)*SUMR - ALFR(I)*SUMI;
      S := S+abs((SUMR3,SUMI3)+(SUMR2,SUMI2)*BETA(I));
    end loop; -- 110
    NORMAB := NORMA*BETA(I)+NORMB*sqrt(ALFR(I)**2+ALFI(I)**2);
    if NORMAB = 0.0 then NORMAB := 1.0; end if;
    NORM(I) := SUMZ;
    if SUMZ = 0.0 then goto L150; end if;
       --  THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL
       --  ALWAYS EXIST AN ELEMENT IN THE VECTOR (Z(I),ZI(I))
       --  LARGER THAN ..(Z(I),ZI(I)../N
    for L in 1..N loop -- 120
       if sqrt(ZR(L,I)**2+ZI(L,I)**2) >= NORM(I)/Long_float(N) then
         XR := NORM(I)*ZR(L,I)/sqrt(ZR(L,I)**2+ZI(L,I)**2);
         XI := NORM(I)*ZI(L,I)/sqrt(ZR(L,I)**2+ZI(L,I)**2);
         C := (XR,XI);
         exit;
       end if;
    end loop; -- 120

    for L in 1..N loop -- 140
      C1 := (ZR(L,I),ZI(L,I))/C;
      ZR(L,I) := Re(C1);
      ZI(L,I) := Im(C1);
    end loop; -- 140
    NORM(I) := S/(NORM(I)*NORMAB);
    <<L150>>
    RESDUL := max(NORM(I),RESDUL);
  end loop; -- 160
end cxcheck;