File:  [CENS] / python / pyGiNaC / wrappers3 / matrix_py.cpp
Revision 1.3: download - view: text, annotated - select for diffs - revision graph
Thu Apr 26 19:44:48 2001 UTC (16 years, 7 months ago) by pearu
Branches: MAIN
CVS tags: HEAD
Impl. pseries,ncmul. Preparing for a release.

/*
# This file is part of the PyGiNaC package.
# http://cens.ioc.ee/projects/pyginac/
#
# $Revision: 1.3 $
# $Id: matrix_py.cpp,v 1.3 2001-04-26 19:44:48 pearu Exp $
#
# Copyright 2001 Pearu Peterson all rights reserved,
# Pearu Peterson <pearu@cens.ioc.ee>
# Permission to use, modify, and distribute this software is given under the
# terms of the LGPL.  See http://www.fsf.org
#
# NO WARRANTY IS EXPRESSED OR IMPLIED.  USE AT YOUR OWN RISK.
#
*/
/*DT
  Matrix
  ------

  >>> from ginac import matrix,is_matrix,numeric

  Column vector:
  >>> print matrix([1,2])
  [[1],
   [2]]

  Row vector:
  >>> print matrix([[1,2]])
  [[1, 2]]

  Matrix:
  >>> a = matrix([[1,2],[3,4]])
  >>> print a
  [[1, 2],
   [3, 4]]
  >>> a
  matrix([[numeric('1'), numeric('2')],
          [numeric('3'), numeric('4')]])

  Take a row:
  >>> print a[1]
  [[3, 4]]

  Take a column:
  >>> print a[:,1]
  [[2],
   [4]]

  Take an element:
  >>> print a[0,1]
  2
  
  Auxiliary functions:
  >>> a.is_matrix(),is_matrix(a),numeric().is_matrix(),is_matrix(numeric())
  (1, 1, 0, 0)

  Matrix constructors:
  >>> print matrix(2,3)
  [[0, 0, 0],
   [0, 0, 0]]
  >>> print matrix(2,3,range(6))
  [[0, 1, 2],
   [3, 4, 5]]
  >>> print matrix(2,3,[1,2,3])
  [[1, 2, 3],
   [0, 0, 0]]
 
  Slicing:
  >>> a = matrix(3,4,range(12))
  >>> print a
  [[0, 1, 2, 3],
   [4, 5, 6, 7],
   [8, 9, 10, 11]]
  >>> print a[1:]
  [[4, 5, 6, 7],
   [8, 9, 10, 11]]
  >>> print a[:,1:]
  [[1, 2, 3],
   [5, 6, 7],
   [9, 10, 11]]
  >>> print a[0:2, 1:3]
  [[1, 2],
   [5, 6]]
  >>> print a[:,::-1]
  [[3, 2, 1, 0],
   [7, 6, 5, 4],
   [11, 10, 9, 8]]
  >>> print a[:,::2]
  [[0, 2],
   [4, 6],
   [8, 10]]
  >>> print a[1,1:-1]
  [[5, 6]]

  Modifying in-place:
  >>> a = matrix(3,3)
  >>> a[1,1] = 8        # middle element
  >>> a[-1] = [3,2,1]   # last row
  >>> a[:,-1] = [5,4,1] # last column
  >>> print a
  [[0, 0, 5],
   [0, 8, 4],
   [3, 2, 1]]
  >>> a = matrix(4,4)   # 4x4 matrix of zeros
  >>> a[1:3,1:3] = [[1,2],[3,4]]
  >>> print a
  [[0, 0, 0, 0],
   [0, 1, 2, 0],
   [0, 3, 4, 0],
   [0, 0, 0, 0]]

  Functions:
  >>> from ginac import rows, cols, transpose
  >>> rows(a), cols(a), a.rows(), a.cols()
  (4, 4, 4, 4)
  >>> print transpose(a)
  [[0, 0, 0, 0],
   [0, 1, 3, 0],
   [0, 2, 4, 0],
   [0, 0, 0, 0]]
  >>> print a.transpose()
  [[0, 0, 0, 0],
   [0, 1, 3, 0],
   [0, 2, 4, 0],
   [0, 0, 0, 0]]

  >>> from ginac import determinant,trace,inverse,charpoly
  >>> a = matrix(2,2,[1,2,3,4])
  >>> print a
  [[1, 2],
   [3, 4]]
  >>> print a.determinant(), a.determinant('gauss')
  -2 -2
  >>> print determinant(a), determinant(a,'divfree')
  -2 -2
  >>> print a.trace(), trace(a)
  5 5
  >>> print a.inverse()
  [[-2, 1],
   [3/2, -1/2]]
  >>> print a*inverse(a)
  [[1, 0],
   [0, 1]]
  >>> print charpoly(a,'lambda')
  -2 - 5 * lambda + lambda ** 2
  >>> from ginac import symbol
  >>> print a.charpoly(symbol('lambda'))
  -2 - 5 * lambda + lambda ** 2
  >>> syms = [symbol('a1'),symbol('a2')]
  >>> print a.solve(syms,[5,-6])
  [[-16],
   [21/2]]
  >>> print a*a.solve(syms,[5,-6])
  [[5],
   [-6]]

  >>> from ginac import diag_matrix
  >>> print diag_matrix([1,2,3])
  [[1, 0, 0],
   [0, 2, 0],
   [0, 0, 3]]
*/


