[Rcpp-commits] r3163 - in pkg/RcppEigen/inst: include/unsupported/Eigen/src/SparseExtra unitTests

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

Author: dmbates
Date: 2011-07-29 04:46:53 +0200 (Fri, 29 Jul 2011)
New Revision: 3163

Modified:
   pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
   pkg/RcppEigen/inst/unitTests/runit.sparse.R
Log:
Allow rectangular matrix in CholmodDecomposition's factorize method.  Test rectangular version of KNex example.


Modified: pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
===================================================================
--- pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h	2011-07-28 22:55:48 UTC (rev 3162)
+++ pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h	2011-07-29 02:46:53 UTC (rev 3163)
@@ -311,8 +311,11 @@
 	{
 	    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
 	    
-	    cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
-	    if (M_cholmod_factorize(&A, m_cholmodFactor, &m_cholmod) != CHOLMOD_OK)
+	    cholmod_sparse    A = (matrix.rows() == matrix.cols()) ?
+		viewAsCholmod(matrix.template selfadjointView<UpLo>()) :
+		viewAsCholmod(matrix);
+	    
+	    if (!M_cholmod_factorize(&A, m_cholmodFactor, &m_cholmod))
 		throw std::runtime_error("CholmodDecomposition.factorize() error");
 	    
 	    this->m_info = Success;

Modified: pkg/RcppEigen/inst/unitTests/runit.sparse.R
===================================================================
--- pkg/RcppEigen/inst/unitTests/runit.sparse.R	2011-07-28 22:55:48 UTC (rev 3162)
+++ pkg/RcppEigen/inst/unitTests/runit.sparse.R	2011-07-29 02:46:53 UTC (rev 3163)
@@ -57,26 +57,58 @@
     using Eigen::CholmodAuto;
     using Eigen::Success;
 
-    List input(input_) ;
-    const SparseMatrix<double>      m1 = input[0];
-    const Map<VectorXd>             v1 = input[1];
-    SparseMatrix<double>            m2(m1.cols(), m1.cols());
-
+    List input(input_);
+    const MappedSparseMatrix<double> m1 = input[0];
+    const Map<VectorXd>         v1 = input[1];
+    SparseMatrix<double>        m2(m1.cols(), m1.cols());
     m2.selfadjointView<Lower>().rankUpdate(m1.adjoint());
 
-    CholmodDecomposition<SparseMatrix<double> > ff;
-    ff.setMode(CholmodAuto);
-    ff.compute(m2);
-    VectorXd                       res = ff.solve(m1.adjoint() * v1);
+    CholmodDecomposition<SparseMatrix<double> > ff(m2);
+    VectorXd                   res = ff.solve(m1.adjoint() * v1);
     
     return List::create(_["res"]   = res,
                         _["rows"]  = ff.rows(),
                         _["cols"]  = ff.cols());
 ',
-                      plugin = "RcppEigen",
-                      inc='#include "RcppEigenStubs.cpp"')
+                      plugin = "RcppEigen")
+
     rr <- fx(KNex)
     checkEquals(rr[[1]], as.vector(solve(crossprod(KNex[[1]]),
                                          crossprod(KNex[[1]], KNex[[2]]))),
                                        "Cholmod solution")
 }
+
+test.solveCholmodRect.R <- function() {
+    suppressMessages(require("Matrix", character.only=TRUE))
+    data("KNex", package = "Matrix")
+
+    fx <- cxxfunction( signature(input_ = "list"), '
+    using Eigen::VectorXd;
+    using Eigen::MatrixXd;
+    using Eigen::Lower;
+    using Eigen::Map;
+    using Eigen::MappedSparseMatrix;
+    using Eigen::SparseMatrix;
+    using Eigen::CholmodDecomposition;
+    using Eigen::CholmodAuto;
+    using Eigen::Success;
+
+    List input(input_);
+    const MappedSparseMatrix<double> m1 = input[0];
+    const Map<VectorXd>         v1 = input[1];
+    SparseMatrix<double>        m2(m1.adjoint());
+
+    CholmodDecomposition<SparseMatrix<double> > ff(m2);
+    VectorXd                   res = ff.solve(m2 * v1);
+    
+    return List::create(_["res"]   = res,
+                        _["rows"]  = ff.rows(),
+                        _["cols"]  = ff.cols());
+',
+                      plugin = "RcppEigen")
+
+    rr <- fx(KNex)
+    checkEquals(rr[[1]], as.vector(solve(crossprod(KNex[[1]]),
+                                         crossprod(KNex[[1]], KNex[[2]]))),
+                                       "Cholmod solution")
+}