
package body EisPack1 is

   function dsign(a,b: Real) return Real is
---  ancient fortran routine
   begin
      if (b<Zero_Real) then return -abs(a); else return abs(a); end if;
   end dsign;

   procedure tred1(n: in Integer; a: in out Matrix; d,e,e2: in out Vector) is
-- translated from eispack routine tred1, http://www.netlib.org/eispack/
      l: Integer;
      f,g,h,scale: Real;

begin

      for i in 1 .. n loop
         d(i) := a(n,i);
         a(n,i) := a(i,i);
      end loop;

      for i in reverse 1 .. n loop
         l := i - 1;
         h := Zero_Real;
         scale := Zero_Real;

         if (l < 1) then
            e(i) := Zero_Real;
            e2(i) := Zero_Real;
            goto L300;
         else

--     .......... scale row (algol tol then not needed) ..........

            for k in 1 .. l loop
               scale := scale + abs(d(k));
            end loop;
            if (scale = Zero_Real) then
               for j in 1 .. l loop
                  d(j) := a(l,j);
                  a(l,j) := a(i,j);
                  a(i,j) := Zero_Real;
               end loop;
               e(i) := Zero_Real;
               e2(i) := Zero_Real;
               goto L300;
            end if;
         end if;

         for k in 1 .. l loop
            d(k) := d(k) / scale;
            h := h + d(k) * d(k);
         end loop;
--
         e2(i) := scale * scale * h;
         f := d(l);
         g := -dsign(Sqrt(h),f);
         e(i) := scale * g;
         h := h - f * g;
         d(l) := f - g;
         if (l /= 1) then

--     .......... form a*u ..........

            for j in 1 .. l loop
               e(j) := Zero_Real;
            end loop;

            for j in 1 .. l loop
               f := d(j);
               g := e(j) + a(j,j) * f;
               if (l >= j+1) then
                  for k in j+1 .. l loop
                     g := g + a(k,j) * d(k);
                     e(k) := e(k) + a(k,j) * f;
                  end loop;
               end if;
               e(j) := g;
            end loop;

--     .......... form p ..........

            f := Zero_Real;
            for j in 1 .. l loop
               e(j) := e(j) / h;
               f := f + e(j) * d(j);
            end loop;
            h := f / (h + h);

--     .......... form q ..........

            for j in 1 .. l loop
               e(j) := e(j) - h * d(j);
            end loop;

--     .......... form reduced a ..........

            for j in 1 .. l loop
               f := d(j);
               g := e(j);
               for k in j .. l loop
                  a(k,j) := a(k,j) - f * e(k) - g * d(k);
               end loop;
            end loop;
         end if;

         for j in 1 .. l loop
            f := d(j);
            d(j) := a(l,j);
            a(l,j) := a(i,j);
            a(i,j) := f * scale;
         end loop;

  <<L300>> null;

      end loop;
   end tred1;

   procedure tred2(n: in Integer; a: in Matrix; d,e: in out Vector; z: in out Matrix) is
-- translated from eispack routine tql1, http://www.netlib.org/eispack/
      l: integer;
      f,g,h,hh,scale: Real;

   begin
      if (n = 1) then
         d(1) := a(1,1);
         z(1,1) := One_Real;
         e(1) := Zero_Real;
         return;
      end if;

      for i in 1 .. n loop
         for j in i .. n loop
            z(j,i) := a(j,i);
         end loop;
         d(i) := a(n,i);
      end loop;

      for i in reverse 2 .. n loop
         l := i - 1;
         h := Zero_Real;
         scale := Zero_Real;

         if (l < 2) then
            e(i) := d(l);
            for j in 1 .. l loop
               d(j) := z(l,j);
               z(i,j) := Zero_Real;
               z(j,i) := Zero_Real;
            end loop;
            goto L290;
         else

--     .......... scale row (algol tol then not needed) ..........

            for k in 1 .. l loop
               scale := scale + abs(d(k));
            end loop;
            if (scale = Zero_Real) then
               e(i) := d(l);
               for j in 1 .. l loop
                  d(j) := z(l,j);
                  z(i,j) := Zero_Real;
                  z(j,i) := Zero_Real;
               end loop;
               goto L290;
            end if;
         end if;

         for k in 1 .. l loop
            d(k) := d(k) / scale;
            h := h + d(k) * d(k);
         end loop;

         f := d(l);
         g := -dsign(Sqrt(h),f);
         e(i) := scale * g;
         h := h - f * g;
         d(l) := f - g;

