* Add functions sub_matrix() and reduced_matrix().

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 9c28e7e7b0c6adc5cef8d8b762c43fd3b819b6ac..924ad5d7e375db0fb7c7a6ccd01748ce5bcdcf24 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -1974,6 +1974,38 @@ by @samp{c}) unit matrix. And finally, @code{symbolic_matrix} constructs a
 matrix filled with newly generated symbols made of the specified base name
 and the position of each element in the matrix.
 
 matrix filled with newly generated symbols made of the specified base name
 and the position of each element in the matrix.
 

+Matrices often arise by omitting elements of another matrix. For
+instance, the submatrix @code{S} of a matrix @code{M} takes a
+rectangular block from @code{M}. The reduced matrix @code{R} is defined
+by removing one row and one column from a matrix @code{M}. (The
+determinant of a reduced matrix is called a @emph{Minor} of @code{M} and
+can be used for computing the inverse using Cramer's rule.)
+
+@cindex @code{sub_matrix()}
+@cindex @code{reduced_matrix()}
+@example
+ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc);
+ex reduced_matrix(const matrix& m, unsigned r, unsigned c);
+@end example
+
+The function @code{sub_matrix()} takes a row offset @code{r} and a
+column offset @code{c} and takes a block of @code{nr} rows and @code{nc}
+columns. The function @code{reduced_matrix()} has two integer arguments
+that specify which row and column to remove:
+
+@example
+@{
+    matrix m(3,3);
+    m = 11, 12, 13,
+        21, 22, 23,
+        31, 32, 33;
+    cout << reduced_matrix(m, 1, 1) << endl;
+    // -> [[11,13],[31,33]]
+    cout << sub_matrix(m, 1, 2, 1, 2) << endl;
+    // -> [[22,23],[32,33]]
+@}
+@end example
+

 Matrix elements can be accessed and set using the parenthesis (function call)
 operator:
 
 Matrix elements can be accessed and set using the parenthesis (function call)
 operator:
 

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 623537febbb965f5ac1324e77f85e07d59a00768..710190376a1c28b59c787d902655aa24d97ccf73 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1598,4 +1598,52 @@ ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const
        return M;
 }
 
        return M;
 }
 

+ex reduced_matrix(const matrix& m, unsigned r, unsigned c)
+{
+       if (r+1>m.rows() || c+1>m.cols() || m.cols()<2 || m.rows()<2)
+               throw std::runtime_error("minor_matrix(): index out of bounds");
+
+       const unsigned rows = m.rows()-1;
+       const unsigned cols = m.cols()-1;
+       matrix &M = *new matrix(rows, cols);
+       M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+       unsigned ro = 0;
+       unsigned ro2 = 0;
+       while (ro2<rows) {
+               if (ro==r)
+                       ++ro;
+               unsigned co = 0;
+               unsigned co2 = 0;
+               while (co2<cols) {
+                       if (co==c)
+                               ++co;
+                       M(ro2,co2) = m(ro, co);
+                       ++co;
+                       ++co2;
+               }
+               ++ro;
+               ++ro2;
+       }
+
+       return M;
+}
+
+ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc)
+{
+       if (r+nr>m.rows() || c+nc>m.cols())
+               throw std::runtime_error("sub_matrix(): index out of bounds");
+
+       matrix &M = *new matrix(nr, nc);
+       M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+       for (unsigned ro=0; ro<nr; ++ro) {
+               for (unsigned co=0; co<nc; ++co) {
+                       M(ro,co) = m(ro+r,co+c);
+               }
+       }
+
+       return M;
+}
+

 } // namespace GiNaC
 } // namespace GiNaC

diff --git a/ginac/matrix.h b/ginac/matrix.h
index ca4b3c8fd4509423e63f6291937856e598749458..3941d3c5bd21630dcb9ca58493fbe9c990ef07fd 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -233,6 +233,13 @@ inline ex unit_matrix(unsigned x)
  *  The base name for LaTeX output is specified separately. */
 extern ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name);
 
  *  The base name for LaTeX output is specified separately. */
 extern ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name);
 

+/** Return the reduced matrix that is formed by deleting the rth row and cth
+ *  column of matrix m. The determinant of the result is the Minor r, c. */
+extern ex reduced_matrix(const matrix& m, unsigned r, unsigned c);
+
+/** Return the nr times nc submatrix starting at position r, c of matrix m. */
+extern ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc);
+

 /** Create an r times c matrix of newly generated symbols consisting of the
  *  given base name plus the numeric row/column position of each element. */
 inline ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name)
 /** Create an r times c matrix of newly generated symbols consisting of the
  *  given base name plus the numeric row/column position of each element. */
 inline ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name)