File Coverage

lib/Math/MatrixDecomposition/Eigen.pm

Criterion	Covered	Total	%
statement	440	622	70.7
branch	85	210	40.4
condition	27	82	32.9
subroutine	15	17	88.2
pod	7	7	100.0
total	574	938	61.1

## 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 Math::Complex qw();
use Scalar::Util qw(looks_like_number);

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

BEGIN
{
  our $VERSION = '1.04';
  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 vectors).
  $$self{vector} //= [];

  # Matrix $Z contains the eigenvectors.
  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 (an array of row vectors).
      my @H = map ([], 0 .. $end);

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

      # 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.
      #
      # See EiSPACK procedure 'balanc' and LAPACK procedure 'dgebal'.
      if ($prop{balance})
	{
	  @perm = (0 .. $end);
	  @scale = map (1, 0 .. $end);

	  # Work variables.
	  my ($i, $j, $k, $l,
	      $b, $b2, $c, $f, $g, $r, $s, $no_conv);

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

	  $k = 0;
	  $l = $end;

	  if ($l > 0)
	    {
	      # Search for rows isolating an eigenvalue
	      # and push them down.
	    L:
	      {
		for $j (reverse (0 .. $l))
		  {
		    $r = 0;

		    for $i (0 .. $l)
		      {
			$r = 1 if $i != $j && $H[$j][$i] != 0;
		      }

		    next if $r != 0;

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

			for $i (0 .. $l)
			  {
			    ($H[$i][$j], $H[$i][$l])
			      = ($H[$i][$l], $H[$i][$j]);
			  }

			for $i ($k .. $end)
			  {
			    ($H[$j][$i], $H[$l][$i])
			      = ($H[$l][$i], $H[$j][$i]);
			  }
		      }

		    $l -= 1;
		    next L;
		  }
	      }

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

		    for $i ($k .. $l)
		      {
			$c = 1 if $i != $j && $H[$i][$j] != 0;
		      }

		    next if $c != 0;

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

			for $i (0 .. $l)
			  {
			    ($H[$i][$j], $H[$i][$k])
			      = ($H[$i][$k], $H[$i][$j]);
			  }

			for $i ($k .. $end)
			  {
			    ($H[$j][$i], $H[$k][$i])
			      = ($H[$k][$i], $H[$j][$i]);
			  }
		      }

		    $k += 1;
		    next K;
		  }
	      }

	      ## Now balance the sub-matrix in rows k to l.

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

		  for $i ($k .. $l)
		    {
		      $c = 0;
		      $r = 0;

		      for $j ($k .. $l)
			{
			  next if $j == $i;

			  $c += abs ($H[$j][$i]);
			  $r += abs ($H[$i][$j]);
			}

		      # 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)
			{
			  $f *= $b;
			  $c *= $b2;
			}

		      $g = $r * $b;
		      while ($c >= $g)
			{
			  $f /= $b;
			  $c /= $b2;
			}

		      # Now balance.
		      if (($c + $r) / $f < 0.95 * $s)
			{
			  $g = 1 / $f;

			  for $j ($k .. $end)
			    {
			      $H[$i][$j] *= $g;
			    }

			  for $j (0 .. $l)
			    {
			      $H[$j][$i] *= $f;
			    }

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

		  last unless $no_conv;
		}
	    }

	  $lo = $k;
	  $hi = $l;
	}

      # 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][$m - 1]);
		}

	      next if $scale == 0;

	      $h = 0;

	      for $i (reverse ($m .. $hi))
		{
		  $ort[$i] = $H[$i][$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][$j];
		    }

		  $f /= $h;

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

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

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

		  $f /= $h;

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

	      $ort[$m] *= $scale;
	      $H[$m][$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][$m - 1] == 0;

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

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

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

		  $g = ($g / $ort[$m]) / $H[$m][$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][$i];
		  $$e[$i] = 0;
		}
	    }

	  # Compute matrix norm.
	  $norm = 0;

	  for $i (0 .. $end)
	    {
	      for $j ($i .. $end)
		{
		  $norm += abs ($H[$i][$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][$l - 1]) + abs ($H[$l][$l]);
		  $s = $norm
		    if $s == 0;

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

	      $x = $H[$n][$n];

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

		  $$d[$n] = $H[$n][$n];
		  $$e[$n] = 0;

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

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

		  $H[$n][$n] = $x + $t;
		  $H[$n - 1][$n - 1] = $y + $t;
		  $x = $H[$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 - 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][$j];
			  $H[$n - 1][$j] = $q * $z + $p * $H[$n][$j];
			  $H[$n][$j] = $q * $H[$n][$j] - $p * $z;
			}

		      # Column modification.
		      for $i (0 .. $n)
			{
			  $z = $H[$i][$n - 1];
			  $H[$i][$n - 1] = $q * $z + $p * $H[$i][$n];
			  $H[$i][$n] = $q * $H[$i][$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 - 1];
		  $w = $H[$n][$n - 1] * $H[$n - 1][$n];

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

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

		      $s = abs ($H[$n][$n - 1]) + abs ($H[$n - 1][$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][$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][$m];
		      $r = $x - $z;
		      $s = $y - $z;
		      $p = ($r * $s - $w) / $H[$m + 1][$m] + $H[$m][$m + 1];
		      $q = $H[$m + 1][$m + 1] - $z - $r - $s;
		      $r = $H[$m + 2][$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][$m - 1]) * (abs ($q) + abs ($r)) < eps * (abs ($p) * (abs ($H[$m - 1][$m - 1]) + abs ($z) + abs ($H[$m + 1][$m + 1])));
		    }

		  for $i ($m + 2 .. $n)
		    {
		      $H[$i][$i - 2] = 0;
		      $H[$i][$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][$k - 1];
			  $q = $H[$k + 1][$k - 1];
			  $r = $not_last ? $H[$k + 2][$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][$k - 1] = 0 - $s * $x;
			    }
			  elsif ($l != $m)
			    {
			      $H[$k][$k - 1] = - $H[$k][$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][$j] + $q * $H[$k + 1][$j] + $r * $H[$k + 2][$j];

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

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

				  $H[$i][$k] -= $p;
				  $H[$i][$k + 1] -= $p * $q;
				  $H[$i][$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][$j] + $q * $H[$k + 1][$j];

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

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

				  $H[$i][$k] -= $p;
				  $H[$i][$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] = 1;

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

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

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

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

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

			  # Overflow control.
			  $t = abs ($H[$i][$n]);
			  if ((eps * $t) * $t > 1)
			    {
			      for $j ($i .. $n)
				{
				  $H[$j][$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 - 1]) > abs ($H[$n - 1][$n]))
		    {
		      $H[$n - 1][$n - 1] = $q / $H[$n][$n - 1];
		      $H[$n - 1][$n] = ($p - $H[$n][$n]) / $H[$n][$n - 1];
		    }
		  else
		    {
		      ($H[$n - 1][$n - 1], $H[$n - 1][$n])
			= cdiv (0, - $H[$n - 1][$n],
				$H[$n - 1][$n - 1] - $p, $q);
		    }

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

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

		      $ra = 0;
		      $sa = 0;

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

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

			  if ($$e[$i] == 0)
			    {
			      ($H[$i][$n - 1], $H[$i][$n])
				= cdiv (- $ra, - $sa, $w, $q);
			    }
			  else
			    {
			      # Solve complex equations.
			      $x = $H[$i][$i + 1];
			      $y = $H[$i + 1][$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 - 1], $H[$i][$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 - 1] = (0 - $ra - $w * $H[$i][$n - 1] + $q * $H[$i][$n]) / $x;
				  $H[$i + 1][$n] = (0 - $sa - $w * $H[$i][$n] - $q * $H[$i][$n - 1]) / $x;
				}
			      else
				{
				  ($H[$i + 1][$n - 1], $H[$i + 1][$n])
				    = cdiv (0 - $r - $y * $H[$i][$n - 1],
					    0 - $s - $y * $H[$i][$n],
					    $z,
					    $q);
				}
			    }

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

	  # Vectors of isolated roots.
	  for $i (0 .. $end)
	    {
	      if ($i < $lo || $i > $hi)
		{
		  for $j ($i .. $end)
		    {
		      $$Z[$j][$i] = $H[$i][$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][$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 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;
}

# 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__