--     .......... form a*u ..........

         for j in 1 .. l loop
            e(j) := Zero_Real;
         end loop;

         for j in 1 .. l loop
            f := d(j);
            z(j,i) := f;
            g := e(j) + z(j,j) * f;
            if (l > j) then
               for k in j+1 .. l loop
                  g := g + z(k,j) * d(k);
                  e(k) := e(k) + z(k,j) * f;
               end loop;
            end if;
            e(j) := g;
         end loop;

--     .......... form p ..........

         f := Zero_Real;

         for j in 1 .. l loop
            e(j) := e(j) / h;
            f := f + e(j) * d(j);
         end loop;
         hh := f / (h + h);

--     .......... form q ..........

         for j in 1 .. l loop
            e(j) := e(j) - hh * d(j);
         end loop;

--     .......... form reduced a ..........

         for j in 1 .. l loop
            f := d(j);
            g := e(j);
            for k in j .. l loop
               z(k,j) := z(k,j) - f * e(k) - g * d(k);
            end loop;
            d(j) := z(l,j);
            z(i,j) := Zero_Real;
         end loop;

  <<L290>> null;

         d(i) := h;
      end loop;

--     .......... accumulation of transformation matrices ..........

      for i in 2 .. n loop
         l := i - 1;
         z(n,l) := z(l,l);
         z(l,l) := One_Real;
         h := d(i);
          if (h /= Zero_Real) then
            for k in 1 .. l loop
               d(k) := z(k,i) / h;
            end loop;

            for j in 1 .. l loop
               g := Zero_Real;
               for k in 1 .. l loop
                  g := g + z(k,i) * z(k,j);
               end loop;
               for k in 1 .. l loop
                  z(k,j) := z(k,j) - g * d(k);
               end loop;
            end loop;
         end if;
         for k in 1 .. l loop
            z(k,i) := Zero_Real;
         end loop;
      end loop;

      for i in 1 .. n loop
         d(i) := z(n,i);
         z(n,i) := Zero_Real;
      end loop;
      z(n,n) := One_Real;
      e(1) := Zero_Real;
   end tred2;

   procedure tql1(n: in Integer; d,e: in out Vector; ierr: out Integer) is
-- translated from eispack routine tql1, http://www.netlib.org/eispack/
     j,m,l1,l2: integer;
     c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2: Real;

   begin
      ierr := 0;
      if (n = 1) then return; end if;
      for i in 2 .. n loop
          e(i-1) := e(i);
      end loop;
      f := Zero_Real;
      tst1 := Zero_Real;
      e(n) := Zero_Real;

      for l in 1 .. n loop
         j := 0;
         h := abs(d(l)) + abs(e(l));
         if (tst1 < h) then tst1 := h; end if;

--     .......... look for small sub-diagonal element ..........

         m := l;
         loop
            tst2 := tst1 + abs(e(m));
            exit when (tst2 = tst1);
            m := m+1;
         end loop;

         if (m /= l) then
            loop
               if (j = 30) then

--     .......... set error -- no convergence to an
--                eigenvalue after 30 iterations ..........

                  ierr := l;
                  return;
               end if;

               j := j + 1;

--     .......... form shift ..........

               l1 := l + 1;
               l2 := l1 + 1;
               g := d(l);
               p := (d(l1) - g) / (Two_Real * e(l));
               r := Sqrt(One_Real+p*p);
               d(l) := e(l) / (p + dsign(r,p));
               d(l1) := e(l) * (p + dsign(r,p));
               dl1 := d(l1);
               h := g - d(l);
               for i in l2 .. n loop
                  d(i) := d(i) - h;
               end loop;
               f := f + h;

