[Rcpp-commits] r3230 - pkg/RcppEigen/src

[noreply at r-forge.r-project.org]

Author: dmbates
Date: 2011-10-26 22:16:52 +0200 (Wed, 26 Oct 2011)
New Revision: 3230

Modified:
   pkg/RcppEigen/src/fastLm.cpp
   pkg/RcppEigen/src/fastLm.h
Log:
Use a method for XtX but returning a MatrixXd, not a view.


Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp	2011-10-26 19:01:46 UTC (rev 3229)
+++ pkg/RcppEigen/src/fastLm.cpp	2011-10-26 20:16:52 UTC (rev 3230)
@@ -55,13 +55,11 @@
 	return *this;
     }
 
-    SelfAdjointView<MatrixXd,Lower> lm::XtX() const {
+    MatrixXd lm::XtX() const {
 	return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
 	    rankUpdate(m_X.adjoint());
     }
-// For some reason that function returning a SelfAdjointView encounters a bad_alloc error
-// Use a macro for the time being
-#define XtX MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().rankUpdate(m_X.adjoint())
+
     /** Returns the threshold that will be used by certain methods such as rank().
      * 
      *  The default value comes from experimenting (see "LU precision
@@ -112,7 +110,7 @@
     
     
     Llt::Llt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	LLT<MatrixXd>  Ch(XtX);
+	LLT<MatrixXd>  Ch(XtX().selfadjointView<Lower>());
 	m_coef            = Ch.solve(X.adjoint() * y);
 	m_fitted          = X * m_coef;
 	m_se              = Ch.matrixL().solve(I_p()).colwise().norm();
@@ -127,7 +125,7 @@
     }
 
     Ldlt::Ldlt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	LDLT<MatrixXd> Ch(XtX);
+	LDLT<MatrixXd> Ch(XtX().selfadjointView<Lower>());
 	ArrayXd         D(Ch.vectorD());
 	m_r               = (D > D.maxCoeff() * threshold()).count();
 	// FIXME: work out how to use Dplus with elements of D unsorted.
@@ -149,7 +147,7 @@
 
     SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y)
 	: lm(X, y) {
-	SelfAdjointEigenSolver<MatrixXd> eig(XtX);
+	SelfAdjointEigenSolver<MatrixXd> eig(XtX().selfadjointView<Lower>());
 	ArrayXd                      D(eig.eigenvalues());
 	m_r                            = (D > D[m_p - 1] * threshold()).count();
 	D                              = D.sqrt();

Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h	2011-10-26 19:01:46 UTC (rev 3229)
+++ pkg/RcppEigen/src/fastLm.h	2011-10-26 20:16:52 UTC (rev 3230)
@@ -71,7 +71,7 @@
 	    return MatrixXd::Identity(m_p, m_p);
 	}
 	RealScalar                threshold() const;
-	SelfAdjointView<MatrixXd,Lower> XtX() const;
+	MatrixXd                        XtX() const;
 	const VectorXd&                  se() const {return m_se;}
 	const VectorXd&                coef() const {return m_coef;}
 	const VectorXd&              fitted() const {return m_fitted;}