matrix.hpp

Go to the documentation of this file.

/*****************************************************************************
    
    matrix.hpp -- îïðåäåëåíèå ìàòðèöû.

    Ýòîò ôàéë ÿâëÿåòñÿ ÷àñòüþ áèáëèîòåêè Arageli.

    Copyright (C) Nikolai Yu. Zolotykh, 1999--2006
    Copyright (C) Sergey S. Lyalin, 2005
    University of Nizhni Novgorod, Russia

*****************************************************************************/

#ifndef _ARAGELI_matrix_hpp_
#define _ARAGELI_matrix_hpp_

#include "config.hpp"

#include <cstddef>
#include <iostream>
#include <vector>
#include <utility>

#include "frwrddecl.hpp"

#include "type_opers.hpp"
#include "type_traits.hpp"
#include "type_pair_traits.hpp"
#include "exception.hpp"
#include "refcntr.hpp"
#include "factory.hpp"
#include "cmp.hpp"
#include "io.hpp"
#include "vector.hpp"

#include "std_import.hpp"


namespace Arageli
{



struct nonsquare_t {};

//const nonsquare_t fromsize = nonsquare_t();



struct colwise_t {};


const colwise_t colwise = colwise_t();



struct diag_t {};


const diag_t diag = diag_t();



struct eye_t {};


const eye_t eye = eye_t();


struct frommat_t {};
const frommat_t frommat = frommat_t();


namespace _Internal
{

template <typename T>
struct Rep_matrix   // internal representation for matrix template
{
    // std::vector -- âðåìåííî, íóæíî ñâîé îáëåã÷¸ííûé âåêòîð
    typedef std::vector<T> Mem;
    typedef typename Mem::size_type size_type;
    size_type rows, cols;
    Mem mem;

    Rep_matrix () : rows(0), cols(0) {}

    void assign (size_type rows_a, size_type cols_a, const T& val)
    { rows = rows_a; cols = cols_a; mem.assign(rows*cols, val); }

    template <typename In>
    void assign (size_type rows_a, size_type cols_a, In b, In e)
    {
        rows = rows_a; cols = cols_a; mem.assign(b, e);
        ARAGELI_ASSERT_1(rows*cols == mem.size());
    }
};

} // namespace _Internal
} // namespace Arageli


namespace std
{

template <typename T>
inline void swap
(
    Arageli::_Internal::Rep_matrix<T>& a,
    Arageli::_Internal::Rep_matrix<T>& b
)
{
    swap(a.mem, b.mem);
    swap(a.rows, b.rows);
    swap(a.cols, b.cols);
}

} // namespace std


namespace Arageli
{



template
<
    typename T,
    bool REFCNT
>
class matrix
{
    template <typename T1, bool REFCNT1>
    friend class matrix;
    
    //struct Rep    // âíóòðåííåå ïðåäñòàâëåíèå
    //{
    //  // std::vector -- âðåìåííî, íóæíî ñâîé îáëåã÷¸ííûé âåêòîð
    //  typedef std::vector<T> Mem;
    //  typedef typename Mem::size_type size_type;
    //  size_type rows, cols;
    //  Mem mem;

    //  Rep () : nrows(0), ncols(0) {}

    //  void assign (size_type rows_a, size_type cols_a, const T& val)
    //  { rows = rows_a; cols = cols_a; mem.assign(rows*cols, val); }

    //  template <typename In>
    //  void assign (size_type rows_a, size_type cols_a, In b, In e)
    //  {
    //      rows = rows_a; cols = cols_a; mem.assign(b, e);
    //      ARAGELI_ASSERT_1(rows*cols == mem.size());
    //  }
    //};

    typedef _Internal::Rep_matrix<T> Rep;
    typedef typename Rep::Mem Mem;

public:


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    typedef typename Mem::value_type value_type;
    
    typedef typename Mem::value_type element_type;

    typedef typename Mem::reference reference;
    
    typedef typename Mem::const_reference const_reference;
    
    typedef typename Mem::pointer pointer;
    
    typedef typename Mem::const_pointer const_pointer;
    

    typedef typename Mem::size_type size_type;
    
    typedef typename Mem::difference_type difference_type;
    
    typedef typename Rep::Mem::iterator iterator;

    typedef typename Rep::Mem::const_iterator const_iterator;

    // @}


    template <typename T1> struct other_element_type
    { typedef matrix<T1, REFCNT> type; };

    template <typename T1, bool REFCNT1> struct other_element_type_refcnt
    { typedef matrix<T1, REFCNT1> type; };


    static const bool refcounting = REFCNT;


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    


    matrix () {}

    // BE CAREFUL! There is one more constructor. It is automaticaly generated
    // copy-constructor. Worry about its meaning is being belonged to other
    // constructors during modifications.


    // ------------------------------------------------------------------


    matrix (std::size_t size_a, const fromsize_t& select)
    { assign_fromsize(size_a); }


    matrix (const char* s, const fromstr_t& select)
    { assign_fromstr(s); }


    template <typename Vec>
    matrix (const Vec& vec, const fromvec_t& select)
    { assign_fromvec(vec); }


    template <typename Mat>
    matrix (const Mat& mat, const frommat_t& select)
    { assign_frommat(mat); }


    template <typename X>
    matrix (const X& x)
    { assign(x); }


    // ------------------------------------------------------------------



    template <typename Val>
    matrix (size_type rowscols_a, const Val& val, const fromval_t& select)
    { assign_fromval(rowscols_a, val); }



    template <typename Seq>
    matrix (size_type rowscols_a, const Seq& seq, const fromseq_t& select)
    { assign_fromseq(rowscols_a, seq); }



    matrix (size_type rows_a, size_type cols_a, const fromsize_t& select)
    { assign_fromsize(rows_a, cols_a); }


    template <typename Vec>
    matrix (size_type rows_a, const Vec& vec, const fromvec_t& select)
    { assign_fromvec(rows_a, vec); }


    template <typename X>
    matrix (size_type rowscols_a, const X& x)
    { assign(rowscols_a, x); }


    // ------------------------------------------------------------------


    template <typename Vec>
    matrix (const Vec& vec, const colwise_t& select)
    { assign_colwise(vec); }


    template <typename Vec>
    matrix
    (
        size_type cols_a, const Vec& vec,
        const colwise_t& select_1, const fromvec_t& select_2
    )
    { assign_colwise_fromvec(cols_a, vec); }


    template <typename Seq>
    matrix
    (
        size_type rowcols_a, const Seq& seq,
        const colwise_t& select_1, const fromseq_t& select_2
    )
    { assign_colwise_fromseq(rowcols_a, seq); }


    template <typename X>
    matrix (size_type rowcols_a, const X& x, const colwise_t& select)
    { assign_colwise(rowcols_a, x); }


    template <typename Seq>
    matrix
    (
        size_type rows_a, size_type cols_a,
        const Seq& seq, const colwise_t& select
    )
    { assign_colwise(rows_a, cols_a, seq); }


    // ------------------------------------------------------------------



    template <typename Val>
    matrix
    (
        size_type rows_a, size_type cols_a,
        const Val& val,
        const fromval_t& select
    )
    { assign_fromval(rows_a, cols_a, val); }



    template <typename Seq>
    matrix
    (
        size_type rows_a, size_type cols_a,
        const Seq& seq,
        const fromseq_t& select
    )
    { assign_fromseq(rows_a, cols_a, seq); }



    template <typename X>
    matrix (size_type rows_a, size_type cols_a, const X& x)
    { assign(rows_a, cols_a, x); }


    // ------------------------------------------------------------------



    template <typename Vec>
    matrix (const Vec& vec, const diag_t& select)
    { assign_diag(vec); }



    template <typename Val>
    matrix
    (
        size_type rowscols_a, const Val& val,
        const diag_t& select_1,
        const fromval_t& select_2
    )
    { assign_diag_fromval(rowscols_a, val); }



    template <typename Seq>
    matrix
    (
        size_type rowscols_a, const Seq& seq,
        const diag_t& select_1,
        const fromseq_t& select_2
    )
    { assign_diag_fromseq(rowscols_a, seq); }



    template <typename Vec>
    matrix
    (
        size_type rowscols_a, const Vec& vec,
        const diag_t& select_1,
        const fromvec_t& select_2
    )
    { assign_diag_fromvec(rowscols_a, vec); }


    template <typename X>
    matrix
    (
        size_type rowscols_a, const X& x,
        const diag_t& select
    )
    { assign_diag(rowscols_a, x); }



    template <typename Val>
    matrix
    (
        size_type rows_a, size_type cols_a,
        const Val& val,
        const diag_t& select_1,
        const fromval_t& select_2
    )
    { assign_diag_fromval(rows_a, cols_a, val); }



    template <typename Seq>
    matrix
    (
        size_type rows_a, size_type cols_a, const Seq& seq,
        const diag_t& select_1,
        const fromseq_t& select_2
    )
    { assign_diag_fromseq(rows_a, cols_a, seq); }



    template <typename Vec>
    matrix
    (
        size_type rows_a, size_type cols_a, const Vec& vec,
        const diag_t& select_1,
        const fromvec_t& select_2
    )
    { assign_diag_fromvec(rows_a, cols_a, vec); }



    template <typename X>
    matrix
    (
        size_type rows_a, size_type cols_a,
        const X& x,
        const diag_t& select
    )
    { assign_diag(rows_a, cols_a, x); }


    // ------------------------------------------------------------------



    matrix (size_type rowscols_a, const eye_t& select)
    { assign_eye(rowscols_a); }



    matrix
    (
        size_type rows_a, size_type cols_a,
        const eye_t& select
    )
    { assign_eye(rows_a, cols_a); }

    
    // @}




    //matrix (size_type rows_a, size_type cols_a, const T& val, const diag_t& select)
    //{ assign(rows_a, cols_a, val, diag); }

    //  square_matrix_t parameter if types size_type and T are correlated
    //  and overloading is not worked.
    //  @param rowscols_a size of a matrix, it is number of rows
    //      (or number of columns)
    //  @param val value for initialization of each item
    //  @param select helper argument for choosing this constructor,
    //      use square_matrix  */
    //matrix (size_type rowscols_a, T val, const square_matrix_t& select)
    //{ rep().assign(rowscols_a, rowscols_a, val); }