--     .......... ql transformation ..........

               p := d(m);
               c := One_Real;
               c2 := c;
               el1 := e(l1);
               s := Zero_Real;
               for i in reverse l .. m-1 loop
                  c3 := c2;
                  c2 := c;
                  s2 := s;
                  g := c * e(i);
                  h := c * p;
                  r := Sqrt(p*p+e(i)*e(i));
                  e(i+1) := s * r;
                  s := e(i) / r;
                  c := p / r;
                  p := c * d(i) - s * g;
                  d(i+1) := h + s * (c * g + s * d(i));
               end loop;
--
               p := -s * s2 * c3 * el1 * e(l) / dl1;
               e(l) := s * p;
               d(l) := c * p;
               tst2 := tst1 + abs(e(l));
               exit when (tst2 = tst1);
            end loop;
         end if;

         p := d(l) + f;

--     .......... order eigenvalues ..........

         if (l = 1) then d(1) := p; end if;

         for i in reverse 2 .. l loop
            if not (p < d(i-1)) then
               d(i) := p;
               exit;
             end if;
            d(i) := d(i-1);
            if (i = 2) then d(1) := p; end if;
         end loop;
      end loop;
   end tql1;

   procedure tql2(n: in Integer; d,e: in out Vector; z: in out Matrix; ierr: out Integer) is
-- translated from eispack routine tql2, http://www.netlib.org/eispack/
      j,m,l1,l2: Integer;
      c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2: Real;

   begin
      ierr := 0;
      if (n = 1) then return; end if;
      for i in 2 .. n loop
         e(i-1) := e(i);
      end loop;
      f := Zero_Real;
      tst1 := Zero_Real;
      e(n) := Zero_Real;

      for l in 1 .. n loop
         j := 0;
         h := abs(d(l)) + abs(e(l));
         if (tst1 < h) then tst1 := h; end if;

--     .......... look for small sub-diagonal element ..........

         m := l;
         loop
            tst2 := tst1 + abs(e(m));
            exit when (tst2 = tst1);
            m := m+1;
         end loop;

         if (m /= l) then
            loop
               if (j = 30) then

--     .......... set error -- no convergence to an
--                eigenvalue after 30 iterations ..........

                  ierr := l;
                  return;
               end if;

               j := j + 1;

--     .......... form shift ..........

               l1 := l + 1;
               l2 := l1 + 1;
               g := d(l);
               p := (d(l1) - g) / (Two_Real * e(l));
               r := Sqrt(One_Real+p*p);
               d(l) := e(l) / (p + dsign(r,p));
               d(l1) := e(l) * (p + dsign(r,p));
               dl1 := d(l1);
               h := g - d(l);
               for i in l2 .. n loop
                  d(i) := d(i) - h;
               end loop;
               f := f + h;

--     .......... ql transformation ..........

               p := d(m);
               c := One_Real;
               c2 := c;
               el1 := e(l1);
               s := Zero_Real;
               for i in reverse l .. m-1 loop
                  c3 := c2;
                  c2 := c;
                  s2 := s;
                  g := c * e(i);
                  h := c * p;
                  r := Sqrt(p*p+e(i)*e(i));
                  e(i+1) := s * r;
                  s := e(i) / r;
                  c := p / r;
                  p := c * d(i) - s * g;
                  d(i+1) := h + s * (c * g + s * d(i));

--     .......... form vector ..........
                  for k in 1 .. n loop
                     h := z(k,i+1);
                     z(k,i+1) := s * z(k,i) + c * h;
                     z(k,i) := c * z(k,i) - s * h;
                  end loop;
               end loop;

               p := -s * s2 * c3 * el1 * e(l) / dl1;
               e(l) := s * p;
               d(l) := c * p;
               tst2 := tst1 + abs(e(l));
               exit when (tst2 = tst1);
            end loop;
         end if;

         d(l) := d(l) + f;
      end loop;

--     .......... order eigenvalues and eigenvectors ..........

      for i in 1 .. n-1 loop
         m := i;
         p := d(i);
         for l in i+1 .. n loop
            if (d(l) < p) then
               m := l;
               p := d(l);
            end if;
         end loop;
         if (m /= i) then
            d(m) := d(i);
            d(i) := p;
            for l in 1 .. n loop
               p := z(l,i);
               z(l,i) := z(l,m);
               z(l,m) := p;
            end loop;
         end if;
      end loop;
   end tql2;

end EisPack1;
