| 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 | |