projects / OpenSPARC-T2-DV / blame
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Initial commit of OpenSPARC T2 design and verification files.
[OpenSPARC-T2-DV] / tools / perl-5.8.0 / lib / site_perl / 5.8.0 / Math / MatrixReal.pm
CommitLineData
86530b38
AT		 1# Copyright (c) 1996, 1997 by Steffen Beyer. All rights reserved.
2# Copyright (c) 1999 by Rodolphe Ortalo. All rights reserved.
3# Copyright (c) 2001,2002 by Jonathan Leto. All rights reserved.
4# This package is free software; you can redistribute it and/or
5# modify it under the same terms as Perl itself.
6
7package Math::MatrixReal;
8
9use strict;
10use Carp;
11use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);
12require Exporter;
13
14@ISA = qw(Exporter);
15@EXPORT = qw();
16@EXPORT_OK = qw(min max);
17%EXPORT_TAGS = (all => [@EXPORT_OK]);
18$VERSION = '1.9';
19
20use overload
21 'neg' => '_negate',
22 '~' => '_transpose',
23 'bool' => '_boolean',
24 '!' => '_not_boolean',
25 '""' => '_stringify',
26 'abs' => '_norm',
27 '+' => '_add',
28 '-' => '_subtract',
29 '*' => '_multiply',
30 '**' => '_exponent',
31 '+=' => '_assign_add',
32 '-=' => '_assign_subtract',
33 '*=' => '_assign_multiply',
34 '**=' => '_assign_exponent',
35 '==' => '_equal',
36 '!=' => '_not_equal',
37 '<' => '_less_than',
38 '<=' => '_less_than_or_equal',
39 '>' => '_greater_than',
40 '>=' => '_greater_than_or_equal',
41 'eq' => '_equal',
42 'ne' => '_not_equal',
43 'lt' => '_less_than',
44 'le' => '_less_than_or_equal',
45 'gt' => '_greater_than',
46 'ge' => '_greater_than_or_equal',
47 '=' => '_clone',
48 'exp' => '_exp',
49 'sin' => '_sin',
50 'cos' => '_cos',
51'fallback' => undef;
52
53sub new
54{
55 croak "Usage: \$new_matrix = Math::MatrixReal->new(\$rows,\$columns);"
56 if (@_ != 3);
57
58 my ($proto,$rows,$cols) = @_;
59 my $class = ref($proto) || $proto || 'Math::MatrixReal';
60 my($i,$j,$this);
61
62 croak "Math::MatrixReal::new(): number of rows must be integer > 0"
63 unless ($rows > 0 and $rows == int($rows) );
64
65 croak "Math::MatrixReal::new(): number of columns must be integer > 0"
66 unless ($cols > 0 and $cols == int($cols) );
67
68 $this = [ [ ], $rows, $cols ];
69
70 # Creates first empty row
71 my $empty = [ ];
72 $#$empty = $cols - 1; # Lengthens the array
73 for (my $j = 0; $j < $cols; $j++)
74 {
75 $empty->[$j] = 0.0;
76 }
77 $this->[0][0] = $empty;
78 # Creates other rows (by copying)
79 for (my $i = 1; $i < $rows; $i++)
80 {
81 my $arow = [ ];
82 @$arow = @$empty;
83 $this->[0][$i] = $arow;
84 }
85 bless($this, $class);
86 return($this);
87}
88sub new_diag {
89 croak "Usage: \$new_matrix = Math::MatrixReal->new_diag( [ 1, 2, 3] );"
90 unless (@_ == 2 );
91 my ($proto,$diag) = @_;
92 my $class = ref($proto) || $proto || 'Math::MatrixReal';
93 my $matrix;
94 my $n = scalar(@$diag);
95
96 croak "Math::MatrixReal::new_diag(): Third argument must be an arrayref"
97 unless (ref($diag) eq "ARRAY");
98
99 $matrix = Math::MatrixReal->new($n,$n);
100 $matrix = $matrix->each_diag( sub { shift @$diag } );
101 return $matrix;
102}
103
104sub new_from_string
105{
106 croak "Usage: \$new_matrix = Math::MatrixReal->new_from_string(\$string);"
107 if (@_ != 2);
108
109 my ($proto,$string) = @_;
110 my $class = ref($proto) || $proto || 'Math::MatrixReal';
111 my($line,$values);
112 my($rows,$cols);
113 my($row,$col);
114 my($warn,$this);
115
116 $warn = $rows = $cols = 0;
117
118 $values = [ ];
119 while ($string =~ m!^\s*
120 \[ \s+ ( (?: [+-]? \d+ (?: \. \d* )? (?: E [+-]? \d+ )? \s+ )+ ) \] \s*? \n
121 !ix)
122 {
123 $line = $1;
124 $string = $';
125 $values->[$rows] = [ ];
126 @{$values->[$rows]} = split(' ', $line);
127 $col = @{$values->[$rows]};
128 if ($col != $cols)
129 {
130 unless ($cols == 0) { $warn = 1; }
131 if ($col > $cols) { $cols = $col; }
132 }
133 $rows++;
134 }
135 if ($string !~ m!^\s*$!)
136 {
137 #croak
138 print "Math::MatrixReal::new_from_string(): syntax error in input string";
139 print "String is\n$string\n---\n";
140 croak;
141 }
142 if ($rows == 0)
143 {
144 croak "Math::MatrixReal::new_from_string(): empty input string";
145 }
146 if ($warn)
147 {
148 warn "Math::MatrixReal::new_from_string(): missing elements will be set to zero!\n";
149 }
150 $this = Math::MatrixReal::new($class,$rows,$cols);
151 for ( $row = 0; $row < $rows; $row++ )
152 {
153 for ( $col = 0; $col < @{$values->[$row]}; $col++ )
154 {
155 $this->[0][$row][$col] = $values->[$row][$col];
156 }
157 }
158 return($this);
159}
160# from Math::MatrixReal::Ext1
161sub new_from_cols {
162 my $proto = shift;
163 my $class = ref($proto) || $proto;
164 my $ref_to_cols = shift;
165 my @cols = @{$ref_to_cols};
166
167 my $cols = scalar( @cols );
168
169 my $matrix = 0;
170 my $rows = 0;
171
172 # each arg is a column, but we don't know what form they're
173 # in yet
174
175 my $col_index = 0;
176 foreach my $col (@cols) {
177 # it's one-based
178 $col_index ++;
179 my $ref = ref( $col ) ;
180
181 if ( $ref =~ /^Math::MatrixReal/ ) {
182 # it's already a Math::MatrixReal something
183 } elsif ( $ref eq 'ARRAY' ) {
184 my @array = @$col;
185 my $length = scalar( @array );
186 $col = $class->new_from_string( '[ '. join( " ]\n[ ", @array) ." ]\n" );
187 } elsif ( $ref eq '' ) {
188 # we hope this is a string
189 $col = $class->new_from_string( $col );
190 } else {
191 # we have no idea, error time!
192 croak __PACKAGE__."::new_from_cols(): sorry, I have no clue what you sent me! I only know how to deal with array refs, strings, and things that are already in the Math::MatrixReal hierarchy \n";
193 }
194 my ($length, $one) = $col->dim;
195 croak __PACKAGE__."::new_from_cols(): This isn't a column vector"
196 unless ($one == 1) ;
197
198 # if we already have a height, check that this is the same
199 if ($rows) {
200 croak __PACKAGE__."::new_from_cols(): This column vector has $length elements and an earlier one had $rows"
201 unless ($length == $rows) ;
202 }
203 # else, we have a new height
204 # TODO: maybe this should check for zero
205 else {
206 $rows = $length;
207 }
208
209 # create the matrix the first time through
210 unless ($matrix) {
211 $matrix = $class->new($rows, $cols);
212 }
213
214 foreach my $row_index (1..$rows){
215 my $value = $col->element($row_index, 1);
216 $matrix->assign($row_index, $col_index, $value);
217 }
218
219 }
220 return $matrix;
221}
222#from Math::MatrixReal::Ext1
223sub new_from_rows {
224 my $proto = shift;
225 my $class = ref($proto) || $proto;
226 my $ref_to_rows = shift;
227 my @rows = @{$ref_to_rows};
228
229 my $rows = scalar( @rows );
230
231 my $matrix = 0;
232 my $cols = 0;
233
234 # each arg is a column, but we don't know what form they're
235 # in yet
236
237 my $row_index = 0;
238 foreach my $row (@rows) {
239 # it's one-based
240 $row_index ++;
241 my $ref = ref( $row ) ;
242
243 if ( $ref =~ /^Math::MatrixReal/ ) {
244 # it's already a Math::MatrixReal something
245 }
246 elsif ( $ref eq 'ARRAY' ) {
247 my @array = @$row;
248 my $length = scalar( @array );
249 $row = $class->new_from_string( '[ '. join( " ", @array) ." ]\n" );
250 }
251 elsif ( $ref eq '' ) {
252 # we hope this is a string
253 $row = $class->new_from_string( $row );
254 }
255 else {
256 # we have no idea, error time!
257 croak "Math::MatrixReal::new_from_rows(): Argument must be an arrayref,string or Math::MatrixReal object";
258 }
259 my ($one, $length) = $row->dim;
260 croak __PACKAGE__."::new_from_rows(): This isn't a column vector"
261 unless ($one == 1) ;
262
263 # if we already have a height, check that this is the same
264 if ($cols) {
265 croak "Math::MatrixReal::new_from_rows(): Column mismatch $length != $cols"
266 unless ($length == $cols) ;
267 }
268 # else, we have a new width FIXME maybe this should check for zero
269 else {
270 $cols = $length;
271 }
272
273 # create the matrix the first time through
274 unless ($matrix) {
275 $matrix = $class->new($rows, $cols);
276 }
277
278 foreach my $col_index (1..$cols){
279 my $value = $row->element(1, $col_index);
280 $matrix->assign($row_index, $col_index, $value);
281 }
282
283 }
284 return $matrix;
285}
286
287sub shadow
288{
289 croak "Usage: \$new_matrix = \$some_matrix->shadow();"
290 if (@_ != 1);
291
292 my($matrix) = @_;
293 my($temp);
294
295 $temp = $matrix->new($matrix->[1],$matrix->[2]);
296 return($temp);
297}
298
299
300sub copy
301{
302 croak "Usage: \$matrix1->copy(\$matrix2);"
303 if (@_ != 2);
304
305 my($matrix1,$matrix2) = @_;
306 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
307 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
308 my($i,$j);
309
310 croak "Math::MatrixReal::copy(): matrix size mismatch"
311 unless (($rows1 == $rows2) && ($cols1 == $cols2));
312
313 for ( $i = 0; $i < $rows1; $i++ )
314 {
315 my $r1 = []; # New array ref
316 my $r2 = $matrix2->[0][$i];
317 @$r1 = @$r2; # Copy whole array directly
318 $matrix1->[0][$i] = $r1;
319 }
320 if (defined $matrix2->[3]) # is an LR decomposition matrix!
321 {
322 $matrix1->[3] = $matrix2->[3]; # $sign
323 $matrix1->[4] = $matrix2->[4]; # $perm_row
324 $matrix1->[5] = $matrix2->[5]; # $perm_col
325 }
326}
327
328sub clone
329{
330 croak "Usage: \$twin_matrix = \$some_matrix->clone();"
331 if (@_ != 1);
332
333 my($matrix) = @_;
334 my($temp);
335
336 $temp = $matrix->new($matrix->[1],$matrix->[2]);
337 $temp->copy($matrix);
338 return($temp);
339}
340
341## trace() : return the sum of the diagonal elements
342sub trace {
343 croak "Usage: \$trace = \$matrix->trace();" if (@_ != 1);
344
345 my $matrix = shift;
346 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
347 my $trace = 0;
348 my($j);
349
350 croak "Math::MatrixReal::trace(): matrix is not quadratic"
351 unless ($rows == $cols);
352
353
354 for ( $j = 0; $j < $cols; $j++ )
355 {
356 $trace += $matrix->[0][$j][$j];
357 }
358 return($trace);
359}
360
361## return the minor corresponding to $r and $c
362## eliminate row $r and col $c , and return the $r-1 by $c-1 matrix
363sub minor {
364 croak "Usage: \$minor = \$matrix->minor(\$r,\$c);" unless (@_ == 3);
365 my ($matrix,$r,$c) = @_;
366 my ($rows,$cols) = $matrix->dim();
367
368 croak "Math::MatrixReal::minor(): \$matrix must be at least 2x2"
369 unless ($rows > 1 and $cols > 1);
370 croak "Math::MatrixReal::minor(): $r and $c must be positive"
371 unless ($r > 0 and $c > 0 );
372 croak "Math::MatrixReal::minor(): matrix has no $r,$c element"
373 unless ($r <= $rows and $c <= $cols );
374
375 my $minor = new Math::MatrixReal($rows-1,$cols-1);
376 my ($i,$j) = (0,0);
377
378 ## assign() might have been easier, but this should be faster
379 ## the element can be in any of 4 regions compared to the eliminated
380 ## row and col:
381 ## above and to the left, above and to the right
382 ## below and to the left, below and to the right
383
384 for(; $i < $rows; $i++){
385 for(;$j < $rows; $j++ ){
386 if( $i >= $r && $j >= $c ){
387 $minor->[0][$i-1][$j-1] = $matrix->[0][$i][$j];
388 } elsif ( $i >= $r && $j < $c ){
389 $minor->[0][$i-1][$j] = $matrix->[0][$i][$j];
390 } elsif ( $i < $r && $j < $c ){
391 $minor->[0][$i][$j] = $matrix->[0][$i][$j];
392 } elsif ( $i < $r && $j >= $c ){
393 $minor->[0][$i][$j-1] = $matrix->[0][$i][$j];
394 } else {
395 croak "Very bad things";
396 }
397 }
398 $j = 0;
399 }
400 return ($minor);
401}
402sub swap_col {
403 croak "Usage: \$matrix->swap_col(\$col1,\$col2); " unless (@_ == 3);
404 my ($matrix,$col1,$col2) = @_;
405 my ($rows,$cols) = $matrix->dim();
406 my (@temp);
407
408 croak "Math::MatrixReal::swap_col(): col index is not valid"
409 unless ( $col1 <= $cols && $col2 <= $cols &&
410 $col1 == int($col1) &&
411 $col2 == int($col2) );
412 $col1--;$col2--;
413 for(my $i=0;$i < $rows;$i++){
414 $temp[$i] = $matrix->[0][$i][$col1];
415 $matrix->[0][$i][$col1] = $matrix->[0][$i][$col2];
416 $matrix->[0][$i][$col2] = $temp[$i];
417 }
418}
419sub swap_row {
420 croak "Usage: \$matrix->swap_row(\$row1,\$row2); " unless (@_ == 3);
421 my ($matrix,$row1,$row2) = @_;
422 my ($rows,$cols) = $matrix->dim();
423 my (@temp);
424
425 croak "Math::MatrixReal::swap_row(): row index is not valid"
426 unless ( $row1 <= $rows && $row2 <= $rows &&
427 $row1 == int($row1) &&
428 $row2 == int($row2) );
429 $row1--;$row2--;
430 for(my $j=0;$j < $cols;$j++){
431 $temp[$j] = $matrix->[0][$row1][$j];
432 $matrix->[0][$row1][$j] = $matrix->[0][$row2][$j];
433 $matrix->[0][$row2][$j] = $temp[$j];
434 }
435}
436# TODO: docs
437# no worky
438sub assign_row {
439 croak "Usage: \$matrix->assign_row(\$row,\$row_vec);" unless (@_ == 3);
440 my ($matrix,$row,$row_vec) = @_;
441 my ($rows1,$cols1) = $matrix->dim();
442 my ($rows2,$cols2) = $row_vec->dim();
443
444 croak "Math::MatrixReal::assign_row(): row mismatch" if ($rows1 != $rows2);
445 croak "Math::MatrixReal::assign_row(): not a row vector" unless( $cols2 == 1);
446
447 @{$matrix->[0][--$row]} = @{$row_vec->[0][0]};
448
449}
450# returns the number of zeroes in a row
451sub _count_zeroes_row {
452 my ($matrix) = @_;
453 my ($rows,$cols) = $matrix->dim();
454 my $count = 0;
455 croak "_count_zeroes_row(): only 1 row, buddy" unless ($rows == 1);
456
457 for(my $i=0;$i < $cols;$i++){
458 $count++ unless $matrix->[0][0][$i];
459 }
460 return $count;
461}
462## divide a row by it's largest abs() element
463sub _normalize_row {
464 my ($matrix) = @_;
465 my ($rows,$cols) = $matrix->dim();
466 my $new_row = Math::MatrixReal->new(1,$cols);
467
468 my $big = abs($matrix->[0][0][0]);
469 for(my $j=0;$j < $cols; $j++ ){
470 $big = $big < abs($matrix->[0][0][$j])
471 ? abs($matrix->[0][0][$j]) : $big;
472 }
473 next unless $big;
474 # now $big is biggest element in row
475 for(my $j = 0;$j < $cols; $j++ ){
476 $new_row->[0][0][$j] = $matrix->[0][0][$j] / $big;
477 }
478 return $new_row;
479
480}
481#### doesn't work
482=crap
483sub row_echelon {
484 my ($matrix) = @_;
485 my ($rows,$cols) = $matrix->dim();
486 my $big;
487 my $tmprow;
488
489 for(my $i = 0;$i < $rows; $i++ ){
490 $big = abs($matrix->[0][$i][0]);
491 for(my $j=0;$j < $cols; $j++ ){
492 $big = $big < abs($matrix->[0][$i][$j])
493 ? abs($matrix->[0][$i][$j]) : $big;
494 }
495 next unless $big;
496 # now $big is biggest element in row
497 for(my $j = 0;$j < $cols; $j++ ){
498 $matrix->[0][$i][$j] /= $big;
499 }
500 # now all elements are between [-1,1] and
501 # dependence can be easily checked
502 }
503
504 for(my $i = 0;$i < $rows; $i++ ){
505 for(my $k = 0; $k < $rows; $k++ ){
506 next if ($i == $k);
507 print "i,k = $i,$k\n";
508 $tmprow = $matrix->row($i+1) + $matrix->row($k+1);
509 print "tmprow= $tmprow\n";
510 print "tmprow has " . $tmprow->_count_zeroes_row . " zeroes.\n";
511 print "i+1 has " . $matrix->row($i+1)->_count_zeroes_row . " zeroes\n";
512 print "k+1 has " . $matrix->row($k+1)->_count_zeroes_row . " zeroes\n";
513
514 $tmprow = $tmprow->_normalize_row;
515 if( $tmprow->norm_sum == 0 ){
516 for(my $j=0;$j<$cols;$j++){
517 $matrix->[0][$k][$j] = 0;
518 }
519 } elsif ( $tmprow->_count_zeroes_row > $matrix->row($i+1)->_count_zeroes_row ){
520 print "tmprow has more zeroes, replacing row " . ($i+1) . "\n";
521 ## assign row
522 $matrix->assign_row($i+1,$tmprow);
523 } elsif ( $tmprow->_count_zeroes_row > $matrix->row($k+1)->_count_zeroes_row ){
524 print "tmprow has more zeroes, replacing row " . ($k+1) . "\n";
525 $matrix->assign_row($k+1,$tmprow);
526 }
527
528 }
529 }
530 return $matrix;
531}
532=cut
533
534sub cofactor {
535 my ($matrix) = @_;
536 my ($rows,$cols) = $matrix->dim();
537
538 croak "Math::MatrixReal::cofactor(): Matrix is not quadratic"
539 unless ($rows == $cols);
540
541 # black magic ahead
542 my $cofactor = $matrix->each( sub { my($v,$i,$j) = @_;
543 ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
544 : -1*$matrix->minor($i,$j)->det();
545 } );
546 return ($cofactor);
547}
548
549sub adjoint {
550 my ($matrix) = @_;
551 return ~($matrix->cofactor);
552}
553sub row
554{
555 croak "Usage: \$row_vector = \$matrix->row(\$row);"
556 if (@_ != 2);
557
558 my($matrix,$row) = @_;
559 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
560 my($temp);
561 my($j);
562
563 croak "Math::MatrixReal::row(): row index out of range"
564 if (($row < 1) || ($row > $rows));
565
566 $row--;
567 $temp = $matrix->new(1,$cols);
568 for ( $j = 0; $j < $cols; $j++ )
569 {
570 $temp->[0][0][$j] = $matrix->[0][$row][$j];
571 }
572 return($temp);
573}
574
575sub column
576{
577 croak "Usage: \$column_vector = \$matrix->column(\$column);"
578 if (@_ != 2);
579
580 my($matrix,$col) = @_;
581 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
582 my($temp);
583 my($i);
584
585 croak "Math::MatrixReal::column(): column index out of range"
586 if (($col < 1) || ($col > $cols));
587
588 $col--;
589 $temp = $matrix->new($rows,1);
590 for ( $i = 0; $i < $rows; $i++ )
591 {
592 $temp->[0][$i][0] = $matrix->[0][$i][$col];
593 }
594 return($temp);
595}
596
597sub _undo_LR
598{
599 croak "Usage: \$matrix->_undo_LR();"
600 if (@_ != 1);
601
602 my($this) = @_;
603
604 undef $this->[3];
605 undef $this->[4];
606 undef $this->[5];
607}
608
609sub zero
610{
611 croak "Usage: \$matrix->zero();"
612 if (@_ != 1);
613
614 my($this) = @_;
615 my($rows,$cols) = ($this->[1],$this->[2]);
616 my($i,$j);
617
618 $this->_undo_LR();
619
620 # Zero first row
621 for (my $j = 0; $j < $cols; $j++ )
622 {
623 $this->[0][0][$j] = 0.0;
624 }
625 # Then propagate to other rows
626 for (my $i = 0; $i < $rows; $i++)
627 {
628 @{$this->[0][$i]} = @{$this->[0][0]};
629 }
630}
631
632sub one
633{
634 croak "Usage: \$matrix->one();"
635 if (@_ != 1);
636
637 my($this) = @_;
638 my($rows,$cols) = ($this->[1],$this->[2]);
639 my($i,$j);
640
641# No need for this: done by the 'zero()'
642# $this->_undo_LR();
643 $this->zero(); # We rely on zero() efficiency
644 for (my $i = 0; $i < $rows; $i++ )
645 {
646 $this->[0][$i][$i] = 1.0;
647 }
648}
649
650sub assign
651{
652 croak "Usage: \$matrix->assign(\$row,\$column,\$value);"
653 if (@_ != 4);
654
655 my($this,$row,$col,$value) = @_;
656 my($rows,$cols) = ($this->[1],$this->[2]);
657
658 croak "Math::MatrixReal::assign(): row index out of range"
659 if (($row < 1) || ($row > $rows));
660
661 croak "Math::MatrixReal::assign(): column index out of range"
662 if (($col < 1) || ($col > $cols));
663
664 $this->_undo_LR();
665
666 $this->[0][--$row][--$col] = $value;
667}
668
669sub element
670{
671 croak "Usage: \$value = \$matrix->element(\$row,\$column);"
672 if (@_ != 3);
673
674 my($this,$row,$col) = @_;
675 my($rows,$cols) = ($this->[1],$this->[2]);
676
677 croak "Math::MatrixReal::element(): row index out of range"
678 if (($row < 1) || ($row > $rows));
679
680 croak "Math::MatrixReal::element(): column index out of range"
681 if (($col < 1) || ($col > $cols));
682
683 return( $this->[0][--$row][--$col] );
684}
685
686sub dim # returns dimensions of a matrix
687{
688 croak "Usage: (\$rows,\$columns) = \$matrix->dim();"
689 if (@_ != 1);
690
691 my($matrix) = @_;
692
693 return( $matrix->[1], $matrix->[2] );
694}
695
696sub norm_one # maximum of sums of each column
697{
698 croak "Usage: \$norm_one = \$matrix->norm_one();"
699 if (@_ != 1);
700
701 my($this) = @_;
702 my($rows,$cols) = ($this->[1],$this->[2]);
703
704 my $max = 0.0;
705 for (my $j = 0; $j < $cols; $j++)
706 {
707 my $sum = 0.0;
708 for (my $i = 0; $i < $rows; $i++)
709 {
710 $sum += abs( $this->[0][$i][$j] );
711 }
712 $max = $sum if ($sum > $max);
713 }
714 return($max);
715}
716## sum of absolute value of every element
717sub norm_sum {
718 croak "Usage: \$norm_sum = \$matrix->norm_sum();"
719 unless (@_ == 1);
720 my ($matrix) = @_;
721 my $norm = 0;
722 $matrix->each( sub { $norm+=abs(shift); } );
723 return $norm;
724}
725sub norm_frobenius {
726 my ($m) = @_;
727 my ($r,$c) = $m->dim;
728 my $s=0;
729
730 $m->each( sub { $s+=abs(shift)**2 } );
731 return sqrt($s);
732}
733
734# Vector Norm
735sub norm_p {
736 my ($v,$p) = @_;
737 # sanity check on $p
738 croak "Math::MatrixReal:norm_p: argument must be a row or column vector"
739 unless ( $v->is_row_vector || $v->is_col_vector );
740 croak "Math::MatrixReal::norm_p: $p must be >= 1"
741 unless ($p =~ m/Inf(inity)?/i || $p >= 1);
742
743 if( $p =~ m/^(Inf|Infinity)$/i ){
744 my $max = $v->element(1,1);
745 $v->each ( sub { my $x=abs(shift); $max = $x if( $x > $max ); } );
746 return $max;
747 }
748
749 my $s=0;
750 $v->each( sub { $s+= (abs(shift))**$p; } );
751 return $s ** (1/$p);
752
753}
754sub norm_max # maximum of sums of each row
755{
756 croak "Usage: \$norm_max = \$matrix->norm_max();"
757 if (@_ != 1);
758
759 my($this) = @_;
760 my($rows,$cols) = ($this->[1],$this->[2]);
761
762 my $max = 0.0;
763 for (my $i = 0; $i < $rows; $i++)
764 {
765 my $sum = 0.0;
766 for (my $j = 0; $j < $cols; $j++)
767 {
768 $sum += abs( $this->[0][$i][$j] );
769 }
770 $max = $sum if ($sum > $max);
771 }
772 return($max);
773}
774
775sub negate
776{
777 croak "Usage: \$matrix1->negate(\$matrix2);"
778 if (@_ != 2);
779
780 my($matrix1,$matrix2) = @_;
781 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
782 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
783
784 croak "Math::MatrixReal::negate(): matrix size mismatch"
785 unless (($rows1 == $rows2) && ($cols1 == $cols2));
786
787 $matrix1->_undo_LR();
788
789 for (my $i = 0; $i < $rows1; $i++ )
790 {
791 for (my $j = 0; $j < $cols1; $j++ )
792 {
793 $matrix1->[0][$i][$j] = -($matrix2->[0][$i][$j]);
794 }
795 }
796}
797## each(): evaluate a coderef on each element and return a new matrix
798## of said
799sub each {
800 croak "Usage: \$new_matrix = \$matrix->each( \&sub );" unless (@_ == 2 );
801 my($matrix,$function) = @_;
802 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
803 my($new_matrix) = $matrix->clone();
804
805 croak "Math::MatrixReal::each(): argument is not a sub reference" unless ref($function);
806 $new_matrix->_undo_LR();
807
808 for (my $i = 0; $i < $rows; $i++ ) {
809 for (my $j = 0; $j < $cols; $j++ ) {
810 no strict 'refs';
811 # $i,$j are 1-based as of 1.7
812 $new_matrix->[0][$i][$j] = &{ $function }($matrix->[0][$i][$j],$i+1,$j+1) ;
813 }
814 }
815 return ($new_matrix);
816}
817## each_diag(): same as each() but only diag elements
818sub each_diag {
819 croak "Usage: \$new_matrix = \$matrix->each_diag( \&sub );" unless (@_ == 2 );
820 my($matrix,$function) = @_;
821 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
822 my($new_matrix) = $matrix->clone();
823
824 croak "Math::MatrixReal::each(): argument is not a sub reference" unless ref($function);
825 croak "Matrix is not quadratic" unless ($rows == $cols);
826
827 $new_matrix->_undo_LR();
828
829 for (my $i = 0; $i < $rows; $i++ ) {
830 for (my $j = 0; $j < $cols; $j++ ) {
831 next unless ($i == $j);
832 no strict 'refs';
833 # $i,$j are 1-based as of 1.7
834 $new_matrix->[0][$i][$j] = &{ $function }($matrix->[0][$i][$j],$i+1,$j+1) ;
835 }
836 }
837 return ($new_matrix);
838}
839
840## Make computing the inverse more user friendly
841sub inverse {
842 croak "Usage: \$inverse = \$matrix->inverse();" unless (@_ == 1);
843 my ($matrix) = @_;
844 return $matrix->decompose_LR->invert_LR;
845}
846
847sub transpose {
848
849 croak "Usage: \$matrix1->transpose(\$matrix2);"
850 if (@_ != 2);
851
852 my($matrix1,$matrix2) = @_;
853 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
854 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
855
856 croak "Math::MatrixReal::transpose(): matrix size mismatch"
857 unless (($rows1 == $cols2) && ($cols1 == $rows2));
858
859 $matrix1->_undo_LR();
860
861 if ($rows1 == $cols1)
862 {
863 # more complicated to make in-place possible!
864
865 for (my $i = 0; $i < $rows1; $i++)
866 {
867 for (my $j = ($i + 1); $j < $cols1; $j++)
868 {
869 my $swap = $matrix2->[0][$i][$j];
870 $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
871 $matrix1->[0][$j][$i] = $swap;
872 }
873 $matrix1->[0][$i][$i] = $matrix2->[0][$i][$i];
874 }
875 }
876 else # ($rows1 != $cols1)
877 {
878 for (my $i = 0; $i < $rows1; $i++)
879 {
880 for (my $j = 0; $j < $cols1; $j++)
881 {
882 $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
883 }
884 }
885 }
886}
887
888sub add
889{
890 croak "Usage: \$matrix1->add(\$matrix2,\$matrix3);"
891 if (@_ != 3);
892
893 my($matrix1,$matrix2,$matrix3) = @_;
894 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
895 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
896 my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
897 my($i,$j);
898
899 croak "Math::MatrixReal::add(): matrix size mismatch"
900 unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
901 ($cols1 == $cols2) && ($cols1 == $cols3));
902
903 $matrix1->_undo_LR();
904
905 for ( $i = 0; $i < $rows1; $i++ )
906 {
907 for ( $j = 0; $j < $cols1; $j++ )
908 {
909 $matrix1->[0][$i][$j] =
910 $matrix2->[0][$i][$j] + $matrix3->[0][$i][$j];
911 }
912 }
913}
914
915sub subtract
916{
917 croak "Usage: \$matrix1->subtract(\$matrix2,\$matrix3);"
918 if (@_ != 3);
919
920 my($matrix1,$matrix2,$matrix3) = @_;
921 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
922 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
923 my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
924 my($i,$j);
925
926 croak "Math::MatrixReal::subtract(): matrix size mismatch"
927 unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
928 ($cols1 == $cols2) && ($cols1 == $cols3));
929
930 $matrix1->_undo_LR();
931
932 for ( $i = 0; $i < $rows1; $i++ )
933 {
934 for ( $j = 0; $j < $cols1; $j++ )
935 {
936 $matrix1->[0][$i][$j] =
937 $matrix2->[0][$i][$j] - $matrix3->[0][$i][$j];
938 }
939 }
940}
941
942sub multiply_scalar
943{
944 croak "Usage: \$matrix1->multiply_scalar(\$matrix2,\$scalar);"
945 if (@_ != 3);
946
947 my($matrix1,$matrix2,$scalar) = @_;
948 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
949 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
950 my($i,$j);
951
952 croak "Math::MatrixReal::multiply_scalar(): matrix size mismatch"
953 unless (($rows1 == $rows2) && ($cols1 == $cols2));
954
955 $matrix1->_undo_LR();
956
957 for ( $i = 0; $i < $rows1; $i++ )
958 {
959 for ( $j = 0; $j < $cols1; $j++ )
960 {
961 $matrix1->[0][$i][$j] = $matrix2->[0][$i][$j] * $scalar;
962 }
963 }
964}
965
966sub multiply
967{
968 croak "Usage: \$product_matrix = \$matrix1->multiply(\$matrix2);"
969 if (@_ != 2);
970
971 my($matrix1,$matrix2) = @_;
972 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
973 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
974 my($temp);
975
976 croak "Math::MatrixReal::multiply(): matrix size mismatch"
977 unless ($cols1 == $rows2);
978
979 $temp = $matrix1->new($rows1,$cols2);
980 for (my $i = 0; $i < $rows1; $i++ )
981 {
982 for (my $j = 0; $j < $cols2; $j++ )
983 {
984 my $sum = 0.0;
985 for (my $k = 0; $k < $cols1; $k++ )
986 {
987 $sum += ( $matrix1->[0][$i][$k] * $matrix2->[0][$k][$j] );
988 }
989 $temp->[0][$i][$j] = $sum;
990 }
991 }
992 return($temp);
993}
994
995sub exponent {
996 croak "Usage: \$matrix_exp = \$matrix1->exponent(\$integer);" if(@_ != 2 );
997 my($matrix,$argument) = @_;
998 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
999 my($name) = "'**'";
1000 my($temp) = $matrix->clone();
1001
1002 croak "Matrix is not quadratic" unless ($rows == $cols);
1003 croak "Exponent must be integer" unless ($argument =~ m/^[+-]?\d+$/ );
1004
1005 return($matrix) if ($argument == 1);
1006
1007 $temp->_undo_LR();
1008
1009 # negative exponent is (A^-1)^n
1010 if( $argument < 0 ){
1011 my $LR = $matrix->decompose_LR();
1012 my $inverse = $LR->invert_LR();
1013 unless (defined $inverse){
1014 carp "Matrix has no inverse";
1015 return undef;
1016 }
1017 $temp = $inverse->clone();
1018 if( $inverse ){
1019 return($inverse) if ($argument == -1);
1020 for( 2 .. abs($argument) ){
1021 $temp = multiply($inverse,$temp);
1022 }
1023 return($temp);
1024 } else {
1025 # TODO: is this the right behaviour?
1026 carp "Cannot compute negative exponent, inverse does not exist";
1027 return undef;
1028 }
1029 # matrix to zero power is identity matrix
1030 } elsif( $argument == 0 ){
1031 $temp->one();
1032 return ($temp);
1033 }
1034
1035 # if it is diagonal, just raise diagonal entries to power
1036 if( $matrix->is_diagonal() ){
1037 $temp = $temp->each_diag( sub { (shift)**$argument } );
1038 return ($temp);
1039
1040 } else {
1041 for( 2 .. $argument ){
1042 $temp = multiply($matrix,$temp);
1043 }
1044 return ($temp);
1045 }
1046}
1047
1048sub min
1049{
1050 croak "Usage: \$minimum = Math::MatrixReal::min(\$number1,\$number2);"
1051 if (@_ != 2);
1052
1053 return( $_[0] < $_[1] ? $_[0] : $_[1] );
1054}
1055
1056sub max
1057{
1058 croak "Usage: \$maximum = Math::MatrixReal::max(\$number1,\$number2);"
1059 if (@_ != 2);
1060
1061 return( $_[0] > $_[1] ? $_[0] : $_[1] );
1062}
1063
1064sub kleene
1065{
1066 croak "Usage: \$minimal_cost_matrix = \$cost_matrix->kleene();"
1067 if (@_ != 1);
1068
1069 my($matrix) = @_;
1070 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1071 my($i,$j,$k,$n);
1072 my($temp);
1073
1074 croak "Math::MatrixReal::kleene(): matrix is not quadratic"
1075 unless ($rows == $cols);
1076
1077 $temp = $matrix->new($rows,$cols);
1078 $temp->copy($matrix);
1079 $temp->_undo_LR();
1080 $n = $rows;
1081 for ( $i = 0; $i < $n; $i++ )
1082 {
1083 $temp->[0][$i][$i] = min( $temp->[0][$i][$i] , 0 );
1084 }
1085 for ( $k = 0; $k < $n; $k++ )
1086 {
1087 for ( $i = 0; $i < $n; $i++ )
1088 {
1089 for ( $j = 0; $j < $n; $j++ )
1090 {
1091 $temp->[0][$i][$j] = min( $temp->[0][$i][$j] ,
1092 ( $temp->[0][$i][$k] +
1093 $temp->[0][$k][$j] ) );
1094 }
1095 }
1096 }
1097 return($temp);
1098}
1099
1100sub normalize
1101{
1102 croak "Usage: (\$norm_matrix,\$norm_vector) = \$matrix->normalize(\$vector);"
1103 if (@_ != 2);
1104
1105 my($matrix,$vector) = @_;
1106 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1107 my($norm_matrix,$norm_vector);
1108 my($max,$val);
1109 my($i,$j,$n);
1110
1111 croak "Math::MatrixReal::normalize(): matrix is not quadratic"
1112 unless ($rows == $cols);
1113
1114 $n = $rows;
1115
1116 croak "Math::MatrixReal::normalize(): vector is not a column vector"
1117 unless ($vector->[2] == 1);
1118
1119 croak "Math::MatrixReal::normalize(): matrix and vector size mismatch"
1120 unless ($vector->[1] == $n);
1121
1122 $norm_matrix = $matrix->new($n,$n);
1123 $norm_vector = $vector->new($n,1);
1124
1125 $norm_matrix->copy($matrix);
1126 $norm_vector->copy($vector);
1127
1128 $norm_matrix->_undo_LR();
1129
1130 for ( $i = 0; $i < $n; $i++ )
1131 {
1132 $max = abs($norm_vector->[0][$i][0]);
1133 for ( $j = 0; $j < $n; $j++ )
1134 {
1135 $val = abs($norm_matrix->[0][$i][$j]);
1136 if ($val > $max) { $max = $val; }
1137 }
1138 if ($max != 0)
1139 {
1140 $norm_vector->[0][$i][0] /= $max;
1141 for ( $j = 0; $j < $n; $j++ )
1142 {
1143 $norm_matrix->[0][$i][$j] /= $max;
1144 }
1145 }
1146 }
1147 return($norm_matrix,$norm_vector);
1148}
1149
1150sub decompose_LR
1151{
1152 croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
1153 if (@_ != 1);
1154
1155 my($matrix) = @_;
1156 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1157 my($perm_row,$perm_col);
1158 my($row,$col,$max);
1159 my($i,$j,$k,$n);
1160 my($sign) = 1;
1161 my($swap);
1162 my($temp);
1163
1164 croak "Math::MatrixReal::decompose_LR(): matrix is not quadratic"
1165 unless ($rows == $cols);
1166
1167 $temp = $matrix->new($rows,$cols);
1168 $temp->copy($matrix);
1169 $n = $rows;
1170 $perm_row = [ ];
1171 $perm_col = [ ];
1172 for ( $i = 0; $i < $n; $i++ )
1173 {
1174 $perm_row->[$i] = $i;
1175 $perm_col->[$i] = $i;
1176 }
1177 NONZERO:
1178 for ( $k = 0; $k < $n; $k++ ) # use Gauss's algorithm:
1179 {
1180 # complete pivot-search:
1181
1182 $max = 0;
1183 for ( $i = $k; $i < $n; $i++ )
1184 {
1185 for ( $j = $k; $j < $n; $j++ )
1186 {
1187 if (($swap = abs($temp->[0][$i][$j])) > $max)
1188 {
1189 $max = $swap;
1190 $row = $i;
1191 $col = $j;
1192 }
1193 }
1194 }
1195 last NONZERO if ($max == 0); # (all remaining elements are zero)
1196 if ($k != $row) # swap row $k and row $row:
1197 {
1198 $sign = -$sign;
1199 $swap = $perm_row->[$k];
1200 $perm_row->[$k] = $perm_row->[$row];
1201 $perm_row->[$row] = $swap;
1202 for ( $j = 0; $j < $n; $j++ )
1203 {
1204 # (must run from 0 since L has to be swapped too!)
1205
1206 $swap = $temp->[0][$k][$j];
1207 $temp->[0][$k][$j] = $temp->[0][$row][$j];
1208 $temp->[0][$row][$j] = $swap;
1209 }
1210 }
1211 if ($k != $col) # swap column $k and column $col:
1212 {
1213 $sign = -$sign;
1214 $swap = $perm_col->[$k];
1215 $perm_col->[$k] = $perm_col->[$col];
1216 $perm_col->[$col] = $swap;
1217 for ( $i = 0; $i < $n; $i++ )
1218 {
1219 $swap = $temp->[0][$i][$k];
1220 $temp->[0][$i][$k] = $temp->[0][$i][$col];
1221 $temp->[0][$i][$col] = $swap;
1222 }
1223 }
1224 for ( $i = ($k + 1); $i < $n; $i++ )
1225 {
1226 # scan the remaining rows, add multiples of row $k to row $i:
1227
1228 $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
1229 if ($swap != 0)
1230 {
1231 # calculate a row of matrix R:
1232
1233 for ( $j = ($k + 1); $j < $n; $j++ )
1234 {
1235 $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
1236 }
1237
1238 # store matrix L in same matrix as R:
1239
1240 $temp->[0][$i][$k] = $swap;
1241 }
1242 }
1243 }
1244 $temp->[3] = $sign;
1245 $temp->[4] = $perm_row;
1246 $temp->[5] = $perm_col;
1247 return($temp);
1248}
1249
1250sub solve_LR
1251{
1252 croak "Usage: (\$dimension,\$x_vector,\$base_matrix) = \$LR_matrix->solve_LR(\$b_vector);"
1253 if (@_ != 2);
1254
1255 my($LR_matrix,$b_vector) = @_;
1256 my($rows,$cols) = ($LR_matrix->[1],$LR_matrix->[2]);
1257 my($dimension,$x_vector,$base_matrix);
1258 my($perm_row,$perm_col);
1259 my($y_vector,$sum);
1260 my($i,$j,$k,$n);
1261
1262 croak "Math::MatrixReal::solve_LR(): not an LR decomposition matrix"
1263 unless ((defined $LR_matrix->[3]) && ($rows == $cols));
1264
1265 $n = $rows;
1266
1267 croak "Math::MatrixReal::solve_LR(): vector is not a column vector"
1268 unless ($b_vector->[2] == 1);
1269
1270 croak "Math::MatrixReal::solve_LR(): matrix and vector size mismatch"
1271 unless ($b_vector->[1] == $n);
1272
1273 $perm_row = $LR_matrix->[4];
1274 $perm_col = $LR_matrix->[5];
1275
1276 $x_vector = $b_vector->new($n,1);
1277 $y_vector = $b_vector->new($n,1);
1278 $base_matrix = $LR_matrix->new($n,$n);
1279
1280 # calculate "x" so that LRx = b ==> calculate Ly = b, Rx = y:
1281
1282 for ( $i = 0; $i < $n; $i++ ) # calculate $y_vector:
1283 {
1284 $sum = $b_vector->[0][($perm_row->[$i])][0];
1285 for ( $j = 0; $j < $i; $j++ )
1286 {
1287 $sum -= $LR_matrix->[0][$i][$j] * $y_vector->[0][$j][0];
1288 }
1289 $y_vector->[0][$i][0] = $sum;
1290 }
1291
1292 $dimension = 0;
1293 for ( $i = ($n - 1); $i >= 0; $i-- ) # calculate $x_vector:
1294 {
1295 if ($LR_matrix->[0][$i][$i] == 0)
1296 {
1297 if ($y_vector->[0][$i][0] != 0)
1298 {
1299 return(); # a solution does not exist!
1300 }
1301 else
1302 {
1303 $dimension++;
1304 $x_vector->[0][($perm_col->[$i])][0] = 0;
1305 }
1306 }
1307 else
1308 {
1309 $sum = $y_vector->[0][$i][0];
1310 for ( $j = ($i + 1); $j < $n; $j++ )
1311 {
1312 $sum -= $LR_matrix->[0][$i][$j] *
1313 $x_vector->[0][($perm_col->[$j])][0];
1314 }
1315 $x_vector->[0][($perm_col->[$i])][0] =
1316 $sum / $LR_matrix->[0][$i][$i];
1317 }
1318 }
1319 if ($dimension)
1320 {
1321 if ($dimension == $n)
1322 {
1323 $base_matrix->one();
1324 }
1325 else
1326 {
1327 for ( $k = 0; $k < $dimension; $k++ )
1328 {
1329 $base_matrix->[0][($perm_col->[($n-$k-1)])][$k] = 1;
1330 for ( $i = ($n-$dimension-1); $i >= 0; $i-- )
1331 {
1332 $sum = 0;
1333 for ( $j = ($i + 1); $j < $n; $j++ )
1334 {
1335 $sum -= $LR_matrix->[0][$i][$j] *
1336 $base_matrix->[0][($perm_col->[$j])][$k];
1337 }
1338 $base_matrix->[0][($perm_col->[$i])][$k] =
1339 $sum / $LR_matrix->[0][$i][$i];
1340 }
1341 }
1342 }
1343 }
1344 return( $dimension, $x_vector, $base_matrix );
1345}
1346
1347sub invert_LR
1348{
1349 croak "Usage: \$inverse_matrix = \$LR_matrix->invert_LR();"
1350 if (@_ != 1);
1351
1352 my($matrix) = @_;
1353 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1354 my($inv_matrix,$x_vector,$y_vector);
1355 my($i,$j,$n);
1356
1357 croak "Math::MatrixReal::invert_LR(): not an LR decomposition matrix"
1358 unless ((defined $matrix->[3]) && ($rows == $cols));
1359
1360 $n = $rows;
1361 if ($matrix->[0][$n-1][$n-1] != 0)
1362 {
1363 $inv_matrix = $matrix->new($n,$n);
1364 $y_vector = $matrix->new($n,1);
1365 for ( $j = 0; $j < $n; $j++ )
1366 {
1367 if ($j > 0)
1368 {
1369 $y_vector->[0][$j-1][0] = 0;
1370 }
1371 $y_vector->[0][$j][0] = 1;
1372 if (($rows,$x_vector,$cols) = $matrix->solve_LR($y_vector))
1373 {
1374 for ( $i = 0; $i < $n; $i++ )
1375 {
1376 $inv_matrix->[0][$i][$j] = $x_vector->[0][$i][0];
1377 }
1378 }
1379 else
1380 {
1381 die "Math::MatrixReal::invert_LR(): unexpected error - please inform author!\n";
1382 }
1383 }
1384 return($inv_matrix);
1385 }
1386 else { return(); } # matrix is not invertible!
1387}
1388
1389sub condition
1390{
1391 # 1st matrix MUST be the inverse of 2nd matrix (or vice-versa)
1392 # for a meaningful result!
1393
1394 # make this work when given no args
1395
1396 croak "Usage: \$condition = \$matrix->condition(\$inverse_matrix);"
1397 if (@_ != 2);
1398
1399 my($matrix1,$matrix2) = @_;
1400 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
1401 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
1402
1403 croak "Math::MatrixReal::condition(): 1st matrix is not quadratic"
1404 unless ($rows1 == $cols1);
1405
1406 croak "Math::MatrixReal::condition(): 2nd matrix is not quadratic"
1407 unless ($rows2 == $cols2);
1408
1409 croak "Math::MatrixReal::condition(): matrix size mismatch"
1410 unless (($rows1 == $rows2) && ($cols1 == $cols2));
1411
1412 return( $matrix1->norm_one() * $matrix2->norm_one() );
1413}
1414
1415## easy to use determinant
1416## very fast if matrix is diagonal or triangular
1417
1418sub det {
1419 croak "Usage: \$determinant = \$matrix->det_LR();" unless (@_ == 1);
1420 my ($matrix) = @_;
1421 my ($rows,$cols) = $matrix->dim();
1422 my $det = 1;
1423
1424 croak "Math::MatrixReal::det(): Matrix is not quadratic"
1425 unless ($rows == $cols);
1426
1427 # diagonal will match too
1428 if( $matrix->is_upper_triangular() ){
1429 $matrix->each_diag( sub { $det*=shift; } );
1430 } elsif ( $matrix->is_lower_triangular() ){
1431 $matrix->each_diag( sub { $det*=shift; } );
1432 } else {
1433 return $matrix->decompose_LR->det_LR();
1434 }
1435 return $det;
1436}
1437
1438sub det_LR # determinant of LR decomposition matrix
1439{
1440 croak "Usage: \$determinant = \$LR_matrix->det_LR();"
1441 if (@_ != 1);
1442
1443 my($matrix) = @_;
1444 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1445 my($k,$det);
1446
1447 croak "Math::MatrixReal::det_LR(): not an LR decomposition matrix"
1448 unless ((defined $matrix->[3]) && ($rows == $cols));
1449
1450 $det = $matrix->[3];
1451 for ( $k = 0; $k < $rows; $k++ )
1452 {
1453 $det *= $matrix->[0][$k][$k];
1454 }
1455 return($det);
1456}
1457
1458sub rank_LR {
1459 return (shift)->order_LR;
1460}
1461
1462sub order_LR # order of LR decomposition matrix (number of non-zero equations)
1463{
1464 croak "Usage: \$order = \$LR_matrix->order_LR();"
1465 if (@_ != 1);
1466
1467 my($matrix) = @_;
1468 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1469 my($order);
1470
1471 croak "Math::MatrixReal::order_LR(): not an LR decomposition matrix"
1472 unless ((defined $matrix->[3]) && ($rows == $cols));
1473
1474 ZERO:
1475 for ( $order = ($rows - 1); $order >= 0; $order-- )
1476 {
1477 last ZERO if ($matrix->[0][$order][$order] != 0);
1478 }
1479 return(++$order);
1480}
1481
1482sub scalar_product
1483{
1484 croak "Usage: \$scalar_product = \$vector1->scalar_product(\$vector2);"
1485 if (@_ != 2);
1486
1487 my($vector1,$vector2) = @_;
1488 my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
1489 my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
1490 my($k,$sum);
1491
1492 croak "Math::MatrixReal::scalar_product(): 1st vector is not a column vector"
1493 unless ($cols1 == 1);
1494
1495 croak "Math::MatrixReal::scalar_product(): 2nd vector is not a column vector"
1496 unless ($cols2 == 1);
1497
1498 croak "Math::MatrixReal::scalar_product(): vector size mismatch"
1499 unless ($rows1 == $rows2);
1500
1501 $sum = 0;
1502 for ( $k = 0; $k < $rows1; $k++ )
1503 {
1504 $sum += $vector1->[0][$k][0] * $vector2->[0][$k][0];
1505 }
1506 return($sum);
1507}
1508
1509sub vector_product
1510{
1511 croak "Usage: \$vector_product = \$vector1->vector_product(\$vector2);"
1512 if (@_ != 2);
1513
1514 my($vector1,$vector2) = @_;
1515 my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
1516 my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
1517 my($temp);
1518 my($n);
1519
1520 croak "Math::MatrixReal::vector_product(): 1st vector is not a column vector"
1521 unless ($cols1 == 1);
1522
1523 croak "Math::MatrixReal::vector_product(): 2nd vector is not a column vector"
1524 unless ($cols2 == 1);
1525
1526 croak "Math::MatrixReal::vector_product(): vector size mismatch"
1527 unless ($rows1 == $rows2);
1528
1529 $n = $rows1;
1530
1531 croak "Math::MatrixReal::vector_product(): only defined for 3 dimensions"
1532 unless ($n == 3);
1533
1534 $temp = $vector1->new($n,1);
1535 $temp->[0][0][0] = $vector1->[0][1][0] * $vector2->[0][2][0] -
1536 $vector1->[0][2][0] * $vector2->[0][1][0];
1537 $temp->[0][1][0] = $vector1->[0][2][0] * $vector2->[0][0][0] -
1538 $vector1->[0][0][0] * $vector2->[0][2][0];
1539 $temp->[0][2][0] = $vector1->[0][0][0] * $vector2->[0][1][0] -
1540 $vector1->[0][1][0] * $vector2->[0][0][0];
1541 return($temp);
1542}
1543
1544sub length
1545{
1546 croak "Usage: \$length = \$vector->length();"
1547 if (@_ != 1);
1548
1549 my($vector) = @_;
1550 my($rows,$cols) = ($vector->[1],$vector->[2]);
1551 my($k,$comp,$sum);
1552
1553 croak "Math::MatrixReal::length(): vector is not a column vector"
1554 unless ($cols == 1);
1555
1556 $sum = 0;
1557 for ( $k = 0; $k < $rows; $k++ )
1558 {
1559 $comp = $vector->[0][$k][0];
1560 $sum += $comp * $comp;
1561 }
1562 return( sqrt( $sum ) );
1563}
1564
1565sub _init_iteration
1566{
1567 croak "Usage: \$which_norm = \$matrix->_init_iteration();"
1568 if (@_ != 1);
1569
1570 my($matrix) = @_;
1571 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1572 my($ok,$max,$sum,$norm);
1573 my($i,$j,$n);
1574
1575 croak "Math::MatrixReal::_init_iteration(): matrix is not quadratic"
1576 unless ($rows == $cols);
1577
1578 $ok = 1;
1579 $n = $rows;
1580 for ( $i = 0; $i < $n; $i++ )
1581 {
1582 if ($matrix->[0][$i][$i] == 0) { $ok = 0; }
1583 }
1584 if ($ok)
1585 {
1586 $norm = 1; # norm_one
1587 $max = 0;
1588 for ( $j = 0; $j < $n; $j++ )
1589 {
1590 $sum = 0;
1591 for ( $i = 0; $i < $j; $i++ )
1592 {
1593 $sum += abs($matrix->[0][$i][$j]);
1594 }
1595 for ( $i = ($j + 1); $i < $n; $i++ )
1596 {
1597 $sum += abs($matrix->[0][$i][$j]);
1598 }
1599 $sum /= abs($matrix->[0][$j][$j]);
1600 if ($sum > $max) { $max = $sum; }
1601 }
1602 $ok = ($max < 1);
1603 unless ($ok)
1604 {
1605 $norm = -1; # norm_max
1606 $max = 0;
1607 for ( $i = 0; $i < $n; $i++ )
1608 {
1609 $sum = 0;
1610 for ( $j = 0; $j < $i; $j++ )
1611 {
1612 $sum += abs($matrix->[0][$i][$j]);
1613 }
1614 for ( $j = ($i + 1); $j < $n; $j++ )
1615 {
1616 $sum += abs($matrix->[0][$i][$j]);
1617 }
1618 $sum /= abs($matrix->[0][$i][$i]);
1619 if ($sum > $max) { $max = $sum; }
1620 }
1621 $ok = ($max < 1)
1622 }
1623 }
1624 if ($ok) { return($norm); }
1625 else { return(0); }
1626}
1627
1628sub solve_GSM # Global Step Method
1629{
1630 croak "Usage: \$xn_vector = \$matrix->solve_GSM(\$x0_vector,\$b_vector,\$epsilon);"
1631 if (@_ != 4);
1632
1633 my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
1634 my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1635 my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1636 my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1637 my($norm,$sum,$diff);
1638 my($xn_vector);
1639 my($i,$j,$n);
1640
1641 croak "Math::MatrixReal::solve_GSM(): matrix is not quadratic"
1642 unless ($rows1 == $cols1);
1643
1644 $n = $rows1;
1645
1646 croak "Math::MatrixReal::solve_GSM(): 1st vector is not a column vector"
1647 unless ($cols2 == 1);
1648
1649 croak "Math::MatrixReal::solve_GSM(): 2nd vector is not a column vector"
1650 unless ($cols3 == 1);
1651
1652 croak "Math::MatrixReal::solve_GSM(): matrix and vector size mismatch"
1653 unless (($rows2 == $n) && ($rows3 == $n));
1654
1655 return() unless ($norm = $matrix->_init_iteration());
1656
1657 $xn_vector = $x0_vector->new($n,1);
1658
1659 $diff = $epsilon + 1;
1660 while ($diff >= $epsilon)
1661 {
1662 for ( $i = 0; $i < $n; $i++ )
1663 {
1664 $sum = $b_vector->[0][$i][0];
1665 for ( $j = 0; $j < $i; $j++ )
1666 {
1667 $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
1668 }
1669 for ( $j = ($i + 1); $j < $n; $j++ )
1670 {
1671 $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
1672 }
1673 $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
1674 }
1675 $x0_vector->subtract($x0_vector,$xn_vector);
1676 if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1677 else { $diff = $x0_vector->norm_max(); }
1678 for ( $i = 0; $i < $n; $i++ )
1679 {
1680 $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1681 }
1682 }
1683 return($xn_vector);
1684}
1685
1686sub solve_SSM # Single Step Method
1687{
1688 croak "Usage: \$xn_vector = \$matrix->solve_SSM(\$x0_vector,\$b_vector,\$epsilon);"
1689 if (@_ != 4);
1690
1691 my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
1692 my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1693 my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1694 my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1695 my($norm,$sum,$diff);
1696 my($xn_vector);
1697 my($i,$j,$n);
1698
1699 croak "Math::MatrixReal::solve_SSM(): matrix is not quadratic"
1700 unless ($rows1 == $cols1);
1701
1702 $n = $rows1;
1703
1704 croak "Math::MatrixReal::solve_SSM(): 1st vector is not a column vector"
1705 unless ($cols2 == 1);
1706
1707 croak "Math::MatrixReal::solve_SSM(): 2nd vector is not a column vector"
1708 unless ($cols3 == 1);
1709
1710 croak "Math::MatrixReal::solve_SSM(): matrix and vector size mismatch"
1711 unless (($rows2 == $n) && ($rows3 == $n));
1712
1713 return() unless ($norm = $matrix->_init_iteration());
1714
1715 $xn_vector = $x0_vector->new($n,1);
1716 $xn_vector->copy($x0_vector);
1717
1718 $diff = $epsilon + 1;
1719 while ($diff >= $epsilon)
1720 {
1721 for ( $i = 0; $i < $n; $i++ )
1722 {
1723 $sum = $b_vector->[0][$i][0];
1724 for ( $j = 0; $j < $i; $j++ )
1725 {
1726 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1727 }
1728 for ( $j = ($i + 1); $j < $n; $j++ )
1729 {
1730 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1731 }
1732 $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
1733 }
1734 $x0_vector->subtract($x0_vector,$xn_vector);
1735 if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1736 else { $diff = $x0_vector->norm_max(); }
1737 for ( $i = 0; $i < $n; $i++ )
1738 {
1739 $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1740 }
1741 }
1742 return($xn_vector);
1743}
1744
1745sub solve_RM # Relaxation Method
1746{
1747 croak "Usage: \$xn_vector = \$matrix->solve_RM(\$x0_vector,\$b_vector,\$weight,\$epsilon);"
1748 if (@_ != 5);
1749
1750 my($matrix,$x0_vector,$b_vector,$weight,$epsilon) = @_;
1751 my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1752 my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1753 my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1754 my($norm,$sum,$diff);
1755 my($xn_vector);
1756 my($i,$j,$n);
1757
1758 croak "Math::MatrixReal::solve_RM(): matrix is not quadratic"
1759 unless ($rows1 == $cols1);
1760
1761 $n = $rows1;
1762
1763 croak "Math::MatrixReal::solve_RM(): 1st vector is not a column vector"
1764 unless ($cols2 == 1);
1765
1766 croak "Math::MatrixReal::solve_RM(): 2nd vector is not a column vector"
1767 unless ($cols3 == 1);
1768
1769 croak "Math::MatrixReal::solve_RM(): matrix and vector size mismatch"
1770 unless (($rows2 == $n) && ($rows3 == $n));
1771
1772 return() unless ($norm = $matrix->_init_iteration());
1773
1774 $xn_vector = $x0_vector->new($n,1);
1775 $xn_vector->copy($x0_vector);
1776
1777 $diff = $epsilon + 1;
1778 while ($diff >= $epsilon)
1779 {
1780 for ( $i = 0; $i < $n; $i++ )
1781 {
1782 $sum = $b_vector->[0][$i][0];
1783 for ( $j = 0; $j < $i; $j++ )
1784 {
1785 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1786 }
1787 for ( $j = ($i + 1); $j < $n; $j++ )
1788 {
1789 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1790 }
1791 $xn_vector->[0][$i][0] = $weight * ( $sum / $matrix->[0][$i][$i] )
1792 + (1 - $weight) * $xn_vector->[0][$i][0];
1793 }
1794 $x0_vector->subtract($x0_vector,$xn_vector);
1795 if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1796 else { $diff = $x0_vector->norm_max(); }
1797 for ( $i = 0; $i < $n; $i++ )
1798 {
1799 $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1800 }
1801 }
1802 return($xn_vector);
1803}
1804
1805# Core householder reduction routine (when eigenvector
1806# are wanted).
1807# Adapted from: Numerical Recipes, 2nd edition.
1808sub _householder_vectors ($)
1809{
1810 my ($Q) = @_;
1811 my ($rows, $cols) = ($Q->[1], $Q->[2]);
1812
1813 # Creates tridiagonal
1814 # Set up tridiagonal needed elements
1815 my $d = []; # N Diagonal elements 0...n-1
1816 my $e = []; # N-1 Off-Diagonal elements 0...n-2
1817
1818 my @p = ();
1819 for (my $i = ($rows-1); $i > 1; $i--)
1820 {
1821 my $scale = 0.0;
1822 # Computes norm of one column (below diagonal)
1823 for (my $k = 0; $k < $i; $k++)
1824 {
1825 $scale += abs($Q->[0][$i][$k]);
1826 }
1827 if ($scale == 0.0)
1828 { # skip the transformation
1829 $e->[$i-1] = $Q->[0][$i][$i-1];
1830 }
1831 else
1832 {
1833 my $h = 0.0;
1834 for (my $k = 0; $k < $i; $k++)
1835 { # Used scaled Q for transformation
1836 $Q->[0][$i][$k] /= $scale;
1837 # Form sigma in h
1838 $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
1839 }
1840 my $t1 = $Q->[0][$i][$i-1];
1841 my $t2 = (($t1 >= 0.0) ? -sqrt($h) : sqrt($h));
1842 $e->[$i-1] = $scale * $t2; # Update off-diagonals
1843 $h -= $t1 * $t2;
1844 $Q->[0][$i][$i-1] -= $t2;
1845 my $f = 0.0;
1846 for (my $j = 0; $j < $i; $j++)
1847 {
1848 $Q->[0][$j][$i] = $Q->[0][$i][$j] / $h;
1849 my $g = 0.0;
1850 for (my $k = 0; $k <= $j; $k++)
1851 {
1852 $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
1853 }
1854 for (my $k = $j+1; $k < $i; $k++)
1855 {
1856 $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
1857 }
1858 # Form elements of P
1859 $p[$j] = $g / $h;
1860 $f += $p[$j] * $Q->[0][$i][$j];
1861 }
1862 my $hh = $f / ($h + $h);
1863 for (my $j = 0; $j < $i; $j++)
1864 {
1865 my $t3 = $Q->[0][$i][$j];
1866 my $t4 = $p[$j] - $hh * $t3;
1867 $p[$j] = $t4;
1868 for (my $k = 0; $k <= $j; $k++)
1869 {
1870 $Q->[0][$j][$k] -= $t3 * $p[$k]
1871 + $t4 * $Q->[0][$i][$k];
1872 }
1873 }
1874 }
1875 }
1876 # Updates for i == 0,1
1877 $e->[0] = $Q->[0][1][0];
1878 $d->[0] = $Q->[0][0][0]; # i==0
1879 $Q->[0][0][0] = 1.0;
1880 $d->[1] = $Q->[0][1][1]; # i==1
1881 $Q->[0][1][1] = 1.0;
1882 $Q->[0][1][0] = $Q->[0][0][1] = 0.0;
1883 for (my $i = 2; $i < $rows; $i++)
1884 {
1885 for (my $j = 0; $j < $i; $j++)
1886 {
1887 my $g = 0.0;
1888 for (my $k = 0; $k < $i; $k++)
1889 {
1890 $g += $Q->[0][$i][$k] * $Q->[0][$k][$j];
1891 }
1892 for (my $k = 0; $k < $i; $k++)
1893 {
1894 $Q->[0][$k][$j] -= $g * $Q->[0][$k][$i];
1895 }
1896 }
1897 $d->[$i] = $Q->[0][$i][$i];
1898 # Reset row and column of Q for next iteration
1899 $Q->[0][$i][$i] = 1.0;
1900 for (my $j = 0; $j < $i; $j++)
1901 {
1902 $Q->[0][$i][$j] = $Q->[0][$j][$i] = 0.0;
1903 }
1904 }
1905 return ($d, $e);
1906}
1907
1908# Computes sqrt(a*a + b*b), but more carefully...
1909sub _pythag ($$)
1910{
1911 my ($a, $b) = @_;
1912 my $aa = abs($a);
1913 my $ab = abs($b);
1914 if ($aa > $ab)
1915 {
1916 # NB: Not needed!: return 0.0 if ($aa == 0.0);
1917 my $t = $ab / $aa;
1918 return ($aa * sqrt(1.0 + $t*$t));
1919 }
1920 else
1921 {
1922 return 0.0 if ($ab == 0.0);
1923 my $t = $aa / $ab;
1924 return ($ab * sqrt(1.0 + $t*$t));
1925 }
1926}
1927
1928# QL algorithm with implicit shifts to determine the eigenvalues
1929# of a tridiagonal matrix. Internal routine.
1930sub _tridiagonal_QLimplicit
1931{
1932 my ($EV, $d, $e) = @_;
1933 my ($rows, $cols) = ($EV->[1], $EV->[2]);
1934
1935 $e->[$rows-1] = 0.0;
1936 # Start real computation
1937 for (my $l = 0; $l < $rows; $l++)
1938 {
1939 my $iter = 0;
1940 my $m;
1941 OUTER:
1942 do {
1943 for ($m = $l; $m < ($rows - 1); $m++)
1944 {
1945 my $dd = abs($d->[$m]) + abs($d->[$m+1]);
1946 last if ((abs($e->[$m]) + $dd) == $dd);
1947 }
1948 if ($m != $l)
1949 {
1950 ## why only allow 30 iterations?
1951 croak("Too many iterations!") if ($iter++ >= 30);
1952 my $g = ($d->[$l+1] - $d->[$l])
1953 / (2.0 * $e->[$l]);
1954 my $r = _pythag($g, 1.0);
1955 $g = $d->[$m] - $d->[$l]
1956 + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
1957 my ($p,$s,$c) = (0.0, 1.0,1.0);
1958 for (my $i = ($m-1); $i >= $l; $i--)
1959 {
1960 my $ii = $i + 1;
1961 my $f = $s * $e->[$i];
1962 my $t = _pythag($f, $g);
1963 $e->[$ii] = $t;
1964 if ($t == 0.0)
1965 {
1966 $d->[$ii] -= $p;
1967 $e->[$m] = 0.0;
1968 next OUTER;
1969 }
1970 my $b = $c * $e->[$i];
1971 $s = $f / $t;
1972 $c = $g / $t;
1973 $g = $d->[$ii] - $p;
1974 my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
1975 $p = $s * $t2;
1976 $d->[$ii] = $g + $p;
1977 $g = $c * $t2 - $b;
1978 for (my $k = 0; $k < $rows; $k++)
1979 {
1980 my $t1 = $EV->[0][$k][$ii];
1981 my $t2 = $EV->[0][$k][$i];
1982 $EV->[0][$k][$ii] = $s * $t2 + $c * $t1;
1983 $EV->[0][$k][$i] = $c * $t2 - $s * $t1;
1984 }
1985 }
1986 $d->[$l] -= $p;
1987 $e->[$l] = $g;
1988 $e->[$m] = 0.0;
1989 }
1990 } while ($m != $l);
1991 }
1992 return;
1993}
1994
1995# Core householder reduction routine (when eigenvector
1996# are NOT wanted).
1997sub _householder_values ($)
1998{
1999 my ($Q) = @_; # NB: Q is destroyed on output...
2000 my ($rows, $cols) = ($Q->[1], $Q->[2]);
2001
2002 # Creates tridiagonal
2003 # Set up tridiagonal needed elements
2004 my $d = []; # N Diagonal elements 0...n-1
2005 my $e = []; # N-1 Off-Diagonal elements 0...n-2
2006
2007 my @p = ();
2008 for (my $i = ($rows - 1); $i > 1; $i--)
2009 {
2010 my $scale = 0.0;
2011 for (my $k = 0; $k < $i; $k++)
2012 {
2013 $scale += abs($Q->[0][$i][$k]);
2014 }
2015 if ($scale == 0.0)
2016 { # skip the transformation
2017 $e->[$i-1] = $Q->[0][$i][$i-1];
2018 }
2019 else
2020 {
2021 my $h = 0.0;
2022 for (my $k = 0; $k < $i; $k++)
2023 { # Used scaled Q for transformation
2024 $Q->[0][$i][$k] /= $scale;
2025 # Form sigma in h
2026 $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
2027 }
2028 my $t = $Q->[0][$i][$i-1];
2029 my $t2 = (($t >= 0.0) ? -sqrt($h) : sqrt($h));
2030 $e->[$i-1] = $scale * $t2; # Updates off-diagonal
2031 $h -= $t * $t2;
2032 $Q->[0][$i][$i-1] -= $t2;
2033 my $f = 0.0;
2034 for (my $j = 0; $j < $i; $j++)
2035 {
2036 my $g = 0.0;
2037 for (my $k = 0; $k <= $j; $k++)
2038 {
2039 $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
2040 }
2041 for (my $k = $j+1; $k < $i; $k++)
2042 {
2043 $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
2044 }
2045 # Form elements of P
2046 $p[$j] = $g / $h;
2047 $f += $p[$j] * $Q->[0][$i][$j];
2048 }
2049 my $hh = $f / ($h + $h);
2050 for (my $j = 0; $j < $i; $j++)
2051 {
2052 my $t = $Q->[0][$i][$j];
2053 my $g = $p[$j] - $hh * $t;
2054 $p[$j] = $g;
2055 for (my $k = 0; $k <= $j; $k++)
2056 {
2057 $Q->[0][$j][$k] -= $t * $p[$k]
2058 + $g * $Q->[0][$i][$k];
2059 }
2060 }
2061 }
2062 }
2063 # Updates for i==1
2064 $e->[0] = $Q->[0][1][0];
2065 # Updates diagonal elements
2066 for (my $i = 0; $i < $rows; $i++)
2067 {
2068 $d->[$i] = $Q->[0][$i][$i];
2069 }
2070 return ($d, $e);
2071}
2072
2073# QL algorithm with implicit shifts to determine the
2074# eigenvalues ONLY. This is O(N^2) only...
2075sub _tridiagonal_QLimplicit_values
2076{
2077 my ($M, $d, $e) = @_; # NB: M is not touched...
2078 my ($rows, $cols) = ($M->[1], $M->[2]);
2079
2080 $e->[$rows-1] = 0.0;
2081 # Start real computation
2082 for (my $l = 0; $l < $rows; $l++)
2083 {
2084 my $iter = 0;
2085 my $m;
2086 OUTER:
2087 do {
2088 for ($m = $l; $m < ($rows - 1); $m++)
2089 {
2090 my $dd = abs($d->[$m]) + abs($d->[$m+1]);
2091 last if ((abs($e->[$m]) + $dd) == $dd);
2092 }
2093 if ($m != $l)
2094 {
2095 croak("Too many iterations!") if ($iter++ >= 30);
2096 my $g = ($d->[$l+1] - $d->[$l])
2097 / (2.0 * $e->[$l]);
2098 my $r = _pythag($g, 1.0);
2099 $g = $d->[$m] - $d->[$l]
2100 + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
2101 my ($p,$s,$c) = (0.0, 1.0,1.0);
2102 for (my $i = ($m-1); $i >= $l; $i--)
2103 {
2104 my $ii = $i + 1;
2105 my $f = $s * $e->[$i];
2106 my $t = _pythag($f, $g);
2107 $e->[$ii] = $t;
2108 if ($t == 0.0)
2109 {
2110 $d->[$ii] -= $p;
2111 $e->[$m] = 0.0;
2112 next OUTER;
2113 }
2114 my $b = $c * $e->[$i];
2115 $s = $f / $t;
2116 $c = $g / $t;
2117 $g = $d->[$ii] - $p;
2118 my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
2119 $p = $s * $t2;
2120 $d->[$ii] = $g + $p;
2121 $g = $c * $t2 - $b;
2122 }
2123 $d->[$l] -= $p;
2124 $e->[$l] = $g;
2125 $e->[$m] = 0.0;
2126 }
2127 } while ($m != $l);
2128 }
2129 return;
2130}
2131
2132# Householder reduction of a real, symmetric matrix A.
2133# Returns a tridiagonal matrix T and the orthogonal matrix
2134# Q effecting the transformation between A and T.
2135sub householder ($)
2136{
2137 my ($A) = @_;
2138 my ($rows, $cols) = ($A->[1], $A->[2]);
2139
2140 croak "Matrix is not quadratic"
2141 unless ($rows = $cols);
2142 croak "Matrix is not symmetric"
2143 unless ($A->is_symmetric());
2144
2145 # Copy given matrix TODO: study if we should do in-place modification
2146 my $Q = $A->clone();
2147
2148 # Do the computation of tridiagonal elements and of
2149 # transformation matrix
2150 my ($diag, $offdiag) = $Q->_householder_vectors();
2151
2152 # Creates the tridiagonal matrix
2153 my $T = $A->shadow();
2154 for (my $i = 0; $i < $rows; $i++)
2155 { # Set diagonal
2156 $T->[0][$i][$i] = $diag->[$i];
2157 }
2158 for (my $i = 0; $i < ($rows-1); $i++)
2159 { # Set off diagonals
2160 $T->[0][$i+1][$i] = $offdiag->[$i];
2161 $T->[0][$i][$i+1] = $offdiag->[$i];
2162 }
2163 return ($T, $Q);
2164}
2165
2166# QL algorithm with implicit shifts to determine the eigenvalues
2167# and eigenvectors of a real tridiagonal matrix - or of a matrix
2168# previously reduced to tridiagonal form.
2169sub tri_diagonalize ($;$)
2170{
2171 my ($T,$Q) = @_; # Q may be 0 if the original matrix is really tridiagonal
2172
2173 my ($rows, $cols) = ($T->[1], $T->[2]);
2174
2175 croak "Matrix is not quadratic"
2176 unless ($rows = $cols);
2177 croak "Matrix is not tridiagonal"
2178 unless ($T->is_tridiagonal()); # DONE
2179
2180 my $EV;
2181 # Obtain/Creates the todo eigenvectors matrix
2182 if ($Q)
2183 {
2184 $EV = $Q->clone();
2185 }
2186 else
2187 {
2188 $EV = $T->shadow();
2189 $EV->one();
2190 }
2191 # Allocates diagonal vector
2192 my $diag = [ ];
2193 # Initializes it with T
2194 for (my $i = 0; $i < $rows; $i++)
2195 {
2196 $diag->[$i] = $T->[0][$i][$i];
2197 }
2198 # Allocate temporary vector for off-diagonal elements
2199 my $offdiag = [ ];
2200 for (my $i = 1; $i < $rows; $i++)
2201 {
2202 $offdiag->[$i-1] = $T->[0][$i][$i-1];
2203 }
2204
2205 # Calls the calculus routine
2206 $EV->_tridiagonal_QLimplicit($diag, $offdiag);
2207
2208 # Allocate eigenvalues vector
2209 my $v = Math::MatrixReal->new($rows,1);
2210 # Fills it
2211 for (my $i = 0; $i < $rows; $i++)
2212 {
2213 $v->[0][$i][0] = $diag->[$i];
2214 }
2215 return ($v, $EV);
2216}
2217
2218# Main routine for diagonalization of a real symmetric
2219# matrix M. Operates by transforming M into a tridiagonal
2220# matrix and then obtaining the eigenvalues and eigenvectors
2221# for that matrix (taking into account the transformation to
2222# tridiagonal).
2223sub sym_diagonalize ($)
2224{
2225 my ($M) = @_;
2226 my ($rows, $cols) = ($M->[1], $M->[2]);
2227
2228 croak "Matrix is not quadratic"
2229 unless ($rows = $cols);
2230 croak "Matrix is not symmetric"
2231 unless ($M->is_symmetric());
2232
2233 # Copy initial matrix
2234 # TODO: study if we should allow in-place modification
2235 my $VEC = $M->clone();
2236
2237 # Do the computation of tridiagonal elements and of
2238 # transformation matrix
2239 my ($diag, $offdiag) = $VEC->_householder_vectors();
2240 # Calls the calculus routine for diagonalization
2241 $VEC->_tridiagonal_QLimplicit($diag, $offdiag);
2242
2243 # Allocate eigenvalues vector
2244 my $val = Math::MatrixReal->new($rows,1);
2245 # Fills it
2246 for (my $i = 0; $i < $rows; $i++)
2247 {
2248 $val->[0][$i][0] = $diag->[$i];
2249 }
2250 return ($val, $VEC);
2251}
2252
2253# Householder reduction of a real, symmetric matrix A.
2254# Returns a tridiagonal matrix T equivalent to A.
2255sub householder_tridiagonal ($)
2256{
2257 my ($A) = @_;
2258 my ($rows, $cols) = ($A->[1], $A->[2]);
2259
2260 croak "Matrix is not quadratic"
2261 unless ($rows = $cols);
2262 croak "Matrix is not symmetric"
2263 unless ($A->is_symmetric());
2264
2265 # Copy given matrix
2266 my $Q = $A->clone();
2267
2268 # Do the computation of tridiagonal elements and of
2269 # transformation matrix
2270 # Q is destroyed after reduction
2271 my ($diag, $offdiag) = $Q->_householder_values();
2272
2273 # Creates the tridiagonal matrix in Q (avoid allocation)
2274 my $T = $Q;
2275 $T->zero();
2276 for (my $i = 0; $i < $rows; $i++)
2277 { # Set diagonal
2278 $T->[0][$i][$i] = $diag->[$i];
2279 }
2280 for (my $i = 0; $i < ($rows-1); $i++)
2281 { # Set off diagonals
2282 $T->[0][$i+1][$i] = $offdiag->[$i];
2283 $T->[0][$i][$i+1] = $offdiag->[$i];
2284 }
2285 return $T;
2286}
2287
2288# QL algorithm with implicit shifts to determine ONLY
2289# the eigenvalues a real tridiagonal matrix - or of a
2290# matrix previously reduced to tridiagonal form.
2291sub tri_eigenvalues ($;$)
2292{
2293 my ($T) = @_;
2294 my ($rows, $cols) = ($T->[1], $T->[2]);
2295
2296 croak "Matrix is not quadratic"
2297 unless ($rows = $cols);
2298 croak "Matrix is not tridiagonal"
2299 unless ($T->is_tridiagonal() ); # DONE
2300
2301 # Allocates diagonal vector
2302 my $diag = [ ];
2303 # Initializes it with T
2304 for (my $i = 0; $i < $rows; $i++)
2305 {
2306 $diag->[$i] = $T->[0][$i][$i];
2307 }
2308 # Allocate temporary vector for off-diagonal elements
2309 my $offdiag = [ ];
2310 for (my $i = 1; $i < $rows; $i++)
2311 {
2312 $offdiag->[$i-1] = $T->[0][$i][$i-1];
2313 }
2314
2315 # Calls the calculus routine (T is not touched)
2316 $T->_tridiagonal_QLimplicit_values($diag, $offdiag);
2317
2318 # Allocate eigenvalues vector
2319 my $v = Math::MatrixReal->new($rows,1);
2320 # Fills it
2321 for (my $i = 0; $i < $rows; $i++)
2322 {
2323 $v->[0][$i][0] = $diag->[$i];
2324 }
2325 return $v;
2326}
2327
2328## more general routine than sym_eigenvalues
2329sub eigenvalues ($){
2330 my ($matrix) = @_;
2331 my ($rows,$cols) = $matrix->dim();
2332 my $i=0;
2333
2334 #return sym_eigenvalues($matrix) if $matrix->is_symmetric();
2335
2336 croak "Matrix is not quadratic" unless ($rows == $cols);
2337
2338 if($matrix->is_upper_triangular() || $matrix->is_lower_triangular() ){
2339 my $l = Math::MatrixReal->new($rows,1);
2340 for(;$i < $rows; $i++ ){
2341 $l->[0][$i][0] = $matrix->[0][$i][$i];
2342 }
2343 return $l;
2344 }
2345
2346 return sym_eigenvalues($matrix) if $matrix->is_symmetric();
2347
2348 carp "Math::MatrixReal::eigenvalues(): Matrix is not symmetric or triangular";
2349 return undef;
2350
2351}
2352# Main routine for diagonalization of a real symmetric
2353# matrix M. Operates by transforming M into a tridiagonal
2354# matrix and then obtaining the eigenvalues and eigenvectors
2355# for that matrix (taking into account the transformation to
2356# tridiagonal).
2357sub sym_eigenvalues ($)
2358{
2359 my ($M) = @_;
2360 my ($rows, $cols) = ($M->[1], $M->[2]);
2361
2362 croak "Matrix is not quadratic"
2363 unless ($rows == $cols);
2364 croak "Matrix is not symmetric"
2365 unless ($M->is_symmetric());
2366
2367 # Copy matrix in temporary
2368 my $A = $M->clone();
2369 # Do the computation of tridiagonal elements and of
2370 # transformation matrix. A is destroyed
2371 my ($diag, $offdiag) = $A->_householder_values();
2372 # Calls the calculus routine for diagonalization
2373 # (M is not touched)
2374 $M->_tridiagonal_QLimplicit_values($diag, $offdiag);
2375
2376 # Allocate eigenvalues vector
2377 my $val = Math::MatrixReal->new($rows,1);
2378 # Fills it
2379 for (my $i = 0; $i < $rows; $i++)
2380 {
2381 $val->[0][$i][0] = $diag->[$i];
2382 }
2383 return $val;
2384}
2385#TODO: docs+test
2386sub is_positive_definite {
2387 my ($matrix) = @_;
2388 my ($r,$c) = $matrix->dim;
2389 my $ev = $matrix->eigenvalues;
2390 my $pos = 1;
2391 $ev->each(sub { my $x = shift; if ($x <= 0){ $pos=0;return; } } );
2392 return $pos;
2393}
2394
2395sub is_row_vector {
2396 my ($m) = @_;
2397 my ($r,$c) = $m->dim;
2398 return 1 if ($r == 1);
2399}
2400sub is_col_vector {
2401 my ($m) = @_;
2402 my ($r,$c) = $m->dim;
2403 return 1 if ($c == 1);
2404}
2405
2406sub is_orthogonal($) {
2407 my ($matrix) = @_;
2408 ##croak "Math::MatrixReal::is_orthogonal(): Matrix is not quadratic"
2409
2410 return 0 unless( $matrix->is_quadratic() );
2411
2412 my $one = $matrix->shadow();
2413 $one->one;
2414
2415 abs(~$matrix * $matrix - $one) < 1e-12 ? return 1 : return 0;
2416
2417}
2418
2419sub is_positive($) {
2420 my ($m) = @_;
2421 my $pos = 1;
2422 $m->each( sub { if( (shift) <= 0){ $pos = 0;return;} } );
2423 return $pos;
2424}
2425sub is_negative($) {
2426 my ($m) = @_;
2427 my $neg = 1;
2428 $m->each( sub { if( (shift) >= 0){ $neg = 0;return;} } );
2429 return $neg;
2430}
2431
2432
2433sub is_periodic($$) {
2434 my ($m,$k) = @_;
2435 return 0 unless $m->is_quadratic();
2436 abs($m**(int($k)+1) - $m) < 1e-12 ? return 1 : return 0;
2437}
2438sub is_idempotent($) {
2439 return (shift)->is_periodic(1);
2440}
2441
2442# Boolean check routine to see if a matrix is
2443# symmetric
2444sub is_symmetric ($)
2445{
2446 my ($M) = @_;
2447 my ($rows, $cols) = ($M->[1], $M->[2]);
2448 # if it is not quadratic it cannot be symmetric...
2449 return 0 unless ($rows == $cols);
2450 # skip when $i=$j?
2451 for (my $i = 1; $i < $rows; $i++)
2452 {
2453 for (my $j = 0; $j < $i; $j++)
2454 {
2455 return 0 unless ($M->[0][$i][$j] == $M->[0][$j][$i]);
2456 }
2457 }
2458 return 1;
2459}
2460# Boolean check to see if matrix is tridiagonal
2461sub is_tridiagonal ($) {
2462 my ($M) = @_;
2463 my ($rows,$cols) = ($M->[1],$M->[2]);
2464 my ($i,$j) = (0,0);
2465 # if it is not quadratic it cannot be tridiag
2466 return 0 unless ($rows == $cols);
2467
2468 for(;$i < $rows; $i++ ){
2469 for(;$j < $cols; $j++ ){
2470 #print "debug: testing $i,$j = " . $M->[0][$i][$j] . "\n";
2471 # skip diag and diag+-1
2472 next if ($i == $j);
2473 next if ($i+1 == $j);
2474 next if ($i-1 == $j);
2475 return 0 if $M->[0][$i][$j];
2476 }
2477 $j = 0;
2478 }
2479 return 1;
2480}
2481# Boolean check to see if matrix is upper triangular
2482# i.e all nonzero elements are above main diagonal
2483sub is_upper_triangular {
2484 my ($M) = @_;
2485 my ($rows,$cols) = $M->dim();
2486 my ($i,$j) = (1,0);
2487 return 0 unless ($rows == $cols);
2488
2489 for(;$i < $rows; $i++ ){
2490 for(;$j < $cols;$j++ ){
2491 next if ($i <= $j);
2492 return 0 if $M->[0][$i][$j];
2493 }
2494 $j = 0;
2495 }
2496 return 1;
2497}
2498# Boolean check to see if matrix is lower triangular
2499# i.e all nonzero elements are lower main diagonal
2500sub is_lower_triangular {
2501 my ($M) = @_;
2502 my ($rows,$cols) = $M->dim();
2503 my ($i,$j) = (0,1);
2504 return 0 unless ($rows == $cols);
2505
2506 for(;$i < $rows; $i++ ){
2507 for(;$j < $cols;$j++ ){
2508 next if ($i >= $j);
2509 return 0 if $M->[0][$i][$j];
2510 }
2511 $j = 0;
2512 }
2513 return 1;
2514}
2515
2516# Boolean check to see if matrix is diagonal
2517sub is_diagonal ($) {
2518 my ($M) = @_;
2519 my ($rows,$cols) = ($M->[1],$M->[2]);
2520 my ($i,$j) = (0,0);
2521 return 0 unless ($rows == $cols );
2522 for(;$i < $rows; $i++ ){
2523 for(;$j < $cols; $j++ ){
2524 # skip diag elements
2525 next if ($i == $j);
2526 return 0 if $M->[0][$i][$j];
2527 }
2528 $j = 0;
2529 }
2530 return 1;
2531}
2532sub is_quadratic ($) {
2533 croak "Usage: \$matrix->is_quadratic()" unless (@_ == 1);
2534 my ($matrix) = @_;
2535 $matrix->[1] == $matrix->[2] ? return 1 : return 0;
2536}
2537
2538sub is_square($) {
2539 croak "Usage: \$matrix->is_square()" unless (@_ == 1);
2540 return (shift)->is_quadratic();
2541}
2542
2543sub is_LR($) {
2544 croak "Usage: \$matrix->is_LR()" unless (@_ == 1);
2545 return (shift)->[3] ? 1 : 0;
2546}
2547###
2548sub is_normal{
2549 my ($matrix) = @_;
2550 my ($rows,$cols) = $matrix->dim;
2551 return 0 unless ($rows == $cols);
2552
2553 return 1 if ( ~$matrix * $matrix - $matrix * ~$matrix < 1e-8 );
2554 return 0;
2555
2556}
2557sub is_skew_symmetric{
2558 my ($m) = @_;
2559 my ($rows, $cols) = $m->dim;
2560 # if it is not quadratic it cannot be skew symmetric...
2561 return 0 unless ($rows == $cols);
2562 for (my $i = 1; $i < $rows; $i++) {
2563 for (my $j = 0; $j < $i; $j++) {
2564 return 0 unless ($m->[0][$i][$j] == -$m->[0][$j][$i]);
2565 }
2566 }
2567 return 1;
2568
2569}
2570####
2571sub is_gramian{
2572 my ($m) = @_;
2573 my ($rows,$cols) = $m->dim;
2574 my $neg=0;
2575 # gramian matrix must be symmetric
2576 return 0 unless $m->is_symmetric;
2577
2578 # must have all non-negative eigenvalues
2579 my $ev = $m->eigenvalues;
2580 $ev->each(sub { $neg++ if ((shift)<0) } );
2581
2582 return $neg ? 0 : 1;
2583}
2584sub is_binary{
2585 my ($m) = @_;
2586 my ($rows, $cols) = $m->dim;
2587 for (my $i = 0; $i < $rows; $i++) {
2588 for (my $j = 0; $j < $cols; $j++) {
2589 return 0 unless ($m->[0][$i][$j] == 1 || $m->[0][$i][$j] == 0);
2590 }
2591 }
2592 return 1;
2593
2594}
2595sub as_scilab {
2596 return (shift)->as_matlab;
2597}
2598
2599sub as_matlab {
2600 my ($m) = shift;
2601 my %args = (
2602 format => "%s",
2603 name => "",
2604 semi => 0,
2605 @_);
2606 my ($row,$col) = $m->dim;
2607 my $s = "";
2608
2609 if( $args{name} ){
2610 $s = "$args{name} = ";
2611 }
2612 $s .= "[";
2613 $m->each(
2614 sub { my($x,$i,$j) = @_;
2615 $s .= sprintf(" $args{format}",$x);
2616 $s .= ";\n" if( $j == $col && $i != $row);
2617 }
2618 );
2619 $s .= "]";
2620 $s .= ";" if $args{semi};
2621 return $s;
2622}
2623#TODO: docs+test
2624sub as_yacas{
2625 my ($m) = shift;
2626 my %args = (
2627 format => "%s",
2628 name => "",
2629 semi => 0,
2630 @_);
2631 my ($row,$col) = $m->dim;
2632 my $s = "";
2633
2634 if( $args{name} ){
2635 $s = "$args{name} := ";
2636 }
2637 $s .= "{";
2638
2639 $m->each(
2640 sub { my($x,$i,$j) = @_;
2641 $s .= "{" if ($j == 1);
2642 $s .= sprintf("$args{format}",$x);
2643 $s .= "," if( $j != $col );
2644 $s .= "}," if ($j == $col && $i != $row);
2645 }
2646 );
2647 $s .= "}}";
2648
2649 return $s;
2650}
2651#TODO: any ideas for a good test?
2652sub as_latex{
2653 my ($m) = shift;
2654 my %args = (
2655 format => "%s",
2656 name => "",
2657 align => "c",
2658 @_);
2659
2660
2661 my ($row,$col) = $m->dim;
2662 my $inside;
2663 my $s = <<LATEX;
2664\$
2665\\left( \\begin{array}{%COLS%}
2666%INSIDE%\\end{array} \\right)
2667\$
2668LATEX
2669 $args{align} = lc $args{align};
2670 if( $args{align} !~ m/^(c|l|r)$/ ){
2671 croak "Math::MatrixReal::as_latex(): Invalid alignment '$args{align}'";
2672 }
2673
2674 $s =~ s/%COLS%/$args{align} x $col/em;
2675
2676 if( $args{name} ){
2677 $s = "\$$args{name} = \$ " . $s;
2678 }
2679 $m->each(
2680 sub { my($x,$i,$j)=@_;
2681
2682 $x = sprintf($args{format},$x);
2683
2684 # last element in each row gets a \\
2685 if ($j == $col && $i != $row){
2686 $inside .= "$x \\\\"."\n";
2687 # the annoying last line has neither
2688 } elsif( $j == $col && $i == $row){
2689 $inside .= "$x\n";
2690 } else {
2691 $inside .= "$x&";
2692 }
2693 }
2694 );
2695 $s =~ s/%INSIDE%/$inside/gm;
2696 return $s;
2697}
2698####
2699sub spectral_radius {
2700 my ($matrix) = @_;
2701 my ($r,$c) = $matrix->dim;
2702 my $ev = $matrix->eigenvalues;
2703 my $radius=0;
2704 $ev->each(sub { my $x = shift; $radius = $x if (abs($x) > $radius); } );
2705 return $radius;
2706}
2707
2708
2709 ########################################
2710 # #
2711 # define overloaded operators section: #
2712 # #
2713 ########################################
2714
2715sub _negate
2716{
2717 my($object,$argument,$flag) = @_;
2718# my($name) = "neg";
2719 my($temp);
2720
2721 $temp = $object->new($object->[1],$object->[2]);
2722 $temp->negate($object);
2723 return($temp);
2724}
2725
2726sub _transpose
2727{
2728 my($object,$argument,$flag) = @_;
2729# my($name) = "'~'";
2730 my($temp);
2731
2732 $temp = $object->new($object->[2],$object->[1]);
2733 $temp->transpose($object);
2734 return($temp);
2735
2736}
2737
2738sub _boolean
2739{
2740 my($object,$argument,$flag) = @_;
2741# my($name) = "bool";
2742 my($rows,$cols) = ($object->[1],$object->[2]);
2743 my($i,$j,$result);
2744
2745 #TODO: ugly...
2746 $result = 0;
2747 BOOL:
2748 for ( $i = 0; $i < $rows; $i++ )
2749 {
2750 for ( $j = 0; $j < $cols; $j++ )
2751 {
2752 if ($object->[0][$i][$j] != 0)
2753 {
2754 $result = 1;
2755 last BOOL;
2756 }
2757 }
2758 }
2759 return($result);
2760}
2761#TODO: ugly copy+paste
2762sub _not_boolean
2763{
2764 my($object,$argument,$flag) = @_;
2765# my($name) = "'!'";
2766 my($rows,$cols) = ($object->[1],$object->[2]);
2767 my($i,$j,$result);
2768
2769 $result = 1;
2770 NOTBOOL:
2771 for ( $i = 0; $i < $rows; $i++ )
2772 {
2773 for ( $j = 0; $j < $cols; $j++ )
2774 {
2775 if ($object->[0][$i][$j] != 0)
2776 {
2777 $result = 0;
2778 last NOTBOOL;
2779 }
2780 }
2781 }
2782 return($result);
2783}
2784
2785sub _stringify
2786{
2787 my($object,$argument,$flag) = @_;
2788# my($name) = '""';
2789 my($rows,$cols) = ($object->[1],$object->[2]);
2790 my($i,$j,$s);
2791
2792 $s = '';
2793 for ( $i = 0; $i < $rows; $i++ )
2794 {
2795 $s .= "[ ";
2796 for ( $j = 0; $j < $cols; $j++ )
2797 {
2798 $s .= sprintf("% #-19.12E ", $object->[0][$i][$j]);
2799 }
2800 $s .= "]\n";
2801 }
2802 return($s);
2803}
2804
2805sub _norm
2806{
2807 my($object,$argument,$flag) = @_;
2808# my($name) = "abs";
2809
2810 return( $object->norm_one() );
2811}
2812
2813sub _add
2814{
2815 my($object,$argument,$flag) = @_;
2816 my($name) = "'+'";
2817 my($temp);
2818
2819 if ((defined $argument) && ref($argument) &&
2820 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2821 {
2822 if (defined $flag)
2823 {
2824 $temp = $object->new($object->[1],$object->[2]);
2825 $temp->add($object,$argument);
2826 return($temp);
2827 }
2828 else
2829 {
2830 $object->add($object,$argument);
2831 return($object);
2832 }
2833 }
2834 else
2835 {
2836 croak "Math::MatrixReal $name: wrong argument type";
2837 }
2838}
2839
2840sub _subtract
2841{
2842 my($object,$argument,$flag) = @_;
2843 my($name) = "'-'";
2844 my($temp);
2845
2846 if ((defined $argument) && ref($argument) &&
2847 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2848 {
2849 if (defined $flag)
2850 {
2851 $temp = $object->new($object->[1],$object->[2]);
2852 if ($flag) { $temp->subtract($argument,$object); }
2853 else { $temp->subtract($object,$argument); }
2854 return($temp);
2855 }
2856 else
2857 {
2858 $object->subtract($object,$argument);
2859 return($object);
2860 }
2861 }
2862 else
2863 {
2864 croak "Math::MatrixReal $name: wrong argument type";
2865 }
2866}
2867
2868sub _exponent
2869{
2870 my($matrix,$argument,$flag) = @_;
2871 my($rows,$cols) = ($matrix->[1],$matrix->[2]);
2872 my($name) = "'**'";
2873
2874 return $matrix->exponent( $argument );
2875}
2876
2877
2878sub _multiply
2879{
2880 my($object,$argument,$flag) = @_;
2881 my($name) = "'*'";
2882 my($temp);
2883
2884 if ((defined $argument) && ref($argument) &&
2885 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2886 {
2887 if ((defined $flag) && $flag)
2888 {
2889 return( multiply($argument,$object) );
2890 }
2891 else
2892 {
2893 return( multiply($object,$argument) );
2894 }
2895 }
2896 elsif ((defined $argument) && !(ref($argument)))
2897 {
2898 if (defined $flag)
2899 {
2900 $temp = $object->new($object->[1],$object->[2]);
2901 $temp->multiply_scalar($object,$argument);
2902 return($temp);
2903 }
2904 else
2905 {
2906 $object->multiply_scalar($object,$argument);
2907 return($object);
2908 }
2909 }
2910 else
2911 {
2912 croak "Math::MatrixReal $name: wrong argument type";
2913 }
2914}
2915
2916sub _assign_add
2917{
2918 my($object,$argument,$flag) = @_;
2919# my($name) = "'+='";
2920
2921 return( &_add($object,$argument,undef) );
2922}
2923
2924sub _assign_subtract
2925{
2926 my($object,$argument,$flag) = @_;
2927# my($name) = "'-='";
2928
2929 return( &_subtract($object,$argument,undef) );
2930}
2931
2932sub _assign_multiply
2933{
2934 my($object,$argument,$flag) = @_;
2935# my($name) = "'*='";
2936
2937 return( &_multiply($object,$argument,undef) );
2938}
2939
2940sub _assign_exponent {
2941 my($object,$arg,$flag) = @_;
2942 return ( &_exponent($object,$arg,undef) );
2943}
2944
2945sub _equal
2946{
2947 my($object,$argument,$flag) = @_;
2948 my($name) = "'=='";
2949 my($rows,$cols) = ($object->[1],$object->[2]);
2950 my($i,$j,$result);
2951
2952 if ((defined $argument) && ref($argument) &&
2953 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2954 {
2955 $result = 1;
2956 EQUAL:
2957 for ( $i = 0; $i < $rows; $i++ )
2958 {
2959 for ( $j = 0; $j < $cols; $j++ )
2960 {
2961 if ($object->[0][$i][$j] != $argument->[0][$i][$j])
2962 {
2963 $result = 0;
2964 last EQUAL;
2965 }
2966 }
2967 }
2968 return($result);
2969 }
2970 else
2971 {
2972 croak "Math::MatrixReal $name: wrong argument type";
2973 }
2974}
2975
2976sub _not_equal
2977{
2978 my($object,$argument,$flag) = @_;
2979 my($name) = "'!='";
2980 my($rows,$cols) = ($object->[1],$object->[2]);
2981 my($i,$j,$result);
2982
2983 if ((defined $argument) && ref($argument) &&
2984 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2985 {
2986 $result = 0;
2987 NOTEQUAL:
2988 for ( $i = 0; $i < $rows; $i++ )
2989 {
2990 for ( $j = 0; $j < $cols; $j++ )
2991 {
2992 if ($object->[0][$i][$j] != $argument->[0][$i][$j])
2993 {
2994 $result = 1;
2995 last NOTEQUAL;
2996 }
2997 }
2998 }
2999 return($result);
3000 }
3001 else
3002 {
3003 croak "Math::MatrixReal $name: wrong argument type";
3004 }
3005}
3006
3007sub _less_than
3008{
3009 my($object,$argument,$flag) = @_;
3010 my($name) = "'<'";
3011
3012 if ((defined $argument) && ref($argument) &&
3013 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
3014 {
3015 if ((defined $flag) && $flag)
3016 {
3017 return( $argument->norm_one() < $object->norm_one() );
3018 }
3019 else
3020 {
3021 return( $object->norm_one() < $argument->norm_one() );
3022 }
3023 }
3024 elsif ((defined $argument) && !(ref($argument)))
3025 {
3026 if ((defined $flag) && $flag)
3027 {
3028 return( abs($argument) < $object->norm_one() );
3029 }
3030 else
3031 {
3032 return( $object->norm_one() < abs($argument) );
3033 }
3034 }
3035 else
3036 {
3037 croak "Math::MatrixReal $name: wrong argument type";
3038 }
3039}
3040
3041sub _less_than_or_equal
3042{
3043 my($object,$argument,$flag) = @_;
3044 my($name) = "'<='";
3045
3046 if ((defined $argument) && ref($argument) &&
3047 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
3048 {
3049 if ((defined $flag) && $flag)
3050 {
3051 return( $argument->norm_one() <= $object->norm_one() );
3052 }
3053 else
3054 {
3055 return( $object->norm_one() <= $argument->norm_one() );
3056 }
3057 }
3058 elsif ((defined $argument) && !(ref($argument)))
3059 {
3060 if ((defined $flag) && $flag)
3061 {
3062 return( abs($argument) <= $object->norm_one() );
3063 }
3064 else
3065 {
3066 return( $object->norm_one() <= abs($argument) );
3067 }
3068 }
3069 else
3070 {
3071 croak "Math::MatrixReal $name: wrong argument type";
3072 }
3073}
3074
3075sub _greater_than
3076{
3077 my($object,$argument,$flag) = @_;
3078 my($name) = "'>'";
3079
3080 if ((defined $argument) && ref($argument) &&
3081 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
3082 {
3083 if ((defined $flag) && $flag)
3084 {
3085 return( $argument->norm_one() > $object->norm_one() );
3086 }
3087 else
3088 {
3089 return( $object->norm_one() > $argument->norm_one() );
3090 }
3091 }
3092 elsif ((defined $argument) && !(ref($argument)))
3093 {
3094 if ((defined $flag) && $flag)
3095 {
3096 return( abs($argument) > $object->norm_one() );
3097 }
3098 else
3099 {
3100 return( $object->norm_one() > abs($argument) );
3101 }
3102 }
3103 else
3104 {
3105 croak "Math::MatrixReal $name: wrong argument type";
3106 }
3107}
3108
3109sub _greater_than_or_equal
3110{
3111 my($object,$argument,$flag) = @_;
3112 my($name) = "'>='";
3113
3114 if ((defined $argument) && ref($argument) &&
3115 (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
3116 {
3117 if ((defined $flag) && $flag)
3118 {
3119 return( $argument->norm_one() >= $object->norm_one() );
3120 }
3121 else
3122 {
3123 return( $object->norm_one() >= $argument->norm_one() );
3124 }
3125 }
3126 elsif ((defined $argument) && !(ref($argument)))
3127 {
3128 if ((defined $flag) && $flag)
3129 {
3130 return( abs($argument) >= $object->norm_one() );
3131 }
3132 else
3133 {
3134 return( $object->norm_one() >= abs($argument) );
3135 }
3136 }
3137 else
3138 {
3139 croak "Math::MatrixReal $name: wrong argument type";
3140 }
3141}
3142sub _exp {
3143 my ($matrix,$arg,$flag) = @_;
3144 my $new_matrix = $matrix->clone();
3145 my ($rows,$cols) = $matrix->dim();
3146
3147 $new_matrix->_undo_LR();
3148
3149 croak "Math::MatrixReal::exp(): Matrix is not quadratic" unless ($rows == $cols);
3150 croak "Math::MatrixReal::exp(): Only diagonal matrices supported" unless ( $matrix->is_diagonal() );
3151
3152 $new_matrix = $matrix->each_diag( sub { exp(shift) } );
3153 return $new_matrix;
3154
3155}
3156sub _cos {
3157 my ($matrix,$arg,$flag) = @_;
3158 my $new_matrix = $matrix->clone();
3159 my ($rows,$cols) = $matrix->dim();
3160
3161 $new_matrix->_undo_LR();
3162
3163 croak "Math::MatrixReal::cos(): Matrix is not quadratic" unless ($rows == $cols);
3164 croak "Math::MatrixReal::cos(): Only diagonal matrices supported" unless ( $matrix->is_diagonal() );
3165
3166 $new_matrix = $matrix->each_diag( sub { cos(shift) } );
3167 return $new_matrix;
3168}
3169sub _sin {
3170 my ($matrix,$arg,$flag) = @_;
3171 my $new_matrix = $matrix->clone();
3172 my ($rows,$cols) = $matrix->dim();
3173
3174 $new_matrix->_undo_LR();
3175
3176 croak "Math::MatrixReal::sin(): Matrix is not quadratic" unless ($rows == $cols);
3177 croak "Math::MatrixReal::sin(): Only diagonal matrices supported" unless ( $matrix->is_diagonal() );
3178
3179 $new_matrix = $matrix->each_diag( sub { sin(shift) } );
3180 return $new_matrix;
3181
3182}
3183
3184sub _clone
3185{
3186 my($object,$argument,$flag) = @_;
3187# my($name) = "'='";
3188 my($temp);
3189
3190 $temp = $object->new($object->[1],$object->[2]);
3191 $temp->copy($object);
3192 $temp->_undo_LR();
3193 return($temp);
3194}
3195
3196sub _trace
3197{
3198 my($text,$object,$argument,$flag) = @_;
3199
3200 unless (defined $object) { $object = 'undef'; };
3201 unless (defined $argument) { $argument = 'undef'; };
3202 unless (defined $flag) { $flag = 'undef'; };
3203 if (ref($object)) { $object = ref($object); }
3204 if (ref($argument)) { $argument = ref($argument); }
3205 print "$text: \$obj='$object' \$arg='$argument' \$flag='$flag'\n";
3206}
3207
32081;
3209
3210__END__
3211
3212=head1 NAME
3213
3214Math::MatrixReal - Matrix of Reals
3215
3216Implements the data type "matrix of reals" (and consequently also
3217"vector of reals").
3218
3219=head1 DESCRIPTION
3220
3221Implements the data type "matrix of reals", which can be used almost
3222like any other basic Perl type thanks to B<OPERATOR OVERLOADING>, i.e.,
3223
3224 $product = $matrix1 * $matrix2;
3225
3226does what you would like it to do (a matrix multiplication).
3227
3228Also features many important operations and methods: matrix norm,
3229matrix transposition, matrix inverse, determinant of a matrix, order
3230and numerical condition of a matrix, scalar product of vectors, vector
3231product of vectors, vector length, projection of row and column vectors,
3232a comfortable way for reading in a matrix from a file, the keyboard or
3233your code, and many more.
3234
3235Allows to solve linear equation systems using an efficient algorithm
3236known as "L-R-decomposition" and several approximative (iterative) methods.
3237
3238Features an implementation of Kleene's algorithm to compute the minimal
3239costs for all paths in a graph with weighted edges (the "weights" being
3240the costs associated with each edge).
3241
3242=head1 SYNOPSIS
3243
3244=head2 Constructor Methods And Such
3245
3246=item *
3247
3248C<use Math::MatrixReal;>
3249
3250Makes the methods and overloaded operators of this module available
3251to your program.
3252
3253=item *
3254
3255C<$new_matrix = new Math::MatrixReal($rows,$columns);>
3256
3257The matrix object constructor method. A new matrix of size $rows by $columns
3258will be created, with the value C<0.0> for all elements.
3259
3260Note that this method is implicitly called by many of the other methods
3261in this module.
3262
3263=item *
3264
3265C<$new_matrix = $some_matrix-E<gt>>C<new($rows,$columns);>
3266
3267Another way of calling the matrix object constructor method.
3268
3269Matrix "C<$some_matrix>" is not changed by this in any way.
3270
3271=item *
3272
3273C<$new_matrix = $matrix-E<gt>new_from_cols( [ $column_vector|$array_ref|$string, ... ] )>
3274
3275Creates a new matrix given a reference to an array of any of the following:
3276
3277=over 4
3278
3279=item * column vectors ( n by 1 Math::MatrixReal matrices )
3280
3281=item * references to arrays
3282
3283=item * strings properly formatted to create a column with Math::MatrixReal's
3284C<new_from_string> command
3285
3286=back
3287
3288You may mix and match these as you wish. However, all must be of the
3289same dimension--no padding happens automatically. Example:
3290
3291 my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
3292 print $matrix;
3293
3294will print
3295
3296 [ 1.000000000000E+00 3.000000000000E+00 ]
3297 [ 2.000000000000E+00 4.000000000000E+00 ]
3298
3299
3300=item *
3301
3302C<new_from_rows( [ $row_vector|$array_ref|$string, ... ] )>
3303
3304Creates a new matrix given a reference to an array of any of the following:
3305
3306=over 4
3307
3308=item * row vectors ( 1 by n Math::MatrixReal matrices )
3309
3310=item * references to arrays
3311
3312=item * strings properly formatted to create a row with Math::MatrixReal's C<new_from_string> command
3313
3314=back
3315
3316You may mix and match these as you wish. However, all must be of the
3317same dimension--no padding happens automatically. Example:
3318
3319 my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
3320 print $matrix;
3321
3322will print
3323
3324 [ 1.000000000000E+00 2.000000000000E+00 ]
3325 [ 3.000000000000E+00 4.000000000000E+00 ]
3326
3327
3328=item *
3329
3330C<$new_matrix = Math::MatrixReal-E<gt>new_diag( $array_ref );>
3331
3332This method allows you to create a diagonal matrix by only specifying
3333the diagonal elements. Example:
3334
3335 $matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
3336 print $matrix;
3337
3338will print
3339
3340 [ 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
3341 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
3342 [ 0.000000000000E+00 0.000000000000E+00 3.000000000000E+00 0.000000000000E+00 ]
3343 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 4.000000000000E+00 ]
3344
3345
3346=item *
3347
3348C<$new_matrix = Math::MatrixReal-E<gt>>C<new_from_string($string);>
3349
3350This method allows you to read in a matrix from a string (for
3351instance, from the keyboard, from a file or from your code).
3352
3353The syntax is simple: each row must start with "C<[ >" and end with
3354"C< ]\n>" ("C<\n>" being the newline character and "C< >" a space or
3355tab) and contain one or more numbers, all separated from each other
3356by spaces or tabs.
3357
3358Additional spaces or tabs can be added at will, but no comments.
3359
3360Examples:
3361
3362 $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
3363 $matrix = Math::MatrixReal->new_from_string($string);
3364 print "$matrix";
3365
3366By the way, this prints
3367
3368 [ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
3369 [ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
3370 [ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
3371
3372But you can also do this in a much more comfortable way using the
3373shell-like "here-document" syntax:
3374
3375 $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
3376 [ 1 0 0 0 0 0 1 ]
3377 [ 0 1 0 0 0 0 0 ]
3378 [ 0 0 1 0 0 0 0 ]
3379 [ 0 0 0 1 0 0 0 ]
3380 [ 0 0 0 0 1 0 0 ]
3381 [ 0 0 0 0 0 1 0 ]
3382 [ 1 0 0 0 0 0 -1 ]
3383 MATRIX
3384
3385You can even use variables in the matrix:
3386
3387 $c1 = 2 / 3;
3388 $c2 = -2 / 5;
3389 $c3 = 26 / 9;
3390
3391 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
3392
3393 [ 3 2 0 ]
3394 [ 0 3 2 ]
3395 [ $c1 $c2 $c3 ]
3396
3397 MATRIX
3398
3399(Remember that you may use spaces and tabs to format the matrix to
3400your taste)
3401
3402Note that this method uses exactly the same representation for a
3403matrix as the "stringify" operator "": this means that you can convert
3404any matrix into a string with C<$string = "$matrix";> and read it back
3405in later (for instance from a file!).
3406
3407Note however that you may suffer a precision loss in this process
3408because only 13 digits are supported in the mantissa when printed!!
3409
3410If the string you supply (or someone else supplies) does not obey
3411the syntax mentioned above, an exception is raised, which can be
3412caught by "eval" as follows:
3413
3414 print "Please enter your matrix (in one line): ";
3415 $string = <STDIN>;
3416 $string =~ s/\\n/\n/g;
3417 eval { $matrix = Math::MatrixReal->new_from_string($string); };
3418 if ($@)
3419 {
3420 print "$@";
3421 # ...
3422 # (error handling)
3423 }
3424 else
3425 {
3426 # continue...
3427 }
3428
3429or as follows:
3430
3431 eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
3432 [ 3 2 0 ]
3433 [ 0 3 2 ]
3434 [ $c1 $c2 $c3 ]
3435 MATRIX
3436 if ($@)
3437 # ...
3438
3439Actually, the method shown above for reading a matrix from the keyboard
3440is a little awkward, since you have to enter a lot of "\n"'s for the
3441newlines.
3442
3443A better way is shown in this piece of code:
3444
3445 while (1)
3446 {
3447 print "\nPlease enter your matrix ";
3448 print "(multiple lines, <ctrl-D> = done):\n";
3449 eval { $new_matrix =
3450 Math::MatrixReal->new_from_string(join('',<STDIN>)); };
3451 if ($@)
3452 {
3453 $@ =~ s/\s+at\b.*?$//;
3454 print "${@}Please try again.\n";
3455 }
3456 else { last; }
3457 }
3458
3459Possible error messages of the "new_from_string()" method are:
3460
3461 Math::MatrixReal::new_from_string(): syntax error in input string
3462 Math::MatrixReal::new_from_string(): empty input string
3463
3464If the input string has rows with varying numbers of columns,
3465the following warning will be printed to STDERR:
3466
3467 Math::MatrixReal::new_from_string(): missing elements will be set to zero!
3468
3469If everything is okay, the method returns an object reference to the
3470(newly allocated) matrix containing the elements you specified.
3471
3472=item *
3473
3474C<$new_matrix = $some_matrix-E<gt>shadow();>
3475
3476Returns an object reference to a B<NEW> but B<EMPTY> matrix
3477(filled with zero's) of the B<SAME SIZE> as matrix "C<$some_matrix>".
3478
3479Matrix "C<$some_matrix>" is not changed by this in any way.
3480
3481=item *
3482
3483C<$matrix1-E<gt>copy($matrix2);>
3484
3485Copies the contents of matrix "C<$matrix2>" to an B<ALREADY EXISTING>
3486matrix "C<$matrix1>" (which must have the same size as matrix "C<$matrix2>"!).
3487
3488Matrix "C<$matrix2>" is not changed by this in any way.
3489
3490=item *
3491
3492C<$twin_matrix = $some_matrix-E<gt>clone();>
3493
3494Returns an object reference to a B<NEW> matrix of the B<SAME SIZE> as
3495matrix "C<$some_matrix>". The contents of matrix "C<$some_matrix>" have
3496B<ALREADY BEEN COPIED> to the new matrix "C<$twin_matrix>". This
3497is the method that the operator "=" is overloaded to when you type
3498C<$a = $b>, when C<$a> and C<$b> are matrices.
3499
3500Matrix "C<$some_matrix>" is not changed by this in any way.
3501
3502
3503=head2 Matrix Row, Column and Element operations
3504
3505=item *
3506
3507C<$row_vector = $matrix-E<gt>row($row);>
3508
3509This is a projection method which returns an object reference to
3510a B<NEW> matrix (which in fact is a (row) vector since it has only
3511one row) to which row number "C<$row>" of matrix "C<$matrix>" has
3512already been copied.
3513
3514Matrix "C<$matrix>" is not changed by this in any way.
3515
3516=item *
3517
3518C<$column_vector = $matrix-E<gt>column($column);>
3519
3520This is a projection method which returns an object reference to
3521a B<NEW> matrix (which in fact is a (column) vector since it has
3522only one column) to which column number "C<$column>" of matrix
3523"C<$matrix>" has already been copied.
3524
3525Matrix "C<$matrix>" is not changed by this in any way.
3526
3527=item *
3528
3529C<$matrix-E<gt>assign($row,$column,$value);>
3530
3531Explicitly assigns a value "C<$value>" to a single element of the
3532matrix "C<$matrix>", located in row "C<$row>" and column "C<$column>",
3533thereby replacing the value previously stored there.
3534
3535=item *
3536
3537C<$value = $matrix-E<gt>>C<element($row,$column);>
3538
3539Returns the value of a specific element of the matrix "C<$matrix>",
3540located in row "C<$row>" and column "C<$column>".
3541
3542=item *
3543
3544C<$new_matrix = $matrix-E<gt>each( \&function )>;
3545
3546Creates a new matrix by evaluating a code reference on each element of the
3547given matrix. The function is passed the element, the row index and the column
3548index, in that order. The value the function returns ( or the value of the last
3549executed statement ) is the value given to the corresponding element in $new_matrix.
3550
3551Example:
3552
3553 # add 1 to every element in the matrix
3554 $matrix = $matrix->each ( sub { (shift) + 1 } );
3555
3556
3557Example:
3558
3559 my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
3560 ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
3561 : -1*$matrix->minor($i,$j)->det();
3562 } );
3563
3564This code needs some explanation. For each element of $matrix, it throws away the actual value
3565and stores the row and column indexes in $i and $j. Then it sets element [$i,$j] in $cofactor
3566to the determinant of C<$matrix-E<gt>minor($i,$j)> if it is an "even" element, or C<-1*$matrix-E<gt>minor($i,$j)>
3567if it is an "odd" element.
3568
3569=item *
3570
3571C<$new_matrix = $matrix-E<gt>each_diag( \&function )>;
3572
3573Creates a new matrix by evaluating a code reference on each diagonal element of the
3574given matrix. The function is passed the element, the row index and the column
3575index, in that order. The value the function returns ( or the value of the last
3576executed statement ) is the value given to the corresponding element in $new_matrix.
3577
3578
3579=item *
3580
3581C<$matrix-E<gt>swap_col( $col1, $col2 );>
3582
3583This method takes two one-based column numbers and swaps the values of each element in each column.
3584C<$matrix-E<gt>swap_col(2,3)> would replace column 2 in $matrix with column 3, and replace column
35853 with column 2.
3586
3587=item *
3588
3589C<$matrix-E<gt>swap_row( $row1, $row2 );>
3590
3591This method takes two one-based row numbers and swaps the values of each element in each row.
3592C<$matrix-E<gt>swap_row(2,3)> would replace row 2 in $matrix with row 3, and replace row
35933 with row 2.
3594
3595
3596=head2 Matrix Operations
3597
3598=item *
3599
3600C<$det = $matrix-E<gt>det();>
3601
3602Returns the determinant of the matrix, without going through
3603the rigamarole of computing a LR decomposition. This method should
3604be much faster than LR decomposition if the matrix is diagonal or
3605triangular. Otherwise, it is just a wrapper for
3606C<$matrix-E<gt>decompose_LR-E<gt>det_LR>. If the determinant is zero,
3607there is no inverse and vice-versa. Only quadratic matrices have
3608determinants.
3609
3610=item *
3611
3612C<$inverse = $matrix-E<gt>inverse();>
3613
3614Returns the inverse of a matrix, without going through the
3615rigamarole of computing a LR decomposition. If no inverse exists,
3616undef is returned and an error is printed via C<carp()>.
3617This is nothing but a wrapper for C<$matrix-E<gt>decompose_LR-E<gt>invert_LR>.
3618
3619=item *
3620
3621C<($rows,$columns) = $matrix-E<gt>dim();>
3622
3623Returns a list of two items, representing the number of rows
3624and columns the given matrix "C<$matrix>" contains.
3625
3626=item *
3627
3628C<$norm_one = $matrix-E<gt>norm_one();>
3629
3630Returns the "one"-norm of the given matrix "C<$matrix>".
3631
3632The "one"-norm is defined as follows:
3633
3634For each column, the sum of the absolute values of the elements in the
3635different rows of that column is calculated. Finally, the maximum
3636of these sums is returned.
3637
3638Note that the "one"-norm and the "maximum"-norm are mathematically
3639equivalent, although for the same matrix they usually yield a different
3640value.
3641
3642Therefore, you should only compare values that have been calculated
3643using the same norm!
3644
3645Throughout this package, the "one"-norm is (arbitrarily) used
3646for all comparisons, for the sake of uniformity and comparability,
3647except for the iterative methods "solve_GSM()", "solve_SSM()" and
3648"solve_RM()" which use either norm depending on the matrix itself.
3649
3650=item *
3651
3652C<$norm_max = $matrix-E<gt>norm_max();>
3653
3654Returns the "maximum"-norm of the given matrix $matrix.
3655
3656The "maximum"-norm is defined as follows:
3657
3658For each row, the sum of the absolute values of the elements in the
3659different columns of that row is calculated. Finally, the maximum
3660of these sums is returned.
3661
3662Note that the "maximum"-norm and the "one"-norm are mathematically
3663equivalent, although for the same matrix they usually yield a different
3664value.
3665
3666Therefore, you should only compare values that have been calculated
3667using the same norm!
3668
3669Throughout this package, the "one"-norm is (arbitrarily) used
3670for all comparisons, for the sake of uniformity and comparability,
3671except for the iterative methods "solve_GSM()", "solve_SSM()" and
3672"solve_RM()" which use either norm depending on the matrix itself.
3673
3674=item *
3675
3676C<$norm_sum = $matrix-E<gt>norm_sum();>
3677
3678This is a very simple norm which is defined as the sum of the
3679absolute values of every element.
3680
3681=item *
3682
3683C<$p_norm> = $matrix-E<gt>norm_p($n);>
3684
3685This function returns the "p-norm" of a vector. The argument $n
3686must be a number greater than or equal to 1 or the string "Inf".
3687The p-norm is defined as (sum(x_i^p))^(1/p). In words, it raised
3688each element to the p-th power, adds them up, and then takes the
3689p-th root of that number. If the string "Inf" is passed, the
3690"infinity-norm" is computed, which is really the limit of the
3691p-norm as p goes to infinity. It is defined as the maximum element
3692of the vector. Also, note that the familiar Euclidean distance
3693between two vectors is just a special case of a p-norm, when p is
3694equal to 2.
3695
3696Example:
3697 $a = Math::MatrixReal->new_from_cols([[1,2,3]]);
3698 $p1 = $a->norm_p(1);
3699 $p2 = $a->norm_p(2);
3700 $p3 = $a->norm_p(3);
3701 $pinf = $a->norm_p("Inf");
3702
3703 print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";
3704
3705 $i1 = $a->new_from_rows([[1,0]]);
3706 $i2 = $a->new_from_rows([[0,1]]);
3707
3708 # this should be sqrt(2) since it is the same as the
3709 # hypotenuse of a 1 by 1 right triangle
3710
3711 $dist = ($i1-$i2)->norm_p(2);
3712 print "Distance is $dist, which should be " . sqrt(2) . "\n";
3713
3714Output:
3715
3716 (1,2,3,Inf) norm:
3717 6
3718 3.74165738677394139
3719 3.30192724889462668
3720 3
3721
3722 Distance is 1.41421356237309505, which should be 1.41421356237309505
3723
3724
3725
3726=item *
3727
3728C<$frob_norm> = C<$matrix-E<gt>norm_frobenius();>
3729
3730This norm is similar to that of a p-norm where p is 2, except it
3731acts on a B<matrix>, not a vector. Each element of the matrix is
3732squared, this is added up, and then a square root is taken.
3733
3734=item *
3735
3736C<$matrix-E<gt>spectral_radius();>
3737
3738Returns the maximum value of the absolute value of all eigenvalues.
3739Currently this computes B<all> eigenvalues, then sifts through them
3740to find the largest in absolute value. Needless to say, this is very
3741inefficient, and in the future an algorithm that computes only the
3742largest eigenvalue may be implemented.
3743
3744=item *
3745
3746C<$matrix1-E<gt>transpose($matrix2);>
3747
3748Calculates the transposed matrix of matrix $matrix2 and stores
3749the result in matrix "C<$matrix1>" (which must already exist and have
3750the same size as matrix "C<$matrix2>"!).
3751
3752This operation can also be carried out "in-place", i.e., input and
3753output matrix may be identical.
3754
3755Transposition is a symmetry operation: imagine you rotate the matrix
3756along the axis of its main diagonal (going through elements (1,1),
3757(2,2), (3,3) and so on) by 180 degrees.
3758
3759Another way of looking at it is to say that rows and columns are
3760swapped. In fact the contents of element C<(i,j)> are swapped
3761with those of element C<(j,i)>.
3762
3763Note that (especially for vectors) it makes a big difference if you
3764have a row vector, like this:
3765
3766 [ -1 0 1 ]
3767
3768or a column vector, like this:
3769
3770 [ -1 ]
3771 [ 0 ]
3772 [ 1 ]
3773
3774the one vector being the transposed of the other!
3775
3776This is especially true for the matrix product of two vectors:
3777
3778 [ -1 ]
3779 [ -1 0 1 ] * [ 0 ] = [ 2 ] , whereas
3780 [ 1 ]
3781
3782 * [ -1 0 1 ]
3783 [ -1 ] [ 1 0 -1 ]
3784 [ 0 ] * [ -1 0 1 ] = [ -1 ] [ 1 0 -1 ] = [ 0 0 0 ]
3785 [ 1 ] [ 0 ] [ 0 0 0 ] [ -1 0 1 ]
3786 [ 1 ] [ -1 0 1 ]
3787
3788So be careful about what you really mean!
3789
3790Hint: throughout this module, whenever a vector is explicitly required
3791for input, a B<COLUMN> vector is expected!
3792
3793=item *
3794
3795C<$trace = $matrix-E<gt>trace();>
3796
3797This returns the trace of the matrix, which is defined as
3798the sum of the diagonal elements. The matrix must be
3799quadratic.
3800
3801=item *
3802
3803C<$minor = $matrix-E<gt>minor($row,$col);>
3804
3805Returns the minor matrix corresponding to $row and $col. $matrix must be quadratic.
3806If $matrix is n rows by n cols, the minor of $row and $col will be an (n-1) by (n-1)
3807matrix. The minor is defined as crossing out the row and the col specified and returning
3808the remaining rows and columns as a matrix. This method is used by C<cofactor()>.
3809
3810=item *
3811
3812C<$cofactor = $matrix-E<gt>cofactor();>
3813
3814The cofactor matrix is constructed as follows:
3815
3816For each element, cross out the row and column that it sits in.
3817Now, take the determinant of the matrix that is left in the other
3818rows and columns.
3819Multiply the determinant by (-1)^(i+j), where i is the row index,
3820and j is the column index.
3821Replace the given element with this value.
3822
3823The cofactor matrix can be used to find the inverse of the matrix. One formula for the
3824inverse of a matrix is the cofactor matrix transposed divided by the original
3825determinant of the matrix.
3826
3827The following two inverses should be exactly the same:
3828
3829 my $inverse1 = $matrix->inverse;
3830 my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
3831
3832Caveat: Although the cofactor matrix is simple algorithm to compute the inverse of a matrix, and
3833can be used with pencil and paper for small matrices, it is comically slower than
3834the native C<inverse()> function. Here is a small benchmark:
3835
3836 # $matrix1 is 15x15
3837 $det = $matrix1->det;
3838 timethese( 10,
3839 {'inverse' => sub { $matrix1->inverse(); },
3840 'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
3841 } );
3842
3843
3844 Benchmark: timing 10 iterations of LR, cofactor, inverse...
3845 inverse: 1 wallclock secs ( 0.56 usr + 0.00 sys = 0.56 CPU) @ 17.86/s (n=10)
3846 cofactor: 36 wallclock secs (36.62 usr + 0.01 sys = 36.63 CPU) @ 0.27/s (n=10)
3847
3848=item *
3849
3850C<$adjoint = $matrix-E<gt>adjoint();>
3851
3852The adjoint is just the transpose of the cofactor matrix. This method is
3853just an alias for C< ~($matrix-E<gt>cofactor)>.
3854
3855
3856=head2 Arithmetic Operations
3857
3858=item *
3859
3860C<$matrix1-E<gt>add($matrix2,$matrix3);>
3861
3862Calculates the sum of matrix "C<$matrix2>" and matrix "C<$matrix3>"
3863and stores the result in matrix "C<$matrix1>" (which must already exist
3864and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
3865
3866This operation can also be carried out "in-place", i.e., the output and
3867one (or both) of the input matrices may be identical.
3868
3869=item *
3870
3871C<$matrix1-E<gt>subtract($matrix2,$matrix3);>
3872
3873Calculates the difference of matrix "C<$matrix2>" minus matrix "C<$matrix3>"
3874and stores the result in matrix "C<$matrix1>" (which must already exist
3875and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
3876
3877This operation can also be carried out "in-place", i.e., the output and
3878one (or both) of the input matrices may be identical.
3879
3880Note that this operation is the same as
3881C<$matrix1-E<gt>add($matrix2,-$matrix3);>, although the latter is
3882a little less efficient.
3883
3884=item *
3885
3886C<$matrix1-E<gt>multiply_scalar($matrix2,$scalar);>
3887
3888Calculates the product of matrix "C<$matrix2>" and the number "C<$scalar>"
3889(i.e., multiplies each element of matrix "C<$matrix2>" with the factor
3890"C<$scalar>") and stores the result in matrix "C<$matrix1>" (which must
3891already exist and have the same size as matrix "C<$matrix2>"!).
3892
3893This operation can also be carried out "in-place", i.e., input and
3894output matrix may be identical.
3895
3896=item *
3897
3898C<$product_matrix = $matrix1-E<gt>multiply($matrix2);>
3899
3900Calculates the product of matrix "C<$matrix1>" and matrix "C<$matrix2>"
3901and returns an object reference to a new matrix "C<$product_matrix>" in
3902which the result of this operation has been stored.
3903
3904Note that the dimensions of the two matrices "C<$matrix1>" and "C<$matrix2>"
3905(i.e., their numbers of rows and columns) must harmonize in the following
3906way (example):
3907
3908 [ 2 2 ]
3909 [ 2 2 ]
3910 [ 2 2 ]
3911
3912 [ 1 1 1 ] [ * * ]
3913 [ 1 1 1 ] [ * * ]
3914 [ 1 1 1 ] [ * * ]
3915 [ 1 1 1 ] [ * * ]
3916
3917I.e., the number of columns of matrix "C<$matrix1>" has to be the same
3918as the number of rows of matrix "C<$matrix2>".
3919
3920The number of rows and columns of the resulting matrix "C<$product_matrix>"
3921is determined by the number of rows of matrix "C<$matrix1>" and the number
3922of columns of matrix "C<$matrix2>", respectively.
3923
3924=item *
3925
3926C<$matrix1-E<gt>negate($matrix2);>
3927
3928Calculates the negative of matrix "C<$matrix2>" (i.e., multiplies
3929all elements with "-1") and stores the result in matrix "C<$matrix1>"
3930(which must already exist and have the same size as matrix "C<$matrix2>"!).
3931
3932This operation can also be carried out "in-place", i.e., input and
3933output matrix may be identical.
3934
3935
3936=item *
3937
3938C<$matrix_to_power = $matrix1-E<gt>exponent($integer);>
3939
3940Raises the matrix to the C<$integer> power. Obviously, C<$integer> must
3941be an integer. If it is zero, the identity matrix is returned. If a negative
3942integer is given, the inverse will be computed (if it exists) and then raised
3943the the absolute value of C<$integer>. The matrix must be quadratic.
3944
3945
3946=head2 Boolean Matrix Operations
3947
3948=item *
3949
3950C<$matrix-E<gt>is_quadratic();>
3951
3952Returns a boolean value indicating if the given matrix is
3953quadratic (also know as "square" or "n by n"). A matrix is
3954quadratic if it has the same number of rows as it does columns.
3955
3956=item *
3957
3958C<$matrix-E<gt>is_square();>
3959
3960This is an alias for C<is_quadratic()>.
3961
3962
3963=item *
3964
3965C<$matrix-E<gt>is_symmetric();>
3966
3967Returns a boolean value indicating if the given matrix is
3968symmetric. By definition, a matrix is symmetric if and only
3969if (B<M>[I<i>,I<j>]=B<M>[I<j>,I<i>]). This is equivalent to
3970C<($matrix == ~$matrix)> but without memory allocation.
3971Only quadratic matrices can be symmetric.
3972
3973Notes: A symmetric matrix always has real eigenvalues/eigenvectors.
3974A matrix plus its transpose is always symmetric.
3975
3976=item *
3977
3978C<$matrix-E<gt>is_skew_symmetric();>
3979
3980Returns a boolean value indicating if the given matrix is
3981skew symmetric. By definition, a matrix is symmetric if and only
3982if (B<M>[I<i>,I<j>]=B<-M>[I<j>,I<i>]). This is equivalent to
3983C<($matrix == -(~$matrix))> but without memory allocation.
3984Only quadratic matrices can be skew symmetric.
3985
3986
3987=item *
3988
3989C<$matrix-E<gt>is_diagonal();>
3990
3991Returns a boolean value indicating if the given matrix is
3992diagonal, i.e. all of the nonzero elements are on the main diagonal.
3993Only quadratic matrices can be diagonal.
3994
3995=item *
3996
3997C<$matrix-E<gt>is_tridiagonal();>
3998
3999Returns a boolean value indicating if the given matrix is
4000tridiagonal, i.e. all of the nonzero elements are on the main diagonal
4001or the diagonals above and below the main diagonal.
4002Only quadratic matrices can be tridiagonal.
4003
4004
4005=item *
4006
4007C<$matrix-E<gt>is_upper_triangular();>
4008
4009Returns a boolean value indicating if the given matrix is upper triangular,
4010i.e. all of the nonzero elements not on the main diagonal are above it.
4011Only quadratic matrices can be upper triangular.
4012Note: diagonal matrices are both upper and lower triangular.
4013
4014
4015=item *
4016
4017C<$matrix-E<gt>is_lower_triangular();>
4018
4019Returns a boolean value indicating if the given matrix is lower triangular,
4020i.e. all of the nonzero elements not on the main diagonal are below it.
4021Only quadratic matrices can be lower triangular.
4022Note: diagonal matrices are both upper and lower triangular.
4023
4024=item *
4025
4026C<$matrix-E<gt>is_orthogonal();>
4027
4028Returns a boolean value indicating if the given matrix is orthogonal.
4029An orthogonal matrix is has the property that the transpose equals the
4030inverse of the matrix. Instead of computing each and comparing them, this
4031method multiplies the matrix by it's transpose, and returns true if this
4032turns out to be the identity matrix, false otherwise.
4033Only quadratic matrices can orthogonal.
4034
4035=item *
4036
4037C<$matrix-E<gt>is_binary();>
4038
4039Returns a boolean value indicating if the given matrix is binary.
4040A matrix is binary if it contains only zeroes or ones.
4041
4042=item *
4043
4044C<$matrix-E<gt>is_gramian();>
4045
4046Returns a boolean value indicating if the give matrix is Gramian.
4047A matrix C<$A> is Gramian if and only if there exists a
4048square matrix C<$B> such that C<$A = ~$B*$B>. This is equivalent to
4049checking if C<$A> is symmetric and has all nonnegative eigenvalues, which
4050is what Math::MatrixReal uses to check for this property.
4051
4052=item *
4053
4054C<$matrix-E<gt>is_LR();>
4055
4056Returns a boolean value indicating if the matrix is an LR decomposition
4057matrix.
4058
4059=item *
4060
4061C<$matrix-E<gt>is_positive();>
4062
4063Returns a boolean value indicating if the matrix contains only
4064positive entries. Note that a zero entry is not positive and
4065will cause C<is_positive()> to return false.
4066
4067=item *
4068
4069C<$matrix-E<gt>is_negative();>
4070
4071Returns a boolean value indicating if the matrix contains only
4072negative entries. Note that a zero entry is not negative and
4073will cause C<is_negative()> to return false.
4074
4075=item *
4076
4077C<$matrix-E<gt>is_periodic($k);>
4078
4079Returns a boolean value indicating if the matrix is periodic
4080with period $k. This is true if C<$matrix ** ($k+1) == $matrix>.
4081When C<$k == 1>, this reduces down to the C<is_idempotent()>
4082function.
4083
4084=item *
4085
4086C<$matrix-E<gt>is_idempotent();>
4087
4088Returns a boolean value indicating if the matrix is idempotent,
4089which is defined as the square of the matrix being equal to
4090the original matrix, i.e C<$matrix ** 2 == $matrix>.
4091
4092=item *
4093
4094C<$matrix-E<gt>is_row_vector();>
4095
4096Returns a boolean value indicating if the matrix is a row vector.
4097A row vector is a matrix which is 1xn. Note that the 1x1 matrix is
4098both a row and column vector.
4099
4100=item *
4101
4102C<$matrix-E<gt>is_col_vector();>
4103
4104Returns a boolean value indicating if the matrix is a col vector.
4105A col vector is a matrix which is nx1. Note that the 1x1 matrix is
4106both a row and column vector.
4107
4108=head2 Eigensystems
4109
4110=over 2
4111
4112=item *
4113
4114C<($l, $V) = $matrix-E<gt>sym_diagonalize();>
4115
4116This method performs the diagonalization of the quadratic
4117I<symmetric> matrix B<M> stored in $matrix.
4118On output, B<l> is a column vector containing all the eigenvalues
4119of B<M> and B<V> is an orthogonal matrix which columns are the
4120corresponding normalized eigenvectors.
4121The primary property of an eigenvalue I<l> and an eigenvector
4122B<x> is of course that: B<M> * B<x> = I<l> * B<x>.
4123
4124The method uses a Householder reduction to tridiagonal form
4125followed by a QL algoritm with implicit shifts on this
4126tridiagonal. (The tridiagonal matrix is kept internally
4127in a compact form in this routine to save memory.)
4128In fact, this routine wraps the householder() and
4129tri_diagonalize() methods described below when their
4130intermediate results are not desired.
4131The overall algorithmic complexity of this technique
4132is O(N^3). According to several books, the coefficient
4133hidden by the 'O' is one of the best possible for general
4134(symmetric) matrixes.
4135
4136=item *
4137
4138C<($T, $Q) = $matrix-E<gt>householder();>
4139
4140This method performs the Householder algorithm which reduces
4141the I<n> by I<n> real I<symmetric> matrix B<M> contained
4142in $matrix to tridiagonal form.
4143On output, B<T> is a symmetric tridiagonal matrix (only
4144diagonal and off-diagonal elements are non-zero) and B<Q>
4145is an I<orthogonal> matrix performing the tranformation
4146between B<M> and B<T> (C<$M == $Q * $T * ~$Q>).
4147
4148=item *
4149
4150C<($l, $V) = $T-E<gt>tri_diagonalize([$Q]);>
4151
4152This method diagonalizes the symmetric tridiagonal
4153matrix B<T>. On output, $l and $V are similar to the
4154output values described for sym_diagonalize().
4155
4156The optional argument $Q corresponds to an orthogonal
4157transformation matrix B<Q> that should be used additionally
4158during B<V> (eigenvectors) computation. It should be supplied
4159if the desired eigenvectors correspond to a more general
4160symmetric matrix B<M> previously reduced by the
4161householder() method, not a mere tridiagonal. If B<T> is
4162really a tridiagonal matrix, B<Q> can be omitted (it
4163will be internally created in fact as an identity matrix).
4164The method uses a QL algorithm (with implicit shifts).
4165
4166=item *
4167
4168C<$l = $matrix-E<gt>sym_eigenvalues();>
4169
4170This method computes the eigenvalues of the quadratic
4171I<symmetric> matrix B<M> stored in $matrix.
4172On output, B<l> is a column vector containing all the eigenvalues
4173of B<M>. Eigenvectors are not computed (on the contrary of
4174C<sym_diagonalize()>) and this method is more efficient
4175(even though it uses a similar algorithm with two phases).
4176However, understand that the algorithmic complexity of this
4177technique is still also O(N^3). But the coefficient hidden
4178by the 'O' is better by a factor of..., well, see your
4179benchmark, it's wiser.
4180
4181This routine wraps the householder_tridiagonal() and
4182tri_eigenvalues() methods described below when the
4183intermediate tridiagonal matrix is not needed.
4184
4185=item *
4186
4187C<$T = $matrix-E<gt>householder_tridiagonal();>
4188
4189This method performs the Householder algorithm which reduces
4190the I<n> by I<n> real I<symmetric> matrix B<M> contained
4191in $matrix to tridiagonal form.
4192On output, B<T> is the obtained symmetric tridiagonal matrix
4193(only diagonal and off-diagonal elements are non-zero). The
4194operation is similar to the householder() method, but potentially
4195a little more efficient as the transformation matrix is not
4196computed.
4197
4198=item *
4199
4200C<$l = $T-E<gt>tri_eigenvalues();>
4201
4202This method computesthe eigenvalues of the symmetric
4203tridiagonal matrix B<T>. On output, $l is a vector
4204containing the eigenvalues (similar to C<sym_eigenvalues()>).
4205This method is much more efficient than tri_diagonalize()
4206when eigenvectors are not needed.
4207
4208=back
4209
4210=head2 Miscellaneous
4211
4212=item *
4213
4214C<$matrix-E<gt>zero();>
4215
4216Assigns a zero to every element of the matrix "C<$matrix>", i.e.,
4217erases all values previously stored there, thereby effectively
4218transforming the matrix into a "zero"-matrix or "null"-matrix,
4219the neutral element of the addition operation in a Ring.
4220
4221(For instance the (quadratic) matrices with "n" rows and columns
4222and matrix addition and multiplication form a Ring. Most prominent
4223characteristic of a Ring is that multiplication is not commutative,
4224i.e., in general, "C<matrix1 * matrix2>" is not the same as
4225"C<matrix2 * matrix1>"!)
4226
4227=item *
4228
4229C<$matrix-E<gt>one();>
4230
4231Assigns one's to the elements on the main diagonal (elements (1,1),
4232(2,2), (3,3) and so on) of matrix "C<$matrix>" and zero's to all others,
4233thereby erasing all values previously stored there and transforming the
4234matrix into a "one"-matrix, the neutral element of the multiplication
4235operation in a Ring.
4236
4237(If the matrix is quadratic (which this method doesn't require, though),
4238then multiplying this matrix with itself yields this same matrix again,
4239and multiplying it with some other matrix leaves that other matrix
4240unchanged!)
4241
4242=item *
4243
4244C<$latex_string = $matrix-E<gt>as_latex( align=E<gt> "c", format =E<gt> "%s", name =E<gt> "" );>
4245
4246This function returns the matrix as a LaTeX string. It takes a hash as an
4247argument which is used to control the style of the output. The hash element C<align>
4248may be "c","l" or "r", corresponding to center, left and right, respectively. The
4249C<format> element is a format string that is given to C<sprintf> to control the
4250style of number format, such a floating point or scientific notation. The C<name>
4251element can be used so that a LaTeX string of "$name = " is prepended to the string.
4252
4253Example:
4254
4255 my $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
4256 print $a->as_latex( ( format => "%.2f", align => "l",name => "A" ) );
4257
4258Output:
4259 $A = $ $
4260 \left( \begin{array}{ll}
4261 1.23&1.00 \\
4262 5.68&2.00 \\
4263 9.10&3.00
4264 \end{array} \right)
4265 $
4266
4267=item *
4268
4269C<$yacas_string = $matrix-E<gt>as_yacas( format =E<gt> "%s", name =E<gt> "", semi =E<gt> 0 );>
4270
4271This function returns the matrix as a string that can be read by Yacas.
4272It takes a hash as
4273an an argument which controls the style of the output. The
4274C<format> element is a format string that is given to C<sprintf> to control the
4275style of number format, such a floating point or scientific notation. The C<name>
4276element can be used so that "$name = " is prepended to the string. The <semi> element can
4277be set to 1 to that a semicolon is appended (so Matlab does not print out the matrix.)
4278
4279Example:
4280
4281 $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
4282 print $a->as_yacas( ( format => "%.2f", align => "l",name => "A" ) );
4283
4284Output:
4285
4286 A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
4287
4288=item *
4289
4290C<$matlab_string = $matrix-E<gt>as_matlab( format =E<gt> "%s", name =E<gt> "", semi =E<gt> 0 );>
4291
4292This function returns the matrix as a string that can be read by Matlab. It takes a hash as
4293an an argument which controls the style of the output. The
4294C<format> element is a format string that is given to C<sprintf> to control the
4295style of number format, such a floating point or scientific notation. The C<name>
4296element can be used so that "$name = " is prepended to the string. The <semi> element can
4297be set to 1 to that a semicolon is appended (so Matlab does not print out the matrix.)
4298
4299Example:
4300
4301 my $a = Math::MatrixReal->new_from_rows([[ 1.234, 5.678, 9.1011],[1,2,3]] );
4302 print $a->as_matlab( ( format => "%.3f", name => "A",semi => 1 ) );
4303
4304Output:
4305 A = [ 1.234 5.678 9.101;
4306 1.000 2.000 3.000];
4307
4308
4309=item *
4310
4311C<$scilab_string = $matrix-E<gt>as_scilab( format =E<gt> "%s", name =E<gt> "", semi =E<gt> 0 );>
4312
4313This function is just an alias for C<as_matlab()>, since both Scilab and Matlab have the
4314same matrix format.
4315
4316=item *
4317
4318C<$minimum = Math::MatrixReal::min($number1,$number2);>
4319
4320Returns the minimum of the two numbers "C<number1>" and "C<number2>".
4321
4322=item *
4323
4324C<$maximum = Math::MatrixReal::max($number1,$number2);>
4325
4326Returns the maximum of the two numbers "C<number1>" and "C<number2>".
4327
4328=item *
4329
4330C<$minimal_cost_matrix = $cost_matrix-E<gt>kleene();>
4331
4332Copies the matrix "C<$cost_matrix>" (which has to be quadratic!) to
4333a new matrix of the same size (i.e., "clones" the input matrix) and
4334applies Kleene's algorithm to it.
4335
4336See L<Math::Kleene(3)> for more details about this algorithm!
4337
4338The method returns an object reference to the new matrix.
4339
4340Matrix "C<$cost_matrix>" is not changed by this method in any way.
4341
4342=item *
4343
4344C<($norm_matrix,$norm_vector) = $matrix-E<gt>normalize($vector);>
4345
4346This method is used to improve the numerical stability when solving
4347linear equation systems.
4348
4349Suppose you have a matrix "A" and a vector "b" and you want to find
4350out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
4351solves the equation system represented by the matrix "A" and the
4352vector "b".
4353
4354Applying this method to the pair (A,b) yields a pair (A',b') where
4355each row has been divided by (the absolute value of) the greatest
4356coefficient appearing in that row. So this coefficient becomes equal
4357to "1" (or "-1") in the new pair (A',b') (all others become smaller
4358than one and greater than minus one).
4359
4360Note that this operation does not change the equation system itself
4361because the same division is carried out on either side of the equation
4362sign!
4363
4364The method requires a quadratic (!) matrix "C<$matrix>" and a vector
4365"C<$vector>" for input (the vector must be a column vector with the same
4366number of rows as the input matrix) and returns a list of two items
4367which are object references to a new matrix and a new vector, in this
4368order.
4369
4370The output matrix and vector are clones of the input matrix and vector
4371to which the operation explained above has been applied.
4372
4373The input matrix and vector are not changed by this in any way.
4374
4375Example of how this method can affect the result of the methods to solve
4376equation systems (explained immediately below following this method):
4377
4378Consider the following little program:
4379
4380 #!perl -w
4381
4382 use Math::MatrixReal qw(new_from_string);
4383
4384 $A = Math::MatrixReal->new_from_string(<<"MATRIX");
4385 [ 1 2 3 ]
4386 [ 5 7 11 ]
4387 [ 23 19 13 ]
4388 MATRIX
4389
4390 $b = Math::MatrixReal->new_from_string(<<"MATRIX");
4391 [ 0 ]
4392 [ 1 ]
4393 [ 29 ]
4394 MATRIX
4395
4396 $LR = $A->decompose_LR();
4397 if (($dim,$x,$B) = $LR->solve_LR($b))
4398 {
4399 $test = $A * $x;
4400 print "x = \n$x";
4401 print "A * x = \n$test";
4402 }
4403
4404 ($A_,$b_) = $A->normalize($b);
4405
4406 $LR = $A_->decompose_LR();
4407 if (($dim,$x,$B) = $LR->solve_LR($b_))
4408 {
4409 $test = $A * $x;
4410 print "x = \n$x";
4411 print "A * x = \n$test";
4412 }
4413
4414This will print:
4415
4416 x =
4417 [ 1.000000000000E+00 ]
4418 [ 1.000000000000E+00 ]
4419 [ -1.000000000000E+00 ]
4420 A * x =
4421 [ 4.440892098501E-16 ]
4422 [ 1.000000000000E+00 ]
4423 [ 2.900000000000E+01 ]
4424 x =
4425 [ 1.000000000000E+00 ]
4426 [ 1.000000000000E+00 ]
4427 [ -1.000000000000E+00 ]
4428 A * x =
4429 [ 0.000000000000E+00 ]
4430 [ 1.000000000000E+00 ]
4431 [ 2.900000000000E+01 ]
4432
4433You can see that in the second example (where "normalize()" has been used),
4434the result is "better", i.e., more accurate!
4435
4436=item *
4437
4438C<$LR_matrix = $matrix-E<gt>decompose_LR();>
4439
4440This method is needed to solve linear equation systems.
4441
4442Suppose you have a matrix "A" and a vector "b" and you want to find
4443out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
4444solves the equation system represented by the matrix "A" and the
4445vector "b".
4446
4447You might also have a matrix "A" and a whole bunch of different
4448vectors "b1".."bk" for which you need to find vectors "x1".."xk"
4449so that C<A * xi = bi>, for C<i=1..k>.
4450
4451Using Gaussian transformations (multiplying a row or column with
4452a factor, swapping two rows or two columns and adding a multiple
4453of one row or column to another), it is possible to decompose any
4454matrix "A" into two triangular matrices, called "L" and "R" (for
4455"Left" and "Right").
4456
4457"L" has one's on the main diagonal (the elements (1,1), (2,2), (3,3)
4458and so so), non-zero values to the left and below of the main diagonal
4459and all zero's in the upper right half of the matrix.
4460
4461"R" has non-zero values on the main diagonal as well as to the right
4462and above of the main diagonal and all zero's in the lower left half
4463of the matrix, as follows:
4464
4465 [ 1 0 0 0 0 ] [ x x x x x ]
4466 [ x 1 0 0 0 ] [ 0 x x x x ]
4467 L = [ x x 1 0 0 ] R = [ 0 0 x x x ]
4468 [ x x x 1 0 ] [ 0 0 0 x x ]
4469 [ x x x x 1 ] [ 0 0 0 0 x ]
4470
4471Note that "C<L * R>" is equivalent to matrix "A" in the sense that
4472C<L * R * x = b E<lt>==E<gt> A * x = b> for all vectors "x", leaving
4473out of account permutations of the rows and columns (these are taken
4474care of "magically" by this module!) and numerical errors.
4475
4476Trick:
4477
4478Because we know that "L" has one's on its main diagonal, we can
4479store both matrices together in the same array without information
4480loss! I.e.,
4481
4482 [ R R R R R ]
4483 [ L R R R R ]
4484 LR = [ L L R R R ]
4485 [ L L L R R ]
4486 [ L L L L R ]
4487
4488Beware, though, that "LR" and "C<L * R>" are not the same!!!
4489
4490Note also that for the same reason, you cannot apply the method "normalize()"
4491to an "LR" decomposition matrix. Trying to do so will yield meaningless
4492rubbish!
4493
4494(You need to apply "normalize()" to each pair (Ai,bi) B<BEFORE> decomposing
4495the matrix "Ai'"!)
4496
4497Now what does all this help us in solving linear equation systems?
4498
4499It helps us because a triangular matrix is the next best thing
4500that can happen to us besides a diagonal matrix (a matrix that
4501has non-zero values only on its main diagonal - in which case
4502the solution is trivial, simply divide "C<b[i]>" by "C<A[i,i]>"
4503to get "C<x[i]>"!).
4504
4505To find the solution to our problem "C<A * x = b>", we divide this
4506problem in parts: instead of solving C<A * x = b> directly, we first
4507decompose "A" into "L" and "R" and then solve "C<L * y = b>" and
4508finally "C<R * x = y>" (motto: divide and rule!).
4509
4510From the illustration above it is clear that solving "C<L * y = b>"
4511and "C<R * x = y>" is straightforward: we immediately know that
4512C<y[1] = b[1]>. We then deduce swiftly that
4513
4514 y[2] = b[2] - L[2,1] * y[1]
4515
4516(and we know "C<y[1]>" by now!), that
4517
4518 y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
4519
4520and so on.
4521
4522Having effortlessly calculated the vector "y", we now proceed to
4523calculate the vector "x" in a similar fashion: we see immediately
4524that C<x[n] = y[n] / R[n,n]>. It follows that
4525
4526 x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
4527
4528and
4529
4530 x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
4531 / R[n-2,n-2]
4532
4533and so on.
4534
4535You can see that - especially when you have many vectors "b1".."bk"
4536for which you are searching solutions to C<A * xi = bi> - this scheme
4537is much more efficient than a straightforward, "brute force" approach.
4538
4539This method requires a quadratic matrix as its input matrix.
4540
4541If you don't have that many equations, fill up with zero's (i.e., do
4542nothing to fill the superfluous rows if it's a "fresh" matrix, i.e.,
4543a matrix that has been created with "new()" or "shadow()").
4544
4545The method returns an object reference to a new matrix containing the
4546matrices "L" and "R".
4547
4548The input matrix is not changed by this method in any way.
4549
4550Note that you can "copy()" or "clone()" the result of this method without
4551losing its "magical" properties (for instance concerning the hidden
4552permutations of its rows and columns).
4553
4554However, as soon as you are applying any method that alters the contents
4555of the matrix, its "magical" properties are stripped off, and the matrix
4556immediately reverts to an "ordinary" matrix (with the values it just happens
4557to contain at that moment, be they meaningful as an ordinary matrix or not!).
4558
4559=item *
4560
4561C<($dimension,$x_vector,$base_matrix) = $LR_matrix>C<-E<gt>>C<solve_LR($b_vector);>
4562
4563Use this method to actually solve an equation system.
4564
4565Matrix "C<$LR_matrix>" must be a (quadratic) matrix returned by the
4566method "decompose_LR()", the LR decomposition matrix of the matrix
4567"A" of your equation system C<A * x = b>.
4568
4569The input vector "C<$b_vector>" is the vector "b" in your equation system
4570C<A * x = b>, which must be a column vector and have the same number of
4571rows as the input matrix "C<$LR_matrix>".
4572
4573The method returns a list of three items if a solution exists or an
4574empty list otherwise (!).
4575
4576Therefore, you should always use this method like this:
4577
4578 if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
4579 {
4580 # do something with the solution...
4581 }
4582 else
4583 {
4584 # do something with the fact that there is no solution...
4585 }
4586
4587The three items returned are: the dimension "C<$dimension>" of the solution
4588space (which is zero if only one solution exists, one if the solution is
4589a straight line, two if the solution is a plane, and so on), the solution
4590vector "C<$x_vector>" (which is the vector "x" of your equation system
4591C<A * x = b>) and a matrix "C<$base_matrix>" representing a base of the
4592solution space (a set of vectors which put up the solution space like
4593the spokes of an umbrella).
4594
4595Only the first "C<$dimension>" columns of this base matrix actually
4596contain entries, the remaining columns are all zero.
4597
4598Now what is all this stuff with that "base" good for?
4599
4600The output vector "x" is B<ALWAYS> a solution of your equation system
4601C<A * x = b>.
4602
4603But also any vector "C<$vector>"
4604
4605 $vector = $x_vector->clone();
4606
4607 $machine_infinity = 1E+99; # or something like that
4608
4609 for ( $i = 1; $i <= $dimension; $i++ )
4610 {
4611 $vector += rand($machine_infinity) * $base_matrix->column($i);
4612 }
4613
4614is a solution to your problem C<A * x = b>, i.e., if "C<$A_matrix>" contains
4615your matrix "A", then
4616
4617 print abs( $A_matrix * $vector - $b_vector ), "\n";
4618
4619should print a number around 1E-16 or so!
4620
4621By the way, note that you can actually calculate those vectors "C<$vector>"
4622a little more efficient as follows:
4623
4624 $rand_vector = $x_vector->shadow();
4625
4626 $machine_infinity = 1E+99; # or something like that
4627
4628 for ( $i = 1; $i <= $dimension; $i++ )
4629 {
4630 $rand_vector->assign($i,1, rand($machine_infinity) );
4631 }
4632
4633 $vector = $x_vector + ( $base_matrix * $rand_vector );
4634
4635Note that the input matrix and vector are not changed by this method
4636in any way.
4637
4638=item *
4639
4640C<$inverse_matrix = $LR_matrix-E<gt>invert_LR();>
4641
4642Use this method to calculate the inverse of a given matrix "C<$LR_matrix>",
4643which must be a (quadratic) matrix returned by the method "decompose_LR()".
4644
4645The method returns an object reference to a new matrix of the same size as
4646the input matrix containing the inverse of the matrix that you initially
4647fed into "decompose_LR()" B<IF THE INVERSE EXISTS>, or an empty list
4648otherwise.
4649
4650Therefore, you should always use this method in the following way:
4651
4652 if ( $inverse_matrix = $LR->invert_LR() )
4653 {
4654 # do something with the inverse matrix...
4655 }
4656 else
4657 {
4658 # do something with the fact that there is no inverse matrix...
4659 }
4660
4661Note that by definition (disregarding numerical errors), the product
4662of the initial matrix and its inverse (or vice-versa) is always a matrix
4663containing one's on the main diagonal (elements (1,1), (2,2), (3,3) and
4664so on) and zero's elsewhere.
4665
4666The input matrix is not changed by this method in any way.
4667
4668=item *
4669
4670C<$condition = $matrix-E<gt>condition($inverse_matrix);>
4671
4672In fact this method is just a shortcut for
4673
4674 abs($matrix) * abs($inverse_matrix)
4675
4676Both input matrices must be quadratic and have the same size, and the result
4677is meaningful only if one of them is the inverse of the other (for instance,
4678as returned by the method "invert_LR()").
4679
4680The number returned is a measure of the "condition" of the given matrix
4681"C<$matrix>", i.e., a measure of the numerical stability of the matrix.
4682
4683This number is always positive, and the smaller its value, the better the
4684condition of the matrix (the better the stability of all subsequent
4685computations carried out using this matrix).
4686
4687Numerical stability means for example that if
4688
4689 abs( $vec_correct - $vec_with_error ) < $epsilon
4690
4691holds, there must be a "C<$delta>" which doesn't depend on the vector
4692"C<$vec_correct>" (nor "C<$vec_with_error>", by the way) so that
4693
4694 abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
4695
4696also holds.
4697
4698=item *
4699
4700C<$determinant = $LR_matrix-E<gt>det_LR();>
4701
4702Calculates the determinant of a matrix, whose LR decomposition matrix
4703"C<$LR_matrix>" must be given (which must be a (quadratic) matrix
4704returned by the method "decompose_LR()").
4705
4706In fact the determinant is a by-product of the LR decomposition: It is
4707(in principle, that is, except for the sign) simply the product of the
4708elements on the main diagonal (elements (1,1), (2,2), (3,3) and so on)
4709of the LR decomposition matrix.
4710
4711(The sign is taken care of "magically" by this module)
4712
4713=item *
4714
4715C<$order = $LR_matrix-E<gt>order_LR();>
4716
4717Calculates the order (called "Rang" in German) of a matrix, whose
4718LR decomposition matrix "C<$LR_matrix>" must be given (which must
4719be a (quadratic) matrix returned by the method "decompose_LR()").
4720
4721This number is a measure of the number of linear independent row
4722and column vectors (= number of linear independent equations in
4723the case of a matrix representing an equation system) of the
4724matrix that was initially fed into "decompose_LR()".
4725
4726If "n" is the number of rows and columns of the (quadratic!) matrix,
4727then "n - order" is the dimension of the solution space of the
4728associated equation system.
4729
4730=item *
4731
4732C<$rank = $LR_matrix-E<gt>rank_LR();>
4733
4734This is an alias for the C<order_LR()> function. The "order"
4735is usually called the "rank" in the United States.
4736
4737=item *
4738
4739C<$scalar_product = $vector1-E<gt>scalar_product($vector2);>
4740
4741Returns the scalar product of vector "C<$vector1>" and vector "C<$vector2>".
4742
4743Both vectors must be column vectors (i.e., a matrix having
4744several rows but only one column).
4745
4746This is a (more efficient!) shortcut for
4747
4748 $temp = ~$vector1 * $vector2;
4749 $scalar_product = $temp->element(1,1);
4750
4751or the sum C<i=1..n> of the products C<vector1[i] * vector2[i]>.
4752
4753Provided none of the two input vectors is the null vector, then
4754the two vectors are orthogonal, i.e., have an angle of 90 degrees
4755between them, exactly when their scalar product is zero, and
4756vice-versa.
4757
4758=item *
4759
4760C<$vector_product = $vector1-E<gt>vector_product($vector2);>
4761
4762Returns the vector product of vector "C<$vector1>" and vector "C<$vector2>".
4763
4764Both vectors must be column vectors (i.e., a matrix having several rows
4765but only one column).
4766
4767Currently, the vector product is only defined for 3 dimensions (i.e.,
4768vectors with 3 rows); all other vectors trigger an error message.
4769
4770In 3 dimensions, the vector product of two vectors "x" and "y"
4771is defined as
4772
4773 | x[1] y[1] e[1] |
4774 determinant | x[2] y[2] e[2] |
4775 | x[3] y[3] e[3] |
4776
4777where the "C<x[i]>" and "C<y[i]>" are the components of the two vectors
4778"x" and "y", respectively, and the "C<e[i]>" are unity vectors (i.e.,
4779vectors with a length equal to one) with a one in row "i" and zero's
4780elsewhere (this means that you have numbers and vectors as elements
4781in this matrix!).
4782
4783This determinant evaluates to the rather simple formula
4784
4785 z[1] = x[2] * y[3] - x[3] * y[2]
4786 z[2] = x[3] * y[1] - x[1] * y[3]
4787 z[3] = x[1] * y[2] - x[2] * y[1]
4788
4789A characteristic property of the vector product is that the resulting
4790vector is orthogonal to both of the input vectors (if neither of both
4791is the null vector, otherwise this is trivial), i.e., the scalar product
4792of each of the input vectors with the resulting vector is always zero.
4793
4794=item *
4795
4796C<$length = $vector-E<gt>length();>
4797
4798This is actually a shortcut for
4799
4800 $length = sqrt( $vector->scalar_product($vector) );
4801
4802and returns the length of a given (column!) vector "C<$vector>".
4803
4804Note that the "length" calculated by this method is in fact the
4805"two"-norm of a vector "C<$vector>"!
4806
4807The general definition for norms of vectors is the following:
4808
4809 sub vector_norm
4810 {
4811 croak "Usage: \$norm = \$vector->vector_norm(\$n);"
4812 if (@_ != 2);
4813
4814 my($vector,$n) = @_;
4815 my($rows,$cols) = ($vector->[1],$vector->[2]);
4816 my($k,$comp,$sum);
4817
4818 croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
4819 unless ($cols == 1);
4820
4821 croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
4822 unless ($n > 0);
4823
4824 croak "Math::MatrixReal::vector_norm(): norm index must be integer"
4825 unless ($n == int($n));
4826
4827 $sum = 0;
4828 for ( $k = 0; $k < $rows; $k++ )
4829 {
4830 $comp = abs( $vector->[0][$k][0] );
4831 $sum += $comp ** $n;
4832 }
4833 return( $sum ** (1 / $n) );
4834 }
4835
4836Note that the case "n = 1" is the "one"-norm for matrices applied to a
4837vector, the case "n = 2" is the euclidian norm or length of a vector,
4838and if "n" goes to infinity, you have the "infinity"- or "maximum"-norm
4839for matrices applied to a vector!
4840
4841=item *
4842
4843C<$xn_vector = $matrix-E<gt>>C<solve_GSM($x0_vector,$b_vector,$epsilon);>
4844
4845=item *
4846
4847C<$xn_vector = $matrix-E<gt>>C<solve_SSM($x0_vector,$b_vector,$epsilon);>
4848
4849=item *
4850
4851C<$xn_vector = $matrix-E<gt>>C<solve_RM($x0_vector,$b_vector,$weight,$epsilon);>
4852
4853In some cases it might not be practical or desirable to solve an
4854equation system "C<A * x = b>" using an analytical algorithm like
4855the "decompose_LR()" and "solve_LR()" method pair.
4856
4857In fact in some cases, due to the numerical properties (the "condition")
4858of the matrix "A", the numerical error of the obtained result can be
4859greater than by using an approximative (iterative) algorithm like one
4860of the three implemented here.
4861
4862All three methods, GSM ("Global Step Method" or "Gesamtschrittverfahren"),
4863SSM ("Single Step Method" or "Einzelschrittverfahren") and RM ("Relaxation
4864Method" or "Relaxationsverfahren"), are fix-point iterations, that is, can
4865be described by an iteration function "C<x(t+1) = Phi( x(t) )>" which has
4866the property:
4867
4868 Phi(x) = x <==> A * x = b
4869
4870We can define "C<Phi(x)>" as follows:
4871
4872 Phi(x) := ( En - A ) * x + b
4873
4874where "En" is a matrix of the same size as "A" ("n" rows and columns)
4875with one's on its main diagonal and zero's elsewhere.
4876
4877This function has the required property.
4878
4879Proof:
4880
4881 A * x = b
4882
4883 <==> -( A * x ) = -b
4884
4885 <==> -( A * x ) + x = -b + x
4886
4887 <==> -( A * x ) + x + b = x
4888
4889 <==> x - ( A * x ) + b = x
4890
4891 <==> ( En - A ) * x + b = x
4892
4893This last step is true because
4894
4895 x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
4896
4897is the same as
4898
4899 ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
4900
4901qed
4902
4903Note that actually solving the equation system "C<A * x = b>" means
4904to calculate
4905
4906 a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] = b[i]
4907
4908 <==> a[i,i] x[i] =
4909 b[i]
4910 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
4911 + a[i,i] x[i]
4912
4913 <==> x[i] =
4914 ( b[i]
4915 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
4916 + a[i,i] x[i]
4917 ) / a[i,i]
4918
4919 <==> x[i] =
4920 ( b[i] -
4921 ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
4922 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
4923 ) / a[i,i]
4924
4925There is one major restriction, though: a fix-point iteration is
4926guaranteed to converge only if the first derivative of the iteration
4927function has an absolute value less than one in an area around the
4928point "C<x(*)>" for which "C<Phi( x(*) ) = x(*)>" is to be true, and
4929if the start vector "C<x(0)>" lies within that area!
4930
4931This is best verified grafically, which unfortunately is impossible
4932to do in this textual documentation!
4933
4934See literature on Numerical Analysis for details!
4935
4936In our case, this restriction translates to the following three conditions:
4937
4938There must exist a norm so that the norm of the matrix of the iteration
4939function, C<( En - A )>, has a value less than one, the matrix "A" may
4940not have any zero value on its main diagonal and the initial vector
4941"C<x(0)>" must be "good enough", i.e., "close enough" to the solution
4942"C<x(*)>".
4943
4944(Remember school math: the first derivative of a straight line given by
4945"C<y = a * x + b>" is "a"!)
4946
4947The three methods expect a (quadratic!) matrix "C<$matrix>" as their
4948first argument, a start vector "C<$x0_vector>", a vector "C<$b_vector>"
4949(which is the vector "b" in your equation system "C<A * x = b>"), in the
4950case of the "Relaxation Method" ("RM"), a real number "C<$weight>" best
4951between zero and two, and finally an error limit (real number) "C<$epsilon>".
4952
4953(Note that the weight "C<$weight>" used by the "Relaxation Method" ("RM")
4954is B<NOT> checked to lie within any reasonable range!)
4955
4956The three methods first test the first two conditions of the three
4957conditions listed above and return an empty list if these conditions
4958are not fulfilled.
4959
4960Therefore, you should always test their return value using some
4961code like:
4962
4963 if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
4964 {
4965 # do something with the solution...
4966 }
4967 else
4968 {
4969 # do something with the fact that there is no solution...
4970 }
4971
4972Otherwise, they iterate until C<abs( Phi(x) - x ) E<lt> epsilon>.
4973
4974(Beware that theoretically, infinite loops might result if the starting
4975vector is too far "off" the solution! In practice, this shouldn't be
4976a problem. Anyway, you can always press <ctrl-C> if you think that the
4977iteration takes too long!)
4978
4979The difference between the three methods is the following:
4980
4981In the "Global Step Method" ("GSM"), the new vector "C<x(t+1)>"
4982(called "y" here) is calculated from the vector "C<x(t)>"
4983(called "x" here) according to the formula:
4984
4985 y[i] =
4986 ( b[i]
4987 - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
4988 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
4989 ) / a[i,i]
4990
4991In the "Single Step Method" ("SSM"), the components of the vector
4992"C<x(t+1)>" which have already been calculated are used to calculate
4993the remaining components, i.e.
4994
4995 y[i] =
4996 ( b[i]
4997 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
4998 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
4999 ) / a[i,i]
5000
5001In the "Relaxation method" ("RM"), the components of the vector
5002"C<x(t+1)>" are calculated by "mixing" old and new value (like
5003cold and hot water), and the weight "C<$weight>" determines the
5004"aperture" of both the "hot water tap" as well as of the "cold
5005water tap", according to the formula:
5006
5007 y[i] =
5008 ( b[i]
5009 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
5010 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
5011 ) / a[i,i]
5012 y[i] = weight * y[i] + (1 - weight) * x[i]
5013
5014Note that the weight "C<$weight>" should be greater than zero and
5015less than two (!).
5016
5017The three methods are supposed to be of different efficiency.
5018Experiment!
5019
5020Remember that in most cases, it is probably advantageous to first
5021"normalize()" your equation system prior to solving it!
5022
5023=head1 OVERLOADED OPERATORS
5024
5025=head2 SYNOPSIS
5026
5027=over 2
5028
5029=item *
5030
5031Unary operators:
5032
5033"C<->", "C<~>", "C<abs>", C<test>, "C<!>", 'C<"">'
5034
5035=item *
5036
5037Binary (arithmetic) operators:
5038
5039"C<+>", "C<->", "C<*>", "C<**>",
5040"C<+=>", "C<-=>", "C<*=>", "C<**=>"
5041
5042=item *
5043
5044Binary (relational) operators:
5045
5046"C<==>", "C<!=>", "C<E<lt>>", "C<E<lt>=>", "C<E<gt>>", "C<E<gt>=>"
5047
5048"C<eq>", "C<ne>", "C<lt>", "C<le>", "C<gt>", "C<ge>"
5049
5050Note that the latter ("C<eq>", "C<ne>", ... ) are just synonyms
5051of the former ("C<==>", "C<!=>", ... ), defined for convenience
5052only.
5053
5054=back
5055
5056=head2 DESCRIPTION
5057
5058=over 5
5059
5060=item '-'
5061
5062Unary minus
5063
5064Returns the negative of the given matrix, i.e., the matrix with
5065all elements multiplied with the factor "-1".
5066
5067Example:
5068
5069 $matrix = -$matrix;
5070
5071=item '~'
5072
5073Transposition
5074
5075Returns the transposed of the given matrix.
5076
5077Examples:
5078
5079 $temp = ~$vector * $vector;
5080 $length = sqrt( $temp->element(1,1) );
5081
5082 if (~$matrix == $matrix) { # matrix is symmetric ... }
5083
5084=item abs
5085
5086Norm
5087
5088Returns the "one"-Norm of the given matrix.
5089
5090Example:
5091
5092 $error = abs( $A * $x - $b );
5093
5094=item test
5095
5096Boolean test
5097
5098Tests wether there is at least one non-zero element in the matrix.
5099
5100Example:
5101
5102 if ($xn_vector) { # result of iteration is not zero ... }
5103
5104=item '!'
5105
5106Negated boolean test
5107
5108Tests wether the matrix contains only zero's.
5109
5110Examples:
5111
5112 if (! $b_vector) { # heterogenous equation system ... }
5113 else { # homogenous equation system ... }
5114
5115 unless ($x_vector) { # not the null-vector! }
5116
5117=item '""""'
5118
5119"Stringify" operator
5120
5121Converts the given matrix into a string.
5122
5123Uses scientific representation to keep precision loss to a minimum in case
5124you want to read this string back in again later with "new_from_string()".
5125
5126Uses a 13-digit mantissa and a 20-character field for each element so that
5127lines will wrap nicely on an 80-column screen.
5128
5129Examples:
5130
5131 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
5132 [ 1 0 ]
5133 [ 0 -1 ]
5134 MATRIX
5135 print "$matrix";
5136
5137 [ 1.000000000000E+00 0.000000000000E+00 ]
5138 [ 0.000000000000E+00 -1.000000000000E+00 ]
5139
5140 $string = "$matrix";
5141 $test = Math::MatrixReal->new_from_string($string);
5142 if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
5143
5144=item '+'
5145
5146Addition
5147
5148Returns the sum of the two given matrices.
5149
5150Examples:
5151
5152 $matrix_S = $matrix_A + $matrix_B;
5153
5154 $matrix_A += $matrix_B;
5155
5156=item '-'
5157
5158Subtraction
5159
5160Returns the difference of the two given matrices.
5161
5162Examples:
5163
5164 $matrix_D = $matrix_A - $matrix_B;
5165
5166 $matrix_A -= $matrix_B;
5167
5168Note that this is the same as:
5169
5170 $matrix_S = $matrix_A + -$matrix_B;
5171
5172 $matrix_A += -$matrix_B;
5173
5174(The latter are less efficient, though)
5175
5176=item '*'
5177
5178Multiplication
5179
5180Returns the matrix product of the two given matrices or
5181the product of the given matrix and scalar factor.
5182
5183Examples:
5184
5185 $matrix_P = $matrix_A * $matrix_B;
5186
5187 $matrix_A *= $matrix_B;
5188
5189 $vector_b = $matrix_A * $vector_x;
5190
5191 $matrix_B = -1 * $matrix_A;
5192
5193 $matrix_B = $matrix_A * -1;
5194
5195 $matrix_A *= -1;
5196
5197=item '**'
5198
5199Exponentiation
5200
5201Returns the matrix raised to an integer power. If 0 is passed,
5202the identity matrix is returned. If a negative integer is passed,
5203it computes the inverse (if it exists) and then raised the inverse
5204to the absolute value of the integer. The matrix must be quadratic.
5205
5206Examples:
5207
5208 $matrix2 = $matrix ** 2;
5209
5210 $matrix **= 2;
5211
5212 $inv2 = $matrix ** -2;
5213
5214 $ident = $matrix ** 0;
5215
5216
5217
5218=item '=='
5219
5220Equality
5221
5222Tests two matrices for equality.
5223
5224Example:
5225
5226 if ( $A * $x == $b ) { print "EUREKA!\n"; }
5227
5228Note that in most cases, due to numerical errors (due to the finite
5229precision of computer arithmetics), it is a bad idea to compare two
5230matrices or vectors this way.
5231
5232Better use the norm of the difference of the two matrices you want
5233to compare and compare that norm with a small number, like this:
5234
5235 if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
5236
5237=item '!='
5238
5239Inequality
5240
5241Tests two matrices for inequality.
5242
5243Example:
5244
5245 while ($x0_vector != $xn_vector) { # proceed with iteration ... }
5246
5247(Stops when the iteration becomes stationary)
5248
5249Note that (just like with the '==' operator), it is usually a bad idea
5250to compare matrices or vectors this way. Compare the norm of the difference
5251of the two matrices with a small number instead.
5252
5253=item 'E<lt>'
5254
5255Less than
5256
5257Examples:
5258
5259 if ( $matrix1 < $matrix2 ) { # ... }
5260
5261 if ( $vector < $epsilon ) { # ... }
5262
5263 if ( 1E-12 < $vector ) { # ... }
5264
5265 if ( $A * $x - $b < 1E-12 ) { # ... }
5266
5267These are just shortcuts for saying:
5268
5269 if ( abs($matrix1) < abs($matrix2) ) { # ... }
5270
5271 if ( abs($vector) < abs($epsilon) ) { # ... }
5272
5273 if ( abs(1E-12) < abs($vector) ) { # ... }
5274
5275 if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
5276
5277Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
5278
5279=item 'E<lt>='
5280
5281Less than or equal
5282
5283As with the '<' operator, this is just a shortcut for the same expression
5284with "abs()" around all arguments.
5285
5286Example:
5287
5288 if ( $A * $x - $b <= 1E-12 ) { # ... }
5289
5290which in fact is the same as:
5291
5292 if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
5293
5294Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
5295
5296=item 'E<gt>'
5297
5298Greater than
5299
5300As with the '<' and '<=' operator, this
5301
5302 if ( $xn - $x0 > 1E-12 ) { # ... }
5303
5304is just a shortcut for:
5305
5306 if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
5307
5308Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
5309
5310=item 'E<gt>='
5311
5312Greater than or equal
5313
5314As with the '<', '<=' and '>' operator, the following
5315
5316 if ( $LR >= $A ) { # ... }
5317
5318is simply a shortcut for:
5319
5320 if ( abs($LR) >= abs($A) ) { # ... }
5321
5322Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
5323
5324=back
5325
5326=head1 SEE ALSO
5327
5328Math::VectorReal, Math::PARI, Math::MatrixBool,
5329DFA::Kleene, Math::Kleene,
5330Set::IntegerRange, Set::IntegerFast .
5331
5332=head1 VERSION
5333
5334This man page documents Math::MatrixReal version 1.9.
5335
5336The latest version can always be found at
5337http://leto.net/code/Math-MatrixReal/
5338
5339=head1 AUTHORS
5340
5341Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>,
5342Jonathan Leto <jonathan@leto.net>.
5343
5344Currently maintained by Jonathan Leto, send all bugs/patches to me.
5345
5346=head1 CREDITS
5347
5348Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
5349Algebra and Linear Algebra at the university (RWTH Aachen, Germany), and
5350to Prof. Esser and his assistant, Mr. Jarausch, for their fascinating
5351lectures in Numerical Analysis!
5352
5353=head1 COPYRIGHT
5354
5355Copyright (c) 1996-2002 by Steffen Beyer, Rodolphe Ortalo, Jonathan Leto.
5356All rights reserved.
5357
5358=head1 LICENSE AGREEMENT
5359
5360This package is free software; you can redistribute it and/or
5361modify it under the same terms as Perl itself.
5362
Full OpenSPARC design & verification tarball including full HDL.