Parent Directory
| 
Revision Log
Copied PG modules over from old system. -sam
1 # 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(a*a + b*b), 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 |