
package body EisPack2 is

   procedure cdiv(ar,ai,br,bi: in Real; cr,ci: in out Real) is
-- complex division, eispack style
     s: Real := abs(br)+abs(bi);
     ars: constant Real := ar/s;
     ais: constant Real := ai/s;
     brs: constant Real := br/s;
     bis: constant Real := bi/s;
   begin
     s := brs*brs+bis*bis;
     cr := (ars*brs+ais*bis)/s;
     ci := (ais*brs-ars*bis)/s;
   end cdiv;

   procedure csroot(xr,xi: in Real; yr,yi: in out Real) is
-- complex square root, eispack style
-- branch chosen so that yr>=0 and sign(yi)=sign(xi)
      s: Real;
   begin
      s := Sqrt(Half_Real*(Sqrt(xr*xr + xi*xi)+abs(xr)));
      if not (xr < Zero_Real) then yr := s; end if;
      if (xi < Zero_Real) then s := -s; end if;
      if not (Zero_Real < xr) then yi := s; end if;
      if (xr <  Zero_Real) then yr := Half_Real*(xi/yi); end if;
      if (Zero_Real < xr) then yi := Half_Real*(xi/yr); end if;
   end csroot;

   procedure cbal(n: in Integer; ar,ai: in out Matrix; low,igh: out Integer; scale: in out Vector) is
-- translated from eispack routine cbal, http://www.netlib.org/eispack/
     allzero,noconv: Boolean;
     k,l,m: Integer;
     c,f,g,r,s,b2,radix: Real;

--     .......... row and column exchange ..........
      procedure exc(m,j: in integer) is
      begin
         scale(m) := Scal(Flt(j));
         if (j = m) then return; end if;
         for i in 1 .. l loop
            Swap(ar(i,j),ar(i,m));
            Swap(ai(i,j),ai(i,m));
         end loop;
         for i in k .. n loop
            Swap(ar(j,i),ar(m,i));
            Swap(ai(j,i),ai(m,i));
         end loop;
      end exc;

   begin
      radix := Two_Real;
      b2 := radix * radix;
      k := 1;
      l := n;

--     .......... search for rows isolating an eigenvalue
--                and push them down ..........
  <<L100>> null;

      for j in reverse 1 .. l loop
         allzero := True;
         for i in 1 .. l loop
            if (i /= j) then
               allzero := (ar(j,i)=Zero_Real) and (ai(j,i)=Zero_Real);
               exit when not allzero;
            end if;
         end loop;
         if allzero then
            m := l;
            exc(m,j);
            if (l = 1) then
               low := k;
               igh := l;
               return;
            end if;
            l := l - 1;
            goto L100;
         end if;
      end loop;

--     .......... search for columns isolating an eigenvalue
--                and push them left ..........
  <<L140>> null;

      for j in k .. l loop
         allzero := True;
         for i in k .. L loop
            if (i /= j) then
               allzero := (ar(i,j)=Zero_Real) and (ai(i,j)=Zero_Real);
               exit when not allzero;
            end if;
         end loop;
         if allzero then
            m := k;
            exc(m,j);
            k := k + 1;
            goto L140;
         end if;
      end loop;

--     .......... now balance the submatrix in rows k to l ..........
      for i in k .. l loop
         scale(i) := One_Real;
      end loop;

--     .......... iterative loop for norm reduction ..........
  <<L190>> noconv := false;

      for i in k .. l loop
         c := Zero_Real;
         r := Zero_Real;
--
         for j in k .. l loop
            if (j /= i) then
               c := c + abs(ar(j,i)) + abs(ai(j,i));
               r := r + abs(ar(i,j)) + abs(ai(i,j));
            end if;
         end loop;

--     .......... guard against zero c or r due to underflow ..........
         if (c /= Zero_Real) and (r /= Zero_Real) then
            g := r / radix;
            f := One_Real;
            s := c + r;
            while (c < g) loop
               f := f * radix;
               c := c * b2;
            end loop;
            g := r * radix;
            while not (c < g) loop
               f := f / radix;
               c := c / b2;
            end loop;

