File Coverage

blib/lib/Math/Matrix.pm
Criterion Covered Total %
statement 2055 2194 93.6
branch 919 1472 62.4
condition 120 255 47.0
subroutine 228 233 97.8
pod 187 187 100.0
total 3509 4341 80.8


line stmt bran cond sub pod time code
1             # -*- mode: perl; coding: utf-8-unix -*-
2              
3             package Math::Matrix;
4              
5 124     124   1274448 use strict;
  124         973  
  124         3588  
6 124     124   730 use warnings;
  124         236  
  124         2996  
7              
8 124     124   632 use Carp;
  124         250  
  124         12830  
9 124     124   829 use Scalar::Util 'blessed';
  124         295  
  124         63969  
10              
11             our $VERSION = '0.94';
12             our $eps = 0.00001;
13              
14             use overload
15              
16             '+' => sub {
17 3     3   19 my ($x, $y, $swap) = @_;
18 3         9 $x -> add($y);
19             },
20              
21             '-' => sub {
22 6     6   40 my ($x, $y, $swap) = @_;
23 6 100       28 if ($swap) {
24 3 100 66     40 return $x -> neg() if !ref($y) && $y == 0;
25              
26 1         4 my $class = ref $x;
27 1         4 return $class -> new($y) -> sub($x);
28             }
29 3         8 $x -> sub($y);
30             },
31              
32             '*' => sub {
33 23     23   112 my ($x, $y, $swap) = @_;
34 23         51 $x -> mul($y);
35             },
36              
37             '**' => sub {
38 3     3   21 my ($x, $y, $swap) = @_;
39 3 100       14 if ($swap) {
40 1         3 my $class = ref $x;
41 1         5 return $class -> new($y) -> pow($x);
42             }
43 2         6 $x -> pow($y);
44             },
45              
46             '==' => sub {
47 2     2   19 my ($x, $y, $swap) = @_;
48 2         7 $x -> meq($y);
49             },
50              
51             '!=' => sub {
52 2     2   5661 my ($x, $y, $swap) = @_;
53 2         8 $x -> mne($y);
54             },
55              
56             'int' => sub {
57 2     2   12 my ($x, $y, $swap) = @_;
58 2         6 $x -> int();
59             },
60              
61             'abs' => sub {
62 1     1   8 my ($x, $y, $swap) = @_;
63 1         4 $x -> abs();
64             },
65              
66 124         2860 '~' => 'transpose',
67             '""' => 'as_string',
68 124     124   157010 '=' => 'clone';
  124         133553  
