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
63
64
65
66
67
68
69
70
71
72
|
$DEBUG = 0;
my $eps = 1e-8;
######### help funcs
sub ok_matrix ($$$)
{
my ($a, $b, $msg) = @_;
my $res = abs($a-$b);
ok( similar($a,$b) , $msg);
print " (|Delta| = $res)\n" if $DEBUG;
}
sub ok_matrix_orthogonal ($)
{
my ($M) = @_;
my $tmp = $M->shadow();
$tmp->one();
my $transp = $M->shadow();
$transp->transpose($M);
$tmp->subtract($M->multiply($transp), $tmp);
my $v = $tmp->norm_one();
ok(($v < $eps), 'matrix is orthogonal');
print " (|M * ~M - I| = $v)\n" if $DEBUG;
}
sub ok_eigenvectors ($$$;$)
{
my ($M, $L, $V, $msg) = @_;
$msg ||= 'eigenvectors computed correctly';
# Now check that all of them correspond to eigenvalue * eigenvector
my ($rows, $columns) = $M->dim();
unless ($rows == $columns) {
ok(0,'matrix should be square to compute eigenvalues');
return;
}
# Computes the result of all eigenvectors...
my $test = $M * $V;
my $test2 = $V->clone();
for (my $i = 1; $i <= $columns; $i++)
{
my $lambda = $L->element($i,1);
for (my $j = 1; $j <= $rows; $j++)
{ # Compute new vector via lambda * x
$test2->assign($j, $i, $lambda * $test2->element($j, $i));
}
}
ok_matrix($test,$test2, $msg );
return;
}
sub similar($$;$) {
my ($x,$y, $eps) = @_;
$eps ||= 1e-8;
abs($x-$y) < $eps ? 1 : 0;
}
sub _debug_info
{
my($text,$object,$argument,$flag) = @_;
unless (defined $object) { $object = 'undef'; };
unless (defined $argument) { $argument = 'undef'; };
unless (defined $flag) { $flag = 'undef'; };
if (ref($object)) { $object = ref($object); }
if (ref($argument)) { $argument = ref($argument); }
print "$text: \$obj='$object' \$arg='$argument' \$flag='$flag'\n";
}
sub assert_dies($;$)
{
my ($code,$msg) = @_;
eval { &$code };
ok($@, $msg);
}
1;
|