--     .......... now balance ..........
            if (((c + r) / F) < (Near_One * S)) then
               g := One_Real / f;
               scale(i) := scale(i) * f;
               noconv := true;
               for j in k .. n loop
                  ar(i,j) := ar(i,j) * g;
                  ai(i,j) := ai(i,j) * g;
               end loop;
               for j in 1 .. l loop
                  ar(j,i) := ar(j,i) * f;
                  ai(j,i) := ai(j,i) * f;
               end loop;
            end if;
         end if;
      end loop;

      if noconv then goto L190; end if;
      low := k;
      igh := l;
   end cbal;

   procedure cbabk2(n,low,igh,m: in Integer; scale: in Vector; zr,zi: in out Matrix) is
     k: Integer;
     s: Real;

   begin
      if (m = 0) then return; end if;

      if (igh /= low) then
         for i in low .. igh loop
            s := scale(i);
            for j in 1 .. m loop
               zr(i,j) := zr(i,j) * s;
               zi(i,j) := zi(i,j) * s;
            end loop;
         end loop;
      end if;

      for i in reverse 1 .. low-1 loop
         k := integer(Approx(scale(i)));
         if (k /= i) then
            for j in 1 .. m loop
               Swap(zr(i,j),zr(k,j));
               Swap(zi(i,j),zi(k,j));
            end loop;
         end if;
      end loop;

      for i in igh+1 .. n loop
         k := integer(Approx(scale(i)));
         if (k /= i) then
            for j in 1 .. m loop
               Swap(zr(i,j),zr(k,j));
               Swap(zi(i,j),zi(k,j));
            end loop;
         end if;
      end loop;
   end cbabk2;

   procedure comhes(n,low,igh: in Integer; ar,ai: in out Matrix; int: out Pivot) is
-- translated from eispack routine comhes, http://www.netlib.org/eispack/
      i,la,kp1,mm1,mp1: Integer;
      xr,xi,yr,yi: Real;

   begin
      la := igh - 1;
      kp1 := low + 1;
--
      for m in kp1 .. la loop
         mm1 := m - 1;
         xr := Zero_Real;
         xi := Zero_Real;
         i := m;
--
         for j in m .. igh loop
            if ((abs(xr) + abs(xi)) < (abs(ar(j,mm1)) + abs(ai(j,mm1)))) then
               xr := ar(j,mm1);
               xi := ai(j,mm1);
               i := j;
            end if;
         end loop;
--
         int(m) := i;
         if (i /= m) then

--     .......... interchange rows and columns of ar and ai ..........
            for j in mm1 .. n loop
               Swap(ar(i,j),ar(m,j));
               Swap(ai(i,j),ai(m,j));
            end loop;
--
            for j in 1 .. igh loop
               Swap(ar(j,i),ar(j,m));
               Swap(ai(j,i),ai(j,m));
            end loop;
         end if;

--     .......... end interchange ..........
         if (xr /= Zero_Real) or (xi /= Zero_Real) then
            mp1 := m + 1;
            for i in mp1 .. igh loop
               yr := ar(i,mm1);
               yi := ai(i,mm1);
               if (yr /= Zero_Real) or (yi /= Zero_Real) then
                  cdiv(yr,yi,xr,xi,yr,yi);
                  ar(i,mm1) := yr;
                  ai(i,mm1) := yi;
                  for j in m .. n loop
                     ar(i,j) := ar(i,j) - yr * ar(m,j) + yi * ai(m,j);
                     ai(i,j) := ai(i,j) - yr * ai(m,j) - yi * ar(m,j);
                  end loop;
                  for j in 1 .. igh loop
                     ar(j,m) := ar(j,m) + yr * ar(j,i) - yi * ai(j,i);
                     ai(j,m) := ai(j,m) + yr * ai(j,i) + yi * ar(j,i);
                  end loop;
               end if;
            end loop;
         end if;
      end loop;
   end comhes;

   procedure comlr2(n,low,igh: in Integer; int: in Pivot; hr,hi: in out Matrix;
                    wr,wi: in out Vector; zr,zi: in out Matrix; ierr: out Integer) is

-- translated from eispack routine comlr2, http://www.netlib.org/eispack/
      ii,l,m,en,im1,itn,its,mp1,enm1: Integer;
      si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,orm,tst1,tst2: Real;
   begin
      ierr := 0;
--     .......... initialize eigenvector matrix ..........
      for i in  1 .. n loop
         for j in 1 .. n loop
            zr(i,j) := Zero_Real;
            zi(i,j) := Zero_Real;
         end loop;
         zr(i,i) := One_Real;
      end loop;
