File Coverage

lib/Math/MatrixDecomposition/Eigen.pm

Criterion	Covered	Total	%
statement	530	649	81.6
branch	104	220	47.2
condition	54	120	45.0
subroutine	17	19	89.4
pod	7	8	87.5
total	712	1016	70.0

## Math/MatrixDecomposition/Eigen.pm --- eigenvalues and eigenvectors.

# Copyright (C) 2010 Ralph Schleicher.  All rights reserved.

# This program is free software; you can redistribute it and/or
# modify it under the same terms as Perl itself.

# Commentary:
#
# Code derived from EiSPACK procedures and
# Jama's 'EigenvalueDecomposition' class.

## Code:

package Math::MatrixDecomposition::Eigen;

use strict;
use warnings;
use Carp;
use Exporter qw(import);
use POSIX qw(:float_h);
use Math::Complex qw();
use Scalar::Util qw(looks_like_number);

use Math::BLAS;
use Math::MatrixDecomposition::Util qw(:all);

BEGIN
{
  our $VERSION = '1.06';
  our @EXPORT_OK = qw(eig);
}

# Calculate eigenvalues and eigenvectors (convenience function).
sub eig
{
  __PACKAGE__->new (@_);
}

# Standard constructor.
sub new
{
  my $class = shift;
  my $self =
    {
     # Eigenvalues (a vector).
     value => undef,

     # Eigenvectors (an array of vectors).
     vector => undef,
    };

  bless ($self, ref ($class) || $class);

  # Process arguments.
  $self->decompose (@_)
    if @_ > 0;

  # Return object.
  $self;
}

