form_concat.h

Go to the documentation of this file.

#ifndef _RHEOLEF_FORM_CONCAT_H
#define _RHEOLEF_FORM_CONCAT_H
#include "rheolef/form.h"
#include "rheolef/csr_concat.h"

namespace rheolef {

//  A = {a, b};                 // matrix & vector
template <class T, class M>
struct form_concat_value {
 typedef enum { scalar, field, field_transpose, form} variant_type;
 form_concat_value (const T& x)               : s(x), v(), m(),  variant(scalar) {}
 form_concat_value (const form_basic<T,M>& x) : s(),  v(), m(x), variant(form) {}
// io/debug:
 friend std::ostream& operator<< (std::ostream& o, const form_concat_value<T,M>& x) {
       if (x.variant == scalar) return o << "s";
  else if (x.variant == field)  return o << "f";
  else if (x.variant == field_transpose) return o << "ft";
  else return o << "m";
 }
public:
 T                 s;
 field_basic<T,M>  v;
 form_basic<T,M>   m;
 variant_type      variant;
};
template <class T, class M>
struct form_concat_line {
 typedef typename form_basic<T,M>::size_type   size_type;
 typedef form_concat_value<T,M>                value_type;
 typedef typename std::list<value_type>::const_iterator const_iterator;
 form_concat_line () : _l() {}

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 form_concat_line (const std::initializer_list<value_type>& il) : _l() {
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef typename std::initializer_list<value_type>::const_iterator const_iterator;
#else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef const value_type* const_iterator;
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    for(const_iterator iter = il.begin(); iter != il.end(); ++iter) {
        _l.push_back(*iter);
    }
 }
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

 const_iterator begin() const { return _l.begin(); }
 const_iterator end()   const { return _l.end(); }

 friend std::ostream& operator<< (std::ostream& o, const form_concat_line<T,M>& x) {
    std::cout << "{";
    for(typename std::list<value_type>::const_iterator iter = x._l.begin(); iter != x._l.end(); ++iter) {
        std::cout << *iter << " ";
    }
    return std::cout << "}";
 }

 void build_form_pass0 (std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Yh, size_t i_comp = 0) const;
 static void build_first_space (const std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Xh);
 void build_form_pass1 (space_basic<T,M>& Xh, space_basic<T,M>& Yh) const;
 form_basic<T,M> build_form_pass2 (const space_basic<T,M>& Xh, const space_basic<T,M>& Yh) const;
 form_basic<T,M> build_form () const;

protected:
 std::list<value_type> _l;
};
template <class T, class M>
void
form_concat_line<T,M>::build_form_pass0 (std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Yh, size_t i_comp) const
{
  size_t j_comp = 0;
  bool have_Yh = false;
  typename std::vector<std::pair<bool,space_basic<T,M> > >::iterator xh_iter = l_Xh.begin();
  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, xh_iter++, j_comp++) {
    const value_type& x = *iter;
    switch (x.variant) {
      case form_concat_value<T,M>::scalar: {
        check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation"
                   << " at ("<<i_comp<<","<<j_comp<<")");
        break;
      }
      case form_concat_value<T,M>::form: {
        if (!(*xh_iter).first) {
          (*xh_iter).first  = true;
          (*xh_iter).second = x.m.get_first_space();
        } else {
          check_macro (x.m.get_first_space() == (*xh_iter).second, "form initializer: invalid second space `"
                << x.m.get_first_space().stamp() << "': expect `" << (*xh_iter).second.stamp() << "'"
                << " at ("<<i_comp<<","<<j_comp<<")");
        }
        if (!have_Yh) {
          have_Yh = true;
          Yh      = x.m.get_second_space();
        } else {
          check_macro (x.m.get_second_space() == Yh, "form initializer: invalid second space `"
                << x.m.get_second_space().stamp() << "': expect `" << Yh.stamp() << "'"
                << " at ("<<i_comp<<","<<j_comp<<")");
        }
        break;
      }
      default: error_macro ("non-form or scalar concatenation not yet supported"
                << " at ("<<i_comp<<","<<j_comp<<")");
    }
  }
  check_macro (have_Yh, "form concatenation: "<<i_comp<<"th row space remains undefined");
}
template <class T, class M>
void
form_concat_line<T,M>::build_first_space (const std::vector<std::pair<bool,space_basic<T,M> > >& l_Xh, space_basic<T,M>& Xh)
{
  space_mult_list<T,M> sml_X;
  size_t j_comp = 0;
  for (typename std::vector<std::pair<bool,space_basic<T,M> > >::const_iterator
        xh_iter = l_Xh.begin(),
        xh_last = l_Xh.end(); xh_iter != xh_last; xh_iter++, j_comp++) {
    check_macro ((*xh_iter).first, "form concatenation: "<<j_comp<<"th column space remains undefined");
    sml_X *= (*xh_iter).second;
  }
  Xh = space_basic<T,M>(sml_X);
}
template <class T, class M>
void
form_concat_line<T,M>::build_form_pass1 (space_basic<T,M>& Xh, space_basic<T,M>& Yh) const
{
  std::vector<std::pair<bool,space_basic<T,M> > > l_Xh (_l.size(), std::pair<bool,space_basic<T,M> >(false, space_basic<T,M>()));
  build_form_pass0 (l_Xh, Yh);
  build_first_space (l_Xh, Xh);
}
template <class T, class M>
form_basic<T,M>
form_concat_line<T,M>::build_form_pass2 (const space_basic<T,M>& Xh, const space_basic<T,M>& Yh) const
{
  form_basic<T,M> a (Xh, Yh);
  csr_concat_line<T,M> uu, ub, bu, bb;
  for(typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
    const value_type& x = *iter;
    switch (x.variant) {
      case form_concat_value<T,M>::form: {
        uu.push_back (x.m.uu());
        ub.push_back (x.m.ub());
        bu.push_back (x.m.bu());
        bb.push_back (x.m.bb());
        break;
      }
      default: error_macro ("non-form concatenation not yet supported");
    }
  }
  a.set_uu() = uu.build_csr();
  a.set_ub() = ub.build_csr();
  a.set_bu() = bu.build_csr();
  a.set_bb() = bb.build_csr();
  return a;
}
template <class T, class M>
form_basic<T,M>
form_concat_line<T,M>::build_form() const
{
  space_basic<T,M> Xh, Yh;
  build_form_pass1 (Xh, Yh);
  return build_form_pass2 (Xh, Yh);
}
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
template <class T, class M>
inline
form_basic<T,M>::form_basic (const std::initializer_list<form_concat_value<T,M> >& init_list)
{
  form_concat_line<T,M> cc (init_list);
  form_basic<T,M>::operator= (cc.build_form());
}
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
template <class T, class M>
struct form_concat {
 typedef typename form_basic<T,M>::size_type   size_type;
 typedef form_concat_line<T,M>           line_type;
 typedef form_concat_value<T,M>          value_type;
 form_concat () : _l() {}

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 form_concat (const std::initializer_list<line_type>& il) : _l() {
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef typename std::initializer_list<line_type>::const_iterator const_iterator;
#else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef const line_type* const_iterator;
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    for(const_iterator iter = il.begin(); iter != il.end(); ++iter) {
        _l.push_back(*iter);
    }
 }
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

