diff --git a/DESCRIPTION b/DESCRIPTION index 05138634..85f123ae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,7 +24,7 @@ License: GPL (>= 2) | file LICENSE LazyLoad: yes Depends: R (>= 3.6.0) LinkingTo: Rcpp -Imports: Matrix (>= 1.1-0), Rcpp (>= 0.11.0), stats, utils -Suggests: inline, tinytest, pkgKitten, microbenchmark +Imports: Rcpp (>= 0.11.0), stats, utils +Suggests: Matrix, inline, microbenchmark, pkgKitten, tinytest URL: https://github.com/RcppCore/RcppEigen, https://dirk.eddelbuettel.com/code/rcpp.eigen.html BugReports: https://github.com/RcppCore/RcppEigen/issues diff --git a/NAMESPACE b/NAMESPACE index fe46f9b4..740f024a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,5 @@ useDynLib("RcppEigen", .registration=TRUE) -importClassesFrom("Matrix", "dgCMatrix", "dgeMatrix", "dsCMatrix", "dtCMatrix") importFrom("Rcpp", "evalCpp") importFrom("utils", "packageDescription", "package.skeleton") importFrom("stats", "model.frame", "model.matrix", "model.response", "fitted", "coef", "printCoefmat", "pt") diff --git a/inst/include/Eigen/CholmodSupport b/inst/include/Eigen/CholmodSupport index 573b8221..fea39040 100644 --- a/inst/include/Eigen/CholmodSupport +++ b/inst/include/Eigen/CholmodSupport @@ -12,10 +12,6 @@ #include "src/Core/util/DisableStupidWarnings.h" -extern "C" { - #include -} - /** \ingroup Support_modules * \defgroup CholmodSupport_Module CholmodSupport module * diff --git a/inst/include/Eigen/src/CholmodSupport/CholmodSupport.h b/inst/include/Eigen/src/CholmodSupport/CholmodSupport.h index 57197202..5a7056a5 100644 --- a/inst/include/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/inst/include/Eigen/src/CholmodSupport/CholmodSupport.h @@ -10,6 +10,10 @@ #ifndef EIGEN_CHOLMODSUPPORT_H #define EIGEN_CHOLMODSUPPORT_H +#ifndef R_MATRIX_CHOLMOD +# define R_MATRIX_CHOLMOD(_NAME_) cholmod_ ## _NAME_ +#endif + namespace Eigen { namespace internal { @@ -195,7 +199,7 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + R_MATRIX_CHOLMOD(start)(&m_cholmod); } explicit CholmodBase(const MatrixType& matrix) @@ -203,15 +207,15 @@ class CholmodBase : public SparseSolverBase { EIGEN_STATIC_ASSERT((internal::is_same::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); m_shiftOffset[0] = m_shiftOffset[1] = 0.0; - cholmod_start(&m_cholmod); + R_MATRIX_CHOLMOD(start)(&m_cholmod); compute(matrix); } ~CholmodBase() { if(m_cholmodFactor) - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); - cholmod_finish(&m_cholmod); + R_MATRIX_CHOLMOD(free_factor)(&m_cholmodFactor, &m_cholmod); + R_MATRIX_CHOLMOD(finish)(&m_cholmod); } inline StorageIndex cols() const { return internal::convert_index(m_cholmodFactor->n); } @@ -246,11 +250,11 @@ class CholmodBase : public SparseSolverBase { if(m_cholmodFactor) { - cholmod_free_factor(&m_cholmodFactor, &m_cholmod); + R_MATRIX_CHOLMOD(free_factor)(&m_cholmodFactor, &m_cholmod); m_cholmodFactor = 0; } cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - m_cholmodFactor = cholmod_analyze(&A, &m_cholmod); + m_cholmodFactor = R_MATRIX_CHOLMOD(analyze)(&A, &m_cholmod); this->m_isInitialized = true; this->m_info = Success; @@ -268,7 +272,7 @@ class CholmodBase : public SparseSolverBase { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView()); - cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod); + R_MATRIX_CHOLMOD(factorize_p)(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod); // If the factorization failed, minor is the column at which it did. On success minor == n. this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue); @@ -293,7 +297,7 @@ class CholmodBase : public SparseSolverBase Ref > b_ref(b.derived()); cholmod_dense b_cd = viewAsCholmod(b_ref); - cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod); + cholmod_dense* x_cd = R_MATRIX_CHOLMOD(solve)(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod); if(!x_cd) { this->m_info = NumericalIssue; @@ -301,7 +305,7 @@ class CholmodBase : public SparseSolverBase } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) dest = Matrix::Map(reinterpret_cast(x_cd->x),b.rows(),b.cols()); - cholmod_free_dense(&x_cd, &m_cholmod); + R_MATRIX_CHOLMOD(free_dense)(&x_cd, &m_cholmod); } /** \internal */ @@ -316,7 +320,7 @@ class CholmodBase : public SparseSolverBase // note: cs stands for Cholmod Sparse Ref > b_ref(b.const_cast_derived()); cholmod_sparse b_cs = viewAsCholmod(b_ref); - cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); + cholmod_sparse* x_cs = R_MATRIX_CHOLMOD(spsolve)(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); if(!x_cs) { this->m_info = NumericalIssue; @@ -324,7 +328,7 @@ class CholmodBase : public SparseSolverBase } // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) dest.derived() = viewAsEigen(*x_cs); - cholmod_free_sparse(&x_cs, &m_cholmod); + R_MATRIX_CHOLMOD(free_sparse)(&x_cs, &m_cholmod); } #endif // EIGEN_PARSED_BY_DOXYGEN diff --git a/inst/include/RcppEigenCholmod.h b/inst/include/RcppEigenCholmod.h index 5f169aa3..147a2553 100644 --- a/inst/include/RcppEigenCholmod.h +++ b/inst/include/RcppEigenCholmod.h @@ -1,6 +1,8 @@ // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*- // -// cholmod.h: selected headers from Tim Davis's CHOLMOD package +// RcppEigenCholmod.h: Provide access to the Matrix API and in turn +// to Eigen's CholmodSupport module. Use of this header relies on +// compilation of ../../src/RcppEigenStubs.cpp and LinkingTo: Matrix. // // Copyright (C) 2011 Douglas Bates, Martin Maechler, Dirk Eddelbuettel and Romain Francois // @@ -19,1072 +21,15 @@ // You should have received a copy of the GNU General Public License // along with RcppEigen. If not, see . -#ifndef RcppEigen_CHOLMOD_H -#define RcppEigen_CHOLMOD_H +#ifndef RcppEigen__RcppEigenCholmod__h +#define RcppEigen__RcppEigenCholmod__h -/* Original of this is Matrix/inst/include/cholmod.h --- sync with it (e.g for lme4) ! - ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - */ -#include -#include - -// Rather use C99 -- which we require in R anyway -#include - - -#ifdef __cplusplus -extern "C" { +#include +#ifndef R_MATRIX_CHOLMOD /* Matrix <= 1.6-1.1 */ +# define R_MATRIX_CHOLMOD(_NAME_) M_cholmod_ ## _NAME_ +# define M_cholmod_start M_R_cholmod_start /* sigh */ #endif -// from Matrix/src/SuiteSparse_config/SuiteSparse_config.h - line 51 : -#ifndef SuiteSparse_long - -/* #ifdef _WIN64 */ - -/* #define SuiteSparse_long __int64 */ -/* #define SuiteSparse_long_max _I64_MAX */ -/* #define SuiteSparse_long_idd "I64d" */ - -/* #else */ - -/* #define SuiteSparse_long long */ -/* #define SuiteSparse_long_max LONG_MAX */ -/* #define SuiteSparse_long_idd "ld" */ - -/* #endif */ - -#define SuiteSparse_long int64_t - // typically long (but on WIN64) -#define SuiteSparse_ulong uint64_t - // only needed for ../COLAMD/Source/colamd.c (original has 'unsigned Int' which fails!!) -#define SuiteSparse_long_max 9223372036854775801 - // typically LONG_MAX (but ..) -#define SuiteSparse_long_idd PRId64 - // typically "ld" - -#define SuiteSparse_long_id "%" SuiteSparse_long_idd -#endif - -/* For backward compatibility with prior versions of SuiteSparse. The UF_* - * macros are deprecated and will be removed in a future version. */ -#ifndef UF_long -#define UF_long SuiteSparse_long -#define UF_long_max SuiteSparse_long_max -#define UF_long_idd SuiteSparse_long_idd -#define UF_long_id SuiteSparse_long_id -#endif - - -#define CHOLMOD_HAS_VERSION_FUNCTION - -#define CHOLMOD_DATE "April 25, 2013" -#define CHOLMOD_VER_CODE(main,sub) ((main) * 1000 + (sub)) -#define CHOLMOD_MAIN_VERSION 2 -#define CHOLMOD_SUB_VERSION 1 -#define CHOLMOD_SUBSUB_VERSION 2 -#define CHOLMOD_VERSION \ - CHOLMOD_VER_CODE(CHOLMOD_MAIN_VERSION,CHOLMOD_SUB_VERSION) -// from ../../src/CHOLMOD/Include/cholmod_core.h - line 275 : -/* Each CHOLMOD object has its own type code. */ - -#define CHOLMOD_COMMON 0 -#define CHOLMOD_SPARSE 1 -#define CHOLMOD_FACTOR 2 -#define CHOLMOD_DENSE 3 -#define CHOLMOD_TRIPLET 4 - -/* ========================================================================== */ -/* === CHOLMOD Common ======================================================= */ -/* ========================================================================== */ - -/* itype defines the types of integer used: */ -#define CHOLMOD_INT 0 /* all integer arrays are int */ -#define CHOLMOD_INTLONG 1 /* most are int, some are SuiteSparse_long */ -#define CHOLMOD_LONG 2 /* all integer arrays are SuiteSparse_long */ - -/* The itype of all parameters for all CHOLMOD routines must match. - * FUTURE WORK: CHOLMOD_INTLONG is not yet supported. - */ - -/* dtype defines what the numerical type is (double or float): */ -#define CHOLMOD_DOUBLE 0 /* all numerical values are double */ -#define CHOLMOD_SINGLE 1 /* all numerical values are float */ - -/* The dtype of all parameters for all CHOLMOD routines must match. - * - * Scalar floating-point values are always passed as double arrays of size 2 - * (for the real and imaginary parts). They are typecast to float as needed. - * FUTURE WORK: the float case is not supported yet. - */ - -/* xtype defines the kind of numerical values used: */ -#define CHOLMOD_PATTERN 0 /* pattern only, no numerical values */ -#define CHOLMOD_REAL 1 /* a real matrix */ -#define CHOLMOD_COMPLEX 2 /* a complex matrix (ANSI C99 compatible) */ -#define CHOLMOD_ZOMPLEX 3 /* a complex matrix (MATLAB compatible) */ - -/* Definitions for cholmod_common: */ -#define CHOLMOD_MAXMETHODS 9 /* maximum number of different methods that */ - /* cholmod_analyze can try. Must be >= 9. */ - -/* Common->status values. zero means success, negative means a fatal error, - * positive is a warning. */ -#define CHOLMOD_OK 0 /* success */ -#define CHOLMOD_NOT_INSTALLED (-1) /* failure: method not installed */ -#define CHOLMOD_OUT_OF_MEMORY (-2) /* failure: out of memory */ -#define CHOLMOD_TOO_LARGE (-3) /* failure: integer overflow occured */ -#define CHOLMOD_INVALID (-4) /* failure: invalid input */ -#define CHOLMOD_GPU_PROBLEM (-5) /* failure: GPU fatal error */ -#define CHOLMOD_NOT_POSDEF (1) /* warning: matrix not pos. def. */ -#define CHOLMOD_DSMALL (2) /* warning: D for LDL' or diag(L) or */ - /* LL' has tiny absolute value */ - -/* ordering method (also used for L->ordering) */ -#define CHOLMOD_NATURAL 0 /* use natural ordering */ -#define CHOLMOD_GIVEN 1 /* use given permutation */ -#define CHOLMOD_AMD 2 /* use minimum degree (AMD) */ -#define CHOLMOD_METIS 3 /* use METIS' nested dissection */ -#define CHOLMOD_NESDIS 4 /* use CHOLMOD's version of nested dissection:*/ - /* node bisector applied recursively, followed - * by constrained minimum degree (CSYMAMD or - * CCOLAMD) */ -#define CHOLMOD_COLAMD 5 /* use AMD for A, COLAMD for A*A' */ - -/* POSTORDERED is not a method, but a result of natural ordering followed by a - * weighted postorder. It is used for L->ordering, not method [ ].ordering. */ -#define CHOLMOD_POSTORDERED 6 /* natural ordering, postordered. */ - -/* supernodal strategy (for Common->supernodal) */ -#define CHOLMOD_SIMPLICIAL 0 /* always do simplicial */ -#define CHOLMOD_AUTO 1 /* select simpl/super depending on matrix */ -#define CHOLMOD_SUPERNODAL 2 /* always do supernodal */ - -typedef struct cholmod_common_struct -{ - /* ---------------------------------------------------------------------- */ - /* parameters for symbolic/numeric factorization and update/downdate */ - /* ---------------------------------------------------------------------- */ - - double dbound ; /* Smallest absolute value of diagonal entries of D - * for LDL' factorization and update/downdate/rowadd/ - * rowdel, or the diagonal of L for an LL' factorization. - * Entries in the range 0 to dbound are replaced with dbound. - * Entries in the range -dbound to 0 are replaced with -dbound. No - * changes are made to the diagonal if dbound <= 0. Default: zero */ - - double grow0 ; /* For a simplicial factorization, L->i and L->x can - * grow if necessary. grow0 is the factor by which - * it grows. For the initial space, L is of size MAX (1,grow0) times - * the required space. If L runs out of space, the new size of L is - * MAX(1.2,grow0) times the new required space. If you do not plan on - * modifying the LDL' factorization in the Modify module, set grow0 to - * zero (or set grow2 to 0, see below). Default: 1.2 */ - - double grow1 ; - - size_t grow2 ; /* For a simplicial factorization, each column j of L - * is initialized with space equal to - * grow1*L->ColCount[j] + grow2. If grow0 < 1, grow1 < 1, or grow2 == 0, - * then the space allocated is exactly equal to L->ColCount[j]. If the - * column j runs out of space, it increases to grow1*need + grow2 in - * size, where need is the total # of nonzeros in that column. If you do - * not plan on modifying the factorization in the Modify module, set - * grow2 to zero. Default: grow1 = 1.2, grow2 = 5. */ - - size_t maxrank ; /* rank of maximum update/downdate. Valid values: - * 2, 4, or 8. A value < 2 is set to 2, and a - * value > 8 is set to 8. It is then rounded up to the next highest - * power of 2, if not already a power of 2. Workspace (Xwork, below) of - * size nrow-by-maxrank double's is allocated for the update/downdate. - * If an update/downdate of rank-k is requested, with k > maxrank, - * it is done in steps of maxrank. Default: 8, which is fastest. - * Memory usage can be reduced by setting maxrank to 2 or 4. - */ - - double supernodal_switch ; /* supernodal vs simplicial factorization */ - int supernodal ; /* If Common->supernodal <= CHOLMOD_SIMPLICIAL - * (0) then cholmod_analyze performs a - * simplicial analysis. If >= CHOLMOD_SUPERNODAL (2), then a supernodal - * analysis is performed. If == CHOLMOD_AUTO (1) and - * flop/nnz(L) < Common->supernodal_switch, then a simplicial analysis - * is done. A supernodal analysis done otherwise. - * Default: CHOLMOD_AUTO. Default supernodal_switch = 40 */ - - int final_asis ; /* If TRUE, then ignore the other final_* parameters - * (except for final_pack). - * The factor is left as-is when done. Default: TRUE.*/ - - int final_super ; /* If TRUE, leave a factor in supernodal form when - * supernodal factorization is finished. If FALSE, - * then convert to a simplicial factor when done. - * Default: TRUE */ - - int final_ll ; /* If TRUE, leave factor in LL' form when done. - * Otherwise, leave in LDL' form. Default: FALSE */ - - int final_pack ; /* If TRUE, pack the columns when done. If TRUE, and - * cholmod_factorize is called with a symbolic L, L is - * allocated with exactly the space required, using L->ColCount. If you - * plan on modifying the factorization, set Common->final_pack to FALSE, - * and each column will be given a little extra slack space for future - * growth in fill-in due to updates. Default: TRUE */ - - int final_monotonic ; /* If TRUE, ensure columns are monotonic when done. - * Default: TRUE */ - - int final_resymbol ;/* if cholmod_factorize performed a supernodal - * factorization, final_resymbol is true, and - * final_super is FALSE (convert a simplicial numeric factorization), - * then numerically zero entries that resulted from relaxed supernodal - * amalgamation are removed. This does not remove entries that are zero - * due to exact numeric cancellation, since doing so would break the - * update/downdate rowadd/rowdel routines. Default: FALSE. */ - - /* supernodal relaxed amalgamation parameters: */ - double zrelax [3] ; - size_t nrelax [3] ; - - /* Let ns be the total number of columns in two adjacent supernodes. - * Let z be the fraction of zero entries in the two supernodes if they - * are merged (z includes zero entries from prior amalgamations). The - * two supernodes are merged if: - * (ns <= nrelax [0]) || (no new zero entries added) || - * (ns <= nrelax [1] && z < zrelax [0]) || - * (ns <= nrelax [2] && z < zrelax [1]) || (z < zrelax [2]) - * - * Default parameters result in the following rule: - * (ns <= 4) || (no new zero entries added) || - * (ns <= 16 && z < 0.8) || (ns <= 48 && z < 0.1) || (z < 0.05) - */ - - int prefer_zomplex ; /* X = cholmod_solve (sys, L, B, Common) computes - * x=A\b or solves a related system. If L and B are - * both real, then X is real. Otherwise, X is returned as - * CHOLMOD_COMPLEX if Common->prefer_zomplex is FALSE, or - * CHOLMOD_ZOMPLEX if Common->prefer_zomplex is TRUE. This parameter - * is needed because there is no supernodal zomplex L. Suppose the - * caller wants all complex matrices to be stored in zomplex form - * (MATLAB, for example). A supernodal L is returned in complex form - * if A is zomplex. B can be real, and thus X = cholmod_solve (L,B) - * should return X as zomplex. This cannot be inferred from the input - * arguments L and B. Default: FALSE, since all data types are - * supported in CHOLMOD_COMPLEX form and since this is the native type - * of LAPACK and the BLAS. Note that the MATLAB/cholmod.c mexFunction - * sets this parameter to TRUE, since MATLAB matrices are in - * CHOLMOD_ZOMPLEX form. - */ - - int prefer_upper ; /* cholmod_analyze and cholmod_factorize work - * fastest when a symmetric matrix is stored in - * upper triangular form when a fill-reducing ordering is used. In - * MATLAB, this corresponds to how x=A\b works. When the matrix is - * ordered as-is, they work fastest when a symmetric matrix is in lower - * triangular form. In MATLAB, R=chol(A) does the opposite. This - * parameter affects only how cholmod_read returns a symmetric matrix. - * If TRUE (the default case), a symmetric matrix is always returned in - * upper-triangular form (A->stype = 1). */ - - int quick_return_if_not_posdef ; /* if TRUE, the supernodal numeric - * factorization will return quickly if - * the matrix is not positive definite. Default: FALSE. */ - - /* ---------------------------------------------------------------------- */ - /* printing and error handling options */ - /* ---------------------------------------------------------------------- */ - - int print ; /* print level. Default: 3 */ - int precise ; /* if TRUE, print 16 digits. Otherwise print 5 */ - int (*print_function) (const char *, ...) ; /* pointer to printf */ - - int try_catch ; /* if TRUE, then ignore errors; CHOLMOD is in the middle - * of a try/catch block. No error message is printed - * and the Common->error_handler function is not called. */ - - void (*error_handler) (int status, const char *file, - int line, const char *message) ; - - /* Common->error_handler is the user's error handling routine. If not - * NULL, this routine is called if an error occurs in CHOLMOD. status - * can be CHOLMOD_OK (0), negative for a fatal error, and positive for - * a warning. file is a string containing the name of the source code - * file where the error occured, and line is the line number in that - * file. message is a string describing the error in more detail. */ - - /* ---------------------------------------------------------------------- */ - /* ordering options */ - /* ---------------------------------------------------------------------- */ - - /* The cholmod_analyze routine can try many different orderings and select - * the best one. It can also try one ordering method multiple times, with - * different parameter settings. The default is to use three orderings, - * the user's permutation (if provided), AMD which is the fastest ordering - * and generally gives good fill-in, and METIS. CHOLMOD's nested dissection - * (METIS with a constrained AMD) usually gives a better ordering than METIS - * alone (by about 5% to 10%) but it takes more time. - * - * If you know the method that is best for your matrix, set Common->nmethods - * to 1 and set Common->method [0] to the set of parameters for that method. - * If you set it to 1 and do not provide a permutation, then only AMD will - * be called. - * - * If METIS is not available, the default # of methods tried is 2 (the user - * permutation, if any, and AMD). - * - * To try other methods, set Common->nmethods to the number of methods you - * want to try. The suite of default methods and their parameters is - * described in the cholmod_defaults routine, and summarized here: - * - * Common->method [i]: - * i = 0: user-provided ordering (cholmod_analyze_p only) - * i = 1: AMD (for both A and A*A') - * i = 2: METIS - * i = 3: CHOLMOD's nested dissection (NESDIS), default parameters - * i = 4: natural - * i = 5: NESDIS with nd_small = 20000 - * i = 6: NESDIS with nd_small = 4, no constrained minimum degree - * i = 7: NESDIS with no dense node removal - * i = 8: AMD for A, COLAMD for A*A' - * - * You can modify the suite of methods you wish to try by modifying - * Common.method [...] after calling cholmod_start or cholmod_defaults. - * - * For example, to use AMD, followed by a weighted postordering: - * - * Common->nmethods = 1 ; - * Common->method [0].ordering = CHOLMOD_AMD ; - * Common->postorder = TRUE ; - * - * To use the natural ordering (with no postordering): - * - * Common->nmethods = 1 ; - * Common->method [0].ordering = CHOLMOD_NATURAL ; - * Common->postorder = FALSE ; - * - * If you are going to factorize hundreds or more matrices with the same - * nonzero pattern, you may wish to spend a great deal of time finding a - * good permutation. In this case, try setting Common->nmethods to 9. - * The time spent in cholmod_analysis will be very high, but you need to - * call it only once. - * - * cholmod_analyze sets Common->current to a value between 0 and nmethods-1. - * Each ordering method uses the set of options defined by this parameter. - */ - - int nmethods ; /* The number of ordering methods to try. Default: 0. - * nmethods = 0 is a special case. cholmod_analyze - * will try the user-provided ordering (if given) and AMD. Let fl and - * lnz be the flop count and nonzeros in L from AMD's ordering. Let - * anz be the number of nonzeros in the upper or lower triangular part - * of the symmetric matrix A. If fl/lnz < 500 or lnz/anz < 5, then this - * is a good ordering, and METIS is not attempted. Otherwise, METIS is - * tried. The best ordering found is used. If nmethods > 0, the - * methods used are given in the method[ ] array, below. The first - * three methods in the default suite of orderings is (1) use the given - * permutation (if provided), (2) use AMD, and (3) use METIS. Maximum - * allowed value is CHOLMOD_MAXMETHODS. */ - - int current ; /* The current method being tried. Default: 0. Valid - * range is 0 to nmethods-1. */ - - int selected ; /* The best method found. */ - - /* The suite of ordering methods and parameters: */ - - struct cholmod_method_struct - { - /* statistics for this method */ - double lnz ; /* nnz(L) excl. zeros from supernodal amalgamation, - * for a "pure" L */ - - double fl ; /* flop count for a "pure", real simplicial LL' - * factorization, with no extra work due to - * amalgamation. Subtract n to get the LDL' flop count. Multiply - * by about 4 if the matrix is complex or zomplex. */ - - /* ordering method parameters */ - double prune_dense ;/* dense row/col control for AMD, SYMAMD, CSYMAMD, - * and NESDIS (cholmod_nested_dissection). For a - * symmetric n-by-n matrix, rows/columns with more than - * MAX (16, prune_dense * sqrt (n)) entries are removed prior to - * ordering. They appear at the end of the re-ordered matrix. - * - * If prune_dense < 0, only completely dense rows/cols are removed. - * - * This paramater is also the dense column control for COLAMD and - * CCOLAMD. For an m-by-n matrix, columns with more than - * MAX (16, prune_dense * sqrt (MIN (m,n))) entries are removed prior - * to ordering. They appear at the end of the re-ordered matrix. - * CHOLMOD factorizes A*A', so it calls COLAMD and CCOLAMD with A', - * not A. Thus, this parameter affects the dense *row* control for - * CHOLMOD's matrix, and the dense *column* control for COLAMD and - * CCOLAMD. - * - * Removing dense rows and columns improves the run-time of the - * ordering methods. It has some impact on ordering quality - * (usually minimal, sometimes good, sometimes bad). - * - * Default: 10. */ - - double prune_dense2 ;/* dense row control for COLAMD and CCOLAMD. - * Rows with more than MAX (16, dense2 * sqrt (n)) - * for an m-by-n matrix are removed prior to ordering. CHOLMOD's - * matrix is transposed before ordering it with COLAMD or CCOLAMD, - * so this controls the dense *columns* of CHOLMOD's matrix, and - * the dense *rows* of COLAMD's or CCOLAMD's matrix. - * - * If prune_dense2 < 0, only completely dense rows/cols are removed. - * - * Default: -1. Note that this is not the default for COLAMD and - * CCOLAMD. -1 is best for Cholesky. 10 is best for LU. */ - - double nd_oksep ; /* in NESDIS, when a node separator is computed, it - * discarded if nsep >= nd_oksep*n, where nsep is - * the number of nodes in the separator, and n is the size of the - * graph being cut. Valid range is 0 to 1. If 1 or greater, the - * separator is discarded if it consists of the entire graph. - * Default: 1 */ - - double other1 [4] ; /* future expansion */ - - size_t nd_small ; /* do not partition graphs with fewer nodes than - * nd_small, in NESDIS. Default: 200 (same as - * METIS) */ - - size_t other2 [4] ; /* future expansion */ - - int aggressive ; /* Aggresive absorption in AMD, COLAMD, SYMAMD, - * CCOLAMD, and CSYMAMD. Default: TRUE */ - - int order_for_lu ; /* CCOLAMD can be optimized to produce an ordering - * for LU or Cholesky factorization. CHOLMOD only - * performs a Cholesky factorization. However, you may wish to use - * CHOLMOD as an interface for CCOLAMD but use it for your own LU - * factorization. In this case, order_for_lu should be set to FALSE. - * When factorizing in CHOLMOD itself, you should *** NEVER *** set - * this parameter FALSE. Default: TRUE. */ - - int nd_compress ; /* If TRUE, compress the graph and subgraphs before - * partitioning them in NESDIS. Default: TRUE */ - - int nd_camd ; /* If 1, follow the nested dissection ordering - * with a constrained minimum degree ordering that - * respects the partitioning just found (using CAMD). If 2, use - * CSYMAMD instead. If you set nd_small very small, you may not need - * this ordering, and can save time by setting it to zero (no - * constrained minimum degree ordering). Default: 1. */ - - int nd_components ; /* The nested dissection ordering finds a node - * separator that splits the graph into two parts, - * which may be unconnected. If nd_components is TRUE, each of - * these connected components is split independently. If FALSE, - * each part is split as a whole, even if it consists of more than - * one connected component. Default: FALSE */ - - /* fill-reducing ordering to use */ - int ordering ; - - size_t other3 [4] ; /* future expansion */ - - } method [CHOLMOD_MAXMETHODS + 1] ; - - int postorder ; /* If TRUE, cholmod_analyze follows the ordering with a - * weighted postorder of the elimination tree. Improves - * supernode amalgamation. Does not affect fundamental nnz(L) and - * flop count. Default: TRUE. */ - - /* ---------------------------------------------------------------------- */ - /* memory management routines */ - /* ---------------------------------------------------------------------- */ - - void *(*malloc_memory) (size_t) ; /* pointer to malloc */ - void *(*realloc_memory) (void *, size_t) ; /* pointer to realloc */ - void (*free_memory) (void *) ; /* pointer to free */ - void *(*calloc_memory) (size_t, size_t) ; /* pointer to calloc */ - - /* ---------------------------------------------------------------------- */ - /* routines for complex arithmetic */ - /* ---------------------------------------------------------------------- */ - - int (*complex_divide) (double ax, double az, double bx, double bz, - double *cx, double *cz) ; - - /* flag = complex_divide (ax, az, bx, bz, &cx, &cz) computes the complex - * division c = a/b, where ax and az hold the real and imaginary part - * of a, and b and c are stored similarly. flag is returned as 1 if - * a divide-by-zero occurs, or 0 otherwise. By default, the function - * pointer Common->complex_divide is set equal to cholmod_divcomplex. - */ - - double (*hypotenuse) (double x, double y) ; - - /* s = hypotenuse (x,y) computes s = sqrt (x*x + y*y), but does so more - * accurately. By default, the function pointer Common->hypotenuse is - * set equal to cholmod_hypot. See also the hypot function in the C99 - * standard, which has an identical syntax and function. If you have - * a C99-compliant compiler, you can set Common->hypotenuse = hypot. */ - - /* ---------------------------------------------------------------------- */ - /* METIS workarounds */ - /* ---------------------------------------------------------------------- */ - - double metis_memory ; /* This is a parameter for CHOLMOD's interface to - * METIS, not a parameter to METIS itself. METIS - * uses an amount of memory that is difficult to estimate precisely - * beforehand. If it runs out of memory, it terminates your program. - * All routines in CHOLMOD except for CHOLMOD's interface to METIS - * return an error status and safely return to your program if they run - * out of memory. To mitigate this problem, the CHOLMOD interface - * can allocate a single block of memory equal in size to an empirical - * upper bound of METIS's memory usage times the Common->metis_memory - * parameter, and then immediately free it. It then calls METIS. If - * this pre-allocation fails, it is possible that METIS will fail as - * well, and so CHOLMOD returns with an out-of-memory condition without - * calling METIS. - * - * METIS_NodeND (used in the CHOLMOD_METIS ordering option) with its - * default parameter settings typically uses about (4*nz+40n+4096) - * times sizeof(int) memory, where nz is equal to the number of entries - * in A for the symmetric case or AA' if an unsymmetric matrix is - * being ordered (where nz includes both the upper and lower parts - * of A or AA'). The observed "upper bound" (with 2 exceptions), - * measured in an instrumented copy of METIS 4.0.1 on thousands of - * matrices, is (10*nz+50*n+4096) * sizeof(int). Two large matrices - * exceeded this bound, one by almost a factor of 2 (Gupta/gupta2). - * - * If your program is terminated by METIS, try setting metis_memory to - * 2.0, or even higher if needed. By default, CHOLMOD assumes that METIS - * does not have this problem (so that CHOLMOD will work correctly when - * this issue is fixed in METIS). Thus, the default value is zero. - * This work-around is not guaranteed anyway. - * - * If a matrix exceeds this predicted memory usage, AMD is attempted - * instead. It, too, may run out of memory, but if it does so it will - * not terminate your program. - */ +#include - double metis_dswitch ; /* METIS_NodeND in METIS 4.0.1 gives a seg */ - size_t metis_nswitch ; /* fault with one matrix of order n = 3005 and - * nz = 6,036,025. This is a very dense graph. - * The workaround is to use AMD instead of METIS for matrices of dimension - * greater than Common->metis_nswitch (default 3000) or more and with - * density of Common->metis_dswitch (default 0.66) or more. - * cholmod_nested_dissection has no problems with the same matrix, even - * though it uses METIS_NodeComputeSeparator on this matrix. If this - * seg fault does not affect you, set metis_nswitch to zero or less, - * and CHOLMOD will not switch to AMD based just on the density of the - * matrix (it will still switch to AMD if the metis_memory parameter - * causes the switch). - */ - - /* ---------------------------------------------------------------------- */ - /* workspace */ - /* ---------------------------------------------------------------------- */ - - /* CHOLMOD has several routines that take less time than the size of - * workspace they require. Allocating and initializing the workspace would - * dominate the run time, unless workspace is allocated and initialized - * just once. CHOLMOD allocates this space when needed, and holds it here - * between calls to CHOLMOD. cholmod_start sets these pointers to NULL - * (which is why it must be the first routine called in CHOLMOD). - * cholmod_finish frees the workspace (which is why it must be the last - * call to CHOLMOD). - */ - - size_t nrow ; /* size of Flag and Head */ - SuiteSparse_long mark ; /* mark value for Flag array */ - size_t iworksize ; /* size of Iwork. Upper bound: 6*nrow+ncol */ - size_t xworksize ; /* size of Xwork, in bytes. - * maxrank*nrow*sizeof(double) for update/downdate. - * 2*nrow*sizeof(double) otherwise */ - - /* initialized workspace: contents needed between calls to CHOLMOD */ - void *Flag ; /* size nrow, an integer array. Kept cleared between - * calls to cholmod rouines (Flag [i] < mark) */ - - void *Head ; /* size nrow+1, an integer array. Kept cleared between - * calls to cholmod routines (Head [i] = EMPTY) */ - - void *Xwork ; /* a double array. Its size varies. It is nrow for - * most routines (cholmod_rowfac, cholmod_add, - * cholmod_aat, cholmod_norm, cholmod_ssmult) for the real case, twice - * that when the input matrices are complex or zomplex. It is of size - * 2*nrow for cholmod_rowadd and cholmod_rowdel. For cholmod_updown, - * its size is maxrank*nrow where maxrank is 2, 4, or 8. Kept cleared - * between calls to cholmod (set to zero). */ - - /* uninitialized workspace, contents not needed between calls to CHOLMOD */ - void *Iwork ; /* size iworksize, 2*nrow+ncol for most routines, - * up to 6*nrow+ncol for cholmod_analyze. */ - - int itype ; /* If CHOLMOD_LONG, Flag, Head, and Iwork are SuiteSparse_long. - * Otherwise all three arrays are int. */ - - int dtype ; /* double or float */ - - /* Common->itype and Common->dtype are used to define the types of all - * sparse matrices, triplet matrices, dense matrices, and factors - * created using this Common struct. The itypes and dtypes of all - * parameters to all CHOLMOD routines must match. */ - - int no_workspace_reallocate ; /* this is an internal flag, used as a - * precaution by cholmod_analyze. It is normally false. If true, - * cholmod_allocate_work is not allowed to reallocate any workspace; - * they must use the existing workspace in Common (Iwork, Flag, Head, - * and Xwork). Added for CHOLMOD v1.1 */ - - /* ---------------------------------------------------------------------- */ - /* statistics */ - /* ---------------------------------------------------------------------- */ - - /* fl and lnz are set only in cholmod_analyze and cholmod_rowcolcounts, - * in the Cholesky modudle. modfl is set only in the Modify module. */ - - int status ; /* error code */ - double fl ; /* LL' flop count from most recent analysis */ - double lnz ; /* fundamental nz in L */ - double anz ; /* nonzeros in tril(A) if A is symmetric/lower, - * triu(A) if symmetric/upper, or tril(A*A') if - * unsymmetric, in last call to cholmod_analyze. */ - double modfl ; /* flop count from most recent update/downdate/ - * rowadd/rowdel (excluding flops to modify the - * solution to Lx=b, if computed) */ - size_t malloc_count ; /* # of objects malloc'ed minus the # free'd*/ - size_t memory_usage ; /* peak memory usage in bytes */ - size_t memory_inuse ; /* current memory usage in bytes */ - - double nrealloc_col ; /* # of column reallocations */ - double nrealloc_factor ;/* # of factor reallocations due to col. reallocs */ - double ndbounds_hit ; /* # of times diagonal modified by dbound */ - - double rowfacfl ; /* # of flops in last call to cholmod_rowfac */ - double aatfl ; /* # of flops to compute A(:,f)*A(:,f)' */ - - /* ---------------------------------------------------------------------- */ - /* future expansion */ - /* ---------------------------------------------------------------------- */ - - /* To allow CHOLMOD to be updated without recompiling the user application, - * additional space is set aside here for future statistics, parameters, - * and workspace. Note: additional entries were added in v1.1 to the - * method array, above, and thus v1.0 and v1.1 are not binary compatible. - * - * v1.1 to the current version are binary compatible. - */ - - /* ---------------------------------------------------------------------- */ - double other1 [10] ; - - double SPQR_xstat [4] ; /* for SuiteSparseQR statistics */ - - /* SuiteSparseQR control parameters: */ - double SPQR_grain ; /* task size is >= max (total flops / grain) */ - double SPQR_small ; /* task size is >= small */ - - /* ---------------------------------------------------------------------- */ - SuiteSparse_long SPQR_istat [10] ; /* for SuiteSparseQR statistics */ - SuiteSparse_long other2 [6] ; /* reduced from size 16 in v1.6 */ - - /* ---------------------------------------------------------------------- */ - int other3 [10] ; /* reduced from size 16 in v1.1. */ - - int prefer_binary ; /* cholmod_read_triplet converts a symmetric - * pattern-only matrix into a real matrix. If - * prefer_binary is FALSE, the diagonal entries are set to 1 + the degree - * of the row/column, and off-diagonal entries are set to -1 (resulting - * in a positive definite matrix if the diagonal is zero-free). Most - * symmetric patterns are the pattern a positive definite matrix. If - * this parameter is TRUE, then the matrix is returned with a 1 in each - * entry, instead. Default: FALSE. Added in v1.3. */ - - /* control parameter (added for v1.2): */ - int default_nesdis ; /* Default: FALSE. If FALSE, then the default - * ordering strategy (when Common->nmethods == 0) - * is to try the given ordering (if present), AMD, and then METIS if AMD - * reports high fill-in. If Common->default_nesdis is TRUE then NESDIS - * is used instead in the default strategy. */ - - /* statistic (added for v1.2): */ - int called_nd ; /* TRUE if the last call to - * cholmod_analyze called NESDIS or METIS. */ - - int blas_ok ; /* FALSE if BLAS int overflow; TRUE otherwise */ - - /* SuiteSparseQR control parameters: */ - int SPQR_shrink ; /* controls stack realloc method */ - int SPQR_nthreads ; /* number of TBB threads, 0 = auto */ - - /* ---------------------------------------------------------------------- */ - size_t other4 [16] ; - - /* ---------------------------------------------------------------------- */ - void *other5 [16] ; - -} cholmod_common ; - -// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1170 : -/* A sparse matrix stored in compressed-column form. */ - -typedef struct cholmod_sparse_struct -{ - size_t nrow ; /* the matrix is nrow-by-ncol */ - size_t ncol ; - size_t nzmax ; /* maximum number of entries in the matrix */ - - /* pointers to int or SuiteSparse_long: */ - void *p ; /* p [0..ncol], the column pointers */ - void *i ; /* i [0..nzmax-1], the row indices */ - - /* for unpacked matrices only: */ - void *nz ; /* nz [0..ncol-1], the # of nonzeros in each col. In - * packed form, the nonzero pattern of column j is in - * A->i [A->p [j] ... A->p [j+1]-1]. In unpacked form, column j is in - * A->i [A->p [j] ... A->p [j]+A->nz[j]-1] instead. In both cases, the - * numerical values (if present) are in the corresponding locations in - * the array x (or z if A->xtype is CHOLMOD_ZOMPLEX). */ - - /* pointers to double or float: */ - void *x ; /* size nzmax or 2*nzmax, if present */ - void *z ; /* size nzmax, if present */ - - int stype ; /* Describes what parts of the matrix are considered: - * - * 0: matrix is "unsymmetric": use both upper and lower triangular parts - * (the matrix may actually be symmetric in pattern and value, but - * both parts are explicitly stored and used). May be square or - * rectangular. - * >0: matrix is square and symmetric, use upper triangular part. - * Entries in the lower triangular part are ignored. - * <0: matrix is square and symmetric, use lower triangular part. - * Entries in the upper triangular part are ignored. - * - * Note that stype>0 and stype<0 are different for cholmod_sparse and - * cholmod_triplet. See the cholmod_triplet data structure for more - * details. - */ - - int itype ; /* CHOLMOD_INT: p, i, and nz are int. - * CHOLMOD_INTLONG: p is SuiteSparse_long, - * i and nz are int. - * CHOLMOD_LONG: p, i, and nz are SuiteSparse_long */ - - int xtype ; /* pattern, real, complex, or zomplex */ - int dtype ; /* x and z are double or float */ - int sorted ; /* TRUE if columns are sorted, FALSE otherwise */ - int packed ; /* TRUE if packed (nz ignored), FALSE if unpacked - * (nz is required) */ - -} cholmod_sparse ; - -// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1552 : -/* A symbolic and numeric factorization, either simplicial or supernodal. - * In all cases, the row indices in the columns of L are kept sorted. */ - -typedef struct cholmod_factor_struct -{ - /* ---------------------------------------------------------------------- */ - /* for both simplicial and supernodal factorizations */ - /* ---------------------------------------------------------------------- */ - - size_t n ; /* L is n-by-n */ - - size_t minor ; /* If the factorization failed, L->minor is the column - * at which it failed (in the range 0 to n-1). A value - * of n means the factorization was successful or - * the matrix has not yet been factorized. */ - - /* ---------------------------------------------------------------------- */ - /* symbolic ordering and analysis */ - /* ---------------------------------------------------------------------- */ - - void *Perm ; /* size n, permutation used */ - void *ColCount ; /* size n, column counts for simplicial L */ - - void *IPerm ; /* size n, inverse permutation. Only created by - * cholmod_solve2 if Bset is used. */ - - /* ---------------------------------------------------------------------- */ - /* simplicial factorization */ - /* ---------------------------------------------------------------------- */ - - size_t nzmax ; /* size of i and x */ - - void *p ; /* p [0..ncol], the column pointers */ - void *i ; /* i [0..nzmax-1], the row indices */ - void *x ; /* x [0..nzmax-1], the numerical values */ - void *z ; - void *nz ; /* nz [0..ncol-1], the # of nonzeros in each column. - * i [p [j] ... p [j]+nz[j]-1] contains the row indices, - * and the numerical values are in the same locatins - * in x. The value of i [p [k]] is always k. */ - - void *next ; /* size ncol+2. next [j] is the next column in i/x */ - void *prev ; /* size ncol+2. prev [j] is the prior column in i/x. - * head of the list is ncol+1, and the tail is ncol. */ - - /* ---------------------------------------------------------------------- */ - /* supernodal factorization */ - /* ---------------------------------------------------------------------- */ - - /* Note that L->x is shared with the simplicial data structure. L->x has - * size L->nzmax for a simplicial factor, and size L->xsize for a supernodal - * factor. */ - - size_t nsuper ; /* number of supernodes */ - size_t ssize ; /* size of s, integer part of supernodes */ - size_t xsize ; /* size of x, real part of supernodes */ - size_t maxcsize ; /* size of largest update matrix */ - size_t maxesize ; /* max # of rows in supernodes, excl. triangular part */ - - void *super ; /* size nsuper+1, first col in each supernode */ - void *pi ; /* size nsuper+1, pointers to integer patterns */ - void *px ; /* size nsuper+1, pointers to real parts */ - void *s ; /* size ssize, integer part of supernodes */ - - /* ---------------------------------------------------------------------- */ - /* factorization type */ - /* ---------------------------------------------------------------------- */ - - int ordering ; /* ordering method used */ - - int is_ll ; /* TRUE if LL', FALSE if LDL' */ - int is_super ; /* TRUE if supernodal, FALSE if simplicial */ - int is_monotonic ; /* TRUE if columns of L appear in order 0..n-1. - * Only applicable to simplicial numeric types. */ - - /* There are 8 types of factor objects that cholmod_factor can represent - * (only 6 are used): - * - * Numeric types (xtype is not CHOLMOD_PATTERN) - * -------------------------------------------- - * - * simplicial LDL': (is_ll FALSE, is_super FALSE). Stored in compressed - * column form, using the simplicial components above (nzmax, p, i, - * x, z, nz, next, and prev). The unit diagonal of L is not stored, - * and D is stored in its place. There are no supernodes. - * - * simplicial LL': (is_ll TRUE, is_super FALSE). Uses the same storage - * scheme as the simplicial LDL', except that D does not appear. - * The first entry of each column of L is the diagonal entry of - * that column of L. - * - * supernodal LDL': (is_ll FALSE, is_super TRUE). Not used. - * FUTURE WORK: add support for supernodal LDL' - * - * supernodal LL': (is_ll TRUE, is_super TRUE). A supernodal factor, - * using the supernodal components described above (nsuper, ssize, - * xsize, maxcsize, maxesize, super, pi, px, s, x, and z). - * - * - * Symbolic types (xtype is CHOLMOD_PATTERN) - * ----------------------------------------- - * - * simplicial LDL': (is_ll FALSE, is_super FALSE). Nothing is present - * except Perm and ColCount. - * - * simplicial LL': (is_ll TRUE, is_super FALSE). Identical to the - * simplicial LDL', except for the is_ll flag. - * - * supernodal LDL': (is_ll FALSE, is_super TRUE). Not used. - * FUTURE WORK: add support for supernodal LDL' - * - * supernodal LL': (is_ll TRUE, is_super TRUE). A supernodal symbolic - * factorization. The simplicial symbolic information is present - * (Perm and ColCount), as is all of the supernodal factorization - * except for the numerical values (x and z). - */ - - int itype ; /* The integer arrays are Perm, ColCount, p, i, nz, - * next, prev, super, pi, px, and s. If itype is - * CHOLMOD_INT, all of these are int arrays. - * CHOLMOD_INTLONG: p, pi, px are SuiteSparse_long, others int. - * CHOLMOD_LONG: all integer arrays are SuiteSparse_long. */ - int xtype ; /* pattern, real, complex, or zomplex */ - int dtype ; /* x and z double or float */ - -} cholmod_factor ; - -// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1836 : -/* A dense matrix in column-oriented form. It has no itype since it contains - * no integers. Entry in row i and column j is located in x [i+j*d]. - */ - -typedef struct cholmod_dense_struct -{ - size_t nrow ; /* the matrix is nrow-by-ncol */ - size_t ncol ; - size_t nzmax ; /* maximum number of entries in the matrix */ - size_t d ; /* leading dimension (d >= nrow must hold) */ - void *x ; /* size nzmax or 2*nzmax, if present */ - void *z ; /* size nzmax, if present */ - int xtype ; /* pattern, real, complex, or zomplex */ - int dtype ; /* x and z double or float */ - -} cholmod_dense ; - -// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 2033 : -/* A sparse matrix stored in triplet form. */ - -typedef struct cholmod_triplet_struct -{ - size_t nrow ; /* the matrix is nrow-by-ncol */ - size_t ncol ; - size_t nzmax ; /* maximum number of entries in the matrix */ - size_t nnz ; /* number of nonzeros in the matrix */ - - void *i ; /* i [0..nzmax-1], the row indices */ - void *j ; /* j [0..nzmax-1], the column indices */ - void *x ; /* size nzmax or 2*nzmax, if present */ - void *z ; /* size nzmax, if present */ - - int stype ; /* Describes what parts of the matrix are considered: - * - * 0: matrix is "unsymmetric": use both upper and lower triangular parts - * (the matrix may actually be symmetric in pattern and value, but - * both parts are explicitly stored and used). May be square or - * rectangular. - * >0: matrix is square and symmetric. Entries in the lower triangular - * part are transposed and added to the upper triangular part when - * the matrix is converted to cholmod_sparse form. - * <0: matrix is square and symmetric. Entries in the upper triangular - * part are transposed and added to the lower triangular part when - * the matrix is converted to cholmod_sparse form. - * - * Note that stype>0 and stype<0 are different for cholmod_sparse and - * cholmod_triplet. The reason is simple. You can permute a symmetric - * triplet matrix by simply replacing a row and column index with their - * new row and column indices, via an inverse permutation. Suppose - * P = L->Perm is your permutation, and Pinv is an array of size n. - * Suppose a symmetric matrix A is represent by a triplet matrix T, with - * entries only in the upper triangular part. Then the following code: - * - * Ti = T->i ; - * Tj = T->j ; - * for (k = 0 ; k < n ; k++) Pinv [P [k]] = k ; - * for (k = 0 ; k < nz ; k++) Ti [k] = Pinv [Ti [k]] ; - * for (k = 0 ; k < nz ; k++) Tj [k] = Pinv [Tj [k]] ; - * - * creates the triplet form of C=P*A*P'. However, if T initially - * contains just the upper triangular entries (T->stype = 1), after - * permutation it has entries in both the upper and lower triangular - * parts. These entries should be transposed when constructing the - * cholmod_sparse form of A, which is what cholmod_triplet_to_sparse - * does. Thus: - * - * C = cholmod_triplet_to_sparse (T, 0, &Common) ; - * - * will return the matrix C = P*A*P'. - * - * Since the triplet matrix T is so simple to generate, it's quite easy - * to remove entries that you do not want, prior to converting T to the - * cholmod_sparse form. So if you include these entries in T, CHOLMOD - * assumes that there must be a reason (such as the one above). Thus, - * no entry in a triplet matrix is ever ignored. - */ - - int itype ; /* CHOLMOD_LONG: i and j are SuiteSparse_long. Otherwise int */ - int xtype ; /* pattern, real, complex, or zomplex */ - int dtype ; /* x and z are double or float */ - -} cholmod_triplet ; - -// -------- our (Matrix) short and const_ forms of of the pointers : -typedef cholmod_common* CHM_CM; -typedef cholmod_dense* CHM_DN; -typedef const cholmod_dense* const_CHM_DN; -typedef cholmod_factor* CHM_FR; -typedef const cholmod_factor* const_CHM_FR; -typedef cholmod_sparse* CHM_SP; -typedef const cholmod_sparse* const_CHM_SP; -typedef cholmod_triplet* CHM_TR; -typedef const cholmod_triplet* const_CHM_TR; - - -// --------- Matrix ("M_") R ("R_") pkg routines "re-exported": --------------- - -int cholmod_start(CHM_CM); -void R_cholmod_error(int status, const char *file, int line, const char *message); -int cholmod_finish(CHM_CM); - -CHM_SP cholmod_allocate_sparse(size_t nrow, size_t ncol, - size_t nzmax, int sorted, - int packed, int stype, int xtype, - CHM_CM); -int cholmod_free_factor(CHM_FR *L, CHM_CM); -int cholmod_free_dense(CHM_DN *A, CHM_CM); -int cholmod_free_sparse(CHM_SP *A, CHM_CM); -int cholmod_free_triplet(CHM_TR *T, CHM_CM); - -long int cholmod_nnz(const_CHM_SP, CHM_CM); -CHM_SP cholmod_speye(size_t nrow, size_t ncol, int xtype, CHM_CM); -CHM_SP cholmod_transpose(const_CHM_SP, int values, CHM_CM); -int cholmod_sort(CHM_SP A, CHM_CM); -CHM_SP cholmod_vertcat(const_CHM_SP, const_CHM_SP, int values, CHM_CM); -CHM_SP cholmod_copy(const_CHM_SP, int stype, int mode, CHM_CM); -CHM_SP cholmod_add(const_CHM_SP, const_CHM_SP, double alpha [2], double beta [2], - int values, int sorted, CHM_CM); - -// from ../../src/CHOLMOD/Include/cholmod_cholesky.h - line 178 : -#define CHOLMOD_A 0 /* solve Ax=b */ -#define CHOLMOD_LDLt 1 /* solve LDL'x=b */ -#define CHOLMOD_LD 2 /* solve LDx=b */ -#define CHOLMOD_DLt 3 /* solve DL'x=b */ -#define CHOLMOD_L 4 /* solve Lx=b */ -#define CHOLMOD_Lt 5 /* solve L'x=b */ -#define CHOLMOD_D 6 /* solve Dx=b */ -#define CHOLMOD_P 7 /* permute x=Px */ -#define CHOLMOD_Pt 8 /* permute x=P'x */ - -CHM_DN cholmod_solve(int, const_CHM_FR, const_CHM_DN, CHM_CM); -CHM_SP cholmod_spsolve(int, const_CHM_FR, const_CHM_SP, CHM_CM); -int cholmod_sdmult(const_CHM_SP, int, const double*, const double*, - const_CHM_DN, CHM_DN Y, CHM_CM); -CHM_SP cholmod_ssmult(const_CHM_SP, const_CHM_SP, int, int, int, - CHM_CM); -int cholmod_factorize(const_CHM_SP, CHM_FR L, CHM_CM); -int cholmod_factorize_p(const_CHM_SP, double *beta, int *fset, - size_t fsize, CHM_FR L, CHM_CM); -CHM_SP cholmod_copy_sparse(const_CHM_SP, CHM_CM); -CHM_DN cholmod_copy_dense(const_CHM_DN, CHM_CM); -CHM_SP cholmod_aat(const_CHM_SP, int *fset, size_t fsize, int mode, - CHM_CM); -int cholmod_band_inplace(CHM_SP A, int k1, int k2, int mode, CHM_CM); -CHM_SP cholmod_add(const_CHM_SP, const_CHM_SP, double alpha[2], double beta[2], - int values, int sorted, CHM_CM); -CHM_DN cholmod_allocate_dense(size_t nrow, size_t ncol, size_t d, - int xtype, CHM_CM); -CHM_FR cholmod_analyze(const_CHM_SP, CHM_CM); -CHM_FR cholmod_analyze_p(const_CHM_SP, int *Perm, int *fset, - size_t fsize, CHM_CM); -int cholmod_change_factor(int to_xtype, int to_ll, int to_super, - int to_packed, int to_monotonic, - CHM_FR L, CHM_CM); -CHM_FR cholmod_copy_factor(const_CHM_FR, CHM_CM); -CHM_SP cholmod_factor_to_sparse(const_CHM_FR, CHM_CM); -CHM_SP cholmod_dense_to_sparse(const_CHM_DN, int values, CHM_CM); -int cholmod_defaults (CHM_CM); -CHM_SP cholmod_triplet_to_sparse(const cholmod_triplet*, int nzmax, CHM_CM); -CHM_SP cholmod_submatrix(const_CHM_SP, int *rset, int rsize, int *cset, - int csize, int values, int sorted, - CHM_CM); -CHM_TR cholmod_sparse_to_triplet(const_CHM_SP, CHM_CM); -CHM_DN cholmod_sparse_to_dense(const_CHM_SP, CHM_CM); -CHM_TR cholmod_allocate_triplet (size_t nrow, size_t ncol, size_t nzmax, - int stype, int xtype, CHM_CM); - -// from ../../src/CHOLMOD/Include/cholmod_matrixops.h - line 107 : -/* scaling modes, selected by the scale input parameter: */ -#define CHOLMOD_SCALAR 0 /* A = s*A */ -#define CHOLMOD_ROW 1 /* A = diag(s)*A */ -#define CHOLMOD_COL 2 /* A = A*diag(s) */ -#define CHOLMOD_SYM 3 /* A = diag(s)*A*diag(s) */ - -int cholmod_scale(const_CHM_DN, int scale, CHM_SP, CHM_CM); - -// added in the Matrix package - the log of the determinant of the matrix that was factored -double chm_factor_ldetL2(const_CHM_FR); - -#ifdef __cplusplus -} #endif - -#endif /* RcppEigen_CHOLMOD_H */ diff --git a/inst/include/RcppEigenForward.h b/inst/include/RcppEigenForward.h index b5ecf289..ef347e4c 100644 --- a/inst/include/RcppEigenForward.h +++ b/inst/include/RcppEigenForward.h @@ -25,11 +25,8 @@ #include #include #include -#include -#include #include #include -#include //#include // causes problems redefining sign #include #include @@ -42,10 +39,6 @@ /* forward declarations */ namespace Rcpp { - /* support for wrap */ - - template - SEXP wrap(const Eigen::CholmodDecomposition >& obj); namespace traits { diff --git a/inst/include/RcppEigenStubs.cpp b/inst/include/RcppEigenStubs.cpp new file mode 100644 index 00000000..ee8b20b0 --- /dev/null +++ b/inst/include/RcppEigenStubs.cpp @@ -0,0 +1,43 @@ +// RcppEigenStubs.cpp: Definitions for CHOLMOD stubs declared +// in RcppEigenCholmod.h, which packages including the header +// must compile +// +// Copyright (C) 2011 Douglas Bates and Martin Maechler +// +// This file is part of RcppEigen. +// +// RcppEigen is free software: 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. +// +// RcppEigen 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with RcppEigen. If not, see . + +/* MJ: Packages could do the following themselves, but that would */ +/* destroy the illusion that they are only linking RcppEigen. */ + +/* MJ: Matrix <= 1.6-1.1 used 'error' and 'warning' unsafely ... sigh ... */ +#include +#ifndef error +# define error Rf_error +# define __ERROR__WAS__UNDEFINED__ +#endif +#ifndef warning +# define warning Rf_warning +# define __WARNING__WAS__UNDEFINED__ +#endif +#include +#ifdef __ERROR__WAS__UNDEFINED__ +# undef error +# undef __ERROR__WAS__UNDEFINED__ +#endif +#ifdef __WARNING__WAS__UNDEFINED__ +# undef warning +# undef __WARNING__WAS__UNDEFINED__ +#endif diff --git a/inst/include/RcppEigenStubs.h b/inst/include/RcppEigenStubs.h deleted file mode 100644 index 04d54184..00000000 --- a/inst/include/RcppEigenStubs.h +++ /dev/null @@ -1,585 +0,0 @@ -// RcppEigenStubs.cpp: Provide access to compiled CHOLMOD functions in -// the Matrix package. -// -// Copyright (C) 2011 Douglas Bates and Martin Maechler -// -// This file is part of RcppEigen. -// -// RcppEigen is free software: 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. -// -// RcppEigen 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 General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with RcppEigen. If not, see . - -// Yes, this really is a C++ source file in an include directory. To -// use the Cholmod support functions in RcppEigen you should create a -// source file, say MyPackage/src/local_stubs.c that contains the -// single line #include "RcppEigenStubs.cpp" - -#include -#include -#include - -#ifdef __cplusplus -extern "C" { -#endif - -inline CHM_DN as_cholmod_dense(CHM_DN ans, SEXP x) -{ - static CHM_DN(*fun)(CHM_DN,SEXP) = NULL; - if(fun == NULL) - fun = (CHM_DN(*)(CHM_DN,SEXP)) - R_GetCCallable("Matrix", "as_cholmod_dense"); - return fun(ans, x); -} - -inline CHM_FR as_cholmod_factor(CHM_FR ans, SEXP x) -{ - static CHM_FR(*fun)(CHM_FR,SEXP) = NULL; - if(fun == NULL) - fun = (CHM_FR(*)(CHM_FR,SEXP)) - R_GetCCallable("Matrix", "as_cholmod_factor"); - return fun(ans, x); -} - -inline CHM_SP as_cholmod_sparse(CHM_SP ans, SEXP x, Rboolean check_Udiag, Rboolean sort_in_place) -{ - static CHM_SP(*fun)(CHM_SP,SEXP,Rboolean,Rboolean)= NULL; - if(fun == NULL) - fun = (CHM_SP(*)(CHM_SP,SEXP,Rboolean,Rboolean)) - R_GetCCallable("Matrix", "as_cholmod_sparse"); - return fun(ans, x, check_Udiag, sort_in_place); -} - -inline CHM_TR as_cholmod_triplet(CHM_TR ans, SEXP x, Rboolean check_Udiag) -{ - static CHM_TR(*fun)(CHM_TR,SEXP,Rboolean)= NULL; - if(fun == NULL) - fun = (CHM_TR(*)(CHM_TR,SEXP,Rboolean)) - R_GetCCallable("Matrix", "as_cholmod_triplet"); - return fun(ans, x, check_Udiag); -} - -inline SEXP Csparse_diagU2N(SEXP x) -{ - static SEXP(*fun)(SEXP) = NULL; - if(fun == NULL) - fun = (SEXP(*)(SEXP)) - R_GetCCallable("Matrix", "Csparse_diagU2N"); - return fun(x); -} - -inline SEXP -M_chm_factor_to_SEXP(const_CHM_FR f, int dofree) -{ - static SEXP(*fun)(const_CHM_FR,int) = NULL; - if(fun == NULL) - fun = (SEXP(*)(const_CHM_FR,int)) - R_GetCCallable("Matrix", "chm_factor_to_SEXP"); - return fun(f, dofree); -} - -inline double -M_chm_factor_ldetL2(const_CHM_FR f) -{ - static double(*fun)(const_CHM_FR) = NULL; - if(fun == NULL) - fun = (double(*)(const_CHM_FR)) - R_GetCCallable("Matrix", "chm_factor_ldetL2"); - return fun(f); -} - -inline CHM_FR -M_chm_factor_update(CHM_FR f, const_CHM_SP A, double mult) -{ - static CHM_FR(*fun)(CHM_FR,const_CHM_SP,double) = NULL; - if(fun == NULL) - fun = (CHM_FR(*)(CHM_FR,const_CHM_SP,double)) - R_GetCCallable("Matrix", "chm_factor_update"); - return fun(f, A, mult); -} - -inline SEXP -M_chm_sparse_to_SEXP(const_CHM_SP a, int dofree, - int uploT, int Rkind, const char *diag, SEXP dn) -{ - static SEXP(*fun)(const_CHM_SP,int,int,int,const char*,SEXP) = NULL; - if(fun == NULL) - fun = (SEXP(*)(const_CHM_SP,int,int,int,const char*,SEXP)) - R_GetCCallable("Matrix", "chm_sparse_to_SEXP"); - return fun(a, dofree, uploT, Rkind, diag, dn); -} - -inline SEXP -M_chm_triplet_to_SEXP(const CHM_TR a, int dofree, - int uploT, int Rkind, const char *diag, SEXP dn) -{ - static SEXP(*fun)(const CHM_TR,int,int,int,const char*,SEXP) = NULL; - if(fun == NULL) - fun = (SEXP(*)(const CHM_TR,int,int,int,const char*,SEXP)) - R_GetCCallable("Matrix", "chm_triplet_to_SEXP"); - return fun(a, dofree, uploT, Rkind, diag, dn); -} - -inline CHM_SP -cholmod_aat(const_CHM_SP A, int *fset, size_t fsize, - int mode, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,int*,size_t, - int,CHM_CM) = NULL; - if(fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,int*,size_t, - int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_aat"); - return fun(A, fset, fsize, mode, Common); -} - -inline int -M_cholmod_band_inplace(CHM_SP A, int k1, int k2, int mode, - CHM_CM Common) -{ - static int(*fun)(CHM_SP,int,int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_SP,int,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_band_inplace"); - return fun(A, k1, k2, mode, Common); -} - -inline CHM_SP -cholmod_add(const_CHM_SP A, const_CHM_SP B, - double alpha[2], double beta[2], int values, - int sorted, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,const_CHM_SP, - double*,double*,int,int, - CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,const_CHM_SP, - double*,double*,int,int, - CHM_CM)) - R_GetCCallable("Matrix", "cholmod_add"); - return fun(A, B, alpha, beta, values, sorted, Common); -} - -inline CHM_DN -cholmod_allocate_dense(size_t nrow, size_t ncol, size_t d, - int xtype, CHM_CM Common) -{ - static CHM_DN(*fun)(size_t,size_t,size_t, - int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_DN(*)(size_t,size_t,size_t, - int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_allocate_dense"); - return fun(nrow, ncol, d, xtype, Common); -} - -inline CHM_SP -cholmod_allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, - int sorted, int packed, int stype, - int xtype, CHM_CM Common) -{ - static CHM_SP(*fun)(size_t,size_t,size_t,int,int, - int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*) - (size_t,size_t,size_t,int,int,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_allocate_sparse"); - return fun(nrow,ncol,nzmax,sorted,packed,stype,xtype,Common); -} - -inline CHM_TR -cholmod_allocate_triplet(size_t nrow, size_t ncol, size_t nzmax, - int stype, int xtype, CHM_CM Common) -{ - static CHM_TR(*fun)(size_t,size_t,size_t, int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_TR(*)(size_t,size_t,size_t,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_allocate_triplet"); - return fun(nrow,ncol,nzmax,stype,xtype,Common); -} - -inline CHM_SP -cholmod_triplet_to_sparse(const cholmod_triplet* T, int nzmax, - CHM_CM Common) -{ - static CHM_SP(*fun)(const cholmod_triplet*,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const cholmod_triplet*,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_triplet_to_sparse"); - return fun(T, nzmax, Common); -} - -inline CHM_TR -cholmod_sparse_to_triplet(const_CHM_SP A, CHM_CM Common) -{ - static CHM_TR(*fun)(const_CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_TR(*)(const_CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_sparse_to_triplet"); - return fun(A, Common); -} - -inline CHM_DN -cholmod_sparse_to_dense(const_CHM_SP A, CHM_CM Common) -{ - static CHM_DN(*fun)(const_CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_DN(*)(const_CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_sparse_to_dense"); - return fun(A, Common); -} - -inline CHM_FR -cholmod_analyze(const_CHM_SP A, CHM_CM Common) -{ - static CHM_FR(*fun)(const_CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_FR(*)(const_CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_analyze"); - return fun(A, Common); -} - -inline CHM_FR -cholmod_analyze_p(const_CHM_SP A, int *Perm, int *fset, - size_t fsize, CHM_CM Common) { - static CHM_FR(*fun)(const_CHM_SP,int*,int*,size_t, - CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_FR(*)(const_CHM_SP,int*,int*, - size_t,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_analyze_p"); - return fun(A, Perm, fset, fsize, Common); -} - -inline CHM_SP -cholmod_copy(const_CHM_SP A, int stype, - int mode, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_copy"); - return fun(A, stype, mode, Common); -} - -inline CHM_DN -cholmod_copy_dense(const_CHM_DN A, CHM_CM Common) -{ - static CHM_DN(*fun)(const_CHM_DN,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_DN(*)(const_CHM_DN,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_copy_dense"); - return fun(A, Common); -} - -inline CHM_FR -cholmod_copy_factor(const_CHM_FR L, CHM_CM Common) -{ - static CHM_FR(*fun)(const_CHM_FR,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_FR(*)(const_CHM_FR,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_copy_factor"); - return fun(L, Common); -} - -inline int -cholmod_change_factor(int to_xtype, int to_ll, int to_super, int to_packed, - int to_monotonic, CHM_FR L, CHM_CM Common) -{ - static int(*fun)(int,int,int,int,int,CHM_FR,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(int,int,int,int,int,CHM_FR,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_change_factor"); - return fun(to_xtype, to_ll, to_super, to_packed, to_monotonic, L, Common); -} - -inline CHM_SP -cholmod_copy_sparse(const_CHM_SP A, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_copy_sparse"); - return fun(A, Common); -} - -inline CHM_SP -cholmod_factor_to_sparse(const_CHM_FR L, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_FR,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_FR,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_factor_to_sparse"); - return fun(L, Common); -} - -inline CHM_SP -cholmod_submatrix(const_CHM_SP A, int *rset, int rsize, int *cset, - int csize, int values, int sorted, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,int*,int,int*,int, - int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,int*,int,int*, - int,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_submatrix"); - return fun(A, rset, rsize, cset, csize, values, sorted, Common); -} - -inline CHM_SP -cholmod_dense_to_sparse(const_CHM_DN X, int values, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_DN,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_DN,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_dense_to_sparse"); - return fun(X, values, Common); -} - -inline int -cholmod_factorize(const_CHM_SP A, CHM_FR L, CHM_CM Common) -{ - static int(*fun)(const_CHM_SP,CHM_FR,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(const_CHM_SP,CHM_FR,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_factorize"); - return fun(A, L, Common); -} - -inline int -cholmod_factorize_p(const_CHM_SP A, double *beta, int *fset, - size_t fsize, CHM_FR L, - CHM_CM Common) -{ - static int(*fun)(const_CHM_SP,double*,int*,size_t, - CHM_FR,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(const_CHM_SP,double*,int*,size_t, - CHM_FR,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_factorize_p"); - return fun(A, beta, fset, fsize, L, Common); -} - -inline int -cholmod_finish(CHM_CM Common) -{ - - static int(*fun)(CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_CM)) - R_GetCCallable("Matrix", "cholmod_finish"); - return fun(Common); -} - -inline int -cholmod_sort(CHM_SP A, CHM_CM Common) -{ - static int(*fun)(CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_sort"); - return fun(A, Common); -} - -inline int -cholmod_free_dense(CHM_DN *A, CHM_CM Common) -{ - static int(*fun)(CHM_DN*,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_DN*,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_free_dense"); - return fun(A, Common); -} - -inline int -cholmod_free_factor(CHM_FR *L, CHM_CM Common) -{ - static int(*fun)(CHM_FR*,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_FR*,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_free_factor"); - return fun(L, Common); -} - -inline int -cholmod_free_sparse(CHM_SP *A, CHM_CM Common) -{ - static int(*fun)(CHM_SP*,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_SP*,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_free_sparse"); - return fun(A, Common); -} - -inline int -cholmod_free_triplet(cholmod_triplet **T, CHM_CM Common) -{ - static int(*fun)(cholmod_triplet**,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(cholmod_triplet**,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_free_triplet"); - return fun(T, Common); -} - -inline long int -cholmod_nnz(const_CHM_SP A, CHM_CM Common) { - static long(*fun)(const_CHM_SP,CHM_CM) = NULL; - if (fun == NULL) - fun = (long(*)(const_CHM_SP,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_nnz"); - return fun(A, Common); -} - -inline int -cholmod_sdmult(const_CHM_SP A, int transpose, - const double *alpha, const double *beta, - const_CHM_DN X, CHM_DN Y, - CHM_CM Common) -{ - static int(*fun)(const_CHM_SP,int,const double*, - const double*,const_CHM_DN, - CHM_DN,CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(const_CHM_SP,int,const double*, - const double*, const_CHM_DN, - CHM_DN,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_sdmult"); - return fun(A, transpose, alpha, beta, X, Y, Common); -} - -inline CHM_SP -cholmod_ssmult(const_CHM_SP A, const_CHM_SP B, int stype, - int values, int sorted, CHM_CM Common) { - static CHM_SP(*fun)(const_CHM_SP,const_CHM_SP, - int,int,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,const_CHM_SP, - int,int,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_ssmult"); - return fun(A, B, stype, values, sorted, Common); -} - -inline CHM_DN -cholmod_solve(int sys, const_CHM_FR L, const_CHM_DN B, CHM_CM Common) { - static CHM_DN(*fun)(int,const_CHM_FR,const_CHM_DN, - CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_DN(*)(int,const_CHM_FR,const_CHM_DN, - CHM_CM)) - R_GetCCallable("Matrix", "cholmod_solve"); - return fun(sys, L, B, Common); -} - -inline CHM_SP -cholmod_speye(size_t nrow, size_t ncol, - int xtype, CHM_CM Common) -{ - static CHM_SP(*fun)(size_t,size_t,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(size_t,size_t,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_speye"); - return fun(nrow, ncol, xtype, Common); -} - -inline CHM_SP -cholmod_spsolve(int sys, const_CHM_FR L, - const_CHM_SP B, CHM_CM Common) -{ - static CHM_SP(*fun)(int,const_CHM_FR, - const_CHM_SP, CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(int,const_CHM_FR, - const_CHM_SP, CHM_CM)) - R_GetCCallable("Matrix", "cholmod_spsolve"); - return fun(sys, L, B, Common); -} - -inline int -cholmod_defaults (CHM_CM Common) -{ - static int(*fun)(CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_CM)) - R_GetCCallable("Matrix", "cholmod_defaults"); - return fun(Common); -} - -/* extern cholmod_common c; */ - -inline void -R_cholmod_error(int status, const char *file, int line, const char *message) -{ -/* NB: keep in sync with R_cholmod_error(), ../../src/chm_common.c */ - - if(status < 0) { -/* Note: Matrix itself uses CHM_set_common_env, CHM_store_common - * and CHM_restore_common to preserve settings through error calls. - * Consider defining your own error handler, *and* possibly restoring - * *your* version of the cholmod_common that *you* use. - */ - Rf_error("Cholmod error '%s' at file:%s, line %d", message, file, line); - } - else - Rf_warning("Cholmod warning '%s' at file:%s, line %d", - message, file, line); -} - -inline int -cholmod_start(CHM_CM Common) -{ - int val; - static int(*fun)(CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(CHM_CM)) - R_GetCCallable("Matrix", "cholmod_start"); - val = fun(Common); -/*-- NB: keep in sync with R_cholmod_start() --> ../../src/chm_common.c */ - /* do not allow CHOLMOD printing - currently */ - Common->print_function = NULL;/* was R_cholmod_printf; /.* Rprintf gives warning */ -/* Consider using your own error handler: */ - Common->error_handler = R_cholmod_error; - return val; -} - -inline CHM_SP -cholmod_transpose(const_CHM_SP A, int values, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,int,CHM_CM)) - R_GetCCallable("Matrix", "cholmod_transpose"); - return fun(A, values, Common); -} - -inline CHM_SP -cholmod_vertcat(const_CHM_SP A, const_CHM_SP B, int values, CHM_CM Common) -{ - static CHM_SP(*fun)(const_CHM_SP,const_CHM_SP,int,CHM_CM) = NULL; - if (fun == NULL) - fun = (CHM_SP(*)(const_CHM_SP,const_CHM_SP, int, CHM_CM)) - R_GetCCallable("Matrix", "cholmod_vertcat"); - return fun(A, B, values, Common); -} - -inline int -cholmod_scale(const_CHM_DN S, int scale, CHM_SP A, - CHM_CM Common) -{ - static int(*fun)(const_CHM_DN,int,CHM_SP, CHM_CM) = NULL; - if (fun == NULL) - fun = (int(*)(const_CHM_DN,int,CHM_SP, CHM_CM)) - R_GetCCallable("Matrix", "cholmod_scale"); - return fun(S, scale, A, Common); -} - -#ifdef __cplusplus -} -#endif diff --git a/inst/include/RcppEigenWrap.h b/inst/include/RcppEigenWrap.h index 595ca831..c9a05447 100644 --- a/inst/include/RcppEigenWrap.h +++ b/inst/include/RcppEigenWrap.h @@ -24,51 +24,6 @@ namespace Rcpp{ - namespace RcppEigen{ - - template - SEXP Eigen_cholmod_wrap(const Eigen::CholmodDecomposition >& obj) { - typedef T* Tpt; - const cholmod_factor* f = obj.factor(); - if (f->minor < f->n) - throw std::runtime_error("CHOLMOD factorization was unsuccessful"); - - //FIXME: Should extend this selection according to T - S4 ans(std::string(f->is_super ? "dCHMsuper" : "dCHMsimpl")); - IntegerVector dd(2); - dd[0] = dd[1] = f->n; - ans.slot("Dim") = dd; - ans.slot("perm") = ::Rcpp::wrap((int*)f->Perm, (int*)f->Perm + f->n); - ans.slot("colcount") = ::Rcpp::wrap((int*)f->ColCount, (int*)f->ColCount + f->n); - IntegerVector tt(f->is_super ? 6 : 4); - tt[0] = f->ordering; tt[1] = f->is_ll; - tt[2] = f->is_super; tt[3] = f->is_monotonic; - ans.slot("type") = tt; - if (f->is_super) { - tt[4] = f->maxcsize; tt[5] = f->maxesize; - ans.slot("super") = ::Rcpp::wrap((int*)f->super, ((int*)f->super) + f->nsuper + 1); - ans.slot("pi") = ::Rcpp::wrap((int*)f->pi, ((int*)f->pi) + f->nsuper + 1); - ans.slot("px") = ::Rcpp::wrap((int*)f->px, ((int*)f->px) + f->nsuper + 1); - ans.slot("s") = ::Rcpp::wrap((int*)f->s, ((int*)f->s) + f->ssize); - ans.slot("x") = ::Rcpp::wrap((Tpt)f->x, ((T*)f->x) + f->xsize); - } else { - ans.slot("i") = ::Rcpp::wrap((int*)f->i, ((int*)f->i) + f->nzmax); - ans.slot("p") = ::Rcpp::wrap((int*)f->p, ((int*)f->p) + f->n + 1); - ans.slot("x") = ::Rcpp::wrap((Tpt)f->x, ((T*)f->x) + f->nzmax); - ans.slot("nz") = ::Rcpp::wrap((int*)f->nz, ((int*)f->nz) + f->n); - ans.slot("nxt") = ::Rcpp::wrap((int*)f->next, ((int*)f->next) + f->n + 2); - ans.slot("prv") = ::Rcpp::wrap((int*)f->prev, ((int*)f->prev) + f->n + 2); - } - return ::Rcpp::wrap(ans); - } - - } /* namespace RcppEigen */ - - template - SEXP wrap(const Eigen::CholmodDecomposition >& obj) { - return RcppEigen::Eigen_cholmod_wrap(obj); - } - namespace RcppEigen{ // helper trait to identify if T is a plain object type