File Coverage

blib/lib/Math/Cephes/Matrix.pm
Criterion Covered Total %
statement 139 147 94.5
branch 19 34 55.8
condition 7 11 63.6
subroutine 17 18 94.4
pod 1 15 6.6
total 183 225 81.3


line stmt bran cond sub pod time code
1             package Math::Cephes::Matrix;
2 1     1   4259 use strict;
  1         3  
  1         32  
3 1     1   5 use warnings;
  1         1  
  1         27  
4 1     1   5 use vars qw(@EXPORT_OK $VERSION);
  1         1  
  1         2082  
5              
6             require Exporter;
7             *import = \&Exporter::import;
8             @EXPORT_OK = qw(mat);
9              
10             $VERSION = '0.5304';
11              
12             require Math::Cephes;
13              
14             sub new {
15 15     15 0 2182 my ($caller, $arr) = @_;
16 15         27 my $refer = ref($caller);
17 15   66     62 my $class = $refer || $caller;
18 15 50 66     59 die "Must supply data for the matrix"
19             unless ($refer or $arr);
20 15 100       32 unless ($refer) {
21 14 50 33     66 die "Please supply an array of arrays for the matrix data"
22             unless (ref($arr) eq 'ARRAY' and ref($arr->[0]) eq 'ARRAY');
23 14         14 my $n = scalar @$arr;
24 14         14 my $m = scalar @{$arr->[0]};
  14         18  
25 14 50       27 die "Matrices must be square" unless $m == $n;
26             }
27 15         16 my ($coef, $n);
28 15 100       28 if ($refer) {
29 1         4 $n = $caller->{n};
30 1         3 my $cdata = $caller->{coef};
31 1         21 foreach (@$cdata) {
32 3         13 push @$coef, [ @$_];
33             }
34             }
35             else {
36 14         23 ($coef, $n) = ($arr, scalar @$arr);
37             }
38 15         104 bless { coef => $coef,
39             n => $n,
40             }, $class;
41             }
42              
43             sub mat {
44 0     0 0 0 return Math::Cephes::Matrix->new(shift);
45             }
46              
47             sub mat_to_vec {
48 4     4 0 6 my $self = shift;
49 4         8 my ($M, $n) = ($self->{coef}, $self->{n});
50 4         5 my $A = [];
51 4         13 for (my $i=0; $i<$n; $i++) {
52 12         33 for (my $j=0; $j<$n; $j++) {
53 36         35 my $index = $i*$n+$j;
54 36         90 $A->[$index] = $M->[$i]->[$j];
55             }
56             }
57 4         9 return $A;
58             }
59              
60             sub vec_to_mat {
61 4     4 0 6 my ($self, $X) = @_;
62 4         9 my $n = $self->{n};
63 4         6 my $I = [];
64 4         13 for (my $i=0; $i<$n; $i++) {
65 12         25 for (my $j=0; $j<$n; $j++) {
66 36         33 my $index = $i*$n+$j;
67 36         99 $I->[$i]->[$j] = $X->[$index];
68             }
69             }
70 4         6 return $I;
71             }
72              
73             sub check {
74 12     12 0 20 my ($self, $B) = @_;
75 12         24 my $na = $self->{n};
76 12         22 my $ref = ref($B);
77 12 100       75 if ($ref eq 'Math::Cephes::Matrix') {
    50          
78 7 50       18 die "Matrices must be of the same size"
79             unless $B->{n} == $na;
80 7         17 return $B->coef;
81             }
82             elsif ($ref eq 'ARRAY') {
83 5         6 my $nb = scalar @$B;
84 5         9 my $ref0 = ref($B->[0]);
85 5 50       13 if ($ref0 eq 'ARRAY') {
    50          
86 0         0 my $m = scalar @{$B->[0]};
  0         0  
87 0 0       0 die "Can only use square matrices" unless $m == $nb;
88 0 0       0 die "Can only use matrices of the same size"
89             unless $na == $nb;
90 0         0 return $B;
91             }
92             elsif (not $ref0) {
93 5 50       525 die "Can only use vectors of the same size" unless $nb == $na;
94 5         13 return $B;
95             }
96             else {
97 0         0 die "Unknown reference '$ref0' for data";
98             }
99             }
100             else {
101 0         0 die "Unknown reference '$ref' for data";
102             }
103             }
104              
105             sub coef {
106 18     18 0 59 return $_[0]->{coef};
107             }
108              
109             sub clr {
110 2     2 0 1740 my $self = shift;
111 2   100     26 my $what = shift || 0;
112 2         5 my $n = $self->{n};
113 2         3 my $B = [];
114 2         7 for (my $i=0; $i<$n; $i++) {
115 6         13 for (my $j=0; $j<$n; $j++) {
116 18         271 $B->[$i]->[$j] = $what;
117             }
118             }
119 2         7 $self->{coef} = $B;
120             }
121              
122             sub simq {
123 1     1 0 7 my ($self, $B) = @_;
124 1         4 $B = $self->check($B);
125 1         3 my ($M, $n) = ($self->{coef}, $self->{n});
126 1 50       4 die "Must supply an array reference for B" unless ref($B) eq 'ARRAY';
127 1         3 my $A = $self->mat_to_vec();
128 1         12 my $X = [split //, 0 x $n];
129 1         4 my $IPS = [split //, 0 x $n];
130 1         2 my $flag = 0;
131 1         38 my $ret = Math::Cephes::simq($A, $B, $X, $n, $flag, $IPS);
132 1 50       5 return $ret ? undef : $X;
133             }
134              
135              
136             sub inv {
137 2     2 0 9 my $self = shift;
138 2         6 my ($M, $n) = ($self->{coef}, $self->{n});
139 2         7 my $A = $self->mat_to_vec();
140 2         16 my $X = [split //, 0 x ($n*$n)];
141 2         9 my $B = [split //, 0 x $n];
142 2         8 my $IPS = [split //, 0 x $n];
143 2         5 my $flag = 0;
144 2         91 my $ret = Math::Cephes::minv($A, $X, $n, $B, $IPS);
145 2 50       10 return undef if $ret;
146 2         7 my $I = $self->vec_to_mat($X);
147 2         7 return Math::Cephes::Matrix->new($I);
148             }
149              
150             sub transp {
151 1     1 0 827 my $self = shift;
152 1         3 my ($M, $n) = ($self->{coef}, $self->{n});
153 1         4 my $A = $self->mat_to_vec();
154 1         8 my $T = [split //, 0 x ($n*$n)];
155 1         13 Math::Cephes::mtransp($n, $A, $T);
156 1         4 my $R = $self->vec_to_mat($T);
157 1         4 return Math::Cephes::Matrix->new($R);
158             }
159              
160             sub add {
161 1     1 1 468 my ($self, $B) = @_;
162 1         3 $B = $self->check($B);
163 1         4 my ($A, $n) = ($self->{coef}, $self->{n});
164 1         2 my $C = [];
165 1         4 for (my $i=0; $i<$n; $i++) {
166 3         5 for (my $j=0; $j<$n; $j++) {
167 9         26 $C->[$i]->[$j] = $A->[$i]->[$j] + $B->[$i]->[$j];
168             }
169             }
170 1         4 return Math::Cephes::Matrix->new($C);
171             }
172              
173             sub sub {
174 1     1 0 878 my ($self, $B) = @_;
175 1         3 $B = $self->check($B);
176 1         13 my ($A, $n) = ($self->{coef}, $self->{n});
177 1         2 my $C = [];
178 1         5 for (my $i=0; $i<$n; $i++) {
179 3         5 for (my $j=0; $j<$n; $j++) {
180 9         25 $C->[$i]->[$j] = $A->[$i]->[$j] - $B->[$i]->[$j];
181             }
182             }
183 1         4 return Math::Cephes::Matrix->new($C);
184             }
185              
186             sub mul {
187 8     8 0 3270 my ($self, $B) = @_;
188 8         19 $B = $self->check($B);
189 8         19 my ($A, $n) = ($self->{coef}, $self->{n});
190 8         23 my $C = [];
191 8 100       28 if (ref($B->[0]) eq 'ARRAY') {
192 4         12 for (my $i=0; $i<$n; $i++) {
193 12         33 for (my $j=0; $j<$n; $j++) {
194 36         57 for (my $m=0; $m<$n; $m++) {
195 108         382 $C->[$i]->[$j] += $A->[$i]->[$m] * $B->[$m]->[$j];
196             }
197             }
198             }
199 4         13 return Math::Cephes::Matrix->new($C);
200             }
201             else {
202 4         20 for (my $i=0; $i<$n; $i++) {
203 12         26 for (my $m=0; $m<$n; $m++) {
204 36         152 $C->[$i] += $A->[$i]->[$m] * $B->[$m];
205             }
206             }
207 4         12 return $C;
208             }
209             }
210              
211             sub div {
212 1     1 0 811 my ($self, $B) = @_;
213 1         10 $B = $self->check($B);
214 1         3 my $C = Math::Cephes::Matrix->new($B)->inv();
215 1         4 my $D = $self->mul($C);
216 1         5 return $D;
217             }
218              
219             sub eigens {
220 1     1 0 5 my $self = shift;
221 1         3 my ($M, $n) = ($self->{coef}, $self->{n});
222 1         2 my $A = [];
223 1         5 for (my $i=0; $i<$n; $i++) {
224 3         5 for (my $j=0; $j<$n; $j++) {
225 9         12 my $index = ($i*$i+$i)/2 + $j;
226 9         21 $A->[$index] = $M->[$i]->[$j];
227             }
228             }
229 1         19 my $EV1 = [split //, 0 x ($n*$n)];
230 1         4 my $E = [split //, 0 x $n];
231 1         4 my $IPS = [split //, 0 x $n];
232 1         36 Math::Cephes::eigens($A, $EV1, $E, $n);
233 1         5 my $EV = $self->vec_to_mat($EV1);
234 1         6 return ($E, Math::Cephes::Matrix->new($EV));
235             }
236              
237             1;
238              
239             __END__