From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Thu, 20 Oct 2005 01:20:15 +0000 (+0000)
Subject: * Add functions sub_matrix() and reduced_matrix().
X-Git-Tag: release_1-4-0~143
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=a8c81ff424cab3ac522a71665b0eda55a8ca2f4d

* Add functions sub_matrix() and reduced_matrix().
---

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 9c28e7e7..924ad5d7 100644
--- 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.
 
+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:
 
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 623537fe..71019037 100644
--- 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;
 }
 
+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
diff --git a/ginac/matrix.h b/ginac/matrix.h
index ca4b3c8f..3941d3c5 100644
--- 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);
 
+/** 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)