[Rcpp-commits] r3083 - in pkg/RcppEigen/inst/include: . unsupported unsupported/Eigen unsupported/Eigen/src unsupported/Eigen/src/SparseExtra
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 15 21:06:44 CEST 2011
Author: dmbates
Date: 2011-06-15 21:06:44 +0200 (Wed, 15 Jun 2011)
New Revision: 3083
Added:
pkg/RcppEigen/inst/include/unsupported/
pkg/RcppEigen/inst/include/unsupported/Eigen/
pkg/RcppEigen/inst/include/unsupported/Eigen/SparseExtra
pkg/RcppEigen/inst/include/unsupported/Eigen/src/
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/Amd.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CMakeLists.txt
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupportLegacy.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/RandomSetter.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/Solve.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/SparseLDLTLegacy.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/SparseLLT.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/SparseLU.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h
pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h
Log:
Add support for sparse Cholesky and LU
Added: pkg/RcppEigen/inst/include/unsupported/Eigen/SparseExtra
===================================================================
--- pkg/RcppEigen/inst/include/unsupported/Eigen/SparseExtra (rev 0)
+++ pkg/RcppEigen/inst/include/unsupported/Eigen/SparseExtra 2011-06-15 19:06:44 UTC (rev 3083)
@@ -0,0 +1,69 @@
+#ifndef EIGEN_SPARSE_EXTRA_MODULE_H
+#define EIGEN_SPARSE_EXTRA_MODULE_H
+
+#include "../../Eigen/Sparse"
+
+#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
+
+#include <vector>
+#include <map>
+#include <cstdlib>
+#include <cstring>
+#include <algorithm>
+
+#ifdef EIGEN_GOOGLEHASH_SUPPORT
+ #include <google/dense_hash_map>
+#endif
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+ * \defgroup SparseExtra_Module SparseExtra module
+ *
+ * This module contains some experimental features extending the sparse module.
+ *
+ * \code
+ * #include <Eigen/SparseExtra>
+ * \endcode
+ */
+
+struct DefaultBackend {};
+
+
+// solver flags
+enum {
+ CompleteFactorization = 0x0000, // the default
+ IncompleteFactorization = 0x0001,
+ MemoryEfficient = 0x0002,
+
+ // For LLT Cholesky:
+ SupernodalMultifrontal = 0x0010,
+ SupernodalLeftLooking = 0x0020,
+
+ // Ordering methods:
+ NaturalOrdering = 0x0100, // the default
+ MinimumDegree_AT_PLUS_A = 0x0200,
+ MinimumDegree_ATA = 0x0300,
+ ColApproxMinimumDegree = 0x0400,
+ Metis = 0x0500,
+ Scotch = 0x0600,
+ Chaco = 0x0700,
+ OrderingMask = 0x0f00
+};
+
+#include "../../Eigen/src/misc/Solve.h"
+
+#include "src/SparseExtra/RandomSetter.h"
+#include "src/SparseExtra/Solve.h"
+#include "src/SparseExtra/Amd.h"
+#include "src/SparseExtra/SimplicialCholesky.h"
+
+#include "src/SparseExtra/SparseLLT.h"
+#include "src/SparseExtra/SparseLDLTLegacy.h"
+#include "src/SparseExtra/SparseLU.h"
+
+} // namespace Eigen
+
+#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
+
+#endif // EIGEN_SPARSE_EXTRA_MODULE_H
Added: pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/Amd.h
===================================================================
--- pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/Amd.h (rev 0)
+++ pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/Amd.h 2011-06-15 19:06:44 UTC (rev 3083)
@@ -0,0 +1,448 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud at inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/*
+
+NOTE: this routine has been adapted from the CSparse library:
+
+Copyright (c) 2006, Timothy A. Davis.
+http://www.cise.ufl.edu/research/sparse/CSparse
+
+CSparse is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+CSparse is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+*/
+
+#ifndef EIGEN_SPARSE_AMD_H
+#define EIGEN_SPARSE_AMD_H
+
+namespace internal {
+
+
+#define CS_FLIP(i) (-(i)-2)
+#define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
+#define CS_MARKED(w,j) (w[j] < 0)
+#define CS_MARK(w,j) { w[j] = CS_FLIP (w[j]); }
+
+/* clear w */
+template<typename Index>
+static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
+{
+ Index k;
+ if(mark < 2 || (mark + lemax < 0))
+ {
+ for(k = 0; k < n; k++)
+ if(w[k] != 0)
+ w[k] = 1;
+ mark = 2;
+ }
+ return (mark); /* at this point, w[0..n-1] < mark holds */
+}
+
+/* depth-first search and postorder of a tree rooted at node j */
+template<typename Index>
+Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
+{
+ int i, p, top = 0;
+ if(!head || !next || !post || !stack) return (-1); /* check inputs */
+ stack[0] = j; /* place j on the stack */
+ while (top >= 0) /* while (stack is not empty) */
+ {
+ p = stack[top]; /* p = top of stack */
+ i = head[p]; /* i = youngest child of p */
+ if(i == -1)
+ {
+ top--; /* p has no unordered children left */
+ post[k++] = p; /* node p is the kth postordered node */
+ }
+ else
+ {
+ head[p] = next[i]; /* remove i from children of p */
+ stack[++top] = i; /* start dfs on child node i */
+ }
+ }
+ return k;
+}
+
+
+/** \internal
+ * Approximate minimum degree ordering algorithm.
+ * \returns the permutation P reducing the fill-in of the input matrix \a C
+ * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.
+ * On exit the values of C are destroyed */
+template<typename Scalar, typename Index>
+void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic>& perm)
+{
+ typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
+
+ int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
+ k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
+ ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
+ unsigned int h;
+
+ Index n = C.cols();
+ dense = std::max<Index> (16, 10 * sqrt ((double) n)); /* find dense threshold */
+ dense = std::min<Index> (n-2, dense);
+
+ Index cnz = C.nonZeros();
+ perm.resize(n+1);
+ t = cnz + cnz/5 + 2*n; /* add elbow room to C */
+ C.resizeNonZeros(t);
+
+ Index* W = new Index[8*(n+1)]; /* get workspace */
+ Index* len = W;
+ Index* nv = W + (n+1);
+ Index* next = W + 2*(n+1);
+ Index* head = W + 3*(n+1);
+ Index* elen = W + 4*(n+1);
+ Index* degree = W + 5*(n+1);
+ Index* w = W + 6*(n+1);
+ Index* hhead = W + 7*(n+1);
+ Index* last = perm.indices().data(); /* use P as workspace for last */
+
+ /* --- Initialize quotient graph ---------------------------------------- */
+ Index* Cp = C._outerIndexPtr();
+ Index* Ci = C._innerIndexPtr();
+ for(k = 0; k < n; k++)
+ len[k] = Cp[k+1] - Cp[k];
+ len[n] = 0;
+ nzmax = t;
+
+ for(i = 0; i <= n; i++)
+ {
+ head[i] = -1; // degree list i is empty
+ last[i] = -1;
+ next[i] = -1;
+ hhead[i] = -1; // hash list i is empty
+ nv[i] = 1; // node i is just one node
+ w[i] = 1; // node i is alive
+ elen[i] = 0; // Ek of node i is empty
+ degree[i] = len[i]; // degree of node i
+ }
+ mark = cs_wclear (0, 0, w, n); /* clear w */
+ elen[n] = -2; /* n is a dead element */
+ Cp[n] = -1; /* n is a root of assembly tree */
+ w[n] = 0; /* n is a dead element */
+
+ /* --- Initialize degree lists ------------------------------------------ */
+ for(i = 0; i < n; i++)
+ {
+ d = degree[i];
+ if(d == 0) /* node i is empty */
+ {
+ elen[i] = -2; /* element i is dead */
+ nel++;
+ Cp[i] = -1; /* i is a root of assembly tree */
+ w[i] = 0;
+ }
+ else if(d > dense) /* node i is dense */
+ {
+ nv[i] = 0; /* absorb i into element n */
+ elen[i] = -1; /* node i is dead */
+ nel++;
+ Cp[i] = CS_FLIP (n);
+ nv[n]++;
+ }
+ else
+ {
+ if(head[d] != -1) last[head[d]] = i;
+ next[i] = head[d]; /* put node i in degree list d */
+ head[d] = i;
+ }
+ }
+
+ while (nel < n) /* while (selecting pivots) do */
+ {
+ /* --- Select node of minimum approximate degree -------------------- */
+ for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
+ if(next[k] != -1) last[next[k]] = -1;
+ head[mindeg] = next[k]; /* remove k from degree list */
+ elenk = elen[k]; /* elenk = |Ek| */
+ nvk = nv[k]; /* # of nodes k represents */
+ nel += nvk; /* nv[k] nodes of A eliminated */
+
+ /* --- Garbage collection ------------------------------------------- */
+ if(elenk > 0 && cnz + mindeg >= nzmax)
+ {
+ for(j = 0; j < n; j++)
+ {
+ if((p = Cp[j]) >= 0) /* j is a live node or element */
+ {
+ Cp[j] = Ci[p]; /* save first entry of object */
+ Ci[p] = CS_FLIP (j); /* first entry is now CS_FLIP(j) */
+ }
+ }
+ for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
+ {
+ if((j = CS_FLIP (Ci[p++])) >= 0) /* found object j */
+ {
+ Ci[q] = Cp[j]; /* restore first entry of object */
+ Cp[j] = q++; /* new pointer to object j */
+ for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
+ }
+ }
+ cnz = q; /* Ci[cnz...nzmax-1] now free */
+ }
+
+ /* --- Construct new element ---------------------------------------- */
+ dk = 0;
+ nv[k] = -nvk; /* flag k as in Lk */
+ p = Cp[k];
+ pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */
+ pk2 = pk1;
+ for(k1 = 1; k1 <= elenk + 1; k1++)
+ {
+ if(k1 > elenk)
+ {
+ e = k; /* search the nodes in k */
+ pj = p; /* list of nodes starts at Ci[pj]*/
+ ln = len[k] - elenk; /* length of list of nodes in k */
+ }
+ else
+ {
+ e = Ci[p++]; /* search the nodes in e */
+ pj = Cp[e];
+ ln = len[e]; /* length of list of nodes in e */
+ }
+ for(k2 = 1; k2 <= ln; k2++)
+ {
+ i = Ci[pj++];
+ if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
+ dk += nvi; /* degree[Lk] += size of node i */
+ nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/
+ Ci[pk2++] = i; /* place i in Lk */
+ if(next[i] != -1) last[next[i]] = last[i];
+ if(last[i] != -1) /* remove i from degree list */
+ {
+ next[last[i]] = next[i];
+ }
+ else
+ {
+ head[degree[i]] = next[i];
+ }
+ }
+ if(e != k)
+ {
+ Cp[e] = CS_FLIP (k); /* absorb e into k */
+ w[e] = 0; /* e is now a dead element */
+ }
+ }
+ if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */
+ degree[k] = dk; /* external degree of k - |Lk\i| */
+ Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */
+ len[k] = pk2 - pk1;
+ elen[k] = -2; /* k is now an element */
+
+ /* --- Find set differences ----------------------------------------- */
+ mark = cs_wclear (mark, lemax, w, n); /* clear w if necessary */
+ for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */
+ {
+ i = Ci[pk];
+ if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
+ nvi = -nv[i]; /* nv[i] was negated */
+ wnvi = mark - nvi;
+ for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */
+ {
+ e = Ci[p];
+ if(w[e] >= mark)
+ {
+ w[e] -= nvi; /* decrement |Le\Lk| */
+ }
+ else if(w[e] != 0) /* ensure e is a live element */
+ {
+ w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
+ }
+ }
+ }
+
+ /* --- Degree update ------------------------------------------------ */
+ for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */
+ {
+ i = Ci[pk]; /* consider node i in Lk */
+ p1 = Cp[i];
+ p2 = p1 + elen[i] - 1;
+ pn = p1;
+ for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */
+ {
+ e = Ci[p];
+ if(w[e] != 0) /* e is an unabsorbed element */
+ {
+ dext = w[e] - mark; /* dext = |Le\Lk| */
+ if(dext > 0)
+ {
+ d += dext; /* sum up the set differences */
+ Ci[pn++] = e; /* keep e in Ei */
+ h += e; /* compute the hash of node i */
+ }
+ else
+ {
+ Cp[e] = CS_FLIP (k); /* aggressive absorb. e->k */
+ w[e] = 0; /* e is a dead element */
+ }
+ }
+ }
+ elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */
+ p3 = pn;
+ p4 = p1 + len[i];
+ for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
+ {
+ j = Ci[p];
+ if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
+ d += nvj; /* degree(i) += |j| */
+ Ci[pn++] = j; /* place j in node list of i */
+ h += j; /* compute hash for node i */
+ }
+ if(d == 0) /* check for mass elimination */
+ {
+ Cp[i] = CS_FLIP (k); /* absorb i into k */
+ nvi = -nv[i];
+ dk -= nvi; /* |Lk| -= |i| */
+ nvk += nvi; /* |k| += nv[i] */
+ nel += nvi;
+ nv[i] = 0;
+ elen[i] = -1; /* node i is dead */
+ }
+ else
+ {
+ degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */
+ Ci[pn] = Ci[p3]; /* move first node to end */
+ Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */
+ Ci[p1] = k; /* add k as 1st element in of Ei */
+ len[i] = pn - p1 + 1; /* new len of adj. list of node i */
+ h %= n; /* finalize hash of i */
+ next[i] = hhead[h]; /* place i in hash bucket */
+ hhead[h] = i;
+ last[i] = h; /* save hash of i in last[i] */
+ }
+ } /* scan2 is done */
+ degree[k] = dk; /* finalize |Lk| */
+ lemax = std::max<Index>(lemax, dk);
+ mark = cs_wclear (mark+lemax, lemax, w, n); /* clear w */
+
+ /* --- Supernode detection ------------------------------------------ */
+ for(pk = pk1; pk < pk2; pk++)
+ {
+ i = Ci[pk];
+ if(nv[i] >= 0) continue; /* skip if i is dead */
+ h = last[i]; /* scan hash bucket of node i */
+ i = hhead[h];
+ hhead[h] = -1; /* hash bucket will be empty */
+ for(; i != -1 && next[i] != -1; i = next[i], mark++)
+ {
+ ln = len[i];
+ eln = elen[i];
+ for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
+ jlast = i;
+ for(j = next[i]; j != -1; ) /* compare i with all j */
+ {
+ ok = (len[j] == ln) && (elen[j] == eln);
+ for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
+ {
+ if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/
+ }
+ if(ok) /* i and j are identical */
+ {
+ Cp[j] = CS_FLIP (i); /* absorb j into i */
+ nv[i] += nv[j];
+ nv[j] = 0;
+ elen[j] = -1; /* node j is dead */
+ j = next[j]; /* delete j from hash bucket */
+ next[jlast] = j;
+ }
+ else
+ {
+ jlast = j; /* j and i are different */
+ j = next[j];
+ }
+ }
+ }
+ }
+
+ /* --- Finalize new element------------------------------------------ */
+ for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */
+ {
+ i = Ci[pk];
+ if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
+ nv[i] = nvi; /* restore nv[i] */
+ d = degree[i] + dk - nvi; /* compute external degree(i) */
+ d = std::min<Index> (d, n - nel - nvi);
+ if(head[d] != -1) last[head[d]] = i;
+ next[i] = head[d]; /* put i back in degree list */
+ last[i] = -1;
+ head[d] = i;
+ mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */
+ degree[i] = d;
+ Ci[p++] = i; /* place i in Lk */
+ }
+ nv[k] = nvk; /* # nodes absorbed into k */
+ if((len[k] = p-pk1) == 0) /* length of adj list of element k*/
+ {
+ Cp[k] = -1; /* k is a root of the tree */
+ w[k] = 0; /* k is now a dead element */
+ }
+ if(elenk != 0) cnz = p; /* free unused space in Lk */
+ }
+
+ /* --- Postordering ----------------------------------------------------- */
+ for(i = 0; i < n; i++) Cp[i] = CS_FLIP (Cp[i]);/* fix assembly tree */
+ for(j = 0; j <= n; j++) head[j] = -1;
+ for(j = n; j >= 0; j--) /* place unordered nodes in lists */
+ {
+ if(nv[j] > 0) continue; /* skip if j is an element */
+ next[j] = head[Cp[j]]; /* place j in list of its parent */
+ head[Cp[j]] = j;
+ }
+ for(e = n; e >= 0; e--) /* place elements in lists */
+ {
+ if(nv[e] <= 0) continue; /* skip unless e is an element */
+ if(Cp[e] != -1)
+ {
+ next[e] = head[Cp[e]]; /* place e in list of its parent */
+ head[Cp[e]] = e;
+ }
+ }
+ for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */
+ {
+ if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, perm.indices().data(), w);
+ }
+
+ perm.indices().conservativeResize(n);
+
+ delete[] W;
+}
+
+} // namespace internal
+
+#endif // EIGEN_SPARSE_AMD_H
Added: pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CMakeLists.txt
===================================================================
--- pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CMakeLists.txt (rev 0)
+++ pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CMakeLists.txt 2011-06-15 19:06:44 UTC (rev 3083)
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_SparseExtra_SRCS "*.h")
+
+INSTALL(FILES
+ ${Eigen_SparseExtra_SRCS}
+ DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/SparseExtra COMPONENT Devel
+ )
Added: pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
===================================================================
--- pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h (rev 0)
+++ pkg/RcppEigen/inst/include/unsupported/Eigen/src/SparseExtra/CholmodSupport.h 2011-06-15 19:06:44 UTC (rev 3083)
@@ -0,0 +1,399 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud at inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_CHOLMODSUPPORT_H
+#define EIGEN_CHOLMODSUPPORT_H
+
+namespace internal {
+
+template<typename Scalar, typename CholmodType>
+void cholmod_configure_matrix(CholmodType& mat)
+{
+ if (internal::is_same<Scalar,float>::value)
+ {
+ mat.xtype = CHOLMOD_REAL;
+ mat.dtype = CHOLMOD_SINGLE;
+ }
+ else if (internal::is_same<Scalar,double>::value)
+ {
+ mat.xtype = CHOLMOD_REAL;
+ mat.dtype = CHOLMOD_DOUBLE;
+ }
+ else if (internal::is_same<Scalar,std::complex<float> >::value)
+ {
+ mat.xtype = CHOLMOD_COMPLEX;
+ mat.dtype = CHOLMOD_SINGLE;
+ }
+ else if (internal::is_same<Scalar,std::complex<double> >::value)
+ {
+ mat.xtype = CHOLMOD_COMPLEX;
+ mat.dtype = CHOLMOD_DOUBLE;
+ }
+ else
+ {
+ eigen_assert(false && "Scalar type not supported by CHOLMOD");
+ }
+}
+
+} // namespace internal
+
+/** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
+ * Note that the data are shared.
+ */
+template<typename _Scalar, int _Options, typename _Index>
+cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
+{
+ typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
+ cholmod_sparse res;
+ res.nzmax = mat.nonZeros();
+ res.nrow = mat.rows();;
+ res.ncol = mat.cols();
+ res.p = mat._outerIndexPtr();
+ res.i = mat._innerIndexPtr();
+ res.x = mat._valuePtr();
+ res.sorted = 1;
+ res.packed = 1;
+ res.dtype = 0;
+ res.stype = -1;
+
+ if (internal::is_same<_Index,int>::value)
+ {
+ res.itype = CHOLMOD_INT;
+ }
+ else
+ {
+ eigen_assert(false && "Index type different than int is not supported yet");
+ }
+
+ // setup res.xtype
+ internal::cholmod_configure_matrix<_Scalar>(res);
+
+ res.stype = 0;
+
+ return res;
+}
+
+template<typename _Scalar, int _Options, typename _Index>
+const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
+{
+ cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
+ return res;
+}
+
+/** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
+ * The data are not copied but shared. */
+template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
+cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
+{
+ cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
+
+ if(UpLo==Upper) res.stype = 1;
+ if(UpLo==Lower) res.stype = -1;
+
+ return res;
+}
+
+/** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
+ * The data are not copied but shared. */
+template<typename Derived>
+cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
+{
+ EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
+ typedef typename Derived::Scalar Scalar;
+
+ cholmod_dense res;
+ res.nrow = mat.rows();
+ res.ncol = mat.cols();
+ res.nzmax = res.nrow * res.ncol;
+ res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
+ res.x = mat.derived().data();
+ res.z = 0;
+
+ internal::cholmod_configure_matrix<Scalar>(res);
+
+ return res;
+}
+
+/** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
+ * The data are not copied but shared. */
+template<typename Scalar, int Flags, typename Index>
+MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
+{
+ return MappedSparseMatrix<Scalar,Flags,Index>
+ (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
+ reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
+}
+
+enum CholmodMode {
+ CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
+};
+
+/** \brief A Cholesky factorization and solver based on Cholmod
+ *
+ * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
+ * using the Cholmod library. The sparse matrix A must be selfajoint and positive definite. The vectors or matrices
+ * X and B can be either dense or sparse.
+ *
+ * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
+ * or Upper. Default is Lower.
+ *
+ */
+template<typename _MatrixType, int _UpLo = Lower>
+class CholmodDecomposition
+{
+ public:
+ typedef _MatrixType MatrixType;
+ enum { UpLo = _UpLo };
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef MatrixType CholMatrixType;
+ typedef typename MatrixType::Index Index;
+
+ public:
+
+ CholmodDecomposition()
+ : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
+ {
+ cholmod_start(&m_cholmod);
+ setMode(CholmodLDLt);
+ }
+
+ CholmodDecomposition(const MatrixType& matrix)
+ : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
+ {
+ cholmod_start(&m_cholmod);
+ compute(matrix);
+ }
+
+ ~CholmodDecomposition()
+ {
+ if(m_cholmodFactor)
+ cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+ cholmod_finish(&m_cholmod);
+ }
+
+ inline Index cols() const { return m_cholmodFactor->n; }
+ inline Index rows() const { return m_cholmodFactor->n; }
+
+ void setMode(CholmodMode mode)
+ {
+ switch(mode)
+ {
+ case CholmodAuto:
+ m_cholmod.final_asis = 1;
+ m_cholmod.supernodal = CHOLMOD_AUTO;
+ break;
+ case CholmodSimplicialLLt:
+ m_cholmod.final_asis = 0;
+ m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+ m_cholmod.final_ll = 1;
+ break;
+ case CholmodSupernodalLLt:
+ m_cholmod.final_asis = 1;
+ m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
+ break;
+ case CholmodLDLt:
+ m_cholmod.final_asis = 1;
+ m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+ break;
+ default:
+ break;
+ }
+ }
+
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful,
+ * \c NumericalIssue if the matrix.appears to be negative.
+ */
+ ComputationInfo info() const
+ {
+ eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+ return m_info;
+ }
+
+ /** Computes the sparse Cholesky decomposition of \a matrix */
+ void compute(const MatrixType& matrix)
+ {
+ analyzePattern(matrix);
+ factorize(matrix);
+ }
+
+ /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+ *
+ * \sa compute()
+ */
+ template<typename Rhs>
+ inline const internal::solve_retval<CholmodDecomposition, Rhs>
+ solve(const MatrixBase<Rhs>& b) const
+ {
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(rows()==b.rows()
+ && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
+ return internal::solve_retval<CholmodDecomposition, Rhs>(*this, b.derived());
+ }
+
+ /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+ *
+ * \sa compute()
+ */
+ template<typename Rhs>
+ inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs>
+ solve(const SparseMatrixBase<Rhs>& b) const
+ {
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(rows()==b.rows()
+ && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
+ return internal::sparse_solve_retval<CholmodDecomposition, Rhs>(*this, b.derived());
+ }
+
+ /** Performs a symbolic decomposition on the sparcity of \a matrix.
+ *
+ * This function is particularly useful when solving for several problems having the same structure.
+ *
+ * \sa factorize()
+ */
+ void analyzePattern(const MatrixType& matrix)
+ {
+ if(m_cholmodFactor)
+ {
+ cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+ m_cholmodFactor = 0;
+ }
+ cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
+ m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
+
+ this->m_isInitialized = true;
+ this->m_info = Success;
+ m_analysisIsOk = true;
+ m_factorizationIsOk = false;
+ }
+
+ /** Performs a numeric decomposition of \a matrix
+ *
+ * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
+ *
+ * \sa analyzePattern()
+ */
+ void factorize(const MatrixType& matrix)
+ {
+ eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+ cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
+ cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
+
+ this->m_info = Success;
+ m_factorizationIsOk = true;
+ }
+
+ /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
+ * See the Cholmod user guide for details. */
+ cholmod_common& cholmod() { return m_cholmod; }
+
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ /** \internal */
+ template<typename Rhs,typename Dest>
+ void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
+ {
+ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
+ const Index size = m_cholmodFactor->n;
+ eigen_assert(size==b.rows());
+
+ // note: cd stands for Cholmod Dense
+ cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
+ cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
+ if(!x_cd)
+ {
+ this->m_info = NumericalIssue;
+ }
+ // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
+ dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
+ cholmod_free_dense(&x_cd, &m_cholmod);
+ }
+
+ /** \internal */
+ template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
+ void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
+ {
+ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
+ const Index size = m_cholmodFactor->n;
+ eigen_assert(size==b.rows());
+
+ // note: cs stands for Cholmod Sparse
+ cholmod_sparse b_cs = viewAsCholmod(b);
+ cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
+ if(!x_cs)
+ {
+ this->m_info = NumericalIssue;
+ }
+ // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
+ dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
+ cholmod_free_sparse(&x_cs, &m_cholmod);
+ }
+ #endif // EIGEN_PARSED_BY_DOXYGEN
+
+ template<typename Stream>
+ void dumpMemory(Stream& s)
+ {}
+
+ protected:
+ mutable cholmod_common m_cholmod;
+ cholmod_factor* m_cholmodFactor;
+ mutable ComputationInfo m_info;
+ bool m_isInitialized;
+ int m_factorizationIsOk;
+ int m_analysisIsOk;
+};
+
+namespace internal {
+
+template<typename _MatrixType, int _UpLo, typename Rhs>
+struct solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
+ : solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
+{
+ typedef CholmodDecomposition<_MatrixType,_UpLo> Dec;
+ EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
+
+ template<typename Dest> void evalTo(Dest& dst) const
+ {
+ dec()._solve(rhs(),dst);
+ }
+};
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 3083
More information about the Rcpp-commits
mailing list