include/algebramix/series_matrix.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : series_matrix.hpp
* DESCRIPTION: Vectors and matrices of univariate power series
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_SERIES_MATRIX_HPP
#define __MMX_SERIES_MATRIX_HPP
#include <basix/vector.hpp>
#include <algebramix/matrix.hpp>
#include <algebramix/series.hpp>
#include <algebramix/series_vector.hpp>

namespace mmx {
#define TMPL template<typename C, typename V, typename U>
#define TMPLWU template<typename C, typename V, typename W, typename U>
#define Series series<C,V>
#define Matrix matrix<C,U>
#define Series_rep series_rep<C,V>
#define Series_matrix series<Matrix,V>
#define Matrix_series matrix<Series,U>
#define Polynomial polynomial<C, typename series_polynomial_helper<C,V>::PV>

/******************************************************************************
* Variant
******************************************************************************/

template<typename C, typename V>
struct series_variant_helper<matrix<C,V> > {
  typedef typename series_variant_helper<C>::SV SV;
};

/******************************************************************************
* Matrices of series
******************************************************************************/

STYPE_TO_TYPE(TMPL,as_matrix_type,Series_matrix,Matrix_series);

TMPL
class matrix_access_series_rep: public Series_rep {
  Series_matrix f;
  nat i, j;
public:
  matrix_access_series_rep (const Series_matrix& f2, nat i2, nat j2):
    Series_rep (get_format1 (CF(f2))), f (f2), i (i2), j (j2) {}
  syntactic expression (const syntactic& z) const {
    return access (flatten (f, z), flatten (i), flatten (j)); }
  void Increase_order (nat l) {
    Series_rep::Increase_order (l);
    increase_order (f, l); }
  C next () { return f[this->n](i,j); }
};

TMPL Series
access (const Series_matrix& f, nat i, nat j) {
  return (Series_rep*) new matrix_access_series_rep<C,V,U> (f, i, j);
}

TMPL Matrix_series
as_matrix (const Series_matrix& f) {
  nat nr= rows (f[0]), nc= cols (f[0]);
  Matrix_series m (Series (0), nr, nc);
  for (nat i=0; i<nr; i++)
    for (nat j=0; j<nc; j++)
      m (i, j)= access (f, i, j);
  return m;
}

TMPL syntactic
flatten (const Matrix_series& m, const syntactic& z) {
  syntactic g;
  nat i, j, nr= rows (m), nc= cols (m);
  matrix<syntactic,U> r (g, nr, nc);
  for (i=0; i<nr; i++)
    for (j=0; j<nc; j++)
      r (i, j)= flatten (m (i, j), z);
  return flatten (r);
}

TMPL format<Matrix >
get_matrix_format (const Matrix_series& m) {
  format<Series > fm1= CF(m);
  format<C> fm2= get_format1 (fm1);
  C zero= promote (0, fm2);
  Matrix r (zero, rows (m), cols (m));
  return get_format (r);
}

TMPL
class matrix_series_rep: public series_rep<Matrix,V> {
  const Matrix_series m;
public:
  matrix_series_rep (const Matrix_series& m2):
    series_rep<Matrix,V> (get_matrix_format (m2)), m (m2) {}
  syntactic expression (const syntactic& z) const {
    return flatten (m, z); }
  void Increase_order (nat l) {
    series_rep<Matrix,V>::Increase_order (l);
    for (nat i=0; i<rows(m); i++)
      for (nat j=0; j<cols(m); j++)
        increase_order (m (i, j), l); }
  Matrix next () {
    nat i, j, nr= rows (m), nc= cols (m);
    Matrix coeff (C (0), nr, nc);
    for (i=0; i<nr; i++)
      for (j=0; j<nc; j++)
        coeff (i, j)= m (i, j) [this->n];
    return coeff; }
};

TMPL Series_matrix
as_series (const Matrix_series& m) {
  return (series_rep<Matrix,V>*) new matrix_series_rep<C,V,U> (m);
}

TMPL Series_matrix
from_matrix (const Matrix_series& m) {
  return (series_rep<Matrix,V>*) new matrix_series_rep<C,V,U> (m);
}

TMPL matrix<Polynomial >
range (const Matrix_series& m, nat start, nat end) {
  C zero= promote (0, get_format1 (CF (m)));
  matrix<Polynomial > r (zero, rows (m), cols (m));
  for (nat i=0; i<rows(m); i++)
    for (nat j=0; j<cols(m); j++)
      r (i, j)= range (m (i, j), start, end);
  return r;
}

/******************************************************************************
* Linear differential matrix equations
******************************************************************************/

TMPL Series_matrix
solve_matrix_lde_init (const Series_matrix& f, const Matrix& c) {
  return unary_recursive_series<solve_matrix_lde_op> (f, c);
}

TMPL Matrix_series
solve_lde (const Matrix_series& f) {
  ASSERT (is_square_matrix (f), "square matrix expected");
  Matrix c (C(1), rows (f), cols (f));
  return as_matrix (solve_matrix_lde_init (as_series (f), c));
}

TMPL Matrix_series
solve_lde_init (const Matrix_series& f, const Matrix& c) {
  ASSERT (is_square_matrix (f), "square matrix expected");
  ASSERT (rows (f) == rows (c), "unequal matrix dimensions");
  return as_matrix (solve_matrix_lde_init (as_series (f), c));
}

#undef Polynomial
#undef Series_matrix
#undef Matrix_series
#undef Matrix
#undef Series
#undef Series_rep
#undef TMPLWU
#undef TMPL
} // namespace mmx
#endif // __MMX_SERIES_MATRIX_HPP

---