    //  @param cols_a number of columns
    //  @param val value for initialization of each items,
    //      by default this is TT::null() */
    //matrix (size_type rows_a, size_type cols_a, T val = factory<T>::null())
    //{ rep().assign(rows_a, cols_a, val); }


    //  nonsquare_t parameter if types size_type and T are correlated
    //  and overloading is not worked.
    //  @param rows_a number of rows
    //  @param cols_a number of columns
    //  @param select helper argument for choosing this constructor,
    //      use fromsize  */
    //matrix (size_type rows_a, size_type cols_a, const nonsquare_t& select)
    //{ rep().assign(rows_a, cols_a, factory<T>::null()); }


    //      and number of columns
    //  @param in   an iterator to begin of the sequence, and elements of matrix
    //      will be columnwise constructed from this sequence. */
    //template <typename In>
    //matrix
    //(
    //  size_type rowscols_a, const In& in,
    //  const colwise_t& select
    //)
    //{ assign(rowscols_a, in, select); }


    //  @param cols_a   number of columns
    //  @param in   an iterator to begin of the sequence, and elements of matrix
    //      will be rowwise constructed from this sequence */
    //template <typename In>
    //matrix
    //(
    //  size_type rows_a, size_type cols_a,
    //  const In& in,
    //  const colwise_t& select
    //)
    //{ assign(rows_a, cols_a, in, select); }

    //  @param cols_a   number of columns */
    //matrix (size_type rows_a, size_type cols_a, const nonsquare_t&)
    //{ rep().assign(rows_a, cols_a, factory<T>::null()); }


    //  you want to call this version of a constructor, you can call more specific
    //  form of this: matrix (rows_a, cols_a, nonsquare_t).
    //  @param rows_a   number of rows
    //  @param cols_a   number of columns */
    //matrix (size_type rows_a, size_type cols_a)
    //{ rep().assign(rows_a, cols_a, factory<T>::null()); }


    //#if 0
    //
    //  /// Initializes matrix by its string notation.
    //  /** Input rules are the same as for operator>> with default behaviour. */
    //  matrix (const char* s);
    //
    //
    //  /// Copy constructor from matrix with any other type of template parameters.
    //  template <typename T1, bool REFCNT1>
    //  matrix (const matrix<T1, REFCNT1>& x)
    //  {
    //      ARAGELI_ASSERT_1
    //      (
    //          reinterpret_cast<const void*>(this) !=
    //          reinterpret_cast<const void*>(&x)
    //      );
    //
    //      rep().assign(x.nrows(), x.ncols(), x.mem().begin(), x.mem().end());
    //  }
    //
    //  /// Creates square matrix, each item is equal to factory<T>::null().
    //  /** @param rowscols_a size of a matrix, it is number of rows
    //          and number of columns */
    //  explicit matrix (size_type rowscols_a)
    //  { rep().assign(rowscols_a, rowscols_a, factory<T>::null()); }
    //
    //#endif


    
    



    void assign_fromsize (std::size_t size_a)
    { assign_fromsize(size_a, size_a); }


    void assign_fromstr (const char* s);


    template <typename Vec>
    void assign_fromvec (const Vec& vec)
    { assign_fromseq(1u, vec.size(), vec.begin()); }
    

    template <typename Mat>
    void assign_frommat (const Mat& mat)
    { assign_fromseq(mat.nrows(), mat.ncols(), mat.begin()); }


    // ------------------------------------------------------------------



    template <typename Val>
    void assign_fromval (size_type rowscols_a, const Val& val)
    { assign_fromval(rowscols_a, rowscols_a, val); }



    template <typename Seq>
    void assign_fromseq (size_type rowscols_a, const Seq& seq)
    { assign_fromseq(rowscols_a, rowscols_a, seq); }



    void assign_fromsize (size_type rows_a, size_type cols_a)
    { assign_fromval(rows_a, cols_a, factory<T>::null()); }


    template <typename Vec>
    void assign_fromvec (size_type rows_a, const Vec& vec); // It isn't implemented yet.


    // ------------------------------------------------------------------


    template <typename Vec>
    void assign_colwise (const Vec& vec);   // It isn't implemented yet.


    template <typename Vec>
    void assign_colwise_fromvec (size_type cols_a, const Vec& vec); // It isn't implemented yet.


    template <typename Seq>
    void assign_colwise_fromseq (size_type rowcols_a, const Seq& seq);  // It isn't implemented yet.


    template <typename X>
    void assign_colwise (size_type rowcols_a, const X& x);  // It isn't implemented yet.


    template <typename Seq>
    void assign_colwise
    (size_type rows_a, size_type cols_a, const Seq& seq);   // It isn't implemented yet.


    // ------------------------------------------------------------------



    template <typename Val>
    void assign_fromval
    (size_type rows_a, size_type cols_a, const Val& val);



    template <typename Seq>
    void assign_fromseq
    (size_type rows_a, size_type cols_a, Seq seq);


    // ------------------------------------------------------------------



    template <typename Vec>
    void assign_diag (const Vec& vec)
    { assign_diag_fromseq(vec.size(), vec.begin()); }



    template <typename Val>
    void assign_diag_fromval (size_type rowscols_a, const Val& val)
    { assign_diag_fromval(rowscols_a, rowscols_a, val); }



    template <typename Seq>
    void assign_diag_fromseq (size_type rowscols_a, const Seq& seq)
    { assign_diag_fromseq(rowscols_a, rowscols_a, seq); }



    template <typename Vec>
    void assign_diag_fromvec (size_type rowscols_a, const Vec& vec)
    { assign_diag_fromseq(rowscols_a, vec.begin()); }


    template <typename X>
    void assign_diag (size_type rowscols_a, const X& x)
    { assign_diag(rowscols_a, rowscols_a, x); }



    template <typename Val>
    void assign_diag_fromval
    (size_type rows_a, size_type cols_a, const Val& val);



    template <typename Seq>
    void assign_diag_fromseq
    (size_type rows_a, size_type cols_a, Seq seq);



    template <typename Vec>
    void assign_diag_fromvec
    (size_type rows_a, size_type cols_a, const Vec& vec)
    { assign_diag_fromseq(rows_a, cols_a, vec.begin()); }



    template <typename X>
    void assign_diag
    (size_type rows_a, size_type cols_a, const X& x)
    { _assign_diag_val_seq_vec(rows_a, cols_a, x, type_traits<X>::category_value); }

private:

    template <typename Val>
    void _assign_diag_val_seq_vec
    (size_type rows_a, size_type cols_a, const Val& val, const type_category::type&)
    { assign_diag_fromval(rows_a, cols_a, val); }

    template <typename Seq>
    void _assign_diag_val_seq_vec
    (size_type rows_a, size_type cols_a, const Seq& seq, const type_category::iterator&)
    { assign_diag_fromseq(rows_a, cols_a, seq); }

    template <typename Vec>
    void _assign_diag_val_seq_vec
    (size_type rows_a, size_type cols_a, const Vec& vec, const type_category::vector&)
    { assign_diag_fromvec(rows_a, cols_a, vec); }

public:

    // ------------------------------------------------------------------



    void assign_eye (size_type rowscols_a)
    { assign_eye(rowscols_a, rowscols_a); }



    void assign_eye (size_type rows_a, size_type cols_a)
    { assign_diag_fromval(rows_a, cols_a, factory<T>::unit()); }
    
    


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    // @{


    void assign ()
    { resize(0, 0); }


    // ------------------------------------------------------------------


    //void assign (std::size_t size_a, const fromsize_t& select)
    //{ assign_fromsize(size_a); }


    //void assign (const char* s, const fromstr_t& select)
    //{ assign_fromstr(s); }


    //template <typename Vec>
    //void assign (const Vec& vec, const fromvec_t& select)
    //{ assign_fromvec(vec); }
    //

    //template <typename Mat>
    //void assign (const Mat& mat, const frommat_t& select)
    //{ assign_frommat(mat); }


    template <typename X>
    void assign (const X& x)
    { _assign_size_str_vec_mat(x, type_traits<X>::category_value); }

private:

    template <typename Size>
    void _assign_size_str_vec_mat
    (const Size& size_a, const type_category::integer&)
    { assign_fromsize(size_a); }

    template <typename Str>
    void _assign_size_str_vec_mat
    (const Str& s, const type_category::string&)
    { assign_fromstr(s); }

    template <typename Vec>
    void _assign_size_str_vec_mat
    (const Vec& vec, const type_category::vector&)
    { assign_fromvec(vec); }

    template <typename Mat>
    void _assign_size_str_vec_mat
    (const Mat& mat, const type_category::matrix&)
    { assign_frommat(mat); }

public:

    // ------------------------------------------------------------------


    //      and number of columns
    //  @param val an initial value for each element */
    //template <typename Val>
    //void assign (size_type rowscols_a, const Val& val, const fromval_t& select)
    //{ assign_fromval(rowscols_a, val); }


    //      and number of columns
    //  @param seq  an iterator to begin of the sequence, and elements of matrix
    //      will be rowwise constructed from this sequence. */
    //template <typename Seq>
    //void assign (size_type rowscols_a, const Seq& seq, const fromseq_t& select)
    //{ assign_fromseq(rowscols_a, seq); }


    //  @param cols_a   number of columns */
    //void assign (size_type rows_a, size_type cols_a, const fromsize_t& select)
    //{ assign_fromsize(rows_a, cols_a); }


    //template <typename Vec>
    //void assign (size_type rows_a, const Vec& vec, const fromvec_t& select)
    //{ assign_fromvec(rows_a, vec); }


    template <typename X>
    void assign (size_type rowscols_a, const X& x)
    {
        _assign_seq_size_vec
        (
            rowscols_a, x,
            type_traits<X>::category_value//,
            //bool_type<type_pair_traits<X, T>::is_assignable>::value
        );
    }


private:

