[system] / trunk / pg / lib / MatrixReal1.pm Repository:
ViewVC logotype

Annotation of /trunk/pg/lib/MatrixReal1.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3372 - (view) (download) (as text)

1 : sh002i 1050 # Copyright (c) 1996, 1997 by Steffen Beyer. All rights reserved.
2 :     # Copyright (c) 1999 by Rodolphe Ortalo. All rights reserved.
3 :     # This package is free software; you can redistribute it and/or
4 :     # modify it under the same terms as Perl itself.
5 :    
6 :     # slightly modified for use in WeBWorK
7 :     # modifications by Michael E Gage -- added a reference to options in the object array ($this)
8 :     # a better approach would be to rewrite this package so that $this is a hash rather than an array
9 :     # grep for MEG to see changes.
10 :    
11 :     # Changed package name to MatrixReal1 throughout.
12 :     package MatrixReal1;
13 :    
14 :     use strict;
15 :     use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);
16 :    
17 :     require Exporter;
18 :    
19 :     @ISA = qw(Exporter);
20 :    
21 :     @EXPORT = qw();
22 :    
23 :     @EXPORT_OK = qw(min max);
24 :    
25 :     %EXPORT_TAGS = (all => [@EXPORT_OK]);
26 :    
27 :     $VERSION = '1.3a5';
28 : gage 3372 $MatrixReal1::OPTION_ENTRY=7;
29 : sh002i 1050
30 :     use Carp;
31 :    
32 :     use overload
33 :     'neg' => '_negate',
34 :     '~' => '_transpose',
35 :     'bool' => '_boolean',
36 :     '!' => '_not_boolean',
37 :     '""' => '_stringify',
38 :     'abs' => '_norm',
39 :     '+' => '_add',
40 :     '-' => '_subtract',
41 :     '*' => '_multiply',
42 :     '+=' => '_assign_add',
43 :     '-=' => '_assign_subtract',
44 :     '*=' => '_assign_multiply',
45 :     '==' => '_equal',
46 :     '!=' => '_not_equal',
47 :     '<' => '_less_than',
48 :     '<=' => '_less_than_or_equal',
49 :     '>' => '_greater_than',
50 :     '>=' => '_greater_than_or_equal',
51 :     'eq' => '_equal',
52 :     'ne' => '_not_equal',
53 :     'lt' => '_less_than',
54 :     'le' => '_less_than_or_equal',
55 :     'gt' => '_greater_than',
56 :     'ge' => '_greater_than_or_equal',
57 :     '=' => '_clone',
58 :     'fallback' => undef;
59 :    
60 :     sub new
61 :     {
62 :     croak "Usage: \$new_matrix = MatrixReal1->new(\$rows,\$columns);"
63 :     if (@_ != 3);
64 :    
65 :     my $proto = shift;
66 :     my $class = ref($proto) || $proto || 'MatrixReal1';
67 :     my $rows = shift;
68 :     my $cols = shift;
69 :     my($i,$j);
70 :     my($this);
71 :    
72 :     croak "MatrixReal1::new(): number of rows must be > 0"
73 :     if ($rows <= 0);
74 :    
75 :     croak "MatrixReal1::new(): number of columns must be > 0"
76 :     if ($cols <= 0);
77 :    
78 :     # $this = [ [ ], $rows, $cols ];
79 : gage 3372 $this = [ [ ], $rows, $cols ];
80 :     $this->[$MatrixReal1::OPTION_ENTRY] = {}; # added a holder for options MEG
81 :     # see also modifications to LR decomposition
82 : sh002i 1050
83 :     # Creates first empty row
84 :     my $empty = [ ];
85 :     $#$empty = $cols - 1; # Lengthens the array
86 :     for (my $j = 0; $j < $cols; $j++)
87 :     {
88 :     $empty->[$j] = 0.0;
89 :     }
90 :     $this->[0][0] = $empty;
91 :     # Creates other rows (by copying)
92 :     for (my $i = 1; $i < $rows; $i++)
93 :     {
94 :     my $arow = [ ];
95 :     @$arow = @$empty;
96 :     $this->[0][$i] = $arow;
97 :     }
98 :     bless($this, $class);
99 :     return($this);
100 :     }
101 :    
102 :     sub new_from_string
103 :     {
104 :     croak "Usage: \$new_matrix = MatrixReal1->new_from_string(\$string);"
105 :     if (@_ != 2);
106 :    
107 :     my $proto = shift;
108 :     my $class = ref($proto) || $proto || 'MatrixReal1';
109 :     my $string = shift;
110 :     my($line,$values);
111 :     my($rows,$cols);
112 :     my($row,$col);
113 :     my($warn);
114 :     my($this);
115 :    
116 :     $warn = 0;
117 :     $rows = 0;
118 :     $cols = 0;
119 :     $values = [ ];
120 :     while ($string =~ m!^\s*
121 :     \[ \s+ ( (?: [+-]? \d+ (?: \. \d* )? (?: E [+-]? \d+ )? \s+ )+ ) \] \s*? \n
122 :     !x)
123 :     {
124 :     $line = $1;
125 :     $string = $';
126 :     $values->[$rows] = [ ];
127 :     @{$values->[$rows]} = split(' ', $line);
128 :     $col = @{$values->[$rows]};
129 :     if ($col != $cols)
130 :     {
131 :     unless ($cols == 0) { $warn = 1; }
132 :     if ($col > $cols) { $cols = $col; }
133 :     }
134 :     $rows++;
135 :     }
136 :     if ($string !~ m!^\s*$!)
137 :     {
138 :     croak "MatrixReal1::new_from_string(): syntax error in input string";
139 :     }
140 :     if ($rows == 0)
141 :     {
142 :     croak "MatrixReal1::new_from_string(): empty input string";
143 :     }
144 :     if ($warn)
145 :     {
146 :     warn "MatrixReal1::new_from_string(): missing elements will be set to zero!\n";
147 :     }
148 :     $this = MatrixReal1::new($class,$rows,$cols);
149 :     for ( $row = 0; $row < $rows; $row++ )
150 :     {
151 :     for ( $col = 0; $col < @{$values->[$row]}; $col++ )
152 :     {
153 :     $this->[0][$row][$col] = $values->[$row][$col];
154 :     }
155 :     }
156 :     return($this);
157 :     }
158 :    
159 :     sub shadow
160 :     {
161 :     croak "Usage: \$new_matrix = \$some_matrix->shadow();"
162 :     if (@_ != 1);
163 :    
164 :     my($matrix) = @_;
165 :     my($temp);
166 :    
167 :     $temp = $matrix->new($matrix->[1],$matrix->[2]);
168 :     return($temp);
169 :     }
170 :    
171 :    
172 :     sub copy
173 :     {
174 :     croak "Usage: \$matrix1->copy(\$matrix2);"
175 :     if (@_ != 2);
176 :    
177 :     my($matrix1,$matrix2) = @_;
178 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
179 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
180 :     my($i,$j);
181 :    
182 :     croak "MatrixReal1::copy(): matrix size mismatch"
183 :     unless (($rows1 == $rows2) && ($cols1 == $cols2));
184 :    
185 :     for ( $i = 0; $i < $rows1; $i++ )
186 :     {
187 :     my $r1 = []; # New array ref
188 :     my $r2 = $matrix2->[0][$i];
189 :     @$r1 = @$r2; # Copy whole array directly
190 :     $matrix1->[0][$i] = $r1;
191 :     }
192 :     $matrix1->[3] = $matrix2->[3]; # sign or option
193 :     if (defined $matrix2->[4]) # is an LR decomposition matrix!
194 :     {
195 :     # $matrix1->[3] = $matrix2->[3]; # $sign
196 :     $matrix1->[4] = $matrix2->[4]; # $perm_row
197 :     $matrix1->[5] = $matrix2->[5]; # $perm_col
198 :     $matrix1->[6] = $matrix2->[6]; # $option
199 :     }
200 :     }
201 :    
202 :     sub clone
203 :     {
204 :     croak "Usage: \$twin_matrix = \$some_matrix->clone();"
205 :     if (@_ != 1);
206 :    
207 :     my($matrix) = @_;
208 :     my($temp);
209 :    
210 :     $temp = $matrix->new($matrix->[1],$matrix->[2]);
211 :     $temp->copy($matrix);
212 :     return($temp);
213 :     }
214 :    
215 :     sub row
216 :     {
217 :     croak "Usage: \$row_vector = \$matrix->row(\$row);"
218 :     if (@_ != 2);
219 :    
220 :     my($matrix,$row) = @_;
221 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
222 :     my($temp);
223 :     my($j);
224 :    
225 :     croak "MatrixReal1::row(): row index out of range"
226 :     if (($row < 1) || ($row > $rows));
227 :    
228 :     $row--;
229 :     $temp = $matrix->new(1,$cols);
230 :     for ( $j = 0; $j < $cols; $j++ )
231 :     {
232 :     $temp->[0][0][$j] = $matrix->[0][$row][$j];
233 :     }
234 :     return($temp);
235 :     }
236 :    
237 :     sub column
238 :     {
239 :     croak "Usage: \$column_vector = \$matrix->column(\$column);"
240 :     if (@_ != 2);
241 :    
242 :     my($matrix,$col) = @_;
243 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
244 :     my($temp);
245 :     my($i);
246 :    
247 :     croak "MatrixReal1::column(): column index out of range"
248 :     if (($col < 1) || ($col > $cols));
249 :    
250 :     $col--;
251 :     $temp = $matrix->new($rows,1);
252 :     for ( $i = 0; $i < $rows; $i++ )
253 :     {
254 :     $temp->[0][$i][0] = $matrix->[0][$i][$col];
255 :     }
256 :     return($temp);
257 :     }
258 :    
259 : gage 3372 sub _undo_LR # I don't think gives original matrix. LR is not the same as the original matrix
260 : sh002i 1050 {
261 :     croak "Usage: \$matrix->_undo_LR();"
262 :     if (@_ != 1);
263 :    
264 :     my($this) = @_;
265 : gage 3372 my $rh_options = $this->[$MatrixReal1::OPTION_ENTRY];
266 : sh002i 1050 undef $this->[3];
267 :     undef $this->[4];
268 :     undef $this->[5];
269 :     undef $this->[6];
270 : gage 3372 $this->[$MatrixReal1::OPTION_ENTRY] = $rh_options;
271 : sh002i 1050 }
272 :    
273 :     sub zero
274 :     {
275 :     croak "Usage: \$matrix->zero();"
276 :     if (@_ != 1);
277 :    
278 :     my($this) = @_;
279 :     my($rows,$cols) = ($this->[1],$this->[2]);
280 :     my($i,$j);
281 :    
282 :     $this->_undo_LR();
283 :    
284 :     # Zero first row
285 :     for (my $j = 0; $j < $cols; $j++ )
286 :     {
287 :     $this->[0][0][$j] = 0.0;
288 :     }
289 :     # Then propagate to other rows
290 :     for (my $i = 0; $i < $rows; $i++)
291 :     {
292 :     @{$this->[0][$i]} = @{$this->[0][0]};
293 :     }
294 :     }
295 :    
296 :     sub one
297 :     {
298 :     croak "Usage: \$matrix->one();"
299 :     if (@_ != 1);
300 :    
301 :     my($this) = @_;
302 :     my($rows,$cols) = ($this->[1],$this->[2]);
303 :     my($i,$j);
304 :    
305 :     # No need for this: done by the 'zero()'
306 :     # $this->_undo_LR();
307 :     $this->zero(); # We rely on zero() efficiency
308 :     for (my $i = 0; $i < $rows; $i++ )
309 :     {
310 :     $this->[0][$i][$i] = 1.0;
311 :     }
312 :     }
313 :    
314 :     sub assign
315 :     {
316 :     croak "Usage: \$matrix->assign(\$row,\$column,\$value);"
317 :     if (@_ != 4);
318 :    
319 :     my($this,$row,$col,$value) = @_;
320 :     my($rows,$cols) = ($this->[1],$this->[2]);
321 :    
322 :     croak "MatrixReal1::assign(): row index out of range"
323 :     if (($row < 1) || ($row > $rows));
324 :    
325 :     croak "MatrixReal1::assign(): column index out of range"
326 :     if (($col < 1) || ($col > $cols));
327 :    
328 :     $this->_undo_LR();
329 :    
330 :     $this->[0][--$row][--$col] = $value;
331 :     }
332 :    
333 :     sub element
334 :     {
335 :     croak "Usage: \$value = \$matrix->element(\$row,\$column);"
336 :     if (@_ != 3);
337 :    
338 :     my($this,$row,$col) = @_;
339 :     my($rows,$cols) = ($this->[1],$this->[2]);
340 :    
341 :     croak "MatrixReal1::element(): row index out of range"
342 :     if (($row < 1) || ($row > $rows));
343 :    
344 :     croak "MatrixReal1::element(): column index out of range"
345 :     if (($col < 1) || ($col > $cols));
346 :    
347 :     return( $this->[0][--$row][--$col] );
348 :     }
349 :    
350 :     sub dim # returns dimensions of a matrix
351 :     {
352 :     croak "Usage: (\$rows,\$columns) = \$matrix->dim();"
353 :     if (@_ != 1);
354 :    
355 :     my($matrix) = @_;
356 :    
357 :     return( $matrix->[1], $matrix->[2] );
358 :     }
359 :    
360 :     sub norm_one # maximum of sums of each column
361 :     {
362 :     croak "Usage: \$norm_one = \$matrix->norm_one();"
363 :     if (@_ != 1);
364 :    
365 :     my($this) = @_;
366 :     my($rows,$cols) = ($this->[1],$this->[2]);
367 :    
368 :     my $max = 0.0;
369 :     for (my $j = 0; $j < $cols; $j++)
370 :     {
371 :     my $sum = 0.0;
372 :     for (my $i = 0; $i < $rows; $i++)
373 :     {
374 :     $sum += abs( $this->[0][$i][$j] );
375 :     }
376 :     $max = $sum if ($sum > $max);
377 :     }
378 :     return($max);
379 :     }
380 :    
381 :     sub norm_max # maximum of sums of each row
382 :     {
383 :     croak "Usage: \$norm_max = \$matrix->norm_max();"
384 :     if (@_ != 1);
385 :    
386 :     my($this) = @_;
387 :     my($rows,$cols) = ($this->[1],$this->[2]);
388 :    
389 :     my $max = 0.0;
390 :     for (my $i = 0; $i < $rows; $i++)
391 :     {
392 :     my $sum = 0.0;
393 :     for (my $j = 0; $j < $cols; $j++)
394 :     {
395 :     $sum += abs( $this->[0][$i][$j] );
396 :     }
397 :     $max = $sum if ($sum > $max);
398 :     }
399 :     return($max);
400 :     }
401 :    
402 :     sub negate
403 :     {
404 :     croak "Usage: \$matrix1->negate(\$matrix2);"
405 :     if (@_ != 2);
406 :    
407 :     my($matrix1,$matrix2) = @_;
408 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
409 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
410 :    
411 :     croak "MatrixReal1::negate(): matrix size mismatch"
412 :     unless (($rows1 == $rows2) && ($cols1 == $cols2));
413 :    
414 :     $matrix1->_undo_LR();
415 :    
416 :     for (my $i = 0; $i < $rows1; $i++ )
417 :     {
418 :     for (my $j = 0; $j < $cols1; $j++ )
419 :     {
420 :     $matrix1->[0][$i][$j] = -($matrix2->[0][$i][$j]);
421 :     }
422 :     }
423 :     }
424 :    
425 :     sub transpose
426 :     {
427 :     croak "Usage: \$matrix1->transpose(\$matrix2);"
428 :     if (@_ != 2);
429 :    
430 :     my($matrix1,$matrix2) = @_;
431 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
432 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
433 :    
434 :     croak "MatrixReal1::transpose(): matrix size mismatch"
435 :     unless (($rows1 == $cols2) && ($cols1 == $rows2));
436 :    
437 :     $matrix1->_undo_LR();
438 :    
439 :     if ($rows1 == $cols1)
440 :     {
441 :     # more complicated to make in-place possible!
442 :    
443 :     for (my $i = 0; $i < $rows1; $i++)
444 :     {
445 :     for (my $j = ($i + 1); $j < $cols1; $j++)
446 :     {
447 :     my $swap = $matrix2->[0][$i][$j];
448 :     $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
449 :     $matrix1->[0][$j][$i] = $swap;
450 :     }
451 :     $matrix1->[0][$i][$i] = $matrix2->[0][$i][$i];
452 :     }
453 :     }
454 :     else # ($rows1 != $cols1)
455 :     {
456 :     for (my $i = 0; $i < $rows1; $i++)
457 :     {
458 :     for (my $j = 0; $j < $cols1; $j++)
459 :     {
460 :     $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
461 :     }
462 :     }
463 :     }
464 :     }
465 :    
466 :     sub add
467 :     {
468 :     croak "Usage: \$matrix1->add(\$matrix2,\$matrix3);"
469 :     if (@_ != 3);
470 :    
471 :     my($matrix1,$matrix2,$matrix3) = @_;
472 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
473 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
474 :     my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
475 :     my($i,$j);
476 :    
477 :     croak "MatrixReal1::add(): matrix size mismatch"
478 :     unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
479 :     ($cols1 == $cols2) && ($cols1 == $cols3));
480 :    
481 :     $matrix1->_undo_LR();
482 :    
483 :     for ( $i = 0; $i < $rows1; $i++ )
484 :     {
485 :     for ( $j = 0; $j < $cols1; $j++ )
486 :     {
487 :     $matrix1->[0][$i][$j] =
488 :     $matrix2->[0][$i][$j] + $matrix3->[0][$i][$j];
489 :     }
490 :     }
491 :     }
492 :    
493 :     sub subtract
494 :     {
495 :     croak "Usage: \$matrix1->subtract(\$matrix2,\$matrix3);"
496 :     if (@_ != 3);
497 :    
498 :     my($matrix1,$matrix2,$matrix3) = @_;
499 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
500 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
501 :     my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
502 :     my($i,$j);
503 :    
504 :     croak "MatrixReal1::subtract(): matrix size mismatch"
505 :     unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
506 :     ($cols1 == $cols2) && ($cols1 == $cols3));
507 :    
508 :     $matrix1->_undo_LR();
509 :    
510 :     for ( $i = 0; $i < $rows1; $i++ )
511 :     {
512 :     for ( $j = 0; $j < $cols1; $j++ )
513 :     {
514 :     $matrix1->[0][$i][$j] =
515 :     $matrix2->[0][$i][$j] - $matrix3->[0][$i][$j];
516 :     }
517 :     }
518 :     }
519 :    
520 :     sub multiply_scalar
521 :     {
522 :     croak "Usage: \$matrix1->multiply_scalar(\$matrix2,\$scalar);"
523 :     if (@_ != 3);
524 :    
525 :     my($matrix1,$matrix2,$scalar) = @_;
526 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
527 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
528 :     my($i,$j);
529 :    
530 :     croak "MatrixReal1::multiply_scalar(): matrix size mismatch"
531 :     unless (($rows1 == $rows2) && ($cols1 == $cols2));
532 :    
533 :     $matrix1->_undo_LR();
534 :    
535 : gage 3372 for ( my $i = 0; $i < $rows1; $i++ )
536 : sh002i 1050 {
537 : gage 3372 for ( my $j = 0; $j < $cols1; $j++ )
538 : sh002i 1050 {
539 :     $matrix1->[0][$i][$j] = $matrix2->[0][$i][$j] * $scalar;
540 :     }
541 :     }
542 :     }
543 :    
544 :     sub multiply
545 :     {
546 :     croak "Usage: \$product_matrix = \$matrix1->multiply(\$matrix2);"
547 :     if (@_ != 2);
548 :    
549 :     my($matrix1,$matrix2) = @_;
550 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
551 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
552 :     my($temp);
553 :    
554 :     croak "MatrixReal1::multiply(): matrix size mismatch"
555 :     unless ($cols1 == $rows2);
556 :    
557 :     $temp = $matrix1->new($rows1,$cols2);
558 :     for (my $i = 0; $i < $rows1; $i++ )
559 :     {
560 :     for (my $j = 0; $j < $cols2; $j++ )
561 :     {
562 :     my $sum = 0.0;
563 :     for (my $k = 0; $k < $cols1; $k++ )
564 :     {
565 :     $sum += ( $matrix1->[0][$i][$k] * $matrix2->[0][$k][$j] );
566 :     }
567 :     $temp->[0][$i][$j] = $sum;
568 :     }
569 :     }
570 :     return($temp);
571 :     }
572 :    
573 :     sub min
574 :     {
575 :     croak "Usage: \$minimum = MatrixReal1::min(\$number1,\$number2);"
576 :     if (@_ != 2);
577 :    
578 :     return( $_[0] < $_[1] ? $_[0] : $_[1] );
579 :     }
580 :    
581 :     sub max
582 :     {
583 :     croak "Usage: \$maximum = MatrixReal1::max(\$number1,\$number2);"
584 :     if (@_ != 2);
585 :    
586 :     return( $_[0] > $_[1] ? $_[0] : $_[1] );
587 :     }
588 :    
589 :     sub kleene
590 :     {
591 :     croak "Usage: \$minimal_cost_matrix = \$cost_matrix->kleene();"
592 :     if (@_ != 1);
593 :    
594 :     my($matrix) = @_;
595 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
596 :     my($i,$j,$k,$n);
597 :     my($temp);
598 :    
599 :     croak "MatrixReal1::kleene(): matrix is not quadratic"
600 :     unless ($rows == $cols);
601 :    
602 :     $temp = $matrix->new($rows,$cols);
603 :     $temp->copy($matrix);
604 :     $temp->_undo_LR();
605 :     $n = $rows;
606 :     for ( $i = 0; $i < $n; $i++ )
607 :     {
608 :     $temp->[0][$i][$i] = min( $temp->[0][$i][$i] , 0 );
609 :     }
610 :     for ( $k = 0; $k < $n; $k++ )
611 :     {
612 :     for ( $i = 0; $i < $n; $i++ )
613 :     {
614 :     for ( $j = 0; $j < $n; $j++ )
615 :     {
616 :     $temp->[0][$i][$j] = min( $temp->[0][$i][$j] ,
617 :     ( $temp->[0][$i][$k] +
618 :     $temp->[0][$k][$j] ) );
619 :     }
620 :     }
621 :     }
622 :     return($temp);
623 :     }
624 :    
625 :     sub normalize
626 :     {
627 :     croak "Usage: (\$norm_matrix,\$norm_vector) = \$matrix->normalize(\$vector);"
628 :     if (@_ != 2);
629 :    
630 :     my($matrix,$vector) = @_;
631 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
632 :     my($norm_matrix,$norm_vector);
633 :     my($max,$val);
634 :     my($i,$j,$n);
635 :    
636 :     croak "MatrixReal1::normalize(): matrix is not quadratic"
637 :     unless ($rows == $cols);
638 :    
639 :     $n = $rows;
640 :    
641 :     croak "MatrixReal1::normalize(): vector is not a column vector"
642 :     unless ($vector->[2] == 1);
643 :    
644 :     croak "MatrixReal1::normalize(): matrix and vector size mismatch"
645 :     unless ($vector->[1] == $n);
646 :    
647 :     $norm_matrix = $matrix->new($n,$n);
648 :     $norm_vector = $vector->new($n,1);
649 :    
650 :     $norm_matrix->copy($matrix);
651 :     $norm_vector->copy($vector);
652 :    
653 :     $norm_matrix->_undo_LR();
654 :    
655 :     for ( $i = 0; $i < $n; $i++ )
656 :     {
657 :     $max = abs($norm_vector->[0][$i][0]);
658 :     for ( $j = 0; $j < $n; $j++ )
659 :     {
660 :     $val = abs($norm_matrix->[0][$i][$j]);
661 :     if ($val > $max) { $max = $val; }
662 :     }
663 :     if ($max != 0)
664 :     {
665 :     $norm_vector->[0][$i][0] /= $max;
666 :     for ( $j = 0; $j < $n; $j++ )
667 :     {
668 :     $norm_matrix->[0][$i][$j] /= $max;
669 :     }
670 :     }
671 :     }
672 :     return($norm_matrix,$norm_vector);
673 :     }
674 :    
675 :     sub decompose_LR
676 :     {
677 :     croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
678 :     if (@_ != 1);
679 :    
680 :     my($matrix) = @_;
681 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
682 :     my($perm_row,$perm_col);
683 :     my($row,$col,$max);
684 :     my($i,$j,$k,$n);
685 :     my($sign) = 1;
686 :     my($swap);
687 :     my($temp);
688 :    
689 :     croak "MatrixReal1::decompose_LR(): matrix is not quadratic"
690 :     unless ($rows == $cols);
691 :    
692 :     $temp = $matrix->new($rows,$cols);
693 :     $temp->copy($matrix);
694 :     $n = $rows;
695 :     $perm_row = [ ];
696 :     $perm_col = [ ];
697 :     for ( $i = 0; $i < $n; $i++ )
698 :     {
699 :     $perm_row->[$i] = $i;
700 :     $perm_col->[$i] = $i;
701 :     }
702 :     NONZERO:
703 :     for ( $k = 0; $k < $n; $k++ ) # use Gauss's algorithm:
704 :     {
705 :     # complete pivot-search:
706 :    
707 :     $max = 0;
708 :     for ( $i = $k; $i < $n; $i++ )
709 :     {
710 :     for ( $j = $k; $j < $n; $j++ )
711 :     {
712 :     if (($swap = abs($temp->[0][$i][$j])) > $max)
713 :     {
714 :     $max = $swap;
715 :     $row = $i;
716 :     $col = $j;
717 :     }
718 :     }
719 :     }
720 :     last NONZERO if ($max == 0); # (all remaining elements are zero)
721 :     if ($k != $row) # swap row $k and row $row:
722 :     {
723 :     $sign = -$sign;
724 :     $swap = $perm_row->[$k];
725 :     $perm_row->[$k] = $perm_row->[$row];
726 :     $perm_row->[$row] = $swap;
727 :     for ( $j = 0; $j < $n; $j++ )
728 :     {
729 :     # (must run from 0 since L has to be swapped too!)
730 :    
731 :     $swap = $temp->[0][$k][$j];
732 :     $temp->[0][$k][$j] = $temp->[0][$row][$j];
733 :     $temp->[0][$row][$j] = $swap;
734 :     }
735 :     }
736 :     if ($k != $col) # swap column $k and column $col:
737 :     {
738 :     $sign = -$sign;
739 :     $swap = $perm_col->[$k];
740 :     $perm_col->[$k] = $perm_col->[$col];
741 :     $perm_col->[$col] = $swap;
742 :     for ( $i = 0; $i < $n; $i++ )
743 :     {
744 :     $swap = $temp->[0][$i][$k];
745 :     $temp->[0][$i][$k] = $temp->[0][$i][$col];
746 :     $temp->[0][$i][$col] = $swap;
747 :     }
748 :     }
749 :     for ( $i = ($k + 1); $i < $n; $i++ )
750 :     {
751 :     # scan the remaining rows, add multiples of row $k to row $i:
752 :    
753 :     $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
754 :     if ($swap != 0)
755 :     {
756 :     # calculate a row of matrix R:
757 :    
758 :     for ( $j = ($k + 1); $j < $n; $j++ )
759 :     {
760 :     $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
761 :     }
762 :    
763 :     # store matrix L in same matrix as R:
764 :    
765 :     $temp->[0][$i][$k] = $swap;
766 :     }
767 :     }
768 :     }
769 :     my $rh_options = $temp->[3];
770 :     $temp->[3] = $sign;
771 :     $temp->[4] = $perm_row;
772 :     $temp->[5] = $perm_col;
773 :     $temp->[6] = $temp->[3];
774 :     return($temp);
775 :     }
776 :    
777 :     sub solve_LR
778 :     {
779 :     croak "Usage: (\$dimension,\$x_vector,\$base_matrix) = \$LR_matrix->solve_LR(\$b_vector);"
780 :     if (@_ != 2);
781 :    
782 :     my($LR_matrix,$b_vector) = @_;
783 :     my($rows,$cols) = ($LR_matrix->[1],$LR_matrix->[2]);
784 :     my($dimension,$x_vector,$base_matrix);
785 :     my($perm_row,$perm_col);
786 :     my($y_vector,$sum);
787 :     my($i,$j,$k,$n);
788 :    
789 :     croak "MatrixReal1::solve_LR(): not an LR decomposition matrix"
790 :     unless ((defined $LR_matrix->[4]) && ($rows == $cols));
791 :    
792 :     $n = $rows;
793 :    
794 :     croak "MatrixReal1::solve_LR(): vector is not a column vector"
795 :     unless ($b_vector->[2] == 1);
796 :    
797 :     croak "MatrixReal1::solve_LR(): matrix and vector size mismatch"
798 :     unless ($b_vector->[1] == $n);
799 :    
800 :     $perm_row = $LR_matrix->[4];
801 :     $perm_col = $LR_matrix->[5];
802 :    
803 :     $x_vector = $b_vector->new($n,1);
804 :     $y_vector = $b_vector->new($n,1);
805 :     $base_matrix = $LR_matrix->new($n,$n);
806 :    
807 :     # calculate "x" so that LRx = b ==> calculate Ly = b, Rx = y:
808 :    
809 :     for ( $i = 0; $i < $n; $i++ ) # calculate $y_vector:
810 :     {
811 :     $sum = $b_vector->[0][($perm_row->[$i])][0];
812 :     for ( $j = 0; $j < $i; $j++ )
813 :     {
814 :     $sum -= $LR_matrix->[0][$i][$j] * $y_vector->[0][$j][0];
815 :     }
816 :     $y_vector->[0][$i][0] = $sum;
817 :     }
818 :    
819 :     $dimension = 0;
820 :     for ( $i = ($n - 1); $i >= 0; $i-- ) # calculate $x_vector:
821 :     {
822 :     if ($LR_matrix->[0][$i][$i] == 0)
823 :     {
824 :     if ($y_vector->[0][$i][0] != 0)
825 :     {
826 :     return(); # a solution does not exist!
827 :     }
828 :     else
829 :     {
830 :     $dimension++;
831 :     $x_vector->[0][($perm_col->[$i])][0] = 0;
832 :     }
833 :     }
834 :     else
835 :     {
836 :     $sum = $y_vector->[0][$i][0];
837 :     for ( $j = ($i + 1); $j < $n; $j++ )
838 :     {
839 :     $sum -= $LR_matrix->[0][$i][$j] *
840 :     $x_vector->[0][($perm_col->[$j])][0];
841 :     }
842 :     $x_vector->[0][($perm_col->[$i])][0] =
843 :     $sum / $LR_matrix->[0][$i][$i];
844 :     }
845 :     }
846 :     if ($dimension)
847 :     {
848 :     if ($dimension == $n)
849 :     {
850 :     $base_matrix->one();
851 :     }
852 :     else
853 :     {
854 :     for ( $k = 0; $k < $dimension; $k++ )
855 :     {
856 :     $base_matrix->[0][($perm_col->[($n-$k-1)])][$k] = 1;
857 :     for ( $i = ($n-$dimension-1); $i >= 0; $i-- )
858 :     {
859 :     $sum = 0;
860 :     for ( $j = ($i + 1); $j < $n; $j++ )
861 :     {
862 :     $sum -= $LR_matrix->[0][$i][$j] *
863 :     $base_matrix->[0][($perm_col->[$j])][$k];
864 :     }
865 :     $base_matrix->[0][($perm_col->[$i])][$k] =
866 :     $sum / $LR_matrix->[0][$i][$i];
867 :     }
868 :     }
869 :     }
870 :     }
871 :     return( $dimension, $x_vector, $base_matrix );
872 :     }
873 :    
874 :     sub invert_LR
875 :     {
876 :     croak "Usage: \$inverse_matrix = \$LR_matrix->invert_LR();"
877 :     if (@_ != 1);
878 :    
879 :     my($matrix) = @_;
880 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
881 :     my($inv_matrix,$x_vector,$y_vector);
882 :     my($i,$j,$n);
883 :    
884 :     croak "MatrixReal1::invert_LR(): not an LR decomposition matrix"
885 :     unless ((defined $matrix->[4]) && ($rows == $cols));
886 :    
887 :     $n = $rows;
888 :     if ($matrix->[0][$n-1][$n-1] != 0)
889 :     {
890 :     $inv_matrix = $matrix->new($n,$n);
891 :     $y_vector = $matrix->new($n,1);
892 :     for ( $j = 0; $j < $n; $j++ )
893 :     {
894 :     if ($j > 0)
895 :     {
896 :     $y_vector->[0][$j-1][0] = 0;
897 :     }
898 :     $y_vector->[0][$j][0] = 1;
899 :     if (($rows,$x_vector,$cols) = $matrix->solve_LR($y_vector))
900 :     {
901 :     for ( $i = 0; $i < $n; $i++ )
902 :     {
903 :     $inv_matrix->[0][$i][$j] = $x_vector->[0][$i][0];
904 :     }
905 :     }
906 :     else
907 :     {
908 :     die "MatrixReal1::invert_LR(): unexpected error - please inform author!\n";
909 :     }
910 :     }
911 :     return($inv_matrix);
912 :     }
913 :     else { return(); } # matrix is not invertible!
914 :     }
915 :    
916 :     sub condition
917 :     {
918 :     # 1st matrix MUST be the inverse of 2nd matrix (or vice-versa)
919 :     # for a meaningful result!
920 :    
921 :     croak "Usage: \$condition = \$matrix->condition(\$inverse_matrix);"
922 :     if (@_ != 2);
923 :    
924 :     my($matrix1,$matrix2) = @_;
925 :     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
926 :     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
927 :    
928 :     croak "MatrixReal1::condition(): 1st matrix is not quadratic"
929 :     unless ($rows1 == $cols1);
930 :    
931 :     croak "MatrixReal1::condition(): 2nd matrix is not quadratic"
932 :     unless ($rows2 == $cols2);
933 :    
934 :     croak "MatrixReal1::condition(): matrix size mismatch"
935 :     unless (($rows1 == $rows2) && ($cols1 == $cols2));
936 :    
937 :     return( $matrix1->norm_one() * $matrix2->norm_one() );
938 :     }
939 :    
940 :     sub det_LR # determinant of LR decomposition matrix
941 :     {
942 :     croak "Usage: \$determinant = \$LR_matrix->det_LR();"
943 :     if (@_ != 1);
944 :    
945 :     my($matrix) = @_;
946 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
947 :     my($k,$det);
948 :    
949 :     croak "MatrixReal1::det_LR(): not an LR decomposition matrix"
950 :     # unless ((defined $matrix->[3]) && ($rows == $cols));
951 :     unless ((defined $matrix->[4]) && ($rows == $cols)); #options might be in [3] position-- MEG
952 :    
953 :     $det = $matrix->[3]; # grab the sign from permutation shifts
954 :     for ( $k = 0; $k < $rows; $k++ )
955 :     {
956 :     $det *= $matrix->[0][$k][$k];
957 :     }
958 :     return($det);
959 :     }
960 :    
961 :     sub order_LR # order of LR decomposition matrix (number of non-zero equations)
962 :     {
963 :     croak "Usage: \$order = \$LR_matrix->order_LR();"
964 :     if (@_ != 1);
965 :    
966 :     my($matrix) = @_;
967 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
968 :     my($order);
969 :    
970 :     croak "MatrixReal1::order_LR(): not an LR decomposition matrix"
971 :     unless ((defined $matrix->[4]) && ($rows == $cols));
972 :    
973 :     ZERO:
974 :     for ( $order = ($rows - 1); $order >= 0; $order-- )
975 :     {
976 :     last ZERO if ($matrix->[0][$order][$order] != 0);
977 :     }
978 :     return(++$order);
979 :     }
980 :    
981 :     sub scalar_product
982 :     {
983 :     croak "Usage: \$scalar_product = \$vector1->scalar_product(\$vector2);"
984 :     if (@_ != 2);
985 :    
986 :     my($vector1,$vector2) = @_;
987 :     my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
988 :     my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
989 :     my($k,$sum);
990 :    
991 :     croak "MatrixReal1::scalar_product(): 1st vector is not a column vector"
992 :     unless ($cols1 == 1);
993 :    
994 :     croak "MatrixReal1::scalar_product(): 2nd vector is not a column vector"
995 :     unless ($cols2 == 1);
996 :    
997 :     croak "MatrixReal1::scalar_product(): vector size mismatch"
998 :     unless ($rows1 == $rows2);
999 :    
1000 :     $sum = 0;
1001 :     for ( $k = 0; $k < $rows1; $k++ )
1002 :     {
1003 :     $sum += $vector1->[0][$k][0] * $vector2->[0][$k][0];
1004 :     }
1005 :     return($sum);
1006 :     }
1007 :    
1008 :     sub vector_product
1009 :     {
1010 :     croak "Usage: \$vector_product = \$vector1->vector_product(\$vector2);"
1011 :     if (@_ != 2);
1012 :    
1013 :     my($vector1,$vector2) = @_;
1014 :     my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
1015 :     my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
1016 :     my($temp);
1017 :     my($n);
1018 :    
1019 :     croak "MatrixReal1::vector_product(): 1st vector is not a column vector"
1020 :     unless ($cols1 == 1);
1021 :    
1022 :     croak "MatrixReal1::vector_product(): 2nd vector is not a column vector"
1023 :     unless ($cols2 == 1);
1024 :    
1025 :     croak "MatrixReal1::vector_product(): vector size mismatch"
1026 :     unless ($rows1 == $rows2);
1027 :    
1028 :     $n = $rows1;
1029 :    
1030 :     croak "MatrixReal1::vector_product(): only defined for 3 dimensions"
1031 :     unless ($n == 3);
1032 :    
1033 :     $temp = $vector1->new($n,1);
1034 :     $temp->[0][0][0] = $vector1->[0][1][0] * $vector2->[0][2][0] -
1035 :     $vector1->[0][2][0] * $vector2->[0][1][0];
1036 :     $temp->[0][1][0] = $vector1->[0][2][0] * $vector2->[0][0][0] -
1037 :     $vector1->[0][0][0] * $vector2->[0][2][0];
1038 :     $temp->[0][2][0] = $vector1->[0][0][0] * $vector2->[0][1][0] -
1039 :     $vector1->[0][1][0] * $vector2->[0][0][0];
1040 :     return($temp);
1041 :     }
1042 :    
1043 :     sub length
1044 :     {
1045 :     croak "Usage: \$length = \$vector->length();"
1046 :     if (@_ != 1);
1047 :    
1048 :     my($vector) = @_;
1049 :     my($rows,$cols) = ($vector->[1],$vector->[2]);
1050 :     my($k,$comp,$sum);
1051 :    
1052 :     croak "MatrixReal1::length(): vector is not a column vector"
1053 :     unless ($cols == 1);
1054 :    
1055 :     $sum = 0;
1056 :     for ( $k = 0; $k < $rows; $k++ )
1057 :     {
1058 :     $comp = $vector->[0][$k][0];
1059 :     $sum += $comp * $comp;
1060 :     }
1061 :     return( sqrt( $sum ) );
1062 :     }
1063 :    
1064 :     sub _init_iteration
1065 :     {
1066 :     croak "Usage: \$which_norm = \$matrix->_init_iteration();"
1067 :     if (@_ != 1);
1068 :    
1069 :     my($matrix) = @_;
1070 :     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
1071 :     my($ok,$max,$sum,$norm);
1072 :     my($i,$j,$n);
1073 :    
1074 :     croak "MatrixReal1::_init_iteration(): matrix is not quadratic"
1075 :     unless ($rows == $cols);
1076 :    
1077 :     $ok = 1;
1078 :     $n = $rows;
1079 :     for ( $i = 0; $i < $n; $i++ )
1080 :     {
1081 :     if ($matrix->[0][$i][$i] == 0) { $ok = 0; }
1082 :     }
1083 :     if ($ok)
1084 :     {
1085 :     $norm = 1; # norm_one
1086 :     $max = 0;
1087 :     for ( $j = 0; $j < $n; $j++ )
1088 :     {
1089 :     $sum = 0;
1090 :     for ( $i = 0; $i < $j; $i++ )
1091 :     {
1092 :     $sum += abs($matrix->[0][$i][$j]);
1093 :     }
1094 :     for ( $i = ($j + 1); $i < $n; $i++ )
1095 :     {
1096 :     $sum += abs($matrix->[0][$i][$j]);
1097 :     }
1098 :     $sum /= abs($matrix->[0][$j][$j]);
1099 :     if ($sum > $max) { $max = $sum; }
1100 :     }
1101 :     $ok = ($max < 1);
1102 :     unless ($ok)
1103 :     {
1104 :     $norm = -1; # norm_max
1105 :     $max = 0;
1106 :     for ( $i = 0; $i < $n; $i++ )
1107 :     {
1108 :     $sum = 0;
1109 :     for ( $j = 0; $j < $i; $j++ )
1110 :     {
1111 :     $sum += abs($matrix->[0][$i][$j]);
1112 :     }
1113 :     for ( $j = ($i + 1); $j < $n; $j++ )
1114 :     {
1115 :     $sum += abs($matrix->[0][$i][$j]);
1116 :     }
1117 :     $sum /= abs($matrix->[0][$i][$i]);
1118 :     if ($sum > $max) { $max = $sum; }
1119 :     }
1120 :     $ok = ($max < 1)
1121 :     }
1122 :     }
1123 :     if ($ok) { return($norm); }
1124 :     else { return(0); }
1125 :     }
1126 :    
1127 :     sub solve_GSM # Global Step Method
1128 :     {
1129 :     croak "Usage: \$xn_vector = \$matrix->solve_GSM(\$x0_vector,\$b_vector,\$epsilon);"
1130 :     if (@_ != 4);
1131 :    
1132 :     my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
1133 :     my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1134 :     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1135 :     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1136 :     my($norm,$sum,$diff);
1137 :     my($xn_vector);
1138 :     my($i,$j,$n);
1139 :    
1140 :     croak "MatrixReal1::solve_GSM(): matrix is not quadratic"
1141 :     unless ($rows1 == $cols1);
1142 :    
1143 :     $n = $rows1;
1144 :    
1145 :     croak "MatrixReal1::solve_GSM(): 1st vector is not a column vector"
1146 :     unless ($cols2 == 1);
1147 :    
1148 :     croak "MatrixReal1::solve_GSM(): 2nd vector is not a column vector"
1149 :     unless ($cols3 == 1);
1150 :    
1151 :     croak "MatrixReal1::solve_GSM(): matrix and vector size mismatch"
1152 :     unless (($rows2 == $n) && ($rows3 == $n));
1153 :    
1154 :     return() unless ($norm = $matrix->_init_iteration());
1155 :    
1156 :     $xn_vector = $x0_vector->new($n,1);
1157 :    
1158 :     $diff = $epsilon + 1;
1159 :     while ($diff >= $epsilon)
1160 :     {
1161 :     for ( $i = 0; $i < $n; $i++ )
1162 :     {
1163 :     $sum = $b_vector->[0][$i][0];
1164 :     for ( $j = 0; $j < $i; $j++ )
1165 :     {
1166 :     $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
1167 :     }
1168 :     for ( $j = ($i + 1); $j < $n; $j++ )
1169 :     {
1170 :     $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
1171 :     }
1172 :     $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
1173 :     }
1174 :     $x0_vector->subtract($x0_vector,$xn_vector);
1175 :     if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1176 :     else { $diff = $x0_vector->norm_max(); }
1177 :     for ( $i = 0; $i < $n; $i++ )
1178 :     {
1179 :     $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1180 :     }
1181 :     }
1182 :     return($xn_vector);
1183 :     }
1184 :    
1185 :     sub solve_SSM # Single Step Method
1186 :     {
1187 :     croak "Usage: \$xn_vector = \$matrix->solve_SSM(\$x0_vector,\$b_vector,\$epsilon);"
1188 :     if (@_ != 4);
1189 :    
1190 :     my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
1191 :     my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1192 :     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1193 :     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1194 :     my($norm,$sum,$diff);
1195 :     my($xn_vector);
1196 :     my($i,$j,$n);
1197 :    
1198 :     croak "MatrixReal1::solve_SSM(): matrix is not quadratic"
1199 :     unless ($rows1 == $cols1);
1200 :    
1201 :     $n = $rows1;
1202 :    
1203 :     croak "MatrixReal1::solve_SSM(): 1st vector is not a column vector"
1204 :     unless ($cols2 == 1);
1205 :    
1206 :     croak "MatrixReal1::solve_SSM(): 2nd vector is not a column vector"
1207 :     unless ($cols3 == 1);
1208 :    
1209 :     croak "MatrixReal1::solve_SSM(): matrix and vector size mismatch"
1210 :     unless (($rows2 == $n) && ($rows3 == $n));
1211 :    
1212 :     return() unless ($norm = $matrix->_init_iteration());
1213 :    
1214 :     $xn_vector = $x0_vector->new($n,1);
1215 :     $xn_vector->copy($x0_vector);
1216 :    
1217 :     $diff = $epsilon + 1;
1218 :     while ($diff >= $epsilon)
1219 :     {
1220 :     for ( $i = 0; $i < $n; $i++ )
1221 :     {
1222 :     $sum = $b_vector->[0][$i][0];
1223 :     for ( $j = 0; $j < $i; $j++ )
1224 :     {
1225 :     $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1226 :     }
1227 :     for ( $j = ($i + 1); $j < $n; $j++ )
1228 :     {
1229 :     $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1230 :     }
1231 :     $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
1232 :     }
1233 :     $x0_vector->subtract($x0_vector,$xn_vector);
1234 :     if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1235 :     else { $diff = $x0_vector->norm_max(); }
1236 :     for ( $i = 0; $i < $n; $i++ )
1237 :     {
1238 :     $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1239 :     }
1240 :     }
1241 :     return($xn_vector);
1242 :     }
1243 :    
1244 :     sub solve_RM # Relaxation Method
1245 :     {
1246 :     croak "Usage: \$xn_vector = \$matrix->solve_RM(\$x0_vector,\$b_vector,\$weight,\$epsilon);"
1247 :     if (@_ != 5);
1248 :    
1249 :     my($matrix,$x0_vector,$b_vector,$weight,$epsilon) = @_;
1250 :     my($rows1,$cols1) = ( $matrix->[1], $matrix->[2]);
1251 :     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
1252 :     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
1253 :     my($norm,$sum,$diff);
1254 :     my($xn_vector);
1255 :     my($i,$j,$n);
1256 :    
1257 :     croak "MatrixReal1::solve_RM(): matrix is not quadratic"
1258 :     unless ($rows1 == $cols1);
1259 :    
1260 :     $n = $rows1;
1261 :    
1262 :     croak "MatrixReal1::solve_RM(): 1st vector is not a column vector"
1263 :     unless ($cols2 == 1);
1264 :    
1265 :     croak "MatrixReal1::solve_RM(): 2nd vector is not a column vector"
1266 :     unless ($cols3 == 1);
1267 :    
1268 :     croak "MatrixReal1::solve_RM(): matrix and vector size mismatch"
1269 :     unless (($rows2 == $n) && ($rows3 == $n));
1270 :    
1271 :     return() unless ($norm = $matrix->_init_iteration());
1272 :    
1273 :     $xn_vector = $x0_vector->new($n,1);
1274 :     $xn_vector->copy($x0_vector);
1275 :    
1276 :     $diff = $epsilon + 1;
1277 :     while ($diff >= $epsilon)
1278 :     {
1279 :     for ( $i = 0; $i < $n; $i++ )
1280 :     {
1281 :     $sum = $b_vector->[0][$i][0];
1282 :     for ( $j = 0; $j < $i; $j++ )
1283 :     {
1284 :     $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1285 :     }
1286 :     for ( $j = ($i + 1); $j < $n; $j++ )
1287 :     {
1288 :     $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
1289 :     }
1290 :     $xn_vector->[0][$i][0] = $weight * ( $sum / $matrix->[0][$i][$i] )
1291 :     + (1 - $weight) * $xn_vector->[0][$i][0];
1292 :     }
1293 :     $x0_vector->subtract($x0_vector,$xn_vector);
1294 :     if ($norm > 0) { $diff = $x0_vector->norm_one(); }
1295 :     else { $diff = $x0_vector->norm_max(); }
1296 :     for ( $i = 0; $i < $n; $i++ )
1297 :     {
1298 :     $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
1299 :     }
1300 :     }
1301 :     return($xn_vector);
1302 :     }
1303 :    
1304 :     # Core householder reduction routine (when eagenvector
1305 :     # are wanted).
1306 :     # Adapted from: Numerical Recipes, 2nd edition.
1307 :     sub _householder_vectors ($)
1308 :     {
1309 :     my ($Q) = @_;
1310 :     my ($rows, $cols) = ($Q->[1], $Q->[2]);
1311 :    
1312 :     # Creates tridiagonal
1313 :     # Set up tridiagonal needed elements
1314 :     my $d = []; # N Diagonal elements 0...n-1
1315 :     my $e = []; # N-1 Off-Diagonal elements 0...n-2
1316 :    
1317 :     my @p = ();
1318 :     for (my $i = ($rows-1); $i > 1; $i--)
1319 :     {
1320 :     my $scale = 0.0;
1321 :     # Computes norm of one column (below diagonal)
1322 :     for (my $k = 0; $k < $i; $k++)
1323 :     {
1324 :     $scale += abs($Q->[0][$i][$k]);
1325 :     }
1326 :     if ($scale == 0.0)
1327 :     { # skip the transformation
1328 :     $e->[$i-1] = $Q->[0][$i][$i-1];
1329 :     }
1330 :     else
1331 :     {
1332 :     my $h = 0.0;
1333 :     for (my $k = 0; $k < $i; $k++)
1334 :     { # Used scaled Q for transformation
1335 :     $Q->[0][$i][$k] /= $scale;
1336 :     # Form sigma in h
1337 :     $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
1338 :     }
1339 :     my $t1 = $Q->[0][$i][$i-1];
1340 :     my $t2 = (($t1 >= 0.0) ? -sqrt($h) : sqrt($h));
1341 :     $e->[$i-1] = $scale * $t2; # Update off-diagonals
1342 :     $h -= $t1 * $t2;
1343 :     $Q->[0][$i][$i-1] -= $t2;
1344 :     my $f = 0.0;
1345 :     for (my $j = 0; $j < $i; $j++)
1346 :     {
1347 :     $Q->[0][$j][$i] = $Q->[0][$i][$j] / $h;
1348 :     my $g = 0.0;
1349 :     for (my $k = 0; $k <= $j; $k++)
1350 :     {
1351 :     $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
1352 :     }
1353 :     for (my $k = $j+1; $k < $i; $k++)
1354 :     {
1355 :     $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
1356 :     }
1357 :     # Form elements of P
1358 :     $p[$j] = $g / $h;
1359 :     $f += $p[$j] * $Q->[0][$i][$j];
1360 :     }
1361 :     my $hh = $f / ($h + $h);
1362 :     for (my $j = 0; $j < $i; $j++)
1363 :     {
1364 :     my $t3 = $Q->[0][$i][$j];
1365 :     my $t4 = $p[$j] - $hh * $t3;
1366 :     $p[$j] = $t4;
1367 :     for (my $k = 0; $k <= $j; $k++)
1368 :     {
1369 :     $Q->[0][$j][$k] -= $t3 * $p[$k]
1370 :     + $t4 * $Q->[0][$i][$k];
1371 :     }
1372 :     }
1373 :     }
1374 :     }
1375 :     # Updates for i == 0,1
1376 :     $e->[0] = $Q->[0][1][0];
1377 :     $d->[0] = $Q->[0][0][0]; # i==0
1378 :     $Q->[0][0][0] = 1.0;
1379 :     $d->[1] = $Q->[0][1][1]; # i==1
1380 :     $Q->[0][1][1] = 1.0;
1381 :     $Q->[0][1][0] = $Q->[0][0][1] = 0.0;
1382 :     for (my $i = 2; $i < $rows; $i++)
1383 :     {
1384 :     for (my $j = 0; $j < $i; $j++)
1385 :     {
1386 :     my $g = 0.0;
1387 :     for (my $k = 0; $k < $i; $k++)
1388 :     {
1389 :     $g += $Q->[0][$i][$k] * $Q->[0][$k][$j];
1390 :     }
1391 :     for (my $k = 0; $k < $i; $k++)
1392 :     {
1393 :     $Q->[0][$k][$j] -= $g * $Q->[0][$k][$i];
1394 :     }
1395 :     }
1396 :     $d->[$i] = $Q->[0][$i][$i];
1397 :     # Reset row and column of Q for next iteration
1398 :     $Q->[0][$i][$i] = 1.0;
1399 :     for (my $j = 0; $j < $i; $j++)
1400 :     {
1401 :     $Q->[0][$i][$j] = $Q->[0][$j][$i] = 0.0;
1402 :     }
1403 :     }
1404 :     return ($d, $e);
1405 :     }
1406 :    
1407 :     # Computes sqrt(a*a + b*b), but more carefully...
1408 :     sub _pythag ($$)
1409 :     {
1410 :     my ($a, $b) = @_;
1411 :     my $aa = abs($a);
1412 :     my $ab = abs($b);
1413 :     if ($aa > $ab)
1414 :     {
1415 :     # NB: Not needed!: return 0.0 if ($aa == 0.0);
1416 :     my $t = $ab / $aa;
1417 :     return ($aa * sqrt(1.0 + $t*$t));
1418 :     }
1419 :     else
1420 :     {
1421 :     return 0.0 if ($ab == 0.0);
1422 :     my $t = $aa / $ab;
1423 :     return ($ab * sqrt(1.0 + $t*$t));
1424 :     }
1425 :     }
1426 :    
1427 :     # QL algorithm with implicit shifts to determine the eigenvalues
1428 :     # of a tridiagonal matrix. Internal routine.
1429 :     sub _tridiagonal_QLimplicit
1430 :     {
1431 :     my ($EV, $d, $e) = @_;
1432 :     my ($rows, $cols) = ($EV->[1], $EV->[2]);
1433 :    
1434 :     $e->[$rows-1] = 0.0;
1435 :     # Start real computation
1436 :     for (my $l = 0; $l < $rows; $l++)
1437 :     {
1438 :     my $iter = 0;
1439 :     my $m;
1440 :     OUTER:
1441 :     do {
1442 :     for ($m = $l; $m < ($rows - 1); $m++)
1443 :     {
1444 :     my $dd = abs($d->[$m]) + abs($d->[$m+1]);
1445 :     last if ((abs($e->[$m]) + $dd) == $dd);
1446 :     }
1447 :     if ($m != $l)
1448 :     {
1449 :     croak("Too many iterations!") if ($iter++ >= 30);
1450 :     my $g = ($d->[$l+1] - $d->[$l])
1451 :     / (2.0 * $e->[$l]);
1452 :     my $r = _pythag($g, 1.0);
1453 :     $g = $d->[$m] - $d->[$l]
1454 :     + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
1455 :     my ($p,$s,$c) = (0.0, 1.0,1.0);
1456 :     for (my $i = ($m-1); $i >= $l; $i--)
1457 :     {
1458 :     my $ii = $i + 1;
1459 :     my $f = $s * $e->[$i];
1460 :     my $t = _pythag($f, $g);
1461 :     $e->[$ii] = $t;
1462 :     if ($t == 0.0)
1463 :     {
1464 :     $d->[$ii] -= $p;
1465 :     $e->[$m] = 0.0;
1466 :     next OUTER;
1467 :     }
1468 :     my $b = $c * $e->[$i];
1469 :     $s = $f / $t;
1470 :     $c = $g / $t;
1471 :     $g = $d->[$ii] - $p;
1472 :     my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
1473 :     $p = $s * $t2;
1474 :     $d->[$ii] = $g + $p;
1475 :     $g = $c * $t2 - $b;
1476 :     for (my $k = 0; $k < $rows; $k++)
1477 :     {
1478 :     my $t1 = $EV->[0][$k][$ii];
1479 :     my $t2 = $EV->[0][$k][$i];
1480 :     $EV->[0][$k][$ii] = $s * $t2 + $c * $t1;
1481 :     $EV->[0][$k][$i] = $c * $t2 - $s * $t1;
1482 :     }
1483 :     }
1484 :     $d->[$l] -= $p;
1485 :     $e->[$l] = $g;
1486 :     $e->[$m] = 0.0;
1487 :     }
1488 :     } while ($m != $l);
1489 :     }
1490 :     return;
1491 :     }
1492 :    
1493 :     # Core householder reduction routine (when eagenvector
1494 :     # are NOT wanted).
1495 :     sub _householder_values ($)
1496 :     {
1497 :     my ($Q) = @_; # NB: Q is destroyed on output...
1498 :     my ($rows, $cols) = ($Q->[1], $Q->[2]);
1499 :    
1500 :     # Creates tridiagonal
1501 :     # Set up tridiagonal needed elements
1502 :     my $d = []; # N Diagonal elements 0...n-1
1503 :     my $e = []; # N-1 Off-Diagonal elements 0...n-2
1504 :    
1505 :     my @p = ();
1506 :     for (my $i = ($rows - 1); $i > 1; $i--)
1507 :     {
1508 :     my $scale = 0.0;
1509 :     for (my $k = 0; $k < $i; $k++)
1510 :     {
1511 :     $scale += abs($Q->[0][$i][$k]);
1512 :     }
1513 :     if ($scale == 0.0)
1514 :     { # skip the transformation
1515 :     $e->[$i-1] = $Q->[0][$i][$i-1];
1516 :     }
1517 :     else
1518 :     {
1519 :     my $h = 0.0;
1520 :     for (my $k = 0; $k < $i; $k++)
1521 :     { # Used scaled Q for transformation
1522 :     $Q->[0][$i][$k] /= $scale;
1523 :     # Form sigma in h
1524 :     $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
1525 :     }
1526 :     my $t = $Q->[0][$i][$i-1];
1527 :     my $t2 = (($t >= 0.0) ? -sqrt($h) : sqrt($h));
1528 :     $e->[$i-1] = $scale * $t2; # Updates off-diagonal
1529 :     $h -= $t * $t2;
1530 :     $Q->[0][$i][$i-1] -= $t2;
1531 :     my $f = 0.0;
1532 :     for (my $j = 0; $j < $i; $j++)
1533 :     {
1534 :     my $g = 0.0;
1535 :     for (my $k = 0; $k <= $j; $k++)
1536 :     {
1537 :     $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
1538 :     }
1539 :     for (my $k = $j+1; $k < $i; $k++)
1540 :     {
1541 :     $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
1542 :     }
1543 :     # Form elements of P
1544 :     $p[$j] = $g / $h;
1545 :     $f += $p[$j] * $Q->[0][$i][$j];
1546 :     }
1547 :     my $hh = $f / ($h + $h);
1548 :     for (my $j = 0; $j < $i; $j++)
1549 :     {
1550 :     my $t = $Q->[0][$i][$j];
1551 :     my $g = $p[$j] - $hh * $t;
1552 :     $p[$j] = $g;
1553 :     for (my $k = 0; $k <= $j; $k++)
1554 :     {
1555 :     $Q->[0][$j][$k] -= $t * $p[$k]
1556 :     + $g * $Q->[0][$i][$k];
1557 :     }
1558 :     }
1559 :     }
1560 :     }
1561 :     # Updates for i==1
1562 :     $e->[0] = $Q->[0][1][0];
1563 :     # Updates diagonal elements
1564 :     for (my $i = 0; $i < $rows; $i++)
1565 :     {
1566 :     $d->[$i] = $Q->[0][$i][$i];
1567 :     }
1568 :     return ($d, $e);
1569 :     }
1570 :    
1571 :     # QL algorithm with implicit shifts to determine the
1572 :     # eigenvalues ONLY. This is O(N^2) only...
1573 :     sub _tridiagonal_QLimplicit_values
1574 :     {
1575 :     my ($M, $d, $e) = @_; # NB: M is not touched...
1576 :     my ($rows, $cols) = ($M->[1], $M->[2]);
1577 :    
1578 :     $e->[$rows-1] = 0.0;
1579 :     # Start real computation
1580 :     for (my $l = 0; $l < $rows; $l++)
1581 :     {
1582 :     my $iter = 0;
1583 :     my $m;
1584 :     OUTER:
1585 :     do {
1586 :     for ($m = $l; $m < ($rows - 1); $m++)
1587 :     {
1588 :     my $dd = abs($d->[$m]) + abs($d->[$m+1]);
1589 :     last if ((abs($e->[$m]) + $dd) == $dd);
1590 :     }
1591 :     if ($m != $l)
1592 :     {
1593 :     croak("Too many iterations!") if ($iter++ >= 30);
1594 :     my $g = ($d->[$l+1] - $d->[$l])
1595 :     / (2.0 * $e->[$l]);
1596 :     my $r = _pythag($g, 1.0);
1597 :     $g = $d->[$m] - $d->[$l]
1598 :     + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
1599 :     my ($p,$s,$c) = (0.0, 1.0,1.0);
1600 :     for (my $i = ($m-1); $i >= $l; $i--)
1601 :     {
1602 :     my $ii = $i + 1;
1603 :     my $f = $s * $e->[$i];
1604 :     my $t = _pythag($f, $g);
1605 :     $e->[$ii] = $t;
1606 :     if ($t == 0.0)
1607 :     {
1608 :     $d->[$ii] -= $p;
1609 :     $e->[$m] = 0.0;
1610 :     next OUTER;
1611 :     }
1612 :     my $b = $c * $e->[$i];
1613 :     $s = $f / $t;
1614 :     $c = $g / $t;
1615 :     $g = $d->[$ii] - $p;
1616 :     my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
1617 :     $p = $s * $t2;
1618 :     $d->[$ii] = $g + $p;
1619 :     $g = $c * $t2 - $b;
1620 :     }
1621 :     $d->[$l] -= $p;
1622 :     $e->[$l] = $g;
1623 :     $e->[$m] = 0.0;
1624 :     }
1625 :     } while ($m != $l);
1626 :     }
1627 :     return;
1628 :     }
1629 :    
1630 :     # Householder reduction of a real, symmetric matrix A.
1631 :     # Returns a tridiagonal matrix T and the orthogonal matrix
1632 :     # Q effecting the transformation between A and T.
1633 :     sub householder ($)
1634 :     {
1635 :     my ($A) = @_;
1636 :     my ($rows, $cols) = ($A->[1], $A->[2]);
1637 :    
1638 :     croak "Matrix is not quadratic"
1639 :     unless ($rows = $cols);
1640 :     croak "Matrix is not symmetric"
1641 :     unless ($A->is_symmetric());
1642 :    
1643 :     # Copy given matrix TODO: study if we should do in-place modification
1644 :     my $Q = $A->clone();
1645 :    
1646 :     # Do the computation of tridiagonal elements and of
1647 :     # transformation matrix
1648 :     my ($diag, $offdiag) = $Q->_householder_vectors();
1649 :    
1650 :     # Creates the tridiagonal matrix
1651 :     my $T = $A->shadow();
1652 :     for (my $i = 0; $i < $rows; $i++)
1653 :     { # Set diagonal
1654 :     $T->[0][$i][$i] = $diag->[$i];
1655 :     }
1656 :     for (my $i = 0; $i < ($rows-1); $i++)
1657 :     { # Set off diagonals
1658 :     $T->[0][$i+1][$i] = $offdiag->[$i];
1659 :     $T->[0][$i][$i+1] = $offdiag->[$i];
1660 :     }
1661 :     return ($T, $Q);
1662 :     }
1663 :    
1664 :     # QL algorithm with implicit shifts to determine the eigenvalues
1665 :     # and eigenvectors of a real tridiagonal matrix - or of a matrix
1666 :     # previously reduced to tridiagonal form.
1667 :     sub tri_diagonalize ($;$)
1668 :     {
1669 :     my ($T,$Q) = @_; # Q may be 0 if the original matrix is really tridiagonal
1670 :    
1671 :     my ($rows, $cols) = ($T->[1], $T->[2]);
1672 :    
1673 :     croak "Matrix is not quadratic"
1674 :     unless ($rows = $cols);
1675 :     croak "Matrix is not tridiagonal"
1676 :     unless ($T->is_symmetric()); # TODO: Do real tridiag check (not symmetric)!
1677 :    
1678 :     my $EV;
1679 :     # Obtain/Creates the todo eigenvectors matrix
1680 :     if ($Q)
1681 :     {
1682 :     $EV = $Q->clone();
1683 :     }
1684 :     else
1685 :     {
1686 :     $EV = $T->shadow();
1687 :     $EV->one();
1688 :     }
1689 :     # Allocates diagonal vector
1690 :     my $diag = [ ];
1691 :     # Initializes it with T
1692 :     for (my $i = 0; $i < $rows; $i++)
1693 :     {
1694 :     $diag->[$i] = $T->[0][$i][$i];
1695 :     }
1696 :     # Allocate temporary vector for off-diagonal elements
1697 :     my $offdiag = [ ];
1698 :     for (my $i = 1; $i < $rows; $i++)
1699 :     {
1700 :     $offdiag->[$i-1] = $T->[0][$i][$i-1];
1701 :     }
1702 :    
1703 :     # Calls the calculus routine
1704 :     $EV->_tridiagonal_QLimplicit($diag, $offdiag);
1705 :    
1706 :     # Allocate eigenvalues vector
1707 :     my $v = MatrixReal1->new($rows,1);
1708 :     # Fills it
1709 :     for (my $i = 0; $i < $rows; $i++)
1710 :     {
1711 :     $v->[0][$i][0] = $diag->[$i];
1712 :     }
1713 :     return ($v, $EV);
1714 :     }
1715 :    
1716 :     # Main routine for diagonalization of a real symmetric
1717 :     # matrix M. Operates by transforming M into a tridiagonal
1718 :     # matrix and then obtaining the eigenvalues and eigenvectors
1719 :     # for that matrix (taking into account the transformation to
1720 :     # tridiagonal).
1721 :     sub sym_diagonalize ($)
1722 :     {
1723 :     my ($M) = @_;
1724 :     my ($rows, $cols) = ($M->[1], $M->[2]);
1725 :    
1726 :     croak "Matrix is not quadratic"
1727 :     unless ($rows = $cols);
1728 :     croak "Matrix is not symmetric"
1729 :     unless ($M->is_symmetric());
1730 :    
1731 :     # Copy initial matrix
1732 :     # TODO: study if we should allow in-place modification
1733 :     my $VEC = $M->clone();
1734 :    
1735 :     # Do the computation of tridiagonal elements and of
1736 :     # transformation matrix
1737 :     my ($diag, $offdiag) = $VEC->_householder_vectors();
1738 :     # Calls the calculus routine for diagonalization
1739 :     $VEC->_tridiagonal_QLimplicit($diag, $offdiag);
1740 :    
1741 :     # Allocate eigenvalues vector
1742 :     my $val = MatrixReal1->new($rows,1);
1743 :     # Fills it
1744 :     for (my $i = 0; $i < $rows; $i++)
1745 :     {
1746 :     $val->[0][$i][0] = $diag->[$i];
1747 :     }
1748 :     return ($val, $VEC);
1749 :     }
1750 :    
1751 :     # Householder reduction of a real, symmetric matrix A.
1752 :     # Returns a tridiagonal matrix T equivalent to A.
1753 :     sub householder_tridiagonal ($)
1754 :     {
1755 :     my ($A) = @_;
1756 :     my ($rows, $cols) = ($A->[1], $A->[2]);
1757 :    
1758 :     croak "Matrix is not quadratic"
1759 :     unless ($rows = $cols);
1760 :     croak "Matrix is not symmetric"
1761 :     unless ($A->is_symmetric());
1762 :    
1763 :     # Copy given matrix
1764 :     my $Q = $A->clone();
1765 :    
1766 :     # Do the computation of tridiagonal elements and of
1767 :     # transformation matrix
1768 :     # Q is destroyed after reduction
1769 :     my ($diag, $offdiag) = $Q->_householder_values();
1770 :    
1771 :     # Creates the tridiagonal matrix in Q (avoid allocation)
1772 :     my $T = $Q;
1773 :     $T->zero();
1774 :     for (my $i = 0; $i < $rows; $i++)
1775 :     { # Set diagonal
1776 :     $T->[0][$i][$i] = $diag->[$i];
1777 :     }
1778 :     for (my $i = 0; $i < ($rows-1); $i++)
1779 :     { # Set off diagonals
1780 :     $T->[0][$i+1][$i] = $offdiag->[$i];
1781 :     $T->[0][$i][$i+1] = $offdiag->[$i];
1782 :     }
1783 :     return $T;
1784 :     }
1785 :    
1786 :     # QL algorithm with implicit shifts to determine ONLY
1787 :     # the eigenvalues a real tridiagonal matrix - or of a
1788 :     # matrix previously reduced to tridiagonal form.
1789 :     sub tri_eigenvalues ($;$)
1790 :     {
1791 :     my ($T) = @_;
1792 :     my ($rows, $cols) = ($T->[1], $T->[2]);
1793 :    
1794 :     croak "Matrix is not quadratic"
1795 :     unless ($rows = $cols);
1796 :     croak "Matrix is not tridiagonal"
1797 :     unless ($T->is_symmetric()); # TODO: Do real tridiag check (not symmetric)!
1798 :    
1799 :     # Allocates diagonal vector
1800 :     my $diag = [ ];
1801 :     # Initializes it with T
1802 :     for (my $i = 0; $i < $rows; $i++)
1803 :     {
1804 :     $diag->[$i] = $T->[0][$i][$i];
1805 :     }
1806 :     # Allocate temporary vector for off-diagonal elements
1807 :     my $offdiag = [ ];
1808 :     for (my $i = 1; $i < $rows; $i++)
1809 :     {
1810 :     $offdiag->[$i-1] = $T->[0][$i][$i-1];
1811 :     }
1812 :    
1813 :     # Calls the calculus routine (T is not touched)
1814 :     $T->_tridiagonal_QLimplicit_values($diag, $offdiag);
1815 :    
1816 :     # Allocate eigenvalues vector
1817 :     my $v = MatrixReal1->new($rows,1);
1818 :     # Fills it
1819 :     for (my $i = 0; $i < $rows; $i++)
1820 :     {
1821 :     $v->[0][$i][0] = $diag->[$i];
1822 :     }
1823 :     return $v;
1824 :     }
1825 :    
1826 :     # Main routine for diagonalization of a real symmetric
1827 :     # matrix M. Operates by transforming M into a tridiagonal
1828 :     # matrix and then obtaining the eigenvalues and eigenvectors
1829 :     # for that matrix (taking into account the transformation to
1830 :     # tridiagonal).
1831 :     sub sym_eigenvalues ($)
1832 :     {
1833 :     my ($M) = @_;
1834 :     my ($rows, $cols) = ($M->[1], $M->[2]);
1835 :    
1836 :     croak "Matrix is not quadratic"
1837 :     unless ($rows = $cols);
1838 :     croak "Matrix is not symmetric"
1839 :     unless ($M->is_symmetric());
1840 :    
1841 :     # Copy matrix in temporary
1842 :     my $A = $M->clone();
1843 :     # Do the computation of tridiagonal elements and of
1844 :     # transformation matrix. A is destroyed
1845 :     my ($diag, $offdiag) = $A->_householder_values();
1846 :     # Calls the calculus routine for diagonalization
1847 :     # (M is not touched)
1848 :     $M->_tridiagonal_QLimplicit_values($diag, $offdiag);
1849 :    
1850 :     # Allocate eigenvalues vector
1851 :     my $val = MatrixReal1->new($rows,1);
1852 :     # Fills it
1853 :     for (my $i = 0; $i < $rows; $i++)
1854 :     {
1855 :     $val->[0][$i][0] = $diag->[$i];
1856 :     }
1857 :     return $val;
1858 :     }
1859 :    
1860 :     # Boolean check routine to see if a matrix is
1861 :     # symmetric
1862 :     sub is_symmetric ($)
1863 :     {
1864 :     my ($M) = @_;
1865 :     my ($rows, $cols) = ($M->[1], $M->[2]);
1866 :     # if it is not quadratic it cannot be symmetric...
1867 :     return 0 unless ($rows == $cols);
1868 :     for (my $i = 1; $i < $rows; $i++)
1869 :     {
1870 :     for (my $j = 0; $j < $i; $j++)
1871 :     {
1872 :     return 0 unless ($M->[0][$i][$j] == $M->[0][$j][$i]);
1873 :     }
1874 :     }
1875 :     return 1;
1876 :     }
1877 :    
1878 :     ########################################
1879 :     # #
1880 :     # define overloaded operators section: #
1881 :     # #
1882 :     ########################################
1883 :    
1884 :     sub _negate
1885 :     {
1886 :     my($object,$argument,$flag) = @_;
1887 :     # my($name) = "neg"; #&_trace($name,$object,$argument,$flag);
1888 :     my($temp);
1889 :    
1890 :     $temp = $object->new($object->[1],$object->[2]);
1891 :     $temp->negate($object);
1892 :     return($temp);
1893 :     }
1894 :    
1895 :     sub _transpose
1896 :     {
1897 :     my($object,$argument,$flag) = @_;
1898 :     # my($name) = "'~'"; #&_trace($name,$object,$argument,$flag);
1899 :     my($temp);
1900 :    
1901 :     $temp = $object->new($object->[2],$object->[1]);
1902 :     $temp->transpose($object);
1903 :     return($temp);
1904 :     }
1905 :    
1906 :     sub _boolean
1907 :     {
1908 :     my($object,$argument,$flag) = @_;
1909 :     # my($name) = "bool"; #&_trace($name,$object,$argument,$flag);
1910 :     my($rows,$cols) = ($object->[1],$object->[2]);
1911 :     my($i,$j,$result);
1912 :    
1913 :     $result = 0;
1914 :     BOOL:
1915 :     for ( $i = 0; $i < $rows; $i++ )
1916 :     {
1917 :     for ( $j = 0; $j < $cols; $j++ )
1918 :     {
1919 :     if ($object->[0][$i][$j] != 0)
1920 :     {
1921 :     $result = 1;
1922 :     last BOOL;
1923 :     }
1924 :     }
1925 :     }
1926 :     return($result);
1927 :     }
1928 :    
1929 :     sub _not_boolean
1930 :     {
1931 :     my($object,$argument,$flag) = @_;
1932 :     # my($name) = "'!'"; #&_trace($name,$object,$argument,$flag);
1933 :     my($rows,$cols) = ($object->[1],$object->[2]);
1934 :     my($i,$j,$result);
1935 :    
1936 :     $result = 1;
1937 :     NOTBOOL:
1938 :     for ( $i = 0; $i < $rows; $i++ )
1939 :     {
1940 :     for ( $j = 0; $j < $cols; $j++ )
1941 :     {
1942 :     if ($object->[0][$i][$j] != 0)
1943 :     {
1944 :     $result = 0;
1945 :     last NOTBOOL;
1946 :     }
1947 :     }
1948 :     }
1949 :     return($result);
1950 :     }
1951 :    
1952 :     sub _stringify
1953 :     {
1954 :     my($object,$argument,$flag) = @_;
1955 :     # my($name) = '""'; #&_trace($name,$object,$argument,$flag);
1956 :     my($rows,$cols) = ($object->[1],$object->[2]);
1957 :     my($i,$j,$s);
1958 :    
1959 :     $s = '';
1960 :     for ( $i = 0; $i < $rows; $i++ )
1961 :     {
1962 :     $s .= "[ ";
1963 :     for ( $j = 0; $j < $cols; $j++ )
1964 :     {
1965 :     $s .= sprintf("% #-19.12E ", $object->[0][$i][$j]);
1966 :     }
1967 :     $s .= "]\n";
1968 :     }
1969 :     return($s);
1970 :     }
1971 :    
1972 :     sub _norm
1973 :     {
1974 :     my($object,$argument,$flag) = @_;
1975 :     # my($name) = "abs"; #&_trace($name,$object,$argument,$flag);
1976 :    
1977 :     return( $object->norm_one() );
1978 :     }
1979 :    
1980 :     sub _add
1981 :     {
1982 :     my($object,$argument,$flag) = @_;
1983 :     my($name) = "'+'"; #&_trace($name,$object,$argument,$flag);
1984 :     my($temp);
1985 :    
1986 :     if ((defined $argument) && ref($argument) &&
1987 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
1988 :     {
1989 :     if (defined $flag)
1990 :     {
1991 :     $temp = $object->new($object->[1],$object->[2]);
1992 :     $temp->add($object,$argument);
1993 :     return($temp);
1994 :     }
1995 :     else
1996 :     {
1997 :     $object->add($object,$argument);
1998 :     return($object);
1999 :     }
2000 :     }
2001 :     else
2002 :     {
2003 :     croak "MatrixReal1 $name: wrong argument type";
2004 :     }
2005 :     }
2006 :    
2007 :     sub _subtract
2008 :     {
2009 :     my($object,$argument,$flag) = @_;
2010 :     my($name) = "'-'"; #&_trace($name,$object,$argument,$flag);
2011 :     my($temp);
2012 :    
2013 :     if ((defined $argument) && ref($argument) &&
2014 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2015 :     {
2016 :     if (defined $flag)
2017 :     {
2018 :     $temp = $object->new($object->[1],$object->[2]);
2019 :     if ($flag) { $temp->subtract($argument,$object); }
2020 :     else { $temp->subtract($object,$argument); }
2021 :     return($temp);
2022 :     }
2023 :     else
2024 :     {
2025 :     $object->subtract($object,$argument);
2026 :     return($object);
2027 :     }
2028 :     }
2029 :     else
2030 :     {
2031 :     croak "MatrixReal1 $name: wrong argument type";
2032 :     }
2033 :     }
2034 :    
2035 :     sub _multiply
2036 :     {
2037 :     my($object,$argument,$flag) = @_;
2038 :     my($name) = "'*'"; #&_trace($name,$object,$argument,$flag);
2039 :     my($temp);
2040 :    
2041 :     if ((defined $argument) && ref($argument) &&
2042 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2043 :     {
2044 :     if ((defined $flag) && $flag)
2045 :     {
2046 :     return( multiply($argument,$object) );
2047 :     }
2048 :     else
2049 :     {
2050 :     return( multiply($object,$argument) );
2051 :     }
2052 :     }
2053 :     elsif ((defined $argument) && !(ref($argument)))
2054 :     {
2055 :     if (defined $flag)
2056 :     {
2057 :     $temp = $object->new($object->[1],$object->[2]);
2058 :     $temp->multiply_scalar($object,$argument);
2059 :     return($temp);
2060 :     }
2061 :     else
2062 :     {
2063 :     $object->multiply_scalar($object,$argument);
2064 :     return($object);
2065 :     }
2066 :     }
2067 :     else
2068 :     {
2069 :     croak "MatrixReal1 $name: wrong argument type";
2070 :     }
2071 :     }
2072 :    
2073 :     sub _assign_add
2074 :     {
2075 :     my($object,$argument,$flag) = @_;
2076 :     # my($name) = "'+='"; #&_trace($name,$object,$argument,$flag);
2077 :    
2078 :     return( &_add($object,$argument,undef) );
2079 :     }
2080 :    
2081 :     sub _assign_subtract
2082 :     {
2083 :     my($object,$argument,$flag) = @_;
2084 :     # my($name) = "'-='"; #&_trace($name,$object,$argument,$flag);
2085 :    
2086 :     return( &_subtract($object,$argument,undef) );
2087 :     }
2088 :    
2089 :     sub _assign_multiply
2090 :     {
2091 :     my($object,$argument,$flag) = @_;
2092 :     # my($name) = "'*='"; #&_trace($name,$object,$argument,$flag);
2093 :    
2094 :     return( &_multiply($object,$argument,undef) );
2095 :     }
2096 :    
2097 :     sub _equal
2098 :     {
2099 :     my($object,$argument,$flag) = @_;
2100 :     my($name) = "'=='"; #&_trace($name,$object,$argument,$flag);
2101 :     my($rows,$cols) = ($object->[1],$object->[2]);
2102 :     my($i,$j,$result);
2103 :    
2104 :     if ((defined $argument) && ref($argument) &&
2105 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2106 :     {
2107 :     $result = 1;
2108 :     EQUAL:
2109 :     for ( $i = 0; $i < $rows; $i++ )
2110 :     {
2111 :     for ( $j = 0; $j < $cols; $j++ )
2112 :     {
2113 :     if ($object->[0][$i][$j] != $argument->[0][$i][$j])
2114 :     {
2115 :     $result = 0;
2116 :     last EQUAL;
2117 :     }
2118 :     }
2119 :     }
2120 :     return($result);
2121 :     }
2122 :     else
2123 :     {
2124 :     croak "MatrixReal1 $name: wrong argument type";
2125 :     }
2126 :     }
2127 :    
2128 :     sub _not_equal
2129 :     {
2130 :     my($object,$argument,$flag) = @_;
2131 :     my($name) = "'!='"; #&_trace($name,$object,$argument,$flag);
2132 :     my($rows,$cols) = ($object->[1],$object->[2]);
2133 :     my($i,$j,$result);
2134 :    
2135 :     if ((defined $argument) && ref($argument) &&
2136 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2137 :     {
2138 :     $result = 0;
2139 :     NOTEQUAL:
2140 :     for ( $i = 0; $i < $rows; $i++ )
2141 :     {
2142 :     for ( $j = 0; $j < $cols; $j++ )
2143 :     {
2144 :     if ($object->[0][$i][$j] != $argument->[0][$i][$j])
2145 :     {
2146 :     $result = 1;
2147 :     last NOTEQUAL;
2148 :     }
2149 :     }
2150 :     }
2151 :     return($result);
2152 :     }
2153 :     else
2154 :     {
2155 :     croak "MatrixReal1 $name: wrong argument type";
2156 :     }
2157 :     }
2158 :    
2159 :     sub _less_than
2160 :     {
2161 :     my($object,$argument,$flag) = @_;
2162 :     my($name) = "'<'"; #&_trace($name,$object,$argument,$flag);
2163 :    
2164 :     if ((defined $argument) && ref($argument) &&
2165 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2166 :     {
2167 :     if ((defined $flag) && $flag)
2168 :     {
2169 :     return( $argument->norm_one() < $object->norm_one() );
2170 :     }
2171 :     else
2172 :     {
2173 :     return( $object->norm_one() < $argument->norm_one() );
2174 :     }
2175 :     }
2176 :     elsif ((defined $argument) && !(ref($argument)))
2177 :     {
2178 :     if ((defined $flag) && $flag)
2179 :     {
2180 :     return( abs($argument) < $object->norm_one() );
2181 :     }
2182 :     else
2183 :     {
2184 :     return( $object->norm_one() < abs($argument) );
2185 :     }
2186 :     }
2187 :     else
2188 :     {
2189 :     croak "MatrixReal1 $name: wrong argument type";
2190 :     }
2191 :     }
2192 :    
2193 :     sub _less_than_or_equal
2194 :     {
2195 :     my($object,$argument,$flag) = @_;
2196 :     my($name) = "'<='"; #&_trace($name,$object,$argument,$flag);
2197 :    
2198 :     if ((defined $argument) && ref($argument) &&
2199 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2200 :     {
2201 :     if ((defined $flag) && $flag)
2202 :     {
2203 :     return( $argument->norm_one() <= $object->norm_one() );
2204 :     }
2205 :     else
2206 :     {
2207 :     return( $object->norm_one() <= $argument->norm_one() );
2208 :     }
2209 :     }
2210 :     elsif ((defined $argument) && !(ref($argument)))
2211 :     {
2212 :     if ((defined $flag) && $flag)
2213 :     {
2214 :     return( abs($argument) <= $object->norm_one() );
2215 :     }
2216 :     else
2217 :     {
2218 :     return( $object->norm_one() <= abs($argument) );
2219 :     }
2220 :     }
2221 :     else
2222 :     {
2223 :     croak "MatrixReal1 $name: wrong argument type";
2224 :     }
2225 :     }
2226 :    
2227 :     sub _greater_than
2228 :     {
2229 :     my($object,$argument,$flag) = @_;
2230 :     my($name) = "'>'"; #&_trace($name,$object,$argument,$flag);
2231 :    
2232 :     if ((defined $argument) && ref($argument) &&
2233 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2234 :     {
2235 :     if ((defined $flag) && $flag)
2236 :     {
2237 :     return( $argument->norm_one() > $object->norm_one() );
2238 :     }
2239 :     else
2240 :     {
2241 :     return( $object->norm_one() > $argument->norm_one() );
2242 :     }
2243 :     }
2244 :     elsif ((defined $argument) && !(ref($argument)))
2245 :     {
2246 :     if ((defined $flag) && $flag)
2247 :     {
2248 :     return( abs($argument) > $object->norm_one() );
2249 :     }
2250 :     else
2251 :     {
2252 :     return( $object->norm_one() > abs($argument) );
2253 :     }
2254 :     }
2255 :     else
2256 :     {
2257 :     croak "MatrixReal1 $name: wrong argument type";
2258 :     }
2259 :     }
2260 :    
2261 :     sub _greater_than_or_equal
2262 :     {
2263 :     my($object,$argument,$flag) = @_;
2264 :     my($name) = "'>='"; #&_trace($name,$object,$argument,$flag);
2265 :    
2266 :     if ((defined $argument) && ref($argument) &&
2267 :     (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
2268 :     {
2269 :     if ((defined $flag) && $flag)
2270 :     {
2271 :     return( $argument->norm_one() >= $object->norm_one() );
2272 :     }
2273 :     else
2274 :     {
2275 :     return( $object->norm_one() >= $argument->norm_one() );
2276 :     }
2277 :     }
2278 :     elsif ((defined $argument) && !(ref($argument)))
2279 :     {
2280 :     if ((defined $flag) && $flag)
2281 :     {
2282 :     return( abs($argument) >= $object->norm_one() );
2283 :     }
2284 :     else
2285 :     {
2286 :     return( $object->norm_one() >= abs($argument) );
2287 :     }
2288 :     }
2289 :     else
2290 :     {
2291 :     croak "MatrixReal1 $name: wrong argument type";
2292 :     }
2293 :     }
2294 :    
2295 :     sub _clone
2296 :     {
2297 :     my($object,$argument,$flag) = @_;
2298 :     # my($name) = "'='"; #&_trace($name,$object,$argument,$flag);
2299 :     my($temp);
2300 :    
2301 :     $temp = $object->new($object->[1],$object->[2]);
2302 :     $temp->copy($object);
2303 :     $temp->_undo_LR();
2304 :     return($temp);
2305 :     }
2306 :    
2307 :     sub _trace
2308 :     {
2309 :     my($text,$object,$argument,$flag) = @_;
2310 :    
2311 :     unless (defined $object) { $object = 'undef'; };
2312 :     unless (defined $argument) { $argument = 'undef'; };
2313 :     unless (defined $flag) { $flag = 'undef'; };
2314 :     if (ref($object)) { $object = ref($object); }
2315 :     if (ref($argument)) { $argument = ref($argument); }
2316 :     print "$text: \$obj='$object' \$arg='$argument' \$flag='$flag'\n";
2317 :     }
2318 :    
2319 :     1;
2320 :    
2321 :     __END__
2322 :    
2323 :     =head1 NAME
2324 :    
2325 :     MatrixReal1 - Matrix of Reals
2326 :    
2327 :     Implements the data type "matrix of reals" (and consequently also
2328 :     "vector of reals")
2329 :    
2330 :     =head1 DESCRIPTION
2331 :    
2332 :     Implements the data type "matrix of reals", which can be used almost
2333 :     like any other basic Perl type thanks to B<OPERATOR OVERLOADING>, i.e.,
2334 :    
2335 :     $product = $matrix1 * $matrix2;
2336 :    
2337 :     does what you would like it to do (a matrix multiplication).
2338 :    
2339 :     Also features many important operations and methods: matrix norm,
2340 :     matrix transposition, matrix inverse, determinant of a matrix, order
2341 :     and numerical condition of a matrix, scalar product of vectors, vector
2342 :     product of vectors, vector length, projection of row and column vectors,
2343 :     a comfortable way for reading in a matrix from a file, the keyboard or
2344 :     your code, and many more.
2345 :    
2346 :     Allows to solve linear equation systems using an efficient algorithm
2347 :     known as "L-R-decomposition" and several approximative (iterative) methods.
2348 :    
2349 :     Features an implementation of Kleene's algorithm to compute the minimal
2350 :     costs for all paths in a graph with weighted edges (the "weights" being
2351 :     the costs associated with each edge).
2352 :    
2353 :     =head1 SYNOPSIS
2354 :    
2355 :     =over 2
2356 :    
2357 :     =item *
2358 :    
2359 :     C<use MatrixReal1;>
2360 :    
2361 :     Makes the methods and overloaded operators of this module available
2362 :     to your program.
2363 :    
2364 :     =item *
2365 :    
2366 :     C<use MatrixReal1 qw(min max);>
2367 :    
2368 :     =item *
2369 :    
2370 :     C<use MatrixReal1 qw(:all);>
2371 :    
2372 :     Use one of these two variants to import (all) the functions which the module
2373 :     offers for export; currently these are "min()" and "max()".
2374 :    
2375 :     =item *
2376 :    
2377 :     C<$new_matrix = new MatrixReal1($rows,$columns);>
2378 :    
2379 :     The matrix object constructor method.
2380 :    
2381 :     Note that this method is implicitly called by many of the other methods
2382 :     in this module!
2383 :    
2384 :     =item *
2385 :    
2386 :     C<$new_matrix = MatrixReal1-E<gt>>C<new($rows,$columns);>
2387 :    
2388 :     An alternate way of calling the matrix object constructor method.
2389 :    
2390 :     =item *
2391 :    
2392 :     C<$new_matrix = $some_matrix-E<gt>>C<new($rows,$columns);>
2393 :    
2394 :     Still another way of calling the matrix object constructor method.
2395 :    
2396 :     Matrix "C<$some_matrix>" is not changed by this in any way.
2397 :    
2398 :     =item *
2399 :    
2400 :     C<$new_matrix = MatrixReal1-E<gt>>C<new_from_string($string);>
2401 :    
2402 :     This method allows you to read in a matrix from a string (for
2403 :     instance, from the keyboard, from a file or from your code).
2404 :    
2405 :     The syntax is simple: each row must start with "C<[ >" and end with
2406 :     "C< ]\n>" ("C<\n>" being the newline character and "C< >" a space or
2407 :     tab) and contain one or more numbers, all separated from each other
2408 :     by spaces or tabs.
2409 :    
2410 :     Additional spaces or tabs can be added at will, but no comments.
2411 :    
2412 :     Examples:
2413 :    
2414 :     $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
2415 :     $matrix = MatrixReal1->new_from_string($string);
2416 :     print "$matrix";
2417 :    
2418 :     By the way, this prints
2419 :    
2420 :     [ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
2421 :     [ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
2422 :     [ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
2423 :    
2424 :     But you can also do this in a much more comfortable way using the
2425 :     shell-like "here-document" syntax:
2426 :    
2427 :     $matrix = MatrixReal1->new_from_string(<<'MATRIX');
2428 :     [ 1 0 0 0 0 0 1 ]
2429 :     [ 0 1 0 0 0 0 0 ]
2430 :     [ 0 0 1 0 0 0 0 ]
2431 :     [ 0 0 0 1 0 0 0 ]
2432 :     [ 0 0 0 0 1 0 0 ]
2433 :     [ 0 0 0 0 0 1 0 ]
2434 :     [ 1 0 0 0 0 0 -1 ]
2435 :     MATRIX
2436 :    
2437 :     You can even use variables in the matrix:
2438 :    
2439 :     $c1 = 2 / 3;
2440 :     $c2 = -2 / 5;
2441 :     $c3 = 26 / 9;
2442 :    
2443 :     $matrix = MatrixReal1->new_from_string(<<"MATRIX");
2444 :    
2445 :     [ 3 2 0 ]
2446 :     [ 0 3 2 ]
2447 :     [ $c1 $c2 $c3 ]
2448 :    
2449 :     MATRIX
2450 :    
2451 :     (Remember that you may use spaces and tabs to format the matrix to
2452 :     your taste)
2453 :    
2454 :     Note that this method uses exactly the same representation for a
2455 :     matrix as the "stringify" operator "": this means that you can convert
2456 :     any matrix into a string with C<$string = "$matrix";> and read it back
2457 :     in later (for instance from a file!).
2458 :    
2459 :     Note however that you may suffer a precision loss in this process
2460 :     because only 13 digits are supported in the mantissa when printed!!
2461 :    
2462 :     If the string you supply (or someone else supplies) does not obey
2463 :     the syntax mentioned above, an exception is raised, which can be
2464 :     caught by "eval" as follows:
2465 :    
2466 :     print "Please enter your matrix (in one line): ";
2467 :     $string = <STDIN>;
2468 :     $string =~ s/\\n/\n/g;
2469 :     eval { $matrix = MatrixReal1->new_from_string($string); };
2470 :     if ($@)
2471 :     {
2472 :     print "$@";
2473 :     # ...
2474 :     # (error handling)
2475 :     }
2476 :     else
2477 :     {
2478 :     # continue...
2479 :     }
2480 :    
2481 :     or as follows:
2482 :    
2483 :     eval { $matrix = MatrixReal1->new_from_string(<<"MATRIX"); };
2484 :     [ 3 2 0 ]
2485 :     [ 0 3 2 ]
2486 :     [ $c1 $c2 $c3 ]
2487 :     MATRIX
2488 :     if ($@)
2489 :     # ...
2490 :    
2491 :     Actually, the method shown above for reading a matrix from the keyboard
2492 :     is a little awkward, since you have to enter a lot of "\n"'s for the
2493 :     newlines.
2494 :    
2495 :     A better way is shown in this piece of code:
2496 :    
2497 :     while (1)
2498 :     {
2499 :     print "\nPlease enter your matrix ";
2500 :     print "(multiple lines, <ctrl-D> = done):\n";
2501 :     eval { $new_matrix =
2502 :     MatrixReal1->new_from_string(join('',<STDIN>)); };
2503 :     if ($@)
2504 :     {
2505 :     $@ =~ s/\s+at\b.*?$//;
2506 :     print "${@}Please try again.\n";
2507 :     }
2508 :     else { last; }
2509 :     }
2510 :    
2511 :     Possible error messages of the "new_from_string()" method are:
2512 :    
2513 :     MatrixReal1::new_from_string(): syntax error in input string
2514 :     MatrixReal1::new_from_string(): empty input string
2515 :    
2516 :     If the input string has rows with varying numbers of columns,
2517 :     the following warning will be printed to STDERR:
2518 :    
2519 :     MatrixReal1::new_from_string(): missing elements will be set to zero!
2520 :    
2521 :     If everything is okay, the method returns an object reference to the
2522 :     (newly allocated) matrix containing the elements you specified.
2523 :    
2524 :     =item *
2525 :    
2526 :     C<$new_matrix = $some_matrix-E<gt>shadow();>
2527 :    
2528 :     Returns an object reference to a B<NEW> but B<EMPTY> matrix
2529 :     (filled with zero's) of the B<SAME SIZE> as matrix "C<$some_matrix>".
2530 :    
2531 :     Matrix "C<$some_matrix>" is not changed by this in any way.
2532 :    
2533 :     =item *
2534 :    
2535 :     C<$matrix1-E<gt>copy($matrix2);>
2536 :    
2537 :     Copies the contents of matrix "C<$matrix2>" to an B<ALREADY EXISTING>
2538 :     matrix "C<$matrix1>" (which must have the same size as matrix "C<$matrix2>"!).
2539 :    
2540 :     Matrix "C<$matrix2>" is not changed by this in any way.
2541 :    
2542 :     =item *
2543 :    
2544 :     C<$twin_matrix = $some_matrix-E<gt>clone();>
2545 :    
2546 :     Returns an object reference to a B<NEW> matrix of the B<SAME SIZE> as
2547 :     matrix "C<$some_matrix>". The contents of matrix "C<$some_matrix>" have
2548 :     B<ALREADY BEEN COPIED> to the new matrix "C<$twin_matrix>".
2549 :    
2550 :     Matrix "C<$some_matrix>" is not changed by this in any way.
2551 :    
2552 :     =item *
2553 :    
2554 :     C<$row_vector = $matrix-E<gt>row($row);>
2555 :    
2556 :     This is a projection method which returns an object reference to
2557 :     a B<NEW> matrix (which in fact is a (row) vector since it has only
2558 :     one row) to which row number "C<$row>" of matrix "C<$matrix>" has
2559 :     already been copied.
2560 :    
2561 :     Matrix "C<$matrix>" is not changed by this in any way.
2562 :    
2563 :     =item *
2564 :    
2565 :     C<$column_vector = $matrix-E<gt>column($column);>
2566 :    
2567 :     This is a projection method which returns an object reference to
2568 :     a B<NEW> matrix (which in fact is a (column) vector since it has
2569 :     only one column) to which column number "C<$column>" of matrix
2570 :     "C<$matrix>" has already been copied.
2571 :    
2572 :     Matrix "C<$matrix>" is not changed by this in any way.
2573 :    
2574 :     =item *
2575 :    
2576 :     C<$matrix-E<gt>zero();>
2577 :    
2578 :     Assigns a zero to every element of the matrix "C<$matrix>", i.e.,
2579 :     erases all values previously stored there, thereby effectively
2580 :     transforming the matrix into a "zero"-matrix or "null"-matrix,
2581 :     the neutral element of the addition operation in a Ring.
2582 :    
2583 :     (For instance the (quadratic) matrices with "n" rows and columns
2584 :     and matrix addition and multiplication form a Ring. Most prominent
2585 :     characteristic of a Ring is that multiplication is not commutative,
2586 :     i.e., in general, "C<matrix1 * matrix2>" is not the same as
2587 :     "C<matrix2 * matrix1>"!)
2588 :    
2589 :     =item *
2590 :    
2591 :     C<$matrix-E<gt>one();>
2592 :    
2593 :     Assigns one's to the elements on the main diagonal (elements (1,1),
2594 :     (2,2), (3,3) and so on) of matrix "C<$matrix>" and zero's to all others,
2595 :     thereby erasing all values previously stored there and transforming the
2596 :     matrix into a "one"-matrix, the neutral element of the multiplication
2597 :     operation in a Ring.
2598 :    
2599 :     (If the matrix is quadratic (which this method doesn't require, though),
2600 :     then multiplying this matrix with itself yields this same matrix again,
2601 :     and multiplying it with some other matrix leaves that other matrix
2602 :     unchanged!)
2603 :    
2604 :     =item *
2605 :    
2606 :     C<$matrix-E<gt>assign($row,$column,$value);>
2607 :    
2608 :     Explicitly assigns a value "C<$value>" to a single element of the
2609 :     matrix "C<$matrix>", located in row "C<$row>" and column "C<$column>",
2610 :     thereby replacing the value previously stored there.
2611 :    
2612 :     =item *
2613 :    
2614 :     C<$value = $matrix-E<gt>>C<element($row,$column);>
2615 :    
2616 :     Returns the value of a specific element of the matrix "C<$matrix>",
2617 :     located in row "C<$row>" and column "C<$column>".
2618 :    
2619 :     =item *
2620 :    
2621 :     C<($rows,$columns) = $matrix-E<gt>dim();>
2622 :    
2623 :     Returns a list of two items, representing the number of rows
2624 :     and columns the given matrix "C<$matrix>" contains.
2625 :    
2626 :     =item *
2627 :    
2628 :     C<$norm_one = $matrix-E<gt>norm_one();>
2629 :    
2630 :     Returns the "one"-norm of the given matrix "C<$matrix>".
2631 :    
2632 :     The "one"-norm is defined as follows:
2633 :    
2634 :     For each column, the sum of the absolute values of the elements in the
2635 :     different rows of that column is calculated. Finally, the maximum
2636 :     of these sums is returned.
2637 :    
2638 :     Note that the "one"-norm and the "maximum"-norm are mathematically
2639 :     equivalent, although for the same matrix they usually yield a different
2640 :     value.
2641 :    
2642 :     Therefore, you should only compare values that have been calculated
2643 :     using the same norm!
2644 :    
2645 :     Throughout this package, the "one"-norm is (arbitrarily) used
2646 :     for all comparisons, for the sake of uniformity and comparability,
2647 :     except for the iterative methods "solve_GSM()", "solve_SSM()" and
2648 :     "solve_RM()" which use either norm depending on the matrix itself.
2649 :    
2650 :     =item *
2651 :    
2652 :     C<$norm_max = $matrix-E<gt>norm_max();>
2653 :    
2654 :     Returns the "maximum"-norm of the given matrix "C<$matrix>".
2655 :    
2656 :     The "maximum"-norm is defined as follows:
2657 :    
2658 :     For each row, the sum of the absolute values of the elements in the
2659 :     different columns of that row is calculated. Finally, the maximum
2660 :     of these sums is returned.
2661 :    
2662 :     Note that the "maximum"-norm and the "one"-norm are mathematically
2663 :     equivalent, although for the same matrix they usually yield a different
2664 :     value.
2665 :    
2666 :     Therefore, you should only compare values that have been calculated
2667 :     using the same norm!
2668 :    
2669 :     Throughout this package, the "one"-norm is (arbitrarily) used
2670 :     for all comparisons, for the sake of uniformity and comparability,
2671 :     except for the iterative methods "solve_GSM()", "solve_SSM()" and
2672 :     "solve_RM()" which use either norm depending on the matrix itself.
2673 :    
2674 :     =item *
2675 :    
2676 :     C<$matrix1-E<gt>negate($matrix2);>
2677 :    
2678 :     Calculates the negative of matrix "C<$matrix2>" (i.e., multiplies
2679 :     all elements with "-1") and stores the result in matrix "C<$matrix1>"
2680 :     (which must already exist and have the same size as matrix "C<$matrix2>"!).
2681 :    
2682 :     This operation can also be carried out "in-place", i.e., input and
2683 :     output matrix may be identical.
2684 :    
2685 :     =item *
2686 :    
2687 :     C<$matrix1-E<gt>transpose($matrix2);>
2688 :    
2689 :     Calculates the transposed matrix of matrix "C<$matrix2>" and stores
2690 :     the result in matrix "C<$matrix1>" (which must already exist and have
2691 :     the same size as matrix "C<$matrix2>"!).
2692 :    
2693 :     This operation can also be carried out "in-place", i.e., input and
2694 :     output matrix may be identical.
2695 :    
2696 :     Transposition is a symmetry operation: imagine you rotate the matrix
2697 :     along the axis of its main diagonal (going through elements (1,1),
2698 :     (2,2), (3,3) and so on) by 180 degrees.
2699 :    
2700 :     Another way of looking at it is to say that rows and columns are
2701 :     swapped. In fact the contents of element C<(i,j)> are swapped
2702 :     with those of element C<(j,i)>.
2703 :    
2704 :     Note that (especially for vectors) it makes a big difference if you
2705 :     have a row vector, like this:
2706 :    
2707 :     [ -1 0 1 ]
2708 :    
2709 :     or a column vector, like this:
2710 :    
2711 :     [ -1 ]
2712 :     [ 0 ]
2713 :     [ 1 ]
2714 :    
2715 :     the one vector being the transposed of the other!
2716 :    
2717 :     This is especially true for the matrix product of two vectors:
2718 :    
2719 :     [ -1 ]
2720 :     [ -1 0 1 ] * [ 0 ] = [ 2 ] , whereas
2721 :     [ 1 ]
2722 :    
2723 :     * [ -1 0 1 ]
2724 :     [ -1 ] [ 1 0 -1 ]
2725 :     [ 0 ] * [ -1 0 1 ] = [ -1 ] [ 1 0 -1 ] = [ 0 0 0 ]
2726 :     [ 1 ] [ 0 ] [ 0 0 0 ] [ -1 0 1 ]
2727 :     [ 1 ] [ -1 0 1 ]
2728 :    
2729 :     So be careful about what you really mean!
2730 :    
2731 :     Hint: throughout this module, whenever a vector is explicitly required
2732 :     for input, a B<COLUMN> vector is expected!
2733 :    
2734 :     =item *
2735 :    
2736 :     C<$matrix1-E<gt>add($matrix2,$matrix3);>
2737 :    
2738 :     Calculates the sum of matrix "C<$matrix2>" and matrix "C<$matrix3>"
2739 :     and stores the result in matrix "C<$matrix1>" (which must already exist
2740 :     and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
2741 :    
2742 :     This operation can also be carried out "in-place", i.e., the output and
2743 :     one (or both) of the input matrices may be identical.
2744 :    
2745 :     =item *
2746 :    
2747 :     C<$matrix1-E<gt>subtract($matrix2,$matrix3);>
2748 :    
2749 :     Calculates the difference of matrix "C<$matrix2>" minus matrix "C<$matrix3>"
2750 :     and stores the result in matrix "C<$matrix1>" (which must already exist
2751 :     and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
2752 :    
2753 :     This operation can also be carried out "in-place", i.e., the output and
2754 :     one (or both) of the input matrices may be identical.
2755 :    
2756 :     Note that this operation is the same as
2757 :     C<$matrix1-E<gt>add($matrix2,-$matrix3);>, although the latter is
2758 :     a little less efficient.
2759 :    
2760 :     =item *
2761 :    
2762 :     C<$matrix1-E<gt>multiply_scalar($matrix2,$scalar);>
2763 :    
2764 :     Calculates the product of matrix "C<$matrix2>" and the number "C<$scalar>"
2765 :     (i.e., multiplies each element of matrix "C<$matrix2>" with the factor
2766 :     "C<$scalar>") and stores the result in matrix "C<$matrix1>" (which must
2767 :     already exist and have the same size as matrix "C<$matrix2>"!).
2768 :    
2769 :     This operation can also be carried out "in-place", i.e., input and
2770 :     output matrix may be identical.
2771 :    
2772 :     =item *
2773 :    
2774 :     C<$product_matrix = $matrix1-E<gt>multiply($matrix2);>
2775 :    
2776 :     Calculates the product of matrix "C<$matrix1>" and matrix "C<$matrix2>"
2777 :     and returns an object reference to a new matrix "C<$product_matrix>" in
2778 :     which the result of this operation has been stored.
2779 :    
2780 :     Note that the dimensions of the two matrices "C<$matrix1>" and "C<$matrix2>"
2781 :     (i.e., their numbers of rows and columns) must harmonize in the following
2782 :     way (example):
2783 :    
2784 :     [ 2 2 ]
2785 :     [ 2 2 ]
2786 :     [ 2 2 ]
2787 :    
2788 :     [ 1 1 1 ] [ * * ]
2789 :     [ 1 1 1 ] [ * * ]
2790 :     [ 1 1 1 ] [ * * ]
2791 :     [ 1 1 1 ] [ * * ]
2792 :    
2793 :     I.e., the number of columns of matrix "C<$matrix1>" has to be the same
2794 :     as the number of rows of matrix "C<$matrix2>".
2795 :    
2796 :     The number of rows and columns of the resulting matrix "C<$product_matrix>"
2797 :     is determined by the number of rows of matrix "C<$matrix1>" and the number
2798 :     of columns of matrix "C<$matrix2>", respectively.
2799 :    
2800 :     =item *
2801 :    
2802 :     C<$minimum = MatrixReal1::min($number1,$number2);>
2803 :    
2804 :     Returns the minimum of the two numbers "C<number1>" and "C<number2>".
2805 :    
2806 :     =item *
2807 :    
2808 :     C<$minimum = MatrixReal1::max($number1,$number2);>
2809 :    
2810 :     Returns the maximum of the two numbers "C<number1>" and "C<number2>".
2811 :    
2812 :     =item *
2813 :    
2814 :     C<$minimal_cost_matrix = $cost_matrix-E<gt>kleene();>
2815 :    
2816 :     Copies the matrix "C<$cost_matrix>" (which has to be quadratic!) to
2817 :     a new matrix of the same size (i.e., "clones" the input matrix) and
2818 :     applies Kleene's algorithm to it.
2819 :    
2820 :     See L<Math::Kleene(3)> for more details about this algorithm!
2821 :    
2822 :     The method returns an object reference to the new matrix.
2823 :    
2824 :     Matrix "C<$cost_matrix>" is not changed by this method in any way.
2825 :    
2826 :     =item *
2827 :    
2828 :     C<($norm_matrix,$norm_vector) = $matrix-E<gt>normalize($vector);>
2829 :    
2830 :     This method is used to improve the numerical stability when solving
2831 :     linear equation systems.
2832 :    
2833 :     Suppose you have a matrix "A" and a vector "b" and you want to find
2834 :     out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
2835 :     solves the equation system represented by the matrix "A" and the
2836 :     vector "b".
2837 :    
2838 :     Applying this method to the pair (A,b) yields a pair (A',b') where
2839 :     each row has been divided by (the absolute value of) the greatest
2840 :     coefficient appearing in that row. So this coefficient becomes equal
2841 :     to "1" (or "-1") in the new pair (A',b') (all others become smaller
2842 :     than one and greater than minus one).
2843 :    
2844 :     Note that this operation does not change the equation system itself
2845 :     because the same division is carried out on either side of the equation
2846 :     sign!
2847 :    
2848 :     The method requires a quadratic (!) matrix "C<$matrix>" and a vector
2849 :     "C<$vector>" for input (the vector must be a column vector with the same
2850 :     number of rows as the input matrix) and returns a list of two items
2851 :     which are object references to a new matrix and a new vector, in this
2852 :     order.
2853 :    
2854 :     The output matrix and vector are clones of the input matrix and vector
2855 :     to which the operation explained above has been applied.
2856 :    
2857 :     The input matrix and vector are not changed by this in any way.
2858 :    
2859 :     Example of how this method can affect the result of the methods to solve
2860 :     equation systems (explained immediately below following this method):
2861 :    
2862 :     Consider the following little program:
2863 :    
2864 :     #!perl -w
2865 :    
2866 :     use MatrixReal1 qw(new_from_string);
2867 :    
2868 :     $A = MatrixReal1->new_from_string(<<"MATRIX");
2869 :     [ 1 2 3 ]
2870 :     [ 5 7 11 ]
2871 :     [ 23 19 13 ]
2872 :     MATRIX
2873 :    
2874 :     $b = MatrixReal1->new_from_string(<<"MATRIX");
2875 :     [ 0 ]
2876 :     [ 1 ]
2877 :     [ 29 ]
2878 :     MATRIX
2879 :    
2880 :     $LR = $A->decompose_LR();
2881 :     if (($dim,$x,$B) = $LR->solve_LR($b))
2882 :     {
2883 :     $test = $A * $x;
2884 :     print "x = \n$x";
2885 :     print "A * x = \n$test";
2886 :     }
2887 :    
2888 :     ($A_,$b_) = $A->normalize($b);
2889 :    
2890 :     $LR = $A_->decompose_LR();
2891 :     if (($dim,$x,$B) = $LR->solve_LR($b_))
2892 :     {
2893 :     $test = $A * $x;
2894 :     print "x = \n$x";
2895 :     print "A * x = \n$test";
2896 :     }
2897 :    
2898 :     This will print:
2899 :    
2900 :     x =
2901 :     [ 1.000000000000E+00 ]
2902 :     [ 1.000000000000E+00 ]
2903 :     [ -1.000000000000E+00 ]
2904 :     A * x =
2905 :     [ 4.440892098501E-16 ]
2906 :     [ 1.000000000000E+00 ]
2907 :     [ 2.900000000000E+01 ]
2908 :     x =
2909 :     [ 1.000000000000E+00 ]
2910 :     [ 1.000000000000E+00 ]
2911 :     [ -1.000000000000E+00 ]
2912 :     A * x =
2913 :     [ 0.000000000000E+00 ]
2914 :     [ 1.000000000000E+00 ]
2915 :     [ 2.900000000000E+01 ]
2916 :    
2917 :     You can see that in the second example (where "normalize()" has been used),
2918 :     the result is "better", i.e., more accurate!
2919 :    
2920 :     =item *
2921 :    
2922 :     C<$LR_matrix = $matrix-E<gt>decompose_LR();>
2923 :    
2924 :     This method is needed to solve linear equation systems.
2925 :    
2926 :     Suppose you have a matrix "A" and a vector "b" and you want to find
2927 :     out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
2928 :     solves the equation system represented by the matrix "A" and the
2929 :     vector "b".
2930 :    
2931 :     You might also have a matrix "A" and a whole bunch of different
2932 :     vectors "b1".."bk" for which you need to find vectors "x1".."xk"
2933 :     so that C<A * xi = bi>, for C<i=1..k>.
2934 :    
2935 :     Using Gaussian transformations (multiplying a row or column with
2936 :     a factor, swapping two rows or two columns and adding a multiple
2937 :     of one row or column to another), it is possible to decompose any
2938 :     matrix "A" into two triangular matrices, called "L" and "R" (for
2939 :     "Left" and "Right").
2940 :    
2941 :     "L" has one's on the main diagonal (the elements (1,1), (2,2), (3,3)
2942 :     and so so), non-zero values to the left and below of the main diagonal
2943 :     and all zero's in the upper right half of the matrix.
2944 :    
2945 :     "R" has non-zero values on the main diagonal as well as to the right
2946 :     and above of the main diagonal and all zero's in the lower left half
2947 :     of the matrix, as follows:
2948 :    
2949 :     [ 1 0 0 0 0 ] [ x x x x x ]
2950 :     [ x 1 0 0 0 ] [ 0 x x x x ]
2951 :     L = [ x x 1 0 0 ] R = [ 0 0 x x x ]
2952 :     [ x x x 1 0 ] [ 0 0 0 x x ]
2953 :     [ x x x x 1 ] [ 0 0 0 0 x ]
2954 :    
2955 :     Note that "C<L * R>" is equivalent to matrix "A" in the sense that
2956 :     C<L * R * x = b E<lt>==E<gt> A * x = b> for all vectors "x", leaving
2957 :     out of account permutations of the rows and columns (these are taken
2958 :     care of "magically" by this module!) and numerical errors.
2959 :    
2960 :     Trick:
2961 :    
2962 :     Because we know that "L" has one's on its main diagonal, we can
2963 :     store both matrices together in the same array without information
2964 :     loss! I.e.,
2965 :    
2966 :     [ R R R R R ]
2967 :     [ L R R R R ]
2968 :     LR = [ L L R R R ]
2969 :     [ L L L R R ]
2970 :     [ L L L L R ]
2971 :    
2972 :     Beware, though, that "LR" and "C<L * R>" are not the same!!!
2973 :    
2974 :     Note also that for the same reason, you cannot apply the method "normalize()"
2975 :     to an "LR" decomposition matrix. Trying to do so will yield meaningless
2976 :     rubbish!
2977 :    
2978 :     (You need to apply "normalize()" to each pair (Ai,bi) B<BEFORE> decomposing
2979 :     the matrix "Ai'"!)
2980 :    
2981 :     Now what does all this help us in solving linear equation systems?
2982 :    
2983 :     It helps us because a triangular matrix is the next best thing
2984 :     that can happen to us besides a diagonal matrix (a matrix that
2985 :     has non-zero values only on its main diagonal - in which case
2986 :     the solution is trivial, simply divide "C<b[i]>" by "C<A[i,i]>"
2987 :     to get "C<x[i]>"!).
2988 :    
2989 :     To find the solution to our problem "C<A * x = b>", we divide this
2990 :     problem in parts: instead of solving C<A * x = b> directly, we first
2991 :     decompose "A" into "L" and "R" and then solve "C<L * y = b>" and
2992 :     finally "C<R * x = y>" (motto: divide and rule!).
2993 :    
2994 :     From the illustration above it is clear that solving "C<L * y = b>"
2995 :     and "C<R * x = y>" is straightforward: we immediately know that
2996 :     C<y[1] = b[1]>. We then deduce swiftly that
2997 :    
2998 :     y[2] = b[2] - L[2,1] * y[1]
2999 :    
3000 :     (and we know "C<y[1]>" by now!), that
3001 :    
3002 :     y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
3003 :    
3004 :     and so on.
3005 :    
3006 :     Having effortlessly calculated the vector "y", we now proceed to
3007 :     calculate the vector "x" in a similar fashion: we see immediately
3008 :     that C<x[n] = y[n] / R[n,n]>. It follows that
3009 :    
3010 :     x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
3011 :    
3012 :     and
3013 :    
3014 :     x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
3015 :     / R[n-2,n-2]
3016 :    
3017 :     and so on.
3018 :    
3019 :     You can see that - especially when you have many vectors "b1".."bk"
3020 :     for which you are searching solutions to C<A * xi = bi> - this scheme
3021 :     is much more efficient than a straightforward, "brute force" approach.
3022 :    
3023 :     This method requires a quadratic matrix as its input matrix.
3024 :    
3025 :     If you don't have that many equations, fill up with zero's (i.e., do
3026 :     nothing to fill the superfluous rows if it's a "fresh" matrix, i.e.,
3027 :     a matrix that has been created with "new()" or "shadow()").
3028 :    
3029 :     The method returns an object reference to a new matrix containing the
3030 :     matrices "L" and "R".
3031 :    
3032 :     The input matrix is not changed by this method in any way.
3033 :    
3034 :     Note that you can "copy()" or "clone()" the result of this method without
3035 :     losing its "magical" properties (for instance concerning the hidden
3036 :     permutations of its rows and columns).
3037 :    
3038 :     However, as soon as you are applying any method that alters the contents
3039 :     of the matrix, its "magical" properties are stripped off, and the matrix
3040 :     immediately reverts to an "ordinary" matrix (with the values it just happens
3041 :     to contain at that moment, be they meaningful as an ordinary matrix or not!).
3042 :    
3043 :     =item *
3044 :    
3045 :     C<($dimension,$x_vector,$base_matrix) = $LR_matrix>C<-E<gt>>C<solve_LR($b_vector);>
3046 :    
3047 :     Use this method to actually solve an equation system.
3048 :    
3049 :     Matrix "C<$LR_matrix>" must be a (quadratic) matrix returned by the
3050 :     method "decompose_LR()", the LR decomposition matrix of the matrix
3051 :     "A" of your equation system C<A * x = b>.
3052 :    
3053 :     The input vector "C<$b_vector>" is the vector "b" in your equation system
3054 :     C<A * x = b>, which must be a column vector and have the same number of
3055 :     rows as the input matrix "C<$LR_matrix>".
3056 :    
3057 :     The method returns a list of three items if a solution exists or an
3058 :     empty list otherwise (!).
3059 :    
3060 :     Therefore, you should always use this method like this:
3061 :    
3062 :     if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
3063 :     {
3064 :     # do something with the solution...
3065 :     }
3066 :     else
3067 :     {
3068 :     # do something with the fact that there is no solution...
3069 :     }
3070 :    
3071 :     The three items returned are: the dimension "C<$dimension>" of the solution
3072 :     space (which is zero if only one solution exists, one if the solution is
3073 :     a straight line, two if the solution is a plane, and so on), the solution
3074 :     vector "C<$x_vector>" (which is the vector "x" of your equation system
3075 :     C<A * x = b>) and a matrix "C<$base_matrix>" representing a base of the
3076 :     solution space (a set of vectors which put up the solution space like
3077 :     the spokes of an umbrella).
3078 :    
3079 :     Only the first "C<$dimension>" columns of this base matrix actually
3080 :     contain entries, the remaining columns are all zero.
3081 :    
3082 :     Now what is all this stuff with that "base" good for?
3083 :    
3084 :     The output vector "x" is B<ALWAYS> a solution of your equation system
3085 :     C<A * x = b>.
3086 :    
3087 :     But also any vector "C<$vector>"
3088 :    
3089 :     $vector = $x_vector->clone();
3090 :    
3091 :     $machine_infinity = 1E+99; # or something like that
3092 :    
3093 :     for ( $i = 1; $i <= $dimension; $i++ )
3094 :     {
3095 :     $vector += rand($machine_infinity) * $base_matrix->column($i);
3096 :     }
3097 :    
3098 :     is a solution to your problem C<A * x = b>, i.e., if "C<$A_matrix>" contains
3099 :     your matrix "A", then
3100 :    
3101 :     print abs( $A_matrix * $vector - $b_vector ), "\n";
3102 :    
3103 :     should print a number around 1E-16 or so!
3104 :    
3105 :     By the way, note that you can actually calculate those vectors "C<$vector>"
3106 :     a little more efficient as follows:
3107 :    
3108 :     $rand_vector = $x_vector->shadow();
3109 :    
3110 :     $machine_infinity = 1E+99; # or something like that
3111 :    
3112 :     for ( $i = 1; $i <= $dimension; $i++ )
3113 :     {
3114 :     $rand_vector->assign($i,1, rand($machine_infinity) );
3115 :     }
3116 :    
3117 :     $vector = $x_vector + ( $base_matrix * $rand_vector );
3118 :    
3119 :     Note that the input matrix and vector are not changed by this method
3120 :     in any way.
3121 :    
3122 :     =item *
3123 :    
3124 :     C<$inverse_matrix = $LR_matrix-E<gt>invert_LR();>
3125 :    
3126 :     Use this method to calculate the inverse of a given matrix "C<$LR_matrix>",
3127 :     which must be a (quadratic) matrix returned by the method "decompose_LR()".
3128 :    
3129 :     The method returns an object reference to a new matrix of the same size as
3130 :     the input matrix containing the inverse of the matrix that you initially
3131 :     fed into "decompose_LR()" B<IF THE INVERSE EXISTS>, or an empty list
3132 :     otherwise.
3133 :    
3134 :     Therefore, you should always use this method in the following way:
3135 :    
3136 :     if ( $inverse_matrix = $LR->invert_LR() )
3137 :     {
3138 :     # do something with the inverse matrix...
3139 :     }
3140 :     else
3141 :     {
3142 :     # do something with the fact that there is no inverse matrix...
3143 :     }
3144 :    
3145 :     Note that by definition (disregarding numerical errors), the product
3146 :     of the initial matrix and its inverse (or vice-versa) is always a matrix
3147 :     containing one's on the main diagonal (elements (1,1), (2,2), (3,3) and
3148 :     so on) and zero's elsewhere.
3149 :    
3150 :     The input matrix is not changed by this method in any way.
3151 :    
3152 :     =item *
3153 :    
3154 :     C<$condition = $matrix-E<gt>condition($inverse_matrix);>
3155 :    
3156 :     In fact this method is just a shortcut for
3157 :    
3158 :     abs($matrix) * abs($inverse_matrix)
3159 :    
3160 :     Both input matrices must be quadratic and have the same size, and the result
3161 :     is meaningful only if one of them is the inverse of the other (for instance,
3162 :     as returned by the method "invert_LR()").
3163 :    
3164 :     The number returned is a measure of the "condition" of the given matrix
3165 :     "C<$matrix>", i.e., a measure of the numerical stability of the matrix.
3166 :    
3167 :     This number is always positive, and the smaller its value, the better the
3168 :     condition of the matrix (the better the stability of all subsequent
3169 :     computations carried out using this matrix).
3170 :    
3171 :     Numerical stability means for example that if
3172 :    
3173 :     abs( $vec_correct - $vec_with_error ) < $epsilon
3174 :    
3175 :     holds, there must be a "C<$delta>" which doesn't depend on the vector
3176 :     "C<$vec_correct>" (nor "C<$vec_with_error>", by the way) so that
3177 :    
3178 :     abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
3179 :    
3180 :     also holds.
3181 :    
3182 :     =item *
3183 :    
3184 :     C<$determinant = $LR_matrix-E<gt>det_LR();>
3185 :    
3186 :     Calculates the determinant of a matrix, whose LR decomposition matrix
3187 :     "C<$LR_matrix>" must be given (which must be a (quadratic) matrix
3188 :     returned by the method "decompose_LR()").
3189 :    
3190 :     In fact the determinant is a by-product of the LR decomposition: It is
3191 :     (in principle, that is, except for the sign) simply the product of the
3192 :     elements on the main diagonal (elements (1,1), (2,2), (3,3) and so on)
3193 :     of the LR decomposition matrix.
3194 :    
3195 :     (The sign is taken care of "magically" by this module)
3196 :    
3197 :     =item *
3198 :    
3199 :     C<$order = $LR_matrix-E<gt>order_LR();>
3200 :    
3201 :     Calculates the order (called "Rang" in German) of a matrix, whose
3202 :     LR decomposition matrix "C<$LR_matrix>" must be given (which must
3203 :     be a (quadratic) matrix returned by the method "decompose_LR()").
3204 :    
3205 :     This number is a measure of the number of linear independent row
3206 :     and column vectors (= number of linear independent equations in
3207 :     the case of a matrix representing an equation system) of the
3208 :     matrix that was initially fed into "decompose_LR()".
3209 :    
3210 :     If "n" is the number of rows and columns of the (quadratic!) matrix,
3211 :     then "n - order" is the dimension of the solution space of the
3212 :     associated equation system.
3213 :    
3214 :     =item *
3215 :    
3216 :     C<$scalar_product = $vector1-E<gt>scalar_product($vector2);>
3217 :    
3218 :     Returns the scalar product of vector "C<$vector1>" and vector "C<$vector2>".
3219 :    
3220 :     Both vectors must be column vectors (i.e., a matrix having
3221 :     several rows but only one column).
3222 :    
3223 :     This is a (more efficient!) shortcut for
3224 :    
3225 :     $temp = ~$vector1 * $vector2;
3226 :     $scalar_product = $temp->element(1,1);
3227 :    
3228 :     or the sum C<i=1..n> of the products C<vector1[i] * vector2[i]>.
3229 :    
3230 :     Provided none of the two input vectors is the null vector, then
3231 :     the two vectors are orthogonal, i.e., have an angle of 90 degrees
3232 :     between them, exactly when their scalar product is zero, and
3233 :     vice-versa.
3234 :    
3235 :     =item *
3236 :    
3237 :     C<$vector_product = $vector1-E<gt>vector_product($vector2);>
3238 :    
3239 :     Returns the vector product of vector "C<$vector1>" and vector "C<$vector2>".
3240 :    
3241 :     Both vectors must be column vectors (i.e., a matrix having several rows
3242 :     but only one column).
3243 :    
3244 :     Currently, the vector product is only defined for 3 dimensions (i.e.,
3245 :     vectors with 3 rows); all other vectors trigger an error message.
3246 :    
3247 :     In 3 dimensions, the vector product of two vectors "x" and "y"
3248 :     is defined as
3249 :    
3250 :     | x[1] y[1] e[1] |
3251 :     determinant | x[2] y[2] e[2] |
3252 :     | x[3] y[3] e[3] |
3253 :    
3254 :     where the "C<x[i]>" and "C<y[i]>" are the components of the two vectors
3255 :     "x" and "y", respectively, and the "C<e[i]>" are unity vectors (i.e.,
3256 :     vectors with a length equal to one) with a one in row "i" and zero's
3257 :     elsewhere (this means that you have numbers and vectors as elements
3258 :     in this matrix!).
3259 :    
3260 :     This determinant evaluates to the rather simple formula
3261 :    
3262 :     z[1] = x[2] * y[3] - x[3] * y[2]
3263 :     z[2] = x[3] * y[1] - x[1] * y[3]
3264 :     z[3] = x[1] * y[2] - x[2] * y[1]
3265 :    
3266 :     A characteristic property of the vector product is that the resulting
3267 :     vector is orthogonal to both of the input vectors (if neither of both
3268 :     is the null vector, otherwise this is trivial), i.e., the scalar product
3269 :     of each of the input vectors with the resulting vector is always zero.
3270 :    
3271 :     =item *
3272 :    
3273 :     C<$length = $vector-E<gt>length();>
3274 :    
3275 :     This is actually a shortcut for
3276 :    
3277 :     $length = sqrt( $vector->scalar_product($vector) );
3278 :    
3279 :     and returns the length of a given (column!) vector "C<$vector>".
3280 :    
3281 :     Note that the "length" calculated by this method is in fact the
3282 :     "two"-norm of a vector "C<$vector>"!
3283 :    
3284 :     The general definition for norms of vectors is the following:
3285 :    
3286 :     sub vector_norm
3287 :     {
3288 :     croak "Usage: \$norm = \$vector->vector_norm(\$n);"
3289 :     if (@_ != 2);
3290 :    
3291 :     my($vector,$n) = @_;
3292 :     my($rows,$cols) = ($vector->[1],$vector->[2]);
3293 :     my($k,$comp,$sum);
3294 :    
3295 :     croak "MatrixReal1::vector_norm(): vector is not a column vector"
3296 :     unless ($cols == 1);
3297 :    
3298 :     croak "MatrixReal1::vector_norm(): norm index must be > 0"
3299 :     unless ($n > 0);
3300 :    
3301 :     croak "MatrixReal1::vector_norm(): norm index must be integer"
3302 :     unless ($n == int($n));
3303 :    
3304 :     $sum = 0;
3305 :     for ( $k = 0; $k < $rows; $k++ )
3306 :     {
3307 :     $comp = abs( $vector->[0][$k][0] );
3308 :     $sum += $comp ** $n;
3309 :     }
3310 :     return( $sum ** (1 / $n) );
3311 :     }
3312 :    
3313 :     Note that the case "n = 1" is the "one"-norm for matrices applied to a
3314 :     vector, the case "n = 2" is the euclidian norm or length of a vector,
3315 :     and if "n" goes to infinity, you have the "infinity"- or "maximum"-norm
3316 :     for matrices applied to a vector!
3317 :    
3318 :     =item *
3319 :    
3320 :     C<$xn_vector = $matrix-E<gt>>C<solve_GSM($x0_vector,$b_vector,$epsilon);>
3321 :    
3322 :     =item *
3323 :    
3324 :     C<$xn_vector = $matrix-E<gt>>C<solve_SSM($x0_vector,$b_vector,$epsilon);>
3325 :    
3326 :     =item *
3327 :    
3328 :     C<$xn_vector = $matrix-E<gt>>C<solve_RM($x0_vector,$b_vector,$weight,$epsilon);>
3329 :    
3330 :     In some cases it might not be practical or desirable to solve an
3331 :     equation system "C<A * x = b>" using an analytical algorithm like
3332 :     the "decompose_LR()" and "solve_LR()" method pair.
3333 :    
3334 :     In fact in some cases, due to the numerical properties (the "condition")
3335 :     of the matrix "A", the numerical error of the obtained result can be
3336 :     greater than by using an approximative (iterative) algorithm like one
3337 :     of the three implemented here.
3338 :    
3339 :     All three methods, GSM ("Global Step Method" or "Gesamtschrittverfahren"),
3340 :     SSM ("Single Step Method" or "Einzelschrittverfahren") and RM ("Relaxation
3341 :     Method" or "Relaxationsverfahren"), are fix-point iterations, that is, can
3342 :     be described by an iteration function "C<x(t+1) = Phi( x(t) )>" which has
3343 :     the property:
3344 :    
3345 :     Phi(x) = x <==> A * x = b
3346 :    
3347 :     We can define "C<Phi(x)>" as follows:
3348 :    
3349 :     Phi(x) := ( En - A ) * x + b
3350 :    
3351 :     where "En" is a matrix of the same size as "A" ("n" rows and columns)
3352 :     with one's on its main diagonal and zero's elsewhere.
3353 :    
3354 :     This function has the required property.
3355 :    
3356 :     Proof:
3357 :    
3358 :     A * x = b
3359 :    
3360 :     <==> -( A * x ) = -b
3361 :    
3362 :     <==> -( A * x ) + x = -b + x
3363 :    
3364 :     <==> -( A * x ) + x + b = x
3365 :    
3366 :     <==> x - ( A * x ) + b = x
3367 :    
3368 :     <==> ( En - A ) * x + b = x
3369 :    
3370 :     This last step is true because
3371 :    
3372 :     x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
3373 :    
3374 :     is the same as
3375 :    
3376 :     ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
3377 :    
3378 :     qed
3379 :    
3380 :     Note that actually solving the equation system "C<A * x = b>" means
3381 :     to calculate
3382 :    
3383 :     a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] = b[i]
3384 :    
3385 :     <==> a[i,i] x[i] =
3386 :     b[i]
3387 :     - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
3388 :     + a[i,i] x[i]
3389 :    
3390 :     <==> x[i] =
3391 :     ( b[i]
3392 :     - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
3393 :     + a[i,i] x[i]
3394 :     ) / a[i,i]
3395 :    
3396 :     <==> x[i] =
3397 :     ( b[i] -
3398 :     ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
3399 :     a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
3400 :     ) / a[i,i]
3401 :    
3402 :     There is one major restriction, though: a fix-point iteration is
3403 :     guaranteed to converge only if the first derivative of the iteration
3404 :     function has an absolute value less than one in an area around the
3405 :     point "C<x(*)>" for which "C<Phi( x(*) ) = x(*)>" is to be true, and
3406 :     if the start vector "C<x(0)>" lies within that area!
3407 :    
3408 :     This is best verified grafically, which unfortunately is impossible
3409 :     to do in this textual documentation!
3410 :    
3411 :     See literature on Numerical Analysis for details!
3412 :    
3413 :     In our case, this restriction translates to the following three conditions:
3414 :    
3415 :     There must exist a norm so that the norm of the matrix of the iteration
3416 :     function, C<( En - A )>, has a value less than one, the matrix "A" may
3417 :     not have any zero value on its main diagonal and the initial vector
3418 :     "C<x(0)>" must be "good enough", i.e., "close enough" to the solution
3419 :     "C<x(*)>".
3420 :    
3421 :     (Remember school math: the first derivative of a straight line given by
3422 :     "C<y = a * x + b>" is "a"!)
3423 :    
3424 :     The three methods expect a (quadratic!) matrix "C<$matrix>" as their
3425 :     first argument, a start vector "C<$x0_vector>", a vector "C<$b_vector>"
3426 :     (which is the vector "b" in your equation system "C<A * x = b>"), in the
3427 :     case of the "Relaxation Method" ("RM"), a real number "C<$weight>" best
3428 :     between zero and two, and finally an error limit (real number) "C<$epsilon>".
3429 :    
3430 :     (Note that the weight "C<$weight>" used by the "Relaxation Method" ("RM")
3431 :     is B<NOT> checked to lie within any reasonable range!)
3432 :    
3433 :     The three methods first test the first two conditions of the three
3434 :     conditions listed above and return an empty list if these conditions
3435 :     are not fulfilled.
3436 :    
3437 :     Therefore, you should always test their return value using some
3438 :     code like:
3439 :    
3440 :     if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
3441 :     {
3442 :     # do something with the solution...
3443 :     }
3444 :     else
3445 :     {
3446 :     # do something with the fact that there is no solution...
3447 :     }
3448 :    
3449 :     Otherwise, they iterate until C<abs( Phi(x) - x ) E<lt> epsilon>.
3450 :    
3451 :     (Beware that theoretically, infinite loops might result if the starting
3452 :     vector is too far "off" the solution! In practice, this shouldn't be
3453 :     a problem. Anyway, you can always press <ctrl-C> if you think that the
3454 :     iteration takes too long!)
3455 :    
3456 :     The difference between the three methods is the following:
3457 :    
3458 :     In the "Global Step Method" ("GSM"), the new vector "C<x(t+1)>"
3459 :     (called "y" here) is calculated from the vector "C<x(t)>"
3460 :     (called "x" here) according to the formula:
3461 :    
3462 :     y[i] =
3463 :     ( b[i]
3464 :     - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
3465 :     a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
3466 :     ) / a[i,i]
3467 :    
3468 :     In the "Single Step Method" ("SSM"), the components of the vector
3469 :     "C<x(t+1)>" which have already been calculated are used to calculate
3470 :     the remaining components, i.e.
3471 :    
3472 :     y[i] =
3473 :     ( b[i]
3474 :     - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
3475 :     a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
3476 :     ) / a[i,i]
3477 :    
3478 :     In the "Relaxation method" ("RM"), the components of the vector
3479 :     "C<x(t+1)>" are calculated by "mixing" old and new value (like
3480 :     cold and hot water), and the weight "C<$weight>" determines the
3481 :     "aperture" of both the "hot water tap" as well as of the "cold
3482 :     water tap", according to the formula:
3483 :    
3484 :     y[i] =
3485 :     ( b[i]
3486 :     - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
3487 :     a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
3488 :     ) / a[i,i]
3489 :     y[i] = weight * y[i] + (1 - weight) * x[i]
3490 :    
3491 :     Note that the weight "C<$weight>" should be greater than zero and
3492 :     less than two (!).
3493 :    
3494 :     The three methods are supposed to be of different efficiency.
3495 :     Experiment!
3496 :    
3497 :     Remember that in most cases, it is probably advantageous to first
3498 :     "normalize()" your equation system prior to solving it!
3499 :    
3500 :     =back
3501 :    
3502 :     =head2 Eigensystems
3503 :    
3504 :     =over 2
3505 :    
3506 :     =item *
3507 :    
3508 :     C<$matrix-E<gt>is_symmetric();>
3509 :    
3510 :     Returns a boolean value indicating if the given matrix is
3511 :     symmetric (B<M>[I<i>,I<j>]=B<M>[I<j>,I<i>]). This is equivalent to
3512 :     C<($matrix == ~$matrix)> but without memory allocation.
3513 :    
3514 :     =item *
3515 :    
3516 :     C<($l, $V) = $matrix-E<gt>sym_diagonalize();>
3517 :    
3518 :     This method performs the diagonalization of the quadratic
3519 :     I<symmetric> matrix B<M> stored in $matrix.
3520 :     On output, B<l> is a column vector containing all the eigenvalues
3521 :     of B<M> and B<V> is an orthogonal matrix which columns are the
3522 :     corresponding normalized eigenvectors.
3523 :     The primary property of an eigenvalue I<l> and an eigenvector
3524 :     B<x> is of course that: B<M> * B<x> = I<l> * B<x>.
3525 :    
3526 :     The method uses a Householder reduction to tridiagonal form
3527 :     followed by a QL algoritm with implicit shifts on this
3528 :     tridiagonal. (The tridiagonal matrix is kept internally
3529 :     in a compact form in this routine to save memory.)
3530 :     In fact, this routine wraps the householder() and
3531 :     tri_diagonalize() methods described below when their
3532 :     intermediate results are not desired.
3533 :     The overall algorithmic complexity of this technique
3534 :     is O(N^3). According to several books, the coefficient
3535 :     hidden by the 'O' is one of the best possible for general
3536 :     (symmetric) matrixes.
3537 :    
3538 :     =item *
3539 :    
3540 :     C<($T, $Q) = $matrix-E<gt>householder();>
3541 :    
3542 :     This method performs the Householder algorithm which reduces
3543 :     the I<n> by I<n> real I<symmetric> matrix B<M> contained
3544 :     in $matrix to tridiagonal form.
3545 :     On output, B<T> is a symmetric tridiagonal matrix (only
3546 :     diagonal and off-diagonal elements are non-zero) and B<Q>
3547 :     is an I<orthogonal> matrix performing the tranformation
3548 :     between B<M> and B<T> (C<$M == $Q * $T * ~$Q>).
3549 :    
3550 :     =item *
3551 :    
3552 :     C<($l, $V) = $T-E<gt>tri_diagonalize([$Q]);>
3553 :    
3554 :     This method diagonalizes the symmetric tridiagonal
3555 :     matrix B<T>. On output, $l and $V are similar to the
3556 :     output values described for sym_diagonalize().
3557 :    
3558 :     The optional argument $Q corresponds to an orthogonal
3559 :     transformation matrix B<Q> that should be used additionally
3560 :     during B<V> (eigenvectors) computation. It should be supplied
3561 :     if the desired eigenvectors correspond to a more general
3562 :     symmetric matrix B<M> previously reduced by the
3563 :     householder() method, not a mere tridiagonal. If B<T> is
3564 :     really a tridiagonal matrix, B<Q> can be omitted (it
3565 :     will be internally created in fact as an identity matrix).
3566 :     The method uses a QL algorithm (with implicit shifts).
3567 :    
3568 :     =item *
3569 :    
3570 :     C<$l = $matrix-E<gt>sym_eigenvalues();>
3571 :    
3572 :     This method computes the eigenvalues of the quadratic
3573 :     I<symmetric> matrix B<M> stored in $matrix.
3574 :     On output, B<l> is a column vector containing all the eigenvalues
3575 :     of B<M>. Eigenvectors are not computed (on the contrary of
3576 :     C<sym_diagonalize()>) and this method is more efficient
3577 :     (even though it uses a similar algorithm with two phases).
3578 :     However, understand that the algorithmic complexity of this
3579 :     technique is still also O(N^3). But the coefficient hidden
3580 :     by the 'O' is better by a factor of..., well, see your
3581 :     benchmark, it's wiser.
3582 :    
3583 :     This routine wraps the householder_tridiagonal() and
3584 :     tri_eigenvalues() methods described below when the
3585 :     intermediate tridiagonal matrix is not needed.
3586 :    
3587 :     =item *
3588 :    
3589 :     C<$T = $matrix-E<gt>householder_tridiagonal();>
3590 :    
3591 :     This method performs the Householder algorithm which reduces
3592 :     the I<n> by I<n> real I<symmetric> matrix B<M> contained
3593 :     in $matrix to tridiagonal form.
3594 :     On output, B<T> is the obtained symmetric tridiagonal matrix
3595 :     (only diagonal and off-diagonal elements are non-zero). The
3596 :     operation is similar to the householder() method, but potentially
3597 :     a little more efficient as the transformation matrix is not
3598 :     computed.
3599 :    
3600 :     =item *
3601 :    
3602 :     C<$l = $T-E<gt>tri_eigenvalues();>
3603 :    
3604 :     This method compute the eigenvalues of the symmetric
3605 :     tridiagonal matrix B<T>. On output, $l is a vector
3606 :     containing the eigenvalues (similar to C<sym_eigenvalues()>).
3607 :     This method is much more efficient than tri_diagonalize()
3608 :     when eigenvectors are not needed.
3609 :    
3610 :     =back
3611 :    
3612 :     =head1 OVERLOADED OPERATORS
3613 :    
3614 :     =head2 SYNOPSIS
3615 :    
3616 :     =over 2
3617 :    
3618 :     =item *
3619 :    
3620 :     Unary operators:
3621 :    
3622 :     "C<->", "C<~>", "C<abs>", C<test>, "C<!>", 'C<"">'
3623 :    
3624 :     =item *
3625 :    
3626 :     Binary (arithmetic) operators:
3627 :    
3628 :     "C<+>", "C<->", "C<*>"
3629 :    
3630 :     =item *
3631 :    
3632 :     Binary (relational) operators:
3633 :    
3634 :     "C<==>", "C<!=>", "C<E<lt>>", "C<E<lt>=>", "C<E<gt>>", "C<E<gt>=>"
3635 :    
3636 :     "C<eq>", "C<ne>", "C<lt>", "C<le>", "C<gt>", "C<ge>"
3637 :    
3638 :     Note that the latter ("C<eq>", "C<ne>", ... ) are just synonyms
3639 :     of the former ("C<==>", "C<!=>", ... ), defined for convenience
3640 :     only.
3641 :    
3642 :     =back
3643 :    
3644 :     =head2 DESCRIPTION
3645 :    
3646 :     =over 5
3647 :    
3648 :     =item '-'
3649 :    
3650 :     Unary minus
3651 :    
3652 :     Returns the negative of the given matrix, i.e., the matrix with
3653 :     all elements multiplied with the factor "-1".
3654 :    
3655 :     Example:
3656 :    
3657 :     $matrix = -$matrix;
3658 :    
3659 :     =item '~'
3660 :    
3661 :     Transposition
3662 :    
3663 :     Returns the transposed of the given matrix.
3664 :    
3665 :     Examples:
3666 :    
3667 :     $temp = ~$vector * $vector;
3668 :     $length = sqrt( $temp->element(1,1) );
3669 :    
3670 :     if (~$matrix == $matrix) { # matrix is symmetric ... }
3671 :    
3672 :     =item abs
3673 :    
3674 :     Norm
3675 :    
3676 :     Returns the "one"-Norm of the given matrix.
3677 :    
3678 :     Example:
3679 :    
3680 :     $error = abs( $A * $x - $b );
3681 :    
3682 :     =item test
3683 :    
3684 :     Boolean test
3685 :    
3686 :     Tests wether there is at least one non-zero element in the matrix.
3687 :    
3688 :     Example:
3689 :    
3690 :     if ($xn_vector) { # result of iteration is not zero ... }
3691 :    
3692 :     =item '!'
3693 :    
3694 :     Negated boolean test
3695 :    
3696 :     Tests wether the matrix contains only zero's.
3697 :    
3698 :     Examples:
3699 :    
3700 :     if (! $b_vector) { # heterogenous equation system ... }
3701 :     else { # homogenous equation system ... }
3702 :    
3703 :     unless ($x_vector) { # not the null-vector! }
3704 :    
3705 :     =item '""""'
3706 :    
3707 :     "Stringify" operator
3708 :    
3709 :     Converts the given matrix into a string.
3710 :    
3711 :     Uses scientific representation to keep precision loss to a minimum in case
3712 :     you want to read this string back in again later with "new_from_string()".
3713 :    
3714 :     Uses a 13-digit mantissa and a 20-character field for each element so that
3715 :     lines will wrap nicely on an 80-column screen.
3716 :    
3717 :     Examples:
3718 :    
3719 :     $matrix = MatrixReal1->new_from_string(<<"MATRIX");
3720 :     [ 1 0 ]
3721 :     [ 0 -1 ]
3722 :     MATRIX
3723 :     print "$matrix";
3724 :    
3725 :     [ 1.000000000000E+00 0.000000000000E+00 ]
3726 :     [ 0.000000000000E+00 -1.000000000000E+00 ]
3727 :    
3728 :     $string = "$matrix";
3729 :     $test = MatrixReal1->new_from_string($string);
3730 :     if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
3731 :    
3732 :     =item '+'
3733 :    
3734 :     Addition
3735 :    
3736 :     Returns the sum of the two given matrices.
3737 :    
3738 :     Examples:
3739 :    
3740 :     $matrix_S = $matrix_A + $matrix_B;
3741 :    
3742 :     $matrix_A += $matrix_B;
3743 :    
3744 :     =item '-'
3745 :    
3746 :     Subtraction
3747 :    
3748 :     Returns the difference of the two given matrices.
3749 :    
3750 :     Examples:
3751 :    
3752 :     $matrix_D = $matrix_A - $matrix_B;
3753 :    
3754 :     $matrix_A -= $matrix_B;
3755 :    
3756 :     Note that this is the same as:
3757 :    
3758 :     $matrix_S = $matrix_A + -$matrix_B;
3759 :    
3760 :     $matrix_A += -$matrix_B;
3761 :    
3762 :     (The latter are less efficient, though)
3763 :    
3764 :     =item '*'
3765 :    
3766 :     Multiplication
3767 :    
3768 :     Returns the matrix product of the two given matrices or
3769 :     the product of the given matrix and scalar factor.
3770 :    
3771 :     Examples:
3772 :    
3773 :     $matrix_P = $matrix_A * $matrix_B;
3774 :    
3775 :     $matrix_A *= $matrix_B;
3776 :    
3777 :     $vector_b = $matrix_A * $vector_x;
3778 :    
3779 :     $matrix_B = -1 * $matrix_A;
3780 :    
3781 :     $matrix_B = $matrix_A * -1;
3782 :    
3783 :     $matrix_A *= -1;
3784 :    
3785 :     =item '=='
3786 :    
3787 :     Equality
3788 :    
3789 :     Tests two matrices for equality.
3790 :    
3791 :     Example:
3792 :    
3793 :     if ( $A * $x == $b ) { print "EUREKA!\n"; }
3794 :    
3795 :     Note that in most cases, due to numerical errors (due to the finite
3796 :     precision of computer arithmetics), it is a bad idea to compare two
3797 :     matrices or vectors this way.
3798 :    
3799 :     Better use the norm of the difference of the two matrices you want
3800 :     to compare and compare that norm with a small number, like this:
3801 :    
3802 :     if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
3803 :    
3804 :     =item '!='
3805 :    
3806 :     Inequality
3807 :    
3808 :     Tests two matrices for inequality.
3809 :    
3810 :     Example:
3811 :    
3812 :     while ($x0_vector != $xn_vector) { # proceed with iteration ... }
3813 :    
3814 :     (Stops when the iteration becomes stationary)
3815 :    
3816 :     Note that (just like with the '==' operator), it is usually a bad idea
3817 :     to compare matrices or vectors this way. Compare the norm of the difference
3818 :     of the two matrices with a small number instead.
3819 :    
3820 :     =item 'E<lt>'
3821 :    
3822 :     Less than
3823 :    
3824 :     Examples:
3825 :    
3826 :     if ( $matrix1 < $matrix2 ) { # ... }
3827 :    
3828 :     if ( $vector < $epsilon ) { # ... }
3829 :    
3830 :     if ( 1E-12 < $vector ) { # ... }
3831 :    
3832 :     if ( $A * $x - $b < 1E-12 ) { # ... }
3833 :    
3834 :     These are just shortcuts for saying:
3835 :    
3836 :     if ( abs($matrix1) < abs($matrix2) ) { # ... }
3837 :    
3838 :     if ( abs($vector) < abs($epsilon) ) { # ... }
3839 :    
3840 :     if ( abs(1E-12) < abs($vector) ) { # ... }
3841 :    
3842 :     if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
3843 :    
3844 :     Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
3845 :    
3846 :     =item 'E<lt>='
3847 :    
3848 :     Less than or equal
3849 :    
3850 :     As with the '<' operator, this is just a shortcut for the same expression
3851 :     with "abs()" around all arguments.
3852 :    
3853 :     Example:
3854 :    
3855 :     if ( $A * $x - $b <= 1E-12 ) { # ... }
3856 :    
3857 :     which in fact is the same as:
3858 :    
3859 :     if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
3860 :    
3861 :     Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
3862 :    
3863 :     =item 'E<gt>'
3864 :    
3865 :     Greater than
3866 :    
3867 :     As with the '<' and '<=' operator, this
3868 :    
3869 :     if ( $xn - $x0 > 1E-12 ) { # ... }
3870 :    
3871 :     is just a shortcut for:
3872 :    
3873 :     if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
3874 :    
3875 :     Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
3876 :    
3877 :     =item 'E<gt>='
3878 :    
3879 :     Greater than or equal
3880 :    
3881 :     As with the '<', '<=' and '>' operator, the following
3882 :    
3883 :     if ( $LR >= $A ) { # ... }
3884 :    
3885 :     is simply a shortcut for:
3886 :    
3887 :     if ( abs($LR) >= abs($A) ) { # ... }
3888 :    
3889 :     Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
3890 :    
3891 :     =back
3892 :    
3893 :     =head1 SEE ALSO
3894 :    
3895 :     Math::MatrixBool(3), DFA::Kleene(3), Math::Kleene(3),
3896 :     Set::IntegerRange(3), Set::IntegerFast(3).
3897 :    
3898 :     =head1 VERSION
3899 :    
3900 :     This man page documents MatrixReal1 version 1.3.
3901 :    
3902 :     =head1 AUTHORS
3903 :    
3904 :     Steffen Beyer <sb@sdm.de>, Rodolphe Ortalo <ortalo@laas.fr>.
3905 :    
3906 :     =head1 CREDITS
3907 :    
3908 :     Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
3909 :     Algebra and Linear Algebra at the university (RWTH Aachen, Germany), and
3910 :     to Prof. Esser and his assistant, Mr. Jarausch, for their fascinating
3911 :     lectures in Numerical Analysis!
3912 :    
3913 :     =head1 COPYRIGHT
3914 :    
3915 :     Copyright (c) 1996, 1997, 1999 by Steffen Beyer and Rodolphe Ortalo.
3916 :     All rights reserved.
3917 :    
3918 :     =head1 LICENSE AGREEMENT
3919 :    
3920 :     This package is free software; you can redistribute it and/or
3921 :     modify it under the same terms as Perl itself.
3922 :    

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9