# Calculate eigenvalues and eigenvectors.
sub decompose
{
  my $self = shift;

  # Check arguments.
  my $a = shift;

  my $m = @_ > 0 && looks_like_number ($_[0]) ? shift : sqrt (@$a);
  my $n = @_ > 0 && looks_like_number ($_[0]) ? shift : $m ? @$a / $m : 0;

  croak ('Invalid argument')
    if (@$a != ($m * $n)
	|| mod ($m, 1) != 0 || $m < 1
	|| mod ($n, 1) != 0 || $n < 1);

  croak ('Matrix has to be square')
    if $m != $n;

  # Get properties.
  my %prop = @_;

  $prop{balance} //= 1;
  $prop{normalize} //= 1;
  $prop{positive} //= 1;

  # Index of last row/column.
  my $end = $n - 1;

  # Eigenvalues (a vector).
  $$self{value} //= [];

  # Vector $d contains the real part of the eigenvalues and vector $e
  # contains the imaginary part of the eigenvalues.
  my $d = $$self{value};
  my $e = [];

  splice (@$d, $n)
    if @$d > $n;

  # Eigenvectors (an array of column vectors).
  $$self{vector} //= [];

  # Matrix $Z contains the eigenvectors, please note
  # that $$Z[$j][$i] denotes matrix element Z(i,j).
  my $Z = $$self{vector};

  splice (@$Z, $n)
    if @$Z > $n;

  for my $v (@$Z)
    {
      splice (@$v, $n)
	if @$v > $n;
    }

  # True if matrix is symmetric.
  my $sym = 1;

 SYM:

  for my $i (0 .. $end)
    {
      for my $j ($i + 1 .. $end)
	{
	  if ($$a[$i * $n + $j] != $$a[$j * $n + $i])
	    {
	      $sym = 0;
	      last SYM;
	    }
	}
    }

  if ($sym)
    {
      # Copy matrix elements.
      for my $j (0 .. $end)
	{
	  $$Z[$j] //= [];

	  for my $i (0 .. $end)
	    {
	      $$Z[$j][$i] = $$a[$i * $n + $j];
	    }
	}

      # Reduce a real symmetric matrix to a symmetric tridiagonal matrix
      # using and accumulating orthogonal similarity transformations.
      #
      # See EiSPACK procedure 'tred2'.
      if (1)
	{
	  my ($i, $j, $k, $l,
	      $f, $g, $h, $hh, $scale);

	  for $i (0 .. $end)
	    {
	      $$d[$i] = $$Z[$i][$end];
	    }

	  for $i (reverse (1 .. $end))
	    {
	      $l = $i - 1;
	      $h = 0;

	      # Scale row.
	      $scale = 0;

	      for $k (0 .. $l)
		{
		  $scale += abs ($$d[$k]);
		}

	      if ($scale == 0)
		{
		  $$e[$i] = $$d[$l];

		  for $j (0 .. $l)
		    {
		      $$d[$j] = $$Z[$j][$l];

		      $$Z[$j][$i] = 0;
		      $$Z[$i][$j] = 0;
		    }
		}
	      else
		{
		  for $k (0 .. $l)
		    {
		      $$d[$k] /= $scale;
		      $h += $$d[$k] ** 2;
		    }

		  $f = $$d[$l];
		  $g = - sign (sqrt ($h), $f);
		  $h -= $f * $g;

		  $$d[$l] = $f - $g;
		  $$e[$i] = $scale * $g;

		  # Form a*u.
		  for $j (0 .. $l)
		    {
		      $$e[$j] = 0;
		    }

		  for $j (0 .. $l)
		    {
		      $f = $$d[$j];
		      $$Z[$i][$j] = $f;
		      $g = $$e[$j] + $$Z[$j][$j] * $f;

		      for $k ($j + 1 .. $l)
			{
			  $g += $$Z[$j][$k] * $$d[$k];
			  $$e[$k] += $$Z[$j][$k] * $f;
			}

		      $$e[$j] = $g;
		    }

		  # Form p.
		  $f = 0;

		  for $j (0 .. $l)
		    {
		      $$e[$j] /= $h;
		      $f += $$e[$j] * $$d[$j];
		    }

		  # Form q.
		  $hh = $f / ($h + $h);

		  for $j (0 .. $l)
		    {
		      $$e[$j] -= $hh * $$d[$j];
		    }

		  # Form reduced a.
		  for $j (0 .. $l)
		    {
		      $f = $$d[$j];
		      $g = $$e[$j];

		      for $k ($j .. $l)
			{
			  $$Z[$j][$k] -= ($f * $$e[$k] + $g * $$d[$k]);
			}

		      $$d[$j] = $$Z[$j][$l];
		      $$Z[$j][$i] = 0;
		    }
		}

	      $$d[$i] = $h;
	    }

	  # Accumulation of transformation matrices.
	  for $i (1 .. $end)
	    {
	      $l = $i - 1;

	      $$Z[$l][$end] = $$Z[$l][$l];
	      $$Z[$l][$l] = 1;

	      $h = $$d[$i];
	      if ($h != 0)
		{
		  for $k (0 .. $l)
		    {
		      $$d[$k] = $$Z[$i][$k] / $h;
		    }

		  for $j (0 .. $l)
		    {
		      $g = 0;

		      for $k (0 .. $l)
			{
			  $g += $$Z[$i][$k] * $$Z[$j][$k];
			}

		      for $k (0 .. $l)
			{
			  $$Z[$j][$k] -= $g * $$d[$k];
			}
		    }
		}

	      for $k (0 .. $l)
		{
		  $$Z[$i][$k] = 0;
		}
	    }

	  for $j (0 .. $end)
	    {
	      $$d[$j] = $$Z[$j][$end];
	      $$Z[$j][$end] = 0;
	    }

	  $$Z[$end][$end] = 1;
	  $$e[0] = 0;
	}

      # Find the eigenvalues and eigenvectors of a symmetric tridiagonal
      # matrix by the QL method.
      #
      # See EiSPACK procedure 'tql2'.
      if (1)
	{
	  my ($i, $j, $k, $l, $m,
	      $c, $c2, $c3, $dl1, $el1, $f, $g, $h, $p, $r, $s, $s2, $t, $t2);

	  for $i (1 .. $end)
	    {
	      $$e[$i - 1] = $$e[$i];
	    }

	  $$e[$end] = 0;

	  $f = 0;
	  $t = 0;

	  for $l (0 .. $end)
	    {
	      $t2 = abs ($$d[$l]) + abs ($$e[$l]);
	      $t = $t2 if $t2 > $t;

              for ($m = $l; $m < $n; ++$m)
		{
		  last if abs ($$e[$m]) <= eps * $t;
		}

	      if ($m > $l)
		{
		  while (1)
		    {
		      $g = $$d[$l];
		      $p = ($$d[$l + 1] - $g) / (2 * $$e[$l]);
		      $r = sign (hypot ($p, 1), $p);

		      $$d[$l] = $$e[$l] / ($p + $r);
		      $$d[$l + 1] = $$e[$l] * ($p + $r);
		      $dl1 = $$d[$l + 1];
		      $h = $g - $$d[$l];

		      for $i ($l + 2 .. $end)
			{
			  $$d[$i] -= $h;
			}

		      $f += $h;

		      $p = $$d[$m];
		      $c = 1;
		      $c2 = $c;
		      $el1 = $$e[$l + 1];
		      $s = 0;

		      for $i (reverse ($l .. $m - 1))
			{
			  $c3 = $c2;
			  $c2 = $c;
			  $s2 = $s;
			  $g = $c * $$e[$i];
			  $h = $c * $p;
			  $r = hypot ($p, $$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]);

			  for $k (0 .. $end)
			    {
			      $h = $$Z[$i + 1][$k];
			      $$Z[$i + 1][$k] = $s * $$Z[$i][$k] + $c * $h;
			      $$Z[$i][$k] = $c * $$Z[$i][$k] - $s * $h;
			    }
			}

		      $p = 0 - $s * $s2 * $c3 * $el1 * $$e[$l] / $dl1;

		      $$e[$l] = $s * $p;
		      $$d[$l] = $c * $p;

		      # Check convergence.
		      last if abs ($$e[$l]) <= eps * $t;
		    }
		}

	      $$d[$l] = $$d[$l] + $f;
	      $$e[$l] = 0;
	    }
	}
    }
  else
    {
      # Hessenberg matrix.
      my @H = @$a;
      my $N = $n;

      # Row and column indices of the beginning and end
      # of the principal sub-matrix.
      my $lo = 0;
      my $hi = $end;

      # Permutation vector.
      my @perm = ();

      # Scaling vector.
      my @scale = ();

      # Balance a real matrix and isolate eigenvalues whenever possible.
      if ($prop{balance})
	{
	  my ($p, $s) = gebal ($N, \@H, low => \$lo, high => \$hi);

	  @perm = @$p;
	  @scale = @$s;
	}

      # Reduce matrix to upper Hessenberg form by orthogonal similarity
      # transformations.
      #
      # See EiSPACK procedure 'orthes'.
      my @ort = ();

      if (1)
	{
	  my ($i, $j, $k, $m,
	      $f, $g, $h, $scale);

	  for $m ($lo + 1 .. $hi - 1)
	    {
	      $scale = 0;

	      for $i ($m .. $hi)
		{
		  $scale += abs ($H[$i * $N + ($m - 1)]);
		}

	      next if $scale == 0;

	      $h = 0;

	      for $i (reverse ($m .. $hi))
		{
		  $ort[$i] = $H[$i * $N + ($m - 1)] / $scale;
		  $h += $ort[$i] ** 2;
		}

	      $g = - sign (sqrt ($h), $ort[$m]);
	      $h -= $ort[$m] * $g;
	      $ort[$m] -= $g;

	      for $j ($m .. $end)
		{
		  $f = 0;

		  for $i (reverse ($m .. $hi))
		    {
		      $f += $ort[$i] * $H[$i * $N + $j];
		    }

		  $f /= $h;

		  for $i ($m .. $hi)
		    {
		      $H[$i * $N + $j] -= $f * $ort[$i];
		    }
		}

	      for $i (0 .. $hi)
		{
		  $f = 0;

		  for $j (reverse ($m .. $hi))
		    {
		      $f += $ort[$j] * $H[$i * $N + $j];
		    }

		  $f /= $h;

		  for $j ($m .. $hi)
		    {
		      $H[$i * $N + $j] -= $f * $ort[$j];
		    }
		}

	      $ort[$m] *= $scale;
	      $H[$m * $N + ($m - 1)] = $scale * $g;
	    }
	}

      # Accumulate the orthogonal similarity transformations.
      #
      # See EiSPACK procedure 'ortran'.
      if (1)
	{
	  my ($i, $j, $k, $m,
	      $g);

	  for $j (0 .. $end)
	    {
	      $$Z[$j] //= [];

	      for $i (0 .. $end)
		{
		  $$Z[$j][$i] = ($i == $j ? 1 : 0);
		}
	    }

	  for $m (reverse ($lo + 1 .. $hi - 1))
	    {
	      next if $H[$m * $N + ($m - 1)] == 0;

	      for $i ($m + 1 .. $hi)
		{
		  $ort[$i] = $H[$i * $N + ($m - 1)];
		}

	      for $j ($m .. $hi)
		{
		  $g = 0;

		  for $i ($m .. $hi)
		    {
		      $g += $ort[$i] * $$Z[$j][$i];
		    }

		  $g = ($g / $ort[$m]) / $H[$m * $N + ($m - 1)];

		  for $i ($m .. $hi)
		    {
		      $$Z[$j][$i] += $g * $ort[$i];
		    }
		}
	    }
	}

      # Find the eigenvalues and eigenvectors of a real upper Hessenberg
      # matrix by the QR method.
      #
      # See EiSPACK procedure 'hqr2'.
      if (1)
	{
	  my ($i, $j, $k, $l, $m,
	      $h, $g, $f, $p, $q, $r, $s, $t, $w, $x, $y, $z,
	      $norm, $iter, $not_last, $vr, $vi, $ra, $sa);

	  $t = 0;

	  # Store isolated roots.
	  for $i (0 .. $end)
	    {
	      if ($i < $lo || $i > $hi)
		{
		  $$d[$i] = $H[$i * $N + $i];
		  $$e[$i] = 0;
		}
	    }

	  # Compute matrix norm.
	  $norm = 0;

	  for $i (0 .. $end)
	    {
	      for $j ($i .. $end)
		{
		  $norm += abs ($H[$i * $N + $j]);
		}
	    }

	  # Search for next eigenvalue.
	  $iter = 0;

	  for ($n = $end; $n >= $lo; )
	    {
	      # Look for single small sub-diagonal element.
	      for ($l = $n; $l > $lo; --$l)
		{
		  $s = abs ($H[($l - 1) * $N + ($l - 1)]) + abs ($H[$l * $N + $l]);
		  $s = $norm
		    if $s == 0;

		  last if abs ($H[$l * $N + ($l - 1)]) < eps * $s;
		}

	      $x = $H[$n * $N + $n];

	      if ($l == $n)
		{
		  # One root found,
		  $H[$n * $N + $n] = $x + $t;

		  $$d[$n] = $H[$n * $N + $n];
		  $$e[$n] = 0;

		  $n -= 1;
		  $iter = 0;
		}
	      elsif ($l == $n - 1)
		{
		  # Two roots found.
		  $y = $H[($n - 1) * $N + ($n - 1)];
		  $w = $H[$n * $N + ($n - 1)] * $H[($n - 1) * $N + $n];

		  $p = ($y - $x) / 2;
		  $q = $p * $p + $w;
		  $z = sqrt (abs ($q));

		  $H[$n * $N + $n] = $x + $t;
		  $H[($n - 1) * $N + ($n - 1)] = $y + $t;
		  $x = $H[$n * $N + $n];

		  if ($q >= 0)
		    {
		      # Real pair.
		      $z = $p + sign ($z, $p);

		      $$d[$n - 1] = $x + $z;
		      $$d[$n] = $$d[$n - 1];
		      $$d[$n] = $x - $w / $z
			if $z != 0;

		      $$e[$n - 1] = 0;
		      $$e[$n] = 0;

		      $x = $H[$n * $N + ($n - 1)];
		      $s = abs ($x) + abs ($z);
		      $p = $x / $s;
		      $q = $z / $s;
		      $r = sqrt ($p * $p + $q * $q);
		      $p = $p / $r;
		      $q = $q / $r;

		      # Row modification.
		      for $j ($n - 1 .. $end)
			{
			  $z = $H[($n - 1) * $N + $j];
			  $H[($n - 1) * $N + $j] = $q * $z + $p * $H[$n * $N + $j];
			  $H[$n * $N + $j] = $q * $H[$n * $N + $j] - $p * $z;
			}

		      # Column modification.
		      for $i (0 .. $n)
			{
			  $z = $H[$i * $N + ($n - 1)];
			  $H[$i * $N + ($n - 1)] = $q * $z + $p * $H[$i * $N + $n];
			  $H[$i * $N + $n] = $q * $H[$i * $N + $n] - $p * $z;
			}

		      # Accumulate transformations.
		      for $i ($lo .. $hi)
			{
			  $z = $$Z[$n - 1][$i];
			  $$Z[$n - 1][$i] = $q * $z + $p * $$Z[$n][$i];
			  $$Z[$n][$i] = $q * $$Z[$n][$i] - $p * $z;
			}
		    }
		  else
		    {
		      # Complex pair.
		      $$d[$n - 1] = $x + $p;
		      $$d[$n] = $x + $p;
		      $$e[$n - 1] = $z;
		      $$e[$n] = 0 - $z;
		    }

		  $n -= 2;
		  $iter = 0;
		}
	      else
		{
		  # Form shift.
		  $y = $H[($n - 1) * $N + ($n - 1)];
		  $w = $H[$n * $N + ($n - 1)] * $H[($n - 1) * $N + $n];

		  # Wilkinson's original ad hoc shift.
		  if ($iter == 10 || $iter == 20)
		    {
		      $t += $x;

		      for $i ($lo .. $n)
			{
			  $H[$i * $N + $i] -= $x;
			}

		      $s = abs ($H[$n * $N + ($n - 1)]) + abs ($H[($n - 1) * $N + ($n - 2)]);
		      $x = 0.75 * $s;
		      $y = $x;
		      $w = -0.4375 * $s * $s;
		    }

		  # Matlab's new ad hoc shift.
		  if ($iter == 30)
		    {
		      $s = ($y - $x) / 2;
		      $s = $s * $s + $w;
		      if ($s > 0)
			{
			  $s = sqrt ($s);
			  $s = - $s if $y < $x;
			  $s = $x - $w / (($y - $x) / 2 + $s);

			  for $i ($lo .. $n)
			    {
			      $H[$i * $N + $i] -= $s;
			    }

			  $t += $s;
			  $x = 0.964;
			  $w = $y = $x;
			}
		    }

		  ++$iter;

		  # Look for two consecutive small sub-diagonal elements.
		  for ($m = $n - 2; $m >= $l; --$m)
		    {
		      $z = $H[$m * $N + $m];
		      $r = $x - $z;
		      $s = $y - $z;
		      $p = ($r * $s - $w) / $H[($m + 1) * $N + $m] + $H[$m * $N + ($m + 1)];
		      $q = $H[($m + 1) * $N + ($m + 1)] - $z - $r - $s;
		      $r = $H[($m + 2) * $N + ($m + 1)];
		      $s = abs ($p) + abs ($q) + abs ($r);
		      $p = $p / $s;
		      $q = $q / $s;
		      $r = $r / $s;

		      last if $m == $l;
		      last if abs ($H[$m * $N + ($m - 1)]) * (abs ($q) + abs ($r)) < eps * (abs ($p) * (abs ($H[($m - 1) * $N + ($m - 1)]) + abs ($z) + abs ($H[($m + 1) * $N + ($m + 1)])));
		    }

		  for $i ($m + 2 .. $n)
		    {
		      $H[$i * $N + ($i - 2)] = 0;
		      $H[$i * $N + ($i - 3)] = 0
			if $i > $m + 2;
		    }

		  # Double QR step.
		  for $k ($m .. $n - 1)
		    {
		      $not_last = ($k != $n - 1);

		      if ($k != $m)
			{
			  $p = $H[$k * $N + ($k - 1)];
			  $q = $H[($k + 1) * $N + ($k - 1)];
			  $r = $not_last ? $H[($k + 2) * $N + ($k - 1)] : 0;
			  $x = abs ($p) + abs ($q) + abs ($r);

			  next if $x == 0;

			  $p = $p / $x;
			  $q = $q / $x;
			  $r = $r / $x;
			}

		      $s = sign (sqrt ($p * $p + $q * $q + $r * $r), $p);
		      if ($s != 0)
			{
			  if ($k != $m)
			    {
			      $H[$k * $N + ($k - 1)] = 0 - $s * $x;
			    }
			  elsif ($l != $m)
			    {
			      $H[$k * $N + ($k - 1)] = - $H[$k * $N + ($k - 1)];
			    }

			  $p = $p + $s;
			  $x = $p / $s;
			  $y = $q / $s;
			  $z = $r / $s;
			  $q = $q / $p;
			  $r = $r / $p;

			  if ($not_last)
			    {
			      # Row modification.
			      for $j ($k .. $end)
				{
				  $p = $H[$k * $N + $j] + $q * $H[($k + 1) * $N + $j] + $r * $H[($k + 2) * $N + $j];

				  $H[$k * $N + $j] -= $p * $x;
				  $H[($k + 1) * $N + $j] -= $p * $y;
				  $H[($k + 2) * $N + $j] -= $p * $z;
				}

			      # Column modification.
			      for $i (0 .. min ($n, $k + 3))
				{
				  $p = $x * $H[$i * $N + $k] + $y * $H[$i * $N + ($k + 1)] + $z * $H[$i * $N + ($k + 2)];

				  $H[$i * $N + $k] -= $p;
				  $H[$i * $N + ($k + 1)] -= $p * $q;
				  $H[$i * $N + ($k + 2)] -= $p * $r;
				}

			      # Accumulate transformations.
			      for $i ($lo .. $hi)
				{
				  $p = $x * $$Z[$k][$i] + $y * $$Z[$k + 1][$i] + $z * $$Z[$k + 2][$i];

				  $$Z[$k][$i] -= $p;
				  $$Z[$k + 1][$i] -= $p * $q;
				  $$Z[$k + 2][$i] -= $p * $r;
				}
			    }
			  else
			    {
			      # Row modification.
			      for $j ($k .. $end)
				{
				  $p = $H[$k * $N + $j] + $q * $H[($k + 1) * $N + $j];

				  $H[$k * $N + $j] -= $p * $x;
				  $H[($k + 1) * $N + $j] -= $p * $y;
				}

			      # Column modification.
			      for $i (0 .. min ($n, $k + 3))
				{
				  $p = $x * $H[$i * $N + $k] + $y * $H[$i * $N + ($k + 1)];

				  $H[$i * $N + $k] -= $p;
				  $H[$i * $N + ($k + 1)] -= $p * $q;
				}

			      # Accumulate transformations.
			      for $i (0 .. $end)
				{
				  $p = $x * $$Z[$k][$i] + $y * $$Z[$k + 1][$i];

				  $$Z[$k][$i] -= $p;
				  $$Z[$k + 1][$i] -= $p * $q;
				}
			    }
			}
		    }
		}
	    }

	  # Backsubstitute to find vectors of upper triangular form.
	  return if $norm == 0;

	  for $n (reverse (0 .. $end))
	    {
	      $p = $$d[$n];
	      $q = $$e[$n];

	      if ($q == 0)
		{
		  # Real vector.
		  $m = $n;
		  $H[$n * $N + $n] = 1;

		  for $i (reverse (0 .. $n - 1))
		    {
		      $w = $H[$i * $N + $i] - $p;
		      $r = 0;

		      for $j ($m .. $n)
			{
			  $r += $H[$i * $N + $j] * $H[$j * $N + $n];
			}

		      if ($$e[$i] < 0)
			{
			  $z = $w;
			  $s = $r;
			}
		      else
			{
			  $m = $i;

			  if ($$e[$i] == 0)
			    {
			      $H[$i * $N + $n] = ($w != 0 ?
						  0 - $r / $w :
						  0 - $r / (eps * $norm));
			    }
			  else
			    {
			      # Solve real equations.
			      $x = $H[$i * $N + ($i + 1)];
			      $y = $H[($i + 1) * $N + $i];
			      $q = ($$d[$i] - $p) ** 2 + $$e[$i] ** 2;
			      $t = ($x * $s - $z * $r) / $q;

			      $H[$i * $N + $n] = $t;
			      $H[($i + 1) * $N + $n] = (abs ($x) > abs ($z) ?
						(0 - $r - $w * $t) / $x :
						(0 - $s - $y * $t) / $z);
			    }

			  # Overflow control.
			  $t = abs ($H[$i * $N + $n]);
			  if ((eps * $t) * $t > 1)
			    {
			      for $j ($i .. $n)
				{
				  $H[$j * $N + $n] /= $t;
				}
			    }
			}
		    }
		}
	      elsif ($q < 0)
		{
		  # Complex vector.
		  $m = $n - 1;

		  # Last vector component chosen imaginary so that
		  # eigenvector matrix is triangular.
		  if (abs ($H[$n * $N + ($n - 1)]) > abs ($H[($n - 1) * $N + $n]))
		    {
		      $H[($n - 1) * $N + ($n - 1)] = $q / $H[$n * $N + ($n - 1)];
		      $H[($n - 1) * $N + $n] = ($p - $H[$n * $N + $n]) / $H[$n * $N + ($n - 1)];
		    }
		  else
		    {
		      ($H[($n - 1) * $N + ($n - 1)], $H[($n - 1) * $N + $n])
			= cdiv (0, - $H[($n - 1) * $N + $n],
				$H[($n - 1) * $N + ($n - 1)] - $p, $q);
		    }

		  $H[$n * $N + ($n - 1)] = 0;
		  $H[$n * $N + $n] = 1;

		  for $i (reverse (0 .. $n - 2))
		    {
		      $w = $H[$i * $N + $i] - $p;

		      $ra = 0;
		      $sa = 0;

		      for $j ($m .. $n)
			{
			  $ra += $H[$i * $N + $j] * $H[$j * $N + ($n - 1)];
			  $sa += $H[$i * $N + $j] * $H[$j * $N + $n];
			}

		      if ($$e[$i] < 0)
			{
			  $z = $w;
			  $r = $ra;
			  $s = $sa;
			}
		      else
			{
			  $m = $i;

			  if ($$e[$i] == 0)
			    {
			      ($H[$i * $N + ($n - 1)], $H[$i * $N + $n])
				= cdiv (- $ra, - $sa, $w, $q);
			    }
			  else
			    {
			      # Solve complex equations.
			      $x = $H[$i * $N + ($i + 1)];
			      $y = $H[($i + 1) * $N + $i];

			      $vr = ($$d[$i] - $p) ** 2 + $$e[$i] ** 2 - $q ** 2;
			      $vi = ($$d[$i] - $p) * 2 * $q;

			      if ($vr == 0 && $vi == 0)
				{
				  $vr = eps * $norm * (abs ($w) + abs ($q) + abs ($x) + abs ($y) + abs ($z));
				}

			      ($H[$i * $N + ($n - 1)], $H[$i * $N + $n])
				= cdiv ($x * $r - $z * $ra + $q * $sa,
					$x * $s - $z * $sa - $q * $ra,
					$vr,
					$vi);

			      if (abs ($x) > (abs ($z) + abs ($q)))
				{
				  $H[($i + 1) * $N + ($n - 1)] = (0 - $ra - $w * $H[$i * $N + ($n - 1)] + $q * $H[$i * $N + $n]) / $x;
				  $H[($i + 1) * $N + $n] = (0 - $sa - $w * $H[$i * $N + $n] - $q * $H[$i * $N + ($n - 1)]) / $x;
				}
			      else
				{
				  ($H[($i + 1) * $N + ($n - 1)], $H[($i + 1) * $N + $n])
				    = cdiv (0 - $r - $y * $H[$i * $N + ($n - 1)],
					    0 - $s - $y * $H[$i * $N + $n],
					    $z,
					    $q);
				}
			    }

			  # Overflow control.
			  $t = max (abs ($H[$i * $N + ($n - 1)]), abs ($H[$i * $N + $n]));
			  if ((eps * $t) * $t > 1)
			    {
			      for $j ($i .. $n)
				{
				  $H[$j * $N + ($n - 1)] /= $t;
				  $H[$j * $N + $n] /= $t;
				}
			    }
			}
		    }
		}
	    }

	  # Vectors of isolated roots.
	  for $i (0 .. $end)
	    {
	      if ($i < $lo || $i > $hi)
		{
		  for $j ($i .. $end)
		    {
		      $$Z[$j][$i] = $H[$i * $N + $j];
		    }
		}
	    }

	  # Multiply by transformation matrix to give
	  # vectors of original full matrix.
	  for $j (reverse ($lo .. $end))
	    {
	      $m = min ($j, $hi);

	      for $i ($lo .. $hi)
		{
		  $z = 0;

		  for $k ($lo .. $m)
		    {
		      $z += $$Z[$k][$i] * $H[$k * $N + $j];
		    }

		  $$Z[$j][$i] = $z;
		}
	    }
	}

      # Form the eigenvectors of a real general matrix by back
      # transforming those of the corresponding balanced matrix
      # determined by 'balance'.
      #
      # See EiSPACK procedure 'balbak'.
      if ($prop{balance})
	{
	  my ($i, $j, $k);

	  # Undo scaling.
	  for $i ($lo .. $hi)
	    {
	      my $s = $scale[$i];

	      for $j (0 .. $end)
		{
		  $$Z[$j][$i] *= $s;
		}
	    }

	  # Undo permutations.
	  for $i (reverse (0 .. $lo - 1))
	    {
	      $k = $perm[$i];
	      if ($k != $i)
		{
		  for $j (0 .. $end)
		    {
		      ($$Z[$j][$i], $$Z[$j][$k])
			= ($$Z[$j][$k], $$Z[$j][$i]);
		    }
		}
	    }

	  for $i ($hi + 1 .. $end)
	    {
	      $k = $perm[$i];
	      if ($k != $i)
		{
		  for $j (0 .. $end)
		    {
		      ($$Z[$j][$i], $$Z[$j][$k])
			= ($$Z[$j][$k], $$Z[$j][$i]);
		    }
		}
	    }
	}
    }

  # Create complex eigenvalues.
  for my $i (0 .. $end)
    {
      $$d[$i] = Math::Complex->make ($$d[$i], $$e[$i])
	if $$e[$i] != 0;
    }

  # Normalize eigenvectors.
  $self->normalize
    if $prop{normalize};

  # Make first non-zero vector element a positive number.
  if ($prop{positive})
    {
      my ($i, $j, $k);

      for $j (0 .. $end)
	{
	  for $i (0 .. $end)
	    {
	      next if $$Z[$j][$i] == 0;

	      if ($$Z[$j][$i] < 0)
		{
		  for $k ($i .. $end)
		    {
		      $$Z[$j][$k] = - $$Z[$j][$k];
		    }
		}

	      last;
	    }
	}
    }

  # Return object.
  $self;
}

