matrix.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

   * a view a deep copy is not made, rather the view simply provides
   * access to the extant block of data underlying the copied object.  
   * Furthermore, when
   * you assign to a view, you overwrite the data the view is
   * currently pointing to, rather than generating a new data block.
   * Together, these behaviors allow
   * for matrices that view portions of other matrices
   * (submatrices) and submatrix assignment.  Views do not guarantee
   * unit stride and may even logically access their elements in a
   * different order than they are stored in memory.  Copying between
   * concretes and views is fully supported and often transparent to
   * the user.
   *
   * There is a fundamental trade-off between concrete matrices and
   * views.  Concrete matrices are simpler to work with, but not
   * as flexible as views.  Because they always have unit stride,
   * concrete matrices
   * have fast iterators and access operators but, because they must
   * always be copied deeply, they provide slow copies (for example,
   * copy construction of a Matrix returned from a function wastes
   * cycles).  Views are more flexible but also somewhat more
   * complicated to program with.  Furthermore, because they cannot
   * guarantee unit stride, their iterators and access operations are
   * somewhat slower than those for concrete matrices.  On the other
   * hand, they provide very fast copies.  The average Scythe user may
   * find that she virtually never works with views directly (although
   * they can be quite useful in certain situations) but they provide
   * a variety of functionality underneath the hood of the library and
   * underpin many common operations.
   *
   * \note
   * The Matrix interface is split between two classes: this Matrix
   * class and Matrix_base, which Matrix extends.  Matrix_base
   * includes a range of accessors that provide the programmer with
   * information about such things as the dimensionality of Matrix
   * objects.
   */

	template <typename T_type = double, matrix_order ORDER = Col, 
            matrix_style STYLE = Concrete>
	class Matrix : public Matrix_base<ORDER, STYLE>,
								 public DataBlockReference<T_type>
	{
		public:
			/**** TYPEDEFS ****/

			/* Iterator types */

      /*! \brief Random Access Iterator type.
       *
       * This typedef for matrix_random_access_iterator provides a
       * convenient shorthand for the default, and most general,
       * Matrix iterator type.
       *
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef matrix_random_access_iterator<T_type, ORDER, ORDER, STYLE>
        iterator;

      /*! \brief Const Random Access Iterator type.
       *
       * This typedef for const_matrix_random_access_iterator provides
       * a convenient shorthand for the default, and most general,
       * Matrix const iterator type.
       *
       * \see iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef const_matrix_random_access_iterator<T_type, ORDER, ORDER,
                                                  STYLE> const_iterator;

      /*! \brief Reverse Random Access Iterator type.
       *
       * This typedef uses std::reverse_iterator to describe a
       * reversed matrix_random_access_iterator type.  This is the
       * reverse of iterator.
       *
       * \see iterator
       * \see const_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef typename 
        std::reverse_iterator<matrix_random_access_iterator<T_type, 
                              ORDER, ORDER, STYLE> > reverse_iterator;

      /*! \brief Reverse Const Random Access Iterator type.
       *
       * This typedef uses std::reverse_iterator to describe a
       * reversed const_matrix_random_access_iterator type.  This is
       * the reverse of const_iterator.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef typename
				std::reverse_iterator<const_matrix_random_access_iterator
                              <T_type, ORDER, ORDER, STYLE> >
				const_reverse_iterator;

      /*! \brief Forward Iterator type.
       *
       * This typedef for matrix_forward_iterator provides
       * a convenient shorthand for a fast (when compared to
       * matrix_random_access_iterator) Matrix iterator type.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
      typedef matrix_forward_iterator<T_type, ORDER, ORDER, STYLE>
        forward_iterator;

      /*! \brief Const Forward Iterator type.
       *
       * This typedef for const_matrix_forward_iterator provides a
       * convenient shorthand for a fast (when compared to
       * const_matrix_random_access_iterator) const Matrix iterator
       * type.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
      typedef const_matrix_forward_iterator<T_type, ORDER, ORDER, STYLE>
        const_forward_iterator;

      /*! \brief Bidirectional Iterator type.
       *
       * This typedef for matrix_bidirectional_iterator provides
       * a convenient shorthand for a compromise (with speed and
       * flexibility between matrix_random_access_iterator and
       * matrix_forward_iterator) Matrix iterator type.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef matrix_bidirectional_iterator<T_type, ORDER, ORDER, STYLE>
        bidirectional_iterator;

      /*! \brief Const Bidirectional Iterator type.
       *
       * This typedef for const_matrix_bidirectional_iterator provides
       * a convenient shorthand for a compromise (with speed and
       * flexibility between const_matrix_random_access_iterator and
       * const_matrix_forward_iterator) const Matrix iterator type.
       * 
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see reverse_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef const_matrix_bidirectional_iterator<T_type, ORDER, ORDER,
                                  STYLE> const_bidirectional_iterator;

      /*! \brief Const Bidirectional Iterator type.
       *
       * This typedef uses std::reverse_iterator to describe a
       * reversed matrix_bidirectional_iterator type.  This is
       * the reverse of bidirectional_iterator.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see const_reverse_bidirectional_iterator
       */
			typedef typename 
        std::reverse_iterator<matrix_bidirectional_iterator<T_type, 
                ORDER, ORDER, STYLE> > reverse_bidirectional_iterator;

      /*! \brief Reverse Const Bidirectional Iterator type.
       *
       * This typedef uses std::reverse_iterator to describe a
       * reversed const_matrix_bidirectional_iterator type.  This is
       * the reverse of const_bidirectional_iterator.
       *
       * \see iterator
       * \see const_iterator
       * \see reverse_iterator
       * \see const_reverse_iterator
       * \see forward_iterator
       * \see const_forward_iterator
       * \see reverse_forward_iterator
       * \see const_reverse_forward_iterator
       * \see bidirectional_iterator
       * \see const_bidirectional_iterator
       * \see reverse_bidirectional_iterator
       */
			typedef typename
				std::reverse_iterator<const_matrix_bidirectional_iterator
                              <T_type, ORDER, ORDER, STYLE> >
				const_reverse_bidirectional_iterator;

      /*!\brief The Matrix' element type.
       *
       * This typedef describes the element type (T_type) of this
       * Matrix.
       */
      typedef T_type ttype;
		
		private:
			/* Some convenience typedefs */
			typedef DataBlockReference<T_type> DBRef;
			typedef Matrix_base<ORDER, STYLE> Base;
			
		public:
			/**** CONSTRUCTORS ****/

			/*! \brief Default constructor.
       *
       * The default constructor creates an empty/null matrix.  Using
       * null matrices in operations will typically cause errors; this
       * constructor exists primarily for initialization within
       * aggregate types.
       *
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       *
       * \b Example:
       * \include example.matrix.constructor.default.cc
       */
			Matrix ()
				:	Base (),
					DBRef ()
			{
			}

			/*! \brief Parameterized type constructor.
       *
       * Creates a 1x1 matrix (scalar).
       *
       * \param element The scalar value of the constructed Matrix.
       *
       * \see Matrix()
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.ptype.cc
			 */
			Matrix (T_type element)
				:	Base (1, 1),
					DBRef (1)
			{
				data_[Base::index(0)] = element;  // ALWAYS use index()
			}

      /*! \brief Standard constructor.
       *
       * The standard constructor creates a rowsXcols Matrix, filled
       * with zeros by default.  Optionally, you can leave the Matrix
       * uninitialized, or choose a different fill value.
       * 
       * \param rows The number of rows in the Matrix.
       * \param cols The number of columns in the Matrix.
       * \param fill Indicates whether or not the Matrix should be
       * initialized.
       * \param fill_value The scalar value to fill the Matrix with
       * when fill == true.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.standard.cc
       */
			Matrix (uint rows, uint cols, bool fill = true,
					T_type fill_value = 0)
				:	Base (rows, cols),
					DBRef (rows * cols)
			{
        // TODO Might use iterator here for abstraction.
				if (fill)
					for (uint i = 0; i < Base::size(); ++i)
						data_[Base::index(i)] = fill_value; // we know data contig
			}

      /*! \brief Iterator constructor.
       *
			 * Creates a \a rows X \a cols matrix, filling it sequentially
			 * (based on this template's matrix_order) with values