1Math::MatrixReal(3)   User Contributed Perl Documentation  Math::MatrixReal(3)
2
3
4

NAME

6       Math::MatrixReal - Matrix of Reals
7
8       Implements the data type "matrix of real numbers" (and consequently
9       also "vector of real numbers").
10

SYNOPSIS

12       my $a = Math::MatrixReal->new_random(5, 5);
13
14       my $b = $a->new_random(10, 30, { symmetric=>1, bounded_by=>[-1,1] });
15
16       my $c = $b * $a ** 3;
17
18       my $d = $b->new_from_rows( [ [ 5, 3 ,4], [3, 4, 5], [ 2, 4, 1 ] ] );
19
20       print $a;
21
22       my $row        = ($a * $b)->row(3);
23
24       my $col        = (5*$c)->col(2);
25
26       my $transpose  = ~$c;
27
28       my $transpose  = $c->transpose;
29
30       my $inverse    = $a->inverse;
31
32       my $inverse    = 1/$a;
33
34       my $inverse    = $a ** -1;
35
36       my $determinant= $a->det;
37
38       •   $matrix->display_precision($integer)
39
40           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

FUNCTIONS

48   Constructor Methods And Such
49       •   use Math::MatrixReal;
50
51           Makes the methods and overloaded operators of this module available
52           to your program.
53
54       •   $new_matrix = new Math::MatrixReal($rows,$columns);
55
56           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
59           Note that this method is implicitly called by many of the other
60           methods in this module.
61
62       •   $new_matrix = $some_matrix->new($rows,$columns);
63
64           Another way of calling the matrix object constructor method.
65
66           Matrix $some_matrix is not changed by this in any way.
67
68       •   $new_matrix = $matrix->new_from_cols( [
69           $column_vector|$array_ref|$string, ... ] )
70
71           Creates a new matrix given a reference to an array of any of the
72           following:
73
74           •   column vectors ( n by 1 Math::MatrixReal matrices )
75
76           •   references to arrays
77
78           •   strings properly formatted to create a column with
79               Math::MatrixReal's new_from_string command
80
81           You may mix and match these as you wish.  However, all must be of
82           the same dimension--no padding happens automatically.  Example:
83
84               my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
85               print $matrix;
86
87           will print
88
89               [  1.000000000000E+00  3.000000000000E+00 ]
90               [  2.000000000000E+00  4.000000000000E+00 ]
91
92       •   new_from_rows( [ $row_vector|$array_ref|$string, ... ] )
93
94           Creates a new matrix given a reference to an array of any of the
95           following:
96
97           •   row vectors ( 1 by n Math::MatrixReal matrices )
98
99           •   references to arrays
100
101           •   strings properly formatted to create a row with
102               Math::MatrixReal's new_from_string command
103
104           You may mix and match these as you wish.  However, all must be of
105           the same dimension--no padding happens automatically. Example:
106
107                   my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
108                   print $matrix;
109
110           will print
111
112                   [  1.000000000000E+00  2.000000000000E+00 ]
113                   [  3.000000000000E+00  4.000000000000E+00 ]
114
115       •   $new_matrix = Math::MatrixReal->new_random($rows, $cols, %options
116           );
117
118           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
124            Example:
125
126               $matrix = Math::MatrixReal->new_random(4, { diagonal => 1, integer => 1 }  );
127               print $matrix;
128
129           will print a 4x4 random diagonal matrix with integer entries
130           between zero and ten, something like
131
132               [  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
137       •   $new_matrix = Math::MatrixReal->new_diag( $array_ref );
138
139           This method allows you to create a diagonal matrix by only
140           specifying the diagonal elements. Example:
141
142               $matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
143               print $matrix;
144
145           will print
146
147               [  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
152       •   $new_matrix = Math::MatrixReal->new_tridiag( $lower, $diag, $upper
153           );
154
155           This method allows you to create a tridiagonal matrix by only
156           specifying the lower diagonal, diagonal and upper diagonal,
157           respectively.
158
159               $matrix = Math::MatrixReal->new_tridiag( [ 6, 4, 2 ], [1,2,3,4], [1, 8, 9] );
160               print $matrix;
161
162           will print
163
164               [  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
169       •   $new_matrix = Math::MatrixReal->new_from_string($string);
170
171           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
174           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
179           Additional spaces or tabs can be added at will, but no comments.
180
181           Examples:
182
183             $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
187           By the way, this prints
188
189             [  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
193           But you can also do this in a much more comfortable way using the
194           shell-like "here-document" syntax:
195
196             $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
206           You can even use variables in the matrix:
207
208             $c1 =   2  /  3;
209             $c2 =  -2  /  5;
210             $c3 =  26  /  9;
211
212             $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
213
214                 [   3    2    0   ]
215                 [   0    3    2   ]
216                 [  $c1  $c2  $c3  ]
217
218             MATRIX
219
220           (Remember that you may use spaces and tabs to format the matrix to
221           your taste)
222
223           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
228           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
231           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
235             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
250           or as follows:
251
252             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
260           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
264           A better way is shown in this piece of code:
265
266             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
280           Possible error messages of the "new_from_string()" method are:
281
282             Math::MatrixReal::new_from_string(): syntax error in input string
283             Math::MatrixReal::new_from_string(): empty input string
284
285           If the input string has rows with varying numbers of columns, the
286           following warning will be printed to STDERR:
287
288             Math::MatrixReal::new_from_string(): missing elements will be set to zero!
289
290           If everything is okay, the method returns an object reference to
291           the (newly allocated) matrix containing the elements you specified.
292
293       •   $new_matrix = $some_matrix->shadow();
294
295           Returns an object reference to a NEW but EMPTY matrix (filled with
296           zero's) of the SAME SIZE as matrix "$some_matrix".
297
298           Matrix "$some_matrix" is not changed by this in any way.
299
300       •   $matrix1->copy($matrix2);
301
302           Copies the contents of matrix "$matrix2" to an ALREADY EXISTING
303           matrix "$matrix1" (which must have the same size as matrix
304           "$matrix2"!).
305
306           Matrix "$matrix2" is not changed by this in any way.
307
308       •   $twin_matrix = $some_matrix->clone();
309
310           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
316           Matrix "$some_matrix" is not changed by this in any way.
317
318       •   $matrix = Math::MatrixReal->reshape($rows, $cols, $array_ref);
319
320           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
325                $matrix = Math::MatrixReal->reshape(4, 3, [1..12]);
326
327           Creates the following matrix:
328
329               [ 1    5    9 ]
330               [ 2    6   10 ]
331               [ 3    7   11 ]
332               [ 4    8   12 ]
333
334   Matrix Row, Column and Element operations
335       •   $value = $matrix->element($row,$column);
336
337           Returns the value of a specific element of the matrix "$matrix",
338           located in row "$row" and column "$column".
339
340           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
344               $elem = $matrix->element(1, 1); # first element of the matrix.
345
346       •   $matrix->assign($row,$column,$value);
347
348           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
352       •   $row_vector = $matrix->row($row);
353
354           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
359           Matrix "$matrix" is not changed by this in any way.
360
361       •   $column_vector = $matrix->column($column);
362
363           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
368           Matrix "$matrix" is not changed by this in any way.
369
370       •   @all_elements = $matrix->as_list;
371
372           Get the contents of a Math::MatrixReal object as a Perl list.
373
374           Example:
375
376              my $matrix = Math::MatrixReal->new_from_rows([ [1, 2], [3, 4] ]);
377              my @list = $matrix->as_list; # 1, 2, 3, 4
378
379           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
383            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
406           See OpenGL, OpenGL::Shader, OpenGL::Array, rotation matrix
407           <https://en.wikipedia.org/wiki/Rotation_matrix>.
408
409       •   $new_matrix = $matrix->each( \&function );
410
411           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
417           Example:
418
419               # add 1 to every element in the matrix
420               $matrix = $matrix->each ( sub { (shift) + 1 } );
421
422           Example:
423
424               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
429           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
435       •   $new_matrix = $matrix->each_diag( \&function );
436
437           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
444       •   $matrix->swap_col( $col1, $col2 );
445
446           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
451       •   $matrix->swap_row( $row1, $row2 );
452
453           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
457       •   $matrix->assign_row( $row_number , $new_row_vector );
458
459           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
464       •   $matrix->maximum();  and  $matrix->minimum();
465
466           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
471           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
476           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
481           Consider the matrix and vector below for the following examples:
482
483                      [ 1 9 4 ]
484                 $A = [ 3 5 2 ]       $B = [ 8 7 9 5 3 ]
485                      [ 8 7 6 ]
486
487           When used in scalar context:
488
489               $max = $A->maximum();    # $max = [ 8, 9, 6 ]
490               $min = $B->minimum();    # $min = 3
491
492           When used in list context:
493
494               ($min, $pos) = $A->minimum(); # $min = [ 1 5 2 ]
495                                             # $pos = [ 1 2 2 ]
496               ($max, $pos) = $B->maximum(); # $max = 9
497                                             # $pos = 3
498
499   Matrix Operations
500       •   "$det = $matrix->det();"
501
502           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
510       •   "$inverse = $matrix->inverse();"
511
512           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
517       •   "($rows,$columns) = $matrix->dim();"
518
519           Returns a list of two items, representing the number of rows and
520           columns the given matrix "$matrix" contains.
521
522       •   "$norm_one = $matrix->norm_one();"
523
524           Returns the "one"-norm of the given matrix "$matrix".
525
526           The "one"-norm is defined as follows:
527
528           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
532           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
536           Therefore, you should only compare values that have been calculated
537           using the same norm!
538
539           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
544       •   "$norm_max = $matrix->norm_max();"
545
546           Returns the "maximum"-norm of the given matrix $matrix.
547
548           The "maximum"-norm is defined as follows:
549
550           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
554           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
558           Therefore, you should only compare values that have been calculated
559           using the same norm!
560
561           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
566       •   "$norm_sum = $matrix->norm_sum();"
567
568           This is a very simple norm which is defined as the sum of the
569           absolute values of every element.
570
571       •   $p_norm = $matrix->norm_p($n);>
572
573           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
584           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
591               print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";
592
593               $i1 = $a->new_from_rows([[1,0]]);
594               $i2 = $a->new_from_rows([[0,1]]);
595
596               # this should be sqrt(2) since it is the same as the
597               # hypotenuse of a 1 by 1 right triangle
598
599               $dist  = ($i1-$i2)->norm_p(2);
600               print "Distance is $dist, which should be " . sqrt(2) . "\n";
601
602           Output:
603
604               (1,2,3,Inf) norm:
605               6
606               3.74165738677394139
607               3.30192724889462668
608               3
609
610               Distance is 1.41421356237309505, which should be 1.41421356237309505
611
612       •   $frob_norm = "$matrix->norm_frobenius();"
613
614           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
618       •   "$matrix->spectral_radius();"
619
620           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
626       •   "$matrix1->transpose($matrix2);"
627
628           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
632           This operation can also be carried out "in-place", i.e., input and
633           output matrix may be identical.
634
635           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
639           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
643           Note that (especially for vectors) it makes a big difference if you
644           have a row vector, like this:
645
646             [ -1 0 1 ]
647
648           or a column vector, like this:
649
650             [ -1 ]
651             [  0 ]
652             [  1 ]
653
654           the one vector being the transposed of the other!
655
656           This is especially true for the matrix product of two vectors:
657
658                          [ -1 ]
659             [ -1 0 1 ] * [  0 ]  =  [ 2 ] ,  whereas
660                          [  1 ]
661
662                                        *     [ -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
668           So be careful about what you really mean!
669
670           Hint: throughout this module, whenever a vector is explicitly
671           required for input, a COLUMN vector is expected!
672
673       •   "$trace = $matrix->trace();"
674
675           This returns the trace of the matrix, which is defined as the sum
676           of the diagonal elements. The matrix must be quadratic.
677
678       •   "$minor = $matrix->minor($row,$col);"
679
680           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
687       •   "$cofactor = $matrix->cofactor();"
688
689           The cofactor matrix is constructed as follows:
690
691           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
697           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
701           The following two inverses should be exactly the same:
702
703               my $inverse1 = $matrix->inverse;
704               my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
705
706           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
711               # $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
718
719               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
723       •   "$adjoint = $matrix->adjoint();"
724
725           The adjoint is just the transpose of the cofactor matrix. This
726           method is just an alias for " ~($matrix->cofactor)".
727
728       •   "$part_of_matrix = $matrix->submatrix(x1,y1,x2,Y2);"
729
730           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
734           - firstly: Coordinate of top left corner to select (x1,y1)
735           - secondly: Coordinate of bottom right corner to select (x2,y2)
736
737           Example:
738
739               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
749               my $submatrix = $matrix->submatrix(5,5,7,7);
750               $submatrix->display_precision(0);
751               print $submatrix;
752
753           Output:
754
755               [  1  0  1  ]
756               [  0  1  0  ]
757               [  1  0  1  ]
758
759   Arithmetic Operations
760       •   "$matrix1->add($matrix2,$matrix3);"
761
762           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
767           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
770       •   "$matrix1->subtract($matrix2,$matrix3);"
771
772           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
777           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
780           Note that this operation is the same as
781           "$matrix1->add($matrix2,-$matrix3);", although the latter is a
782           little less efficient.
783
784       •   "$matrix1->multiply_scalar($matrix2,$scalar);"
785
786           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
792           This operation can also be carried out "in-place", i.e., input and
793           output matrix may be identical.
794
795       •   "$product_matrix = $matrix1->multiply($matrix2);"
796
797           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
801           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
805                                     [ 2 2 ]
806                                     [ 2 2 ]
807                                     [ 2 2 ]
808
809                         [ 1 1 1 ]   [ * * ]
810                         [ 1 1 1 ]   [ * * ]
811                         [ 1 1 1 ]   [ * * ]
812                         [ 1 1 1 ]   [ * * ]
813
814           I.e., the number of columns of matrix "$matrix1" has to be the same
815           as the number of rows of matrix "$matrix2".
816
817           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
822       •   "$matrix1->negate($matrix2);"
823
824           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
829           This operation can also be carried out "in-place", i.e., input and
830           output matrix may be identical.
831
832       •   "$matrix_to_power = $matrix1->exponent($integer);"
833
834           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
840   Boolean Matrix Operations
841       •   "$matrix->is_quadratic();"
842
843           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
847       •   "$matrix->is_square();"
848
849           This is an alias for "is_quadratic()".
850
851       •   "$matrix->is_symmetric();"
852
853           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
859           Notes: A symmetric matrix always has real eigenvalues/eigenvectors.
860           A matrix plus its transpose is always symmetric.
861
862       •   "$matrix->is_skew_symmetric();"
863
864           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
870       •   "$matrix->is_diagonal();"
871
872           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
876       •   "$matrix->is_tridiagonal();"
877
878           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
883       •   "$matrix->is_upper_triangular();"
884
885           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
891       •   "$matrix->is_lower_triangular();"
892
893           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
899       •   "$matrix->is_orthogonal();"
900
901           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
908       •   "$matrix->is_binary();"
909
910           Returns a boolean value indicating if the given matrix is binary.
911           A matrix is binary if it contains only zeroes or ones.
912
913       •   "$matrix->is_gramian();"
914
915           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
921       •   "$matrix->is_LR();"
922
923           Returns a boolean value indicating if the matrix is an LR
924           decomposition matrix.
925
926       •   "$matrix->is_positive();"
927
928           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
932       •   "$matrix->is_negative();"
933
934           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
938       •   "$matrix->is_periodic($k);"
939
940           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
944       •   "$matrix->is_idempotent();"
945
946           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
950       •   "$matrix->is_row_vector();"
951
952           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
956       •   "$matrix->is_col_vector();"
957
958           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
962   Eigensystems
963       • "($l, $V) = $matrix->sym_diagonalize();"
964
965         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
972         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
982       • "($T, $Q) = $matrix->householder();"
983
984         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
990       • "($l, $V) = $T->tri_diagonalize([$Q]);"
991
992         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
996         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
1005       • "$l = $matrix->sym_eigenvalues();"
1006
1007         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
1016         This routine wraps the householder_tridiagonal() and
1017         tri_eigenvalues() methods described below when the intermediate
1018         tridiagonal matrix is not needed.
1019
1020       • "$T = $matrix->householder_tridiagonal();"
1021
1022         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
1029       • $l = $T->tri_eigenvalues();
1030
1031         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
1036   Miscellaneous
1037       •   $matrix->zero();
1038
1039           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
1044           (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
1050       •   $matrix->one();
1051
1052           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
1058           (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
1063       •   "$latex_string = $matrix->as_latex( align=> "c", format => "%s",
1064           name => "" );"
1065
1066           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
1075           Example:
1076
1077               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
1080               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
1089       •   "$yacas_string = $matrix->as_yacas( format => "%s", name => "",
1090           semi => 0 );"
1091
1092           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
1101           Example:
1102
1103               $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
1106           Output:
1107
1108               A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
1109
1110       •   "$matlab_string = $matrix->as_matlab( format => "%s", name => "",
1111           semi => 0 );"
1112
1113           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
1122           Example:
1123
1124                   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
1127           Output:
1128                   A = [ 1.234 5.678 9.101;
1129                    1.000 2.000 3.000];
1130
1131       •   "$scilab_string = $matrix->as_scilab( format => "%s", name => "",
1132           semi => 0 );"
1133
1134           This function is just an alias for "as_matlab()", since both Scilab
1135           and Matlab have the same matrix format.
1136
1137       •   "$minimum = Math::MatrixReal::min($number1,$number2);" "$minimum =
1138           Math::MatrixReal::min($matrix);" "<$minimum = $matrix-"min;>>
1139
1140           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
1145       •   "$maximum = Math::MatrixReal::max($number1,$number2);" "$maximum =
1146           Math::MatrixReal::max($number1,$number2);" "$maximum =
1147           Math::MatrixReal::max($matrix);" "<$maximum = $matrix-"max;>>
1148
1149           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
1154       •   "$minimal_cost_matrix = $cost_matrix->kleene();"
1155
1156           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
1160           See Math::Kleene(3) for more details about this algorithm!
1161
1162           The method returns an object reference to the new matrix.
1163
1164           Matrix "$cost_matrix" is not changed by this method in any way.
1165
1166       •   "($norm_matrix,$norm_vector) = $matrix->normalize($vector);"
1167
1168           This method is used to improve the numerical stability when solving
1169           linear equation systems.
1170
1171           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
1176           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
1182           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
1186           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
1192           The output matrix and vector are clones of the input matrix and
1193           vector to which the operation explained above has been applied.
1194
1195           The input matrix and vector are not changed by this in any way.
1196
1197           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
1201           Consider the following little program:
1202
1203             #!perl -w
1204
1205             use Math::MatrixReal qw(new_from_string);
1206
1207             $A = Math::MatrixReal->new_from_string(<<"MATRIX");
1208             [  1   2   3  ]
1209             [  5   7  11  ]
1210             [ 23  19  13  ]
1211             MATRIX
1212
1213             $b = Math::MatrixReal->new_from_string(<<"MATRIX");
1214             [   0   ]
1215             [   1   ]
1216             [  29   ]
1217             MATRIX
1218
1219             $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
1227             ($A_,$b_) = $A->normalize($b);
1228
1229             $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
1237           This will print:
1238
1239             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
1256           You can see that in the second example (where "normalize()" has
1257           been used), the result is "better", i.e., more accurate!
1258
1259       •   "$LR_matrix = $matrix->decompose_LR();"
1260
1261           This method is needed to solve linear equation systems.
1262
1263           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
1268           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
1272           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
1278           "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
1282           "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
1286                     [ 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
1292           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
1297           Trick:
1298
1299           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
1303                            [ 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
1309           Beware, though, that "LR" and ""L * R"" are not the same!!!
1310
1311           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
1315           (You need to apply "normalize()" to each pair (Ai,bi) BEFORE
1316           decomposing the matrix "Ai'"!)
1317
1318           Now what does all this help us in solving linear equation systems?
1319
1320           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
1325           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
1330           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
1334             y[2] = b[2] - L[2,1] * y[1]
1335
1336           (and we know ""y[1]"" by now!), that
1337
1338             y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
1339
1340           and so on.
1341
1342           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
1346             x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
1347
1348           and
1349
1350             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
1353           and so on.
1354
1355           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
1360           This method requires a quadratic matrix as its input matrix.
1361
1362           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
1366           The method returns an object reference to a new matrix containing
1367           the matrices "L" and "R".
1368
1369           The input matrix is not changed by this method in any way.
1370
1371           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
1375           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
1381       •   "($dimension,$x_vector,$base_matrix) =
1382           $LR_matrix""->""solve_LR($b_vector);"
1383
1384           Use this method to actually solve an equation system.
1385
1386           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
1390           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
1394           The method returns a list of three items if a solution exists or an
1395           empty list otherwise (!).
1396
1397           Therefore, you should always use this method like this:
1398
1399             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
1408           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
1417           Only the first "$dimension" columns of this base matrix actually
1418           contain entries, the remaining columns are all zero.
1419
1420           Now what is all this stuff with that "base" good for?
1421
1422           The output vector "x" is ALWAYS a solution of your equation system
1423           "A * x = b".
1424
1425           But also any vector "$vector"
1426
1427             $vector = $x_vector->clone();
1428
1429             $machine_infinity = 1E+99; # or something like that
1430
1431             for ( $i = 1; $i <= $dimension; $i++ )
1432             {
1433                 $vector += rand($machine_infinity) * $base_matrix->column($i);
1434             }
1435
1436           is a solution to your problem "A * x = b", i.e., if "$A_matrix"
1437           contains your matrix "A", then
1438
1439             print abs( $A_matrix * $vector - $b_vector ), "\n";
1440
1441           should print a number around 1E-16 or so!
1442
1443           By the way, note that you can actually calculate those vectors
1444           "$vector" a little more efficient as follows:
1445
1446             $rand_vector = $x_vector->shadow();
1447
1448             $machine_infinity = 1E+99; # or something like that
1449
1450             for ( $i = 1; $i <= $dimension; $i++ )
1451             {
1452                 $rand_vector->assign($i,1, rand($machine_infinity) );
1453             }
1454
1455             $vector = $x_vector + ( $base_matrix * $rand_vector );
1456
1457           Note that the input matrix and vector are not changed by this
1458           method in any way.
1459
1460       •   "$inverse_matrix = $LR_matrix->invert_LR();"
1461
1462           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
1466           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
1471           Therefore, you should always use this method in the following way:
1472
1473             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
1482           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
1487           The input matrix is not changed by this method in any way.
1488
1489       •   "$condition = $matrix->condition($inverse_matrix);"
1490
1491           In fact this method is just a shortcut for
1492
1493             abs($matrix) * abs($inverse_matrix)
1494
1495           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
1499           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
1503           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
1507           Numerical stability means for example that if
1508
1509             abs( $vec_correct - $vec_with_error ) < $epsilon
1510
1511           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
1514             abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
1515
1516           also holds.
1517
1518       •   "$determinant = $LR_matrix->det_LR();"
1519
1520           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
1524           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
1529           (The sign is taken care of "magically" by this module)
1530
1531       •   "$order = $LR_matrix->order_LR();"
1532
1533           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
1537           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
1542           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
1546       •   "$rank = $LR_matrix->rank_LR();"
1547
1548           This is an alias for the "order_LR()" function. The "order" is
1549           usually called the "rank" in the United States.
1550
1551       •   "$scalar_product = $vector1->scalar_product($vector2);"
1552
1553           Returns the scalar product of vector "$vector1" and vector
1554           "$vector2".
1555
1556           Both vectors must be column vectors (i.e., a matrix having several
1557           rows but only one column).
1558
1559           This is a (more efficient!) shortcut for
1560
1561             $temp           = ~$vector1 * $vector2;
1562             $scalar_product =  $temp->element(1,1);
1563
1564           or the sum "i=1..n" of the products "vector1[i] * vector2[i]".
1565
1566           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
1571       •   "$vector_product = $vector1->vector_product($vector2);"
1572
1573           Returns the vector product of vector "$vector1" and vector
1574           "$vector2".
1575
1576           Both vectors must be column vectors (i.e., a matrix having several
1577           rows but only one column).
1578
1579           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
1583           In 3 dimensions, the vector product of two vectors "x" and "y" is
1584           defined as
1585
1586                         |  x[1]  y[1]  e[1]  |
1587             determinant |  x[2]  y[2]  e[2]  |
1588                         |  x[3]  y[3]  e[3]  |
1589
1590           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
1596           This determinant evaluates to the rather simple formula
1597
1598             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
1602           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
1608       •   "$length = $vector->length();"
1609
1610           This is actually a shortcut for
1611
1612             $length = sqrt( $vector->scalar_product($vector) );
1613
1614           and returns the length of a given column or row vector "$vector".
1615
1616           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
1619           The general definition for norms of vectors is the following:
1620
1621             sub vector_norm
1622             {
1623                 croak "Usage: \$norm = \$vector->vector_norm(\$n);"
1624                   if (@_ != 2);
1625
1626                 my($vector,$n) = @_;
1627                 my($rows,$cols) = ($vector->[1],$vector->[2]);
1628                 my($k,$comp,$sum);
1629
1630                 croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
1631                   unless ($cols == 1);
1632
1633                 croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
1634                   unless ($n > 0);
1635
1636                 croak "Math::MatrixReal::vector_norm(): norm index must be integer"
1637                   unless ($n == int($n));
1638
1639                 $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
1648           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
1653       •   "$xn_vector = $matrix->""solve_GSM($x0_vector,$b_vector,$epsilon);"
1654
1655       •   "$xn_vector = $matrix->""solve_SSM($x0_vector,$b_vector,$epsilon);"
1656
1657       •   "$xn_vector =
1658           $matrix->""solve_RM($x0_vector,$b_vector,$weight,$epsilon);"
1659
1660           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
1664           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
1669           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
1676             Phi(x)  =  x    <==>    A * x  =  b
1677
1678           We can define "Phi(x)" as follows:
1679
1680             Phi(x)  :=  ( En - A ) * x  +  b
1681
1682           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
1685           This function has the required property.
1686
1687           Proof:
1688
1689                      A * x        =   b
1690
1691             <==>  -( A * x )      =  -b
1692
1693             <==>  -( A * x ) + x  =  -b + x
1694
1695             <==>  -( A * x ) + x + b  =  x
1696
1697             <==>  x - ( A * x ) + b  =  x
1698
1699             <==>  ( En - A ) * x + b  =  x
1700
1701           This last step is true because
1702
1703             x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
1704
1705           is the same as
1706
1707             ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
1708
1709           qed
1710
1711           Note that actually solving the equation system ""A * x = b"" means
1712           to calculate
1713
1714                   a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n]  =  b[i]
1715
1716             <==>  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
1721             <==>  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
1727             <==>  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
1733           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
1739           This is best verified graphically, which unfortunately is
1740           impossible to do in this textual documentation!
1741
1742           See literature on Numerical Analysis for details!
1743
1744           In our case, this restriction translates to the following three
1745           conditions:
1746
1747           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
1753           (Remember school math: the first derivative of a straight line
1754           given by ""y = a * x + b"" is "a"!)
1755
1756           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
1763           (Note that the weight "$weight" used by the "Relaxation Method"
1764           ("RM") is NOT checked to lie within any reasonable range!)
1765
1766           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
1770           Therefore, you should always test their return value using some
1771           code like:
1772
1773             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
1782           Otherwise, they iterate until "abs( Phi(x) - x ) < epsilon".
1783
1784           (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
1789           The difference between the three methods is the following:
1790
1791           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
1795             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
1801           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
1805             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
1811           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
1817             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
1824           Note that the weight "$weight" should be greater than zero and less
1825           than two (!).
1826
1827           The three methods are supposed to be of different efficiency.
1828           Experiment!
1829
1830           Remember that in most cases, it is probably advantageous to first
1831           "normalize()" your equation system prior to solving it!
1832

OVERLOADED OPERATORS

1834   SYNOPSIS
1835       • Unary operators:
1836
1837         ""-"", ""~"", ""abs"", "test", ""!"", '""'
1838
1839       • Binary operators:
1840
1841         "".""
1842
1843         Binary (arithmetic) operators:
1844
1845         ""+"", ""-"", ""*"", ""**"", ""+="", ""-="", ""*="", ""/="",""**=""
1846
1847       • Binary (relational) operators:
1848
1849         ""=="", ""!="", ""<"", ""<="", "">"", "">=""
1850
1851         ""eq"", ""ne"", ""lt"", ""le"", ""gt"", ""ge""
1852
1853         Note that the latter (""eq"", ""ne"", ... ) are just synonyms of the
1854         former (""=="", ""!="", ... ), defined for convenience only.
1855
1856   DESCRIPTION
1857       '.'  Concatenation
1858
1859            Returns the two matrices concatenated side by side.
1860
1861            Example:      $c = $a . $b;
1862
1863            For example, if
1864
1865                    $a=[ 1 2 ]   $b=[ 5 6 ]
1866                       [ 3 4 ]      [ 7 8 ]
1867            then
1868
1869                    $c=[ 1 2 5 6 ]
1870                       [ 3 4 7 8 ]
1871
1872            Note that only matrices with the same number of rows may be
1873            concatenated.
1874
1875       '-'  Unary minus
1876
1877            Returns the negative of the given matrix, i.e., the matrix with
1878            all elements multiplied with the factor "-1".
1879
1880            Example:
1881
1882                $matrix = -$matrix;
1883
1884       '~'  Transposition
1885
1886            Returns the transposed of the given matrix.
1887
1888            Examples:
1889
1890                $temp = ~$vector * $vector;
1891                $length = sqrt( $temp->element(1,1) );
1892
1893                if (~$matrix == $matrix) { # matrix is symmetric ... }
1894
1895       abs  Norm
1896
1897            Returns the "one"-Norm of the given matrix.
1898
1899            Example:
1900
1901                $error = abs( $A * $x - $b );
1902
1903       test Boolean test
1904
1905            Tests wether there is at least one non-zero element in the matrix.
1906
1907            Example:
1908
1909                if ($xn_vector) { # result of iteration is not zero ... }
1910
1911       '!'  Negated boolean test
1912
1913            Tests wether the matrix contains only zero's.
1914
1915            Examples:
1916
1917                if (! $b_vector) { # heterogenous equation system ... }
1918                else             { # homogenous equation system ... }
1919
1920                unless ($x_vector) { # not the null-vector! }
1921
1922       '""""'
1923            "Stringify" operator
1924
1925            Converts the given matrix into a string.
1926
1927            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
1931            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
1935            Examples:
1936
1937                $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
1938                [ 1  0 ]
1939                [ 0 -1 ]
1940                MATRIX
1941                print "$matrix";
1942
1943                [  1.000000000000E+00  0.000000000000E+00 ]
1944                [  0.000000000000E+00 -1.000000000000E+00 ]
1945
1946                $string = "$matrix";
1947                $test = Math::MatrixReal->new_from_string($string);
1948                if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
1949
1950       '+'  Addition
1951
1952            Returns the sum of the two given matrices.
1953
1954            Examples:
1955
1956                $matrix_S = $matrix_A + $matrix_B;
1957
1958                $matrix_A += $matrix_B;
1959
1960       '-'  Subtraction
1961
1962            Returns the difference of the two given matrices.
1963
1964            Examples:
1965
1966                $matrix_D = $matrix_A - $matrix_B;
1967
1968                $matrix_A -= $matrix_B;
1969
1970            Note that this is the same as:
1971
1972                $matrix_S = $matrix_A + -$matrix_B;
1973
1974                $matrix_A += -$matrix_B;
1975
1976            (The latter are less efficient, though)
1977
1978       '*'  Multiplication
1979
1980            Returns the matrix product of the two given matrices or the
1981            product of the given matrix and scalar factor.
1982
1983            Examples:
1984
1985                $matrix_P = $matrix_A * $matrix_B;
1986
1987                $matrix_A *= $matrix_B;
1988
1989                $vector_b = $matrix_A * $vector_x;
1990
1991                $matrix_B = -1 * $matrix_A;
1992
1993                $matrix_B = $matrix_A * -1;
1994
1995                $matrix_A *= -1;
1996
1997       '/'  Division
1998
1999            Currently a shortcut for doing $a * $b ** -1 is $a / $b, which
2000            works for square matrices. One can also use 1/$a .
2001
2002       '**' Exponentiation
2003
2004            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
2009            Examples:
2010
2011                $matrix2 = $matrix ** 2;
2012
2013                $matrix **= 2;
2014
2015                $inv2 = $matrix ** -2;
2016
2017                $ident = $matrix ** 0;
2018
2019       '==' Equality
2020
2021            Tests two matrices for equality.
2022
2023            Example:
2024
2025                if ( $A * $x == $b ) { print "EUREKA!\n"; }
2026
2027            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
2031            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
2034                if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
2035
2036       '!=' Inequality
2037
2038            Tests two matrices for inequality.
2039
2040            Example:
2041
2042                while ($x0_vector != $xn_vector) { # proceed with iteration ... }
2043
2044            (Stops when the iteration becomes stationary)
2045
2046            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
2050       '<'  Less than
2051
2052            Examples:
2053
2054                if ( $matrix1 < $matrix2 ) { # ... }
2055
2056                if ( $vector < $epsilon ) { # ... }
2057
2058                if ( 1E-12 < $vector ) { # ... }
2059
2060                if ( $A * $x - $b < 1E-12 ) { # ... }
2061
2062            These are just shortcuts for saying:
2063
2064                if ( abs($matrix1) < abs($matrix2) ) { # ... }
2065
2066                if ( abs($vector) < abs($epsilon) ) { # ... }
2067
2068                if ( abs(1E-12) < abs($vector) ) { # ... }
2069
2070                if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
2071
2072            Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2073            scalars.
2074
2075       '<=' Less than or equal
2076
2077            As with the '<' operator, this is just a shortcut for the same
2078            expression with "abs()" around all arguments.
2079
2080            Example:
2081
2082                if ( $A * $x - $b <= 1E-12 ) { # ... }
2083
2084            which in fact is the same as:
2085
2086                if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
2087
2088            Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2089            scalars.
2090
2091       '>'  Greater than
2092
2093            As with the '<' and '<=' operator, this
2094
2095                if ( $xn - $x0 > 1E-12 ) { # ... }
2096
2097            is just a shortcut for:
2098
2099                if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
2100
2101            Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2102            scalars.
2103
2104       '>=' Greater than or equal
2105
2106            As with the '<', '<=' and '>' operator, the following
2107
2108                if ( $LR >= $A ) { # ... }
2109
2110            is simply a shortcut for:
2111
2112                if ( abs($LR) >= abs($A) ) { # ... }
2113
2114            Uses the "one"-norm for matrices and Perl's built-in "abs()" for
2115            scalars.
2116

SEE ALSO

2118       Math::VectorReal, Math::PARI, Math::MatrixBool, Math::Vec, DFA::Kleene,
2119       Math::Kleene, Set::IntegerRange, Set::IntegerFast .
2120

VERSION

2122       This man page documents Math::MatrixReal version 2.13
2123
2124       The latest code can be found at
2125       https://github.com/leto/math--matrixreal .
2126

AUTHORS

2128       Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>,
2129       Jonathan "Duke" Leto <jonathan@leto.net>.
2130
2131       Currently maintained by Jonathan "Duke" Leto, send all bugs/patches to
2132       Github Issues: https://github.com/leto/math--matrixreal/issues
2133

CREDITS

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

LICENSE AGREEMENT

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

POD ERRORS

2151       Hey! The above document had some coding errors, which are explained
2152       below:
2153
2154       Around line 4032:
2155           '=item' outside of any '=over'
2156
2157
2158
2159perl v5.32.1                      2021-01-27               Math::MatrixReal(3)
Impressum