File Coverage

lib/Math/MatrixDecomposition/LU.pm
Criterion Covered Total %
statement 106 125 84.8
branch 20 38 52.6
condition 10 28 35.7
subroutine 11 13 84.6
pod 5 5 100.0
total 152 209 72.7


line stmt bran cond sub pod time code
1             ## Math/MatrixDecomposition/LU.pm --- LU decomposition.
2              
3             # Copyright (C) 2010 Ralph Schleicher. All rights reserved.
4              
5             # This program is free software; you can redistribute it and/or
6             # modify it under the same terms as Perl itself.
7              
8             ## Code:
9              
10             package Math::MatrixDecomposition::LU;
11              
12 1     1   51945 use strict;
  1         8  
  1         24  
13 1     1   4 use warnings;
  1         1  
  1         18  
14 1     1   5 use Carp;
  1         1  
  1         58  
15 1     1   5 use Exporter qw(import);
  1         1  
  1         31  
16 1     1   5 use Scalar::Util qw(looks_like_number);
  1         1  
  1         38  
17              
18 1     1   548 use Math::BLAS;
  1         13286  
  1         99  
19 1     1   310 use Math::MatrixDecomposition::Util qw(mod min);
  1         1  
  1         64  
20              
21             BEGIN
22             {
23 1     1   3 our $VERSION = '1.06';
24 1         899 our @EXPORT_OK = qw(lu);
25             }
26              
27             # Create a LU decomposition (convenience function).
28             sub lu
29             {
30 0     0 1 0 __PACKAGE__->new (@_);
31             }
32              
33             # Standard constructor.
34             sub new
35             {
36 1     1 1 8 my $class = shift;
37 1         4 my $self =
38             {
39             # LU matrix in row major layout.
40             LU => [],
41              
42             # Number of matrix rows and columns.
43             rows => 0,
44             columns => 0,
45              
46             # Pivot row indices (row permutations).
47             pivot => [],
48              
49             # Sign of determinant. A value of zero means that
50             # the matrix is singular.
51             sign => 0,
52             };
53              
54 1   33     6 bless ($self, ref ($class) || $class);
55              
56             # Process arguments.
57 1 50       3 $self->decompose (@_)
58             if @_ > 0;
59              
60             # Return object.
61 1         2 $self;
62             }
63              
64             # LU decomposition.
65             sub decompose
66             {
67 2     2 1 1006 my $self = shift;
68              
69             # Check arguments.
70 2         3 my $a = shift;
71              
72 2 50 33     9 my $m = @_ > 0 && looks_like_number ($_[0]) ? shift : sqrt (@$a);
73 2 50 33     7 my $n = @_ > 0 && looks_like_number ($_[0]) ? shift : $m ? @$a / $m : 0;
    50          
74              
75 2 50 33     26 croak ('Invalid argument')
      33        
      33        
      33        
76             if (@$a != ($m * $n)
77             || mod ($m, 1) != 0 || $m < 1
78             || mod ($n, 1) != 0 || $n < 1);
79              
80             # Get properties.
81 2         5 my %prop = @_;
82              
83 2   50     7 $prop{capture} //= 0;
84              
85             # Index of last row/column.
86 2         3 my $last_row = $m - 1;
87 2         3 my $last_col = $n - 1;
88              
89             # Matrix elements.
90 2 50       4 if ($prop{capture})
91             {
92             # Work in-place.
93 0         0 $$self{LU} = $a;
94             }
95             else
96             {
97             # Copy matrix elements.
98 2   50     7 $$self{LU} //= [];
99 2         4 @{$$self{LU}} = @$a;
  2         4  
100 2         3 $a = $$self{LU};
101             }
102              
103             # Matrix size.
104 2         2 $$self{rows} = $m;
105 2         3 $$self{columns} = $n;
106              
107             # Pivot row indices (row permutations).
108 2         2 my $pivot = $$self{pivot};
109              
110 2         3 @$pivot = ();
111              
112             # Sign of determinant.
113 2         2 my $sign = 1;
114              
115             # Work variables.
116 2         3 my ($i, $j, $k, $p, $row, $tem);
117              
118 2         6 for $k (0 .. min ($last_row, $last_col))
119             {
120             # Search pivot element.
121 7         93 $p = $k + (blas_amax_val ($m - $k, $a,
122             x_ind => $k * $n + $k,
123             x_incr => $n))[0];
124              
125 7         132 push (@$pivot, $p);
126              
127             # Swap rows.
128 7 100       11 if ($p != $k)
129             {
130 4         11 blas_swap ($n, $a, $a,
131             x_ind => $k * $n,
132             y_ind => $p * $n);
133              
134 4         106 $sign = - $sign;
135             }
136              
137             # Divide by the pivot element.
138 7 50       10 if ($$a[$k * $n + $k] == 0)
139             {
140 0         0 $sign = 0;
141 0         0 next;
142             }
143              
144 7         13 for $i ($k + 1 .. $last_row)
145             {
146 9         100 $$a[$i * $n + $k] /= $$a[$k * $n + $k];
147              
148 9         19 blas_axpby ($last_col - $k, $a, $a,
149             alpha => - $$a[$i * $n + $k],
150             x_ind => $k * $n + $k + 1,
151             y_ind => $i * $n + $k + 1);
152             }
153             }
154              
155             # Save sign of determinant.
156 2         3 $$self{sign} = $sign;
157              
158             # Return object.
159 2         4 $self;
160             }
161              
162             # Solve a system of linear equations.
163             sub solve
164             {
165 2     2 1 8 my $self = shift;
166              
167 2         2 my $b = shift;
168 2         3 my $x = shift;
169              
170             # LU matrix elements.
171 2         2 my $a = $$self{LU};
172              
173             # Number of matrix rows and columns.
174 2         3 my $m = $$self{rows};
175 2         3 my $n = $$self{columns};
176              
177 2 50 33     10 croak ('Invalid argument')
178             if @$b == 0 || @$b % $m != 0;
179              
180 2         33 my $l = @$b / $m;
181              
182             # Copy right-hand side.
183 2 50       6 if (defined ($x))
184             {
185 0         0 @$x = @$b;
186             }
187             else
188             {
189             # Work in-place.
190 2         2 $x = $b;
191             }
192              
193             # Index of last row/column.
194 2         10 my $last_row = $m - 1;
195 2         3 my $last_col = $n - 1;
196 2         3 my $last_vec = $l - 1;
197              
198             # Pivot row indices.
199 2         4 my $pivot = $$self{pivot};
200              
201             # Work variables.
202 2         4 my ($i, $j, $k, $p, $row, $tem);
203              
204 2 100       5 if ($l == 1)
205             {
206             # Apply row permutations.
207 1         3 for $i (0 .. $#$pivot)
208             {
209 4         5 $p = $$pivot[$i];
210 4 100       6 next if $p == $i;
211              
212 2         5 @$x[$i, $p] = @$x[$p, $i];
213             }
214              
215             # Forward substitution, that is solve 'Ly = b' for 'y'.
216             # Elements of 'y' are saved in 'x'.
217 1         2 for $i (1 .. $last_row)
218             {
219 3         4 $row = $i * $n;
220              
221 3         3 for $j (0 .. $i - 1)
222             {
223 6         9 $$x[$i] -= $$x[$j] * $$a[$row + $j];
224             }
225             }
226              
227             # Backward substitution, that is solve 'Ux = y' for 'x'.
228 1         3 for $i (reverse (0 .. $last_row))
229             {
230 4         13 $row = $i * $n;
231              
232 4 50       9 if ($$a[$row + $i] == 0)
233             {
234             # Matrix A is singular.
235 0 0       0 return unless $$x[$i] == 0;
236             }
237             else
238             {
239 4         6 for $j ($i + 1 .. $last_col)
240             {
241 6         7 $$x[$i] -= $$x[$j] * $$a[$row + $j];
242             }
243              
244 4         5 $$x[$i] /= $$a[$row + $i];
245             }
246             }
247             }
248             else
249             {
250             # Matrix A is singular.
251 1 50       4 return if $$self{sign} == 0;
252              
253             # Apply row permutations.
254 1         4 for $i (0 .. $#$pivot)
255             {
256 3         48 $p = $$pivot[$i];
257 3 100       13 next if $p == $i;
258              
259 2         8 blas_swap ($l, $x, $x,
260             x_ind => $i * $l,
261             y_ind => $p * $l);
262             }
263              
264             # Forward substitution.
265 1         3 for $k (0 .. $last_col)
266             {
267 3         41 for $i ($k + 1 .. $last_col)
268             {
269 3         27 blas_axpby ($l, $x, $x,
270             alpha => - $$a[$i * $n + $k],
271             x_ind => $k * $l,
272             y_ind => $i * $l);
273             }
274             }
275              
276             # Backward substitution.
277 1         3 for $k (reverse (0 .. $last_col))
278             {
279 3         44 blas_rscale ($l, $x,
280             alpha => $$a[$k * $n + $k],
281             x_ind => $k * $l);
282              
283 3         111 for $i (0 .. $k - 1)
284             {
285 3         30 blas_axpby ($l, $x, $x,
286             alpha => - $$a[$i * $n + $k],
287             x_ind => $k * $l,
288             y_ind => $i * $l);
289             }
290             }
291             }
292              
293             # Resize result.
294 2 50       4 splice (@$x, $n * $l)
295             if $m > $n;
296              
297             # Fix rounding errors.
298 2         3 for (@$x)
299             {
300 13         32 $_ = 0 + "$_";
301             }
302              
303             # Return value.
304 2         5 $x;
305             }
306              
307             # Determinant.
308             sub det
309             {
310 0     0 1   my $self = shift;
311              
312 0           my $m = $$self{rows};
313 0           my $n = $$self{columns};
314              
315 0 0         croak ('Matrix has to be square')
316             if $m != $n;
317              
318             # Calculate determinant.
319 0           my $det = $$self{sign};
320              
321 0 0         if ($det)
322             {
323 0           my $u = $$self{LU};
324              
325 0           my $ind = 0;
326 0           my $incr = $n + 1;
327              
328 0           for (1 .. $n)
329             {
330 0           $det *= $$u[$ind];
331 0           $ind += $incr;
332             }
333             }
334              
335 0           $det;
336             }
337              
338             1;
339              
340             __END__