view release on metacpan or  search on metacpan

t/eigen_3x3.t  view on Meta::CPAN

my $ap1_2 = $al1 * $aev1;    # lambda *x
print "Original computation of A*ev1:\n$ap1_1 Scaled eigenvector:\n$ap1_2"
    if $DEBUG2;
ok_matrix( $ap1_1, $ap1_2, 'eigenvectors match');

print "Direct diagonalization...\n" if $DEBUG;
my ($L12, $V12) = $symm->sym_diagonalize();
print "eigenvalues L:\n$L12 eigenvectors:\n$V12" if $DEBUG2;
ok_eigenvectors($symm, $L12, $V12);
ok_matrix_orthogonal($V12);
# Double check the equality
ok_matrix( $L12, $L1);
ok_matrix( $V12, $V1);

#
# Now test the eigenvalues only computations...
#
print "Recomputing: Eigenvalues only.\n 3x3\n" if $DEBUG;
my $altT = $symm->householder_tridiagonal();
ok_matrix( $altT, $T,'householder_tridiagonal works');
my $altL1 = $altT->tri_eigenvalues();
ok_matrix( $altL1, $L1,'tri_eigenvalues works');
my $altL12 = $symm->sym_eigenvalues();
ok_matrix( $altL12, $L12, 'sym_eigenvalues works');

__END__

# Attempt:
# Obtain the eigenvectors when eigenvalues are known
# using inverse iteration.
# We solve (M - lambda * I) * b(k+1) = b(k)
#  with b(0) a random unit vector and b(k) is
#  normalized at each step.
# This should converge towards the eigenvector,
# but there are problems:
#  - for some value, there is not convergence
#  (the above system is rather singular, so...)
#  - the solution can oscillate between v and -v (?)
# Rodolphe Ortalo, 99/06/14
sub obtain_eigenvector ($$)
{
  my ($M, $eigenvalue) = @_;
  # Form the linear system A - lamda1 * I
  my $inv_it = $M->shadow();
  $inv_it->one();
  $inv_it->multiply_scalar($inv_it, (-1.0 * $eigenvalue));
  $inv_it->add($M, $inv_it);
  print "Linear system matrix:\n $inv_it" if $DEBUG2;
  # Creates a random vector
  my ($rows, $cols) = $inv_it->dim();
  my $b = Math::MatrixReal->new($rows, 1);
  for (my $i = 1; $i <= $rows; $i++)
    {
      $b->assign($i, 1, rand());
    }
  # Normalize it
  my $l = $b->length();
  $b->multiply_scalar($b, (1.0 / $l));
  # Now do LR decomposition for linear system
  my $inv_it_LR = $inv_it->decompose_LR();
  # Check iterations
  my $iter = 0;
  my $delta;
  do {
    my ($dim, $b_base, $base) = $inv_it_LR->solve_LR($b);
    # Normalize
    my $l = $b_base->length();
    $b_base->multiply_scalar($b_base, (-1.0 / $l));
#    print "b_base=\n$b_base";
    $b->subtract($b_base,$b);
    $delta = $b->norm_one();
    print "delta=$delta\n";
    $b = $b_base;
  } while (($delta >= 1e-10) && ($iter++ <= 10));
  return $b;
}
#
# Now, try to find one eigenvector again...
# (Using Steffen's functions...:-)
#
my $ev = obtain_eigenvector($symm, $al1);
print "Real ev:\n $aev1 Found ev:\n $ev" if $DEBUG;
ok_matrix(15, $ev, $aev1);
ok_matrix(16, obtain_eigenvector($symm, $L1->element(2,1)),
	  $V1->column(2));
ok_matrix(17, obtain_eigenvector($symm, $L1->element(3,1)),
	  $V1->column(3));