    //template <typename Val>
    //void _assign_val_seq_size_vec
    //(
    //  size_type rowscols_a, const Val& val,
    //  const type_category::type&, const true_type&
    //)
    //{ assign_fromval(rowscols_a, val); }

    template <typename Seq>
    void _assign_seq_size_vec
    (
        size_type rowscols_a, const Seq& seq,
        const type_category::iterator&//, const false_type&
    )
    { assign_fromseq(rowscols_a, seq); }

    //template <typename Size>
    //void _assign_val_seq_size_vec
    //(
    //  size_type rows_a, const Size& cols_a,
    //  const type_category::integer&//, const true_type&
    //)
    //{ assign_fromsize(rows_a, cols_a); }

    template <typename Size>
    void _assign_seq_size_vec
    (
        size_type rows_a, const Size& cols_a,
        const type_category::integer&//, const false_type&
    )
    { assign_fromsize(rows_a, cols_a); }

    template <typename Vec>
    void _assign_seq_size_vec
    (
        size_type rows_a, const Vec& vec,
        const type_category::vector&//, const false_type&
    )
    { assign_fromvec(rows_a, vec); }

public:

    // ------------------------------------------------------------------


    //template <typename Vec>
    //void assign (const Vec& vec, const colwise_t& select)
    //{ assign_colwise(vec); }


    //template <typename Vec>
    //void assign
    //(
    //  size_type cols_a, const Vec& vec,
    //  const colwise_t& select_1, const fromvec_t& select_2
    //)
    //{ assign_colwise_fromvec(cols_a, vec); }


    //template <typename Seq>
    //void assign
    //(
    //  size_type rowcols_a, const Seq& seq,
    //  const colwise_t& select_1, const fromseq_t& select_2
    //)
    //{ assign_colwise_fromseq(rowcols_a, seq); }


    //template <typename Seq>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const Seq& seq, const colwise_t& select
    //)
    //{ assign_colwise(rows_a, cols_a, seq); }


    // ------------------------------------------------------------------


    //  @param cols_a   number of columns
    //  @param val  an initial value for each element */
    //template <typename Val>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const Val& val,
    //  const fromval_t& select
    //)
    //{ assign_fromval(rows_a, cols_a, val); }


    //  @param cols_a   number of columns
    //  @param in   an iterator to begin of the sequence, and elements of matrix
    //      will be rowwise constructed from this sequence */
    //template <typename Seq>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const Seq& seq,
    //  const fromseq_t& select
    //)
    //{ assign_fromseq(rows_a, cols_a, seq); }



    template <typename X>
    void assign (size_type rows_a, size_type cols_a, const X& x)
    {
        _assign_val_seq
        (
            rows_a, cols_a, x,
            bool_type<type_pair_traits<X, T>::is_assignable>::value
        );
    }


private:

    template <typename Val>
    void _assign_val_seq
    (
        size_type rows_a, size_type cols_a, const Val& val,
        const true_type&
    )
    { assign_fromval(rows_a, cols_a, val); }

    template <typename Seq>
    void _assign_val_seq
    (
        size_type rows_a, size_type cols_a, const Seq& seq,
        const false_type&
    )
    { assign_fromseq(rows_a, cols_a, seq); }

public:

    // ------------------------------------------------------------------


    //  @param seq sequence of init values */
    //template <typename Vec>
    //void assign (const Vec& vec, const diag_t& select)
    //{ assign_diag(vec); }


    //  @param val value on the diagonal */
    //template <typename Val>
    //void assign
    //(
    //  size_type rowscols_a, const Val& val,
    //  const diag_t& select_1,
    //  const fromval_t& select_2
    //)
    //{ assign_diag_fromval(rowscols_a, val); }


    //  @param seq sequence of init values */
    //template <typename Seq>
    //void assign
    //(
    //  size_type rowscols_a, const Seq& seq,
    //  const diag_t& select_1,
    //  const fromseq_t& select_2
    //)
    //{ assign_diag_fromseq(rowscols_a, seq); }


    //  @param seq sequence of init values */
    //template <typename Vec>
    //void assign
    //(
    //  size_type rowscols_a, const Vec& vec,
    //  const diag_t& select_1,
    //  const fromvec_t& select_2
    //)
    //{ assign_diag_fromvec(rowscols_a, vec); }


    //template <typename X>
    //void assign
    //(
    //  size_type rowscols_a, const X& x,
    //  const diag_t& select
    //)
    //{ assign_diag(rowscols_a, x); }


    //  @param cols_a   number of columns
    //  @param val value on the diagonal */
    //template <typename Val>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const Val& val,
    //  const diag_t& select_1,
    //  const fromval_t& select_2
    //)
    //{ assign_diag_fromval(rows_a, cols_a, val); }


    //  @param cols_a   number of columns
    //  @param seq sequence of init values */
    //template <typename Seq>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a, const Seq& seq,
    //  const diag_t& select_1,
    //  const fromseq_t& select_2
    //)
    //{ assign_diag_fromseq(rows_a, cols_a, seq); }


    //  @param cols_a   number of columns
    //  @param vec sequence of init values */
    //template <typename Vec>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a, const Vec& vec,
    //  const diag_t& select_1,
    //  const fromvec_t& select_2
    //)
    //{ assign_diag_fromvec(rows_a, cols_a, vec); }


    //  @param cols_a   number of columns
    //  @param x    Meaning of this parameter depends on a structure of X. */
    //template <typename X>
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const X& x,
    //  const diag_t& select
    //)
    //{ assign_diag(rows_a, cols_a, x); }




    //  @param select helper argument for choosing this constructor,
    //      use eye  */
    //void assign (size_type rowscols_a, const eye_t& select)
    //{ assign_eye(rowscols_a); }


    //  @param select helper argument for choosing this constructor,
    //      use eye  */
    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const eye_t& select
    //)
    //{ assign_eye(rows_a, cols_a); }


    // ------------------------------------------------------------------

    template <typename X>
    matrix& operator= (const X& x)
    {
        _assign_from_val_str_vec_mat
        (
            x, type_traits<X>::category_value,
            bool_type<type_pair_traits<X, T>::is_assignable>::value
        );
        
        return *this;
    }

private:

    template <typename Val>
    void _assign_from_val_str_vec_mat
    (const Val& val, const type_category::type&, const true_type&)
    { fill_fromval(val); }

    template <typename Str>
    void _assign_from_val_str_vec_mat
    (const Str& str, const type_category::string&, const false_type&)
    { assign_fromstr(str); }

    template <typename Vec>
    void _assign_from_val_str_vec_mat
    (const Vec& vec, const type_category::vector&, const false_type&)
    { assign_fromvec(vec); }

    template <typename Mat>
    void _assign_from_val_str_vec_mat
    (const Mat& mat, const type_category::matrix&, const false_type&)
    { assign_frommat(mat); }

public:


#if 0


    template <typename T1, bool REFCNT1>
    matrix& operator= (const matrix<T1, REFCNT1>& x)
    {
        ARAGELI_ASSERT_1
        (
            reinterpret_cast<const void*>(this) !=
            reinterpret_cast<const void*>(&x)
        );

        unique_clear();
        rep().assign(x.nrows(), x.ncols(), x.mem().begin(), x.mem().end());
        return *this;
    }

    template <typename T1>
    matrix& operator= (const T1& v)
    {
        unique_clear();
        std::fill(mem().begin(), mem().end(), v);
        return *this;
    }

    matrix& operator= (const char* s)
    { return *this = matrix(s); }


    void assign ()
    {
        unique_clear();
        mem().clear();
        rep().rows = 0;
        rep().cols = 0;
    }



    void assign (size_type rowscols_a)
    {
        unique_clear();
        rep().assign(rowscols_a, rowscols_a, factory<T>::null());
    }



    void assign (size_type rows_a, size_type cols_a)
    {
        unique_clear();
        rep().assign(rows_a, cols_a, factory<T>::null());
    }
    


    void assign
    (
        size_type rows_a, size_type cols_a,
        const nonsquare_t& select
    )
    {
        unique_clear();
        rep().assign(rows_a, cols_a, factory<T>::null());
    }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    template <typename T1>
    void assign (size_type rowscols_a, const T1& val, const fromval_t& select)
    { assign_fromval(rowscols_a, val); }



    template <typename In>
    void assign (size_type rowscols_a, const In& in, const fromseq_t& select)
    { assign_fromseq(rowscols_a, in); }



    template <typename T1>
    void assign_fromval (size_type rowscols_a, const T1& val)
    {
        unique_clear();
        rep().assign(rowscols_a, rowscols_a, val);
    }



    template <typename In>
    void assign_fromseq (size_type rowscols_a, const In& in)
    {
        unique_clear();
        In last = in;
        std::advance(last, rowscols_a*rowscols_a);  // WARNING! Temporary solution.
        rep().assign(rowscols_a, rowscols_a, in, last);
    }



    template <typename In_or_T1>
    void assign (size_type rowscols_a, const In_or_T1& x)
    {
        assign_from_seq_or_val
        (
            rowscols_a, x,
            bool_type<type_pair_traits<In_or_T1, T>::is_convertible>::value
        );
    }


    template <typename In_or_T1>
    void assign_from_seq_or_val
    (size_type rowscols_a, const In_or_T1& x, const true_type&)
    { assign_fromval(rowscols_a, x); }


    template <typename In_or_T1>
    void assign_from_seq_or_val
    (size_type rowscols_a, const In_or_T1& x, const false_type&)
    { assign_fromseq(rowscols_a, x); }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    template <typename T1>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const T1& val, const fromval_t& select
    )
    { assign_fromval(rows_a, cols_a, val); }



