use Test::More tests => 15;
use File::Spec;
use lib File::Spec->catfile("..","lib");
use Math::MatrixReal;

do 'funcs.pl';

my $DEBUG2 = 0;
my $bigsize = 30; # Size of big matrix tests (be careful: n^3!)

# test on random bigger matrix
print "Matrix ".$bigsize."x$bigsize for eigenvalues & eigenvectors computation:\n" if $DEBUG;
# Creates a random matrix
my $big = Math::MatrixReal->new_random($bigsize, { symmetric => 1 });

# Tests eigenvalues & eigenvectors computation

print "Householder reduction...\n" if $DEBUG;
my ($Tbig, $Qbig) = $big->householder();
print "Is Qbig orthogonal?\n" if $DEBUG;
ok_matrix_orthogonal( $Qbig);
ok_matrix( $big, $Qbig * $Tbig * ~$Qbig);
print "Diagonalization of tridiagonal...\n" if $DEBUG;
my ($Lbig, $Vbig) = $Tbig->tri_diagonalize($Qbig);
ok_eigenvectors( $big, $Lbig, $Vbig);
ok_matrix_orthogonal( $Vbig);

print "Direct diagonalization...\n" if $DEBUG;
my ($Lbig_2, $Vbig_2) = $big->sym_diagonalize();
print "eigenvalues L:\n$Lbig_2 eigenvectors:\n$Vbig_2" if $DEBUG2;
ok_eigenvectors($big, $Lbig_2, $Vbig_2);
ok_matrix_orthogonal( $Vbig_2);
# Double check the equality
ok_matrix( $Lbig_2, $Lbig);
ok_matrix( $Vbig_2, $Vbig);

#
# Now test the eigenvalues only computations...
#
print "Recomputing: Eigenvalues only.\n ".$bigsize."x".$bigsize."\n" if $DEBUG;
my $altTbig = $big->householder_tridiagonal();
ok_matrix( $altTbig, $Tbig);
my $altLbig = $altTbig->tri_eigenvalues();
ok_matrix( $altLbig, $Lbig);
my $altLbig_2 = $big->sym_eigenvalues();
ok_matrix( $altLbig_2, $Lbig_2);

##############
#### lower tri
my $eigen = Math::MatrixReal->new_from_string(<<MAT);
[  0.000000000000E+00 ]
[  3.000000000000E+00 ]
[  4.000000000000E+00 ]
[  5.000000000000E+00 ]
[  1.000000000000E+00 ]
MAT
$matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
[ 0 0 0 0 0 ]
[ 0 3 0 0 0 ]
[ 0 0 4 0 0 ]
[ 1 0 0 5 0 ]
[ 1 1 1 1 1 ]
MATRIX
ok_matrix( $eigen, $matrix->eigenvalues );

$matrix = $eigen->new_from_rows ( [[1,0,0],[0,2,0],[0,0,3]] );
$eigen = $eigen->new_from_string(<<MAT);
[ 1 ]
[ 2 ]
[ 3 ]
MAT
ok_matrix($eigen, $matrix->eigenvalues );

####################
## upper tri
$matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
[ 1 0 0 0 1 ]
[ 0 2 0 0 2 ]
[ 0 0 3 0 0 ]
[ 0 0 0 4 0 ]
[ 0 0 0 0 5 ]
MATRIX
$eigen = Math::MatrixReal->new_from_string(<<MAT);
[ 1 ]
[ 2 ]
[ 3 ]
[ 4 ]
[ 5 ]
MAT

ok_matrix( $eigen, $matrix->eigenvalues );

######################
#### diag
$matrix = $matrix->new_diag ( [ 10, 20, 30 ] );
$eigen =  $matrix->new_from_cols ( [ [ 10, 20, 30 ] ] );
ok_matrix( $eigen, $matrix->eigenvalues );