X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=bf7f8347c9a1f7b31645595d408ce602cf77081d;hp=50f4631546a5795216287f27c9b9f272731fefa6;hb=d0ff428fb5b7bb565a0aea69e05e5705d840c16d;hpb=c6ead0d1f88d99b0d8cc99b2da8ef9e9c41a073f

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 50f46315..bf7f8347 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -3,7 +3,7 @@
  *  Implementation of symbolic matrices */
 
 /*
- *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2018 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -227,28 +227,6 @@ ex & matrix::let_op(size_t i)
 	return m[i];
 }
 
-/** Evaluate matrix entry by entry. */
-ex matrix::eval(int level) const
-{
-	// check if we have to do anything at all
-	if ((level==1)&&(flags & status_flags::evaluated))
-		return *this;
-	
-	// emergency break
-	if (level == -max_recursion_level)
-		throw (std::runtime_error("matrix::eval(): recursion limit exceeded"));
-	
-	// eval() entry by entry
-	exvector m2(row*col);
-	--level;
-	for (unsigned r=0; r<row; ++r)
-		for (unsigned c=0; c<col; ++c)
-			m2[r*col+c] = m[r*col+c].eval(level);
-	
-	return (new matrix(row, col, std::move(m2)))->setflag(status_flags::dynallocated |
-	                                                      status_flags::evaluated);
-}
-
 ex matrix::subs(const exmap & mp, unsigned options) const
 {
 	exvector m2(row * col);
@@ -961,10 +939,11 @@ ex matrix::charpoly(const ex & lambda) const
 
 /** Inverse of this matrix.
  *
+ *  @param algo selects the algorithm (one of solve_algo)
  *  @return    the inverted matrix
  *  @exception logic_error (matrix not square)
  *  @exception runtime_error (singular matrix) */
-matrix matrix::inverse() const
+matrix matrix::inverse(unsigned algo) const
 {
 	if (row != col)
 		throw (std::logic_error("matrix::inverse(): matrix not square"));
@@ -987,7 +966,7 @@ matrix matrix::inverse() const
 	
 	matrix sol(row,col);
 	try {
-		sol = this->solve(vars,identity);
+		sol = this->solve(vars, identity, algo);
 	} catch (const std::runtime_error & e) {
 	    if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
 			throw (std::runtime_error("matrix::inverse(): singular matrix"));
@@ -1075,11 +1054,11 @@ matrix matrix::solve(const matrix & vars,
 		unsigned last_assigned_sol = n+1;
 		for (int r=m-1; r>=0; --r) {
 			unsigned fnz = 1;    // first non-zero in row
-			while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].is_zero()))
+			while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].normal().is_zero()))
 				++fnz;
 			if (fnz>n) {
 				// row consists only of zeros, corresponding rhs must be 0, too
-				if (!aug.m[r*(n+p)+n+co].is_zero()) {
+				if (!aug.m[r*(n+p)+n+co].normal().is_zero()) {
 					throw (std::runtime_error("matrix::solve(): inconsistent linear system"));
 				}
 			} else {
@@ -1572,8 +1551,7 @@ ex lst_to_matrix(const lst & l)
 	}
 
 	// Allocate and fill matrix
-	matrix &M = *new matrix(rows, cols);
-	M.setflag(status_flags::dynallocated);
+	matrix & M = dynallocate<matrix>(rows, cols);
 
 	unsigned i = 0;
 	for (auto & itr : l) {
@@ -1593,8 +1571,7 @@ ex diag_matrix(const lst & l)
 	size_t dim = l.nops();
 
 	// Allocate and fill matrix
-	matrix &M = *new matrix(dim, dim);
-	M.setflag(status_flags::dynallocated);
+	matrix & M = dynallocate<matrix>(dim, dim);
 
 	unsigned i = 0;
 	for (auto & it : l) {
@@ -1610,8 +1587,7 @@ ex diag_matrix(std::initializer_list<ex> l)
 	size_t dim = l.size();
 
 	// Allocate and fill matrix
-	matrix &M = *new matrix(dim, dim);
-	M.setflag(status_flags::dynallocated);
+	matrix & M = dynallocate<matrix>(dim, dim);
 
 	unsigned i = 0;
 	for (auto & it : l) {
@@ -1624,8 +1600,8 @@ ex diag_matrix(std::initializer_list<ex> l)
 
 ex unit_matrix(unsigned r, unsigned c)
 {
-	matrix &Id = *new matrix(r, c);
-	Id.setflag(status_flags::dynallocated | status_flags::evaluated);
+	matrix & Id = dynallocate<matrix>(r, c);
+	Id.setflag(status_flags::evaluated);
 	for (unsigned i=0; i<r && i<c; i++)
 		Id(i,i) = _ex1;
 
@@ -1634,8 +1610,8 @@ ex unit_matrix(unsigned r, unsigned c)
 
 ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name)
 {
-	matrix &M = *new matrix(r, c);
-	M.setflag(status_flags::dynallocated | status_flags::evaluated);
+	matrix & M = dynallocate<matrix>(r, c);
+	M.setflag(status_flags::evaluated);
 
 	bool long_format = (r > 10 || c > 10);
 	bool single_row = (r == 1 || c == 1);
@@ -1676,8 +1652,8 @@ ex reduced_matrix(const matrix& m, unsigned r, unsigned c)
 
 	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);
+	matrix & M = dynallocate<matrix>(rows, cols);
+	M.setflag(status_flags::evaluated);
 
 	unsigned ro = 0;
 	unsigned ro2 = 0;
@@ -1705,8 +1681,8 @@ 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);
+	matrix & M = dynallocate<matrix>(nr, nc);
+	M.setflag(status_flags::evaluated);
 
 	for (unsigned ro=0; ro<nr; ++ro) {
 		for (unsigned co=0; co<nc; ++co) {