--     .......... form the matrix of accumulated transformations
--                from the information left by comhes ..........
      for i in reverse low+1 .. igh-1 loop
         for k in i+1 .. igh loop
            zr(k,i) := hr(k,i-1);
            zi(k,i) := hi(k,i-1);
         end loop;
         ii := int(i);
         if (i /= ii) then
            for k in i .. igh loop
               zr(i,k) := zr(ii,k);
               zi(i,k) := zi(ii,k);
               zr(ii,k) := Zero_Real;
               zi(ii,k) := Zero_Real;
            end loop;
            zr(ii,i) := One_Real;
         end if;
      end loop;

--     .......... store roots isolated by cbal ..........
      for i in 1 .. n loop
         if (i < low) or (i > igh) then
            wr(i) := hr(i,i);
            wi(i) := hi(i,i);
         end if;
      end loop;

      en := igh;
      tr := Zero_Real;
      ti := Zero_Real;
      itn := 30*n;

--     .......... search for next eigenvalue ..........
  <<L220>> if (en < low) then goto L680; end if;

      its := 0;
      enm1 := en - 1;

--     .......... look for single small sub-diagonal element
      loop
         for ll in reverse low .. en loop
            l := ll;
            exit when (l = low);
            tst1 := abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) + abs(hr(l,l)) + abs(hi(l,l));
            tst2 := tst1 + abs(hr(l,l-1)) + abs(hi(l,l-1));
            exit when (tst2 = tst1);
         end loop;

         if (l = en) then goto L660; end if;

         if (itn = 0) then
--     .......... set error -- all eigenvalues have not
--                converged after 30*n iterations ..........
           ierr := en;
           return;
         end if;

--     .......... form shift ..........
         if (its /= 10) and (its /= 20) then
            sr := hr(en,en);
            si := hi(en,en);                                              --- can get TINY
            xr := hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1);
            xi := hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1);
            if (xr /= Zero_Real) or (xi /= Zero_Real) then
               yr := (hr(enm1,enm1) - sr) / Two_Real;
               yi := (hi(enm1,enm1) - si) / Two_Real;
               csroot(yr*yr-yi*yi+xr,Two_Real*yr*yi+xi,zzr,zzi);
               if ((yr * zzr + yi * zzi) < Zero_Real) then
                  zzr := -zzr;
                  zzi := -zzi;
               end if;

               cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi);
               sr := sr - xr;
               si := si - xi;
            end if;
         else
--     .......... form exceptional shift ..........
            sr := abs(hr(en,enm1)) + abs(hr(enm1,en-2));
            si := abs(hi(en,enm1)) + abs(hi(enm1,en-2));
         end if;

         for i in low .. en loop
            hr(i,i) := hr(i,i) - sr;
            hi(i,i) := hi(i,i) - si;
         end loop;

         tr := tr + sr;
         ti := ti + si;
         its := its + 1;
         itn := itn - 1;

--     .......... look for two consecutive small
--                sub-diagonal elements ..........
         xr := abs(hr(enm1,enm1)) + abs(hi(enm1,enm1));
         yr := abs(hr(en,enm1)) + abs(hi(en,enm1));
         zzr := abs(hr(en,en)) + abs(hi(en,en));

         for mm in reverse l .. enm1 loop
            m := mm;
            exit when (m = l);
            yi := yr;
            yr := abs(hr(m,m-1)) + abs(hi(m,m-1));
            xi := zzr;
            zzr := xr;
            xr := abs(hr(m-1,m-1)) + abs(hi(m-1,m-1));
            tst1 := (zzr / yi) * (zzr + xr + xi);
            tst2 := tst1 + yr;
            exit when (tst2 = tst1);
         end loop;

--     .......... triangular decomposition h:=l*r ..........
         mp1 := m + 1;

         for i in mp1 .. en loop
            im1 := i - 1;
            xr := hr(im1,im1);
            xi := hi(im1,im1);
            yr := hr(i,im1);
            yi := hi(i,im1);

            if ((abs(xr) + abs(xi)) < (abs(yr) + abs(yi))) then