#ifdef PYGINAC_DEFS
this_module.def(&ex::is_matrix, "is_matrix");
this_module.def(&return_false, "is_matrix");
ex_class.def(&ex::is_matrix, "is_matrix");
this_module.def(&matrix_1, "matrix");
this_module.def(&matrix_2, "matrix");
this_module.def(&matrix_3, "matrix");

ex_class.def(&rows_0, "rows");
this_module.def(&rows_0, "rows");
ex_class.def(&cols_0, "cols");
this_module.def(&cols_0, "cols");
ex_class.def(&transpose_0, "transpose");
this_module.def(&transpose_0, "transpose");
ex_class.def(&determinant_0, "determinant");
this_module.def(&determinant_0, "determinant");
ex_class.def(&determinant_1, "determinant");
this_module.def(&determinant_1, "determinant");
if ((determinant_algo_map = PyDict_New())==NULL)
  throw py::error_already_set();
PyDict_SetItemString(determinant_algo_map,"automatic",PyInt_FromLong(0));
PyDict_SetItemString(determinant_algo_map,"gauss",PyInt_FromLong(1));
PyDict_SetItemString(determinant_algo_map,"divfree",PyInt_FromLong(2));
PyDict_SetItemString(determinant_algo_map,"laplace",PyInt_FromLong(3));
PyDict_SetItemString(determinant_algo_map,"bareiss",PyInt_FromLong(4));
ex_class.def(&trace_0, "trace");
this_module.def(&trace_0, "trace");
ex_class.def(&inverse_0, "inverse");
this_module.def(&inverse_0, "inverse");
ex_class.def(&charpoly_1, "charpoly");
this_module.def(&charpoly_1, "charpoly");
if ((solve_algo_map = PyDict_New())==NULL)
  throw py::error_already_set();