# Balance a real matrix.
#
# See LAPACK procedure 'dgebal'.
sub gebal ($$%)
{
  # Arguments.
  my ($n, $a, %opt) = @_;

  # Default options.
  $opt{permute} //= 1;
  $opt{scale} //= 1;
  $opt{low} //= undef;
  $opt{high} //= undef;

  # Index of last row/column.
  my $end = $n - 1;

  # Row and column indices of the beginning and end
  # of the principal sub-matrix.
  my $k = 0;
  my $l = $end;

  # Permutation vector.
  my @perm = (0 .. $end);

  # Scaling vector.
  my @scale = map (1, 0 .. $end);

  if ($n > 1)
    {
      # Isolate eigenvalues.
      if ($opt{permute})
	{
	  # Search for rows isolating an eigenvalue
	  # and push them down.
	L:
	  {
	    for my $j (reverse (0 .. $l))
	      {
		my $swap = 1;

		for my $i (0 .. $l)
		  {
		    $swap = 0 if $i != $j && $$a[$j * $n + $i] != 0;
		  }

		next unless $swap;

		if ($j != $l)
		  {
		    # Exchange row and column.
		    @perm[$j, $l] = @perm[$l, $j];

		    blas_swap ($l + 1, $a, $a,
			       x_ind => $j,
			       x_incr => $n,
			       y_ind => $l,
			       y_incr => $n);
		    blas_swap ($n - $k, $a, $a,
			       x_ind => $j * $n + $k,
			       x_incr => 1,
			       y_ind => $l * $n + $k,
			       y_incr => 1);
		  }

		$l -= 1;
		next L;
	      }
	  }

	  # Search for columns isolating an eigenvalue
	  # and push them left.
	K:
	  {
	    for my $j ($k .. $l)
	      {
		my $swap = 1;

		for my $i ($k .. $l)
		  {
		    $swap = 0 if $i != $j && $$a[$i * $n + $j] != 0;
		  }

		next unless $swap;

		if ($j != $k)
		  {
		    # Exchange row and column.
		    @perm[$j, $k] = @perm[$k, $j];

		    blas_swap ($l + 1, $a, $a,
			       x_ind => $j,
			       x_incr => $n,
			       y_ind => $k,
			       y_incr => $n);
		    blas_swap ($n - $k, $a, $a,
			       x_ind => $j * $n + $k,
			       x_incr => 1,
			       y_ind => $k * $n + $k,
			       y_incr => 1);
		  }

		$k += 1;
		next K;
	      }
	  }
	}

      # Balance the sub-matrix in rows k to l.
      if ($opt{scale})
	{
	  my ($b, $c, $ca, $f, $g, $r, $ra, $s, $no_conv,
	      $sfmin1, $sfmax1, $sfmin2, $sfmax2);

	  # Scale factors are powers of two.
	  $b = 2;

	  if (1)
	    {
	      my $base = 2;
	      my $eps = DBL_EPSILON;
	      my $small = 1 / DBL_MAX;
	      # Safe minimum, such that 1/sfmin does not overflow.
	      my $sfmin = ($small >= DBL_MIN ? $small * (1 + $eps) : DBL_MIN);

	      $sfmin1 = $sfmin / ($eps * $base);
	      $sfmax1 = 1 / $sfmin1;
	      $sfmin2 = $sfmin1 * $b;
	      $sfmax2 = 1 / $sfmin2;
	    }

	  # Iterative loop for norm reduction.
	  while (1)
	    {
	      $no_conv = 0;

	      for my $i ($k .. $l)
		{
		  $c = blas_norm ($l - $k + 1, $a,
				  norm => BLAS_TWO_NORM,
				  x_ind => $k * $n + $i,
				  x_incr => $n);
		  $r = blas_norm ($l - $k + 1, $a,
				  norm => BLAS_TWO_NORM,
				  x_ind => $i * $n + $k,
				  x_incr => 1);
		  (undef, $ca) = blas_amax_val ($l + 1, $a,
						x_ind => $i,
						x_incr => $n);
		  (undef, $ra) = blas_amax_val ($n - $k, $a,
						x_ind => $i * $n + $k,
						x_incr => 1);

		  # Guard against zero c or r due to underflow.
		  next if $c == 0 || $r == 0;

		  $s = $c + $r;
		  $f = 1;

		  $g = $r / $b;
		  while ($c < $g
			 && max ($f, max ($c, $ca)) < $sfmax2
			 && min ($g, min ($r, $ra)) > $sfmin2)
		    {
		      croak ('Infinite loop')
			if isnan ($f + $c + $ca + $g + $r + $ra);

		      $f *= $b;
		      $c *= $b;
		      $ca *= $b;
		      $g /= $b;
		      $r /= $b;
		      $ra /= $b;
		    }

		  $g = $c / $b;
		  while ($g >= $r
			 && max ($r, $ra) < $sfmax2
			 && min ($g, min ($f, min ($c, $ca))) > $sfmin2)
		    {
		      $f /= $b;
		      $c /= $b;
		      $ca /= $b;
		      $g /= $b;
		      $r *= $b;
		      $ra *= $b;
		    }

		  # Now balance.
		  if (! (($c + $r) >= 0.95 * $s
			 || ($f < 1 && $scale[$i] < 1 && $f * $scale[$i] <= $sfmin1)
			 || ($f > 1 && $scale[$i] > 1 && $scale[$i] >= $sfmax1 / $f)))
		    {
		      blas_rscale ($n - $k, $a,
				   alpha => $f,
				   x_ind => $i * $n + $k,
				   x_incr => 1);
		      blas_scale ($l + 1, $a,
				  alpha => $f,
				  x_ind => $i,
				  x_incr => $n);

		      $scale[$i] *= $f;
		      $no_conv = 1;
		    }
		}

	      last unless $no_conv;
	    }
	}
    }

  ${$opt{low}} = $k
    if ref ($opt{low});

  ${$opt{high}} = $l
    if ref ($opt{high});

  return (\@perm, \@scale);
}