    template <typename In>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const In& in, const fromseq_t& select
    )
    { assign_fromseq(rows_a, cols_a, in); }



    template <typename T1>
    void assign_fromval (size_type rows_a, size_type cols_a, const T1& val)
    {
        unique_clear();
        rep().assign(rows_a, cols_a, val);
    }



    template <typename In>
    void assign_fromseq (size_type rows_a, size_type cols_a, const In& in)
    {
        unique_clear();
        In last = in;
        std::advance(last, rows_a*cols_a);  // WARNING! Temporary solution.
        rep().assign(rows_a, cols_a, in, last);
    }



    template <typename In_or_T1>
    void assign (size_type rows_a, size_type cols_a, const In_or_T1& x)
    {
        assign_from_seq_or_val
        (
            rows_a, cols_a, x,
            bool_type<type_pair_traits<In_or_T1, T>::is_convertible>::value
        );
    }


    template <typename In_or_T1>
    void assign_from_seq_or_val
    (size_type rows_a, size_type cols_a, const In_or_T1& x, const true_type&)
    { assign_fromval(rows_a, cols_a, x); }


    template <typename In_or_T1>
    void assign_from_seq_or_val
    (size_type rows_a, size_type cols_a, const In_or_T1& x, const false_type&)
    { assign_fromseq(rows_a, cols_a, x); }

    
    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    template <typename In>
    void assign
    (
        size_type rowscols_a, const In& in,
        const colwise_t& select
    )
    { assign_colwise(rowscols_a, in); }



    template <typename In>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const In& in,
        const colwise_t& select
    )
    { assign_colwise(rows_a, cols_a, in); }



    template <typename In>
    void assign_colwise (size_type rowscols_a, const In& in)
    { assign_colwise(rowscols_a, rowscols_a, in); }



    template <typename In>
    void assign_colwise
    (size_type rows_a, size_type cols_a, In in);


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    template <typename T1>
    void assign
    (
        size_type rowscols_a, const T1& val,
        const diag_t& select_1,
        const fromval_t& select_2
    )
    { assign_diag_fromval(rowscols_a, val); }



    template <typename In>
    void assign
    (
        size_type rowscols_a, const In& in,
        const diag_t& select_1,
        const fromseq_t& select_2
    )
    { assign_diag_fromseq(rowscols_a, in); }



    template <typename T1>
    void assign_diag_fromval (size_type rowscols_a, const T1& val)
    { assign_diag_fromval(rowscols_a, rowscols_a, val); }



    template <typename In>
    void assign_diag_fromseq (size_type rowscols_a, const In& in)
    { assign_diag_fromseq(rowscols_a, rowscols_a, in); }



    template <typename In_or_T1>
    void assign_diag (size_type rowscols_a, const In_or_T1& x)
    {
        assign_diag_from_seq_or_val
        (
            rowscols_a, x,
            bool_type<type_pair_traits<In_or_T1, T>::is_convertible>::value
        );
    }



    template <typename In_or_T1>
    void assign
    (
        size_type rowscols_a, const In_or_T1& x,
        const diag_t& select_1
    )
    { assign_diag(rowscols_a, x); }


    template <typename In_or_T1>
    void assign_diag_from_seq_or_val
    (size_type rowscols_a, const In_or_T1& x, const true_type&)
    { assign_diag_fromval(rowscols_a, x); }


    template <typename In_or_T1>
    void assign_diag_from_seq_or_val
    (size_type rowscols_a, const In_or_T1& x, const false_type&)
    { assign_diag_fromseq(rowscols_a, x); }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    template <typename T1>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const T1& val,
        const diag_t& select_1,
        const fromval_t& select_2
    )
    { assign_diag_fromval(rows_a, cols_a, val); }



    template <typename In>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const In& in,
        const diag_t& select_1,
        const fromseq_t& select_2
    )
    { assign_diag_fromseq(rows_a, cols_a, in); }



    template <typename T1>
    void assign_diag_fromval (size_type rows_a, size_type cols_a, const T1& val);



    template <typename In>
    void assign_diag_fromseq (size_type rows_a, size_type cols_a, In in);



    template <typename In_or_T1>
    void assign_diag
    (
        size_type rows_a, size_type cols_a,
        const In_or_T1& x
    )
    {
        assign_diag_from_seq_or_val
        (
            rows_a, cols_a, x,
            bool_type<type_pair_traits<In_or_T1, T>::is_convertible>::value
        );
    }



    template <typename In_or_T1>
    void assign
    (
        size_type rows_a, size_type cols_a,
        const In_or_T1& x,
        const diag_t& select_1
    )
    { assign_diag(rows_a, cols_a, x); }


    template <typename In_or_T1>
    void assign_diag_from_seq_or_val
    (
        size_type rows_a, size_type cols_a,
        const In_or_T1& x, const true_type&
    )
    { assign_diag_fromval(rows_a, cols_a, x); }


    template <typename In_or_T1>
    void assign_diag_from_seq_or_val
    (
        size_type rows_a, size_type cols_a,
        const In_or_T1& x, const false_type&
    )
    { assign_diag_fromseq(rows_a, cols_a, x); }




    //  square_matrix_t parameter if types size_type and T are correlated
    //  and overloading is not worked.
    //  @param rowscols_a size of a matrix, it is number of rows
    //      (or number of columns)
    //  @param val value for initialization of each item
    //  @param select helper argument for choosing this constructor,
    //      use square_matrix  */
    //void assign (size_type rowscols_a, T val, const square_matrix_t& select)
    //{
    //  unique_clear();
    //  rep().assign(rowscols_a, rowscols_a, val);
    //}


    //  @param cols_a number of columns
    //  @param val value for initialization of each items,
    //      by default this is TT::null() */
    //void assign (size_type rows_a, size_type cols_a, T val = factory<T>::null())
    //{
    //  unique_clear();
    //  rep().assign(rows_a, cols_a, val);
    //}




    //  @param val value on the diagonal
    //  @param select helper argument for choosing this constructor,
    //      use diag  */
    //void assign (size_type rowscols_a, const T val, const diag_t& select);

    //void assign
    //(
    //  size_type rows_a, size_type cols_a,
    //  const T val, const diag_t& select
    //);


    void assign
    (
        size_type rowscols_a,
        const eye_t& select
    )
    {
        unique_clear();
        assign_diag_fromval(rowscols_a, factory<T>::unit());
    }


    void assign
    (
        size_type rows_a, size_type cols_a,
        const eye_t& select
    )
    {
        unique_clear();
        assign_diag_fromval(rows_a, cols_a, factory<T>::unit());
    }

