rewrite spchol using eigen. 64/11464/8
Cedric Delamarre [Fri, 10 May 2013 07:23:34 +0000 (09:23 +0200)]
test_run("sparse","bug_6827",["no_check_error_output" ]);
test_run("arnoldi","bug_12137",["no_check_error_output" ]);

Change-Id: Ibb7c961aad9be62c6550e54a3ce62f31dada4a92

scilab/modules/ast/includes/types/sparse.hxx
scilab/modules/ast/src/cpp/types/sparse.cpp
scilab/modules/sparse/Makefile.am
scilab/modules/sparse/Makefile.in
scilab/modules/sparse/includes/sparse_gw.hxx
scilab/modules/sparse/sci_gateway/cpp/sci_spchol.cpp [new file with mode: 0644]
scilab/modules/sparse/sci_gateway/cpp/sparse_gw.cpp
scilab/modules/sparse/sci_gateway/cpp/sparse_gw.vcxproj
scilab/modules/sparse/sci_gateway/cpp/sparse_gw.vcxproj.filters
scilab/modules/sparse/tests/nonreg_tests/bug_6827.dia.ref
scilab/modules/sparse/tests/nonreg_tests/bug_6827.tst

index 6e00387..735a4a6 100644 (file)
@@ -347,6 +347,8 @@ struct EXTERN_AST Sparse : GenericType
         return true;
     }
 
+    int newCholLLT(Sparse** permut, Sparse** factor) const;
+
     /** create a new sparse matrix containing the non zero values set to 1.
        equivalent but faster than calling one_set() on a new copy of the
        current matrix.
index ddf4238..aff05e5 100644 (file)
 #include <iterator>
 #include <algorithm>
 
+#include <Eigen/Core>
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/SparseCholesky>
+
 #include "sparse.hxx"
 #include "types.hxx"
 #include "tostring_common.hxx"
@@ -1871,6 +1875,38 @@ Sparse* Sparse::dotDivide(Sparse SPARSE_CONST& o) const
     return new Sparse(realSp, cplxSp);
 }
 
+int Sparse::newCholLLT(Sparse** _SpPermut, Sparse** _SpFactor) const
+{
+    typedef Eigen::SparseMatrix<double, Eigen::ColMajor> RealSparseCol_t;
+    RealSparseCol_t spColMajor = RealSparseCol_t((const RealSparse_t&) * matrixReal);
+
+    // Constructs and performs the LLT factorization of sparse
+    Eigen::SimplicialLLT<RealSparseCol_t> pLLT(spColMajor);
+    int iInfo = pLLT.info();
+    if (iInfo != Eigen::Success)
+    {
+        *_SpFactor = NULL;
+        *_SpPermut = NULL;
+        return iInfo;
+    }
+
+    // Get the lower matrix of factorization.
+    // The new RealSparse_t will be setted in Sparse without copy.
+    *_SpFactor = new Sparse(new RealSparse_t(pLLT.matrixL()), NULL);
+
+    // Get the permutation matrix.
+    Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> p = pLLT.permutationP();
+    *_SpPermut = new Sparse(p.rows(), p.cols());
+    for (int i = 0; i < p.rows(); i++)
+    {
+        (*_SpPermut)->set(i, p.indices()[i], 1, false);
+    }
+
+    (*_SpPermut)->finalize();
+
+    return iInfo;
+}
+
 struct BoolCast
 {
     BoolCast(std::complex<double> const& c): b(c.real() || c.imag()) {}
@@ -1992,7 +2028,7 @@ template<typename S> struct GetReal: std::unary_function<typename S::InnerIterat
     }
 };
 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
-    : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
+        : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
 {
     double operator()( Sparse::CplxSparse_t::InnerIterator it) const
     {
index 7d3a6f9..c01f10d 100644 (file)
@@ -34,8 +34,8 @@ GATEWAY_CPP_SOURCES = \
     sci_gateway/cpp/sci_lusolve.cpp \
     sci_gateway/cpp/sci_luget.cpp \
     sci_gateway/cpp/sci_ludel.cpp \
-    sci_gateway/cpp/sci_ordmmd.cpp
-
+    sci_gateway/cpp/sci_ordmmd.cpp \
+    sci_gateway/cpp/sci_spchol.cpp
 
 libscisparse_la_CPPFLAGS = \
     -I$(srcdir)/src/c/ \
index 562b721..dd89ce2 100644 (file)
@@ -197,7 +197,8 @@ am__objects_3 = sci_gateway/cpp/libscisparse_la-sparse_gw.lo \
        sci_gateway/cpp/libscisparse_la-sci_lusolve.lo \
        sci_gateway/cpp/libscisparse_la-sci_luget.lo \
        sci_gateway/cpp/libscisparse_la-sci_ludel.lo \
-       sci_gateway/cpp/libscisparse_la-sci_ordmmd.lo
+       sci_gateway/cpp/libscisparse_la-sci_ordmmd.lo \
+       sci_gateway/cpp/libscisparse_la-sci_spchol.lo
 am_libscisparse_la_OBJECTS = $(am__objects_3)
 libscisparse_la_OBJECTS = $(am_libscisparse_la_OBJECTS)
 @MAINTAINER_MODE_FALSE@am_libscisparse_la_rpath =
@@ -614,7 +615,8 @@ GATEWAY_CPP_SOURCES = \
     sci_gateway/cpp/sci_lusolve.cpp \
     sci_gateway/cpp/sci_luget.cpp \
     sci_gateway/cpp/sci_ludel.cpp \
-    sci_gateway/cpp/sci_ordmmd.cpp
+    sci_gateway/cpp/sci_ordmmd.cpp \
+    sci_gateway/cpp/sci_spchol.cpp
 
 libscisparse_la_CPPFLAGS = \
     -I$(srcdir)/src/c/ \
@@ -912,6 +914,9 @@ sci_gateway/cpp/libscisparse_la-sci_ludel.lo:  \
 sci_gateway/cpp/libscisparse_la-sci_ordmmd.lo:  \
        sci_gateway/cpp/$(am__dirstamp) \
        sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
+sci_gateway/cpp/libscisparse_la-sci_spchol.lo:  \
+       sci_gateway/cpp/$(am__dirstamp) \
+       sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
 
 libscisparse.la: $(libscisparse_la_OBJECTS) $(libscisparse_la_DEPENDENCIES) $(EXTRA_libscisparse_la_DEPENDENCIES) 
        $(AM_V_CXXLD)$(CXXLINK) $(am_libscisparse_la_rpath) $(libscisparse_la_OBJECTS) $(libscisparse_la_LIBADD) $(LIBS)
@@ -938,6 +943,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_ordmmd.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_sp2adj.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_sparse.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spchol.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spcompack.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spget.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spones.Plo@am__quote@
@@ -1153,6 +1159,13 @@ sci_gateway/cpp/libscisparse_la-sci_ordmmd.lo: sci_gateway/cpp/sci_ordmmd.cpp
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscisparse_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o sci_gateway/cpp/libscisparse_la-sci_ordmmd.lo `test -f 'sci_gateway/cpp/sci_ordmmd.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_ordmmd.cpp
 
+sci_gateway/cpp/libscisparse_la-sci_spchol.lo: sci_gateway/cpp/sci_spchol.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscisparse_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT sci_gateway/cpp/libscisparse_la-sci_spchol.lo -MD -MP -MF sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spchol.Tpo -c -o sci_gateway/cpp/libscisparse_la-sci_spchol.lo `test -f 'sci_gateway/cpp/sci_spchol.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_spchol.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spchol.Tpo sci_gateway/cpp/$(DEPDIR)/libscisparse_la-sci_spchol.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='sci_gateway/cpp/sci_spchol.cpp' object='sci_gateway/cpp/libscisparse_la-sci_spchol.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscisparse_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o sci_gateway/cpp/libscisparse_la-sci_spchol.lo `test -f 'sci_gateway/cpp/sci_spchol.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_spchol.cpp
+
 .f.o:
        $(AM_V_F77)$(F77COMPILE) -c -o $@ $<
 
index f156f6f..9793a7a 100644 (file)
@@ -46,5 +46,6 @@ CPP_GATEWAY_PROTOTYPE(sci_lufact);
 CPP_GATEWAY_PROTOTYPE(sci_lusolve);
 CPP_GATEWAY_PROTOTYPE(sci_luget);
 CPP_GATEWAY_PROTOTYPE(sci_ludel);
+CPP_GATEWAY_PROTOTYPE(sci_spchol);
 
 #endif /* !__SPARSE_GW_HXX__ */