# Normalize eigenvectors.
sub normalize
{
  my $self = shift;

  # Work variables.
  my $len;

  for my $v (@{ $$self{vector} })
    {
      $len = blas_norm (@$v, $v, norm => BLAS_TWO_NORM);
      blas_rscale (@$v, $v, alpha => $len) if $len != 0;
    }

  # Return object.
  $self;
}

# Sort eigenvalues and corresponding eigenvectors.
sub sort
{
  my $self = shift;
  my $order = shift // 'abs_desc';

  # Permutation vector.
  my @p = ();

  if ($order =~ m/\Avec_/)
    {
      my $Z = $$self{vector};

      @p = ($order eq 'vec_desc' ?
	    sort { _cmp_vec ($$Z[$b], $$Z[$a]) } 0 .. $#$Z :
	    ($order eq 'vec_asc' ?
	     sort { _cmp_vec ($$Z[$a], $$Z[$b]) } 0 .. $#$Z :
	     croak ("Invalid argument")));
    }
  elsif (grep (ref ($_), @{ $$self{value} }))
    {
      # Consider complex eigenvalues.
      my (@d, @e, @m) = ();

      for (@{ $$self{value} })
	{
	  if (ref ($_))
	    {
	      push (@d, $_->Re);
	      push (@e, $_->Im);
	      push (@m, abs ($_));
	    }
	  else
	    {
	      push (@d, $_);
	      push (@e, 0);
	      push (@m, abs ($_));
	    }
	}

      @p = ($order eq 'abs_desc' ?
	    sort { abs ($d[$b]) <=> abs ($d[$a]) || $m[$b] <=> $m[$a] || $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
	    ($order eq 'abs_asc' ?
	     sort { abs ($d[$a]) <=> abs ($d[$b]) || $m[$a] <=> $m[$b] || $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
	     ($order eq 'norm_desc' ?
	      sort { $m[$b] <=> $m[$a] || $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
	      ($order eq 'norm_asc' ?
	       sort { $m[$a] <=> $m[$b] || $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
	       ($order eq 'desc' ?
		sort { $d[$b] <=> $d[$a] || $e[$b] <=> $e[$a] } 0 .. $#d :
		($order eq 'asc' ?
		 sort { $d[$a] <=> $d[$b] || $e[$a] <=> $e[$b] } 0 .. $#d :
		 croak ("Invalid argument")))))));
    }
  else
    {
      # Only real eigenvalues.
      my $d = $$self{value};

      @p = ($order eq 'abs_desc' || $order eq 'norm_desc' ?
	    sort { abs ($$d[$b]) <=> abs ($$d[$a]) || $$d[$b] <=> $$d[$a] } 0 .. $#$d :
	    ($order eq 'abs_asc' || $order eq 'norm_asc' ?
	     sort { abs ($$d[$a]) <=> abs ($$d[$b]) || $$d[$a] <=> $$d[$b] } 0 .. $#$d :
	     ($order eq 'desc' ?
	      sort { $$d[$b] <=> $$d[$a] } 0 .. $#$d :
	      ($order eq 'asc' ?
	       sort { $$d[$a] <=> $$d[$b] } 0 .. $#$d :
	       croak ("Invalid argument")))));
    }

  # Reorder eigenvalues and corresponding eigenvectors.
  if (@p > 0)
    {
      my $ref;

      $ref = $$self{value};
      @$ref = @$ref[@p];

      $ref = $$self{vector};
      @$ref = @$ref[@p];
    }

  # Return object.
  $self;
}

# Compare two eigenvectors.
sub _cmp_vec
{
  my ($u, $v) = @_;

  croak ("Invalid argument")
    if $#$u != $#$v;

  my $d = 0;

  for my $i (0 .. $#$u)
    {
      last if $d = ($$u[$i] <=> $$v[$i]);
    }

  $d;
}

# Return one or more eigenvalues.
sub value
{
  my $self = shift;

  @_ > 0 ? @{ $$self{value} }[@_] : @{ $$self{value} };
}

*values = \&value;

# Return one or more eigenvectors.
sub vector
{
  my $self = shift;

  @_ > 0 ? @{ $$self{vector} }[@_] : @{ $$self{vector} };
}

*vectors = \&vector;

1;

__END__