1#!/usr/bin/perl
2######################### We start with some black magic to print on failure.
3use lib '../blib/lib','../blib/arch';
4use strict;
5use warnings;
6use vars qw($loaded);
7
8BEGIN {$| = 1; print "1..74\n";}
9END {print "not ok 1\n" unless $loaded;}
10use Math::Cephes::Matrix qw(mat);
11$loaded = 1;
12print "ok 1\n";
13
14######################### End of black magic.
15
16# util
17my $count = 1;
18my $eps = 1e-07;
19sub ok {
20  local($^W) = 0;
21  $count++;
22  my ($package, $file, $line) = caller;
23  my ($value, $true, $skip) = @_;
24  $skip ||= '';
25  $skip = "# skip ($skip)" if $skip;
26  my $error = sprintf( "%12.8f", abs($value - $true));
27  print($error < $eps ? "ok $count $skip\n" :
28	"not ok $count (expected $true: got $value) at $file line $line\n");
29}
30
31my $M = Math::Cephes::Matrix->new([ [1, 2, -1], [2, -3, 1], [1, 0, 3]]);
32my $B = [2, -1, 10];
33my $X = $M->simq($B);
34ok( $X->[0], 1);
35ok( $X->[1], 2);
36ok( $X->[2], 3);
37my $C = Math::Cephes::Matrix->new([ [1, 2, 4], [2, 9, 2], [6, 2, 7]]);
38my $I = $C->inv();
39my $T = $I->mul($C)->coef;
40ok( $T->[0]->[0], 1);
41ok( $T->[1]->[1], 1);
42ok( $T->[2]->[2], 1);
43ok( $T->[0]->[1], 0);
44ok( $T->[1]->[0], 0);
45ok( $T->[2]->[0], 0);
46my $V = $M->mul($X);
47ok( $V->[0], $B->[0]);
48ok( $V->[1], $B->[1]);
49ok( $V->[2], $B->[2]);
50my $D = $M->add($C)->coef;
51ok( $D->[0]->[0], 2);
52ok( $D->[1]->[1], 6);
53ok( $D->[2]->[2], 10);
54ok( $D->[0]->[1], 4);
55ok( $D->[1]->[0], 4);
56ok( $D->[2]->[0], 7);
57$D = $M->sub($C)->coef;
58ok( $D->[0]->[0], 0);
59ok( $D->[1]->[1], -12);
60ok( $D->[2]->[2], -4);
61ok( $D->[0]->[1], 0);
62ok( $D->[1]->[0], 0);
63ok( $D->[2]->[0], -5);
64my $H = $C->transp()->coef;
65ok( $H->[0]->[0], 1);
66ok( $H->[1]->[1], 9);
67ok( $H->[2]->[2], 7);
68ok( $H->[0]->[1], 2);
69ok( $H->[1]->[0], 2);
70ok( $H->[2]->[0], 4);
71my $R = $M->div($C);
72my $Q = $R->mul($C)->coef;
73my $Mc = $M->coef;
74for (my $i=0; $i<3; $i++) {
75    for (my $j=0; $j<3; $j++) {
76	ok($Q->[$i]->[$j], $Mc->[$i]->[$j]);
77    }
78}
79$R = $M->mul($C)->coef;
80ok( $R->[0]->[0], -1);
81ok( $R->[1]->[1], -21);
82ok( $R->[2]->[2], 25);
83ok( $R->[0]->[1], 18);
84ok( $R->[1]->[0], 2);
85ok( $R->[2]->[0], 19);
86
87$C->clr();
88$R = $C->coef;
89ok( $R->[0]->[0], 0);
90ok( $R->[2]->[2], 0);
91ok( $R->[1]->[0], 0);
92ok( $R->[2]->[0], 0);
93
94$C->clr(3);
95$R = $C->coef;
96ok( $R->[0]->[0], 3);
97ok( $R->[2]->[2], 3);
98ok( $R->[1]->[0], 3);
99ok( $R->[2]->[0], 3);
100
101my $S = Math::Cephes::Matrix->new([ [1, 2, 3], [2, 2, 3], [3, 3, 4]]);
102my ($E, $EV1) = $S->eigens();
103my $EV = $EV1->coef;
104for (my $i=0; $i<3; $i++) {
105  my $v = [];
106  for (my $j=0; $j<3; $j++) {
107    $v->[$j] = $EV->[$i]->[$j];
108  }
109  my $sv = $S->mul($v);
110  for (my $j=0; $j<3; $j++) {
111    ok($sv->[$j], $E->[$i]*$v->[$j]);
112  }
113}
114
115my $Z = $M->new()->coef;
116for (my $i=0; $i<3; $i++) {
117    for (my $j=0; $j<3; $j++) {
118	ok($Z->[$i]->[$j], $Mc->[$i]->[$j]);
119    }
120}
121$Z->[0]->[0] = 5;
122ok($Mc->[0]->[0], 1);
123ok($Z->[0]->[0], 5);
124
125