1Math::MatrixReal(3) User Contributed Perl Documentation Math::MatrixReal(3)
2
6 Math::MatrixReal - Matrix of Reals
7 Implements the data type "matrix of real numbers" (and consequently
9 also "vector of real numbers").
10
12 my $a = Math::MatrixReal->new_random(5, 5);
13 my $b = $a->new_random(10, 30, { symmetric=>1, bounded_by=>[-1,1] });
15 my $c = $b * $a ** 3;
17 my $d = $b->new_from_rows( [ [ 5, 3 ,4], [3, 4, 5], [ 2, 4, 1 ] ] );
19 print $a;
21 my $row = ($a * $b)->row(3);
23 my $col = (5*$c)->col(2);
25 my $transpose = ~$c;
27 my $transpose = $c->transpose;
29 my $inverse = $a->inverse;
31 my $inverse = 1/$a;
33 my $inverse = $a ** -1;
35 my $determinant= $a->det;
37 • $matrix->display_precision($integer)
39 Sets the default precision when matrices are printed or
41 stringified. $matrix->display_precision(0) will only show the
42 integer part of all the entries of $matrix and
43 $matrix->display_precision() will return to the default scientific
44 display notation. This method does not effect the precision of the
45 calculations.
46
48 Constructor Methods And Such
49 • use Math::MatrixReal;
50 Makes the methods and overloaded operators of this module available
52 to your program.
53 • $new_matrix = new Math::MatrixReal($rows,$columns);
55 The matrix object constructor method. A new matrix of size $rows by
57 $columns will be created, with the value 0.0 for all elements.
58 Note that this method is implicitly called by many of the other
60 methods in this module.
61 • $new_matrix = $some_matrix->new($rows,$columns);
63 Another way of calling the matrix object constructor method.
65 Matrix $some_matrix is not changed by this in any way.
67 • $new_matrix = $matrix->new_from_cols( [
69 $column_vector|$array_ref|$string, ... ] )
70 Creates a new matrix given a reference to an array of any of the
72 following:
73 • column vectors ( n by 1 Math::MatrixReal matrices )
75 • references to arrays
77 • strings properly formatted to create a column with
79 Math::MatrixReal's new_from_string command
80 You may mix and match these as you wish. However, all must be of
82 the same dimension--no padding happens automatically. Example:
83 my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
85 print $matrix;
86 will print
88 [ 1.000000000000E+00 3.000000000000E+00 ]
90 [ 2.000000000000E+00 4.000000000000E+00 ]
91 • new_from_rows( [ $row_vector|$array_ref|$string, ... ] )
93 Creates a new matrix given a reference to an array of any of the
95 following:
96 • row vectors ( 1 by n Math::MatrixReal matrices )
98 • references to arrays
100 • strings properly formatted to create a row with
102 Math::MatrixReal's new_from_string command
103 You may mix and match these as you wish. However, all must be of
105 the same dimension--no padding happens automatically. Example:
106 my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
108 print $matrix;
109 will print
111 [ 1.000000000000E+00 2.000000000000E+00 ]
113 [ 3.000000000000E+00 4.000000000000E+00 ]
114 • $new_matrix = Math::MatrixReal->new_random($rows, $cols, %options
116 );
117 This method allows you to create a random matrix with various
119 properties controlled by the %options matrix, which is optional.
120 The default values of the %options matrix are { integer => 0,
121 symmetric => 0, tridiagonal => 0, diagonal => 0, bounded_by =>
122 [0,10] } .
123 Example:
125 $matrix = Math::MatrixReal->new_random(4, { diagonal => 1, integer => 1 } );
127 print $matrix;
128 will print a 4x4 random diagonal matrix with integer entries
130 between zero and ten, something like
131 [ 5.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
133 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
134 [ 0.000000000000E+00 0.000000000000E+00 1.000000000000E+00 0.000000000000E+00 ]
135 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 8.000000000000E+00 ]
136 • $new_matrix = Math::MatrixReal->new_diag( $array_ref );
138 This method allows you to create a diagonal matrix by only
140 specifying the diagonal elements. Example:
141 $matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
143 print $matrix;
144 will print
146 [ 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
148 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
149 [ 0.000000000000E+00 0.000000000000E+00 3.000000000000E+00 0.000000000000E+00 ]
150 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 4.000000000000E+00 ]
151 • $new_matrix = Math::MatrixReal->new_tridiag( $lower, $diag, $upper
153 );
154 This method allows you to create a tridiagonal matrix by only
156 specifying the lower diagonal, diagonal and upper diagonal,
157 respectively.
158 $matrix = Math::MatrixReal->new_tridiag( [ 6, 4, 2 ], [1,2,3,4], [1, 8, 9] );
160 print $matrix;
161 will print
163 [ 1.000000000000E+00 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
165 [ 6.000000000000E+00 2.000000000000E+00 8.000000000000E+00 0.000000000000E+00 ]
166 [ 0.000000000000E+00 4.000000000000E+00 3.000000000000E+00 9.000000000000E+00 ]
167 [ 0.000000000000E+00 0.000000000000E+00 2.000000000000E+00 4.000000000000E+00 ]
168 • $new_matrix = Math::MatrixReal->new_from_string($string);
170 This method allows you to read in a matrix from a string (for
172 instance, from the keyboard, from a file or from your code).
173 The syntax is simple: each row must start with ""[ "" and end with
175 "" ]\n"" (""\n"" being the newline character and "" "" a space or
176 tab) and contain one or more numbers, all separated from each other
177 by spaces or tabs.
178 Additional spaces or tabs can be added at will, but no comments.
180 Examples:
182 $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
184 $matrix = Math::MatrixReal->new_from_string($string);
185 print "$matrix";
186 By the way, this prints
188 [ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
190 [ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
191 [ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
192 But you can also do this in a much more comfortable way using the
194 shell-like "here-document" syntax:
195 $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
197 [ 1 0 0 0 0 0 1 ]
198 [ 0 1 0 0 0 0 0 ]
199 [ 0 0 1 0 0 0 0 ]
200 [ 0 0 0 1 0 0 0 ]
201 [ 0 0 0 0 1 0 0 ]
202 [ 0 0 0 0 0 1 0 ]
203 [ 1 0 0 0 0 0 -1 ]
204 MATRIX
205 You can even use variables in the matrix:
207 $c1 = 2 / 3;
209 $c2 = -2 / 5;
210 $c3 = 26 / 9;
211 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
213 [ 3 2 0 ]
215 [ 0 3 2 ]
216 [ $c1 $c2 $c3 ]
217 MATRIX
219 (Remember that you may use spaces and tabs to format the matrix to
221 your taste)
222 Note that this method uses exactly the same representation for a
224 matrix as the "stringify" operator "": this means that you can
225 convert any matrix into a string with "$string = "$matrix";" and
226 read it back in later (for instance from a file!).
227 Note however that you may suffer a precision loss in this process
229 because only 13 digits are supported in the mantissa when printed!!
230 If the string you supply (or someone else supplies) does not obey
232 the syntax mentioned above, an exception is raised, which can be
233 caught by "eval" as follows:
234 print "Please enter your matrix (in one line): ";
236 $string = <STDIN>;
237 $string =~ s/\\n/\n/g;
238 eval { $matrix = Math::MatrixReal->new_from_string($string); };
239 if ($@)
240 {
241 print "$@";
242 # ...
243 # (error handling)
244 }
245 else
246 {
247 # continue...
248 }
249 or as follows:
251 eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
253 [ 3 2 0 ]
254 [ 0 3 2 ]
255 [ $c1 $c2 $c3 ]
256 MATRIX
257 if ($@)
258 # ...
259 Actually, the method shown above for reading a matrix from the
261 keyboard is a little awkward, since you have to enter a lot of
262 "\n"'s for the newlines.
263 A better way is shown in this piece of code:
265 while (1)
267 {
268 print "\nPlease enter your matrix ";
269 print "(multiple lines, <ctrl-D> = done):\n";
270 eval { $new_matrix =
271 Math::MatrixReal->new_from_string(join('',<STDIN>)); };
272 if ($@)
273 {
274 $@ =~ s/\s+at\b.*?$//;
275 print "${@}Please try again.\n";
276 }
277 else { last; }
278 }
279 Possible error messages of the "new_from_string()" method are:
281 Math::MatrixReal::new_from_string(): syntax error in input string
283 Math::MatrixReal::new_from_string(): empty input string
284 If the input string has rows with varying numbers of columns, the
286 following warning will be printed to STDERR:
287 Math::MatrixReal::new_from_string(): missing elements will be set to zero!
289 If everything is okay, the method returns an object reference to
291 the (newly allocated) matrix containing the elements you specified.
292 • $new_matrix = $some_matrix->shadow();
294 Returns an object reference to a NEW but EMPTY matrix (filled with
296 zero's) of the SAME SIZE as matrix "$some_matrix".
297 Matrix "$some_matrix" is not changed by this in any way.
299 • $matrix1->copy($matrix2);
301 Copies the contents of matrix "$matrix2" to an ALREADY EXISTING
303 matrix "$matrix1" (which must have the same size as matrix
304 "$matrix2"!).
305 Matrix "$matrix2" is not changed by this in any way.
307 • $twin_matrix = $some_matrix->clone();
309 Returns an object reference to a NEW matrix of the SAME SIZE as
311 matrix "$some_matrix". The contents of matrix "$some_matrix" have
312 ALREADY BEEN COPIED to the new matrix "$twin_matrix". This is the
313 method that the operator "=" is overloaded to when you type "$a =
314 $b", when $a and $b are matrices.
315 Matrix "$some_matrix" is not changed by this in any way.
317 • $matrix = Math::MatrixReal->reshape($rows, $cols, $array_ref);
319 Return a matrix with the specified dimensions ($rows x $cols)
321 whose elements are taken from the array reference $array_ref. The
322 elements of the matrix are accessed in column-major order (like
323 Fortran arrays are stored).
324 $matrix = Math::MatrixReal->reshape(4, 3, [1..12]);
326 Creates the following matrix:
328 [ 1 5 9 ]
330 [ 2 6 10 ]
331 [ 3 7 11 ]
332 [ 4 8 12 ]
333 Matrix Row, Column and Element operations
335 • $value = $matrix->element($row,$column);
336 Returns the value of a specific element of the matrix "$matrix",
338 located in row "$row" and column "$column".
339 NOTE: Unlike Perl, matrices are indexed with base-one indexes.
341 Thus, the first element of the matrix is placed in the first line,
342 first column:
343 $elem = $matrix->element(1, 1); # first element of the matrix.
345 • $matrix->assign($row,$column,$value);
347 Explicitly assigns a value "$value" to a single element of the
349 matrix "$matrix", located in row "$row" and column "$column",
350 thereby replacing the value previously stored there.
351 • $row_vector = $matrix->row($row);
353 This is a projection method which returns an object reference to a
355 NEW matrix (which in fact is a (row) vector since it has only one
356 row) to which row number "$row" of matrix "$matrix" has already
357 been copied.
358 Matrix "$matrix" is not changed by this in any way.
360 • $column_vector = $matrix->column($column);
362 This is a projection method which returns an object reference to a
364 NEW matrix (which in fact is a (column) vector since it has only
365 one column) to which column number "$column" of matrix "$matrix"
366 has already been copied.
367 Matrix "$matrix" is not changed by this in any way.
369 • @all_elements = $matrix->as_list;
371 Get the contents of a Math::MatrixReal object as a Perl list.
373 Example:
375 my $matrix = Math::MatrixReal->new_from_rows([ [1, 2], [3, 4] ]);
377 my @list = $matrix->as_list; # 1, 2, 3, 4
378 This method is suitable for use with OpenGL. For example, there is
380 need to rotate model around X-axis to 90 degrees clock-wise. That
381 could be achieved via:
382 use Math::Trig;
384 use OpenGL;
385 ...;
386 my $axis = [1, 0, 0];
387 my $angle = 90;
388 ...
389 my ($x, $y, $z) = @$axis;
390 my $f = $angle;
391 my $cos_f = cos(deg2rad($f));
392 my $sin_f = sin(deg2rad($f));
393 my $rotation = Math::MatrixReal->new_from_rows([
394 [$cos_f+(1-$cos_f)*$x**2, (1-$cos_f)*$x*$y-$sin_f*$z, (1-$cos_f)*$x*$z+$sin_f*$y, 0 ],
395 [(1-$cos_f)*$y*$z+$sin_f*$z, $cos_f+(1-$cos_f)*$y**2 , (1-$cos_f)*$y*$z-$sin_f*$x, 0 ],
396 [(1-$cos_f)*$z*$x-$sin_f*$y, (1-$cos_f)*$z*$y+$sin_f*$x, $cos_f+(1-$cos_f)*$z**2 ,0 ],
397 [0, 0, 0, 1 ],
398 ]);
399 ...;
400 my $model_initial = Math::MatrixReal->new_diag( [1, 1, 1, 1] ); # identity matrix
401 my $model = $model_initial * $rotation;
402 $model = ~$model; # OpenGL operates on transposed matrices
403 my $model_oga = OpenGL::Array->new_list(GL_FLOAT, $model->as_list);
404 $shader->SetMatrix(model => $model_oga); # instance of OpenGL::Shader
405 See OpenGL, OpenGL::Shader, OpenGL::Array, rotation matrix
407 <https://en.wikipedia.org/wiki/Rotation_matrix>.
408 • $new_matrix = $matrix->each( \&function );
410 Creates a new matrix by evaluating a code reference on each element
412 of the given matrix. The function is passed the element, the row
413 index and the column index, in that order. The value the function
414 returns ( or the value of the last executed statement ) is the
415 value given to the corresponding element in $new_matrix.
416 Example:
418 # add 1 to every element in the matrix
420 $matrix = $matrix->each ( sub { (shift) + 1 } );
421 Example:
423 my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
425 ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
426 : -1*$matrix->minor($i,$j)->det();
427 } );
428 This code needs some explanation. For each element of $matrix, it
430 throws away the actual value and stores the row and column indexes
431 in $i and $j. Then it sets element [$i,$j] in $cofactor to the
432 determinant of "$matrix->minor($i,$j)" if it is an "even" element,
433 or "-1*$matrix->minor($i,$j)" if it is an "odd" element.
434 • $new_matrix = $matrix->each_diag( \&function );
436 Creates a new matrix by evaluating a code reference on each
438 diagonal element of the given matrix. The function is passed the
439 element, the row index and the column index, in that order. The
440 value the function returns ( or the value of the last executed
441 statement ) is the value given to the corresponding element in
442 $new_matrix.
443 • $matrix->swap_col( $col1, $col2 );
445 This method takes two one-based column numbers and swaps the values
447 of each element in each column. "$matrix->swap_col(2,3)" would
448 replace column 2 in $matrix with column 3, and replace column 3
449 with column 2.
450 • $matrix->swap_row( $row1, $row2 );
452 This method takes two one-based row numbers and swaps the values of
454 each element in each row. "$matrix->swap_row(2,3)" would replace
455 row 2 in $matrix with row 3, and replace row 3 with row 2.
456 • $matrix->assign_row( $row_number , $new_row_vector );
458 This method takes a one-based row number and assigns row
460 $row_number of $matrix with $new_row_vector and returns the
461 resulting matrix. "$matrix->assign_row(5, $x)" would replace row 5
462 in $matrix with the row vector $x.
463 • $matrix->maximum(); and $matrix->minimum();
465 These two methods work similarly, one for computing the maximum
467 element or elements from a matrix, and the minimum element or
468 elements from a matrix. They work in a similar way as
469 Octave/MatLab max/min functions.
470 When computing the maximum or minimum from a vector (vertical or
472 horizontal), only one element is returned. When computing the
473 maximum or minimum from a matrix, the maximum/minimum element for
474 each column is returned in an array reference.
475 When called in list context, the function returns a pair, where the
477 first element is the maximum/minimum element (or elements) and the
478 second is the position of that value in the vector (first
479 occurrence), or the row where it occurs, for matrices.
480 Consider the matrix and vector below for the following examples:
482 [ 1 9 4 ]
484 $A = [ 3 5 2 ] $B = [ 8 7 9 5 3 ]
485 [ 8 7 6 ]
486 When used in scalar context:
488 $max = $A->maximum(); # $max = [ 8, 9, 6 ]
490 $min = $B->minimum(); # $min = 3
491 When used in list context:
493 ($min, $pos) = $A->minimum(); # $min = [ 1 5 2 ]
495 # $pos = [ 1 2 2 ]
496 ($max, $pos) = $B->maximum(); # $max = 9
497 # $pos = 3
498 Matrix Operations
500 • "$det = $matrix->det();"
501 Returns the determinant of the matrix, without going through the
503 rigamarole of computing a LR decomposition. This method should be
504 much faster than LR decomposition if the matrix is diagonal or
505 triangular. Otherwise, it is just a wrapper for
506 "$matrix->decompose_LR->det_LR". If the determinant is zero, there
507 is no inverse and vice-versa. Only quadratic matrices have
508 determinants.
509 • "$inverse = $matrix->inverse();"
511 Returns the inverse of a matrix, without going through the
513 rigamarole of computing a LR decomposition. If no inverse exists,
514 undef is returned and an error is printed via "carp()". This is
515 nothing but a wrapper for "$matrix->decompose_LR->invert_LR".
516 • "($rows,$columns) = $matrix->dim();"
518 Returns a list of two items, representing the number of rows and
520 columns the given matrix "$matrix" contains.
521 • "$norm_one = $matrix->norm_one();"
523 Returns the "one"-norm of the given matrix "$matrix".
525 The "one"-norm is defined as follows:
527 For each column, the sum of the absolute values of the elements in
529 the different rows of that column is calculated. Finally, the
530 maximum of these sums is returned.
531 Note that the "one"-norm and the "maximum"-norm are mathematically
533 equivalent, although for the same matrix they usually yield a
534 different value.
535 Therefore, you should only compare values that have been calculated
537 using the same norm!
538 Throughout this package, the "one"-norm is (arbitrarily) used for
540 all comparisons, for the sake of uniformity and comparability,
541 except for the iterative methods "solve_GSM()", "solve_SSM()" and
542 "solve_RM()" which use either norm depending on the matrix itself.
543 • "$norm_max = $matrix->norm_max();"
545 Returns the "maximum"-norm of the given matrix $matrix.
547 The "maximum"-norm is defined as follows:
549 For each row, the sum of the absolute values of the elements in the
551 different columns of that row is calculated. Finally, the maximum
552 of these sums is returned.
553 Note that the "maximum"-norm and the "one"-norm are mathematically
555 equivalent, although for the same matrix they usually yield a
556 different value.
557 Therefore, you should only compare values that have been calculated
559 using the same norm!
560 Throughout this package, the "one"-norm is (arbitrarily) used for
562 all comparisons, for the sake of uniformity and comparability,
563 except for the iterative methods "solve_GSM()", "solve_SSM()" and
564 "solve_RM()" which use either norm depending on the matrix itself.
565 • "$norm_sum = $matrix->norm_sum();"
567 This is a very simple norm which is defined as the sum of the
569 absolute values of every element.
570 • $p_norm = $matrix->norm_p($n);>
572 This function returns the "p-norm" of a vector. The argument $n
574 must be a number greater than or equal to 1 or the string "Inf".
575 The p-norm is defined as (sum(x_i^p))^(1/p). In words, it raised
576 each element to the p-th power, adds them up, and then takes the
577 p-th root of that number. If the string "Inf" is passed, the
578 "infinity-norm" is computed, which is really the limit of the
579 p-norm as p goes to infinity. It is defined as the maximum element
580 of the vector. Also, note that the familiar Euclidean distance
581 between two vectors is just a special case of a p-norm, when p is
582 equal to 2.
583 Example:
585 $a = Math::MatrixReal->new_from_cols([[1,2,3]]);
586 $p1 = $a->norm_p(1);
587 $p2 = $a->norm_p(2);
588 $p3 = $a->norm_p(3);
589 $pinf = $a->norm_p("Inf");
590 print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";
592 $i1 = $a->new_from_rows([[1,0]]);
594 $i2 = $a->new_from_rows([[0,1]]);
595 # this should be sqrt(2) since it is the same as the
597 # hypotenuse of a 1 by 1 right triangle
598 $dist = ($i1-$i2)->norm_p(2);
600 print "Distance is $dist, which should be " . sqrt(2) . "\n";
601 Output:
603 (1,2,3,Inf) norm:
605 6
606 3.74165738677394139
607 3.30192724889462668
608 3
609 Distance is 1.41421356237309505, which should be 1.41421356237309505
611 • $frob_norm = "$matrix->norm_frobenius();"
613 This norm is similar to that of a p-norm where p is 2, except it
615 acts on a matrix, not a vector. Each element of the matrix is
616 squared, this is added up, and then a square root is taken.
617 • "$matrix->spectral_radius();"
619 Returns the maximum value of the absolute value of all eigenvalues.
621 Currently this computes all eigenvalues, then sifts through them to
622 find the largest in absolute value. Needless to say, this is very
623 inefficient, and in the future an algorithm that computes only the
624 largest eigenvalue may be implemented.
625 • "$matrix1->transpose($matrix2);"
627 Calculates the transposed matrix of matrix $matrix2 and stores the
629 result in matrix "$matrix1" (which must already exist and have the
630 same size as matrix "$matrix2"!).
631 This operation can also be carried out "in-place", i.e., input and
633 output matrix may be identical.
634 Transposition is a symmetry operation: imagine you rotate the
636 matrix along the axis of its main diagonal (going through elements
637 (1,1), (2,2), (3,3) and so on) by 180 degrees.
638 Another way of looking at it is to say that rows and columns are
640 swapped. In fact the contents of element "(i,j)" are swapped with
641 those of element "(j,i)".
642 Note that (especially for vectors) it makes a big difference if you
644 have a row vector, like this:
645 [ -1 0 1 ]
647 or a column vector, like this:
649 [ -1 ]
651 [ 0 ]
652 [ 1 ]
653 the one vector being the transposed of the other!
655 This is especially true for the matrix product of two vectors:
657 [ -1 ]
659 [ -1 0 1 ] * [ 0 ] = [ 2 ] , whereas
660 [ 1 ]
661 * [ -1 0 1 ]
663 [ -1 ] [ 1 0 -1 ]
664 [ 0 ] * [ -1 0 1 ] = [ -1 ] [ 1 0 -1 ] = [ 0 0 0 ]
665 [ 1 ] [ 0 ] [ 0 0 0 ] [ -1 0 1 ]
666 [ 1 ] [ -1 0 1 ]
667 So be careful about what you really mean!
669 Hint: throughout this module, whenever a vector is explicitly
671 required for input, a COLUMN vector is expected!
672 • "$trace = $matrix->trace();"
674 This returns the trace of the matrix, which is defined as the sum
676 of the diagonal elements. The matrix must be quadratic.
677 • "$minor = $matrix->minor($row,$col);"
679 Returns the minor matrix corresponding to $row and $col. $matrix
681 must be quadratic. If $matrix is n rows by n cols, the minor of
682 $row and $col will be an (n-1) by (n-1) matrix. The minor is
683 defined as crossing out the row and the col specified and returning
684 the remaining rows and columns as a matrix. This method is used by
685 "cofactor()".
686 • "$cofactor = $matrix->cofactor();"
688 The cofactor matrix is constructed as follows:
690 For each element, cross out the row and column that it sits in.
692 Now, take the determinant of the matrix that is left in the other
693 rows and columns. Multiply the determinant by (-1)^(i+j), where i
694 is the row index, and j is the column index. Replace the given
695 element with this value.
696 The cofactor matrix can be used to find the inverse of the matrix.
698 One formula for the inverse of a matrix is the cofactor matrix
699 transposed divided by the original determinant of the matrix.
700 The following two inverses should be exactly the same:
702 my $inverse1 = $matrix->inverse;
704 my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
705 Caveat: Although the cofactor matrix is simple algorithm to compute
707 the inverse of a matrix, and can be used with pencil and paper for
708 small matrices, it is comically slower than the native "inverse()"
709 function. Here is a small benchmark:
710 # $matrix1 is 15x15
712 $det = $matrix1->det;
713 timethese( 10,
714 {'inverse' => sub { $matrix1->inverse(); },
715 'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
716 } );
717 Benchmark: timing 10 iterations of LR, cofactor, inverse...
720 inverse: 1 wallclock secs ( 0.56 usr + 0.00 sys = 0.56 CPU) @ 17.86/s (n=10)
721 cofactor: 36 wallclock secs (36.62 usr + 0.01 sys = 36.63 CPU) @ 0.27/s (n=10)
722 • "$adjoint = $matrix->adjoint();"
724 The adjoint is just the transpose of the cofactor matrix. This
726 method is just an alias for " ~($matrix->cofactor)".
727 • "$part_of_matrix = $matrix->submatrix(x1,y1,x2,Y2);"
729 Submatrix permit to select only part of existing matrix in order to
731 produce a new one. This method take four arguments to define a
732 selection area:
733 - firstly: Coordinate of top left corner to select (x1,y1)
735 - secondly: Coordinate of bottom right corner to select (x2,y2)
736 Example:
738 my $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
740 [ 0 0 0 0 0 0 0 ]
741 [ 0 0 0 0 0 0 0 ]
742 [ 0 0 0 0 0 0 0 ]
743 [ 0 0 0 0 0 0 0 ]
744 [ 0 0 0 0 1 0 1 ]
745 [ 0 0 0 0 0 1 0 ]
746 [ 0 0 0 0 1 0 1 ]
747 MATRIX
748 my $submatrix = $matrix->submatrix(5,5,7,7);
750 $submatrix->display_precision(0);
751 print $submatrix;
752 Output:
754 [ 1 0 1 ]
756 [ 0 1 0 ]
757 [ 1 0 1 ]
758 Arithmetic Operations
760 • "$matrix1->add($matrix2,$matrix3);"
761 Calculates the sum of matrix "$matrix2" and matrix "$matrix3" and
763 stores the result in matrix "$matrix1" (which must already exist
764 and have the same size as matrix "$matrix2" and matrix
765 "$matrix3"!).
766 This operation can also be carried out "in-place", i.e., the output
768 and one (or both) of the input matrices may be identical.
769 • "$matrix1->subtract($matrix2,$matrix3);"
771 Calculates the difference of matrix "$matrix2" minus matrix
773 "$matrix3" and stores the result in matrix "$matrix1" (which must
774 already exist and have the same size as matrix "$matrix2" and
775 matrix "$matrix3"!).
776 This operation can also be carried out "in-place", i.e., the output
778 and one (or both) of the input matrices may be identical.
779 Note that this operation is the same as
781 "$matrix1->add($matrix2,-$matrix3);", although the latter is a
782 little less efficient.
783 • "$matrix1->multiply_scalar($matrix2,$scalar);"
785 Calculates the product of matrix "$matrix2" and the number
787 "$scalar" (i.e., multiplies each element of matrix "$matrix2" with
788 the factor "$scalar") and stores the result in matrix "$matrix1"
789 (which must already exist and have the same size as matrix
790 "$matrix2"!).
791 This operation can also be carried out "in-place", i.e., input and
793 output matrix may be identical.
794 • "$product_matrix = $matrix1->multiply($matrix2);"
796 Calculates the product of matrix "$matrix1" and matrix "$matrix2"
798 and returns an object reference to a new matrix "$product_matrix"
799 in which the result of this operation has been stored.
800 Note that the dimensions of the two matrices "$matrix1" and
802 "$matrix2" (i.e., their numbers of rows and columns) must harmonize
803 in the following way (example):
804 [ 2 2 ]
806 [ 2 2 ]
807 [ 2 2 ]
808 [ 1 1 1 ] [ * * ]
810 [ 1 1 1 ] [ * * ]
811 [ 1 1 1 ] [ * * ]
812 [ 1 1 1 ] [ * * ]
813 I.e., the number of columns of matrix "$matrix1" has to be the same
815 as the number of rows of matrix "$matrix2".
816 The number of rows and columns of the resulting matrix
818 "$product_matrix" is determined by the number of rows of matrix
819 "$matrix1" and the number of columns of matrix "$matrix2",
820 respectively.
821 • "$matrix1->negate($matrix2);"
823 Calculates the negative of matrix "$matrix2" (i.e., multiplies all
825 elements with "-1") and stores the result in matrix "$matrix1"
826 (which must already exist and have the same size as matrix
827 "$matrix2"!).
828 This operation can also be carried out "in-place", i.e., input and
830 output matrix may be identical.
831 • "$matrix_to_power = $matrix1->exponent($integer);"
833 Raises the matrix to the $integer power. Obviously, $integer must
835 be an integer. If it is zero, the identity matrix is returned. If a
836 negative integer is given, the inverse will be computed (if it
837 exists) and then raised the the absolute value of $integer. The
838 matrix must be quadratic.
839 Boolean Matrix Operations
841 • "$matrix->is_quadratic();"
842 Returns a boolean value indicating if the given matrix is quadratic
844 (also know as "square" or "n by n"). A matrix is quadratic if it
845 has the same number of rows as it does columns.
846 • "$matrix->is_square();"
848 This is an alias for "is_quadratic()".
850 • "$matrix->is_symmetric();"
852 Returns a boolean value indicating if the given matrix is
854 symmetric. By definition, a matrix is symmetric if and only if
855 (M[i,j]=M[j,i]). This is equivalent to "($matrix == ~$matrix)" but
856 without memory allocation. Only quadratic matrices can be
857 symmetric.
858 Notes: A symmetric matrix always has real eigenvalues/eigenvectors.
860 A matrix plus its transpose is always symmetric.
861 • "$matrix->is_skew_symmetric();"
863 Returns a boolean value indicating if the given matrix is skew
865 symmetric. By definition, a matrix is symmetric if and only if
866 (M[i,j]=-M[j,i]). This is equivalent to "($matrix == -(~$matrix))"
867 but without memory allocation. Only quadratic matrices can be skew
868 symmetric.
869 • "$matrix->is_diagonal();"
871 Returns a boolean value indicating if the given matrix is diagonal,
873 i.e. all of the nonzero elements are on the main diagonal. Only
874 quadratic matrices can be diagonal.
875 • "$matrix->is_tridiagonal();"
877 Returns a boolean value indicating if the given matrix is
879 tridiagonal, i.e. all of the nonzero elements are on the main
880 diagonal or the diagonals above and below the main diagonal. Only
881 quadratic matrices can be tridiagonal.
882 • "$matrix->is_upper_triangular();"
884 Returns a boolean value indicating if the given matrix is upper
886 triangular, i.e. all of the nonzero elements not on the main
887 diagonal are above it. Only quadratic matrices can be upper
888 triangular. Note: diagonal matrices are both upper and lower
889 triangular.
890 • "$matrix->is_lower_triangular();"
892 Returns a boolean value indicating if the given matrix is lower
894 triangular, i.e. all of the nonzero elements not on the main
895 diagonal are below it. Only quadratic matrices can be lower
896 triangular. Note: diagonal matrices are both upper and lower
897 triangular.
898 • "$matrix->is_orthogonal();"
900 Returns a boolean value indicating if the given matrix is
902 orthogonal. An orthogonal matrix is has the property that the
903 transpose equals the inverse of the matrix. Instead of computing
904 each and comparing them, this method multiplies the matrix by it's
905 transpose, and returns true if this turns out to be the identity
906 matrix, false otherwise. Only quadratic matrices can orthogonal.
907 • "$matrix->is_binary();"
909 Returns a boolean value indicating if the given matrix is binary.
911 A matrix is binary if it contains only zeroes or ones.
912 • "$matrix->is_gramian();"
914 Returns a boolean value indicating if the give matrix is Gramian.
916 A matrix $A is Gramian if and only if there exists a square matrix
917 $B such that "$A = ~$B*$B". This is equivalent to checking if $A is
918 symmetric and has all nonnegative eigenvalues, which is what
919 Math::MatrixReal uses to check for this property.
920 • "$matrix->is_LR();"
922 Returns a boolean value indicating if the matrix is an LR
924 decomposition matrix.
925 • "$matrix->is_positive();"
927 Returns a boolean value indicating if the matrix contains only
929 positive entries. Note that a zero entry is not positive and will
930 cause "is_positive()" to return false.
931 • "$matrix->is_negative();"
933 Returns a boolean value indicating if the matrix contains only
935 negative entries. Note that a zero entry is not negative and will
936 cause "is_negative()" to return false.
937 • "$matrix->is_periodic($k);"
939 Returns a boolean value indicating if the matrix is periodic with
941 period $k. This is true if "$matrix ** ($k+1) == $matrix". When
942 "$k == 1", this reduces down to the "is_idempotent()" function.
943 • "$matrix->is_idempotent();"
945 Returns a boolean value indicating if the matrix is idempotent,
947 which is defined as the square of the matrix being equal to the
948 original matrix, i.e "$matrix ** 2 == $matrix".
949 • "$matrix->is_row_vector();"
951 Returns a boolean value indicating if the matrix is a row vector.
953 A row vector is a matrix which is 1xn. Note that the 1x1 matrix is
954 both a row and column vector.
955 • "$matrix->is_col_vector();"
957 Returns a boolean value indicating if the matrix is a col vector.
959 A col vector is a matrix which is nx1. Note that the 1x1 matrix is
960 both a row and column vector.
961 Eigensystems
963 • "($l, $V) = $matrix->sym_diagonalize();"
964 This method performs the diagonalization of the quadratic symmetric
966 matrix M stored in $matrix. On output, l is a column vector
967 containing all the eigenvalues of M and V is an orthogonal matrix
968 which columns are the corresponding normalized eigenvectors. The
969 primary property of an eigenvalue l and an eigenvector x is of course
970 that: M * x = l * x.
971 The method uses a Householder reduction to tridiagonal form followed
973 by a QL algoritm with implicit shifts on this tridiagonal. (The
974 tridiagonal matrix is kept internally in a compact form in this
975 routine to save memory.) In fact, this routine wraps the
976 householder() and tri_diagonalize() methods described below when
977 their intermediate results are not desired. The overall algorithmic
978 complexity of this technique is O(N^3). According to several books,
979 the coefficient hidden by the 'O' is one of the best possible for
980 general (symmetric) matrixes.
981 • "($T, $Q) = $matrix->householder();"
983 This method performs the Householder algorithm which reduces the n by
985 n real symmetric matrix M contained in $matrix to tridiagonal form.
986 On output, T is a symmetric tridiagonal matrix (only diagonal and
987 off-diagonal elements are non-zero) and Q is an orthogonal matrix
988 performing the tranformation between M and T ("$M == $Q * $T * ~$Q").
989 • "($l, $V) = $T->tri_diagonalize([$Q]);"
991 This method diagonalizes the symmetric tridiagonal matrix T. On
993 output, $l and $V are similar to the output values described for
994 sym_diagonalize().
995 The optional argument $Q corresponds to an orthogonal transformation
997 matrix Q that should be used additionally during V (eigenvectors)
998 computation. It should be supplied if the desired eigenvectors
999 correspond to a more general symmetric matrix M previously reduced by
1000 the householder() method, not a mere tridiagonal. If T is really a
1001 tridiagonal matrix, Q can be omitted (it will be internally created
1002 in fact as an identity matrix). The method uses a QL algorithm (with
1003 implicit shifts).
1004 • "$l = $matrix->sym_eigenvalues();"
1006 This method computes the eigenvalues of the quadratic symmetric
1008 matrix M stored in $matrix. On output, l is a column vector
1009 containing all the eigenvalues of M. Eigenvectors are not computed
1010 (on the contrary of "sym_diagonalize()") and this method is more
1011 efficient (even though it uses a similar algorithm with two phases).
1012 However, understand that the algorithmic complexity of this technique
1013 is still also O(N^3). But the coefficient hidden by the 'O' is better
1014 by a factor of..., well, see your benchmark, it's wiser.
1015 This routine wraps the householder_tridiagonal() and
1017 tri_eigenvalues() methods described below when the intermediate
1018 tridiagonal matrix is not needed.
1019 • "$T = $matrix->householder_tridiagonal();"
1021 This method performs the Householder algorithm which reduces the n by
1023 n real symmetric matrix M contained in $matrix to tridiagonal form.
1024 On output, T is the obtained symmetric tridiagonal matrix (only
1025 diagonal and off-diagonal elements are non-zero). The operation is
1026 similar to the householder() method, but potentially a little more
1027 efficient as the transformation matrix is not computed.
1028 • $l = $T->tri_eigenvalues();
1030 This method computesthe eigenvalues of the symmetric tridiagonal
1032 matrix T. On output, $l is a vector containing the eigenvalues
1033 (similar to "sym_eigenvalues()"). This method is much more efficient
1034 than tri_diagonalize() when eigenvectors are not needed.
1035 Miscellaneous
1037 • $matrix->zero();
1038 Assigns a zero to every element of the matrix "$matrix", i.e.,
1040 erases all values previously stored there, thereby effectively
1041 transforming the matrix into a "zero"-matrix or "null"-matrix, the
1042 neutral element of the addition operation in a Ring.
1043 (For instance the (quadratic) matrices with "n" rows and columns
1045 and matrix addition and multiplication form a Ring. Most prominent
1046 characteristic of a Ring is that multiplication is not commutative,
1047 i.e., in general, ""matrix1 * matrix2"" is not the same as
1048 ""matrix2 * matrix1""!)
1049 • $matrix->one();
1051 Assigns one's to the elements on the main diagonal (elements (1,1),
1053 (2,2), (3,3) and so on) of matrix "$matrix" and zero's to all
1054 others, thereby erasing all values previously stored there and
1055 transforming the matrix into a "one"-matrix, the neutral element of
1056 the multiplication operation in a Ring.
1057 (If the matrix is quadratic (which this method doesn't require,
1059 though), then multiplying this matrix with itself yields this same
1060 matrix again, and multiplying it with some other matrix leaves that
1061 other matrix unchanged!)
1062 • "$latex_string = $matrix->as_latex( align=> "c", format => "%s",
1064 name => "" );"
1065 This function returns the matrix as a LaTeX string. It takes a hash
1067 as an argument which is used to control the style of the output.
1068 The hash element "align" may be "c","l" or "r", corresponding to
1069 center, left and right, respectively. The "format" element is a
1070 format string that is given to "sprintf" to control the style of
1071 number format, such a floating point or scientific notation. The
1072 "name" element can be used so that a LaTeX string of "$name = " is
1073 prepended to the string.
1074 Example:
1076 my $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
1078 print $a->as_latex( ( format => "%.2f", align => "l",name => "A" ) );
1079 Output:
1081 $A = $ $
1082 \left( \begin{array}{ll}
1083 1.23&1.00 \\
1084 5.68&2.00 \\
1085 9.10&3.00
1086 \end{array} \right)
1087 $
1088 • "$yacas_string = $matrix->as_yacas( format => "%s", name => "",
1090 semi => 0 );"
1091 This function returns the matrix as a string that can be read by
1093 Yacas. It takes a hash as an an argument which controls the style
1094 of the output. The "format" element is a format string that is
1095 given to "sprintf" to control the style of number format, such a
1096 floating point or scientific notation. The "name" element can be
1097 used so that "$name = " is prepended to the string. The <semi>
1098 element can be set to 1 to that a semicolon is appended (so Matlab
1099 does not print out the matrix.)
1100 Example:
1102 $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
1104 print $a->as_yacas( ( format => "%.2f", align => "l",name => "A" ) );
1105 Output:
1107 A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
1109 • "$matlab_string = $matrix->as_matlab( format => "%s", name => "",
1111 semi => 0 );"
1112 This function returns the matrix as a string that can be read by
1114 Matlab. It takes a hash as an an argument which controls the style
1115 of the output. The "format" element is a format string that is
1116 given to "sprintf" to control the style of number format, such a
1117 floating point or scientific notation. The "name" element can be
1118 used so that "$name = " is prepended to the string. The <semi>
1119 element can be set to 1 to that a semicolon is appended (so Matlab
1120 does not print out the matrix.)
1121 Example:
1123 my $a = Math::MatrixReal->new_from_rows([[ 1.234, 5.678, 9.1011],[1,2,3]] );
1125 print $a->as_matlab( ( format => "%.3f", name => "A",semi => 1 ) );
1126 Output:
1128 A = [ 1.234 5.678 9.101;
1129 1.000 2.000 3.000];
1130 • "$scilab_string = $matrix->as_scilab( format => "%s", name => "",
1132 semi => 0 );"
1133 This function is just an alias for "as_matlab()", since both Scilab
1135 and Matlab have the same matrix format.
1136 • "$minimum = Math::MatrixReal::min($number1,$number2);" "$minimum =
1138 Math::MatrixReal::min($matrix);" "<$minimum = $matrix-"min;>>
1139 Returns the minimum of the two numbers ""number1"" and ""number2""
1141 if called with two arguments, or returns the value of the smallest
1142 element of a matrix if called with one argument or as an object
1143 method.
1144 • "$maximum = Math::MatrixReal::max($number1,$number2);" "$maximum =
1146 Math::MatrixReal::max($number1,$number2);" "$maximum =
1147 Math::MatrixReal::max($matrix);" "<$maximum = $matrix-"max;>>
1148 Returns the maximum of the two numbers ""number1"" and ""number2""
1150 if called with two arguments, or returns the value of the largest
1151 element of a matrix if called with one arguemnt or as on object
1152 method.
1153 • "$minimal_cost_matrix = $cost_matrix->kleene();"
1155 Copies the matrix "$cost_matrix" (which has to be quadratic!) to a
1157 new matrix of the same size (i.e., "clones" the input matrix) and
1158 applies Kleene's algorithm to it.
1159 See Math::Kleene(3) for more details about this algorithm!
1161 The method returns an object reference to the new matrix.
1163 Matrix "$cost_matrix" is not changed by this method in any way.
1165 • "($norm_matrix,$norm_vector) = $matrix->normalize($vector);"
1167 This method is used to improve the numerical stability when solving
1169 linear equation systems.
1170 Suppose you have a matrix "A" and a vector "b" and you want to find
1172 out a vector "x" so that "A * x = b", i.e., the vector "x" which
1173 solves the equation system represented by the matrix "A" and the
1174 vector "b".
1175 Applying this method to the pair (A,b) yields a pair (A',b') where
1177 each row has been divided by (the absolute value of) the greatest
1178 coefficient appearing in that row. So this coefficient becomes
1179 equal to "1" (or "-1") in the new pair (A',b') (all others become
1180 smaller than one and greater than minus one).
1181 Note that this operation does not change the equation system itself
1183 because the same division is carried out on either side of the
1184 equation sign!
1185 The method requires a quadratic (!) matrix "$matrix" and a vector
1187 "$vector" for input (the vector must be a column vector with the
1188 same number of rows as the input matrix) and returns a list of two
1189 items which are object references to a new matrix and a new vector,
1190 in this order.
1191 The output matrix and vector are clones of the input matrix and
1193 vector to which the operation explained above has been applied.
1194 The input matrix and vector are not changed by this in any way.
1196 Example of how this method can affect the result of the methods to
1198 solve equation systems (explained immediately below following this
1199 method):
1200 Consider the following little program:
1202 #!perl -w
1204 use Math::MatrixReal qw(new_from_string);
1206 $A = Math::MatrixReal->new_from_string(<<"MATRIX");
1208 [ 1 2 3 ]
1209 [ 5 7 11 ]
1210 [ 23 19 13 ]
1211 MATRIX
1212 $b = Math::MatrixReal->new_from_string(<<"MATRIX");
1214 [ 0 ]
1215 [ 1 ]
1216 [ 29 ]
1217 MATRIX
1218 $LR = $A->decompose_LR();
1220 if (($dim,$x,$B) = $LR->solve_LR($b))
1221 {
1222 $test = $A * $x;
1223 print "x = \n$x";
1224 print "A * x = \n$test";
1225 }
1226 ($A_,$b_) = $A->normalize($b);
1228 $LR = $A_->decompose_LR();
1230 if (($dim,$x,$B) = $LR->solve_LR($b_))
1231 {
1232 $test = $A * $x;
1233 print "x = \n$x";
1234 print "A * x = \n$test";
1235 }
1236 This will print:
1238 x =
1240 [ 1.000000000000E+00 ]
1241 [ 1.000000000000E+00 ]
1242 [ -1.000000000000E+00 ]
1243 A * x =
1244 [ 4.440892098501E-16 ]
1245 [ 1.000000000000E+00 ]
1246 [ 2.900000000000E+01 ]
1247 x =
1248 [ 1.000000000000E+00 ]
1249 [ 1.000000000000E+00 ]
1250 [ -1.000000000000E+00 ]
1251 A * x =
1252 [ 0.000000000000E+00 ]
1253 [ 1.000000000000E+00 ]
1254 [ 2.900000000000E+01 ]
1255 You can see that in the second example (where "normalize()" has
1257 been used), the result is "better", i.e., more accurate!
1258 • "$LR_matrix = $matrix->decompose_LR();"
1260 This method is needed to solve linear equation systems.
1262 Suppose you have a matrix "A" and a vector "b" and you want to find
1264 out a vector "x" so that "A * x = b", i.e., the vector "x" which
1265 solves the equation system represented by the matrix "A" and the
1266 vector "b".
1267 You might also have a matrix "A" and a whole bunch of different
1269 vectors "b1".."bk" for which you need to find vectors "x1".."xk" so
1270 that "A * xi = bi", for "i=1..k".
1271 Using Gaussian transformations (multiplying a row or column with a
1273 factor, swapping two rows or two columns and adding a multiple of
1274 one row or column to another), it is possible to decompose any
1275 matrix "A" into two triangular matrices, called "L" and "R" (for
1276 "Left" and "Right").
1277 "L" has one's on the main diagonal (the elements (1,1), (2,2),
1279 (3,3) and so so), non-zero values to the left and below of the main
1280 diagonal and all zero's in the upper right half of the matrix.
1281 "R" has non-zero values on the main diagonal as well as to the
1283 right and above of the main diagonal and all zero's in the lower
1284 left half of the matrix, as follows:
1285 [ 1 0 0 0 0 ] [ x x x x x ]
1287 [ x 1 0 0 0 ] [ 0 x x x x ]
1288 L = [ x x 1 0 0 ] R = [ 0 0 x x x ]
1289 [ x x x 1 0 ] [ 0 0 0 x x ]
1290 [ x x x x 1 ] [ 0 0 0 0 x ]
1291 Note that ""L * R"" is equivalent to matrix "A" in the sense that
1293 "L * R * x = b <==> A * x = b" for all vectors "x", leaving out
1294 of account permutations of the rows and columns (these are taken
1295 care of "magically" by this module!) and numerical errors.
1296 Trick:
1298 Because we know that "L" has one's on its main diagonal, we can
1300 store both matrices together in the same array without information
1301 loss! I.e.,
1302 [ R R R R R ]
1304 [ L R R R R ]
1305 LR = [ L L R R R ]
1306 [ L L L R R ]
1307 [ L L L L R ]
1308 Beware, though, that "LR" and ""L * R"" are not the same!!!
1310 Note also that for the same reason, you cannot apply the method
1312 "normalize()" to an "LR" decomposition matrix. Trying to do so will
1313 yield meaningless rubbish!
1314 (You need to apply "normalize()" to each pair (Ai,bi) BEFORE
1316 decomposing the matrix "Ai'"!)
1317 Now what does all this help us in solving linear equation systems?
1319 It helps us because a triangular matrix is the next best thing that
1321 can happen to us besides a diagonal matrix (a matrix that has non-
1322 zero values only on its main diagonal - in which case the solution
1323 is trivial, simply divide ""b[i]"" by ""A[i,i]"" to get ""x[i]""!).
1324 To find the solution to our problem ""A * x = b"", we divide this
1326 problem in parts: instead of solving "A * x = b" directly, we first
1327 decompose "A" into "L" and "R" and then solve ""L * y = b"" and
1328 finally ""R * x = y"" (motto: divide and rule!).
1329 From the illustration above it is clear that solving ""L * y = b""
1331 and ""R * x = y"" is straightforward: we immediately know that
1332 "y[1] = b[1]". We then deduce swiftly that
1333 y[2] = b[2] - L[2,1] * y[1]
1335 (and we know ""y[1]"" by now!), that
1337 y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
1339 and so on.
1341 Having effortlessly calculated the vector "y", we now proceed to
1343 calculate the vector "x" in a similar fashion: we see immediately
1344 that "x[n] = y[n] / R[n,n]". It follows that
1345 x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
1347 and
1349 x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
1351 / R[n-2,n-2]
1352 and so on.
1354 You can see that - especially when you have many vectors "b1".."bk"
1356 for which you are searching solutions to "A * xi = bi" - this
1357 scheme is much more efficient than a straightforward, "brute force"
1358 approach.
1359 This method requires a quadratic matrix as its input matrix.
1361 If you don't have that many equations, fill up with zero's (i.e.,
1363 do nothing to fill the superfluous rows if it's a "fresh" matrix,
1364 i.e., a matrix that has been created with "new()" or "shadow()").
1365 The method returns an object reference to a new matrix containing
1367 the matrices "L" and "R".
1368 The input matrix is not changed by this method in any way.
1370 Note that you can "copy()" or "clone()" the result of this method
1372 without losing its "magical" properties (for instance concerning
1373 the hidden permutations of its rows and columns).
1374 However, as soon as you are applying any method that alters the
1376 contents of the matrix, its "magical" properties are stripped off,
1377 and the matrix immediately reverts to an "ordinary" matrix (with
1378 the values it just happens to contain at that moment, be they
1379 meaningful as an ordinary matrix or not!).
1380 • "($dimension,$x_vector,$base_matrix) =
1382 $LR_matrix""->""solve_LR($b_vector);"
1383 Use this method to actually solve an equation system.
1385 Matrix "$LR_matrix" must be a (quadratic) matrix returned by the
1387 method "decompose_LR()", the LR decomposition matrix of the matrix
1388 "A" of your equation system "A * x = b".
1389 The input vector "$b_vector" is the vector "b" in your equation
1391 system "A * x = b", which must be a column vector and have the same
1392 number of rows as the input matrix "$LR_matrix".
1393 The method returns a list of three items if a solution exists or an
1395 empty list otherwise (!).
1396 Therefore, you should always use this method like this:
1398 if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
1400 {
1401 # do something with the solution...
1402 }
1403 else
1404 {
1405 # do something with the fact that there is no solution...
1406 }
1407 The three items returned are: the dimension "$dimension" of the
1409 solution space (which is zero if only one solution exists, one if
1410 the solution is a straight line, two if the solution is a plane,
1411 and so on), the solution vector "$x_vector" (which is the vector
1412 "x" of your equation system "A * x = b") and a matrix
1413 "$base_matrix" representing a base of the solution space (a set of
1414 vectors which put up the solution space like the spokes of an
1415 umbrella).
1416 Only the first "$dimension" columns of this base matrix actually
1418 contain entries, the remaining columns are all zero.
1419 Now what is all this stuff with that "base" good for?
1421 The output vector "x" is ALWAYS a solution of your equation system
1423 "A * x = b".
1424 But also any vector "$vector"
1426 $vector = $x_vector->clone();
1428 $machine_infinity = 1E+99; # or something like that
1430 for ( $i = 1; $i <= $dimension; $i++ )
1432 {
1433 $vector += rand($machine_infinity) * $base_matrix->column($i);
1434 }
1435 is a solution to your problem "A * x = b", i.e., if "$A_matrix"
1437 contains your matrix "A", then
1438 print abs( $A_matrix * $vector - $b_vector ), "\n";
1440 should print a number around 1E-16 or so!
1442 By the way, note that you can actually calculate those vectors
1444 "$vector" a little more efficient as follows:
1445 $rand_vector = $x_vector->shadow();
1447 $machine_infinity = 1E+99; # or something like that
1449 for ( $i = 1; $i <= $dimension; $i++ )
1451 {
1452 $rand_vector->assign($i,1, rand($machine_infinity) );
1453 }
1454 $vector = $x_vector + ( $base_matrix * $rand_vector );
1456 Note that the input matrix and vector are not changed by this
1458 method in any way.
1459 • "$inverse_matrix = $LR_matrix->invert_LR();"
1461 Use this method to calculate the inverse of a given matrix
1463 "$LR_matrix", which must be a (quadratic) matrix returned by the
1464 method "decompose_LR()".
1465 The method returns an object reference to a new matrix of the same
1467 size as the input matrix containing the inverse of the matrix that
1468 you initially fed into "decompose_LR()" IF THE INVERSE EXISTS, or
1469 an empty list otherwise.
1470 Therefore, you should always use this method in the following way:
1472 if ( $inverse_matrix = $LR->invert_LR() )
1474 {
1475 # do something with the inverse matrix...
1476 }
1477 else
1478 {
1479 # do something with the fact that there is no inverse matrix...
1480 }
1481 Note that by definition (disregarding numerical errors), the
1483 product of the initial matrix and its inverse (or vice-versa) is
1484 always a matrix containing one's on the main diagonal (elements
1485 (1,1), (2,2), (3,3) and so on) and zero's elsewhere.
1486 The input matrix is not changed by this method in any way.
1488 • "$condition = $matrix->condition($inverse_matrix);"
1490 In fact this method is just a shortcut for
1492 abs($matrix) * abs($inverse_matrix)
1494 Both input matrices must be quadratic and have the same size, and
1496 the result is meaningful only if one of them is the inverse of the
1497 other (for instance, as returned by the method "invert_LR()").
1498 The number returned is a measure of the "condition" of the given
1500 matrix "$matrix", i.e., a measure of the numerical stability of the
1501 matrix.
1502 This number is always positive, and the smaller its value, the
1504 better the condition of the matrix (the better the stability of all
1505 subsequent computations carried out using this matrix).
1506 Numerical stability means for example that if
1508 abs( $vec_correct - $vec_with_error ) < $epsilon
1510 holds, there must be a "$delta" which doesn't depend on the vector
1512 "$vec_correct" (nor "$vec_with_error", by the way) so that
1513 abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
1515 also holds.
1517 • "$determinant = $LR_matrix->det_LR();"
1519 Calculates the determinant of a matrix, whose LR decomposition
1521 matrix "$LR_matrix" must be given (which must be a (quadratic)
1522 matrix returned by the method "decompose_LR()").
1523 In fact the determinant is a by-product of the LR decomposition: It
1525 is (in principle, that is, except for the sign) simply the product
1526 of the elements on the main diagonal (elements (1,1), (2,2), (3,3)
1527 and so on) of the LR decomposition matrix.
1528 (The sign is taken care of "magically" by this module)
1530 • "$order = $LR_matrix->order_LR();"
1532 Calculates the order (called "Rang" in German) of a matrix, whose
1534 LR decomposition matrix "$LR_matrix" must be given (which must be a
1535 (quadratic) matrix returned by the method "decompose_LR()").
1536 This number is a measure of the number of linear independent row
1538 and column vectors (= number of linear independent equations in the
1539 case of a matrix representing an equation system) of the matrix
1540 that was initially fed into "decompose_LR()".
1541 If "n" is the number of rows and columns of the (quadratic!)
1543 matrix, then "n - order" is the dimension of the solution space of
1544 the associated equation system.
1545 • "$rank = $LR_matrix->rank_LR();"
1547 This is an alias for the "order_LR()" function. The "order" is
1549 usually called the "rank" in the United States.
1550 • "$scalar_product = $vector1->scalar_product($vector2);"
1552 Returns the scalar product of vector "$vector1" and vector
1554 "$vector2".
1555 Both vectors must be column vectors (i.e., a matrix having several
1557 rows but only one column).
1558 This is a (more efficient!) shortcut for
1560 $temp = ~$vector1 * $vector2;
1562 $scalar_product = $temp->element(1,1);
1563 or the sum "i=1..n" of the products "vector1[i] * vector2[i]".
1565 Provided none of the two input vectors is the null vector, then the
1567 two vectors are orthogonal, i.e., have an angle of 90 degrees
1568 between them, exactly when their scalar product is zero, and vice-
1569 versa.
1570 • "$vector_product = $vector1->vector_product($vector2);"
1572 Returns the vector product of vector "$vector1" and vector
1574 "$vector2".
1575 Both vectors must be column vectors (i.e., a matrix having several
1577 rows but only one column).
1578 Currently, the vector product is only defined for 3 dimensions
1580 (i.e., vectors with 3 rows); all other vectors trigger an error
1581 message.
1582 In 3 dimensions, the vector product of two vectors "x" and "y" is
1584 defined as
1585 | x[1] y[1] e[1] |
1587 determinant | x[2] y[2] e[2] |
1588 | x[3] y[3] e[3] |
1589 where the ""x[i]"" and ""y[i]"" are the components of the two
1591 vectors "x" and "y", respectively, and the ""e[i]"" are unity
1592 vectors (i.e., vectors with a length equal to one) with a one in
1593 row "i" and zero's elsewhere (this means that you have numbers and
1594 vectors as elements in this matrix!).
1595 This determinant evaluates to the rather simple formula
1597 z[1] = x[2] * y[3] - x[3] * y[2]
1599 z[2] = x[3] * y[1] - x[1] * y[3]
1600 z[3] = x[1] * y[2] - x[2] * y[1]
1601 A characteristic property of the vector product is that the
1603 resulting vector is orthogonal to both of the input vectors (if
1604 neither of both is the null vector, otherwise this is trivial),
1605 i.e., the scalar product of each of the input vectors with the
1606 resulting vector is always zero.
1607 • "$length = $vector->length();"
1609 This is actually a shortcut for
1611 $length = sqrt( $vector->scalar_product($vector) );
1613 and returns the length of a given column or row vector "$vector".
1615 Note that the "length" calculated by this method is in fact the
1617 "two"-norm (also know as the Euclidean norm) of a vector "$vector"!
1618 The general definition for norms of vectors is the following:
1620 sub vector_norm
1622 {
1623 croak "Usage: \$norm = \$vector->vector_norm(\$n);"
1624 if (@_ != 2);
1625 my($vector,$n) = @_;
1627 my($rows,$cols) = ($vector->[1],$vector->[2]);
1628 my($k,$comp,$sum);
1629 croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
1631 unless ($cols == 1);
1632 croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
1634 unless ($n > 0);
1635 croak "Math::MatrixReal::vector_norm(): norm index must be integer"
1637 unless ($n == int($n));
1638 $sum = 0;
1640 for ( $k = 0; $k < $rows; $k++ )
1641 {
1642 $comp = abs( $vector->[0][$k][0] );
1643 $sum += $comp ** $n;
1644 }
1645 return( $sum ** (1 / $n) );
1646 }
1647 Note that the case "n = 1" is the "one"-norm for matrices applied
1649 to a vector, the case "n = 2" is the euclidian norm or length of a
1650 vector, and if "n" goes to infinity, you have the "infinity"- or
1651 "maximum"-norm for matrices applied to a vector!
1652 • "$xn_vector = $matrix->""solve_GSM($x0_vector,$b_vector,$epsilon);"
1654 • "$xn_vector = $matrix->""solve_SSM($x0_vector,$b_vector,$epsilon);"
1656 • "$xn_vector =
1658 $matrix->""solve_RM($x0_vector,$b_vector,$weight,$epsilon);"
1659 In some cases it might not be practical or desirable to solve an
1661 equation system ""A * x = b"" using an analytical algorithm like
1662 the "decompose_LR()" and "solve_LR()" method pair.
1663 In fact in some cases, due to the numerical properties (the
1665 "condition") of the matrix "A", the numerical error of the obtained
1666 result can be greater than by using an approximative (iterative)
1667 algorithm like one of the three implemented here.
1668 All three methods, GSM ("Global Step Method" or
1670 "Gesamtschrittverfahren"), SSM ("Single Step Method" or
1671 "Einzelschrittverfahren") and RM ("Relaxation Method" or
1672 "Relaxationsverfahren"), are fix-point iterations, that is, can be
1673 described by an iteration function ""x(t+1) = Phi( x(t) )"" which
1674 has the property:
1675 Phi(x) = x <==> A * x = b
1677 We can define "Phi(x)" as follows:
1679 Phi(x) := ( En - A ) * x + b
1681 where "En" is a matrix of the same size as "A" ("n" rows and
1683 columns) with one's on its main diagonal and zero's elsewhere.
1684 This function has the required property.
1686 Proof:
1688 A * x = b
1690 <==> -( A * x ) = -b
1692 <==> -( A * x ) + x = -b + x
1694 <==> -( A * x ) + x + b = x
1696 <==> x - ( A * x ) + b = x
1698 <==> ( En - A ) * x + b = x
1700 This last step is true because
1702 x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
1704 is the same as
1706 ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
1708 qed
1710 Note that actually solving the equation system ""A * x = b"" means
1712 to calculate
1713 a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] = b[i]
1715 <==> a[i,i] x[i] =
1717 b[i]
1718 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1719 + a[i,i] x[i]
1720 <==> x[i] =
1722 ( b[i]
1723 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1724 + a[i,i] x[i]
1725 ) / a[i,i]
1726 <==> x[i] =
1728 ( b[i] -
1729 ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
1730 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
1731 ) / a[i,i]
1732 There is one major restriction, though: a fix-point iteration is
1734 guaranteed to converge only if the first derivative of the
1735 iteration function has an absolute value less than one in an area
1736 around the point "x(*)" for which ""Phi( x(*) ) = x(*)"" is to be
1737 true, and if the start vector "x(0)" lies within that area!
1738 This is best verified graphically, which unfortunately is
1740 impossible to do in this textual documentation!
1741 See literature on Numerical Analysis for details!
1743 In our case, this restriction translates to the following three
1745 conditions:
1746 There must exist a norm so that the norm of the matrix of the
1748 iteration function, "( En - A )", has a value less than one, the
1749 matrix "A" may not have any zero value on its main diagonal and the
1750 initial vector "x(0)" must be "good enough", i.e., "close enough"
1751 to the solution "x(*)".
1752 (Remember school math: the first derivative of a straight line
1754 given by ""y = a * x + b"" is "a"!)
1755 The three methods expect a (quadratic!) matrix "$matrix" as their
1757 first argument, a start vector "$x0_vector", a vector "$b_vector"
1758 (which is the vector "b" in your equation system ""A * x = b""), in
1759 the case of the "Relaxation Method" ("RM"), a real number "$weight"
1760 best between zero and two, and finally an error limit (real number)
1761 "$epsilon".
1762 (Note that the weight "$weight" used by the "Relaxation Method"
1764 ("RM") is NOT checked to lie within any reasonable range!)
1765 The three methods first test the first two conditions of the three
1767 conditions listed above and return an empty list if these
1768 conditions are not fulfilled.
1769 Therefore, you should always test their return value using some
1771 code like:
1772 if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
1774 {
1775 # do something with the solution...
1776 }
1777 else
1778 {
1779 # do something with the fact that there is no solution...
1780 }
1781 Otherwise, they iterate until "abs( Phi(x) - x ) < epsilon".
1783 (Beware that theoretically, infinite loops might result if the
1785 starting vector is too far "off" the solution! In practice, this
1786 shouldn't be a problem. Anyway, you can always press <ctrl-C> if
1787 you think that the iteration takes too long!)
1788 The difference between the three methods is the following:
1790 In the "Global Step Method" ("GSM"), the new vector ""x(t+1)""
1792 (called "y" here) is calculated from the vector "x(t)" (called "x"
1793 here) according to the formula:
1794 y[i] =
1796 ( b[i]
1797 - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
1798 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
1799 ) / a[i,i]
1800 In the "Single Step Method" ("SSM"), the components of the vector
1802 ""x(t+1)"" which have already been calculated are used to calculate
1803 the remaining components, i.e.
1804 y[i] =
1806 ( b[i]
1807 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1808 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1809 ) / a[i,i]
1810 In the "Relaxation method" ("RM"), the components of the vector
1812 ""x(t+1)"" are calculated by "mixing" old and new value (like cold
1813 and hot water), and the weight "$weight" determines the "aperture"
1814 of both the "hot water tap" as well as of the "cold water tap",
1815 according to the formula:
1816 y[i] =
1818 ( b[i]
1819 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1820 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1821 ) / a[i,i]
1822 y[i] = weight * y[i] + (1 - weight) * x[i]
1823 Note that the weight "$weight" should be greater than zero and less
1825 than two (!).
1826 The three methods are supposed to be of different efficiency.
1828 Experiment!
1829 Remember that in most cases, it is probably advantageous to first
1831 "normalize()" your equation system prior to solving it!
1832
1834 SYNOPSIS
1835 • Unary operators:
1836 ""-"", ""~"", ""abs"", "test", ""!"", '""'
1838 • Binary operators:
1840 "".""
1842 Binary (arithmetic) operators:
1844 ""+"", ""-"", ""*"", ""**"", ""+="", ""-="", ""*="", ""/="",""**=""
1846 • Binary (relational) operators:
1848 ""=="", ""!="", ""<"", ""<="", "">"", "">=""
1850 ""eq"", ""ne"", ""lt"", ""le"", ""gt"", ""ge""
1852 Note that the latter (""eq"", ""ne"", ... ) are just synonyms of the
1854 former (""=="", ""!="", ... ), defined for convenience only.
1855 DESCRIPTION
1857 '.' Concatenation
1858 Returns the two matrices concatenated side by side.
1860 Example: $c = $a . $b;
1862 For example, if
1864 $a=[ 1 2 ] $b=[ 5 6 ]
1866 [ 3 4 ] [ 7 8 ]
1867 then
1868 $c=[ 1 2 5 6 ]
1870 [ 3 4 7 8 ]
1871 Note that only matrices with the same number of rows may be
1873 concatenated.
1874 '-' Unary minus
1876 Returns the negative of the given matrix, i.e., the matrix with
1878 all elements multiplied with the factor "-1".
1879 Example:
1881 $matrix = -$matrix;
1883 '~' Transposition
1885 Returns the transposed of the given matrix.
1887 Examples:
1889 $temp = ~$vector * $vector;
1891 $length = sqrt( $temp->element(1,1) );
1892 if (~$matrix == $matrix) { # matrix is symmetric ... }
1894 abs Norm
1896 Returns the "one"-Norm of the given matrix.
1898 Example:
1900 $error = abs( $A * $x - $b );
1902 test Boolean test
1904 Tests wether there is at least one non-zero element in the matrix.
1906 Example:
1908 if ($xn_vector) { # result of iteration is not zero ... }
1910 '!' Negated boolean test
1912 Tests wether the matrix contains only zero's.
1914 Examples:
1916 if (! $b_vector) { # heterogenous equation system ... }
1918 else { # homogenous equation system ... }
1919 unless ($x_vector) { # not the null-vector! }
1921 '""""'
1923 "Stringify" operator
1924 Converts the given matrix into a string.
1926 Uses scientific representation to keep precision loss to a minimum
1928 in case you want to read this string back in again later with
1929 "new_from_string()".
1930 By default a 13-digit mantissa and a 20-character field for each
1932 element is used so that lines will wrap nicely on an 80-column
1933 screen.
1934 Examples:
1936 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
1938 [ 1 0 ]
1939 [ 0 -1 ]
1940 MATRIX
1941 print "$matrix";
1942 [ 1.000000000000E+00 0.000000000000E+00 ]
1944 [ 0.000000000000E+00 -1.000000000000E+00 ]
1945 $string = "$matrix";
1947 $test = Math::MatrixReal->new_from_string($string);
1948 if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
1949 '+' Addition
1951 Returns the sum of the two given matrices.
1953 Examples:
1955 $matrix_S = $matrix_A + $matrix_B;
1957 $matrix_A += $matrix_B;
1959 '-' Subtraction
1961 Returns the difference of the two given matrices.
1963 Examples:
1965 $matrix_D = $matrix_A - $matrix_B;
1967 $matrix_A -= $matrix_B;
1969 Note that this is the same as:
1971 $matrix_S = $matrix_A + -$matrix_B;
1973 $matrix_A += -$matrix_B;
1975 (The latter are less efficient, though)
1977 '*' Multiplication
1979 Returns the matrix product of the two given matrices or the
1981 product of the given matrix and scalar factor.
1982 Examples:
1984 $matrix_P = $matrix_A * $matrix_B;
1986 $matrix_A *= $matrix_B;
1988 $vector_b = $matrix_A * $vector_x;
1990 $matrix_B = -1 * $matrix_A;
1992 $matrix_B = $matrix_A * -1;
1994 $matrix_A *= -1;
1996 '/' Division
1998 Currently a shortcut for doing $a * $b ** -1 is $a / $b, which
2000 works for square matrices. One can also use 1/$a .
2001 '**' Exponentiation
2003 Returns the matrix raised to an integer power. If 0 is passed, the
2005 identity matrix is returned. If a negative integer is passed, it
2006 computes the inverse (if it exists) and then raised the inverse to
2007 the absolute value of the integer. The matrix must be quadratic.
2008 Examples:
2010 $matrix2 = $matrix ** 2;
2012 $matrix **= 2;
2014 $inv2 = $matrix ** -2;
2016 $ident = $matrix ** 0;
2018 '==' Equality
2020 Tests two matrices for equality.
2022 Example:
2024 if ( $A * $x == $b ) { print "EUREKA!\n"; }
2026 Note that in most cases, due to numerical errors (due to the
2028 finite precision of computer arithmetics), it is a bad idea to
2029 compare two matrices or vectors this way.
2030 Better use the norm of the difference of the two matrices you want
2032 to compare and compare that norm with a small number, like this:
2033 if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
2035 '!=' Inequality
2037 Tests two matrices for inequality.
2039 Example:
2041 while ($x0_vector != $xn_vector) { # proceed with iteration ... }
2043 (Stops when the iteration becomes stationary)
2045 Note that (just like with the '==' operator), it is usually a bad
2047 idea to compare matrices or vectors this way. Compare the norm of
2048 the difference of the two matrices with a small number instead.
2049 '<' Less than
2051 Examples:
2053 if ( $matrix1 < $matrix2 ) { # ... }
2055 if ( $vector < $epsilon ) { # ... }
2057 if ( 1E-12 < $vector ) { # ... }
2059 if ( $A * $x - $b < 1E-12 ) { # ... }
2061 These are just shortcuts for saying:
2063 if ( abs($matrix1) < abs($matrix2) ) { # ... }
2065 if ( abs($vector) < abs($epsilon) ) { # ... }
2067 if ( abs(1E-12) < abs($vector) ) { # ... }
2069 if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
2071 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2073 scalars.
2074 '<=' Less than or equal
2076 As with the '<' operator, this is just a shortcut for the same
2078 expression with "abs()" around all arguments.
2079 Example:
2081 if ( $A * $x - $b <= 1E-12 ) { # ... }
2083 which in fact is the same as:
2085 if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
2087 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2089 scalars.
2090 '>' Greater than
2092 As with the '<' and '<=' operator, this
2094 if ( $xn - $x0 > 1E-12 ) { # ... }
2096 is just a shortcut for:
2098 if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
2100 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2102 scalars.
2103 '>=' Greater than or equal
2105 As with the '<', '<=' and '>' operator, the following
2107 if ( $LR >= $A ) { # ... }
2109 is simply a shortcut for:
2111 if ( abs($LR) >= abs($A) ) { # ... }
2113 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2115 scalars.
2116
2118 Math::VectorReal, Math::PARI, Math::MatrixBool, Math::Vec, DFA::Kleene,
2119 Math::Kleene, Set::IntegerRange, Set::IntegerFast .
2120
2122 This man page documents Math::MatrixReal version 2.13
2123 The latest code can be found at
2125 https://github.com/leto/math--matrixreal .
2126
2128 Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>,
2129 Jonathan "Duke" Leto <jonathan@leto.net>.
2130 Currently maintained by Jonathan "Duke" Leto, send all bugs/patches to
2132 Github Issues: https://github.com/leto/math--matrixreal/issues
2133
2135 Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
2136 Algebra and Linear Algebra at the university (RWTH Aachen, Germany),
2137 and to Prof. Esser and his assistant, Mr. Jarausch, for their
2138 fascinating lectures in Numerical Analysis!
2139
2141 Copyright (c) 1996-2016 by various authors including the original
2142 developer Steffen Beyer, Rodolphe Ortalo, the current maintainer
2143 Jonathan "Duke" Leto and all the wonderful people in the AUTHORS file.
2144 All rights reserved.
2145
2147 This package is free software; you can redistribute it and/or modify it
2148 under the same terms as Perl itself. Fuck yeah.
2149
2151 Hey! The above document had some coding errors, which are explained
2152 below:
2153 Around line 4032:
2155 '=item' outside of any '=over'
2156perl v5.34.0 2022-01-21 Math::MatrixReal(3)