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