matrix.h

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

				return res;
			}

			/* Make this matrix a copy of another. The matrix retains its
       * own order and style in this case, because that can't change
       * at run time.
       */
      /*!\brief Make this Matrix a copy of another.
       *
       * Converts this Matrix into a deep copy of another Matrix.
       * This Matrix retains its own matrix_order and matrix_style but
       * contains copies of M's elements and becomes the same size and
       * shape as M.  Calling this method automatically detaches this
       * Matrix from its previous DataBlock before copying.
       *
       * \param M The Matrix to copy.
       *
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see copy()
       * \see reference(const Matrix<T_type, O, S> &)
       * \see detach()
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.copyother.cc
       */
      template <matrix_order O, matrix_style S>
			inline void copy (const Matrix<T_type, O, S>& M)
			{
				resize2Match(M);
        scythe::copy<ORDER,ORDER> (M, *this);
      }

			/**** INDEXING OPERATORS ****/
			
      /*! \brief Access or modify an element in this Matrix.
       *
       * This indexing operator allows the caller to access or modify
       * the ith (indexed in this Matrix's matrix_order) element of
       * this Matrix, indexed from 0 to n - 1, where n is the number
       * of elements in the Matrix.
       *
       * \param i The index of the element to access/modify.
       *
       * \see operator[](uint) const
       * \see operator()(uint)
       * \see operator()(uint) const
       * \see operator()(uint, uint)
       * \see operator()(uint, uint) const
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator[] (uint i)
			{
				SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
						"Index " << i << " out of range");

				return data_[Base::index(i)];
			}
			
      /*! \brief Access an element in this Matrix.
       *
       * This indexing operator allows the caller to access 
       * the ith (indexed in this Matrix's matrix_order) element of
       * this Matrix, indexed from 0 to n - 1, where n is the number
       * of elements in the Matrix.
       *
       * \param i The index of the element to access.
       *
       * \see operator[](uint)
       * \see operator()(uint)
       * \see operator()(uint) const
       * \see operator()(uint, uint)
       * \see operator()(uint, uint) const
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator[] (uint i) const
			{
				SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
						"Index " << i << " out of range");

				return data_[Base::index(i)];
			}

      /*! \brief Access or modify an element in this Matrix.
       *
       * This indexing operator allows the caller to access or modify
       * the ith (indexed in this Matrix's matrix_order) element of
       * this Matrix, indexed from 0 to n - 1, where n is the number
       * of elements in the Matrix.
       *
       * \param i The index of the element to access/modify.
       *
       * \see operator[](uint)
       * \see operator[](uint) const
       * \see operator()(uint) const
       * \see operator()(uint, uint)
       * \see operator()(uint, uint) const
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator() (uint i)
			{
				SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
						"Index " << i << " out of range");

				return data_[Base::index(i)];
			}
			
      /*! \brief Access an element in this Matrix.
       *
       * This indexing operator allows the caller to access 
       * the ith (indexed in this Matrix's matrix_order) element of
       * this Matrix, indexed from 0 to n - 1, where n is the number
       * of elements in the Matrix.
       *
       * \param i The index of the element to access.
       *
       * \see operator[](uint)
       * \see operator[](uint) const
       * \see operator()(uint)
       * \see operator()(uint, uint)
       * \see operator()(uint, uint) const
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator() (uint i) const
			{
				SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
					"Index " << i << " out of range");

				return data_[Base::index(i)];
			}

      /*! \brief Access or modify an element in this Matrix.
       *
       * This indexing operator allows the caller to access or modify
       * the (i, j)th element of
       * this Matrix, where i is an element of 0, 1, ..., rows - 1 and
       * j is an element of 0, 1, ..., columns - 1.
       *
       * \param i The row index of the element to access/modify.
       * \param j The column index of the element to access/modify.
       *
       * \see operator[](uint)
       * \see operator[](uint) const
       * \see operator()(uint)
       * \see operator()(uint) const
       * \see operator()(uint, uint) const
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator() (uint i, uint j)
			{
				SCYTHE_CHECK_30 (! Base::inRange(i, j), scythe_bounds_error,
						"Index (" << i << ", " << j << ") out of range");

				return data_[Base::index(i, j)];
			}
				
      /*! \brief Access an element in this Matrix.
       *
       * This indexing operator allows the caller to access 
       * the (i, j)th element of
       * this Matrix, where i is an element of 0, 1, ..., rows - 1 and
       * j is an element of 0, 1, ..., columns - 1.
       *
       * \param i The row index of the element to access.
       * \param j The column index of the element to access.
       *
       * \see operator[](uint)
       * \see operator[](uint) const
       * \see operator()(uint)
       * \see operator()(uint) const
       * \see operator() (uint, uint)
       * 
       * \throw scythe_bounds_error (Level 3)
       */
			inline T_type& operator() (uint i, uint j) const
			{
				SCYTHE_CHECK_30 (! Base::inRange(i, j), scythe_bounds_error,
						"Index (" << i << ", " << j << ") out of range");

				return data_[Base::index(i, j)];
			}

      /**** SUBMATRIX OPERATORS ****/


      /* Submatrices are always views.  An extra (but relatively
       * cheap) copy constructor call is made when mixing and matching
       * orders like
       *
       * Matrix<> A;
       * ...
       * Matrix<double, Row> B = A(2, 2, 4, 4);
       *
       * It is technically possible to get around this, by providing
       * templates of each function of the form
       * template <matrix_order O>
       * Matrix<T_type, O, View> operator() (...)
       *
       * but the syntax to call them (crappy return type inference):
       *
       * Matrix<double, Row> B = A.template operator()<Row>(2, 2, 4, 4)
       *
       * is such complete gibberish that I don't think this is worth
       * the slight optimization.
       */
      
      /*! \brief Returns a view of a submatrix.
       *
       * This operator returns a rectangular submatrix view of this
       * Matrix with its upper left corner at (x1, y1) and its lower
       * right corner at (x2, y2).
       *
       * \param x1 The upper row of the submatrix.
       * \param y1 The leftmost column of the submatrix.
       * \param x2 The lowest row of the submatrix.
       * \param y2 The rightmost column of the submatrix.
       *
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint)
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements)
       * \see operator()(uint, all_elements) const
       *
       * \throw scythe_bounds_error (Level 2)
       *
       * \b Example:
       * \include example.matrix.submatrix.cc
       */
      inline Matrix<T_type, ORDER, View> 
			operator() (uint x1, uint y1, uint x2, uint y2)
			{
				SCYTHE_CHECK_20 (! Base::inRange(x1, y1) 
            || ! Base::inRange(x2, y2)
						|| x1 > x2 || y1 > y2,
						scythe_bounds_error,
						"Submatrix (" << x1 << ", " << y1 << ") ; ("
						<< x2 << ", " << y2 << ") out of range or ill-formed");
				
				return (Matrix<T_type, ORDER, View>(*this, x1, y1, x2, y2));
			}
			
      /*! \brief Returns a view of a submatrix.
       *
       * This operator returns a rectangular submatrix view of this
       * Matrix with its upper left corner at (x1, y1) and its lower
       * right corner at (x2, y2).
       *
       * \param x1 The upper row of the submatrix.
       * \param y1 The leftmost column of the submatrix.
       * \param x2 The lowest row of the submatrix.
       * \param y2 The rightmost column of the submatrix.
       *
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(all_elements, uint)
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements)
       * \see operator()(uint, all_elements) const
       *
       * \throw scythe_bounds_error (Level 2)
       */
      inline Matrix<T_type, ORDER, View> 
      operator() (uint x1, uint y1, uint x2, uint y2) const
			{
				SCYTHE_CHECK_20 (! Base::inRange(x1, y1) 
            || ! Base::inRange(x2, y2)
						|| x1 > x2 || y1 > y2,
						scythe_bounds_error,
						"Submatrix (" << x1 << ", " << y1 << ") ; ("
						<< x2 << ", " << y2 << ") out of range or ill-formed");

				return (Matrix<T_type, ORDER, View>(*this, x1, y1, x2, y2));
			}

      /*! \brief Returns a view of a column vector.
       *
       * This operator returns a vector view of column j in this Matrix.
       *
       * \param a An all_elements object signifying whole vector access.
       * \param j The column to view.
       *
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements)
       * \see operator()(uint, all_elements) const
       *
       * \throw scythe_bounds_error (Level 2)
       *
       * \b Example:
       * \include example.matrix.vector.cc
       */
      inline Matrix<T_type, ORDER, View> 
			operator() (const all_elements a, uint j)
			{
				SCYTHE_CHECK_20 (j >= Base::cols(), scythe_bounds_error,
						"Column vector index " << j << " out of range");

				return (Matrix<T_type, ORDER, View>
           (*this, 0, j, Base::rows() - 1, j));
			}
			
      /*! \brief Returns a view of a column vector.
       *
       * This operator returns a vector view of column j in this Matrix.
       *
       * \param a An all_elements object signifying whole vector access.
       * \param j The column to view.
       *
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint)
       * \see operator()(uint, all_elements)
       * \see operator()(uint, all_elements) const
       *
       * \throw scythe_bounds_error (Level 2)
       */
      inline Matrix<T_type, ORDER, View> 
			operator() (const all_elements a, uint j) const
			{
				SCYTHE_CHECK_20 (j >= Base::cols(), scythe_bounds_error,
						"Column vector index " << j << " out of range");

				return (Matrix<T_type, ORDER, View>
           (*this, 0, j, Base::rows() - 1, j));
			}

      /*! \brief Returns a view of a row vector.
       *
       * This operator returns a vector view of row i in this Matrix.
       *
       * \param i The row to view.
       * \param b An all_elements object signifying whole vector access.
       *
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint)
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements) const
       *
       * \throw scythe_bounds_error (Level 2)
       *
       * \b Example:
       * \include example.matrix.vector.cc
       */
      inline Matrix<T_type, ORDER, View> 
			operator() (uint i, const all_elements b)
			{
				SCYTHE_CHECK_20 (i >= Base::rows(), scythe_bounds_error,
						"Row vector index " << i << " out of range");

				return (Matrix<T_type, ORDER, View>
            (*this, i, 0, i, Base::cols() - 1));
			}
			
      /*! \brief Returns a view of a row vector.
       *
       * This operator returns a vector view of row i in this Matrix.
       *
       * \param i The row to view.
       * \param b An all_elements object signifying whole vector access.
       *
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint)
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements)
       *
       * \throw s