#endif

    // @}


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    // @{


    const_reference el (size_type i, size_type j) const
    {
        ARAGELI_ASSERT_0(i < nrows() && j < ncols());
        return mem(i*ncols() + j);
    }
    

    reference el (size_type i, size_type j)
    {
        ARAGELI_ASSERT_0(i < nrows() && j < ncols());
        unique();
        return mem(i*ncols() + j);
    }


    const_reference operator() (size_type i, size_type j) const
    { return el(i, j); }
    

    reference operator() (size_type i, size_type j)
    { return el(i, j); }


    const_reference at (size_type i, size_type j) const
    {
        if(i > nrows() || j > ncols())throw out_of_range();
        return el(i, j);
    }
    

    reference at (size_type i, size_type j)
    {
        if(i > nrows() || j > ncols())throw out_of_range();
        return el(i, j);
    }

    // @}


    void resize (size_type rows_a, size_type cols_a)
    {
        unique_clear();
        mem().resize(rows_a*cols_a);
        rep().rows = rows_a;
        rep().cols = cols_a;
    }


    bool pack ()
    {
        if(size() != capacity())
        {
            do_pack();
            return true;
        }
        else return false;
    }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // Some properties access.


    size_type nrows () const { return rep().rows; }

    size_type ncols () const { return rep().cols; }
    
    size_type size () const { return mem().size(); }


    size_type length () const { return std::max(nrows(), ncols()); }

    size_type capacity () const { return mem().capacity(); }


    bool is_empty () const
    {
        ARAGELI_ASSERT_1(!mem().empty() || (nrows() == 0 || ncols() == 0));
        return mem().empty();
    }


    bool is_null () const;


    bool is_unit () const;


    bool is_opposite_unit () const;


    bool is_square () const { return nrows() == ncols(); }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // Some modificatory operations.



    matrix& opposite ();


    matrix& inverse ();


    matrix& each_inverse ();


    template <typename T1> matrix& add_scalar (const T1& x);
    template <typename T1> matrix& sub_scalar (const T1& x);
    template <typename T1> matrix& mul_scalar (const T1& x);
    template <typename T1> matrix& div_scalar (const T1& x);
    template <typename T1> matrix& mod_scalar (const T1& x);
    template <typename T1> matrix& bitand_scalar (const T1& x);
    template <typename T1> matrix& bitor_scalar (const T1& x);
    template <typename T1> matrix& bitxor_scalar (const T1& x);
    template <typename T1> matrix& shl_scalar (const T1& x);
    template <typename T1> matrix& shr_scalar (const T1& x);
    
    template <typename T1> matrix& add_matrix (const T1& x) { return each_add_matrix(x); }
    template <typename T1> matrix& sub_matrix (const T1& x) { return each_sub_matrix(x); }
    template <typename T1> matrix& mul_matrix (const T1& x);
    template <typename T1> matrix& div_matrix (const T1& x);
    template <typename T1> matrix& mod_matrix (const T1& x) { return each_mod_matrix(x); }
    template <typename T1> matrix& bitand_matrix (const T1& x) { return each_bitand_matrix(x); }
    template <typename T1> matrix& bitor_matrix (const T1& x) { return each_bitor_matrix(x); }
    template <typename T1> matrix& bitxor_matrix (const T1& x) { return each_bitxor_matrix(x); }
    template <typename T1> matrix& shl_matrix (const T1& x) { return each_shl_matrix(x); }
    template <typename T1> matrix& shr_matrix (const T1& x) { return each_shr_matrix(x); }

    template <typename T1> matrix& each_add_matrix (const T1& x);
    template <typename T1> matrix& each_sub_matrix (const T1& x);
    template <typename T1> matrix& each_mul_matrix (const T1& x);
    template <typename T1> matrix& each_div_matrix (const T1& x);
    template <typename T1> matrix& each_mod_matrix (const T1& x);
    template <typename T1> matrix& each_bitand_matrix (const T1& x);
    template <typename T1> matrix& each_bitor_matrix (const T1& x);
    template <typename T1> matrix& each_bitxor_matrix (const T1& x);
    template <typename T1> matrix& each_shl_matrix (const T1& x);
    template <typename T1> matrix& each_shr_matrix (const T1& x);



    template <typename T1> matrix& operator+= (const T1& x) { return add_scalar(x); }
    template <typename T1> matrix& operator-= (const T1& x) { return sub_scalar(x); }
    template <typename T1> matrix& operator*= (const T1& x) { return mul_scalar(x); }
    template <typename T1> matrix& operator/= (const T1& x) { return div_scalar(x); }
    template <typename T1> matrix& operator%= (const T1& x) { return mod_scalar(x); }
    template <typename T1> matrix& operator&= (const T1& x) { return bitand_scalar(x); }
    template <typename T1> matrix& operator|= (const T1& x) { return bitor_scalar(x); }
    template <typename T1> matrix& operator^= (const T1& x) { return bitxor_scalar(x); }
    template <typename T1> matrix& operator<<= (const T1& x) { return shl_scalar(x); }
    template <typename T1> matrix& operator>>= (const T1& x) { return shr_scalar(x); }

    // BEGIN OF TEMPORARY EXTENTION
    
    matrix& operator+= (const T& x) { return add_scalar(x); }
    matrix& operator-= (const T& x) { return sub_scalar(x); }
    matrix& operator*= (const T& x) { return mul_scalar(x); }
    matrix& operator/= (const T& x) { return div_scalar(x); }
    matrix& operator%= (const T& x) { return mod_scalar(x); }
    matrix& operator&= (const T& x) { return bitand_scalar(x); }
    matrix& operator|= (const T& x) { return bitor_scalar(x); }
    matrix& operator^= (const T& x) { return bitxor_scalar(x); }
    matrix& operator<<= (const T& x) { return shl_scalar(x); }
    matrix& operator>>= (const T& x) { return shr_scalar(x); }

    // END OF TEMPORARY EXTENTION


    // @}



    template <typename T1, bool REFCNT1>
    matrix& operator+= (const matrix<T1, REFCNT1>& x) { return add_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& operator-= (const matrix<T1, REFCNT1>& x) { return sub_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_add (const matrix<T1, REFCNT1>& x) { return each_add_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_sub (const matrix<T1, REFCNT1>& x) { return each_sub_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_mul (const matrix<T1, REFCNT1>& x) { return each_mul_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_div (const matrix<T1, REFCNT1>& x) { return each_div_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_mod (const matrix<T1, REFCNT1>& x) { return each_mod_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_bitand (const matrix<T1, REFCNT1>& x) { return each_bitand_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_bitor (const matrix<T1, REFCNT1>& x) { return each_bitor_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_bitxor (const matrix<T1, REFCNT1>& x) { return each_bitxor_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_shl (const matrix<T1, REFCNT1>& x) { return each_shl_matrix(x); }

    template <typename T1, bool REFCNT1>
    matrix& each_shr (const matrix<T1, REFCNT1>& x) { return each_shr_matrix(x); }

    // @}


    template <typename T1, bool REFCNT1>
    matrix& operator*= (const matrix<T1, REFCNT1>& x) { return mul_matrix(x); }


    template <typename T1, bool REFCNT1>
    matrix& operator/= (const matrix<T1, REFCNT1>& x) { return div_matrix(x); }


    matrix operator- () const
    {
        matrix res(*this);
        return res.opposite();
    }
    
    const matrix& operator+ () const { return *this; }

    matrix& operator++ ()
    {
        std::for_each(begin(), end(), func::prefix_increment<T, T&>());
        return *this;
    }

    matrix operator++ (int) { matrix t = *this; operator++(); return t; }

    matrix& operator-- ()
    {
        std::for_each(begin(), end(), func::prefix_decrement<T, T&>());
        return *this;
    }

    matrix operator-- (int) { matrix t = *this; operator--(); return t; }


    //std::pair<size_type, size_type> find_max () const;
    //std::pair<size_type, size_type> find_min () const;

    void transpose ();

    template <typename T1, bool REFCNT1>
    void swap (matrix<T1, REFCNT1>& x)
    {
        _Internal::swap_help_1
            (*this, x, store, x.store, equal_types<T, T1>::value);
    }


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    
    // @{


    iterator begin () { unique(); return mem().begin(); }
    

    iterator end () { unique(); return mem().end(); }

    const_iterator begin () const { return mem().begin(); }

    const_iterator end () const { return mem().end(); }

    // @}


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    // @{


    void insert_row (size_type pos, const T& val);
    

    template <typename T1, bool REFCNT1>
    void insert_row (size_type pos, const vector<T1, REFCNT1>& vals);


    template <typename In>
    void insert_row (size_type pos, In first);


    void insert_rows (size_type pos, size_type n, const T& val);
    

    template <typename T1, bool REFCNT1>
    void insert_rows (size_type pos, size_type n, const vector<T1, REFCNT1>& vals);


    template <typename In>
    void insert_rows (size_type pos, size_type n, In first);

    void erase_row (size_type pos);

    void erase_rows (size_type pos, size_type n);

    void swap_rows (size_type xpos, size_type ypos);

    template <typename T1>
    void mult_row (size_type i, const T1& x);
    
    template <typename T1>
    void div_row (size_type i, const T1& x);
    
    template <typename T1>
    void mod_row (size_type i, const T1& x);
    
    void add_rows (size_type dstpos, size_type srcpos);
    
    void sub_rows (size_type dstpos, size_type srcpos);

    template <typename T2>
    void addmult_rows (size_type dstpos, size_type srcpos, const T2& y);

    template <typename T1, typename T2>
    void addmult_rows
    (
        size_type dstpos, const T1& x,
        size_type srcpos, const T2& y
    );

    // WARNING! TEMPORARY DECLARATION. Conflicts are possible.
    template <typename T1, typename T2>
    void addmult_rows
    (
        size_type dstpos, const T1& x,
        int srcpos, const T2& y
    )
    { addmult_rows(dstpos, x, size_type(srcpos), y); }

    // WARNING! TEMPORARY DECLARATION. Conflicts are possible.
    template <typename T1, typename V, typename T2>
    void addmult_rows
    (
        size_type dstpos, const T1& x,
        const V& srcvec, const T2& y
    );

    vector<T, true> copy_row (size_type i) const
    //{ return vector<T, true>(begin() + i*ncols(), begin() + (i + 1)*ncols()); }
    { return vector<T, true>(ncols(), begin() + i*ncols()); }


    template <typename V>
    V& copy_row (size_type i, V& res) const
    {
        //res.assign(begin() + i*ncols(), begin() + (i + 1)*ncols());
        res.assign(ncols(), begin() + i*ncols());
        return res;
    }

    template <typename SV>
    matrix<T, true> copy_rows (const SV& sv) const
    {
        matrix<T, true> res;
        return copy_rows(sv, res);
    }

    template <typename SV, typename M>
    M& copy_rows (const SV& sv, M& res) const;
    
    template <typename SV>
    void erase_rows (const SV& sv);

    template <typename SV, typename M>
    M& take_rows (const SV& sv, M& res)
    {
        copy_rows(sv, res);
        erase_rows(sv);
        return res;
    }

    template <typename SV>
    matrix<T, true> take_rows (const SV& sv)
    {
        matrix<T, true> res;
        return take_rows(sv, res);
    }


    template <typename V>
    V& take_row (size_type i, V& res)
    {
        copy_row(i, res);
        erase_row(i);
        return res;
    }

    vector<T, true> take_row (size_type i)
    {
        vector<T, true> res;
        return take_row(i, res);
    }


    template <typename V>
    void assign_row (size_type i, const V& v)
    {
        ARAGELI_ASSERT_0(i < nrows());
        ARAGELI_ASSERT_0(v.size() >= ncols());
        std::copy(v.begin(), v.begin() + ncols(), begin() + ncols()*i);
    }


    // @}


    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    
    // @{


    void insert_col (size_type pos, const T& val);
    

    template <typename T1, bool REFCNT1>
    void insert_col (size_type pos, const vector<T1, REFCNT1>& vals);


    template <typename In>
    void insert_col (size_type pos, In first);


    void insert_cols (size_type pos, size_type n, const T& val);
    

    template <typename T1, bool REFCNT1>
    void insert_cols (size_type pos, size_type n, const vector<T1, REFCNT1>& vals);


    template <typename In>
    void insert_cols (size_type pos, size_type n, In first);

    void erase_col (size_type pos);

    void erase_cols (size_type pos, size_type n);

    void swap_cols (size_type xpos, size_type ypos);

    template <typename T1>
    void mult_col (size_type i, const T1& x);
    
    template <typename T1>
    void div_col (size_type i, const T1& x);
    
    template <typename T1>
    void mod_col (size_type i, const T1& x);
    
    void add_cols (size_type dstpos, size_type srcpos);
    
    void sub_cols (size_type dstpos, size_type srcpos);

    template <typename T2>
    void addmult_cols (size_type dstpos, size_type srcpos, const T2& y);

    template <typename T1, typename T2>
    void addmult_cols
    (
        size_type dstpos, const T1& x,
        size_type srcpos, const T2& y
    );

    template <typename T1, typename T2>
    void addmult_cols
    (
        size_type dstpos, const T1& x,
        int srcpos, const T2& y
    )
    { addmult_cols(dstpos, x, size_type(srcpos), y); }

    template <typename T1, typename V, typename T2>
    void addmult_cols
    (
        size_type dstpos, const T1& x,
        const V& srcvec, const T2& y
    );

    vector<T, true> copy_col (size_type i) const
    {
        vector<T, true> res;
        return copy_col(i, res);
    }


    template <typename V>
    V& copy_col (size_type i, V& res) const;

    template <typename SV>
    matrix<T, true> copy_cols (const SV& sv) const
    {
        matrix<T, true> res;
        return copy_cols(sv, res);
    }

    template <typename SV, typename M>
    M& copy_cols (const SV& sv, M& res) const;
    
    template <typename SV>
    void erase_cols (const SV& sv);

    template <typename SV, typename M>
    M& take_cols (const SV& sv, M& res)
    {
        copy_cols(sv, res);
        erase_cols(sv);
        return res;
    }

    template <typename SV>
    matrix<T, true> take_cols (const SV& sv)
    {
        matrix<T, true> res;
        return take_cols(sv, res);
    }

    template <typename V>
    V& take_col (size_type i, V& res)
    {
        copy_col(i, res);
        erase_col(i);
        return res;
    }

    vector<T, true> take_col (size_type i)
    {
        vector<T, true> res;
        return take_col(i, res);
    }

    template <typename V>
    void assign_col (size_type j, const V& v)
    {
        ARAGELI_ASSERT_0(j < ncols());
        ARAGELI_ASSERT_0(v.size() >= nrows());
        
        for(size_type i = 0; i < nrows(); ++i)
            el(i, j) = v[i];
    }

    // @}


    template <typename M>
    void insert_matrix_bottom (const M& m)
    {
        ARAGELI_ASSERT_0(m.ncols() == ncols());
        
        // TODO: Make it faster!
        for(size_type i = 0; i < m.nrows(); ++i)
            insert_row(nrows(), m.copy_row(i));
    }


    template <typename RS, typename CS>
    matrix<T, true> copy_submatrix (const RS& rs, const CS& cs) const
    {
        matrix<T, true> res;
        return copy_submatrix(rs, cs, res);
    }

    template <typename RS, typename CS, typename M>
    M& copy_submatrix (const RS& rs, const CS& cs, M& res) const;

    template <typename RS, typename CS>
    void erase_submatrix (const RS& rs, const CS& cs);

    template <typename RS, typename CS, typename M>
    M& take_submatrix (const RS& rs, const CS& cs, M& res)
    {
        copy_submatrix(rs, cs, res);
        erase_submatrix(rs, cs);
        return res;
    }

    template <typename RS, typename CS>
    matrix<T, true> take_submatrix (const RS& rs, const CS& cs)
    {
        matrix<T, true> res;
        return take_submatrix(rs, cs, res);
    }



    bool unique () { return store.unique(); }


    bool unique_clear () { return store.unique_clear(); }

private:

    Rep& rep () { return store.value(); }
    const Rep& rep () const { return store.value(); }

    Mem& mem () { return rep().mem; }
    const Mem& mem () const { return rep().mem; }

    reference mem (size_type i) { return mem()[i]; }
    const_reference mem (size_type i) const { return mem()[i]; }

    bool do_pack ();

    refcntr_proxy<Rep, REFCNT> store;
};


#if 0

template <typename TE2, bool REFCNTE2>
struct type_pair_traits<big_int, matrix<TE2, REFCNTE2> >
    : public type_pair_traits_for_unknown_and_vector_matrix
        <big_int, matrix<TE2, REFCNTE2> > {};

#endif

template <typename T, bool REFCNT>
matrix<T, true> transpose (const matrix<T, REFCNT>& a)
{
    matrix<T, true> res = a;
    res.transpose();
    return res;
}



template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2
>
inline int cmp
(
    const matrix<T1, REFCNT1>& a,
    const matrix<T2, REFCNT2>& b
)
{
    if(a.nrows() < b.nrows())return -1;
    if(a.nrows() > b.nrows())return +1;
    if(a.ncols() < b.ncols())return -1;
    if(a.ncols() > b.ncols())return +1;
    return _Internal::aggregate_cmp(a.begin(), a.end(), b.begin(), b.end());
}



#define _ARAGELI_MATRIX_CMP(OPER)                                           \
    template                                                                \
    <                                                                       \
        typename T1, bool REFCNT1,                                          \
        typename T2, bool REFCNT2                                           \
    >                                                                       \
    inline bool operator OPER                                               \
    (const matrix<T1, REFCNT1>& a, const matrix<T2, REFCNT2>& b)            \
    { return cmp(a, b) OPER 0; }                                            \
                                                                            \
    template                                                                \
    <                                                                       \
        typename T1, bool REFCNT1,                                          \
        typename T2                                                         \
    >                                                                       \
    inline bool operator OPER                                               \
    (const matrix<T1, REFCNT1>& a, const T2& b)                             \
    {                                                                       \
        return cmp                                                          \
        (                                                                   \
            a,                                                              \
            matrix<T2, false>(a.nrows(), a.ncols, b, fromsize)      \
        ) OPER 0;                                                           \
    }                                                                       \
                                                                            \
    template                                                                \
    <                                                                       \
        typename T1,                                                        \
        typename T2, bool REFCNT2                                           \
    >                                                                       \
    inline bool operator OPER                                               \
    (const T1& a, const matrix<T2, REFCNT2>& b)                             \
    {                                                                       \
        return cmp                                                          \
        (                                                                   \
            matrix<T2, false>(b.nrows(), b.ncols(), a, fromsize),   \
            b                                                               \
        ) OPER 0;                                                           \
    }


// @{

_ARAGELI_MATRIX_CMP(<=)
_ARAGELI_MATRIX_CMP(>=)
_ARAGELI_MATRIX_CMP(<)
_ARAGELI_MATRIX_CMP(>)
_ARAGELI_MATRIX_CMP(==)
_ARAGELI_MATRIX_CMP(!=)

// @}

#undef _ARAGELI_MATRIX_CMP


#define _ARAGELI_GENERIC_MATRIX_OPERATORS1(MNEM, OPERB)             \
template <typename T1, typename T2, typename T3>                    \
T3& each_##MNEM##_matrix_matrix (const T1& a, const T2& b, T3& res) \
{                                                                   \
    ARAGELI_ASSERT_0(a.size() == b.size());                         \
    ARAGELI_ASSERT_0(res.size() == a.size());                       \
                                                                    \
    typename T1::const_iterator ia = a.begin(), iend = a.end();     \
    typename T2::const_iterator ib = b.begin();                     \
    typename T3::iterator ic = res.begin();                         \
                                                                    \
    for(; ia != iend; ++ia, ++ib, ++ic)                             \
        *ic = *ia OPERB *ib;                                        \
                                                                    \
    return res;                                                     \
}                                                                   \
                                                                    \
template <typename T1, typename T2, typename T3>                    \
T3& MNEM##_matrix_scalar (const T1& a, const T2& b, T3& res)        \
{                                                                   \
    ARAGELI_ASSERT_0(res.size() == a.size());                       \
                                                                    \
    typename T1::const_iterator ia = a.begin(), iend = a.end();     \
    typename T3::iterator ic = res.begin();                         \
                                                                    \
    for(; ia != iend; ++ia, ++ic)                                   \
        *ic = *ia OPERB b;                                          \
                                                                    \
    return res;                                                     \
}                                                                   \
                                                                    \
template <typename T1, typename T2, typename T3>                    \
T3& MNEM##_scalar_matrix (const T1& a, const T2& b, T3& res)        \
{                                                                   \
    ARAGELI_ASSERT_0(res.size() == b.size());                       \
                                                                    \
    typename T2::const_iterator ib = b.begin(), iend = b.end();     \
    typename T3::iterator ic = res.begin();                         \
                                                                    \
    for(; ib != iend; ++ib, ++ic)                                   \
        *ic = a OPERB *ib;                                          \
                                                                    \
    return res;                                                     \
}


_ARAGELI_GENERIC_MATRIX_OPERATORS1(add, +)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(sub, -)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(mul, *)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(div, /)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(mod, %)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(bitand, &)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(bitor, |)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(bitxor, ^)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(shl, <<)
_ARAGELI_GENERIC_MATRIX_OPERATORS1(shr, >>)


#define _ARAGELI_GENERIC_MATRIX_OPERATORS2(MNEM)                \
template <typename T1, typename T2, typename T3>                \
T3& MNEM##_matrix_matrix (const T1& a, const T2& b, T3& res)    \
{ return each_##MNEM##_matrix_matrix(a, b, res); }

template <typename T1, typename T2, typename T3>
T3& mul_matrix_matrix (const T1& a, const T2& b, T3& res)
{ res = a; return res.mul_matrix(b); }

template <typename T1, typename T2, typename T3>
T3& div_matrix_matrix (const T1& a, const T2& b, T3& res)
{ res = a; return res.div_matrix(b); }


_ARAGELI_GENERIC_MATRIX_OPERATORS2(add)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(sub)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(mod)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(bitand)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(bitor)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(bitxor)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(shl)
_ARAGELI_GENERIC_MATRIX_OPERATORS2(shr)


#undef _ARAGELI_GENERIC_MATRIX_OPERATORS1
#undef _ARAGELI_GENERIC_MATRIX_OPERATORS2


#define _ARAGELI_MATRIX_BINARY_BY_UNARY(MNEM, OPERB)    \
    template                                            \
    <                                                   \
        typename T1, bool REFCNT1,                      \
        typename T2, bool REFCNT2                       \
    >                                                   \
    inline matrix<T1, REFCNT1> operator OPERB           \
    (                                                   \
        const matrix<T1, REFCNT1>& a,                   \
        const matrix<T2, REFCNT2>& b                    \
    )                                                   \
    {                                                   \
        matrix<T1, REFCNT1> t                           \
            (a.nrows(), a.ncols(), fromsize);   \
        return MNEM##_matrix_matrix(a, b, t);           \
    }                                                   \
                                                        \
    template                                            \
    <                                                   \
        typename T1, bool REFCNT1,                      \
        typename T2                                     \
    >                                                   \
    inline matrix<T1, REFCNT1> operator OPERB           \
    (                                                   \
        const matrix<T1, REFCNT1>& a,                   \
        const T2& b                                     \
    )                                                   \
    {                                                   \
        matrix<T1, REFCNT1> t                           \
            (a.nrows(), a.ncols(), fromsize);   \
        return MNEM##_matrix_scalar(a, b, t);           \
    }                                                   \
                                                        \
    template                                            \
    <                                                   \
        typename T1,                                    \
        typename T2, bool REFCNT2                       \
    >                                                   \
    inline matrix<T2, REFCNT2> operator OPERB           \
    (                                                   \
        const T1& a,                                    \
        const matrix<T2, REFCNT2>& b                    \
    )                                                   \
    {                                                   \
        matrix<T2, REFCNT2> t                           \
            (b.nrows(), b.ncols(), fromsize);   \
        return MNEM##_scalar_matrix(a, b, t);           \
    }                                                   \
                                                        \
    template <typename T1, bool REFCNT1>                \
    inline matrix<T1, REFCNT1> operator OPERB           \
    (                                                   \
        const matrix<T1, REFCNT1>& a,                   \
        const T1& b                                     \
    )                                                   \
    {                                                   \
        matrix<T1, REFCNT1> t                           \
            (a.nrows(), a.ncols(), fromsize);   \
        return MNEM##_matrix_scalar(a, b, t);           \
    }                                                   \
                                                        \
    template <typename T2, bool REFCNT2>                \
    inline matrix<T2, REFCNT2> operator OPERB           \
    (                                                   \
        const T2& a,                                    \
        const matrix<T2, REFCNT2>& b                    \
    )                                                   \
    {                                                   \
        matrix<T2, REFCNT2> t                           \
            (b.nrows(), b.ncols(), fromsize);   \
        return MNEM##_scalar_matrix(a, b, t);           \
    }                                                   \
                                                        \
    /* BEGIN OF TEMPORARY EXTENTION */                  \
    template <typename T1, bool REFCNT1>                \
    inline matrix<rational<T1>, REFCNT1>                \
    operator OPERB                                      \
    (                                                   \
        const matrix<rational<T1>, REFCNT1>& a,         \
        const rational<T1>& b                           \
    )                                                   \
    {                                                   \
        matrix<rational<T1>, REFCNT1> t                 \
            (a.nrows(), a.ncols(), fromsize);   \
        return MNEM##_matrix_scalar(a, b, t);           \
    }                                                   \
                                                        \
    template <typename T2, bool REFCNT2>                \
    inline matrix<rational<T2>, REFCNT2>                \
    operator OPERB                                      \
    (                                                   \
        const rational<T2>& a,                          \
        const matrix<rational<T2>, REFCNT2>& b          \
    )                                                   \
    {                                                   \
        matrix<rational<T2>, REFCNT2> t                 \
            (b.nrows(), b.ncols(), fromsize);   \
        return MNEM##_scalar_matrix(a, b, t);           \
    }                                                   \
                                                        \
    template                                            \
    <                                                   \
            typename T1, typename D1,                   \
            bool REFCNT1, bool REFCNT2                  \
    >                                                   \
    inline matrix<sparse_polynom                        \
        <T1, D1, REFCNT1>, REFCNT2>                     \
    operator OPERB                                      \
    (                                                   \
        const matrix<sparse_polynom                     \
            <T1, D1, REFCNT1>, REFCNT2>& a,             \
        const sparse_polynom<T1, D1, REFCNT1>& b        \
    )                                                   \
    {                                                   \
        matrix<sparse_polynom                           \
            <T1, D1, REFCNT1>, REFCNT2> t               \
            (a.nrows(), a.ncols(), fromsize);   \
        return MNEM##_matrix_scalar(a, b, t);           \
    }                                                   \
                                                        \
    template                                            \
    <                                                   \
            typename T1, typename D1,                   \
            bool REFCNT1, bool REFCNT2                  \
    >                                                   \
    inline matrix<sparse_polynom                        \
        <T1, D1, REFCNT1>, REFCNT2>                     \
    operator OPERB                                      \
    (                                                   \
        const sparse_polynom<T1, D1, REFCNT1>& a,       \
        const matrix<sparse_polynom                     \
            <T1, D1, REFCNT1>, REFCNT2>& b              \
    )                                                   \
    {                                                   \
        matrix<sparse_polynom                           \
            <T1, D1, REFCNT1>, REFCNT2> t               \
            (b.nrows(), b.ncols(), fromsize);   \
        return MNEM##_scalar_matrix(a, b, t);           \
    }
    /* END OF TEMPORARY EXTENTION */


// @{

_ARAGELI_MATRIX_BINARY_BY_UNARY(add, +)
_ARAGELI_MATRIX_BINARY_BY_UNARY(sub, -)
_ARAGELI_MATRIX_BINARY_BY_UNARY(mul, *)
_ARAGELI_MATRIX_BINARY_BY_UNARY(div, /)
_ARAGELI_MATRIX_BINARY_BY_UNARY(mod, %)
_ARAGELI_MATRIX_BINARY_BY_UNARY(bitand, &)
_ARAGELI_MATRIX_BINARY_BY_UNARY(bitor, |)
_ARAGELI_MATRIX_BINARY_BY_UNARY(bitxor, ^)
_ARAGELI_MATRIX_BINARY_BY_UNARY(shl, <<)
_ARAGELI_MATRIX_BINARY_BY_UNARY(shr, >>)

// @}


#undef _ARAGELI_MATRIX_BINARY_BY_UNARY


//#define _ARAGELI_MATRIX_BINARY_BY_UNARY(OPERU, OPERB) \
//  template                                            \
//  <                                                   \
//      typename T1, bool REFCNT1,                      \
//      typename T2, bool REFCNT2                       \
//  >                                                   \
//  inline matrix<T1, REFCNT1> operator OPERB           \
//  (                                                   \
//      const matrix<T1, REFCNT1>& a,                   \
//      const matrix<T2, REFCNT2>& b                    \
//  )                                                   \
//  {                                                   \
//      matrix<T1, REFCNT1> t(a);                       \
//      return t OPERU b;                               \
//  }
//
//  
//#define _ARAGELI_MATRIX_BINARY_SCALAR(OPERU, OPERB)           \
//  template                                                \
//  <                                                       \
//      typename T1, bool REFCNT1,                          \
//      typename T2                                         \
//  >                                                       \
//  inline matrix<T1, REFCNT1> operator OPERB               \
//  (                                                       \
//      const matrix<T1, REFCNT1>& a,                       \
//      const T2& b                                         \
//  )                                                       \
//  {                                                       \
//      matrix<T1, REFCNT1> t(a);                           \
//      return t OPERB##= b;                                \
//  }                                                       \
//                                                          \
//  template                                                \
//  <                                                       \
//      typename T1,                                        \
//      typename T2, bool REFCNT2                           \
//  >                                                       \
//  inline matrix<T2, REFCNT2> operator OPERB               \
//  (                                                       \
//      const T1& a,                                        \
//      const matrix<T2, REFCNT2>& b                        \
//  )                                                       \
//  {                                                       \
//      /* ñîçäàâàòü ìàòðèöó ñêàëÿðîâ íå ðàöèîíàëüíî */     \
//      matrix<T2, REFCNT2> t(b.nrows(), b.ncols(), a);     \
//      return t.OPERU(b);                                  \
//  }
//
//
//
//_ARAGELI_MATRIX_BINARY_BY_UNARY(+=, +)
//_ARAGELI_MATRIX_BINARY_BY_UNARY(-=, -)
//_ARAGELI_MATRIX_BINARY_BY_UNARY(*=, *)
//_ARAGELI_MATRIX_BINARY_BY_UNARY(/=, /)
//
//_ARAGELI_MATRIX_BINARY_SCALAR(each_add, +)
//_ARAGELI_MATRIX_BINARY_SCALAR(each_sub, -)
//_ARAGELI_MATRIX_BINARY_SCALAR(each_mul, *)
//_ARAGELI_MATRIX_BINARY_SCALAR(each_div, /)
//_ARAGELI_MATRIX_BINARY_SCALAR(each_mod, %)
//


//template
//<
//  typename T1, bool REFCNT1,
//  typename T2, bool REFCNT2
//>
//bool operator==
//(
//  const matrix<T1, REFCNT1>& a,
//  const matrix<T2, REFCNT2>& b
//);
//
//
//template
//<
//  typename T1, bool REFCNT1,
//  typename T2, bool REFCNT2
//>
//inline bool operator!=
//(
//  const matrix<T1, REFCNT1>& a,
//  const matrix<T2, REFCNT2>& b
//)
//{ return !(a == b); }



template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2
>
vector<T2, REFCNT2> operator*
(
    const matrix<T1, REFCNT1>& m,
    const vector<T2, REFCNT2>& x
);



template
<
    typename T2, bool REFCNT2,
    typename T1, bool REFCNT1
>
vector<T2, REFCNT2> operator*
(
    const vector<T2, REFCNT2>& x,
    const matrix<T1, REFCNT1>& m
);


extern const char* matrix_output_list_first_bracket_default;
extern const char* matrix_output_list_second_bracket_default;
extern const char* matrix_output_list_row_separator_default;
extern const char* matrix_output_list_first_row_bracket_default;
extern const char* matrix_output_list_second_row_bracket_default;
extern const char* matrix_output_list_col_separator_default;
extern const char* matrix_input_list_first_bracket_default;
extern const char* matrix_input_list_second_bracket_default;
extern const char* matrix_input_list_row_separator_default;
extern const char* matrix_input_list_first_row_bracket_default;
extern const char* matrix_input_list_second_row_bracket_default;
extern const char* matrix_input_list_col_separator_default;
extern const char* matrix_output_aligned_left_col_default;
extern const char* matrix_output_aligned_right_col_default;
extern const char* matrix_output_aligned_inter_col_default;



template <typename T, bool REFCNT>
std::ostream& output_list
(
    std::ostream& out,
    const matrix<T, REFCNT>& x,
    const char* first_bracket       = matrix_output_list_first_bracket_default,
    const char* second_bracket      = matrix_output_list_second_bracket_default,
    const char* row_separator       = matrix_output_list_row_separator_default,
    const char* first_row_bracket   = matrix_output_list_first_row_bracket_default,
    const char* second_row_bracket  = matrix_output_list_second_row_bracket_default,
    const char* col_separator       = matrix_output_list_col_separator_default
);


template <typename T, bool REFCNT>
std::istream& input_list
(
    std::istream& out,
    matrix<T, REFCNT>& x,
    const char* first_bracket       = matrix_input_list_first_bracket_default,
    const char* second_bracket      = matrix_input_list_second_bracket_default,
    const char* row_separator       = matrix_input_list_row_separator_default,
    const char* first_row_bracket   = matrix_input_list_first_row_bracket_default,
    const char* second_row_bracket  = matrix_input_list_second_row_bracket_default,
    const char* col_separator       = matrix_input_list_col_separator_default
);


template <typename T, bool REFCNT>
std::ostream& output_aligned
(
    std::ostream& out,
    const matrix<T, REFCNT>& x,
    const char* left_col    = matrix_output_aligned_left_col_default,
    const char* right_col   = matrix_output_aligned_right_col_default,
    const char* inter_col   = matrix_output_aligned_inter_col_default
);


template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2
>
std::ostream& output_aligned_ver_pair
(
    std::ostream& out,
    const matrix<T1, REFCNT1>& m1,
    const matrix<T2, REFCNT2>& m2,
    const char* left_col = "|| ",
    const char* right_col = " ||",
    const char* inter_col = " ",
    const char* hsplitter = "-"
);


template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2
>
std::ostream& output_aligned_hor_pair
(
    std::ostream& out,
    const matrix<T1, REFCNT1>& m1,
    const matrix<T2, REFCNT2>& m2,
    const char* left_col = "|| ",
    const char* right_col = " ||",
    const char* inter_col = " ",
    const char* vsplitter = " | "
);



template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2,
    typename T3, bool REFCNT3
>
std::ostream& output_aligned_corner_triplet_br
(
    std::ostream& out,
    const matrix<T1, REFCNT1>& m1,
    const matrix<T2, REFCNT2>& m2,
    const matrix<T3, REFCNT3>& m3,
    const char* left_col = "|| ",
    const char* right_col = " ||",
    const char* inter_col = " ",
    const char* vsplitter = " | ",
    const char* hsplitter = "-"
);


template <typename T, bool REFCNT>
std::ostream& operator<< (std::ostream& out, const matrix<T, REFCNT>& x)
{ return output_list(out, x); }

template <typename T, bool REFCNT>
std::istream& operator>> (std::istream& out, matrix<T, REFCNT>& x)
{ return input_list(out, x); }


template <typename T, bool REFCNT>
std::ofstream& operator<< (std::ofstream& out, const matrix<T, REFCNT>& x)
{ return static_cast<std::ofstream&>(output_list(out, x)); }

template <typename T, bool REFCNT>
std::ifstream& operator>> (std::ifstream& out, matrix<T, REFCNT>& x)
{ return static_cast<std::ifstream&>(input_list(out, x)); }


// BEGIN OF TEMPORARY EXTENTION

template <typename T, bool REFCNT>
std::ostream& operator<< (std::ostringstream& out, const matrix<T, REFCNT>& x)
{ return output_list(out, x); }

template <typename T, bool REFCNT>
std::istream& operator>> (std::istringstream& out, matrix<T, REFCNT>& x)
{ return input_list(out, x); }

// END OF TEMPORARY EXTENTION


template <typename T, bool REFCNT>
struct factory<matrix<T, REFCNT> >
{
private:

    typedef matrix<T, REFCNT> TTT;

public:

    static const bool is_specialized = true;
    
    static const TTT& unit ()
    {
        static const TTT unit_s(1, eye);
        return unit_s;
    }
    
    static TTT unit (const TTT& pat)
    { return TTT(pat.nrows(), Arageli::unit<T>(), diag); }
    
    static const TTT& opposite_unit ()
    {
        static const TTT opposite_unit_s
            (1, Arageli::opposite_unit<T>(), diag);
        return opposite_unit_s;
    }
    
    static TTT opposite_unit (const TTT& pat)
    { return TTT(pat.nrows(), Arageli::opposite_unit<T>(), diag); }
    
    static const TTT& null ()
    {
        static const TTT null_s(1, Arageli::null<T>(), fromval);
        return null_s;
    }
    
    static TTT null (const TTT& pat)
    { return TTT(pat.nrows(), Arageli::null<T>(), diag); }

};


template <typename T, bool REFCNT>
inline matrix<T, REFCNT> opposite
(const matrix<T, REFCNT>& x)
{ return -x; }

template <typename T, bool REFCNT>
inline matrix<T, REFCNT>& opposite
(const matrix<T, REFCNT>& x, matrix<T, REFCNT>* y)
{ return (*y = x).opposite(); }

template <typename T, bool REFCNT>
inline matrix<T, REFCNT>& opposite
(matrix<T, REFCNT>* x)
{ return x->opposite(); }

//template <typename T, bool REFCNT>
//inline matrix<T, REFCNT> inverse
//(const matrix<T, REFCNT>& x)
//{
//  matrix<T, REFCNT> t(x);
//  return x.inverse();
//}

template <typename T, bool REFCNT>
inline matrix<T, REFCNT>& inverse
(const matrix<T, REFCNT>& x, matrix<T, REFCNT>* y)
{ return (*y = x).inverse(); }

template <typename T, bool REFCNT>
inline matrix<T, REFCNT>& inverse
(matrix<T, REFCNT>* x)
{ return x->inverse(); }


template <typename T, bool REFCNT>
inline bool is_unit (const matrix<T, REFCNT>& x)
{ return x.is_unit(); }

template <typename T, bool REFCNT>
inline bool is_opposite_unit (const matrix<T, REFCNT>& x)
{ return x.is_opposite_unit(); }

template <typename T, bool REFCNT>
inline bool is_null (const matrix<T, REFCNT>& x)
{ return x.is_null(); }


template <typename T, bool REFCNT>
struct type_traits<matrix<T, REFCNT> >
    : public type_traits_default<matrix<T, REFCNT> >
{
    static const bool is_specialized = type_traits<T>::is_specialized;

    static const bool is_number = false;
    static const bool is_integer_number = false;
    static const bool is_polynom = false;
    static const bool is_real_number = false;
    static const bool is_rational_number = false;
    static const bool is_complex_number = false;
    static const bool is_ring = type_traits<T>::is_ring;    // +, -, *, null, unit
    static const bool is_field = false;
    static const bool is_finite = type_traits<T>::is_finite;
    static const bool is_additive_group = type_traits<T>::is_additive_group;    // +, -, null
    static const bool is_multiplicative_group = false;
    static const bool has_zero_divisor = type_traits<T>::has_zero_divisor;
    static const bool is_integer_modulo_ring = false;
    static const bool is_matrix = true;
    static const bool is_vector = false;
    
    static const bool has_commutative_multiplication =
                        type_traits<T>::has_commutative_multiplication &&
                        type_traits<T>::has_commutative_addition;
    
    static const bool has_commutative_addition =
                        type_traits<T>::has_commutative_addition;

    static const bool has_null = type_traits<T>::has_null;
    
    static const bool has_unit =
                        type_traits<T>::has_unit &&
                        type_traits<T>::has_null;
    
    static const bool has_opposite_unit =
                        type_traits<T>::has_opposite_unit &&
                        type_traits<T>::has_null;

    static const bool is_aggregate = true;
    
    typedef T element_type;

    typedef type_category::dense_matrix category_type;
    static const category_type category_value;
};


template <typename T, bool REFCNT>
const type_category::dense_matrix
    type_traits<matrix<T, REFCNT> >::category_value =
    type_category::dense_matrix();


template <typename T, bool REFCNT>
inline std::ostream& output_polynom_first
(std::ostream& out, const matrix<T, REFCNT>& x)
{ return output_list(out, x); }

template <typename T, bool REFCNT>
inline std::ostream& output_polynom_internal
(std::ostream& out, const matrix<T, REFCNT>& x)
{ return output_list(out << "+", x); }

template <typename T, bool REFCNT>
inline std::ostream& output_pow
(std::ostream& out, const matrix<T, REFCNT>& x)
{ return output_list(out, x); }

template <typename T, bool REFCNT>
inline std::istream& input_polynom_first
(std::istream& in, matrix<T, REFCNT>& x)
{ return input_list(in, x); }

template <typename T, bool REFCNT>
std::istream& input_polynom_internal
(std::istream& in, matrix<T, REFCNT>& x);

template <typename T, bool REFCNT>
inline std::istream& input_pow
(std::istream& in, matrix<T, REFCNT>& x)
{ return input_polynom_first(in, x); }


template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2,
    typename T3, bool REFCNT3
>
void fill_submatrix_col
(
    const matrix<T1, REFCNT1>& orig,
    const vector<T2, REFCNT2>& indexes,
    matrix<T3, REFCNT3>& res
)
{
    res.assign_fromsize(orig.nrows(), indexes.size());

    for(typename vector<T2, REFCNT2>::size_type j = 0; j < indexes.size(); ++j)
    {
        const T2& index = indexes[j];

        ARAGELI_ASSERT_0(index >= 0);
        ARAGELI_ASSERT_0(index < orig.ncols());

        for(typename matrix<T3, REFCNT3>::size_type i = 0; i < res.nrows(); ++i)
            res(i, j) = orig(i, index);
    }
}


template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2,
    typename T3, bool REFCNT3
>
void fill_submatrix_row
(
    const matrix<T1, REFCNT1>& orig,
    const vector<T2, REFCNT2>& indexes,
    matrix<T3, REFCNT3>& res
)
{
    res.assign_fromsize(indexes.size(), orig.ncols());

    for(typename vector<T2, REFCNT2>::size_type i = 0; i < indexes.size(); ++i)
    {
        const T2& index = indexes[i];

        ARAGELI_ASSERT_0(index >= 0);
        ARAGELI_ASSERT_0(index < orig.nrows());

        for(typename matrix<T3, REFCNT3>::size_type j = 0; j < res.ncols(); ++j)
            res(i, j) = orig(index, j);
    }
}


}


namespace std
{

template
<
    typename T1, bool REFCNT1,
    typename T2, bool REFCNT2
>
inline void swap
(
    Arageli::matrix<T1, REFCNT1>& a,
    Arageli::matrix<T2, REFCNT2>& b
)
{ a.swap(b); }


template <typename T1, bool REFCNT1>
inline void swap
(
    Arageli::matrix<T1, REFCNT1>& a,
    Arageli::matrix<T1, REFCNT1>& b
)
{ a.swap(b); }

}


#ifdef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
    #define ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_MATRIX
    #include "matrix.cpp"
    #undef  ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_MATRIX
#endif


#endif

---