matrix.h

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

		public:
			/**** ACCESSORS ****/

      /*! \brief Returns the total number of elements in the Matrix.
       *
       * \see rows()
       * \see cols()
       * \see max_size()
       */
			inline uint size () const
			{
				return (rows() * cols());
			}

			/*! \brief Returns the  number of rows in the Matrix.
       *
       * \see size()
       * \see cols()
       */
			inline uint rows() const
			{
				return rows_;
			}

			/*! \brief Returns the number of columns in the Matrix.
       *
       * \see size()
       * \see rows()
       */
			inline uint cols () const
			{
				return cols_;
			}

      /*! \brief Check matrix ordering.
       *
       * This method returns the matrix_order of this Matrix.
       *
       * \see style()
       * \see storeorder()
       */
      inline matrix_order order() const
      {
        return ORDER;
      }

      /*! \brief Check matrix style.
       *
       * This method returns the matrix_style of this Matrix.
       *
       * \see order()
       * \see storeorder()
       */
      inline matrix_style style() const
      {
        return STYLE;
      }

      /*! \brief Returns the storage order of the underlying
       * DataBlock.  
       *
       * In view matrices, the storage order of the data may not be the
       * same as the template ORDER.
       * 
       *
       * \see rowstride()
       * \see colstride()
       * \see order()
       * \see style()
       */
      inline matrix_order storeorder () const
      {
        return storeorder_;
      }

      /*! \brief Returns the in-memory distance between elements in
       * successive rows of the matrix.
       *
       * \see colstride()
       * \see storeorder()
			 */
			inline uint rowstride () const
			{
				return rowstride_;
			}
			
      /*! \brief Returns the in-memory distance between elements in
       * successive columns of the matrix.
       *
       * \see rowstride()
       * \see storeorder()
			 */
      inline uint colstride () const
			{
				return colstride_;
			}

      /*! \brief Returns the maximum possible matrix size.
       *
       * Maximum matrix size is simply the highest available unsigned
       * int on your system.
       *
       * \see size()
			 */
			inline uint max_size () const
			{
				return UINT_MAX;
			}

			/*! \brief Returns true if this Matrix is 1x1.
       *
       * \see isNull()
       */
			inline bool isScalar () const
			{
				return (rows() == 1 && cols() == 1);
			}

			/*! \brief Returns true if this Matrix is 1xm.
       *
       * \see isColVector()
       * \see isVector()
       */
			inline bool isRowVector () const
			{
				return (rows() == 1);
			}
			
			/*! \brief Returns true if this Matrix is nx1.
       *
       * \see isRowVector()
       * \see isVector()
       */
			inline bool isColVector () const
			{
				return (cols() == 1);
			}

			/*! \brief Returns true if this Matrix is nx1 or 1xn.
       *
       * \see isRowVector()
       * \see isColVector()
       */
			inline bool isVector () const
			{
				return (cols() == 1 || rows() == 1);
			}
			
			/*! \brief Returns true if this Matrix is nxn.
       *
       * Null and scalar matrices are both considered square.
       *
       * \see isNull()
       * \see isScalar()
       */
			inline bool isSquare () const
			{
				return (cols() == rows());
			}

      /*! \brief Returns true if this Matrix has 0 elements.
       *  
       *  \see empty()
       *  \see isScalar()
       */
			inline bool isNull () const
			{
				return (rows() == 0);
			}

      /*! \brief Returns true if this Matrix has 0 elements.
       *
       * This function is identical to isNull() but conforms to STL
       * container class conventions.
       *
       * \see isNull()
       */
			inline bool empty () const
			{
				return (rows() == 0);
			}
			

			/**** HELPERS ****/

			/*! \brief Check if an index is in bounds.
       *
       * This function takes a single-argument index into a Matrix and
       * returns true iff that index is within the bounds of the
       * Matrix.  This function is equivalent to the expression:
       * \code
       * i < M.size()
       * \endcode
       * for a given Matrix M.
       *
       * \param i The element index to check.
       *
       * \see inRange(uint, uint)
       */
			inline bool inRange (uint i) const
			{
				return (i < size());
			}

			/*! \brief Check if an index is in bounds.
       *
       * This function takes a two-argument index into a Matrix and
       * returns true iff that index is within the bounds of the
       * Matrix.  This function is equivalent to the expression:
       * \code
       * i < M.rows() && j < M.cols()
       * \endcode
       * for a given Matrix M.
       *
       * \param i The row value of the index to check.
       * \param j The column value of the index to check.
       *
       * \see inRange(uint)
       */
			inline bool inRange (uint i, uint j) const
			{
				return (i < rows() && j < cols());
			}

    protected:
			/* These methods take offsets into a matrix and convert them
			 * into that actual position in the referenced memory block,
			 * taking stride into account.  Protection is debatable.  They
			 * could be useful to outside functions in the library but most
			 * callers should rely on the public () operator in the derived
			 * class or iterators.
       *
       * Note that these are very fast for concrete matrices but not
       * so great for views.  Of course, the type checks are done at
       * compile time.
			 */
			
			/* Turn an index into a true offset into the data. */
			inline uint index (uint i) const
      {
        if (STYLE == View) {
          if (ORDER == Col) {
            uint col = i / rows();
            uint row = i % rows();
            return (index(row, col));
          } else {
            uint row = i / cols();
            uint col = i % cols();
            return (index(row, col));
          }
        } else
          return(i);
      }

			/* Turn an i, j into an index. */
      inline uint index (uint row, uint col) const
      {
        if (STYLE == Concrete) {
          if (ORDER == Col)
				    return (col * rows() + row);
          else
            return (row * cols() + col);
        } else { // view
          if (storeorder_ == Col)
            return (col * colstride() + row);
          else
            return (row * rowstride() + col);
        }
      }

    
    /**** INSTANCE VARIABLES ****/
    protected:
      uint rows_;   // # of rows
      uint cols_;   // # of cols

    private:
      /* The derived class shouldn't have to worry about this
       * implementation detail.
       */
      uint rowstride_;   // the in-memory number of elements from the
      uint colstride_;   // beginning of one column(row) to the next
      matrix_order storeorder_; // The in-memory storage order of this
                                // matrix.  This will always be the
                                // same as ORDER for concrete
                                // matrices but views can look at
                                // matrices with storage orders that
                                // differ from their own.
                                // TODO storeorder is always known at
                                // compile time, so we could probably
                                // add a third template param to deal
                                // with this.  That would speed up
                                // views a touch.  Bit messy maybe.
  };

	/**** MATRIX CLASS ****/

  /*!  \brief An STL-compliant matrix container class.
   *
   * The Matrix class provides a matrix object with an interface similar
   * to standard mathematical notation.  The class provides a number
   * of unary and binary operators for manipulating matrices.
   * Operators provide such functionality as addition, multiplication,
   * and access to specific elements within the Matrix.  One can test
   * two matrices for equality or use provided methods to test the
   * size, shape, or symmetry of a given Matrix.  In addition, we
   * provide a number of facilities for saving, loading, and printing
   * matrices.  Other portions of the library provide functions for
   * manipulating matrices.  Most notably, la.h provides definitions
   * of common linear algebra routines and ide.h defines functions
   * that perform inversion and decomposition.
   * 
   * This Matrix data structure sits at the core of the library.  In
   * addition to standard matrix operations, this class allows
   * multiple matrices to view and modify the same underlying data.
   * This ability provides an elegant way in which to access and
   * modify submatrices such as isolated row vectors and greatly
   * increases the overall flexibility of the class.  In addition, we
   * provide iterators (defined in matrix_random_access_iterator.h,
   * matrix_forward_iterator.h, and matrix_bidirectional_iterator.h)
   * that allow Matrix objects to interact seamlessly with the generic
   * algorithms provided by the Standard Template Library (STL).
   *
   * The Matrix class uses template parameters to define multiple
   * behaviors.  Matrices are templated on data type, matrix_order,
   * and matrix_style.
   *
   * Matrix objects can contain elements of any type.  For the most
   * part, uses will wish to fill their matrices with single or double
   * precision floating point numbers, but matrices of integers,
   * boolean values, complex numbers, and user-defined types are all
   * possible and useful.  Although the basic book-keeping methods in
   * the Matrix class will support virtually any type, certain
   * operators require that one or more mathematical operator be
   * defined for the given type and many of the functions in the wider
   * Scythe library expect, or even demand, matrices containing floating
   * point numbers.
   *
   * There are two possible Matrix element orderings, row- or
   * column-major.  Differences in matrix ordering will be most
   * noticeable at construction time.  Constructors that build matrices
   * from streams or other list-like structures will place elements
   * into the matrix in its given order.  In general, any method that
   * processes a matrix in order will use the given matrix_order.  For
   * the most part, matrices of both orderings should exhibit the same
   * performance, but when a trade-off must be made, we err on the
   * side of column-major ordering.  In one respect, this bias is very
   * strong.  If you enable LAPACK/BLAS support in with the
   * SCYTHE_LAPACK compiler flag, the library will use these optimized
   * fortran routines to perform a number of operations on column
   * major matrices; we provide no LAPACK/BLAS support for row-major
   * matrices.  Operations on matrices with mismatched ordering are
   * legal and supported, but not guaranteed to be as fast as
   * same-order operations, especially when SCYTHE_LAPACK is enabled.
   *
   * There are also two possible styles of Matrix template, concrete
   * and view.  These two types of matrix provide distinct ways in
   * which to interact with an underlying block of data.
   * 
   * Concrete matrices behave like matrices in previous
   * Scythe releases.  They directly encapsulate a block of data and
   * always process it in the same order as it is stored (their
   * matrix_order always matches the underlying storage order).
   * All copy constructions and assignments on
   * concrete matrices make deep copies and it is not possible to use
   * the reference() method to make a concrete Matrix a view of
   * another Matrix.  Furthermore, concrete matrices are guaranteed to
   * have unit stride (That is, adjacent Matrix elements are stored
   * adjacently in memory).  
   *
   * Views, on the other hand, provide references to data blocks.
   * More than one view can look at the same underlying block of data,
   * possibly at different portions of the data at the same time.
   * Furthermore, a view may look at the data block of a concrete
   * matrix, perhaps accessing a single row vector or a small
   * submatrix of a larger matrix.  When you copy construct