PyDict_SetItemString(solve_algo_map,"automatic",PyInt_FromLong(0));
PyDict_SetItemString(solve_algo_map,"gauss",PyInt_FromLong(1));
PyDict_SetItemString(solve_algo_map,"divfree",PyInt_FromLong(2));
PyDict_SetItemString(solve_algo_map,"bareiss",PyInt_FromLong(3));
ex_class.def(&solve_2, "solve");
ex_class.def(&solve_3, "solve");
this_module.def(&GiNaC::diag_matrix, "diag_matrix");
#else
#ifdef PYGINAC_EX_PROTOS
#define PYGINAC_PROTOS
#endif
#ifdef PYGINAC_PROTOS
#ifdef PYGINAC_EX_PROTOS
bool is_matrix(void) const;
#else // PYGINAC_EX_PROTOS
#include "slice.h"
BOOST_PYTHON_BEGIN_CONVERSION_NAMESPACE
PYGINAC_TOPYTHON1(matrix);
GiNaC::matrix from_python (PyObject* o, py::type<const GiNaC::matrix &>);
//unsigned from_python (PyObject* o, py::type<GiNaC::determinant_algo>);
BOOST_PYTHON_END_CONVERSION_NAMESPACE
GiNaC::ex matrix_1(const GiNaC::matrix & m);
GiNaC::ex matrix_2(long r, long c);
GiNaC::ex matrix_3(long r, long c, const GiNaC::lst & l);
PyObject* matrix_getitem(const GiNaC::matrix & m, py::ref index);
void matrix_setitem(GiNaC::matrix & m, py::ref index, py::ref other);
unsigned rows_0(const GiNaC::matrix & m);
unsigned cols_0(const GiNaC::matrix & m);
GiNaC::matrix transpose_0(const GiNaC::matrix & m);
GiNaC::ex determinant_0(const GiNaC::matrix & m);
GiNaC::ex determinant_1(const GiNaC::matrix & m, py::ref algo);
PyObject * determinant_algo_map = NULL;
GiNaC::ex trace_0(const GiNaC::matrix & m);
GiNaC::matrix inverse_0(const GiNaC::matrix & m);
GiNaC::ex charpoly_1(const GiNaC::matrix & m, const GiNaC::symbol & s);
PyObject * solve_algo_map = NULL;
GiNaC::matrix solve_2(const GiNaC::matrix &m, const GiNaC::matrix &vars, const GiNaC::matrix &rhs);
GiNaC::matrix solve_3(const GiNaC::matrix &m, const GiNaC::matrix &vars, const GiNaC::matrix &rhs, py::ref algo);
#endif // !PYGINAC_EX_PROTOS
#else  // PYGINAC_PROTOS
BOOST_PYTHON_BEGIN_CONVERSION_NAMESPACE
GiNaC::matrix from_python (PyObject* o, py::type<const GiNaC::matrix &>) {
  if (ExInstance_Check(o)) {
    GiNaC::ex e = from_python(o, py::type<const GiNaC::ex &>());
    if (is_ex_exactly_of_type(e, matrix))
      return ex_to_matrix(e);
    if (is_ex_exactly_of_type(e, lst)) {
      GiNaC::ex ne = GiNaC::lst_to_matrix(ex_to_lst(e));
      if (is_ex_exactly_of_type(ne, matrix))
	return ex_to_matrix(ne);
    }
  }
  if (Sequence_Check(o)) {
	unsigned r;
	unsigned c;
	r = PySequence_Length(o);
	PyObject *o2 = PySequence_GetItem(o,0);
	if (Sequence_Check(o2))
	  c = PySequence_Length(o2);
	else
	  c = 1;
	GiNaC::exvector v(0);
	for (unsigned i=0;i<r;++i) {
	  PyObject *o2 = PySequence_GetItem(o,i);
	  if (!Sequence_Check(o2))
		if (c>1) {
		  PyErr_SetString(PyExc_TypeError, "matrix() sequence argument must be a sequence");
		  throw py::error_already_set();
		} else
		  v.push_back(ex_from_ref(o2));
	  else {
		if (!(unsigned(PySequence_Length(o2))==c)) {
		  PyErr_SetString(PyExc_TypeError, "matrix() sequence argument must be a sequence of equal length sequences");
		  throw py::error_already_set();
		}
		for (unsigned j=0;j<c;++j)
		  v.push_back(ex_from_ref(PySequence_GetItem(o2,j)));
	  }
	}
	return GiNaC::matrix(r,c,v);
  }
  PYGINAC_FROMPYTHON_TYPEERROR(matrix,(ex(matrix)|int,int|int,int,sequence));
}
BOOST_PYTHON_END_CONVERSION_NAMESPACE

