Annotation of /trunk/webwork/system/courseScripts/MatrixReal1.pm

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

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

aubreyja at gmail dot com	ViewVC Help
Powered by ViewVC 1.0.9