Commit | Line | Data |
---|---|---|
86530b38 AT |
1 | $DEBUG = 0; |
2 | ||
3 | ######### help funcs | |
4 | ||
5 | sub ok ($$) { | |
6 | my($number, $result) = @_ ; | |
7 | ||
8 | print "ok $number\n" if $result ; | |
9 | print "not ok $number\n" if !$result ; | |
10 | } | |
11 | ||
12 | sub random_matrix ($) | |
13 | { | |
14 | my ($size) = @_; | |
15 | my $M = Math::MatrixReal->new($size, $size); | |
16 | for (my $i=1; $i<=$size; $i++) | |
17 | { | |
18 | for (my $j=1; $j<=$size; $j++) | |
19 | { | |
20 | $M->assign($i,$j,rand()); | |
21 | } | |
22 | } | |
23 | return $M; | |
24 | } | |
25 | ||
26 | sub ok_matrix ($$$) | |
27 | { | |
28 | my ($no, $M1, $M2) = @_; | |
29 | my $tmp = $M1->shadow(); | |
30 | $tmp->subtract($M1,$M2); | |
31 | my $v = $tmp->norm_one(); | |
32 | ok($no, ($v < 1e-8)); | |
33 | print " ($no: |Delta| = $v)\n" if $DEBUG; | |
34 | } | |
35 | ||
36 | sub ok_matrix_orthogonal ($$) | |
37 | { | |
38 | my ($no, $M) = @_; | |
39 | my $tmp = $M->shadow(); | |
40 | $tmp->one(); | |
41 | my $transp = $M->shadow(); | |
42 | $transp->transpose($M); | |
43 | $tmp->subtract($M->multiply($transp), $tmp); | |
44 | my $v = $tmp->norm_one(); | |
45 | ok($no, ($v < 1e-8)); | |
46 | print " ($no: |M * ~M - I| = $v)\n" if $DEBUG; | |
47 | } | |
48 | ||
49 | sub ok_eigenvectors ($$$$) | |
50 | { | |
51 | my ($no, $M, $L, $V) = @_; | |
52 | # Now check that all of them correspond to eigenvalue * eigenvector | |
53 | my ($rows, $columns) = $M->dim(); | |
54 | unless ($rows == $columns) { | |
55 | ok("$no", 0); | |
56 | return; | |
57 | } | |
58 | # Computes the result of all eigenvectors... | |
59 | my $test = $M * $V; | |
60 | my $test2 = $V->clone(); | |
61 | for (my $i = 1; $i <= $columns; $i++) | |
62 | { | |
63 | my $lambda = $L->element($i,1); | |
64 | for (my $j = 1; $j <= $rows; $j++) | |
65 | { # Compute new vector via lambda * x | |
66 | $test2->assign($j, $i, $lambda * $test2->element($j, $i)); | |
67 | } | |
68 | } | |
69 | ok_matrix("$no",$test,$test2); | |
70 | return; | |
71 | } | |
72 | sub similar($$$) { | |
73 | my ($x,$y) = @_; | |
74 | my $eps = shift || 1e-8; | |
75 | abs($x - $y) < $eps ? return 1 : return 0; | |
76 | } | |
77 | 1; | |
78 |