/*F_DT is_matrix(obj)
  Check if `obj' is matrix.
*/
/*M_DT is_matrix(self)
  Check if object is matrix.
*/
bool ex::is_matrix(void) const { return is_ex_exactly_of_type(*this, matrix); }
/*F_DT matrix(*args)
  matrix(seq) - construct matrix from a 1- or 2-sequence `seq'.
    If `seq' is 1-sequence, then returned matrix is formed as a column matrix.
  matrix(r,c) - construct r x c matrix filled with zeros.
  matrix(r,c,seq) - construct r x c matrix filled with elements in a flat
    sequence `seq'
  matrix(matrix(...)) -> matrix(...)
*/
GiNaC::ex matrix_1(const GiNaC::matrix & m) {
  return m;
}
GiNaC::ex matrix_2(long r, long c) {
  if (r > 0 && c > 0)
    return *new GiNaC::matrix(r,c);
  PyErr_SetString(PyExc_ValueError, "matrix() number of rows and columns must be positive");
  throw py::error_already_set();
}
GiNaC::ex matrix_3(long r, long c, const GiNaC::lst & l) {
  if (r > 0 && c > 0)
    return *new GiNaC::matrix(r,c,l);
  PyErr_SetString(PyExc_ValueError, "matrix() number of rows and columns must be positive");
  throw py::error_already_set();
}
PyObject* matrix_getitem(const GiNaC::matrix & m, py::ref index) {
  PyObject * o = index.get();
  int rows = m.rows();
  int cols = m.cols();
  int start_r, step_size_r, n_steps_r;
  int start_c, step_size_c, n_steps_c;  
  if (PyTuple_Check(o))
    if (PyTuple_Size(o)==2) {
      start_r = parse_subindex(PyTuple_GetItem(o,0), &step_size_r, &n_steps_r, rows);
      start_c = parse_subindex(PyTuple_GetItem(o,1), &step_size_c, &n_steps_c, cols);
    } else {
      PyErr_SetString(PyExc_IndexError, "matrix_getitem() index must be 2-int|2-slice|int,slice|slice,int|int|slice");
      start_r = start_c = -1;
    }
  else {
    if (rows==1) {
      start_r = 0;
      step_size_r = 1;
      n_steps_r = 1;
      start_c = parse_subindex(o, &step_size_c, &n_steps_c, cols);
    } else {
      start_r = parse_subindex(o, &step_size_r, &n_steps_r, rows);
      start_c = 0;
      step_size_c = 1;
      n_steps_c = cols;
    }
  }
  if (start_r < 0 || start_c < 0)
    throw py::error_already_set();
  if (n_steps_r==PseudoIndex || n_steps_r==RubberIndex || n_steps_c==PseudoIndex || n_steps_c==RubberIndex) {
    PyErr_SetString(PyExc_NotImplementedError, "matrix_getitem() for PseudoIndex|RubberIndex");
    throw py::error_already_set();
  }
  if ((n_steps_r==SingleIndex && n_steps_c==SingleIndex) ||
	  ((n_steps_r==SingleIndex && cols==1) || (n_steps_c==SingleIndex && rows==1)))
    return BOOST_PYTHON_CONVERSION::to_python(m(start_r,start_c));
  if (n_steps_r==SingleIndex)
    n_steps_r = step_size_r = 1;
  if (n_steps_c==SingleIndex)
    n_steps_c = step_size_c = 1;
  GiNaC::matrix rm(n_steps_r,n_steps_c);
  int kr = start_r;
  for(int i=0;i<n_steps_r;++i,kr += step_size_r)
    for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c)
      rm.set(i,j,m((kr==rows)?0:kr,(kc==cols)?0:kc));
  return BOOST_PYTHON_CONVERSION::to_python(rm);
}
void matrix_setitem(GiNaC::matrix & m, py::ref index, py::ref other) {
  PyObject * o = index.get();
  PyObject * oo = other.get();
  int rows = m.rows();
  int cols = m.cols();
  int start_r, step_size_r, n_steps_r;
  int start_c, step_size_c, n_steps_c;
  if (PyTuple_Check(o))
    if (PyTuple_Size(o)==2) {
      start_r = parse_subindex(PyTuple_GetItem(o,0), &step_size_r, &n_steps_r, rows);
      start_c = parse_subindex(PyTuple_GetItem(o,1), &step_size_c, &n_steps_c, cols);
    } else {
      PyErr_SetString(PyExc_IndexError, "matrix_setitem() index must be 2-int|2-slice|int,slice|slice,int|int|slice");
      start_r = start_c = -1;
    }
  else {
    if (rows==1) {
      start_r = 0;
      step_size_r = 1;
      n_steps_r = 1;
      start_c = parse_subindex(o, &step_size_c, &n_steps_c, cols);
    } else {
      start_r = parse_subindex(o, &step_size_r, &n_steps_r, rows);
      start_c = 0;
      step_size_c = 1;
      n_steps_c = cols;
    }
  }
  if (start_r < 0 || start_c < 0)
    throw py::error_already_set();
  if (n_steps_r==PseudoIndex || n_steps_r==RubberIndex || n_steps_c==PseudoIndex || n_steps_c==RubberIndex) {
    PyErr_SetString(PyExc_NotImplementedError, "matrix_setitem() for PseudoIndex|RubberIndex");
    throw py::error_already_set();
  }
  if ((n_steps_r==SingleIndex && n_steps_c==SingleIndex) ||
      ((n_steps_r==SingleIndex && cols==1) || (n_steps_c==SingleIndex && rows==1)))
    m.set(start_r,start_c,ex_from_ref(other));
  else {
    Py_INCREF(oo);
    GiNaC::matrix om = BOOST_PYTHON_CONVERSION::from_python(oo, py::type<const GiNaC::matrix &>());
    int orows = om.rows();
    int ocols = om.cols();
    if (n_steps_r==SingleIndex)
      n_steps_r = step_size_r = 1;
    if (n_steps_c==SingleIndex)
      n_steps_c = step_size_c = 1;
    if ((orows==1 || ocols==1) && (n_steps_r==1 || n_steps_c==1))
      if (orows*ocols == n_steps_r*n_steps_c ) {
        int kr = start_r;
        for(int i=0,l=0;i<n_steps_r;++i,kr += step_size_r)
          for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c,++l)
            m.set((kr==rows)?0:kr,(kc==cols)?0:kc, om.op(l));
      } else {
        PyErr_SetString(PyExc_IndexError, "matrix_setitem() different sizes");
        throw py::error_already_set();
      }
    else {
      if (orows != n_steps_r || ocols != n_steps_c ) {
        PyErr_SetString(PyExc_IndexError, "matrix_setitem() not aligned");
        throw py::error_already_set();
      }
      int kr = start_r;
      for(int i=0;i<n_steps_r;++i,kr += step_size_r)
        for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c)
          m.set((kr==rows)?0:kr,(kc==cols)?0:kc, om(i,j));
    }
  }
}
/*M_DT rows(self)
  Get number of rows.
*/
/*F_DT rows(m)
  Get number of rows.
*/
unsigned rows_0(const GiNaC::matrix & m) {
  return m.rows();
}
/*M_DT cols(self)
  Get number of columns.
*/
/*F_DT cols(m)
  Get number of columns.
*/
unsigned cols_0(const GiNaC::matrix & m) {
  return m.cols();
}
/*M_DT transpose(self)
  Return transposed matrix.
*/
/*F_DT transpose(m)
  Return transposed matrix.
*/
GiNaC::matrix transpose_0(const GiNaC::matrix & m) {
  return m.transpose();
}
/*M_DT determinant(self, algo='automatic')
  Determinant of a square matrix.
  If all the elements of the matrix are elements of an integral
  domain the determinant is also in that integral domain and the
  result is expanded only. If one or more elements are from a
  quotient field the determinant is usually also in that quotient
  field and the result is normalized before it is returned. This
  implies that the determinant of the symbolic 2x2 matrix
  [[a/(a-b),1],[b/(a-b),1]] is returned as unity. (In this respect,
  it behaves like MapleV and unlike Mathematica.)

  Option algo:
    'automatic' - let the system choose  
    'gauss'     - Gauss elimiation
    'divfree'   - Division-free elimination
    'laplace'   - Laplace (or minor) elimination
    'bareiss'   - Bareiss fraction-free elimination
*/
/*F_DT determinant(m, algo='automatic')
  Determinant of a square matrix.
  If all the elements of the matrix are elements of an integral
  domain the determinant is also in that integral domain and the
  result is expanded only. If one or more elements are from a
  quotient field the determinant is usually also in that quotient
  field and the result is normalized before it is returned. This
  implies that the determinant of the symbolic 2x2 matrix
  [[a/(a-b),1],[b/(a-b),1]] is returned as unity. (In this respect,
  it behaves like MapleV and unlike Mathematica.)

  Option algo:
    'automatic' - let the system choose  
    'gauss'     - Gauss elimiation
    'divfree'   - Division-free elimination
    'laplace'   - Laplace (or minor) elimination
    'bareiss'   - Bareiss fraction-free elimination
*/
GiNaC::ex determinant_0(const GiNaC::matrix & m) {
  return m.determinant();
}
GiNaC::ex determinant_1(const GiNaC::matrix & m, py::ref algo) {
  unsigned algo_v = 999;
  PyObject * o = algo.get();
  if (PyString_Check(o)) {
    PyObject *v = PyDict_GetItem(determinant_algo_map,o);
    if (v!=NULL && PyInt_Check(v))
      switch (PyInt_AS_LONG(v)) {
      case 0: algo_v = GiNaC::determinant_algo::automatic; break;
      case 1: algo_v = GiNaC::determinant_algo::gauss; break;
      case 2: algo_v = GiNaC::determinant_algo::divfree; break;
      case 3: algo_v = GiNaC::determinant_algo::laplace; break;
      case 4: algo_v = GiNaC::determinant_algo::bareiss; break;
      }
  }
  if (algo_v == 999) {
    PyErr_SetString(PyExc_TypeError, "determinant_algo argument must be string: automatic|gauss|divfree|laplace|bareiss");
    throw py::error_already_set();
  }
  return m.determinant(algo_v);
}
/*M_DT trace(self)
  Trace of a square matrix.        
  The result is normalized if it is in some quotient field and
  expanded only otherwise. This implies that the trace of the
  symbolic 2x2 matrix [[a/(a-b),x],[y,b/(b-a)]] is recognized to be
  unity.
  Return the sum of diagonal elements.
*/
/*F_DT trace(m)
  Trace of a square matrix.        
  The result is normalized if it is in some quotient field and
  expanded only otherwise. This implies that the trace of the
  symbolic 2x2 matrix [[a/(a-b),x],[y,b/(b-a)]] is recognized to be
  unity.
  Return the sum of diagonal elements.
*/
GiNaC::ex trace_0(const GiNaC::matrix & m) {
  return m.trace();
}
/*M_DT inverse(self)
  Inverse of a square matrix.
*/
/*F_DT inverse(m)
  Inverse of a square matrix.
 */
