forked from leto/math--matrixreal
/
eigen_7x7.t
62 lines (57 loc) · 1.91 KB
/
eigen_7x7.t
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
use Test::More tests => 10;
use Math::MatrixReal;
use lib 't/lib';
use Test::Matrices qw{ok_eigenvectors ok_matrix ok_matrix_orthogonal};
no lib 't/lib';
my $DEBUG = $Test::Matrices::DEBUG = 0;
my $DEBUG2 = 0;
#
# Trying some matrixes creation extracted from the pod...
#
my $matrix77_b = Math::MatrixReal->new_from_string(<<'MATRIX');
[ 1 7 -12 6 -9 0 1 ]
[ 0 5 0 0 0 0 0 ]
[ 0 0 1 4 0 0 0 ]
[ 0 0 0 1 0 0 0 ]
[ 12 0 0 0 5 0 4 ]
[ 0 3 0 8 0 1 0 ]
[ 1 0 0 0 0 0 -5 ]
MATRIX
print "$matrix77_b" if $DEBUG2;
#
# Tests eigenvalues & eigenvectors computation
#
#
# Redo things with the 7x7 matrix
#
my $symm2 = $matrix77_b + ~$matrix77_b;
print "Matrix 7x7 for eigenvalues & eigenvectors computation:\n" if $DEBUG;
print "$symm2" if $DEBUG2;
print "Householder reduction...\n" if $DEBUG;
my ($T2, $Q2) = $symm2->householder();
print "T2=\n$T2 Q2=\n$Q2" if $DEBUG2;
print "Is Q2 orthogonal?\n" if $DEBUG;
ok_matrix_orthogonal($Q2);
ok_matrix($symm2, $Q2 * $T2 * ~$Q2, 'householder reduction for 7x7');
print "Diagonalization of tridiagonal...\n" if $DEBUG;
my ($L, $V) = $T2->tri_diagonalize($Q2);
print "eigenvalues L:\n$L eigenvectors:\n$V" if $DEBUG2;
ok_eigenvectors($symm2, $L, $V);
print "Direct diagonalization...\n" if $DEBUG;
my ($L_2, $V_2) = $symm2->sym_diagonalize();
print "eigenvalues L:\n$L_2 eigenvectors:\n$V_2" if $DEBUG2;
ok_eigenvectors($symm2, $L_2, $V_2);
ok_matrix_orthogonal( $V_2);
# Double check the equality
ok_matrix( $L_2, $L, 'equality');
ok_matrix( $V_2, $V, 'equality');
#
# Now test the eigenvalues only computations...
#
print "Recomputing: Eigenvalues only.\n 7x7\n" if $DEBUG;
my $altT2 = $symm2->householder_tridiagonal();
ok_matrix( $altT2, $T2, 'householder_tridiagonal for 7x7');
my $altL = $altT2->tri_eigenvalues();
ok_matrix( $altL, $L, 'tri_eigenvalues for 7x7');
my $altL_2 = $symm2->sym_eigenvalues();
ok_matrix( $altL_2, $L_2, 'sym_eigenvalues for 7x7');