-- cxhes.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 cxhes(N:Integer; A:in out Complex_Matrix; B:in out Complex_Matrix; Z:in out Complex_Matrix) is -- A=(AR,AI) CONTAINS A COMPLEX GENERAL MATRIX, -- B=(BR,BI) CONTAINS A COMPLEX GENERAL MATRIX, -- ON OUTPUT- -- A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS -- BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO, AND THE -- SUBDIAGONAL ELEMENTS HAVE BEEN MADE REAL (AND NON-NEGATIVE), -- -- B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS -- BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, -- -- Z=(ZR,ZI) CONTAINS THE PRODUCT OF THE RIGHT HAND -- TRANSFORMATIONS L,K1,L1,NK1,NM1:Integer; R,S,T,TI,U1,U2,XI,XR,YI,YR,RHO,U1I:Real; package Elementary_Functions is new Ada.numerics.Generic_Elementary_Functions(Real); use Elementary_Functions; begin for I in 1..N loop -- 3 for J in 1..N loop -- 2 Z(I,J) := (0.0,0.0); end loop; -- 2 Z(I,I) := (1.0,0.0); end loop; -- 3 -- REDUCE B TO UPPER TRIANGULAR FORM WITH -- TEMPORARILY REAL DIAGONAL ELEMENTS if N <= 1 then Return; end if; NM1 := N - 1; for L in 1..NM1 loop -- 100 L1 := L + 1; S := 0.0; for I in L..N loop -- 20 S := S + abs(Re(B(I,L))) + abs(Im(B(I,L))); end loop; -- 20 if S/=0.0 then -- goto l100 RHO := 0.0; for I in L..N loop -- 25 B(I,L) := B(I,L) / S; RHO := RHO + Re(B(I,L))**2 + Im(B(I,L))**2; end loop; -- 25 R := Sqrt(RHO); XR := abs(B(L,L)); if XR /= 0.0 then RHO := RHO + XR * R; U1 := -Re(B(L,L)) / XR; U1I := -Im(B(L,L)) / XR; YR := R / XR + 1.0; B(L,L) := YR*B(L,L); else -- 27 Set_Re(B(L,L),R); U1 := -1.0; U1I := 0.0; end if; -- 28 for J in L1..N loop -- 50 T := 0.0; TI := 0.0; for I in L..N loop -- 30 T := T + Re(B(I,L)) * Re(B(I,J)) + Im(B(I,L)) * Im(B(I,J)); TI := TI + Re(B(I,L)) * Im(B(I,J)) - Im(B(I,L)) * Re(B(I,J)); end loop; -- 30 T := T / RHO; TI := TI / RHO; for I in L..N loop -- 40 Set_Re(B(I,J), Re(B(I,J)) - T * Re(B(I,L)) + TI * Im(B(I,L))); Set_Im(B(I,J), Im(B(I,J)) - T * Im(B(I,L)) - TI * Re(B(I,L))); end loop; -- 40 XI := U1 * Im(B(L,J)) - U1I * Re(B(L,J)); Set_Re(B(L,J), U1 * Re(B(L,J)) + U1I * Im(B(L,J))); Set_Im(B(L,J), XI); end loop; -- 50 for J in 1..N loop -- 80 T := 0.0; TI := 0.0; for I in L..N loop --60 T := T + Re(B(I,L)) * Re(A(I,J)) + Im(B(I,L)) * Im(A(I,J)); TI := TI + Re(B(I,L)) * Im(A(I,J)) - Im(B(I,L)) * Re(A(I,J)); end loop; -- 60 T := T / RHO; TI := TI / RHO; for I in L..N loop -- 70 Set_Re(A(I,J), Re(A(I,J)) - T * Re(B(I,L)) + TI * Im(B(I,L))); Set_Im(A(I,J), Im(A(I,J)) - T * Im(B(I,L)) - TI * Re(B(I,L))); end loop; -- 70 XI := U1 * Im(A(L,J)) - U1I * Re(A(L,J)); Set_Re(A(L,J), U1 * Re(A(L,J)) + U1I * Im(A(L,J))); Set_Im(A(L,J), XI); end loop; -- 80 B(L,L) := (R*S,0.0); for I in L1..N loop -- 90 B(I,L) := (0.0,0.0); end loop; -- 90 end if; end loop; --100 -- REDUCE A TO UPPER HESSENBERG FORM WITH REAL SUBDIAGONAL -- ELEMENTS, WHILE KEEPING B TRIANGULAR for K in 1..NM1 loop -- 160 K1 := K + 1; -- SET BOTTOM ELEMENT IN K-TH COLUMN OF A REAL if Im(A(N,K)) /= 0.0 then -- 105 R := abs(A(N,K)); U1 := Re(A(N,K)) / R; U1I := Im(A(N,K)) / R; A(N,K) := (R,0.0); for J in K1..N loop --103 XI := U1 * Im(A(N,J)) - U1I * Re(A(N,J)); Set_Re(A(N,J), U1 * Re(A(N,J)) + U1I * Im(A(N,J))); Set_Im(A(N,J), XI); end loop; -- 103 XI := U1 * Im(B(N,N)) - U1I * Re(B(N,N)); Set_Re(B(N,N), U1 * Re(B(N,N)) + U1I * Im(B(N,N))); Set_Im(B(N,N), XI); end if; -- 105 if K = NM1 then return; end if; NK1 := NM1 - K; -- FOR L=N-1 STEP -1 UNTIL K+1 DO for LB in 1..NK1 loop -- 150 L := N - LB; L1 := L + 1; -- ZERO A(L+1,K) S := abs(Re(A(L,K))) + abs(Im(A(L,K))) + Re(A(L1,K)); if S /= 0.0 then -- goto 150 U1 := Re(A(L,K)) / S; U1I := Im(A(L,K)) / S; U2 := Re(A(L1,K)) / S; R := sqrt(U1*U1+U1I*U1I+U2*U2); U1 := U1 / R; U1I := U1I / R; U2 := U2 / R; A(L,K) := (R*S,0.0); Set_Re(A(L1,K), 0.0); for J in K1..N loop -- 110 XR := Re(A(L,J)); XI := Im(A(L,J)); YR := Re(A(L1,J)); YI := Im(A(L1,J)); Set_Re(A(L,J), U1 * XR + U1I * XI + U2 * YR); Set_Im(A(L,J), U1 * XI - U1I * XR + U2 * YI); Set_Re(A(L1,J), U1 * YR - U1I * YI - U2 * XR); Set_Im(A(L1,J), U1 * YI + U1I * YR - U2 * XI); end loop; -- 110 XR := Re(B(L,L)); Set_Re(B(L,L), U1 * XR); Set_Im(B(L,L), -U1I * XR); Set_Re(B(L1,L), -U2 * XR); for J in L1..N loop -- 120 XR := Re(B(L,J)); XI := Im(B(L,J)); YR := Re(B(L1,J)); YI := Im(B(L1,J)); Set_Re(B(L,J), U1 * XR + U1I * XI + U2 * YR); Set_Im(B(L,J), U1 * XI - U1I * XR + U2 * YI); Set_Re(B(L1,J), U1 * YR - U1I * YI - U2 * XR); Set_Im(B(L1,J), U1 * YI + U1I * YR - U2 * XI); end loop; --120 -- ZERO B(L+1,L) S := abs(Re(B(L1,L1))) + abs(Im(B(L1,L1))) + abs(Re(B(L1,L))); if S /= 0.0 then -- goto 150 U1 := Re(B(L1,L1)) / S; U1I := Im(B(L1,L1)) / S; U2 := Re(B(L1,L)) / S; R := sqrt(U1*U1+U1I*U1I+U2*U2); U1 := U1 / R; U1I := U1I / R; U2 := U2 / R; B(L1,L1) := (R*S,0.0); Set_Re(B(L1,L), 0.0); for I in 1..L loop -- 130 XR := Re(B(I,L1)); XI := Im(B(I,L1)); YR := Re(B(I,L)); YI := Im(B(I,L)); Set_Re(B(I,L1), U1 * XR + U1I * XI + U2 * YR); Set_Im(B(I,L1), U1 * XI - U1I * XR + U2 * YI); Set_Re(B(I,L), U1 * YR - U1I * YI - U2 * XR); Set_Im(B(I,L), U1 * YI + U1I * YR - U2 * XI); end loop; -- 130 for I in 1..N loop -- 140 XR := Re(A(I,L1)); XI := Im(A(I,L1)); YR := Re(A(I,L)); YI := Im(A(I,L)); Set_Re(A(I,L1), U1 * XR + U1I * XI + U2 * YR); Set_Im(A(I,L1), U1 * XI - U1I * XR + U2 * YI); Set_Re(A(I,L), U1 * YR - U1I * YI - U2 * XR); Set_Im(A(I,L), U1 * YI + U1I * YR - U2 * XI); end loop; -- 140 for I in 1..N loop -- 145 XR := Re(Z(I,L1)); XI := Im(Z(I,L1)); YR := Re(Z(I,L)); YI := Im(Z(I,L)); Set_Re(Z(I,L1), U1 * XR + U1I * XI + U2 * YR); Set_Im(Z(I,L1), U1 * XI - U1I * XR + U2 * YI); Set_Re(Z(I,L), U1 * YR - U1I * YI - U2 * XR); Set_Im(Z(I,L), U1 * YI + U1I * YR - U2 * XI); end loop; -- 145 end if; -- goto 150 end if; -- other goto 150 end loop; --150 end loop; -- 160 end Cxhes;