[Rcpp-commits] r3318 - in pkg/RcppEigen: inst/examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 10 19:47:58 CET 2011
Author: dmbates
Date: 2011-11-10 19:47:58 +0100 (Thu, 10 Nov 2011)
New Revision: 3318
Modified:
pkg/RcppEigen/inst/examples/lmBenchmark.R
pkg/RcppEigen/src/Makevars
pkg/RcppEigen/src/Makevars.win
pkg/RcppEigen/src/fastLm.cpp
pkg/RcppEigen/src/fastLm.h
Log:
Added another fastLm method using the Lapack subroutine dgesdd for the SVD
Modified: pkg/RcppEigen/inst/examples/lmBenchmark.R
===================================================================
--- pkg/RcppEigen/inst/examples/lmBenchmark.R 2011-11-09 21:54:11 UTC (rev 3317)
+++ pkg/RcppEigen/inst/examples/lmBenchmark.R 2011-11-10 18:47:58 UTC (rev 3318)
@@ -21,7 +21,9 @@
exprs$PivQR <- expression(.Call("fastLm", mm, y, 0L, PACKAGE="RcppEigen"))
## LDLt Cholesky decomposition with rank detection
exprs$LDLt <- expression(.Call("fastLm", mm, y, 2L, PACKAGE="RcppEigen"))
-## SVD
+## SVD using the Lapack subroutine dgesdd and Eigen support
+exprs$GESDD <- expression(.Call("fastLm", mm, y, 6L, PACKAGE="RcppEigen"))
+## SVD (the JacobiSVD class from Eigen)
exprs$SVD <- expression(.Call("fastLm", mm, y, 4L, PACKAGE="RcppEigen"))
## eigenvalues and eigenvectors of X'X
exprs$SymmEig <- expression(.Call("fastLm", mm, y, 5L, PACKAGE="RcppEigen"))
Modified: pkg/RcppEigen/src/Makevars
===================================================================
--- pkg/RcppEigen/src/Makevars 2011-11-09 21:54:11 UTC (rev 3317)
+++ pkg/RcppEigen/src/Makevars 2011-11-10 18:47:58 UTC (rev 3318)
@@ -7,5 +7,4 @@
## to enable error checking. However, checking for errors is flagged
## in R CMD check as an undesirable practice.
-PKG_LIBS=`$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"`
-
+PKG_LIBS=`$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
\ No newline at end of file
Modified: pkg/RcppEigen/src/Makevars.win
===================================================================
--- pkg/RcppEigen/src/Makevars.win 2011-11-09 21:54:11 UTC (rev 3317)
+++ pkg/RcppEigen/src/Makevars.win 2011-11-10 18:47:58 UTC (rev 3318)
@@ -1,6 +1,6 @@
## This assumes that we can call Rscript to ask Rcpp about its locations
## Use the R_HOME indirection to support installations of multiple R version
-PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
+PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
PKG_CPPFLAGS = -I../inst/include -I. -DNDEBUG
Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp 2011-11-09 21:54:11 UTC (rev 3317)
+++ pkg/RcppEigen/src/fastLm.cpp 2011-11-10 18:47:58 UTC (rev 3318)
@@ -21,6 +21,7 @@
// <http://www.gnu.org/licenses/>.
#include "fastLm.h"
+#include <R_ext/Lapack.h>
namespace lmsol {
using Rcpp::_;
@@ -54,18 +55,14 @@
return *this;
}
- ArrayXd lm::Dplus(const ArrayXd& D) {
- ArrayXd Di(m_p);
- double comp(D.maxCoeff() * threshold());
- for (int j = 0; j < m_p; ++j) Di[j] = (D[j] < comp) ? 0. : 1./D[j];
- m_r = (Di != 0.).count();
- return Di;
+ inline ArrayXd lm::Dplus(const ArrayXd& d) {
+ ArrayXd di(d.size());
+ double comp(d.maxCoeff() * threshold());
+ for (int j = 0; j < d.size(); ++j) di[j] = (d[j] < comp) ? 0. : 1./d[j];
+ m_r = (di != 0.).count();
+ return di;
}
- MatrixXd lm::I_p() const {
- return MatrixXd::Identity(m_p, m_p);
- }
-
MatrixXd lm::XtX() const {
return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
rankUpdate(m_X.adjoint());
@@ -138,6 +135,31 @@
m_se = Ch.solve(I_p()).diagonal().array().sqrt();
}
+ int gesdd(MatrixXd& A, ArrayXd& S, MatrixXd& Vt) {
+ int info, mone = -1, m = A.rows(), n = A.cols();
+ std::vector<int> iwork(8 * n);
+ double wrk;
+ if (m < n || S.size() != n || Vt.rows() != n || Vt.cols() != n)
+ throw std::invalid_argument("dimension mismatch in gesvd");
+ F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
+ &m, Vt.data(), &n, &wrk, &mone, &iwork[0], &info);
+ int lwork(wrk);
+ std::vector<double> work(lwork);
+ F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
+ &m, Vt.data(), &n, &work[0], &lwork, &iwork[0], &info);
+ return info;
+ }
+
+ GESDD::GESDD(const Map<MatrixXd>& X, const Map<VectorXd> &y) : lm(X, y) {
+ MatrixXd U(X), Vt(m_p, m_p);
+ ArrayXd S(m_p);
+ if (gesdd(U, S, Vt)) throw std::runtime_error("error in gesdd");
+ MatrixXd VDi(Vt.adjoint() * Dplus(S).matrix().asDiagonal());
+ m_coef = VDi * U.adjoint() * y;
+ m_fitted = X * m_coef;
+ m_se = VDi.rowwise().norm();
+ }
+
SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
JacobiSVD<MatrixXd> UDV(X.jacobiSvd(ComputeThinU|ComputeThinV));
MatrixXd VDi(UDV.matrixV() *
@@ -157,7 +179,7 @@
m_se = VDi.rowwise().norm();
}
- enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t};
+ enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t, GESDD_t};
static inline lm do_lm(const Map<MatrixXd> &X, const Map<VectorXd> &y, int type) {
switch(type) {
@@ -173,6 +195,8 @@
return SVD(X, y);
case SymmEigen_t:
return SymmEigen(X, y);
+ case GESDD_t:
+ return GESDD(X, y);
}
throw invalid_argument("invalid type");
return ColPivQR(X, y); // -Wall
Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h 2011-11-09 21:54:11 UTC (rev 3317)
+++ pkg/RcppEigen/src/fastLm.h 2011-11-10 18:47:58 UTC (rev 3318)
@@ -29,8 +29,6 @@
using Eigen::ColPivHouseholderQR;
using Eigen::ComputeThinU;
using Eigen::ComputeThinV;
- using Eigen::DiagonalMatrix;
- using Eigen::Dynamic;
using Eigen::HouseholderQR;
using Eigen::JacobiSVD;
using Eigen::LDLT;
@@ -66,7 +64,7 @@
lm(const Map<MatrixXd>&, const Map<VectorXd>&);
ArrayXd Dplus(const ArrayXd& D);
- MatrixXd I_p() const;
+ MatrixXd I_p() const {return MatrixXd::Identity(m_p, m_p);}
MatrixXd XtX() const;
// setThreshold and threshold are based on ColPivHouseholderQR methods
@@ -98,6 +96,11 @@
QR(const Map<MatrixXd>&, const Map<VectorXd>&);
};
+ class GESDD : public lm {
+ public:
+ GESDD(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
+
class SVD : public lm {
public:
SVD(const Map<MatrixXd>&, const Map<VectorXd>&);
More information about the Rcpp-commits
mailing list