File:  [CENS] / python / pyGiNaC / wrappers4 / matrix.cpp
Revision 1.3: download - view: text, annotated - select for diffs - revision graph
Mon Dec 9 12:14:35 2002 UTC (14 years, 11 months ago) by pearu
Branches: MAIN
CVS tags: HEAD
*** empty log message ***

/*
  This file is part of the PyGiNaC package.
  http://cens.ioc.ee/projects/pyginac/

  $Revision: 1.3 $
  $Id: matrix.cpp,v 1.3 2002-12-09 12:14:35 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.
*/

#ifndef PYGINAC_MATRIX_CPP
#define PYGINAC_MATRIX_CPP
/* prototypes */

GiNaC::ex matrix0(unsigned r,
		  unsigned c,
		  const GiNaC::ex & init);
GiNaC::ex matrix(unsigned r,
		 unsigned c, 
		 const GiNaC::exvector & v);
GiNaC::ex matrix_mul(const GiNaC::matrix & self,
		     const GiNaC::matrix & other);
GiNaC::ex matrix_mul_numeric(const GiNaC::matrix & self,
			     const GiNaC::numeric & other);
GiNaC::ex matrix_determinant(const GiNaC::matrix & self,
			     const std::string & algo);
GiNaC::ex matrix_solve(const GiNaC::matrix & self,
		       const GiNaC::matrix & vars,
		       const GiNaC::matrix & rhs,
		       const std::string & algo);
GiNaC::ex matrix_get_slice(const GiNaC::matrix & self,
			   PyObject *index);
void matrix_set_slice(GiNaC::matrix & self,
		      PyObject *index,
		      const GiNaC::ex &other);

BOOST_PYTHON_BEGIN_CONVERSION_NAMESPACE

GiNaC::matrix & from_python(PyObject* o,
			    py::type<GiNaC::matrix &>);
const GiNaC::matrix & from_python(PyObject* o,
				  py::type<const GiNaC::matrix &>);
PyObject * to_python(const GiNaC::matrix & obj);

BOOST_PYTHON_END_CONVERSION_NAMESPACE

#else
#ifndef PYGINAC_MATRIX_CPP_1
#define PYGINAC_MATRIX_CPP_1
/* definitions */

this_module.def(matrix, "matrix");
this_module.def(matrix0, "matrix0");
ex_class.def(&GiNaC::matrix::rows, "matrix_rows");
ex_class.def(&GiNaC::matrix::cols, "matrix_cols");
ex_class.def(&GiNaC::matrix::trace, "matrix_trace");
ex_class.def(&GiNaC::matrix::transpose, "matrix_transpose");
ex_class.def(&GiNaC::matrix::inverse, "matrix_inverse");
ex_class.def(&GiNaC::matrix::charpoly, "matrix_charpoly");
ex_class.def(&GiNaC::matrix::add, "matrix_add");
ex_class.def(&GiNaC::matrix::sub, "matrix_sub");
ex_class.def(matrix_mul, "matrix_mul");
ex_class.def(matrix_mul_numeric, "matrix_mul_numeric");
ex_class.def(&GiNaC::matrix::mul_scalar, "matrix_mul_scalar");
ex_class.def(&GiNaC::matrix::pow, "matrix_pow");
ex_class.def(matrix_determinant, "matrix_determinant");
ex_class.def(matrix_solve, "matrix_solve");
ex_class.def(matrix_get_slice, "matrix_get_slice");
ex_class.def(matrix_set_slice, "matrix_set_slice");

#else
/* implementation */

GiNaC::ex matrix(unsigned r,
		 unsigned c,
		 const GiNaC::exvector & v) {
  return GiNaC::matrix(r,c,v);
}

GiNaC::ex matrix0(unsigned r,
		  unsigned c,
		  const GiNaC::ex & init) {
  GiNaC::exvector v;
  for (int i=r*c; i; --i) {
    v.push_back(init);
  }
  return GiNaC::matrix(r,c,v);
}

GiNaC::ex matrix_mul(const GiNaC::matrix & self,
		     const GiNaC::matrix & other) {
  return self.mul(other);
}

GiNaC::ex matrix_mul_numeric(const GiNaC::matrix & self,
			     const GiNaC::numeric & other) {
  return self.mul(other);
}

unsigned get_determinant_algo(const std::string & algo) {
  if (algo=="automatic")
    return GiNaC::determinant_algo::automatic;
  if (algo=="gauss")
    return GiNaC::determinant_algo::gauss;
  if (algo=="divfree")
    return GiNaC::determinant_algo::divfree;
  if (algo=="laplace")
    return GiNaC::determinant_algo::laplace;
  if (algo=="bareiss")
    return GiNaC::determinant_algo::bareiss;
  PyErr_SetString(PyExc_ValueError,(std::string("determinant_algo must be automatic|gauss|divfree|laplace|bareiss but got ")+algo).c_str());
  throw py::error_already_set();
}

unsigned get_solve_algo(const std::string & algo) {
  if (algo=="automatic")
    return GiNaC::solve_algo::automatic;
  if (algo=="gauss")
    return GiNaC::solve_algo::gauss;
  if (algo=="divfree")
    return GiNaC::solve_algo::divfree;
  if (algo=="bareiss")
    return GiNaC::solve_algo::bareiss;
  PyErr_SetString(PyExc_ValueError,(std::string("solve_algo must be automatic|gauss|divfree|bareiss but got ")+algo).c_str());
  throw py::error_already_set();
}

GiNaC::ex matrix_determinant(const GiNaC::matrix & self,
			     const std::string & algo) {
  return self.determinant(get_determinant_algo(algo));
}