diff --git a/scilab/modules/sparse/sci_gateway/cpp/sci_spchol.cpp b/scilab/modules/sparse/sci_gateway/cpp/sci_spchol.cpp
new file mode 100644 (file)
index 0000000..50e04ec
--- /dev/null
@@ -0,0 +1,85 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2013 - Scilab Enterprises - Cedric Delamarre
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include "sparse_gw.hxx"
+#include "function.hxx"
+#include "sparse.hxx"
+
+extern "C"
+{
+#include "Scierror.h"
+#include "localization.h"
+#include "sciprint.h"
+}
+
+types::Function::ReturnValue sci_spchol(types::typed_list &in, int _iRetCount, types::typed_list &out)
+{
+    types::Sparse* pSpIn       = NULL;
+    types::Sparse* pSpPermut   = NULL;
+    types::Sparse* pSpFact     = NULL;
+
+    if (in.size() != 1)
+    {
+        Scierror(999, _("%s: Wrong number of input argument(s): %d expected.\n"), "spchol", 1);
+        return types::Function::Error;
+    }
+
+    if (_iRetCount != 2)
+    {
+        Scierror(999, _("%s: Wrong number of output argument(s): %d expected.\n"), "spchol", 2);
+        return types::Function::Error;
+    }
+
+    if (in[0]->isSparse() == false && in[0]->isSparseBool() == false)
+    {
+        Scierror(999, _("%s: Wrong type for argument #%d: Sparse matrix expected.\n"), "spchol", 1);
+        return types::Function::Error;
+    }
+
+    pSpIn = in[0]->getAs<types::Sparse>();
+
+    if (pSpIn->isComplex())
+    {
+        Scierror(999, _("%s: Wrong type for argument #%d: Real matrix expected.\n"), "spchol", 1);
+        return types::Function::Error;
+    }
+
+    if (pSpIn->getRows() != pSpIn->getCols())
+    {
+        Scierror(999, _("%s: Wrong size for argument #%d: Square sparse matrix expected.\n"), "spchol", 1);
+        return types::Function::Error;
+    }
+
+    int ierr = pSpIn->newCholLLT(&pSpPermut, &pSpFact);
+
+    switch (ierr)
+    {
+        case Eigen::Success :
+            break;
+        case Eigen::NumericalIssue :
+            Scierror(999, _("%s: The provided data did not satisfy the prerequisites.\n"), "spchol");
+            return types::Function::Error;
+        case Eigen::NoConvergence :
+            Scierror(999, _("%s: Iterative procedure did not converge.\n"), "spchol");
+            return types::Function::Error;
+        case Eigen::InvalidInput :
+            Scierror(999, _("%s: The inputs are invalid, or the algorithm has been improperly called.\nWhen assertions are enabled, such errors trigger an assert.\n"), "spchol");
+            return types::Function::Error;
+        default :
+            break;
+    }
+
+    out.push_back(pSpFact);
+    out.push_back(pSpPermut);
+
+    return types::Function::OK;
+}
index 28ca824..5519900 100644 (file)
@@ -30,5 +30,6 @@ int SparseModule::Load()
     symbol::Context::getInstance()->addFunction(types::Function::createFunction(L"lusolve", &sci_lusolve, MODULE_NAME));
     symbol::Context::getInstance()->addFunction(types::Function::createFunction(L"luget", &sci_luget, MODULE_NAME));
     symbol::Context::getInstance()->addFunction(types::Function::createFunction(L"ludel", &sci_ludel, MODULE_NAME));
