// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
//
// cholmod.h: selected headers from Tim Davis's CHOLMOD package
//
// Copyright (C) 2011 Douglas Bates, Martin Maechler, Dirk Eddelbuettel and Romain Francois
//
// 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 .
#ifndef RcppEigen_CHOLMOD_H
#define RcppEigen_CHOLMOD_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" {
#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.
*/
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 */