ide.h

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

    // template declaration of the more general template.
    Matrix<> AA = A;

    // Get a pointer to the internal array and set up some vars
    double* Aarray = AA.getArray();  // internal array pointer
    int rows = (int) AA.rows(); // the dim of the matrix
    int err = 0; // The output error condition

    // Cholesky decomposition step
    lapack::dpotrf_("L", &rows, Aarray, &rows, &err);
    SCYTHE_CHECK_10(err > 0, scythe_type_error,
        "Matrix is not positive definite")
    SCYTHE_CHECK_10(err < 0, scythe_invalid_arg,
        "The " << err << "th value of the matrix had an illegal value")

    // Inversion step
    lapack::dpotri_("L", &rows, Aarray, &rows, &err);
    SCYTHE_CHECK_10(err > 0, scythe_type_error,
        "The (" << err << ", " << err << ") element of the matrix is zero"
        << " and the inverse could not be computed")
    SCYTHE_CHECK_10(err < 0, scythe_invalid_arg,
        "The " << err << "th value of the matrix had an illegal value")
    lapack::make_symmetric(Aarray, rows);

    return AA;
  }

  template<>
  inline Matrix<>
  invpd (const Matrix<>& A, const Matrix<>& M)
  {
    SCYTHE_DEBUG_MSG("Using lapack/blas for invpd");
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10 (! A.isSquare(), scythe_dimension_error,
        "A is not square");
    SCYTHE_CHECK_10(A.rows() != M.cols() || A.cols() != M.rows(), 
        scythe_conformation_error, "A and M do not conform");

    // We have to do an explicit copy within the func to match the
    // template declaration of the more general template.
    Matrix<> MM = M;

    // Get pointer and set up some vars
    double* Marray = MM.getArray();
    int rows = (int) MM.rows();
    int err = 0;

    // Inversion step
    lapack::dpotri_("L", &rows, Marray, &rows, &err);
    SCYTHE_CHECK_10(err > 0, scythe_type_error,
        "The (" << err << ", " << err << ") element of the matrix is zero"
        << " and the inverse could not be computed")
    SCYTHE_CHECK_10(err < 0, scythe_invalid_arg,
        "The " << err << "th value of the matrix had an illegal value")
    lapack::make_symmetric(Marray, rows);

    return MM;
  }

  template <>
  inline Matrix<>
  inv(const Matrix<>& A)
  {
    SCYTHE_DEBUG_MSG("Using lapack/blas for inv");
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10 (! A.isSquare(), scythe_dimension_error,
        "A is not square");

    // We have to do an explicit copy within the func to match the
    // template declaration of the more general template.
    Matrix<> AA = A;

    // Get a pointer to the internal array and set up some vars
    double* Aarray = AA.getArray();  // internal array pointer
    int rows = (int) AA.rows(); // the dim of the matrix
    int* ipiv = new int[rows];  // Holds the lu decomp pivot array
    int err = 0; // The output error condition

    // LU decomposition step
    lapack::dgetrf_(&rows, &rows, Aarray, &rows, ipiv, &err);

    SCYTHE_CHECK_10(err > 0, scythe_type_error, "Matrix is singular");
    SCYTHE_CHECK_10(err < 0, scythe_invalid_arg, 
        "The " << err << "th value of the matrix had an illegal value");

    // Inversion step; first do a workspace query, then the actual
    // inversion
    double work_query = 0;
    int work_size = -1;
    lapack::dgetri_(&rows, Aarray, &rows, ipiv, &work_query, 
                    &work_size, &err);
    double* workspace = new double[(work_size = (int) work_query)];
    lapack::dgetri_(&rows, Aarray, &rows, ipiv, workspace, &work_size,
                    &err);
    delete[] ipiv;
    delete[] workspace;

    SCYTHE_CHECK_10(err > 0, scythe_type_error, "Matrix is singular");
    SCYTHE_CHECK_10(err < 0, scythe_invalid_arg, 
        "Internal error in LAPACK routine dgetri");

    return AA;
  }

  /*!\brief The result of a singular value decomposition.
   *
   * Objects of this type hold the results of a singular value
   * decomposition (SVD) of an \f$m \times n\f$ matrix \f$A\f$, as
   * returned by svd().  The SVD takes the form: \f$A = U
   * \cdot \Sigma \cdot V'\f$.  SVD objects contain \a d, which
   * holds the singular values of \f$A\f$ (the diagonal of
   * \f$\Sigma\f$) in descending order.  Furthermore, depending on the
   * options passed to svd(), they may hold some or all of the
   * left singular vectors of \f$A\f$ in \a U and some or all of the
   * right singular vectors of \f$A\f$ in \a Vt.
   *
   * \see svd(const Matrix<>& A, int nu, int nv);
   */
   
  struct SVD {
    Matrix<> d;  // singular values
    Matrix<> U;  // left singular vectors
    Matrix<> Vt; // transpose of right singular vectors
  };

 /*!\brief Calculates the singular value decomposition of a matrix,
  * optionally computing the left and right singular vectors.
  *
  * This function returns the singular value decomposition (SVD) of a
  * \f$m \times n\f$ matrix \a A, optionally computing the left and right
  * singular vectors.  It returns the singular values and vectors in
  * a SVD object.
  *
  * \note This function requires BLAS/LAPACK functionality and is
  * only available on machines that provide these libraries.  Make
  * sure you enable the SCYTHE_LAPACK preprocessor flag if you wish
  * to use this function.  Furthermore, note that this function takes
  * and returns only column-major concrete matrices.  Future versions
  * of Scythe will provide a native C++ implementation of this
  * function with support for general matrix templates.
  *
  * \param A The matrix to decompose.
  * \param nu The number of left singular vectors to compute and return.  Values less than zero are equivalent to min(\f$m\f$, \f$n\f$).
  * \param nv The number of right singular vectors to compute and return.  Values less than zero are equivalent to min(\f$m\f$, \f$n\f$).
  *
  * \throw scythe_null_error (Level 1)
  * \throw scythe_convergence_error (Level 1)
  * \throw scythe_lapack_internal_error (Level 1)
  *
  * \see SVD
  * \see eigen(const Matrix<>& A, bool vectors)
  */
 
 inline SVD
 svd (const Matrix<>& A, int nu = -1, int nv = -1)
 {
   SCYTHE_DEBUG_MSG("Using lapack/blas for eigen");
   SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
       "Matrix is NULL");

   char* jobz;
   int m = (int) A.rows();
   int n = (int) A.cols();
   int mn = (int) std::min(A.rows(), A.cols());
   Matrix<> U;
   Matrix<> V;
   if (nu < 0) nu = mn;
   if (nv < 0) nv = mn;
   if (nu <= mn && nv<= mn) {
     jobz = "S";
     U = Matrix<>(m, mn, false);
     V = Matrix<>(mn, n, false);
   } else if (nu == 0 && nv == 0) {
     jobz = "N";
   } else {
     jobz = "A";
     U = Matrix<>(m, m, false);
     V = Matrix<>(n, n, false);
   }
   double* Uarray = U.getArray();
   double* Varray = V.getArray();

   int ldu = (int) U.rows();
   int ldvt = (int) V.rows();
   Matrix<> X = A;
   double* Xarray = X.getArray();
   Matrix<> d(mn, 1, false);
   double* darray = d.getArray();

   double tmp, *work;
   int lwork, info;
   int *iwork = new int[8 * mn];

   // get optimal workspace
   lwork = -1;
   lapack::dgesdd_(jobz, &m, &n, Xarray, &m, darray, Uarray, &ldu,
                   Varray, &ldvt, &tmp, &lwork, iwork, &info);
   SCYTHE_CHECK_10(info < 0, scythe_lapack_internal_error,
       "Internal error in LAPACK routine dgessd");
   SCYTHE_CHECK_10(info > 0, scythe_convergence_error, "Did not converge");

   lwork = (int) tmp;
   work = new double[lwork];

   // Now for real
   lapack::dgesdd_(jobz, &m, &n, Xarray, &m, darray, Uarray, &ldu,
                   Varray, &ldvt, work, &lwork, iwork, &info);
   SCYTHE_CHECK_10(info < 0, scythe_lapack_internal_error,
       "Internal error in LAPACK routine dgessd");
   SCYTHE_CHECK_10(info > 0, scythe_convergence_error, "Did not converge");
   delete[] work;
  
   if (nu < mn && nu > 0)
     U = U(0, 0, U.rows() - 1, (unsigned int) std::min(m, nu) - 1);
   if (nv < mn && nv > 0)
     V = V(0, 0, (unsigned int) std::min(n, nv) - 1, V.cols() - 1);
   SVD result;
   result.d = d;
   result.U = U;
   result.Vt = V;

   return result;
 }

 /*!\brief The result of an eigenvalue/vector decomposition.
  *
  * Objects of this type hold the results of the eigen() function.
  * That is the eigenvalues and, optionally, the eigenvectors of a
  * symmetric matrix of order \f$n\f$.  The eigenvalues are stored in
  * ascending order in the member column vector \a values.  The
  * vectors are stored in the \f$n \times n\f$ matrix \a vectors.
  *
  * \see eigen(const Matrix<>& A, bool vectors)
  */
 
 struct Eigen {
   Matrix<> values;
   Matrix<> vectors;
 };

 /*!\brief Calculates the eigenvalues and eigenvectors of a symmetric
  * matrix.
  *
  * This function returns the eigenvalues and, optionally,
  * eigenvectors of a symmetric matrix \a A of order \f$n\f$.  It
  * returns an Eigen object containing the vector of values, in
  * ascending order, and, optionally, a matrix holding the vectors.
  *
  * \note This function requires BLAS/LAPACK functionality and is
  * only available on machines that provide these libraries.  Make
  * sure you enable the SCYTHE_LAPACK preprocessor flag if you wish
  * to use this function.  Furthermore, note that this function takes
  * and returns only column-major concrete matrices.  Future versions
  * of Scythe will provide a native C++ implementation of this
  * function with support for general matrix templates.
  *
  * \param A The Matrix to be decomposed.
  * \param vectors This boolean value indicates whether or not to
  * return eigenvectors in addition to eigenvalues.  It is set to true
  * by default.
  *
  * \throw scythe_null_error (Level 1)
  * \throw scythe_dimension_error (Level 1)
  * \throw scythe_lapack_internal_error (Level 1)
  *
  * \see Eigen
  * \see svd(const Matrix<>& A, int nu, int nv);
  */
  inline Eigen
  eigen (const Matrix<>& A, bool vectors=true)
  {
    SCYTHE_DEBUG_MSG("Using lapack/blas for eigen");
    SCYTHE_CHECK_10(! A.isSquare(), scythe_dimension_error,
        "Matrix not square");
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "Matrix is NULL");
    // Should be symmetric but rounding errors make checking for this
    // difficult.

    // Make a copy of A
    Matrix<> AA = A;

    // Get a point to the internal array and set up some vars
    double* Aarray = AA.getArray(); // internal array points
    int order = (int) AA.rows();    // input matrix is order x order
    double dignored = 0;            // we don't use this option
    int iignored = 0;               // or this one
    double abstol = 0.0;            // tolerance (default)
    int m;                          // output value
    Matrix<> result;                // result matrix
    char getvecs[1];                // are we getting eigenvectors?
    if (vectors) {
      getvecs[0] = 'V';
      result = Matrix<>(order, order + 1, false);
    } else {
      result = Matrix<>(order, 1, false);
      getvecs[0] = 'N';
    }
    double* eigenvalues = result.getArray(); // pointer to result array
    int* isuppz = new int[2 * order];        // indices of nonzero eigvecs
    double tmp;   // inital temporary value for getting work-space info
    int lwork, liwork, *iwork, itmp; // stuff for workspace
    double *work; // and more stuff for workspace
    int info = 0;  // error code holder

    // get optimal size for work arrays
    lwork = -1;
    liwork = -1;
    lapack::dsyevr_(getvecs, "A", "L", &order, Aarray, &order, &dignored,
        &dignored, &iignored, &iignored, &abstol, &m, eigenvalues, 
        eigenvalues + order, &order, isuppz, &tmp, &lwork, &itmp,
        &liwork, &info);
    SCYTHE_CHECK_10(info != 0, scythe_lapack_internal_error,
        "Internal error in LAPACK routine dsyevr");
    lwork = (int) tmp;
    liwork = itmp;
    work = new double[lwork];
    iwork = new int[liwork];

    // do the actual operation
    lapack::dsyevr_(getvecs, "A", "L", &order, Aarray, &order, &dignored,
        &dignored, &iignored, &iignored, &abstol, &m, eigenvalues, 
        eigenvalues + order, &order, isuppz, work, &lwork, iwork,
        &liwork, &info);
    SCYTHE_CHECK_10(info != 0, scythe_lapack_internal_error,
        "Internal error in LAPACK routine dsyevr");

    delete[] isuppz;
    delete[] work;
    delete[] iwork;

    Eigen resobj;
    if (vectors) {
      resobj.values = result(_, 0);
      resobj.vectors = result(0, 1, result.rows() -1, result.cols() - 1);
    } else {
      resobj.values = result;
    }

    return resobj;
  }

#endif

} // end namespace scythe

#endif /* SCYTHE_IDE_H */