+    symbol::Context::getInstance()->addFunction(types::Function::createFunction(L"spchol", &sci_spchol, MODULE_NAME));
     return 1;
 }
index 58ccd3c..a9c8617 100644 (file)
     <ClCompile Include="sci_ordmmd.cpp" />
     <ClCompile Include="sci_sp2adj.cpp" />
     <ClCompile Include="sci_sparse.cpp" />
+    <ClCompile Include="sci_spchol.cpp" />
     <ClCompile Include="sci_spcompack.cpp" />
     <ClCompile Include="sci_spget.cpp" />
     <ClCompile Include="sci_spones.cpp" />
index e7e47bf..98d52e0 100644 (file)
@@ -63,6 +63,9 @@
     <ClCompile Include="sci_ludel.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="sci_spchol.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\includes\sparse_gw.hxx">
index 8b8a94e..3458c14 100644 (file)
@@ -23,8 +23,8 @@ X=[3.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ; ...
    0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ; ...
    2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ; ...
    0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ; ...
-   0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.]; 
+   0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.];
 X(1,1) = X(1,1) + %i;
 X      = sparse(X);
-msgerr = msprintf(gettext("Wrong type for argument #%d: Real matrix expected.\n"), 1);
-assert_checkerror('[R,P] = spchol(X);', msgerr);
+msgerr = msprintf(gettext("%s: Wrong type for argument #%d: Real matrix expected.\n"), "spchol", 1);
+assert_checkerror("[R,P] = spchol(X);", msgerr);
index 815765d..87faa56 100644 (file)
 // the function spchol was returning a wrong error message if X is complex
 
 X=[3.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ; ...
-   0.,  5.,  4.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ; ...
-   0.,  4.,  5.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ; ...
-   2.,  0.,  0.,  3.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ; ...
-   0.,  0.,  0.,  0. , 5.,  0.,  0.,  0.,  0.,  0.,  4. ; ...
-   0.,  0.,  0.,  0.,  0.,  4.,  0.,  3.,  0.,  3.,  0. ; ...
-   2.,  0.,  0.,  2.,  0.,  0.,  3.,  0.,  2.,  0.,  0. ; ...
-   0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ; ...
-   2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ; ...
-   0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ; ...
-   0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.]; 
+0.,  5.,  4.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ; ...
+0.,  4.,  5.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ; ...
+2.,  0.,  0.,  3.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ; ...
+0.,  0.,  0.,  0. , 5.,  0.,  0.,  0.,  0.,  0.,  4. ; ...
+0.,  0.,  0.,  0.,  0.,  4.,  0.,  3.,  0.,  3.,  0. ; ...
+2.,  0.,  0.,  2.,  0.,  0.,  3.,  0.,  2.,  0.,  0. ; ...
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ; ...
+2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ; ...
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ; ...
+0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.];
 
 X(1,1) = X(1,1) + %i;
 X      = sparse(X);
 
-msgerr = msprintf(gettext("Wrong type for argument #%d: Real matrix expected.\n"), 1);
-assert_checkerror('[R,P] = spchol(X);', msgerr);
+msgerr = msprintf(gettext("%s: Wrong type for argument #%d: Real matrix expected.\n"), "spchol", 1);
+assert_checkerror("[R,P] = spchol(X);", msgerr);