GiNaC::matrix inverse_0(const GiNaC::matrix & m) {
  return m.inverse();
}
/*M_DT charpoly(self, sym)
  Characteristic Polynomial.
  Following mathematica notation the characteristic polynomial of a
  matrix M is defined as the determiant of `(M - sym * 1)' where 1
  stands for the unit matrix of the same dimension as M. Note that
  some CASs define it with a sign inside the determinant which gives
  rise to an overall sign if the dimension is odd. This method
  returns the characteristic polynomial collected in powers of
  `sym' as a new expression.
  `sym' must be a ex(symbol)|string.
*/
/*F_DT charpoly(m, sym)
  Characteristic Polynomial.
  Following mathematica notation the characteristic polynomial of a
  matrix M is defined as the determiant of `(M - sym * 1)' where 1
  stands for the unit matrix of the same dimension as M. Note that
  some CASs define it with a sign inside the determinant which gives
  rise to an overall sign if the dimension is odd. This method
  returns the characteristic polynomial collected in powers of
  `sym' as a new expression.
  `sym' must be a ex(symbol)|string.
 */
GiNaC::ex charpoly_1(const GiNaC::matrix & m, const GiNaC::symbol & s) {
  return m.charpoly(s);
}

/*M_DT solve(self,vars,rhs,algo='automatic')
  Solve a linear system consisting of a m x n matrix and
  a m x p right hand side by applying an elimination scheme to the
  augmented matrix.
  vars - n x p matrix, all elements must be symbols
  rhs  - m x p matrix
  Option algo:
    automatic - let the system choose  
    gauss     - Gauss elimiation
    divfree   - Division-free elimination
    bareiss   - Bareiss fraction-free elimination
  Return n x p solution matrix.
*/
GiNaC::matrix solve_2(const GiNaC::matrix &m, const GiNaC::matrix &vars, const GiNaC::matrix &rhs) {
  return m.solve(vars,rhs);
}
GiNaC::matrix solve_3(const GiNaC::matrix &m, const GiNaC::matrix &vars, const GiNaC::matrix &rhs, py::ref algo) {
  unsigned algo_v = 999;
  PyObject * o = algo.get();
  if (PyString_Check(o)) {
    PyObject *v = PyDict_GetItem(solve_algo_map,o);
    if (v!=NULL && PyInt_Check(v))
      switch (PyInt_AS_LONG(v)) {
      case 0: algo_v = GiNaC::solve_algo::automatic; break;
      case 1: algo_v = GiNaC::solve_algo::gauss; break;
      case 2: algo_v = GiNaC::solve_algo::divfree; break;
      case 3: algo_v = GiNaC::solve_algo::bareiss; break;
      }
  }
  if (algo_v == 999) {
    PyErr_SetString(PyExc_TypeError, "solve_algo argument must be string: automatic|gauss|divfree|bareiss");
    throw py::error_already_set();
  }  
  return m.solve(vars,rhs,algo_v);
}
/*F_DT diag_matrix(seq)
  Return a new matrix with a diagonal `seq'.
 */
#endif // !PYGINAC_PROTOS
#endif // !PYGINAC_DEFS





FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>