Annotation of /branches/rel-2-0-b1-opt/pg/lib/MatrixReal1.pm

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

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

aubreyja at gmail dot com	ViewVC Help
Powered by ViewVC 1.0.9