69              
70             =pod
71              
72             =encoding utf8
73              
74             =head1 NAME
75              
76             Math::Matrix - multiply and invert matrices
77              
78             =head1 SYNOPSIS
79              
80             use Math::Matrix;
81              
82             # Generate a random 3-by-3 matrix.
83             srand(time);
84             my $A = Math::Matrix -> new([rand, rand, rand],
85             [rand, rand, rand],
86             [rand, rand, rand]);
87             $A -> print("A\n");
88              
89             # Append a fourth column to $A.
90             my $x = Math::Matrix -> new([rand, rand, rand]);
91             my $E = $A -> concat($x -> transpose);
92             $E -> print("Equation system\n");
93              
94             # Compute the solution.
95             my $s = $E -> solve;
96             $s -> print("Solutions s\n");
97              
98             # Verify that the solution equals $x.
99             $A -> multiply($s) -> print("A*s\n");
100              
101             =head1 DESCRIPTION
102              
103             This module implements various constructors and methods for creating and
104             manipulating matrices.
105              
106             All methods return new objects, so, for example, C<$X-Eadd($Y)> does not
107             modify C<$X>.
108              
109             $X -> add($Y); # $X not modified; output is lost
110             $X = $X -> add($Y); # this works
111              
112             Some operators are overloaded (see L) and allow the operand to be
113             modified directly.
114              
115             $X = $X + $Y; # this works
116             $X += $Y; # so does this
117              
118             =head1 METHODS
119              
120             =head2 Constructors
121              
122             =over 4
123              
124             =item new()
125              
126             Creates a new object from the input arguments and returns it.
127              
128             If a single input argument is given, and that argument is a reference to array
129             whose first element is itself a reference to an array, it is assumed that the
130             argument contains the whole matrix, like this:
131              
132             $x = Math::Matrix->new([[1, 2, 3], [4, 5, 6]]); # 2-by-3 matrix
133             $x = Math::Matrix->new([[1, 2, 3]]); # 1-by-3 matrix
134             $x = Math::Matrix->new([[1], [2], [3]]); # 3-by-1 matrix
135              
136             If a single input argument is given, and that argument is not a reference to an
137             array, a 1-by-1 matrix is returned.
138              
139             $x = Math::Matrix->new(1); # 1-by-1 matrix
140              
141             Note that all the folling cases result in an empty matrix:
142              
143             $x = Math::Matrix->new([[], [], []]);
144             $x = Math::Matrix->new([[]]);
145             $x = Math::Matrix->new([]);
146              
147             If C> is called as an instance method with no input arguments, a zero
148             filled matrix with identical dimensions is returned:
149              
150             $b = $a->new(); # $b is a zero matrix with the size of $a
151              
152             Each row must contain the same number of elements.
153              
154             =cut
155              
156             sub new {
157 947     947 1 964889 my $that = shift;
158 947   66     4060 my $class = ref($that) || $that;
159 947         1607 my $self = [];
160              
161             # If called as an instance method and no arguments are given, return a
162             # zero matrix of the same size as the invocand.
163              
164 947 100 100     2841 if (ref($that) && (@_ == 0)) {
165 11         81 @$self = map [ (0) x @$_ ], @$that;
166             }
167              
168             # Otherwise return a new matrix based on the input arguments. The object
169             # data is a blessed reference to an array containing the matrix data. This
170             # array contains a list of arrays, one for each row, which in turn contains
171             # a list of elements. An empty matrix has no rows.
172             #
173             # [[ 1, 2, 3 ], [ 4, 5, 6 ]] 2-by-3 matrix
174             # [[ 1, 2, 3 ]] 1-by-3 matrix
175             # [[ 1 ], [ 2 ], [ 3 ]] 3-by-1 matrix
176             # [[ 1 ]] 1-by-1 matrix
177             # [] empty matrix
178              
179             else {
180              
181 936         1381 my $data;
182              
183             # If there is a single argument, and that is not a reference,
184             # assume new() has been called as, e.g., $class -> new(3).
185              
186 936 100 100     6244 if (@_ == 1 && !ref($_[0])) {
    100 66        
      100        
      100        
187 16         34 $data = [[ $_[0] ]];
188             }
189              
190             # If there is a single argument, and that is a reference to an array,
191             # and that array contains at least one element, and that element is
192             # itself a reference to an array, then assume new() has been called
193             # with the matrix as one argument, i.e., a reference to an array of
194             # arrays, e.g., $class -> new([ [1, 2], [3, 4] ]) ...
195              
196             elsif (@_ == 1 && ref($_[0]) eq 'ARRAY'
197 908         3936 && @{$_[0]} > 0 && ref($_[0][0]) eq 'ARRAY')
198             {
199 643         1055 $data = $_[0];
200             }
201              
202             # ... otherwise assume that each argument to new() is a row. Note that
203             # new() called with no arguments results in an empty matrix.
204              
205             else {
206 277         619 $data = [ @_ ];
207             }
208              
209             # Sanity checking.
210              
211 936 100       2192 if (@$data) {
212 935         1566 my $nrow = @$data;
213 935         1205 my $ncol;
214              
215 935         2378 for my $i (0 .. $nrow - 1) {
216 1828         2640 my $row = $data -> [$i];
217              
218             # Verify that the row is a reference to an array.
219              
220 1828 50       3514 croak "row with index $i is not a reference to an array"
221             unless ref($row) eq 'ARRAY';
222              
223             # In the first round, get the number of elements, i.e., the
224             # number of columns in the matrix. In the successive
225             # rounds, verify that each row has the same number of
226             # elements.
227              
228 1828 100       3106 if ($i == 0) {
229 935         1776 $ncol = @$row;
230             } else {
231 893 50       1961 croak "each row must have the same number of elements"
232             unless @$row == $ncol;
233             }
234             }
235              
236             # Copy the data into $self only if the matrix is non-emtpy.
237              
238 935 100       4200 @$self = map [ @$_ ], @$data if $ncol;
239             }
240             }
241              
242 947         2802 bless $self, $class;
243             }
244              
245             =pod
246              
247             =item new_from_sub()
248              
249             Creates a new matrix object by doing a subroutine call to create each element.
250              
251             $sub = sub { ... };
252             $x = Math::Matrix -> new_from_sub($sub); # 1-by-1
253             $x = Math::Matrix -> new_from_sub($sub, $m); # $m-by-$m
254             $x = Math::Matrix -> new_from_sub($sub, $m, $n); # $m-by-$n
255              
256             The subroutine is called in scalar context with two input arguments, the row and
257             column indices of the element to be created. Note that no checks are performed
258             on the output of the subroutine.
259              
260             Example 1, a 4-by-4 identity matrix can be created with
261              
262             $sub = sub { $_[0] == $_[1] ? 1 : 0 };
263             $x = Math::Matrix -> new_from_sub($sub, 4);
264              
265             Example 2, the code
266              
267             $x = Math::Matrix -> new_from_sub(sub { 2**$_[1] }, 1, 11);
268              
269             creates the following 1-by-11 vector with powers of two
270              
271             [ 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 ]
272              
273             Example 3, the code, using C<$i> and C<$j> for increased readability
274              
275             $sub = sub {
276             ($i, $j) = @_;
277             $d = $j - $i;
278             return $d == -1 ? 5
279             : $d == 0 ? 6
280             : $d == 1 ? 7
281             : 0;
282             };
283             $x = Math::Matrix -> new_from_sub($sub, 5);
284              
285             creates the tridiagonal matrix
286              
287             [ 6 7 0 0 0 ]
288             [ 5 6 7 0 0 ]
289             [ 0 5 6 7 0 ]
290             [ 0 0 5 6 7 ]
291             [ 0 0 0 5 6 ]
292              
293             =cut
294              
295             sub new_from_sub {
296 3 50   3 1 3475 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
297 3 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 4;
298 3         6 my $class = shift;
299              
300 3 50       7 croak +(caller(0))[3], " is a class method, not an instance method"
301             if ref $class;
302              
303 3         6 my $sub = shift;
304 3 50       8 croak "The first input argument must be a code reference"
305             unless ref($sub) eq 'CODE';
306              
307 3 100       15 my ($nrow, $ncol) = @_ == 0 ? (1, 1)
    50          
308             : @_ == 1 ? (@_, @_)
309             : (@_);
310              
311 3         8 my $x = bless [], $class;
312 3         8 for my $i (0 .. $nrow - 1) {
313 10         39 for my $j (0 .. $ncol - 1) {
314 52         229 $x -> [$i][$j] = $sub -> ($i, $j);
315             }
316             }
317              
318 3         22 return $x;
319             }
320              
321             =pod
322              
323             =item new_from_rows()
324              
325             Creates a new matrix by assuming each argument is a row vector.
326              
327             $x = Math::Matrix -> new_from_rows($y, $z, ...);
328              
329             For example
330              
331             $x = Math::Matrix -> new_from_rows([1, 2, 3],[4, 5, 6]);
332              
333             returns the matrix
334              
335             [ 1 2 3 ]
336             [ 4 5 6 ]
337              
338             =cut
339              
340             sub new_from_rows {
341 2 50   2 1 103 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
342 2         5 my $class = shift;
343              
344 2 50       7 croak +(caller(0))[3], " is a class method, not an instance method"
345             if ref $class;
346              
347 2         5 my @args = ();
348 2         10 for (my $i = 0 ; $i <= $#_ ; ++$i) {
349 8         13 my $x = $_[$i];
350 8 50 33     48 $x = $class -> new($x)
351             unless defined(blessed($x)) && $x -> isa($class);
352 8 100       20 if ($x -> is_vector()) {
353 4         14 push @args, $x -> to_row();
354             } else {
355 4         12 push @args, $x;
356             }
357             }
358              
359 2         5 $class -> new([]) -> catrow(@args);
360             }
361              
362             =pod
363              
364             =item new_from_cols()
365              
366             Creates a matrix by assuming each argument is a column vector.
367              
368             $x = Math::Matrix -> new_from_cols($y, $z, ...);
369              
370             For example,
371              
372             $x = Math::Matrix -> new_from_cols([1, 2, 3],[4, 5, 6]);
373              
374             returns the matrix
375              
376             [ 1 4 ]
377             [ 2 5 ]
378             [ 3 6 ]
379              
380             =cut
381              
382             sub new_from_cols {
383 1 50   1 1 105 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
384 1         3 my $class = shift;
385              
386 1 50       3 croak +(caller(0))[3], " is a class method, not an instance method"
387             if ref $class;
388              
389 1         5 $class -> new_from_rows(@_) -> swaprc();
390             }
391              
392             =pod
393              
394             =item id()
395              
396             Returns a new identity matrix.
397              
398             $I = Math::Matrix -> id($n); # $n-by-$n identity matrix
399             $I = $x -> id($n); # $n-by-$n identity matrix
400             $I = $x -> id(); # identity matrix with size of $x
401              
402             =cut
403              
404             sub id {
405 91     91 1 10516 my $self = shift;
406 91         150 my $ref = ref $self;
407 91   66     294 my $class = $ref || $self;
408              
409 91         128 my $n;
410 91 100       170 if (@_) {
411 88         125 $n = shift;
412             } else {
413 3 50       8 if ($ref) {
414 3         9 my ($mx, $nx) = $self -> size();
415 3 50       13 croak "When id() is called as an instance method, the invocand",
416             " must be a square matrix" unless $mx == $nx;
417 3         6 $n = $mx;
418             } else {
419 0         0 croak "When id() is called as a class method, the size must be",
420             " given as an input argument";
421             }
422             }
423              
424 91         561 bless [ map [ (0) x ($_ - 1), 1, (0) x ($n - $_) ], 1 .. $n ], $class;
425             }
426              
427             =pod
428              
429             =item new_identity()
430              
431             This is an alias for C>.
432              
433             =cut
434              
435             sub new_identity {
436 14     14 1 11390 id(@_);
437             }
438              
439             =pod
440              
441             =item eye()
442              
443             This is an alias for C>.
444              
445             =cut
446              
447             sub eye {
448 6     6 1 10607 new_identity(@_);
449             }
450              
451             =pod
452              
453             =item exchg()
454              
455             Exchange matrix.
456              
457             $x = Math::Matrix -> exchg($n); # $n-by-$n exchange matrix
458              
459             =cut
460              
461             sub exchg {
462 5 50   5 1 6983 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
463 5 50       13 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
464 5         10 my $class = shift;
465              
466 5         7 my $n = shift;
467 5         41 bless [ map [ (0) x ($n - $_), 1, (0) x ($_ - 1) ], 1 .. $n ], $class;
468             }
469              
470             =pod
471              
472             =item scalar()
473              
474             Returns a scalar matrix, i.e., a diagonal matrix with all the diagonal elements
475             set to the same value.
476              
477             # Create an $m-by-$m scalar matrix where each element is $c.
478             $x = Math::Matrix -> scalar($c, $m);
479              
480             # Create an $m-by-$n scalar matrix where each element is $c.
481             $x = Math::Matrix -> scalar($c, $m, $n);
482              
483             Multiplying a matrix A by a scalar matrix B is effectively the same as multiply
484             each element in A by the constant on the diagonal of B.
485              
486             =cut
487              
488             sub scalar {
489 6 50   6 1 6768 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
490 6 50       18 croak "Too many arguments for ", (caller(0))[3] if @_ > 4;
491 6         8 my $class = shift;
492              
493 6 50       13 croak +(caller(0))[3], " is a class method, not an instance method"
494             if ref $class;
495              
496 6         11 my $c = shift;
497 6 100       23 my ($m, $n) = @_ == 0 ? (1, 1)
    100          
498             : @_ == 1 ? (@_, @_)
499             : (@_);
500 6 50       14 croak "The number of rows must be equal to the number of columns"
501             unless $m == $n;
502              
503 6         47 bless [ map [ (0) x ($_ - 1), $c, (0) x ($n - $_) ], 1 .. $m ], $class;
504             }
505              
506             =pod
507              
508             =item zeros()
509              
510             Create a zero matrix.
511              
512             # Create an $m-by-$m matrix where each element is 0.
513             $x = Math::Matrix -> zeros($m);
514              
515             # Create an $m-by-$n matrix where each element is 0.
516             $x = Math::Matrix -> zeros($m, $n);
517              
518             =cut
519              
520             sub zeros {
521 13 50   13 1 3887 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
522 13 50       26 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
523 13         76 my $self = shift;
524 13         39 $self -> constant(0, @_);
525             };
526              
527             =pod
528              
529             =item ones()
530              
531             Create a matrix of ones.
532              
533             # Create an $m-by-$m matrix where each element is 1.
534             $x = Math::Matrix -> ones($m);
535              
536             # Create an $m-by-$n matrix where each element is 1.
537             $x = Math::Matrix -> ones($m, $n);
538              
539             =cut
540              
541             sub ones {
542 5 50   5 1 3855 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
543 5 50       12 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
544 5         8 my $self = shift;
545 5         12 $self -> constant(1, @_);
546             };
547              
548             =pod
549              
550             =item inf()
551              
552             Create a matrix of positive infinities.
553              
554             # Create an $m-by-$m matrix where each element is Inf.
555             $x = Math::Matrix -> inf($m);
556              
557             # Create an $m-by-$n matrix where each element is Inf.
558             $x = Math::Matrix -> inf($m, $n);
559              
560             =cut
561              
562             sub inf {
563 5 50   5 1 16254 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
564 5 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
565 5         9 my $self = shift;
566              
567 5         26 require Math::Trig;
568 5         31 my $inf = Math::Trig::Inf();
569 5         21 $self -> constant($inf, @_);
570             };
571              
572             =pod
573              
574             =item nan()
575              
576             Create a matrix of NaNs (Not-a-Number).
577              
578             # Create an $m-by-$m matrix where each element is NaN.
579             $x = Math::Matrix -> nan($m);
580              
581             # Create an $m-by-$n matrix where each element is NaN.
582             $x = Math::Matrix -> nan($m, $n);
583              
584             =cut
585              
586             sub nan {
587 5 50   5 1 3847 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
588 5 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
589 5         9 my $self = shift;
590              
591 5         23 require Math::Trig;
592 5         11 my $inf = Math::Trig::Inf();
593 5         15 my $nan = $inf - $inf;
594 5         12 $self -> constant($nan, @_);
595             };
596              
597             =pod
598              
599             =item constant()
600              
601             Returns a constant matrix, i.e., a matrix whose elements all have the same
602             value.
603              
604             # Create an $m-by-$m matrix where each element is $c.
605             $x = Math::Matrix -> constant($c, $m);
606              
607             # Create an $m-by-$n matrix where each element is $c.
608             $x = Math::Matrix -> constant($c, $m, $n);
609              
610             =cut
611              
612             sub constant {
613 33 50   33 1 3939 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
614 33 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 4;
615 33         45 my $class = shift;
616              
617 33 50       56 croak +(caller(0))[3], " is a class method, not an instance method"
618             if ref $class;
619              
620 33         39 my $c = shift;
621 33 100       92 my ($m, $n) = @_ == 0 ? (1, 1)
    100          
622             : @_ == 1 ? (@_, @_)
623             : (@_);
624              
625 33         210 bless [ map [ ($c) x $n ], 1 .. $m ], $class;
626             }
627              
628             =pod
629              
630             =item rand()
631              
632             Returns a matrix of uniformly distributed random numbers in the range [0,1).
633              
634             $x = Math::Matrix -> rand($m); # $m-by-$m matrix
635             $x = Math::Matrix -> rand($m, $n); # $m-by-$n matrix
636              
637             To generate an C<$m>-by-C<$n> matrix of uniformly distributed random numbers in
638             the range [0,C<$a>), use
639              
640             $x = $a * Math::Matrix -> rand($m, $n);
641              
642             To generate an C<$m>-by-C<$n> matrix of uniformly distributed random numbers in
643             the range [C<$a>,C<$b>), use
644              
645             $x = $a + ($b - $a) * Math::Matrix -> rand($m, $n);
646              
647             See also C> and C>.
648              
649             =cut
650              
651             sub rand {
652 6 50   6 1 12337 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
653 6 50       16 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
654 6         11 my $class = shift;
655              
656 6 50       11 croak +(caller(0))[3], " is a class method, not an instance method"
657             if ref $class;
658              
659 6 100       23 my ($nrow, $ncol) = @_ == 0 ? (1, 1)
    100          
660             : @_ == 1 ? (@_, @_)
661             : (@_);
662              
663 6         13 my $x = bless [], $class;
664 6         20 for my $i (0 .. $nrow - 1) {
665 10         19 for my $j (0 .. $ncol - 1) {
666 22         114 $x -> [$i][$j] = CORE::rand;
667             }
668             }
669              
670 6         16 return $x;
671             }
672              
673             =pod
674              
675             =item randi()
676              
677             Returns a matrix of uniformly distributed random integers.
678              
679             $x = Math::Matrix -> randi($max); # 1-by-1 matrix
680             $x = Math::Matrix -> randi($max, $n); # $n-by-$n matrix
681             $x = Math::Matrix -> randi($max, $m, $n); # $m-by-$n matrix
682              
683             $x = Math::Matrix -> randi([$min, $max]); # 1-by-1 matrix
684             $x = Math::Matrix -> randi([$min, $max], $n); # $n-by-$n matrix
685             $x = Math::Matrix -> randi([$min, $max], $m, $n); # $m-by-$n matrix
686              
687             See also C> and C>.
688              
689             =cut
690              
691             sub randi {
692 6 50   6 1 17959 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
693 6 50       18 croak "Too many arguments for ", (caller(0))[3] if @_ > 4;
694 6         13 my $class = shift;
695              
696 6 50       12 croak +(caller(0))[3], " is a class method, not an instance method"
697             if ref $class;
698              
699 6         8 my ($min, $max);
700 6         8 my $lim = shift;
701 6 100       18 if (ref($lim) eq 'ARRAY') {
702 3         7 ($min, $max) = @$lim;
703             } else {
704 3         4 $min = 0;
705 3         6 $max = $lim;
706             }
707              
708 6 100       23 my ($nrow, $ncol) = @_ == 0 ? (1, 1)
    100          
709             : @_ == 1 ? (@_, @_)
710             : (@_);
711              
712 6         11 my $c = $max - $min + 1;
713 6         12 my $x = bless [], $class;
714 6         17 for my $i (0 .. $nrow - 1) {
715 14         23 for my $j (0 .. $ncol - 1) {
716 50         167 $x -> [$i][$j] = $min + CORE::int(CORE::rand($c));
717             }
718             }
719              
720 6         19 return $x;
721             }
722              
723             =pod
724              
725             =item randn()
726              
727             Returns a matrix of random numbers from the standard normal distribution.
728              
729             $x = Math::Matrix -> randn($m); # $m-by-$m matrix
730             $x = Math::Matrix -> randn($m, $n); # $m-by-$n matrix
731              
732             To generate an C<$m>-by-C<$n> matrix with mean C<$mu> and standard deviation
733             C<$sigma>, use
734              
735             $x = $mu + $sigma * Math::Matrix -> randn($m, $n);
736              
737             See also C> and C>.
738              
739             =cut
740              
741             sub randn {
742 6 50   6 1 6679 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
743 6 50       14 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
744 6         9 my $class = shift;
745              
746 6 50       13 croak +(caller(0))[3], " is a class method, not an instance method"
747             if ref $class;
748              
749 6 100       21 my ($nrow, $ncol) = @_ == 0 ? (1, 1)
    100          
750             : @_ == 1 ? (@_, @_)
751             : (@_);
752              
753 6         13 my $nelm = $nrow * $ncol;
754 6         8 my $twopi = 2 * atan2 0, -1;
755              
756             # The following might generate one value more than we need.
757              
758 6         12 my $x = [];
759 6         15 for (my $k = 0 ; $k < $nelm ; $k += 2) {
760 13         73 my $c1 = sqrt(-2 * log(CORE::rand));
761 13         23 my $c2 = $twopi * CORE::rand;
762 13         70 push @$x, $c1 * cos($c2), $c1 * sin($c2);
763             }
764 6 100       15 pop @$x if @$x > $nelm;
765              
766 6         15 $x = bless [ $x ], $class;
767 6         15 $x -> reshape($nrow, $ncol);
768             }
769              
770             =pod
771              
772             =item clone()
773              
774             Clones a matrix and returns the clone.
775              
776             $b = $a->clone;
777              
778             =cut
779              
780             sub clone {
781 127 50   127 1 2142 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
782 127         183 my $x = shift;
783 127         215 my $class = ref $x;
784              
785 127 50       262 croak +(caller(0))[3], " is an instance method, not a class method"
786             unless $class;
787              
788 127         696 my $y = [ map [ @$_ ], @$x ];
789 127         302 bless $y, $class;
790             }
791              
792             =pod
793              
794             =item diagonal()
795              
796             A constructor method that creates a diagonal matrix from a single list or array
797             of numbers.
798              
799             $p = Math::Matrix->diagonal(1, 4, 4, 8);
800             $q = Math::Matrix->diagonal([1, 4, 4, 8]);
801              
802             The matrix is zero filled except for the diagonal members, which take the
803             values of the vector.
804              
805             The method returns B in case of error.
806              
807             =cut
808              
809             #
810             # Either class or object call, create a square matrix with the same
811             # dimensions as the passed-in list or array.
812             #
813             sub diagonal {
814 2     2 1 2707 my $that = shift;
815 2   33     12 my $class = ref($that) || $that;
816 2         7 my @diag = @_;
817 2         4 my $self = [];
818              
819             # diagonal([2,3]) -> diagonal(2,3)
820 2 100       7 @diag = @{$diag[0]} if (ref $diag[0] eq "ARRAY");
  1         4  
821              
822 2         3 my $len = scalar @diag;
823 2 50       5 return undef if ($len == 0);
824              
825 2         8 for my $idx (0..$len-1) {
826 8         15 my @r = (0) x $len;
827 8         10 $r[$idx] = $diag[$idx];
828 8         12 push(@{$self}, [@r]);
  8         19  
829             }
830 2         7 bless $self, $class;
831             }
832              
833             =pod
834              
835             =item tridiagonal()
836              
837             A constructor method that creates a matrix from vectors of numbers.
838              
839             $p = Math::Matrix->tridiagonal([1, 4, 4, 8]);
840             $q = Math::Matrix->tridiagonal([1, 4, 4, 8], [9, 12, 15]);
841             $r = Math::Matrix->tridiagonal([1, 4, 4, 8], [9, 12, 15], [4, 3, 2]);
842              
843             In the first case, the main diagonal takes the values of the vector, while both
844             of the upper and lower diagonals's values are all set to one.
845              
846             In the second case, the main diagonal takes the values of the first vector,
847             while the upper and lower diagonals are each set to the values of the second
848             vector.
849              
850             In the third case, the main diagonal takes the values of the first vector,
851             while the upper diagonal is set to the values of the second vector, and the
852             lower diagonal is set to the values of the third vector.
853              
854             The method returns B in case of error.
855              
856             =cut
857              
858             #
859             # Either class or object call, create a square matrix with the same
860             # dimensions as the passed-in list or array.
861             #
862             sub tridiagonal {
863 4     4 1 7513 my $that = shift;
864 4   33     21 my $class = ref($that) || $that;
865 4         11 my(@up_d, @main_d, @low_d);
866 4         6 my $self = [];
867              
868             #
869             # Handle the different ways the tridiagonal vectors could
870             # be passed in.
871             #
872 4 50       22 if (ref $_[0] eq "ARRAY") {
873 4         12 @main_d = @{$_[0]};
  4         10  
874              
875 4 100       11 if (ref $_[1] eq "ARRAY") {
876 3         4 @up_d = @{$_[1]};
  3         7  
877              
878 3 100       9 if (ref $_[2] eq "ARRAY") {
879 2         3 @low_d = @{$_[2]};
  2         3  
880             }
881             }
882             } else {
883 0         0 @main_d = @_;
884             }
885              
886 4         7 my $len = scalar @main_d;
887 4 50       11 return undef if ($len == 0);
888              
889             #
890             # Default the upper and lower diagonals if no vector
891             # was passed in for them.
892             #
893 4 100       10 @up_d = (1) x ($len -1) if (scalar @up_d == 0);
894 4 100       12 @low_d = @up_d if (scalar @low_d == 0);
895              
896             #
897             # First row...
898             #
899 4         9 my @arow = (0) x $len;
900 4         9 @arow[0..1] = ($main_d[0], $up_d[0]);
901 4         6 push (@{$self}, [@arow]);
  4         10  
902              
903             #
904             # Bulk of the matrix...
905             #
906 4         14 for my $idx (1 .. $#main_d - 1) {
907 8         14 my @r = (0) x $len;
908 8         23 @r[$idx-1 .. $idx+1] = ($low_d[$idx-1], $main_d[$idx], $up_d[$idx]);
909 8         11 push (@{$self}, [@r]);
  8         18  
910             }
911              
912             #
913             # Last row.
914             #
915 4         9 my @zrow = (0) x $len;
916 4         10 @zrow[$len-2..$len-1] = ($low_d[$#main_d -1], $main_d[$#main_d]);
917 4         6 push (@{$self}, [@zrow]);
  4         8  
918              
919 4         16 bless $self, $class;
920             }
921              
922             =pod
923              
924             =item blkdiag()
925              
926             Create block diagonal matrix. Returns a block diagonal matrix given a list of
927             matrices.
928              
929             $z = Math::Matrix -> blkdiag($x, $y, ...);
930              
931             =cut
932              
933             sub blkdiag {
934 4 50   4 1 3138 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
935             #croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
936 4         8 my $class = shift;
937              
938 4         6 my $y = [];
939 4         5 my $nrowy = 0;
940 4         5 my $ncoly = 0;
941              
942 4         12 for my $i (0 .. $#_) {
943 9         12 my $x = $_[$i];
944              
945 9 100 66     41 $x = $class -> new($x)
946             unless defined(blessed($x)) && $x -> isa($class);
947              
948 9         18 my ($nrowx, $ncolx) = $x -> size();
949              
950             # Upper right submatrix.
951              
952 9         17 for my $i (0 .. $nrowy - 1) {
953 9         14 for my $j (0 .. $ncolx - 1) {
954 11         16 $y -> [$i][$ncoly + $j] = 0;
955             }
956             }
957              
958             # Lower left submatrix.
959              
960 9         14 for my $i (0 .. $nrowx - 1) {
961 11         15 for my $j (0 .. $ncoly - 1) {
962 17         25 $y -> [$nrowy + $i][$j] = 0;
963             }
964             }
965              
966             # Lower right submatrix.
967              
968 9         13 for my $i (0 .. $nrowx - 1) {
969 11         14 for my $j (0 .. $ncolx - 1) {
970 13         22 $y -> [$nrowy + $i][$ncoly + $j] = $x -> [$i][$j];
971             }
972             }
973              
974 9         10 $nrowy += $nrowx;
975 9         18 $ncoly += $ncolx;
976             }
977              
978 4         9 bless $y, $class;
979             }
980              
981             =pod
982              
983             =back
984              
985             =head2 Identify matrices
986              
987             =over 4
988              
989             =item is_empty()
990              
991             Returns 1 is the invocand is empty, i.e., it has no elements.
992              
993             $bool = $x -> is_empty();
994              
995             =cut
996              
997             sub is_empty {
998 56 50   56 1 255 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
999 56 50       119 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1000 56         89 my $x = shift;
1001 56         138 return $x -> nelm() == 0;
1002             }
1003              
1004             =pod
1005              
1006             =item is_scalar()
1007              
1008             Returns 1 is the invocand is a scalar, i.e., it has one element.
1009              
1010             $bool = $x -> is_scalar();
1011              
1012             =cut
1013              
1014             sub is_scalar {
1015 111 50   111 1 254 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1016 111 50       223 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1017 111         158 my $x = shift;
1018 111 100       223 return $x -> nelm() == 1 ? 1 : 0;
1019             }
1020              
1021             =pod
1022              
1023             =item is_vector()
1024              
1025             Returns 1 is the invocand is a vector, i.e., a row vector or a column vector.
1026              
1027             $bool = $x -> is_vector();
1028              
1029             =cut
1030              
1031             sub is_vector {
1032 52 50   52 1 155 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1033 52 50       113 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1034 52         82 my $x = shift;
1035 52 100 100     128 return $x -> is_col() || $x -> is_row() ? 1 : 0;
1036             }
1037              
1038             =pod
1039              
1040             =item is_row()
1041              
1042             Returns 1 if the invocand has exactly one row, and 0 otherwise.
1043              
1044             $bool = $x -> is_row();
1045              
1046             =cut
1047              
1048             sub is_row {
1049 73 50   73 1 185 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1050 73 50       150 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1051 73         126 my $x = shift;
1052 73 100       160 return $x -> nrow() == 1 ? 1 : 0;
1053             }
1054              
1055             =pod
1056              
1057             =item is_col()
1058              
1059             Returns 1 if the invocand has exactly one column, and 0 otherwise.
1060              
1061             $bool = $x -> is_col();
1062              
1063             =cut
1064              
1065             sub is_col {
1066 81 50   81 1 199 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1067 81 50       173 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1068 81         122 my $x = shift;
1069 81 100       190 return $x -> ncol() == 1 ? 1 : 0;
1070             }
1071              
1072             =pod
1073              
1074             =item is_square()
1075              
1076             Returns 1 is the invocand is square, and 0 otherwise.
1077              
1078             $bool = $x -> is_square();
1079              
1080             =cut
1081              
1082             sub is_square {
1083 78 50   78 1 182 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1084 78 50       154 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1085 78         110 my $x = shift;
1086 78         211 my ($nrow, $ncol) = $x -> size();
1087 78 100       316 return $nrow == $ncol ? 1 : 0;
1088             }
1089              
1090             =pod
1091              
1092             =item is_symmetric()
1093              
1094             Returns 1 is the invocand is symmetric, and 0 otherwise.
1095              
1096             $bool = $x -> is_symmetric();
1097              
1098             An symmetric matrix satisfies x(i,j) = x(j,i) for all i and j, for example
1099              
1100             [ 1 2 -3 ]
1101             [ 2 -4 5 ]
1102             [ -3 5 6 ]
1103              
1104             =cut
1105              
1106             sub is_symmetric {
1107 50 50   50 1 160 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1108 50 50       104 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1109 50         77 my $x = shift;
1110              
1111 50         142 my ($nrow, $ncol) = $x -> size();
1112 50 100       132 return 0 unless $nrow == $ncol;
1113              
1114 44         122 for my $i (1 .. $nrow - 1) {
1115 100         170 for my $j (0 .. $i - 1) {
1116 232 100       883 return 0 unless $x -> [$i][$j] == $x -> [$j][$i];
1117             }
1118             }
1119              
1120 25         135 return 1;
1121             }
1122              
1123             =pod
1124              
1125             =item is_antisymmetric()
1126              
1127             Returns 1 is the invocand is antisymmetric a.k.a. skew-symmetric, and 0
1128             otherwise.
1129              
1130             $bool = $x -> is_antisymmetric();
1131              
1132             An antisymmetric matrix satisfies x(i,j) = -x(j,i) for all i and j, for
1133             example
1134              
1135             [ 0 2 -3 ]
1136             [ -2 0 4 ]
1137             [ 3 -4 0 ]
1138              
1139             =cut
1140              
1141             sub is_antisymmetric {
1142 25 50   25 1 85 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1143 25 50       58 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1144 25         50 my $x = shift;
1145              
1146 25         64 my ($nrow, $ncol) = $x -> size();
1147 25 100       73 return 0 unless $nrow == $ncol;
1148              
1149             # Check the diagonal.
1150              
1151 22         58 for my $i (0 .. $nrow - 1) {
1152 38 100       156 return 0 unless $x -> [$i][$i] == 0;
1153             }
1154              
1155             # Check the off-diagonal.
1156              
1157 5         18 for my $i (1 .. $nrow - 1) {
1158 5         13 for my $j (0 .. $i - 1) {
1159 6 100       30 return 0 unless $x -> [$i][$j] == -$x -> [$j][$i];
1160             }
1161             }
1162              
1163 2         11 return 1;
1164             }
1165              
1166             =pod
1167              
1168             =item is_persymmetric()
1169              
1170             Returns 1 is the invocand is persymmetric, and 0 otherwise.
1171              
1172             $bool = $x -> is_persymmetric();
1173              
1174             A persymmetric matrix is a square matrix which is symmetric with respect to the
1175             anti-diagonal, e.g.:
1176              
1177             [ f h j k ]
1178             [ c g i j ]
1179             [ b d g h ]
1180             [ a b c f ]
1181              
1182             =cut
1183              
1184             sub is_persymmetric {
1185 23 50   23 1 74 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1186 23 50       52 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1187 23         40 my $x = shift;
1188              
1189 23         62 $x -> fliplr() -> is_symmetric();
1190             }
1191              
1192             =pod
1193              
1194             =item is_hankel()
1195              
1196             Returns 1 is the invocand is a Hankel matric a.k.a. a catalecticant matrix, and
1197             0 otherwise.
1198              
1199             $bool = $x -> is_hankel();
1200              
1201             A Hankel matrix is a square matrix in which each ascending skew-diagonal from
1202             left to right is constant, e.g.:
1203              
1204             [ e f g h i ]
1205             [ d e f g h ]
1206             [ c d e f g ]
1207             [ b c d e f ]
1208             [ a b c d e ]
1209              
1210             =cut
1211              
1212             sub is_hankel {
1213 23 50   23 1 79 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1214 23 50       48 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1215 23         45 my $x = shift;
1216              
1217 23         57 my ($nrow, $ncol) = $x -> size();
1218 23 100       80 return 0 unless $nrow == $ncol;
1219              
1220             # Check the lower triangular part.
1221              
1222 20         54 for my $i (0 .. $nrow - 2) {
1223 34         49 my $first = $x -> [$i][0];
1224 34         59 for my $k (1 .. $nrow - $i - 1) {
1225 71 100       236 return 0 unless $x -> [$i + $k][$k] == $first;
1226             }
1227             }
1228              
1229             # Check the strictly upper triangular part.
1230              
1231 8         20 for my $j (1 .. $ncol - 2) {
1232 15         27 my $first = $x -> [0][$j];
1233 15         30 for my $k (1 .. $nrow - $j - 1) {
1234 33 100       75 return 0 unless $x -> [$k][$j + $k] == $first;
1235             }
1236             }
1237              
1238 7         33 return 1;
1239             }
1240              
1241             =pod
1242              
1243             =item is_zero()
1244              
1245             Returns 1 is the invocand is a zero matrix, and 0 otherwise. A zero matrix
1246             contains no element whose value is different from zero.
1247              
1248             $bool = $x -> is_zero();
1249              
1250             =cut
1251              
1252             sub is_zero {
1253 25 50   25 1 86 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1254 25 50       52 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1255 25         42 my $x = shift;
1256 25         61 return $x -> is_constant(0);
1257             }
1258              
1259             =pod
1260              
1261             =item is_one()
1262              
1263             Returns 1 is the invocand is a matrix of ones, and 0 otherwise. A matrix of
1264             ones contains no element whose value is different from one.
1265              
1266             $bool = $x -> is_one();
1267              
1268             =cut
1269              
1270             sub is_one {
1271 25 50   25 1 79 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1272 25 50       55 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1273 25         40 my $x = shift;
1274 25         57 return $x -> is_constant(1);
1275             }
1276              
1277             =pod
1278              
1279             =item is_constant()
1280              
1281             Returns 1 is the invocand is a constant matrix, and 0 otherwise. A constant
1282             matrix is a matrix where no two elements have different values.
1283              
1284             $bool = $x -> is_constant();
1285              
1286             =cut
1287              
1288             sub is_constant {
1289 75 50   75 1 189 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1290 75 50       158 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
1291 75         113 my $x = shift;
1292              
1293 75         158 my ($nrow, $ncol) = $x -> size();
1294              
1295             # An empty matrix contains no elements that are different from $c.
1296              
1297 75 100       221 return 1 if $nrow * $ncol == 0;
1298              
1299 72 100       139 my $c = @_ ? shift() : $x -> [0][0];
1300 72         217 for my $i (0 .. $nrow - 1) {
1301 78         142 for my $j (0 .. $ncol - 1) {
1302 147 100       785 return 0 if $x -> [$i][$j] != $c;
1303             }
1304             }
1305              
1306 3         16 return 1;
1307             }
1308              
1309             =pod
1310              
1311             =item is_identity()
1312              
1313             Returns 1 is the invocand is an identity matrix, and 0 otherwise. An
1314             identity matrix contains ones on the main diagonal and zeros elsewhere.
1315              
1316             $bool = $x -> is_identity();
1317              
1318             =cut
1319              
1320             sub is_identity {
1321 25 50   25 1 84 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1322 25 50       53 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1323 25         47 my $x = shift;
1324              
1325 25         57 my ($nrow, $ncol) = $x -> size();
1326 25 100       73 return 0 unless $nrow == $ncol;
1327              
1328 22         65 for my $i (0 .. $nrow - 1) {
1329 24         45 for my $j (0 .. $ncol - 1) {
1330 36 100       164 return 0 if $x -> [$i][$j] != ($i == $j ? 1 : 0);
    100          
1331             }
1332             }
1333              
1334 3         17 return 1;
1335             }
1336              
1337             =pod
1338              
1339             =item is_exchg()
1340              
1341             Returns 1 is the invocand is an exchange matrix, and 0 otherwise.
1342              
1343             $bool = $x -> is_exchg();
1344              
1345             An exchange matrix contains ones on the main anti-diagonal and zeros elsewhere,
1346             for example
1347              
1348             [ 0 0 1 ]
1349             [ 0 1 0 ]
1350             [ 1 0 0 ]
1351              
1352             =cut
1353              
1354             sub is_exchg {
1355 25 50   25 1 81 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1356 25 50       60 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1357 25         42 my $x = shift;
1358              
1359 25         59 my ($nrow, $ncol) = $x -> size();
1360 25 100       73 return 0 unless $nrow == $ncol;
1361              
1362 22         41 my $imax = $nrow - 1;
1363 22         58 for my $i (0 .. $nrow - 1) {
1364 23         37 for my $j (0 .. $ncol - 1) {
1365 47 100       222 return 0 if $x -> [$i][$j] != ($i + $j == $imax ? 1 : 0);
    100          
1366             }
1367             }
1368              
1369 3         18 return 1;
1370             }
1371              
1372             =pod
1373              
1374             =item is_bool()
1375              
1376             Returns 1 is the invocand is a boolean matrix, and 0 otherwise.
1377              
1378             $bool = $x -> is_bool();
1379              
1380             A boolean matrix is a matrix is a matrix whose entries are either 0 or 1, for
1381             example
1382              
1383             [ 0 1 1 ]
1384             [ 1 0 0 ]
1385             [ 0 1 0 ]
1386              
1387             =cut
1388              
1389             sub is_bool {
1390 25 50   25 1 83 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1391 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1392 25         37 my $x = shift;
1393              
1394 25         65 my ($nrow, $ncol) = $x -> size();
1395              
1396 25         73 for my $i (0 .. $nrow - 1) {
1397 34         65 for my $j (0 .. $ncol - 1) {
1398 80         117 my $val = $x -> [$i][$j];
1399 80 100 100     315 return 0 if $val != 0 && $val != 1;
1400             }
1401             }
1402              
1403 5         30 return 1;
1404             }
1405              
1406             =pod
1407              
1408             =item is_perm()
1409              
1410             Returns 1 is the invocand is an permutation matrix, and 0 otherwise.
1411              
1412             $bool = $x -> is_perm();
1413              
1414             A permutation matrix is a square matrix with exactly one 1 in each row and
1415             column, and all other elements 0, for example
1416              
1417             [ 0 1 0 ]
1418             [ 1 0 0 ]
1419             [ 0 0 1 ]
1420              
1421             =cut
1422              
1423             sub is_perm {
1424 25 50   25 1 82 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1425 25 50       53 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1426 25         44 my $x = shift;
1427              
1428 25         63 my ($nrow, $ncol) = $x -> size();
1429 25 100       71 return 0 unless $nrow == $ncol;
1430              
1431 22         57 my $rowsum = [ (0) x $nrow ];
1432 22         43 my $colsum = [ (0) x $ncol ];
1433              
1434 22         55 for my $i (0 .. $nrow - 1) {
1435 30         50 for my $j (0 .. $ncol - 1) {
1436 74         94 my $val = $x -> [$i][$j];
1437 74 100 100     285 return 0 if $val != 0 && $val != 1;
1438 57 100       110 if ($val == 1) {
1439 16 50       32 return 0 if ++$rowsum -> [$i] > 1;
1440 16 50       34 return 0 if ++$colsum -> [$j] > 1;
1441             }
1442             }
1443             }
1444              
1445 5         11 for my $i (0 .. $nrow - 1) {
1446 10 50       18 return 0 if $rowsum -> [$i] != 1;
1447 10 50       21 return 0 if $colsum -> [$i] != 1;
1448             }
1449              
1450 5         26 return 1;
1451             }
1452              
1453             =pod
1454              
1455             =item is_int()
1456              
1457             Returns 1 is the invocand is an integer matrix, i.e., a matrix of integers, and
1458             0 otherwise.
1459              
1460             $bool = $x -> is_int();
1461              
1462             =cut
1463              
1464             sub is_int {
1465 25 50   25 1 91 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1466 25 50       106 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1467 25         45 my $x = shift;
1468              
1469 25         58 my ($nrow, $ncol) = $x -> size();
1470              
1471 25         72 for my $i (0 .. $nrow - 1) {
1472 83         132 for my $j (0 .. $ncol - 1) {
1473 347 50       718 return 0 unless $x -> [$i][$j] == int $x -> [$i][$j];
1474             }
1475             }
1476              
1477 25         116 return 1;
1478             }
1479              
1480             =pod
1481              
1482             =item is_diag()
1483              
1484             Returns 1 is the invocand is diagonal, and 0 otherwise.
1485              
1486             $bool = $x -> is_diag();
1487              
1488             A diagonal matrix is a square matrix where all non-zero elements, if any, are on
1489             the main diagonal. It has the following pattern, where only the elements marked
1490             as C can be non-zero,
1491              
1492             [ x 0 0 0 0 ]
1493             [ 0 x 0 0 0 ]
1494             [ 0 0 x 0 0 ]
1495             [ 0 0 0 x 0 ]
1496             [ 0 0 0 0 x ]
1497              
1498             =cut
1499              
1500             sub is_diag {
1501 25 50   25 1 84 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1502 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1503 25         44 my $x = shift;
1504 25         67 $x -> is_band(0);
1505             }
1506              
1507             =pod
1508              
1509             =item is_adiag()
1510              
1511             Returns 1 is the invocand is anti-diagonal, and 0 otherwise.
1512              
1513             $bool = $x -> is_adiag();
1514              
1515             A diagonal matrix is a square matrix where all non-zero elements, if any, are on
1516             the main antidiagonal. It has the following pattern, where only the elements
1517             marked as C can be non-zero,
1518              
1519             [ 0 0 0 0 x ]
1520             [ 0 0 0 x 0 ]
1521             [ 0 0 x 0 0 ]
1522             [ 0 x 0 0 0 ]
1523             [ x 0 0 0 0 ]
1524              
1525             =cut
1526              
1527             sub is_adiag {
1528 25 50   25 1 83 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1529 25 50       56 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1530 25         45 my $x = shift;
1531 25         62 $x -> is_aband(0);
1532             }
1533              
1534             =pod
1535              
1536             =item is_tridiag()
1537              
1538             Returns 1 is the invocand is tridiagonal, and 0 otherwise.
1539              
1540             $bool = $x -> is_tridiag();
1541              
1542             A tridiagonal matrix is a square matrix with nonzero elements only on the
1543             diagonal and slots horizontally or vertically adjacent the diagonal (i.e., along
1544             the subdiagonal and superdiagonal). It has the following pattern, where only the
1545             elements marked as C can be non-zero,
1546              
1547             [ x x 0 0 0 ]
1548             [ x x x 0 0 ]
1549             [ 0 x x x 0 ]
1550             [ 0 0 x x x ]
1551             [ 0 0 0 x x ]
1552              
1553             =cut
1554              
1555             sub is_tridiag {
1556 25 50   25 1 83 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1557 25 50       56 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1558 25         42 my $x = shift;
1559 25         61 $x -> is_band(1);
1560             }
1561              
1562             =pod
1563              
1564             =item is_atridiag()
1565              
1566             Returns 1 is the invocand is anti-tridiagonal, and 0 otherwise.
1567              
1568             $bool = $x -> is_tridiag();
1569              
1570             A anti-tridiagonal matrix is a square matrix with nonzero elements only on the
1571             anti-diagonal and slots horizontally or vertically adjacent the diagonal (i.e.,
1572             along the anti-subdiagonal and anti-superdiagonal). It has the following
1573             pattern, where only the elements marked as C can be non-zero,
1574              
1575             [ 0 0 0 x x ]
1576             [ 0 0 x x x ]
1577             [ 0 x x x 0 ]
1578             [ x x x 0 0 ]
1579             [ x x 0 0 0 ]
1580              
1581             =cut
1582              
1583             sub is_atridiag {
1584 25 50   25 1 77 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1585 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1586 25         38 my $x = shift;
1587 25         59 $x -> is_aband(1);
1588             }
1589              
1590             =pod
1591              
1592             =item is_pentadiag()
1593              
1594             Returns 1 is the invocand is pentadiagonal, and 0 otherwise.
1595              
1596             $bool = $x -> is_pentadiag();
1597              
1598             A pentadiagonal matrix is a square matrix with nonzero elements only on the
1599             diagonal and the two diagonals above and below the main diagonal. It has the
1600             following pattern, where only the elements marked as C can be non-zero,
1601              
1602             [ x x x 0 0 0 ]
1603             [ x x x x 0 0 ]
1604             [ x x x x x 0 ]
1605             [ 0 x x x x x ]
1606             [ 0 0 x x x x ]
1607             [ 0 0 0 x x x ]
1608              
1609             =cut
1610              
1611             sub is_pentadiag {
1612 25 50   25 1 82 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1613 25 50       69 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1614 25         41 my $x = shift;
1615 25         57 $x -> is_band(2);
1616             }
1617              
1618             =pod
1619              
1620             =item is_apentadiag()
1621              
1622             Returns 1 is the invocand is anti-pentadiagonal, and 0 otherwise.
1623              
1624             $bool = $x -> is_pentadiag();
1625              
1626             A anti-pentadiagonal matrix is a square matrix with nonzero elements only on the
1627             anti-diagonal and two anti-diagonals above and below the main anti-diagonal. It
1628             has the following pattern, where only the elements marked as C can be
1629             non-zero,
1630              
1631             [ 0 0 0 x x x ]
1632             [ 0 0 x x x x ]
1633             [ 0 x x x x x ]
1634             [ x x x x x 0 ]
1635             [ x x x x 0 0 ]
1636             [ x x x 0 0 0 ]
1637              
1638             =cut
1639              
1640             sub is_apentadiag {
1641 25 50   25 1 85 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1642 25 50       86 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1643 25         47 my $x = shift;
1644 25         61 $x -> is_aband(2);
1645             }
1646              
1647             =pod
1648              
1649             =item is_heptadiag()
1650              
1651             Returns 1 is the invocand is heptadiagonal, and 0 otherwise.
1652              
1653             $bool = $x -> is_heptadiag();
1654              
1655             A heptadiagonal matrix is a square matrix with nonzero elements only on the
1656             diagonal and the two diagonals above and below the main diagonal. It has the
1657             following pattern, where only the elements marked as C can be non-zero,
1658              
1659             [ x x x x 0 0 ]
1660             [ x x x x x 0 ]
1661             [ x x x x x x ]
1662             [ x x x x x x ]
1663             [ 0 x x x x x ]
1664             [ 0 0 x x x x ]
1665              
1666             =cut
1667              
1668             sub is_heptadiag {
1669 25 50   25 1 84 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1670 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1671 25         41 my $x = shift;
1672 25         54 $x -> is_band(3);
1673             }
1674              
1675             =pod
1676              
1677             =item is_aheptadiag()
1678              
1679             Returns 1 is the invocand is anti-heptadiagonal, and 0 otherwise.
1680              
1681             $bool = $x -> is_heptadiag();
1682              
1683             A anti-heptadiagonal matrix is a square matrix with nonzero elements only on the
1684             anti-diagonal and two anti-diagonals above and below the main anti-diagonal. It
1685             has the following pattern, where only the elements marked as C can be
1686             non-zero,
1687              
1688             [ 0 0 x x x x ]
1689             [ 0 x x x x x ]
1690             [ x x x x x x ]
1691             [ x x x x x x ]
1692             [ x x x x x 0 ]
1693             [ x x x x 0 0 ]
1694              
1695             =cut
1696              
1697             sub is_aheptadiag {
1698 25 50   25 1 85 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1699 25 50       55 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1700 25         37 my $x = shift;
1701 25         58 $x -> is_aband(3);
1702             }
1703              
1704             =pod
1705              
1706             =item is_band()
1707              
1708             Returns 1 is the invocand is a band matrix with a specified bandwidth, and 0
1709             otherwise.
1710              
1711             $bool = $x -> is_band($k);
1712              
1713             A band matrix is a square matrix with nonzero elements only on the diagonal and
1714             on the C<$k> diagonals above and below the main diagonal. The bandwidth C<$k>
1715             must be non-negative.
1716              
1717             $bool = $x -> is_band(0); # is $x diagonal?
1718             $bool = $x -> is_band(1); # is $x tridiagonal?
1719             $bool = $x -> is_band(2); # is $x pentadiagonal?
1720             $bool = $x -> is_band(3); # is $x heptadiagonal?
1721              
1722             See also C> and C>.
1723              
1724             =cut
1725              
1726             sub is_band {
1727 225 50   225 1 553 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
1728 225 50       435 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
1729 225         310 my $x = shift;
1730 225         434 my $class = ref $x;
1731              
1732 225         509 my ($nrow, $ncol) = $x -> size();
1733 225 100       558 return 0 unless $nrow == $ncol; # must be square
1734              
1735 198         368 my $k = shift; # bandwidth
1736 198 50       393 croak "Bandwidth can not be undefined" unless defined $k;
1737 198 50       361 if (ref $k) {
1738 0 0 0     0 $k = $class -> new($k)
1739             unless defined(blessed($k)) && $k -> isa($class);
1740 0 0       0 croak "Bandwidth must be a scalar" unless $k -> is_scalar();
1741 0         0 $k = $k -> [0][0];
1742             }
1743              
1744 198 100       692 return 0 if $nrow <= $k; # if the band doesn't fit inside
1745 136 100       441 return 1 if $nrow == $k + 1; # if the whole band fits exactly
1746              
1747 106         318 for my $i (0 .. $nrow - $k - 2) {
1748 130         235 for my $j ($k + 1 + $i .. $ncol - 1) {
1749 188 100 100     995 return 0 if ($x -> [$i][$j] != 0 ||
1750             $x -> [$j][$i] != 0);
1751             }
1752             }
1753              
1754 23         127 return 1;
1755             }
1756              
1757             =pod
1758              
1759             =item is_aband()
1760              
1761             Returns 1 is the invocand is "anti-banded" with a specified bandwidth, and 0
1762             otherwise.
1763              
1764             $bool = $x -> is_aband($k);
1765              
1766             Some examples
1767              
1768             $bool = $x -> is_aband(0); # is $x anti-diagonal?
1769             $bool = $x -> is_aband(1); # is $x anti-tridiagonal?
1770             $bool = $x -> is_aband(2); # is $x anti-pentadiagonal?
1771             $bool = $x -> is_aband(3); # is $x anti-heptadiagonal?
1772              
1773             A band matrix is a square matrix with nonzero elements only on the diagonal and
1774             on the C<$k> diagonals above and below the main diagonal. The bandwidth C<$k>
1775             must be non-negative.
1776              
1777             A "anti-banded" matrix is a square matrix with nonzero elements only on the
1778             anti-diagonal and C<$k> anti-diagonals above and below the main anti-diagonal.
1779              
1780             See also C> and C>.
1781              
1782             =cut
1783              
1784             sub is_aband {
1785 225 50   225 1 626 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
1786 225 50       415 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
1787 225         326 my $x = shift;
1788 225         389 my $class = ref $x;
1789              
1790 225         517 my ($nrow, $ncol) = $x -> size();
1791 225 100       614 return 0 unless $nrow == $ncol; # must be square
1792              
1793 198         308 my $k = shift; # bandwidth
1794 198 50       430 croak "Bandwidth can not be undefined" unless defined $k;
1795 198 50       364 if (ref $k) {
1796 0 0 0     0 $k = $class -> new($k)
1797             unless defined(blessed($k)) && $k -> isa($class);
1798 0 0       0 croak "Bandwidth must be a scalar" unless $k -> is_scalar();
1799 0         0 $k = $k -> [0][0];
1800             }
1801              
1802 198 100       716 return 0 if $nrow <= $k; # if the band doesn't fit inside
1803 136 100       428 return 1 if $nrow == $k + 1; # if the whole band fits exactly
1804              
1805             # Check upper part.
1806              
1807 106         278 for my $i (0 .. $nrow - $k - 2) {
1808 134         219 for my $j (0 .. $nrow - $k - 2 - $i) {
1809 210 100       751 return 0 if $x -> [$i][$j] != 0;
1810             }
1811             }
1812              
1813             # Check lower part.
1814              
1815 35         82 for my $i ($k + 1 .. $nrow - 1) {
1816 57         97 for my $j ($nrow - $i + $k .. $nrow - 1) {
1817 89 100       254 return 0 if $x -> [$i][$j] != 0;
1818             }
1819             }
1820              
1821 27         137 return 1;
1822             }
1823              
1824             =pod
1825              
1826             =item is_triu()
1827              
1828             Returns 1 is the invocand is upper triangular, and 0 otherwise.
1829              
1830             $bool = $x -> is_triu();
1831              
1832             An upper triangular matrix is a square matrix where all non-zero elements are on
1833             or above the main diagonal. It has the following pattern, where only the
1834             elements marked as C can be non-zero. It has the following pattern, where
1835             only the elements marked as C can be non-zero,
1836              
1837             [ x x x x ]
1838             [ 0 x x x ]
1839             [ 0 0 x x ]
1840             [ 0 0 0 x ]
1841              
1842             =cut
1843              
1844             sub is_triu {
1845 25 50   25 1 137 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1846 25 50       55 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1847 25         43 my $x = shift;
1848              
1849 25         64 my $nrow = $x -> nrow();
1850 25         56 my $ncol = $x -> ncol();
1851              
1852 25 100       72 return 0 unless $nrow == $ncol;
1853              
1854 22         63 for my $i (1 .. $nrow - 1) {
1855 30         55 for my $j (0 .. $i - 1) {
1856 37 100       157 return 0 unless $x -> [$i][$j] == 0;
1857             }
1858             }
1859              
1860 6         29 return 1;
1861             }
1862              
1863             =pod
1864              
1865             =item is_striu()
1866              
1867             Returns 1 is the invocand is strictly upper triangular, and 0 otherwise.
1868              
1869             $bool = $x -> is_striu();
1870              
1871             A strictly upper triangular matrix is a square matrix where all non-zero
1872             elements are strictly above the main diagonal. It has the following pattern,
1873             where only the elements marked as C can be non-zero,
1874              
1875             [ 0 x x x ]
1876             [ 0 0 x x ]
1877             [ 0 0 0 x ]
1878             [ 0 0 0 0 ]
1879              
1880             =cut
1881              
1882             sub is_striu {
1883 25 50   25 1 84 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1884 25 50       58 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1885 25         49 my $x = shift;
1886              
1887 25         56 my $nrow = $x -> nrow();
1888 25         52 my $ncol = $x -> ncol();
1889              
1890 25 100       67 return 0 unless $nrow == $ncol;
1891              
1892 22         61 for my $i (0 .. $nrow - 1) {
1893 36         66 for my $j (0 .. $i) {
1894 50 100       175 return 0 unless $x -> [$i][$j] == 0;
1895             }
1896             }
1897              
1898 2         9 return 1;
1899             }
1900              
1901             =pod
1902              
1903             =item is_tril()
1904              
1905             Returns 1 is the invocand is lower triangular, and 0 otherwise.
1906              
1907             $bool = $x -> is_tril();
1908              
1909             A lower triangular matrix is a square matrix where all non-zero elements are on
1910             or below the main diagonal. It has the following pattern, where only the
1911             elements marked as C can be non-zero,
1912              
1913             [ x 0 0 0 ]
1914             [ x x 0 0 ]
1915             [ x x x 0 ]
1916             [ x x x x ]
1917              
1918             =cut
1919              
1920             sub is_tril {
1921 25 50   25 1 92 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1922 25 50       58 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1923 25         38 my $x = shift;
1924              
1925 25         55 my $nrow = $x -> nrow();
1926 25         51 my $ncol = $x -> ncol();
1927              
1928 25 100       64 return 0 unless $nrow == $ncol;
1929              
1930 22         61 for my $i (0 .. $nrow - 1) {
1931 28         66 for my $j ($i + 1 .. $ncol - 1) {
1932 35 100       151 return 0 unless $x -> [$i][$j] == 0;
1933             }
1934             }
1935              
1936 6         38 return 1;
1937             }
1938              
1939             =pod
1940              
1941             =item is_stril()
1942              
1943             Returns 1 is the invocand is strictly lower triangular, and 0 otherwise.
1944              
1945             $bool = $x -> is_stril();
1946              
1947             A strictly lower triangular matrix is a square matrix where all non-zero
1948             elements are strictly below the main diagonal. It has the following pattern,
1949             where only the elements marked as C can be non-zero,
1950              
1951             [ 0 0 0 0 ]
1952             [ x 0 0 0 ]
1953             [ x x 0 0 ]
1954             [ x x x 0 ]
1955              
1956             =cut
1957              
1958             sub is_stril {
1959 25 50   25 1 80 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1960 25 50       55 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1961 25         42 my $x = shift;
1962              
1963 25         57 my $nrow = $x -> nrow();
1964 25         54 my $ncol = $x -> ncol();
1965              
1966 25 100       65 return 0 unless $nrow == $ncol;
1967              
1968 22         64 for my $i (0 .. $nrow - 1) {
1969 24         47 for my $j ($i .. $ncol - 1) {
1970 46 100       191 return 0 unless $x -> [$i][$j] == 0;
1971             }
1972             }
1973              
1974 2         23 return 1;
1975             }
1976              
1977             =pod
1978              
1979             =item is_atriu()
1980              
1981             Returns 1 is the invocand is upper anti-triangular, and 0 otherwise.
1982              
1983             $bool = $x -> is_atriu();
1984              
1985             An upper anti-triangular matrix is a square matrix where all non-zero elements
1986             are on or above the main anti-diagonal. It has the following pattern, where only
1987             the elements marked as C can be non-zero,
1988              
1989             [ x x x x ]
1990             [ x x x 0 ]
1991             [ x x 0 0 ]
1992             [ x 0 0 0 ]
1993              
1994             =cut
1995              
1996             sub is_atriu {
1997 25 50   25 1 78 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
1998 25 50       53 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
1999 25         39 my $x = shift;
2000              
2001 25         61 my $nrow = $x -> nrow();
2002 25         56 my $ncol = $x -> ncol();
2003              
2004 25 100       67 return 0 unless $nrow == $ncol;
2005              
2006 22         57 for my $i (1 .. $nrow - 1) {
2007 29         56 for my $j ($ncol - $i .. $ncol - 1) {
2008 34 100       149 return 0 unless $x -> [$i][$j] == 0;
2009             }
2010             }
2011              
2012 6         53 return 1;
2013             }
2014              
2015             =pod
2016              
2017             =item is_satriu()
2018              
2019             Returns 1 is the invocand is strictly upper anti-triangular, and 0 otherwise.
2020              
2021             $bool = $x -> is_satriu();
2022              
2023             A strictly anti-triangular matrix is a square matrix where all non-zero elements
2024             are strictly above the main diagonal. It has the following pattern, where only
2025             the elements marked as C can be non-zero,
2026              
2027             [ x x x 0 ]
2028             [ x x 0 0 ]
2029             [ x 0 0 0 ]
2030             [ 0 0 0 0 ]
2031              
2032             =cut
2033              
2034             sub is_satriu {
2035 25 50   25 1 87 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2036 25 50       50 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2037 25         42 my $x = shift;
2038              
2039 25         61 my $nrow = $x -> nrow();
2040 25         55 my $ncol = $x -> ncol();
2041              
2042 25 100       68 return 0 unless $nrow == $ncol;
2043              
2044 22         58 for my $i (0 .. $nrow - 1) {
2045 34         64 for my $j ($ncol - $i - 1 .. $ncol - 1) {
2046 42 100       158 return 0 unless $x -> [$i][$j] == 0;
2047             }
2048             }
2049              
2050 2         10 return 1;
2051             }
2052              
2053             =pod
2054              
2055             =item is_atril()
2056              
2057             Returns 1 is the invocand is lower anti-triangular, and 0 otherwise.
2058              
2059             $bool = $x -> is_atril();
2060              
2061             A lower anti-triangular matrix is a square matrix where all non-zero elements
2062             are on or below the main anti-diagonal. It has the following pattern, where only
2063             the elements marked as C can be non-zero,
2064              
2065             [ 0 0 0 x ]
2066             [ 0 0 x x ]
2067             [ 0 x x x ]
2068             [ x x x x ]
2069              
2070             =cut
2071              
2072             sub is_atril {
2073 25 50   25 1 80 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2074 25 50       54 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2075 25         43 my $x = shift;
2076              
2077 25         56 my $nrow = $x -> nrow();
2078 25         61 my $ncol = $x -> ncol();
2079              
2080 25 100       72 return 0 unless $nrow == $ncol;
2081              
2082 22         65 for my $i (0 .. $nrow - 1) {
2083 28         66 for my $j (0 .. $ncol - $i - 2) {
2084 39 100       192 return 0 unless $x -> [$i][$j] == 0;
2085             }
2086             }
2087              
2088 6         27 return 1;
2089             }
2090              
2091             =pod
2092              
2093             =item is_satril()
2094              
2095             Returns 1 is the invocand is strictly lower anti-triangular, and 0 otherwise.
2096              
2097             $bool = $x -> is_satril();
2098              
2099             A strictly lower anti-triangular matrix is a square matrix where all non-zero
2100             elements are strictly below the main anti-diagonal. It has the following
2101             pattern, where only the elements marked as C can be non-zero,
2102              
2103             [ 0 0 0 0 ]
2104             [ 0 0 0 x ]
2105             [ 0 0 x x ]
2106             [ 0 x x x ]
2107              
2108             =cut
2109              
2110             sub is_satril {
2111 25 50   25 1 89 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2112 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2113 25         44 my $x = shift;
2114              
2115 25         55 my $nrow = $x -> nrow();
2116 25         60 my $ncol = $x -> ncol();
2117              
2118 25 100       70 return 0 unless $nrow == $ncol;
2119              
2120 22         63 for my $i (0 .. $nrow - 1) {
2121 24         40 for my $j (0 .. $ncol - $i - 1) {
2122 45 100       178 return 0 unless $x -> [$i][$j] == 0;
2123             }
2124             }
2125              
2126 2         9 return 1;
2127             }
2128              
2129             =pod
2130              
2131             =back
2132              
2133             =head2 Identify elements
2134              
2135             This section contains methods for identifying and locating elements within an
2136             array. See also C>.
2137              
2138             =over 4
2139              
2140             =item find()
2141              
2142             Returns the location of each non-zero element.
2143              
2144             $K = $x -> find(); # linear index
2145             ($I, $J) = $x -> find(); # subscripts
2146              
2147             For example, to find the linear index of each element that is greater than or
2148             equal to 1, use
2149              
2150             $K = $x -> sge(1) -> find();
2151              
2152             =cut
2153              
2154             sub find {
2155 4 50   4 1 39 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2156 4 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2157 4         7 my $x = shift;
2158              
2159 4         12 my ($m, $n) = $x -> size();
2160              
2161 4         8 my $I = [];
2162 4         7 my $J = [];
2163 4         10 for my $j (0 .. $n - 1) {
2164 6         23 for my $i (0 .. $m - 1) {
2165 12 100       26 next unless $x->[$i][$j];
2166 10         14 push @$I, $i;
2167 10         26 push @$J, $j;
2168             }
2169             }
2170              
2171 4 100       21 return $I, $J if wantarray;
2172              
2173 2         6 my $K = [ map { $m * $J -> [$_] + $I -> [$_] } 0 .. $#$I ];
  5         12  
2174 2         9 return $K;
2175             }
2176              
2177             =pod
2178              
2179             =item is_finite()
2180              
2181             Returns a matrix of ones and zeros. The element is one if the corresponding
2182             element in the invocand matrix is finite, and zero otherwise.
2183              
2184             $y = $x -> is_finite();
2185              
2186             =cut
2187              
2188             sub is_finite {
2189 2 50   2 1 14 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2190 2 50       5 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2191 2         3 my $x = shift;
2192              
2193 2         12 require Math::Trig;
2194 2         6 my $pinf = Math::Trig::Inf(); # positiv infinity
2195 2         8 my $ninf = -$pinf; # negative infinity
2196              
2197 2         7 bless [ map { [
2198             map {
2199 2 100 100     4 $ninf < $_ && $_ < $pinf ? 1 : 0
  6         39  
2200             } @$_
2201             ] } @$x ], ref $x;
2202             }
2203              
2204             =pod
2205              
2206             =item is_inf()
2207              
2208             Returns a matrix of ones and zeros. The element is one if the corresponding
2209             element in the invocand matrix is positive or negative infinity, and zero
2210             otherwise.
2211              
2212             $y = $x -> is_inf();
2213              
2214             =cut
2215              
2216             sub is_inf {
2217 2 50   2 1 17 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2218 2 50       5 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2219 2         3 my $x = shift;
2220              
2221 2         12 require Math::Trig;
2222 2         9 my $pinf = Math::Trig::Inf(); # positiv infinity
2223 2         13 my $ninf = -$pinf; # negative infinity
2224              
2225 2         5 bless [ map { [
2226             map {
2227 2 100 100     5 $_ == $pinf || $_ == $ninf ? 1 : 0;
  6         27  
2228             } @$_
2229             ] } @$x ], ref $x;
2230             }
2231              
2232             =pod
2233              
2234             =item is_nan()
2235              
2236             Returns a matrix of ones and zeros. The element is one if the corresponding
2237             element in the invocand matrix is a NaN (Not-a-Number), and zero otherwise.
2238              
2239             $y = $x -> is_nan();
2240              
2241             =cut
2242              
2243             sub is_nan {
2244 2 50   2 1 13 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2245 2 50       4 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2246 2         5 my $x = shift;
2247              
2248 2 100       6 bless [ map [ map { $_ != $_ ? 1 : 0 } @$_ ], @$x ], ref $x;
  6         20  
2249             }
2250              
2251             =pod
2252              
2253             =item all()
2254              
2255             Tests whether all of the elements along various dimensions of a matrix are
2256             non-zero. If the dimension argument is not given, the first non-singleton
2257             dimension is used.
2258              
2259             $y = $x -> all($dim);
2260             $y = $x -> all();
2261              
2262             =cut
2263              
2264             sub all {
2265 12 50   12 1 92 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2266 12 50       27 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2267 12         17 my $x = shift;
2268 12         39 $x -> apply(\&_all, @_);
2269             }
2270              
2271             =pod
2272              
2273             =item any()
2274              
2275             Tests whether any of the elements along various dimensions of a matrix are
2276             non-zero. If the dimension argument is not given, the first non-singleton
2277             dimension is used.
2278              
2279             $y = $x -> any($dim);
2280             $y = $x -> any();
2281              
2282             =cut
2283              
2284             sub any {
2285 12 50   12 1 86 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2286 12 50       28 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2287 12         16 my $x = shift;
2288 12         37 $x -> apply(\&_any, @_);
2289             }
2290              
2291             =pod
2292              
2293             =item cumall()
2294              
2295             A cumulative variant of C>. If the dimension argument is not given,
2296             the first non-singleton dimension is used.
2297              
2298             $y = $x -> cumall($dim);
2299             $y = $x -> cumall();
2300              
2301             =cut
2302              
2303             sub cumall {
2304 12 50   12 1 92 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2305 12 50       32 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2306 12         19 my $x = shift;
2307 12         40 $x -> apply(\&_cumall, @_);
2308             }
2309              
2310             =pod
2311              
2312             =item cumany()
2313              
2314             A cumulative variant of C>. If the dimension argument is not given,
2315             the first non-singleton dimension is used.
2316              
2317             $y = $x -> cumany($dim);
2318             $y = $x -> cumany();
2319              
2320             =cut
2321              
2322             sub cumany {
2323 12 50   12 1 80 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2324 12 50       25 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2325 12         17 my $x = shift;
2326 12         39 $x -> apply(\&_cumany, @_);
2327             }
2328              
2329             =pod
2330              
2331             =back
2332              
2333             =head2 Basic properties
2334              
2335             =over 4
2336              
2337             =item size()
2338              
2339             You can determine the dimensions of a matrix by calling:
2340              
2341             ($m, $n) = $a -> size;
2342              
2343             =cut
2344              
2345             sub size {
2346 2087     2087 1 724442 my $self = shift;
2347 2087         2830 my $m = @{$self};
  2087         6091  
2348 2087 100       4253 my $n = $m ? @{$self->[0]} : 0;
  1918         3269  
2349 2087         5160 ($m, $n);
2350             }
2351              
2352             =pod
2353              
2354             =item nelm()
2355              
2356             Returns the number of elements in the matrix.
2357              
2358             $n = $x -> nelm();
2359              
2360             =cut
2361              
2362             sub nelm {
2363 193 50   193 1 1374 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2364 193 50       363 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2365 193         275 my $x = shift;
2366 193 100       665 return @$x ? @$x * @{$x->[0]} : 0;
  164         792  
2367             }
2368              
2369             =pod
2370              
2371             =item nrow()
2372              
2373             Returns the number of rows.
2374              
2375             $m = $x -> nrow();
2376              
2377             =cut
2378              
2379             sub nrow {
2380 657 50   657 1 7366 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2381 657 50       1155 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2382 657         900 my $x = shift;
2383 657         1610 return scalar @$x;
2384             }
2385              
2386             =pod
2387              
2388             =item ncol()
2389              
2390             Returns the number of columns.
2391              
2392             $n = $x -> ncol();
2393              
2394             =cut
2395              
2396             sub ncol {
2397 590 50   590 1 2229 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2398 590 50       1114 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2399 590         802 my $x = shift;
2400 590 100       1655 return @$x ? scalar(@{$x->[0]}) : 0;
  543         1349  
2401             }
2402              
2403             =pod
2404              
2405             =item npag()
2406              
2407             Returns the number of pages. A non-matrix has one page.
2408              
2409             $n = $x -> pag();
2410              
2411             =cut
2412              
2413             sub npag {
2414 6 50   6 1 941 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2415 6 50       14 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2416 6         12 my $x = shift;
2417 6 100       31 return @$x ? 1 : 0;
2418             }
2419              
2420             =pod
2421              
2422             =item ndim()
2423              
2424             Returns the number of dimensions. This is the number of dimensions along which
2425             the length is different from one.
2426              
2427             $n = $x -> ndim();
2428              
2429             =cut
2430              
2431             sub ndim {
2432 6 50   6 1 922 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2433 6 50       14 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2434 6         11 my $x = shift;
2435 6         13 my ($nrow, $ncol) = $x -> size();
2436 6         12 my $ndim = 0;
2437 6 100       12 ++$ndim if $nrow != 1;
2438 6 100       11 ++$ndim if $ncol != 1;
2439 6         26 return $ndim;
2440             }
2441              
2442             =pod
2443              
2444             =item bandwidth()
2445              
2446             Returns the bandwidth of a matrix. In scalar context, returns the number of the
2447             non-zero diagonal furthest away from the main diagonal. In list context,
2448             separate values are returned for the lower and upper bandwidth.
2449              
2450             $n = $x -> bandwidth();
2451             ($l, $u) = $x -> bandwidth();
2452              
2453             The bandwidth is a non-negative integer. If the bandwidth is 0, the matrix is
2454             diagonal or zero. If the bandwidth is 1, the matrix is tridiagonal. If the
2455             bandwidth is 2, the matrix is pentadiagonal etc.
2456              
2457             A matrix with the following pattern, where C denotes a non-zero value, would
2458             return 2 in scalar context, and (1,2) in list context.
2459              
2460             [ x x x 0 0 0 ]
2461             [ x x x x 0 0 ]
2462             [ 0 x x x x 0 ]
2463             [ 0 0 x x x x ]
2464             [ 0 0 0 x x x ]
2465             [ 0 0 0 0 x x ]
2466              
2467             See also C> and C>.
2468              
2469             =cut
2470              
2471             sub bandwidth {
2472 18 50   18 1 3800 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
2473 18 50       40 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
2474 18         26 my $x = shift;
2475              
2476 18         31 my ($nrow, $ncol) = $x -> size();
2477              
2478 18         24 my $upper = 0;
2479 18         25 my $lower = 0;
2480              
2481 18         30 for my $i (0 .. $nrow - 1) {
2482 52         61 for my $j (0 .. $ncol - 1) {
2483 192 100       253 next if $x -> [$i][$j] == 0;
2484 146         154 my $dist = $j - $i;
2485 146 100       164 if ($dist > 0) {
2486 52 100       84 $upper = $dist if $dist > $upper;
2487             } else {
2488 94 100       139 $lower = $dist if $dist < $lower;
2489             }
2490             }
2491             }
2492              
2493 18         24 $lower = -$lower;
2494 18 100       40 return $lower, $upper if wantarray;
2495 9 50       26 return $lower > $upper ? $lower : $upper;
2496             }
2497              
2498             =pod
2499              
2500             =back
2501              
2502             =head2 Manipulate matrices
2503              
2504             These methods are for combining matrices, splitting matrices, extracing parts of
2505             a matrix, inserting new parts into a matrix, deleting parts of a matrix etc.
2506             There are also methods for shuffling elements around (relocating elements)
2507             inside a matrix.
2508              
2509             These methods are not concerned with the values of the elements.
2510              
2511             =over 4
2512              
2513             =item catrow()
2514              
2515             Concatenate rows, i.e., concatenate matrices vertically. Any number of arguments
2516             is allowed. All non-empty matrices must have the same number or columns. The
2517             result is a new matrix.
2518              
2519             $x = Math::Matrix -> new([1, 2], [4, 5]); # 2-by-2 matrix
2520             $y = Math::Matrix -> new([3, 6]); # 1-by-2 matrix
2521             $z = $x -> catrow($y); # 3-by-2 matrix
2522              
2523             =cut
2524              
2525             sub catrow {
2526 25     25 1 174 my $x = shift;
2527 25         51 my $class = ref $x;
2528              
2529 25         32 my $ncol;
2530 25         47 my $z = bless [], $class; # initialize output
2531              
2532 25         56 for my $y ($x, @_) {
2533 63         145 my $ncoly = $y -> ncol();
2534 63 100       135 next if $ncoly == 0; # ignore empty $y
2535              
2536 44 100       83 if (defined $ncol) {
2537 22 50       62 croak "All operands must have the same number of columns in ",
2538             (caller(0))[3] unless $ncoly == $ncol;
2539             } else {
2540 22         28 $ncol = $ncoly;
2541             }
2542              
2543 44         142 push @$z, map [ @$_ ], @$y;
2544             }
2545              
2546 25         62 return $z;
2547             }
2548              
2549             =pod
2550              
2551             =item catcol()
2552              
2553             Concatenate columns, i.e., matrices horizontally. Any number of arguments is
2554             allowed. All non-empty matrices must have the same number or rows. The result is
2555             a new matrix.
2556              
2557             $x = Math::Matrix -> new([1, 2], [4, 5]); # 2-by-2 matrix
2558             $y = Math::Matrix -> new([3], [6]); # 2-by-1 matrix
2559             $z = $x -> catcol($y); # 2-by-3 matrix
2560              
2561             =cut
2562              
2563             sub catcol {
2564 92     92 1 243 my $x = shift;
2565 92         141 my $class = ref $x;
2566              
2567 92         115 my $nrow;
2568 92         180 my $z = bless [], $class; # initialize output
2569              
2570 92         199 for my $y ($x, @_) {
2571 191         356 my $nrowy = $y -> nrow();
2572 191 100       361 next if $nrowy == 0; # ignore empty $y
2573              
2574 176 100       308 if (defined $nrow) {
2575 87 50       174 croak "All operands must have the same number of rows in ",
2576             (caller(0))[3] unless $nrowy == $nrow;
2577             } else {
2578 89         126 $nrow = $nrowy;
2579             }
2580              
2581 176         316 for my $i (0 .. $nrow - 1) {
2582 419         501 push @{ $z -> [$i] }, @{ $y -> [$i] };
  419         613  
  419         799  
2583             }
2584             }
2585              
2586 92         166 return $z;
2587             }
2588              
2589             =pod
2590              
2591             =item getrow()
2592              
2593             Get the specified row(s). Returns a new matrix with the specified rows. The
2594             number of rows in the output is identical to the number of elements in the
2595             input.
2596              
2597             $y = $x -> getrow($i); # get one
2598             $y = $x -> getrow([$i0, $i1, $i2]); # get multiple
2599              
2600             =cut
2601              
2602             sub getrow {
2603 4 50   4 1 49 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
2604 4 50       13 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2605 4         7 my $x = shift;
2606 4         12 my $class = ref $x;
2607              
2608 4         8 my $idx = shift;
2609 4 50       9 croak "Row index can not be undefined" unless defined $idx;
2610 4 100       13 if (ref $idx) {
2611 3 50 33     21 $idx = __PACKAGE__ -> new($idx)
2612             unless defined(blessed($idx)) && $idx -> isa($class);
2613 3         12 $idx = $idx -> to_row();
2614 3         7 $idx = $idx -> [0];
2615             } else {
2616 1         5 $idx = [ $idx ];
2617             }
2618              
2619 4         10 my ($nrowx, $ncolx) = $x -> size();
2620              
2621 4         7 my $y = [];
2622 4         11 for my $iy (0 .. $#$idx) {
2623 5         8 my $ix = $idx -> [$iy];
2624 5 50       11 croak "Row index value $ix too large for $nrowx-by-$ncolx matrix in ",
2625             (caller(0))[3] if $ix >= $nrowx;
2626 5         8 $y -> [$iy] = [ @{ $x -> [$ix] } ];
  5         13  
2627             }
2628              
2629 4         12 bless $y, $class;
2630             }
2631              
2632             =pod
2633              
2634             =item getcol()
2635              
2636             Get the specified column(s). Returns a new matrix with the specified columns.
2637             The number of columns in the output is identical to the number of elements in
2638             the input.
2639              
2640             $y = $x -> getcol($j); # get one
2641             $y = $x -> getcol([$j0, $j1, $j2]); # get multiple
2642              
2643             =cut
2644              
2645             sub getcol {
2646 4 50   4 1 42 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
2647 4 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2648 4         6 my $x = shift;
2649 4         7 my $class = ref $x;
2650              
2651 4         6 my $idx = shift;
2652 4 50       7 croak "Column index can not be undefined" unless defined $idx;
2653 4 100       11 if (ref $idx) {
2654 3 50 33     17 $idx = __PACKAGE__ -> new($idx)
2655             unless defined(blessed($idx)) && $idx -> isa($class);
2656 3         9 $idx = $idx -> to_row();
2657 3         7 $idx = $idx -> [0];
2658             } else {
2659 1         2 $idx = [ $idx ];
2660             }
2661              
2662 4         11 my ($nrowx, $ncolx) = $x -> size();
2663              
2664 4         8 my $y = [];
2665 4         11 for my $jy (0 .. $#$idx) {
2666 5         7 my $jx = $idx -> [$jy];
2667 5 50       13 croak "Column index value $jx too large for $nrowx-by-$ncolx matrix in ",
2668             (caller(0))[3] if $jx >= $ncolx;
2669 5         11 for my $i (0 .. $nrowx - 1) {
2670 20         35 $y -> [$i][$jy] = $x -> [$i][$jx];
2671             }
2672             }
2673              
2674 4         14 bless $y, $class;
2675             }
2676              
2677             =pod
2678              
2679             =item delrow()
2680              
2681             Delete row(s). Returns a new matrix identical to the invocand but with the
2682             specified row(s) deleted.
2683              
2684             $y = $x -> delrow($i); # delete one
2685             $y = $x -> delrow([$i0, $i1, $i2]); # delete multiple
2686              
2687             =cut
2688              
2689             sub delrow {
2690 5 50   5 1 43 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
2691 5 50       12 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2692 5         6 my $x = shift;
2693 5         10 my $class = ref $x;
2694              
2695 5         6 my $idxdel = shift;
2696 5 50       10 croak "Row index can not be undefined" unless defined $idxdel;
2697 5 100       11 if (ref $idxdel) {
2698 4 50 33     19 $idxdel = __PACKAGE__ -> new($idxdel)
2699             unless defined(blessed($idxdel)) && $idxdel -> isa($class);
2700 4         12 $idxdel = $idxdel -> to_row();
2701 4         9 $idxdel = $idxdel -> [0];
2702             } else {
2703 1         3 $idxdel = [ $idxdel ];
2704             }
2705              
2706 5         11 my $nrowx = $x -> nrow();
2707              
2708             # This should be made faster.
2709              
2710 5         10 my $idxget = [];
2711 5         8 for my $i (0 .. $nrowx - 1) {
2712 20         27 my $seen = 0;
2713 20         26 for my $idx (@$idxdel) {
2714 28 100       55 if ($i == int $idx) {
2715 9         10 $seen = 1;
2716 9         11 last;
2717             }
2718             }
2719 20 100       42 push @$idxget, $i unless $seen;
2720             }
2721              
2722 5         7 my $y = [];
2723 5         19 @$y = map [ @$_ ], @$x[ @$idxget ];
2724 5         16 bless $y, $class;
2725             }
2726              
2727             =pod
2728              
2729             =item delcol()
2730              
2731             Delete column(s). Returns a new matrix identical to the invocand but with the
2732             specified column(s) deleted.
2733              
2734             $y = $x -> delcol($j); # delete one
2735             $y = $x -> delcol([$j0, $j1, $j2]); # delete multiple
2736              
2737             =cut
2738              
2739             sub delcol {
2740 5 50   5 1 42 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
2741 5 50       10 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
2742 5         9 my $x = shift;
2743 5         8 my $class = ref $x;
2744              
2745 5         8 my $idxdel = shift;
2746 5 50       9 croak "Column index can not be undefined" unless defined $idxdel;
2747 5 100       11 if (ref $idxdel) {
2748 4 50 33     21 $idxdel = __PACKAGE__ -> new($idxdel)
2749             unless defined(blessed($idxdel)) && $idxdel -> isa($class);
2750 4         12 $idxdel = $idxdel -> to_row();
2751 4         8 $idxdel = $idxdel -> [0];
2752             } else {
2753 1         2 $idxdel = [ $idxdel ];
2754             }
2755              
2756 5         13 my ($nrowx, $ncolx) = $x -> size();
2757              
2758             # This should be made faster.
2759              
2760 5         10 my $idxget = [];
2761 5         11 for my $j (0 .. $ncolx - 1) {
2762 20         28 my $seen = 0;
2763 20         27 for my $idx (@$idxdel) {
2764 28 100       54 if ($j == int $idx) {
2765 9         13 $seen = 1;
2766 9         11 last;
2767             }
2768             }
2769 20 100       39 push @$idxget, $j unless $seen;
2770             }
2771              
2772 5         6 my $y = [];
2773 5 100       12 if (@$idxget) {
2774 4         6 for my $row (@$x) {
2775 16         21 push @$y, [ @{$row}[ @$idxget ] ];
  16         29  
2776             }
2777             }
2778 5         18 bless $y, $class;
2779             }
2780              
2781             =pod
2782              
2783             =item concat()
2784              
2785             Concatenate two matrices horizontally. The matrices must have the same number of
2786             rows. The result is a new matrix or B in case of error.
2787              
2788             $x = Math::Matrix -> new([1, 2], [4, 5]); # 2-by-2 matrix
2789             $y = Math::Matrix -> new([3], [6]); # 2-by-1 matrix
2790             $z = $x -> concat($y); # 2-by-3 matrix
2791              
2792             =cut
2793              
2794             sub concat {
2795 11     11 1 63 my $self = shift;
2796 11         14 my $other = shift;
2797 11         36 my $result = $self->clone();
2798              
2799 11 50       17 return undef if scalar(@{$self}) != scalar(@{$other});
  11         16  
  11         38  
2800 11         19 for my $i (0 .. $#{$self}) {
  11         29  
2801 27         33 push @{$result->[$i]}, @{$other->[$i]};
  27         52  
  27         63  
2802             }
2803 11         31 $result;
2804             }
2805              
2806             =pod
2807              
2808             =item splicerow()
2809              
2810             Row splicing. This is like Perl's built-in splice() function, except that it
2811             works on the rows of a matrix.
2812              
2813             $y = $x -> splicerow($offset);
2814             $y = $x -> splicerow($offset, $length);
2815             $y = $x -> splicerow($offset, $length, $a, $b, ...);
2816              
2817             The built-in splice() function modifies the first argument and returns the
2818             removed elements, if any. However, since splicerow() does not modify the
2819             invocand, it returns the modified version as the first output argument and the
2820             removed part as a (possibly empty) second output argument.
2821              
2822             $x = Math::Matrix -> new([[ 1, 2],
2823             [ 3, 4],
2824             [ 5, 6],
2825             [ 7, 8]]);
2826             $a = Math::Matrix -> new([[11, 12],
2827             [13, 14]]);
2828             ($y, $z) = $x -> splicerow(1, 2, $a);
2829              
2830             Gives C<$y>
2831              
2832             [ 1 2 ]
2833             [ 11 12 ]
2834             [ 13 14 ]
2835             [ 7 8 ]
2836              
2837             and C<$z>
2838              
2839             [ 3 4 ]
2840             [ 5 6 ]
2841              
2842             =cut
2843              
2844             sub splicerow {
2845 12 50   12 1 100 croak "Not enough input arguments" if @_ < 1;
2846 12         18 my $x = shift;
2847 12         20 my $class = ref $x;
2848              
2849 12         17 my $offs = 0;
2850 12         30 my $len = $x -> nrow();
2851 12         29 my $repl = $class -> new([]);
2852              
2853 12 100       27 if (@_) {
2854 10         14 $offs = shift;
2855 10 50       26 croak "Offset can not be undefined" unless defined $offs;
2856 10 50       23 if (ref $offs) {
2857 0 0 0     0 $offs = $class -> new($offs)
2858             unless defined(blessed($offs)) && $offs -> isa($class);
2859 0 0       0 croak "Offset must be a scalar" unless $offs -> is_scalar();
2860 0         0 $offs = $offs -> [0][0];
2861             }
2862              
2863 10 100       21 if (@_) {
2864 4         9 $len = shift;
2865 4 50       10 croak "Length can not be undefined" unless defined $len;
2866 4 50       6 if (ref $len) {
2867 0 0 0     0 $len = $class -> new($len)
2868             unless defined(blessed($len)) && $len -> isa($class);
2869 0 0       0 croak "length must be a scalar" unless $len -> is_scalar();
2870 0         0 $len = $len -> [0][0];
2871             }
2872              
2873 4 100       9 if (@_) {
2874 2         7 $repl = $repl -> catrow(@_);
2875             }
2876             }
2877             }
2878              
2879 12         30 my $y = $x -> clone();
2880 12         24 my $z = $class -> new([]);
2881              
2882 12         48 @$z = splice @$y, $offs, $len, @$repl;
2883 12 100       75 return wantarray ? ($y, $z) : $y;
2884             }
2885              
2886             =pod
2887              
2888             =item splicecol()
2889              
2890             Column splicing. This is like Perl's built-in splice() function, except that it
2891             works on the columns of a matrix.
2892              
2893             $y = $x -> splicecol($offset);
2894             $y = $x -> splicecol($offset, $length);
2895             $y = $x -> splicecol($offset, $length, $a, $b, ...);
2896              
2897             The built-in splice() function modifies the first argument and returns the
2898             removed elements, if any. However, since splicecol() does not modify the
2899             invocand, it returns the modified version as the first output argument and the
2900             removed part as a (possibly empty) second output argument.
2901              
2902             $x = Math::Matrix -> new([[ 1, 3, 5, 7 ],
2903             [ 2, 4, 6, 8 ]]);
2904             $a = Math::Matrix -> new([[11, 13],
2905             [12, 14]]);
2906             ($y, $z) = $x -> splicerow(1, 2, $a);
2907              
2908             Gives C<$y>
2909              
2910             [ 1 11 13 7 ]
2911             [ 2 12 14 8 ]
2912              
2913             and C<$z>
2914              
2915             [ 3 5 ]
2916             [ 4 6 ]
2917              
2918             =cut
2919              
2920             sub splicecol {
2921 22 50   22 1 110 croak "Not enough input arguments" if @_ < 1;
2922 22         41 my $x = shift;
2923 22         40 my $class = ref $x;
2924              
2925 22         53 my ($nrowx, $ncolx) = $x -> size();
2926              
2927 22         45 my $offs = 0;
2928 22         43 my $len = $ncolx;
2929 22         60 my $repl = $class -> new([]);
2930              
2931 22 100       65 if (@_) {
2932 20         31 $offs = shift;
2933 20 50       49 croak "Offset can not be undefined" unless defined $offs;
2934 20 50       43 if (ref $offs) {
2935 0 0 0     0 $offs = $class -> new($offs)
2936             unless defined(blessed($offs)) && $offs -> isa($class);
2937 0 0       0 croak "Offset must be a scalar" unless $offs -> is_scalar();
2938 0         0 $offs = $offs -> [0][0];
2939             }
2940              
2941 20 100       57 if (@_) {
2942 14         27 $len = shift;
2943 14 50       29 croak "Length can not be undefined" unless defined $len;
2944 14 50       33 if (ref $len) {
2945 0 0 0     0 $len = $class -> new($len)
2946             unless defined(blessed($len)) && $len -> isa($class);
2947 0 0       0 croak "length must be a scalar" unless $len -> is_scalar();
2948 0         0 $len = $len -> [0][0];
2949             }
2950              
2951 14 100       37 if (@_) {
2952 2         7 $repl = $repl -> catcol(@_);
2953             }
2954             }
2955             }
2956              
2957 22         55 my $y = $x -> clone();
2958 22         52 my $z = $class -> new([]);
2959              
2960 22 50       60 if ($offs > $len) {
2961 0         0 carp "splicecol() offset past end of array";
2962 0         0 $offs = $len;
2963             }
2964              
2965             # The case when we are not removing anything from the invocand matrix: If
2966             # the offset is identical to the number of columns in the invocand matrix,
2967             # just appending the replacement matrix to the invocand matrix.
2968              
2969 22 100 100     102 if ($offs == $len) {
    100          
2970 2 50       5 unless ($repl -> is_empty()) {
2971 0         0 for my $i (0 .. $nrowx - 1) {
2972 0         0 push @{ $y -> [$i] }, @{ $repl -> [$i] };
  0         0  
  0         0  
2973             }
2974             }
2975             }
2976              
2977             # The case when we are removing everything from the invocand matrix: If the
2978             # offset is zero, and the length is identical to the number of columns in
2979             # the invocand matrix, replace the whole invocand matrix with the
2980             # replacement matrix.
2981              
2982             elsif ($offs == 0 && $len == $ncolx) {
2983 4         10 @$z = @$y;
2984 4         9 @$y = @$repl;
2985             }
2986              
2987             # The case when we are removing parts of the invocand matrix.
2988              
2989             else {
2990 16 100       52 if ($repl -> is_empty()) {
2991 14         35 for my $i (0 .. $nrowx - 1) {
2992 47         65 @{ $z -> [$i] } = splice @{ $y -> [$i] }, $offs, $len;
  47         111  
  47         79  
2993             }
2994             } else {
2995 2         5 for my $i (0 .. $nrowx - 1) {
2996 4         6 @{ $z -> [$i] } = splice @{ $y -> [$i] }, $offs, $len, @{ $repl -> [$i] };
  4         10  
  4         6  
  4         10  
2997             }
2998             }
2999             }
3000              
3001 22 100       106 return wantarray ? ($y, $z) : $y;
3002             }
3003              
3004             =pod
3005              
3006             =item swaprc()
3007              
3008             Swap rows and columns. This method does nothing but shuffle elements around. For
3009             real numbers, swaprc() is identical to the transpose() method.
3010              
3011             A subclass implementing a matrix of complex numbers should provide a transpose()
3012             method that also takes the complex conjugate of each elements. The swaprc()
3013             method, on the other hand, should only shuffle elements around.
3014              
3015             =cut
3016              
3017             sub swaprc {
3018 5     5 1 35 my $x = shift;
3019 5         10 my $class = ref $x;
3020              
3021 5         11 my $y = bless [], $class;
3022 5         12 my $ncolx = $x -> ncol();
3023 5 100       15 return $y if $ncolx == 0;
3024              
3025 4         13 for my $j (0 .. $ncolx - 1) {
3026 9         30 push @$y, [ map $_->[$j], @$x ];
3027             }
3028 4         11 return $y;
3029             }
3030              
3031             =pod
3032              
3033             =item flipud()
3034              
3035             Flip upside-down, i.e., flip along dimension 1.
3036              
3037             $y = $x -> flipud();
3038              
3039             =cut
3040              
3041             sub flipud {
3042 2 50   2 1 28 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3043 2 50       6 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3044 2         4 my $x = shift;
3045 2         5 my $class = ref $x;
3046              
3047 2         51 my $y = [ reverse map [ @$_ ], @$x ];
3048 2         6 bless $y, $class;;
3049             }
3050              
3051             =pod
3052              
3053             =item fliplr()
3054              
3055             Flip left-to-right, i.e., flip along dimension 2.
3056              
3057             $y = $x -> fliplr();
3058              
3059             =cut
3060              
3061             sub fliplr {
3062 25 50   25 1 105 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3063 25 50       57 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3064 25         42 my $x = shift;
3065 25         47 my $class = ref $x;
3066              
3067 25         173 my $y = [ map [ reverse @$_ ], @$x ];
3068 25         85 bless $y, $class;
3069             }
3070              
3071             =pod
3072              
3073             =item flip()
3074              
3075             Flip along various dimensions of a matrix. If the dimension argument is not
3076             given, the first non-singleton dimension is used.
3077              
3078             $y = $x -> flip($dim);
3079             $y = $x -> flip();
3080              
3081             See also C> and C>.
3082              
3083             =cut
3084              
3085             sub flip {
3086 12 50   12 1 74 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3087 12 50       21 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3088 12         19 my $x = shift;
3089 12     20   57 $x -> apply(sub { reverse @_ }, @_);
  20         37  
3090             }
3091              
3092             =pod
3093              
3094             =item rot90()
3095              
3096             Rotate 90 degrees counterclockwise.
3097              
3098             $y = $x -> rot90(); # rotate 90 degrees counterclockwise
3099             $y = $x -> rot90($n); # rotate 90*$n degrees counterclockwise
3100              
3101             =cut
3102              
3103             sub rot90 {
3104 12 50   12 1 63 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3105 12 50       24 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3106 12         20 my $x = shift;
3107 12         20 my $class = ref $x;
3108              
3109 12         22 my $n = 1;
3110 12 100       30 if (@_) {
3111 10         14 $n = shift;
3112 10 100       23 if (ref $n) {
3113 2 100 66     24 $n = $class -> new($n)
3114             unless defined(blessed($n)) && $n -> isa($class);
3115 2 50       7 croak "Argument must be a scalar" unless $n -> is_scalar();
3116 2         6 $n = $n -> [0][0];
3117             }
3118 10 50       25 croak "Argument must be an integer" unless $n == int $n;
3119             }
3120              
3121 12         20 my $y = [];
3122              
3123             # Rotate 0 degrees, i.e., clone.
3124              
3125 12         27 $n %= 4;
3126 12 100       46 if ($n == 0) {
    100          
    100          
    50          
3127 1         5 $y = [ map [ @$_ ], @$x ];
3128             }
3129              
3130             # Rotate 90 degrees.
3131              
3132             elsif ($n == 1) {
3133 5         12 my ($nrowx, $ncolx) = $x -> size();
3134 5         11 my $jmax = $ncolx - 1;
3135 5         10 for my $i (0 .. $nrowx - 1) {
3136 8         15 for my $j (0 .. $ncolx - 1) {
3137 24         46 $y -> [$jmax - $j][$i] = $x -> [$i][$j];
3138             }
3139             }
3140             }
3141              
3142             # Rotate 180 degrees.
3143              
3144             elsif ($n == 2) {
3145 3         43 $y = [ map [ reverse @$_ ], reverse @$x ];
3146             }
3147              
3148             # Rotate 270 degrees.
3149              
3150             elsif ($n == 3) {
3151 3         11 my ($nrowx, $ncolx) = $x -> size();
3152 3         6 my $imax = $nrowx - 1;
3153 3         8 for my $i (0 .. $nrowx - 1) {
3154 4         7 for my $j (0 .. $ncolx - 1) {
3155 12         29 $y -> [$j][$imax - $i] = $x -> [$i][$j];
3156             }
3157             }
3158             }
3159              
3160 12         40 bless $y, $class;
3161             }
3162              
3163             =pod
3164              
3165             =item rot180()
3166              
3167             Rotate 180 degrees.
3168              
3169             $y = $x -> rot180();
3170              
3171             =cut
3172              
3173             sub rot180 {
3174 2 50   2 1 23 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3175 2 50       6 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3176 2         4 my $x = shift;
3177 2         5 $x -> rot90(2);
3178             }
3179              
3180             =pod
3181              
3182             =item rot270()
3183              
3184             Rotate 270 degrees counterclockwise, i.e., 90 degrees clockwise.
3185              
3186             $y = $x -> rot270();
3187              
3188             =cut
3189              
3190             sub rot270 {
3191 2 50   2 1 26 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3192 2 50       7 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3193 2         3 my $x = shift;
3194 2         7 $x -> rot90(3);
3195             }
3196              
3197             =pod
3198              
3199             =item repelm()
3200              
3201             Repeat elements.
3202              
3203             $x -> repelm($y);
3204              
3205             Repeats each element in $x the number of times specified in $y.
3206              
3207             If $x is the matrix
3208              
3209             [ 4 5 6 ]
3210             [ 7 8 9 ]
3211              
3212             and $y is
3213              
3214             [ 3 2 ]
3215              
3216             the returned matrix is
3217              
3218             [ 4 4 5 5 6 6 ]
3219             [ 4 4 5 5 6 6 ]
3220             [ 4 4 5 5 6 6 ]
3221             [ 7 7 8 8 9 9 ]
3222             [ 7 7 8 8 9 9 ]
3223             [ 7 7 8 8 9 9 ]
3224              
3225             =cut
3226              
3227             sub repelm {
3228 5 50   5 1 36 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3229 5 50       11 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3230 5         7 my $x = shift;
3231 5         8 my $class = ref $x;
3232              
3233 5         8 my $y = shift;
3234 5 50 33     36 $y = __PACKAGE__ -> new($y)
3235             unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
3236 5 50       17 croak "Input argument must contain two elements"
3237             unless $y -> nelm() == 2;
3238              
3239 5         12 my ($nrowx, $ncolx) = $x -> size();
3240              
3241 5         12 $y = $y -> to_col();
3242 5         10 my $nrowrep = $y -> [0][0];
3243 5         7 my $ncolrep = $y -> [1][0];
3244              
3245 5         8 my $z = [];
3246 5         10 for my $ix (0 .. $nrowx - 1) {
3247 8         14 for my $jx (0 .. $ncolx - 1) {
3248 24         33 for my $iy (0 .. $nrowrep - 1) {
3249 30         46 for my $jy (0 .. $ncolrep - 1) {
3250 36         51 my $iz = $ix * $nrowrep + $iy;
3251 36         38 my $jz = $jx * $ncolrep + $jy;
3252 36         58 $z -> [$iz][$jz] = $x -> [$ix][$jx];
3253             }
3254             }
3255             }
3256             }
3257              
3258 5         20 bless $z, $class;
3259             }
3260              
3261             =pod
3262              
3263             =item repmat()
3264              
3265             Repeat elements.
3266              
3267             $x -> repmat($y);
3268              
3269             Repeats the matrix $x the number of times specified in $y.
3270              
3271             If $x is the matrix
3272              
3273             [ 4 5 6 ]
3274             [ 7 8 9 ]
3275              
3276             and $y is
3277              
3278             [ 3 2 ]
3279              
3280             the returned matrix is
3281              
3282             [ 4 5 6 4 5 6 ]
3283             [ 7 8 9 7 8 9 ]
3284             [ 4 5 6 4 5 6 ]
3285             [ 7 8 9 7 8 9 ]
3286             [ 4 5 6 4 5 6 ]
3287             [ 7 8 9 7 8 9 ]
3288              
3289             =cut
3290              
3291             sub repmat {
3292 5 50   5 1 50 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3293 5 50       13 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3294 5         9 my $x = shift;
3295 5         7 my $class = ref $x;
3296              
3297 5         9 my $y = shift;
3298 5 50 33     38 $y = __PACKAGE__ -> new($y)
3299             unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
3300 5 50       20 croak "Input argument must contain two elements"
3301             unless $y -> nelm() == 2;
3302              
3303 5         11 my ($nrowx, $ncolx) = $x -> size();
3304              
3305 5         13 $y = $y -> to_col();
3306 5         9 my $nrowrep = $y -> [0][0];
3307 5         6 my $ncolrep = $y -> [1][0];
3308              
3309 5         8 my $z = [];
3310 5         8 for my $ix (0 .. $nrowx - 1) {
3311 8         16 for my $jx (0 .. $ncolx - 1) {
3312 24         40 for my $iy (0 .. $nrowrep - 1) {
3313 30         44 for my $jy (0 .. $ncolrep - 1) {
3314 36         47 my $iz = $iy * $nrowx + $ix;
3315 36         43 my $jz = $jy * $ncolx + $jx;
3316 36         77 $z -> [$iz][$jz] = $x -> [$ix][$jx];
3317             }
3318             }
3319             }
3320             }
3321              
3322 5         23 bless $z, $class;
3323             }
3324              
3325             =pod
3326              
3327             =item reshape()
3328              
3329             Returns a reshaped copy of a matrix. The reshaping is done by creating a new
3330             matrix and looping over the elements in column major order. The new matrix must
3331             have the same number of elements as the invocand matrix. The following returns
3332             an C<$m>-by-C<$n> matrix,
3333              
3334             $y = $x -> reshape($m, $n);
3335              
3336             The code
3337              
3338             $x = Math::Matrix -> new([[1, 3, 5, 7], [2, 4, 6, 8]]);
3339             $y = $x -> reshape(4, 2);
3340              
3341             creates the matrix $x
3342              
3343             [ 1 3 5 7 ]
3344             [ 2 4 6 8 ]
3345              
3346             and returns a reshaped copy $y
3347              
3348             [ 1 5 ]
3349             [ 2 6 ]
3350             [ 3 7 ]
3351             [ 4 8 ]
3352              
3353             =cut
3354              
3355             sub reshape {
3356 30 50   30 1 184 croak "Not enough arguments for ", (caller(0))[3] if @_ < 3;
3357 30 50       65 croak "Too many arguments for ", (caller(0))[3] if @_ > 3;
3358 30         45 my $x = shift;
3359 30         54 my $class = ref $x;
3360              
3361 30         85 my ($nrowx, $ncolx) = $x -> size();
3362 30         57 my $nelmx = $nrowx * $ncolx;
3363              
3364 30         48 my ($nrowy, $ncoly) = @_;
3365 30         47 my $nelmy = $nrowy * $ncoly;
3366              
3367 30 50       60 croak "when reshaping, the number of elements can not change in ",
3368             (caller(0))[3] unless $nelmx == $nelmy;
3369              
3370 30         48 my $y = [];
3371              
3372             # No reshaping; just clone.
3373              
3374 30 100 66     125 if ($nrowx == $nrowy && $ncolx == $ncoly) {
    100          
    100          
3375 5         21 $y = [ map [ @$_ ], @$x ];
3376             }
3377              
3378             elsif ($nrowx == 1) {
3379              
3380             # Reshape from a row vector to a column vector.
3381              
3382 10 100       27 if ($ncoly == 1) {
3383 4         14 $y = [ map [ $_ ], @{ $x -> [0] } ];
  4         32  
3384             }
3385              
3386             # Reshape from a row vector to a matrix.
3387              
3388             else {
3389 6         14 my $k = 0;
3390 6         19 for my $j (0 .. $ncoly - 1) {
3391 15         32 for my $i (0 .. $nrowy - 1) {
3392 33         62 $y -> [$i][$j] = $x -> [0][$k++];
3393             }
3394             }
3395             }
3396             }
3397              
3398             elsif ($ncolx == 1) {
3399              
3400             # Reshape from a column vector to a row vector.
3401              
3402 6 100       23 if ($nrowy == 1) {
3403 3         12 $y = [[ map { @$_ } @$x ]];
  18         31  
3404             }
3405              
3406             # Reshape from a column vector to a matrix.
3407              
3408             else {
3409 3         10 my $k = 0;
3410 3         13 for my $j (0 .. $ncoly - 1) {
3411 9         21 for my $i (0 .. $nrowy - 1) {
3412 18         40 $y -> [$i][$j] = $x -> [$k++][0];
3413             }
3414             }
3415             }
3416             }
3417              
3418             # The invocand is a matrix. This code works in all cases, but is somewhat
3419             # slower than the specialized code above.
3420              
3421             else {
3422 9         25 for my $k (0 .. $nelmx - 1) {
3423 54         74 my $ix = $k % $nrowx;
3424 54         73 my $jx = ($k - $ix) / $nrowx;
3425 54         79 my $iy = $k % $nrowy;
3426 54         78 my $jy = ($k - $iy) / $nrowy;
3427 54         113 $y -> [$iy][$jy] = $x -> [$ix][$jx];
3428             }
3429             }
3430              
3431 30         104 bless $y, $class;
3432             }
3433              
3434             =pod
3435              
3436             =item to_row()
3437              
3438             Reshape to a row.
3439              
3440             $x -> to_row();
3441              
3442             This method reshapes the matrix into a single row. It is essentially the same
3443             as, but faster than,
3444              
3445             $x -> reshape(1, $x -> nelm());
3446              
3447             =cut
3448              
3449             sub to_row {
3450 30 50   30 1 127 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3451 30 50       76 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3452 30         47 my $x = shift;
3453 30         60 my $class = ref $x;
3454              
3455 30         60 my $y = bless [], $class;
3456              
3457 30         76 my $ncolx = $x -> ncol();
3458 30 100       81 return $y if $ncolx == 0;
3459              
3460 25         69 for my $j (0 .. $ncolx - 1) {
3461 61         86 push @{ $y -> [0] }, map $_->[$j], @$x;
  61         210  
3462             }
3463 25         109 return $y;
3464             }
3465              
3466             =pod
3467              
3468             =item to_col()
3469              
3470             Reshape to a column.
3471              
3472             $y = $x -> to_col();
3473              
3474             This method reshapes the matrix into a single column. It is essentially the same
3475             as, but faster than,
3476              
3477             $x -> reshape($x -> nelm(), 1);
3478              
3479             =cut
3480              
3481             sub to_col {
3482 19 50   19 1 69 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3483 19 50       45 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3484 19         34 my $x = shift;
3485              
3486 19         47 my $class = ref $x;
3487              
3488 19         41 my $y = bless [], $class;
3489              
3490 19         47 my $ncolx = $x -> ncol();
3491 19 100       54 return $y if $ncolx == 0;
3492              
3493 18         49 for my $j (0 .. $ncolx - 1) {
3494 43         129 push @$y, map [ $_->[$j] ], @$x;
3495             }
3496 18         41 return $y;
3497             }
3498              
3499             =pod
3500              
3501             =item to_permmat()
3502              
3503             Permutation vector to permutation matrix. Converts a vector of zero-based
3504             permutation indices to a permutation matrix.
3505              
3506             $P = $v -> to_permmat();
3507              
3508             For example
3509              
3510             $v = Math::Matrix -> new([[0, 3, 1, 4, 2]]);
3511             $m = $v -> to_permmat();
3512              
3513             gives the permutation matrix C<$m>
3514              
3515             [ 1 0 0 0 0 ]
3516             [ 0 0 0 1 0 ]
3517             [ 0 1 0 0 0 ]
3518             [ 0 0 0 0 1 ]
3519             [ 0 0 1 0 0 ]
3520              
3521             =cut
3522              
3523             sub to_permmat {
3524 4 50   4 1 22 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3525 4 50       20 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3526 4         6 my $v = shift;
3527 4         8 my $class = ref $v;
3528              
3529 4         9 my $n = $v -> nelm();
3530 4         10 my $P = $class -> zeros($n, $n); # initialize output
3531 4 100       11 return $P if $n == 0; # if emtpy $v
3532              
3533 3 50       7 croak "Invocand must be a vector" unless $v -> is_vector();
3534 3         7 $v = $v -> to_col();
3535              
3536 3         5 for my $i (0 .. $n - 1) {
3537 9         14 my $j = $v -> [$i][0];
3538 9 50 33     27 croak "index out of range" unless 0 <= $j && $j < $n;
3539 9         16 $P -> [$i][$j] = 1;
3540             }
3541              
3542 3         16 return $P;
3543             }
3544              
3545             =pod
3546              
3547             =item to_permvec()
3548              
3549             Permutation matrix to permutation vector. Converts a permutation matrix to a
3550             vector of zero-based permutation indices.
3551              
3552             $v = $P -> to_permvec();
3553              
3554             $v = Math::Matrix -> new([[0, 3, 1, 4, 2]]);
3555             $m = $v -> to_permmat();
3556              
3557             Gives the permutation matrix C<$m>
3558              
3559             [ 1 0 0 0 0 ]
3560             [ 0 0 0 1 0 ]
3561             [ 0 1 0 0 0 ]
3562             [ 0 0 0 0 1 ]
3563             [ 0 0 1 0 0 ]
3564              
3565             See also C>.
3566              
3567             =cut
3568              
3569             sub to_permvec {
3570 4 50   4 1 37 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3571 4 50       13 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3572 4         8 my $P = shift;
3573 4         7 my $class = ref $P;
3574              
3575 4 50       12 croak "Invocand matrix must be square" unless $P -> is_square();
3576 4         14 my $n = $P -> nrow();
3577              
3578 4         14 my $v = $class -> zeros($n, 1); # initialize output
3579              
3580 4         11 my $seen = [ (0) x $n ]; # keep track of the ones
3581              
3582 4         11 for my $i (0 .. $n - 1) {
3583 9         11 my $k;
3584 9         16 for my $j (0 .. $n - 1) {
3585 27 100       52 next if $P -> [$i][$j] == 0;
3586 9 50       18 if ($P -> [$i][$j] == 1) {
3587 9 50       16 croak "invalid permutation matrix; more than one row has",
3588             " an element with value 1 in column $j" if $seen->[$j]++;
3589 9         12 $k = $j;
3590 9         15 next;
3591             }
3592 0         0 croak "invalid permutation matrix; element ($i,$j)",
3593             " is neither 0 nor 1";
3594             }
3595 9 50       17 croak "invalid permutation matrix; row $i has no element with value 1"
3596             unless defined $k;
3597 9         14 $v->[$i][0] = $k;
3598             }
3599              
3600 4         13 return $v;
3601             }
3602              
3603             =pod
3604              
3605             =item triu()
3606              
3607             Upper triangular part. Extract the upper triangular part of a matrix and set all
3608             other elements to zero.
3609              
3610             $y = $x -> triu();
3611             $y = $x -> triu($n);
3612              
3613             The optional second argument specifies how many diagonals above or below the
3614             main diagonal should also be set to zero. The default value of C<$n> is zero
3615             which includes the main diagonal.
3616              
3617             =cut
3618              
3619             sub triu {
3620 12 50   12 1 89 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3621 12 50       26 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3622 12         20 my $x = shift;
3623 12         20 my $class = ref $x;
3624              
3625 12         20 my $n = 0;
3626 12 100       26 if (@_) {
3627 10         14 $n = shift;
3628 10 50       19 if (ref $n) {
3629 0 0 0     0 $n = $class -> new($n)
3630             unless defined(blessed($n)) && $n -> isa($class);
3631 0 0       0 croak "Argument must be a scalar" unless $n -> is_scalar();
3632 0         0 $n = $n -> [0][0];
3633             }
3634 10 50       23 croak "Argument must be an integer" unless $n == int $n;
3635             }
3636              
3637 12         29 my ($nrowx, $ncolx) = $x -> size();
3638              
3639 12         21 my $y = [];
3640 12         25 for my $i (0 .. $nrowx - 1) {
3641 30         50 for my $j (0 .. $ncolx - 1) {
3642 72 100       149 $y -> [$i][$j] = $j - $i >= $n ? $x -> [$i][$j] : 0;
3643             }
3644             }
3645              
3646 12         31 bless $y, $class;
3647             }
3648              
3649             =pod
3650              
3651             =item tril()
3652              
3653             Lower triangular part. Extract the lower triangular part of a matrix and set all
3654             other elements to zero.
3655              
3656             $y = $x -> tril();
3657             $y = $x -> tril($n);
3658              
3659             The optional second argument specifies how many diagonals above or below the
3660             main diagonal should also be set to zero. The default value of C<$n> is zero
3661             which includes the main diagonal.
3662              
3663             =cut
3664              
3665             sub tril {
3666 12 50   12 1 74 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3667 12 50       27 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3668 12         20 my $x = shift;
3669 12         20 my $class = ref $x;
3670              
3671 12         19 my $n = 0;
3672 12 100       22 if (@_) {
3673 10         14 $n = shift;
3674 10 50       21 if (ref $n) {
3675 0 0 0     0 $n = $class -> new($n)
3676             unless defined(blessed($n)) && $n -> isa($class);
3677 0 0       0 croak "Argument must be a scalar" unless $n -> is_scalar();
3678 0         0 $n = $n -> [0][0];
3679             }
3680 10 50       21 croak "Argument must be an integer" unless $n == int $n;
3681             }
3682              
3683 12         27 my ($nrowx, $ncolx) = $x -> size();
3684              
3685 12         22 my $y = [];
3686 12         26 for my $i (0 .. $nrowx - 1) {
3687 30         44 for my $j (0 .. $ncolx - 1) {
3688 72 100       154 $y -> [$i][$j] = $j - $i <= $n ? $x -> [$i][$j] : 0;
3689             }
3690             }
3691              
3692 12         36 bless $y, $class;
3693             }
3694              
3695             =pod
3696              
3697             =item slice()
3698              
3699             Extract columns:
3700              
3701             a->slice(1,3,5);
3702              
3703             =cut
3704              
3705             sub slice {
3706 16     16 1 66 my $self = shift;
3707 16         31 my $class = ref($self);
3708 16         37 my $result = [];
3709              
3710 16         75 for my $i (0 .. $#$self) {
3711 39         59 push @$result, [ @{$self->[$i]}[@_] ];
  39         90  
3712             }
3713              
3714 16         58 bless $result, $class;
3715             }
3716              
3717             =pod
3718              
3719             =item diagonal_vector()
3720              
3721             Extract the diagonal as an array:
3722              
3723             $diag = $a->diagonal_vector;
3724              
3725             =cut
3726              
3727             sub diagonal_vector {
3728 1     1 1 2029 my $self = shift;
3729 1         3 my @diag;
3730 1         2 my $idx = 0;
3731 1         6 my($m, $n) = $self->size();
3732              
3733 1 50       4 croak "Not a square matrix" if $m != $n;
3734              
3735 1         2 foreach my $r (@{$self}) {
  1         3  
3736 4         9 push @diag, $r->[$idx++];
3737             }
3738 1         4 return \@diag;
3739             }
3740              
3741             =pod
3742              
3743             =item tridiagonal_vector()
3744              
3745             Extract the diagonals that make up a tridiagonal matrix:
3746              
3747             ($main_d, $upper_d, $lower_d) = $a->tridiagonal_vector;
3748              
3749             =cut
3750              
3751             sub tridiagonal_vector {
3752 1     1 1 1455 my $self = shift;
3753 1         3 my(@main_d, @up_d, @low_d);
3754 1         4 my($m, $n) = $self->size();
3755 1         2 my $idx = 0;
3756              
3757 1 50       5 croak "Not a square matrix" if $m != $n;
3758              
3759 1         2 foreach my $r (@{$self}) {
  1         3  
3760 4 100       10 push @low_d, $r->[$idx - 1] if ($idx > 0);
3761 4         7 push @main_d, $r->[$idx++];
3762 4 100       11 push @up_d, $r->[$idx] if ($idx < $m);
3763             }
3764 1         6 return ([@main_d],[@up_d],[@low_d]);
3765             }
3766              
3767             =pod
3768              
3769             =back
3770              
3771             =head2 Mathematical functions
3772              
3773             =head3 Addition
3774              
3775             =over 4
3776              
3777             =item add()
3778              
3779             Addition. If one operands is a scalar, it is treated like a constant matrix with
3780             the same size as the other operand. Otherwise ordinary matrix addition is
3781             performed.
3782              
3783             $z = $x -> add($y);
3784              
3785             See also C> and C>.
3786              
3787             =cut
3788              
3789             sub add {
3790 8 50   8 1 66 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3791 8 50       23 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3792 8         12 my $x = shift;
3793 8         15 my $class = ref $x;
3794              
3795 8         10 my $y = shift;
3796 8 100 66     75 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
3797              
3798 8 100 100     23 $x -> is_scalar() || $y -> is_scalar() ? $x -> sadd($y) : $x -> madd($y);
3799             }
3800              
3801             =pod
3802              
3803             =item madd()
3804              
3805             Matrix addition. Add two matrices of the same dimensions. An error is thrown if
3806             the matrices don't have the same size.
3807              
3808             $z = $x -> madd($y);
3809              
3810             See also C> and C>.
3811              
3812             =cut
3813              
3814             sub madd {
3815 5 50   5 1 32 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3816 5 50       14 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3817 5         7 my $x = shift;
3818 5         11 my $class = ref $x;
3819              
3820 5         8 my $y = shift;
3821 5 50 33     50 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
3822              
3823 5         20 my ($nrowx, $ncolx) = $x -> size();
3824 5         10 my ($nrowy, $ncoly) = $y -> size();
3825              
3826 5 50 33     26 croak "Can't add $nrowx-by-$ncolx matrix to $nrowy-by-$ncoly matrix"
3827             unless $nrowx == $nrowy && $ncolx == $ncoly;
3828              
3829 5         8 my $z = [];
3830 5         18 for my $i (0 .. $nrowx - 1) {
3831 7         13 for my $j (0 .. $ncolx - 1) {
3832 21         45 $z->[$i][$j] = $x->[$i][$j] + $y->[$i][$j];
3833             }
3834             }
3835              
3836 5         15 bless $z, $class;
3837             }
3838              
3839             =pod
3840              
3841             =item sadd()
3842              
3843             Scalar (element by element) addition with scalar expansion. This method places
3844             no requirements on the size of the input matrices.
3845              
3846             $z = $x -> sadd($y);
3847              
3848             See also C> and C>.
3849              
3850             =cut
3851              
3852             sub sadd {
3853 9 50   9 1 46 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3854 9 50       20 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3855 9         13 my $x = shift;
3856              
3857 9     33   37 my $sub = sub { $_[0] + $_[1] };
  33         67  
3858 9         29 $x -> sapply($sub, @_);
3859             }
3860              
3861             =pod
3862              
3863             =back
3864              
3865             =head3 Subtraction
3866              
3867             =over 4
3868              
3869             =item sub()
3870              
3871             Subtraction. If one operands is a scalar, it is treated as a constant matrix
3872             with the same size as the other operand. Otherwise, ordinarly matrix subtraction
3873             is performed.
3874              
3875             $z = $x -> sub($y);
3876              
3877             See also C> and C>.
3878              
3879             =cut
3880              
3881             sub sub {
3882 9 50   9 1 53 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3883 9 50       22 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3884 9         14 my $x = shift;
3885 9         17 my $class = ref $x;
3886              
3887 9         22 my $y = shift;
3888 9 100 66     72 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
3889              
3890 9 100 100     28 $x -> is_scalar() || $y -> is_scalar() ? $x -> ssub($y) : $x -> msub($y);
3891             }
3892              
3893             =pod
3894              
3895             =item msub()
3896              
3897             Matrix subtraction. Subtract two matrices of the same size. An error is thrown
3898             if the matrices don't have the same size.
3899              
3900             $z = $x -> msub($y);
3901              
3902             See also C> and C>.
3903              
3904             =cut
3905              
3906             sub msub {
3907 6 50   6 1 40 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3908 6 50       29 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3909 6         12 my $x = shift;
3910 6         21 my $class = ref $x;
3911              
3912 6         42 my $y = shift;
3913 6 50 33     66 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
3914              
3915 6         32 my ($nrowx, $ncolx) = $x -> size();
3916 6         16 my ($nrowy, $ncoly) = $y -> size();
3917              
3918 6 50 33     30 croak "Can't subtract $nrowy-by-$ncoly matrix from $nrowx-by-$ncolx matrix"
3919             unless $nrowx == $nrowy && $ncolx == $ncoly;
3920              
3921 6         13 my $z = [];
3922 6         20 for my $i (0 .. $nrowx - 1) {
3923 8         17 for my $j (0 .. $ncolx - 1) {
3924 24         51 $z->[$i][$j] = $x->[$i][$j] - $y->[$i][$j];
3925             }
3926             }
3927              
3928 6         31 bless $z, $class;
3929             }
3930              
3931             =pod
3932              
3933             =item ssub()
3934              
3935             Scalar (element by element) subtraction with scalar expansion. This method
3936             places no requirements on the size of the input matrices.
3937              
3938             $z = $x -> ssub($y);
3939              
3940             See also C> and C>.
3941              
3942             =cut
3943              
3944             sub ssub {
3945 9 50   9 1 44 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
3946 9 50       20 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
3947 9         15 my $x = shift;
3948              
3949 9     33   34 my $sub = sub { $_[0] - $_[1] };
  33         73  
3950 9         29 $x -> sapply($sub, @_);
3951             }
3952              
3953             =pod
3954              
3955             =item subtract()
3956              
3957             This is an alias for C>.
3958              
3959             =cut
3960              
3961             sub subtract {
3962 1     1 1 7 my $x = shift;
3963 1         5 $x -> sub(@_);
3964             }
3965              
3966             =pod
3967              
3968             =back
3969              
3970             =head3 Negation
3971              
3972             =over 4
3973              
3974             =item neg()
3975              
3976             Negation. Negate a matrix.
3977              
3978             $y = $x -> neg();
3979              
3980             It is effectively equivalent to
3981              
3982             $y = $x -> map(sub { -$_ });
3983              
3984             =cut
3985              
3986             sub neg {
3987 4 50   4 1 34 croak "Not enough arguments for ", (caller(0))[3] if @_ < 1;
3988 4 50       11 croak "Too many arguments for ", (caller(0))[3] if @_ > 1;
3989 4         8 my $x = shift;
3990 4         71 bless [ map [ map -$_, @$_ ], @$x ], ref $x;
3991             }
3992              
3993             =pod
3994              
3995             =item negative()
3996              
3997             This is an alias for C>.
3998              
3999             =cut
4000              
4001             sub negative {
4002 0     0 1 0 my $x = shift;
4003 0         0 $x -> neg(@_);
4004             }
4005              
4006             =pod
4007              
4008             =back
4009              
4010             =head3 Multiplication
4011              
4012             =over 4
4013              
4014             =item mul()
4015              
4016             Multiplication. If one operands is a scalar, it is treated as a constant matrix
4017             with the same size as the other operand. Otherwise, ordinary matrix
4018             multiplication is performed.
4019              
4020             $z = $x -> mul($y);
4021              
4022             =cut
4023              
4024             sub mul {
4025 27 50   27 1 93 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
4026 27 50       51 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
4027 27         37 my $x = shift;
4028 27         43 my $class = ref $x;
4029              
4030 27         37 my $y = shift;
4031 27 100 66     145 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
4032              
4033 27 100 100     72 $x -> is_scalar() || $y -> is_scalar() ? $x -> smul($y) : $x -> mmul($y);
4034             }
4035              
4036             =pod
4037              
4038             =item mmul()
4039              
4040             Matrix multiplication. An error is thrown if the sizes don't match; the number
4041             of columns in the first operand must be equal to the number of rows in the
4042             second operand.
4043              
4044             $z = $x -> mmul($y);
4045              
4046             =cut
4047              
4048             sub mmul {
4049 34 50   34 1 107 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
4050 34 50       72 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
4051 34         48 my $x = shift;
4052 34         55 my $class = ref $x;
4053              
4054 34         65 my $y = shift;
4055 34 50 33     195 $y = $class -> new($y) unless defined(blessed($y)) && $y -> isa($class);
4056              
4057 34         113 my $mx = $x -> nrow();
4058 34         76 my $nx = $x -> ncol();
4059              
4060 34         61 my $my = $y -> nrow();
4061 34         60 my $ny = $y -> ncol();
4062              
4063 34 50       70 croak "Can't multiply $mx-by-$nx matrix with $my-by-$ny matrix"
4064             unless $nx == $my;
4065              
4066 34         60 my $z = [];
4067 34         50 my $l = $nx - 1; # "inner length"
4068 34         88 for my $i (0 .. $mx - 1) {
4069 74         127 for my $j (0 .. $ny - 1) {
4070 192         534 $z -> [$i][$j] = _sum(map $x -> [$i][$_] * $y -> [$_][$j], 0 .. $l);
4071             }
4072             }
4073              
4074 34         183 bless $z, $class;
4075             }
4076              
4077             =pod
4078              
4079             =item smul()
4080              
4081             Scalar (element by element) multiplication with scalar expansion. This method
4082             places no requirements on the size of the input matrices.
4083              
4084             $z = $x -> smul($y);
4085              
4086             =cut
4087              
4088             sub smul {
4089 12 50   12 1 76 croak "Not enough arguments for ", (caller(0))[3] if @_ < 2;
4090 12 50       31 croak "Too many arguments for ", (caller(0))[3] if @_ > 2;
4091 12         20 my $x = shift;
4092              
4093 12     36   43 my $sub = sub { $_[0] * $_[1] };
  36         84  
4094 12         39 $x -> sapply($sub, @_);
4095             }
4096              
4097             =pod
4098              
4099             =item mmuladd()
4100              
4101             Matrix fused multiply and add. If C<$x> is a C<$p>-by-C<$q> matrix, then C<$y>
4102             must be a C<$q>-by-C<$r> matrix and C<$z> must be a C<$p>-by-C<$r> matrix. An
4103             error is thrown if the sizes don't match.
4104              
4105             $w = $x -> mmuladd($y, $z);
4106              
4107             The fused multiply and add is equivalent to, but computed with higher accuracy
4108             than
4109              
4110             $w = $x -> mmul($y) -> madd($z);
4111              
4112             This method can be used to improve the solution of linear systems.
4113              
4114             =cut
4115              
4116             sub mmuladd {