GiNaC::ex matrix_solve(const GiNaC::matrix & self,
		       const GiNaC::matrix & vars,
		       const GiNaC::matrix & rhs,
		       const std::string & algo) {
  return self.solve(vars,rhs,get_solve_algo(algo));
}

BOOST_PYTHON_BEGIN_CONVERSION_NAMESPACE

GiNaC::matrix & from_python(PyObject* o,
			    py::type<GiNaC::matrix &>) {
  if (PYOBJ_IS_EX(o)) {
    GiNaC::ex & e = pyobj_get_ex_nonconst(o);
    if (GiNaC::is_a<GiNaC::matrix>(e))
      return GiNaC::ex_to_nonconst<GiNaC::matrix>(e);
  }
  PyErr_SetString(PyExc_TypeError,(std::string("from_python(obj, matrix) expected matrix but got ")+PyString_AsString(PyObject_Repr(PyObject_Type(o)))).c_str());
  throw py::error_already_set();
}

const GiNaC::matrix & from_python(PyObject* o,
				  py::type<const GiNaC::matrix &>) {
  return from_python(o, py::type<GiNaC::matrix &>());
}

PyObject * to_python(const GiNaC::matrix & obj) {
  return BOOST_PYTHON_CONVERSION::to_python(GiNaC::ex(obj));
}

BOOST_PYTHON_END_CONVERSION_NAMESPACE

GiNaC::ex matrix_get_slice(const GiNaC::matrix & self,PyObject *index) {
  int rows = self.rows();
  int cols = self.cols();
  int start_r, step_size_r, n_steps_r;
  int start_c, step_size_c, n_steps_c;
  if (PyTuple_Check(index)) {
    if (PyTuple_Size(index)==2) {
      start_r = parse_subindex(PyTuple_GetItem(index,0), 
			       &step_size_r, &n_steps_r, rows);
      start_c = parse_subindex(PyTuple_GetItem(index,1), 
			       &step_size_c, &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_get_slice() for PseudoIndex|RubberIndex");
	throw py::error_already_set();
      }
      if (n_steps_r==SingleIndex && n_steps_c==SingleIndex)
	return self(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::exvector v;
      for(int i=0,kr=start_r; 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)
	  v.push_back(self((kr==rows)?0:kr,(kc==cols)?0:kc));
      return GiNaC::matrix(n_steps_r,n_steps_c,v);
    } else {
      PyErr_SetString(PyExc_IndexError, "matrix_get_slice() index must be 2-tuple of int|slice or int|slice");
      throw py::error_already_set();
    }
  } else {
    start_c = parse_subindex(index, &step_size_c, &n_steps_c, cols);
    if (start_c < 0)
      throw py::error_already_set();
    if (n_steps_c==PseudoIndex || n_steps_c==RubberIndex) {
      PyErr_SetString(PyExc_NotImplementedError, "matric_get_slice() for PseudoIndex|RubberIndex");
      throw py::error_already_set();
    }
    if (n_steps_c==SingleIndex)
      n_steps_c=1;
    GiNaC::exvector v;
    for(int i=0;i<rows;++i)
      for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c)
	v.push_back(self(i,(kc==cols)?0:kc));
    return GiNaC::matrix(rows,n_steps_c,v);
  }
}

void matrix_set_slice(GiNaC::matrix & self,
		      PyObject *index,
		      const GiNaC::ex &other) {
  int rows = self.rows();
  int cols = self.cols();
  int start_r, step_size_r, n_steps_r;
  int start_c, step_size_c, n_steps_c;
  int n = other.nops();
  if (PyTuple_Check(index)) {
    if (PyTuple_Size(index)==2) {
      start_r = parse_subindex(PyTuple_GetItem(index,0), 
			       &step_size_r, &n_steps_r, rows);
      start_c = parse_subindex(PyTuple_GetItem(index,1), 
			       &step_size_c, &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_set_slice() for PseudoIndex|RubberIndex");
	throw py::error_already_set();
      }
      if (n_steps_r==SingleIndex && n_steps_c==SingleIndex) {
	self(start_r,start_c) = other;
	return;
      }
      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 (n==0) {
	for(int i=0,kr=start_r; 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)
	    self((kr==rows)?0:kr,(kc==cols)?0:kc) = other;
	return;
      }
      if (n!=n_steps_r*n_steps_c) {
	PyErr_SetString(PyExc_IndexError, "matrix_set_slice(): not matching index vectors");
	throw py::error_already_set();
      }
      for(int i=0,l=0,kr=start_r; 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)
	    self((kr==rows)?0:kr,(kc==cols)?0:kc) = other.op(l);
      return;
    } else {
      PyErr_SetString(PyExc_IndexError, "matrix_set_slice() index must be 2-tuple of int|slice or int|slice");
      throw py::error_already_set();
    }
  } else {
    start_c = parse_subindex(index, &step_size_c, &n_steps_c, cols);
    if (start_c < 0)
      throw py::error_already_set();
    if (n_steps_c==PseudoIndex || n_steps_c==RubberIndex) {
      PyErr_SetString(PyExc_NotImplementedError, "matric_set_slice() for PseudoIndex|RubberIndex");
      throw py::error_already_set();
    }
    if (n_steps_c==SingleIndex)
      n_steps_c=1;
    if (n==0) {
      for(int i=0;i<rows;++i)
	for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c)
	  self(i,(kc==cols)?0:kc) = other;
      return;
    }
    if (n!=rows*n_steps_c) {
      PyErr_SetString(PyExc_IndexError, "matrix_set_slice(): not matching index vectors");
      throw py::error_already_set();
    }
    for(int i=0,l=0;i<rows;++i)
      for(int j=0, kc=start_c;j<n_steps_c;++j,kc += step_size_c,++l)
	self(i,(kc==cols)?0:kc) = other.op(l);
    return;
  }
}

#endif
#endif

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