--     .......... interchange rows of hr and hi ..........
               for j in im1 .. n loop
                  Swap(hr(im1,j),hr(i,j));
                  Swap(hi(im1,j),hi(i,j));
               end loop;

               cdiv(xr,xi,yr,yi,zzr,zzi);
               wr(i) := One_Real;
            else
               cdiv(yr,yi,xr,xi,zzr,zzi);
               wr(i) := -One_Real;
            end if;

            hr(i,im1) := zzr;
            hi(i,im1) := zzi;

            for j in i .. n loop
               hr(i,j) := hr(i,j) - zzr * hr(im1,j) + zzi * hi(im1,j);
               hi(i,j) := hi(i,j) - zzr * hi(im1,j) - zzi * hr(im1,j);
            end loop;
         end loop;

--     .......... composition r*l:=h ..........
         for j in mp1 .. en loop
            xr := hr(j,j-1);
            xi := hi(j,j-1);
            hr(j,j-1) := Zero_Real;
            hi(j,j-1) := Zero_Real;

--     .......... interchange columns of hr, hi, zr, and zi, if necessary ..........
            if (Zero_Real < wr(j)) then
               for i in 1 .. j loop
                  Swap(hr(i,j-1),hr(i,j));
                  Swap(hi(i,j-1),hi(i,j));
               end loop;

               for i in low .. igh loop
                  Swap(zr(i,j-1),zr(i,j));
                  Swap(zi(i,j-1),zi(i,j));
               end loop;
            end if;

            for i in 1 .. j loop
               hr(i,j-1) := hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j);
               hi(i,j-1) := hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j);
            end loop;

--     .......... accumulate transformations ..........
            for i in low .. igh loop
               zr(i,j-1) := zr(i,j-1) + xr * zr(i,j) - xi * zi(i,j);
               zi(i,j-1) := zi(i,j-1) + xr * zi(i,j) + xi * zr(i,j);
            end loop;
         end loop;
      end loop;

--     .......... a root found ..........
  <<L660>> null;
      hr(en,en) := hr(en,en) + tr;
      wr(en) := hr(en,en);
      hi(en,en) := hi(en,en) + ti;
      wi(en) := hi(en,en);
      en := enm1;
      goto L220;

--     .......... all roots found.  backsubstitute to find
--                vectors of upper triangular form ..........
  <<L680>> null;
      orm := Zero_Real;
      for i in 1 .. n loop
         for j in i .. n loop
            tr := abs(hr(i,j)) + abs(hi(i,j));
            if (orm < tr) then orm := tr; end if;
         end loop;
      end loop;

      hr(1,1) := orm;
      if (n = 1) or (orm = Zero_Real) then return; end if;

      for nn in reverse 2 .. n loop
         en := nn;
         xr := wr(en);
         xi := wi(en);
         hr(en,en) := One_Real;
         hi(en,en) := Zero_Real;
         enm1 := en - 1;

         for i in reverse 1 .. enm1 loop
            zzr := Zero_Real;
            zzi := Zero_Real;

            for j in i+1 .. en loop
               zzr := zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en);
               zzi := zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en);
            end loop;

            yr := xr - wr(i);
            yi := xi - wi(i);
            if (yr = Zero_Real) and (yi = Zero_Real) then
               tst1 := orm;
               yr := tst1;
               loop
                  yr := Small * yr;
                  tst2 := orm + yr;
                  exit when not (tst1 < tst2);
               end loop;
            end if;
            cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en));

--     .......... overflow control ..........
            tr := abs(hr(i,en)) + abs(hi(i,en));
            if (tr /= Zero_Real) then
               tst1 := tr;
               tst2 := tst1 + One_Real/tst1;
               if not (tst1 < tst2) then
                  for j in i .. en loop
                     hr(j,en) := hr(j,en)/tr;
                     hi(j,en) := hi(j,en)/tr;
                  end loop;
               end if;
            end if;

         end loop;
      end loop;

--     .......... end backsubstitution ..........
--     .......... vectors of isolated roots ..........
      for i in 1 .. n loop
         if (i < low) or (i > igh) then
            for j in i .. n loop
               zr(i,j) := hr(i,j);
               zi(i,j) := hi(i,j);
            end loop;
         end if;
      end loop;

--     .......... multiply by transformation matrix to give
--                vectors of original full matrix.
      for j in reverse low .. n loop
         m := Integer'Min(j,igh);
--
         for i in low .. igh loop
            zzr := Zero_Real;
            zzi := Zero_Real;
            for k in low .. m loop
               zzr := zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j);
               zzi := zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j);
            end loop;
            zr(i,j) := zzr;
            zi(i,j) := zzi;
         end loop;
      end loop;
   end comlr2;

end Eispack2;