 friend std::ostream& operator<< (std::ostream& o, const form_concat<T,M>& x) {
    std::cout << "{";
    for(typename std::list<line_type>::const_iterator iter = x._l.begin(); iter != x._l.end(); ++iter) {
        std::cout << *iter << " ";
    }
    return std::cout << "}";
 }
 form_basic<T,M> build_form () const;

protected:
 std::list<line_type> _l;
};
template <class T, class M>
form_basic<T,M>
form_concat<T,M>::build_form() const
{
  size_t i_comp = 0;
  space_mult_list<T,M> sml_Y;
  std::vector<std::pair<bool,space_basic<T,M> > > l_Xh (_l.size(), std::pair<bool,space_basic<T,M> >(false, space_basic<T,M>()));
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
    const line_type& line = *iter;
    space_basic<T,M> Yih;
    line.build_form_pass0 (l_Xh, Yih, i_comp);
    sml_Y *= Yih;
  }
  space_basic<T,M> Xh;
  form_concat_line<T,M>::build_first_space (l_Xh, Xh);
  space_basic<T,M> Yh (sml_Y);
  csr_concat<T,M> uu, ub, bu, bb;
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter) {
    const line_type& line = *iter;
    csr_concat_line<T,M> uu_i, ub_i, bu_i, bb_i;
    for (typename std::list<value_type>::const_iterator jter = line.begin(); jter != line.end(); ++jter) {
      const value_type& x = *jter;
      switch (x.variant) {
        case form_concat_value<T,M>::scalar: {
          check_macro (x.s == 0, "unsupported non-nul scalar `"<<x.s<<"' in form concatenation");
          // zero: no values to insert in the sparse structure => nothing more to do
          uu_i.push_back (x.s);
          ub_i.push_back (x.s);
          bu_i.push_back (x.s);
          bb_i.push_back (x.s);
          break;
        }
        case form_concat_value<T,M>::form: {
          uu_i.push_back (x.m.uu());
          ub_i.push_back (x.m.ub());
          bu_i.push_back (x.m.bu());
          bb_i.push_back (x.m.bb());
          break;
        }
        default: error_macro ("non-form or scalar concatenation not yet supported");
      }
    }
    uu.push_back (uu_i);
    ub.push_back (ub_i);
    bu.push_back (bu_i);
    bb.push_back (bb_i);
  }
  form_basic<T,M> a(Xh, Yh);
  a.set_uu() = uu.build_csr();
  a.set_ub() = ub.build_csr();
  a.set_bu() = bu.build_csr();
  a.set_bb() = bb.build_csr();
  return a;
}

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
template <class T, class M>
inline
form_basic<T,M>::form_basic (const std::initializer_list<form_concat_line<T,M> >& init_list)
{
  form_concat<T,M> cc (init_list);
  form_basic<T,M>::operator= (cc.build_form());
}
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

} // namespace rheolef
#endif // _RHEOLEF_FORM_CONCAT_H