12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

move colamd and spqr_mini out of gtsam

release/4.3a0
Kai Ni 2010-10-14 02:41:08 +00:00
parent fd794e0da2
commit 5cbf67eeba
38 changed files with 50 additions and 15705 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Makefile.am
View File

@ -11,15 +11,14 @@ ACLOCAL_AMFLAGS = -I m4
AUTOMAKE_OPTIONS = foreign nostdinc
# All the sub-directories that need to be built
SUBDIRS = CppUnitLite colamd spqr_mini base geometry inference linear nonlinear slam . tests wrap examples
SUBDIRS = CppUnitLite base geometry inference linear nonlinear slam . tests wrap examples
# And the corresponding libraries produced
SUBLIBS = colamd/libcolamd.la \
base/libbase.la geometry/libgeometry.la inference/libinference.la \
SUBLIBS = base/libbase.la geometry/libgeometry.la inference/libinference.la \
linear/liblinear.la nonlinear/libnonlinear.la slam/libslam.la
if USE_LAPACK
SUBLIBS += spqr_mini/libspqr_mini.la
SUBLIBS += -L$(SparseLib) -lcholmod -lspqr
endif
# TODO: UFconfig, CCOLAMD, and LDL automake magic without adding or touching any file
@ -31,7 +30,7 @@ endif
lib_LTLIBRARIES = libgtsam.la
libgtsam_la_SOURCES =
nodist_EXTRA_libgtsam_la_SOURCES = dummy.cxx
libgtsam_la_LIBADD = $(SUBLIBS)
libgtsam_la_LIBADD = $(SUBLIBS) -L$(SparseLib) -lcolamd -lccolamd
libgtsam_la_LDFLAGS = -no-undefined -version-info 0:0:0
if USE_ACCELERATE_MACOS

4
base/Makefile.am
View File

@ -50,7 +50,7 @@ base_HEADERS = $(headers)
noinst_LTLIBRARIES = libbase.la
libbase_la_SOURCES = $(sources)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
if USE_BLAS
AM_CPPFLAGS += -DGT_USE_CBLAS
@ -78,7 +78,7 @@ AM_LDFLAGS += -lcblas -latlas
endif
if USE_LAPACK
LDADD += ../spqr_mini/libspqr_mini.la
AM_LDFLAGS += -L$(SparseLib) -lcholmod -lspqr
endif
if USE_LAPACK_LINUX

5
base/SPQRUtil.h
View File

@ -11,7 +11,10 @@
#include <gtsam/base/Matrix.h>
#ifdef GT_USE_LAPACK
#include <gtsam/spqr_mini/spqr_subset.hpp>
extern "C" {
#include <cholmod.h>
}
#include <spqr.hpp>
namespace gtsam {

16
colamd/Makefile.am
View File

@ -1,16 +0,0 @@
#----------------------------------------------------------------------------------------------------
# colamd
# replaced Makefile with automake for easy linking
#----------------------------------------------------------------------------------------------------
# Create a libtool library that is not installed
# It will be packaged in the toplevel libgtsam.la as specfied in ../Makefile.am
noinst_LTLIBRARIES = libcolamd.la
# We normally would not install these headers
# but they are included in the templated class FactorGraph-inl.h so we need them
colamddir = $(pkgincludedir)/colamd
colamd_HEADERS = colamd.h ccolamd.h UFconfig.h
# These are the sources for the library:
libcolamd_la_SOURCES = colamd.c colamd_global.c ccolamd.c ccolamd_global.c

118
colamd/UFconfig.h
View File

@ -1,118 +0,0 @@
/* ========================================================================== */
/* === UFconfig.h =========================================================== */
/* ========================================================================== */
/* Configuration file for SuiteSparse: a Suite of Sparse matrix packages
* (AMD, COLAMD, CCOLAMD, CAMD, CHOLMOD, UMFPACK, CXSparse, and others).
*
* UFconfig.h provides the definition of the long integer. On most systems,
* a C program can be compiled in LP64 mode, in which long's and pointers are
* both 64-bits, and int's are 32-bits. Windows 64, however, uses the LLP64
* model, in which int's and long's are 32-bits, and long long's and pointers
* are 64-bits.
*
* SuiteSparse packages that include long integer versions are
* intended for the LP64 mode. However, as a workaround for Windows 64
* (and perhaps other systems), the long integer can be redefined.
*
* If _WIN64 is defined, then the __int64 type is used instead of long.
*
* The long integer can also be defined at compile time. For example, this
* could be added to UFconfig.mk:
*
* CFLAGS = -O -D'UF_long=long long' -D'UF_long_max=9223372036854775801' \
* -D'UF_long_id="%lld"'
*
* This file defines UF_long as either long (on all but _WIN64) or
* __int64 on Windows 64. The intent is that a UF_long is always a 64-bit
* integer in a 64-bit code. ptrdiff_t might be a better choice than long;
* it is always the same size as a pointer.
*
* This file also defines the SUITESPARSE_VERSION and related definitions.
*
* Copyright (c) 2007, University of Florida. No licensing restrictions
* apply to this file or to the UFconfig directory. Author: Timothy A. Davis.
*/
#ifndef _UFCONFIG_H
#define _UFCONFIG_H
#ifdef __cplusplus
extern "C" {
#endif
#include <limits.h>
/* ========================================================================== */
/* === UF_long ============================================================== */
/* ========================================================================== */
#ifndef UF_long
#ifdef _WIN64
#define UF_long __int64
#define UF_long_max _I64_MAX
#define UF_long_id "%I64d"
#else
#define UF_long long
#define UF_long_max LONG_MAX
#define UF_long_id "%ld"
#endif
#endif
/* ========================================================================== */
/* === SuiteSparse version ================================================== */
/* ========================================================================== */
/* SuiteSparse is not a package itself, but a collection of packages, some of
* which must be used together (UMFPACK requires AMD, CHOLMOD requires AMD,
* COLAMD, CAMD, and CCOLAMD, etc). A version number is provided here for the
* collection itself. The versions of packages within each version of
* SuiteSparse are meant to work together. Combining one packge from one
* version of SuiteSparse, with another package from another version of
* SuiteSparse, may or may not work.
*
* SuiteSparse Version 3.2.0 contains the following packages:
*
* AMD version 2.2.0
* CAMD version 2.2.0
* COLAMD version 2.7.1
* CCOLAMD version 2.7.1
* CHOLMOD version 1.7.0
* CSparse version 2.2.1
* CXSparse version 2.2.1
* KLU version 1.0.1
* BTF version 1.0.1
* LDL version 2.0.1
* UFconfig version number is the same as SuiteSparse
* UMFPACK version 5.2.0
* RBio version 1.1.1
* UFcollection version 1.1.1
* LINFACTOR version 1.1.0
* MESHND version 1.1.0
* SSMULT version 1.1.0
* MATLAB_Tools no specific version number
* SuiteSparseQR version 1.1.0
*
* Other package dependencies:
* BLAS required by CHOLMOD and UMFPACK
* LAPACK required by CHOLMOD
* METIS 4.0.1 required by CHOLMOD (optional) and KLU (optional)
*/
#define SUITESPARSE_DATE "Sept 20, 2008"
#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
#define SUITESPARSE_MAIN_VERSION 3
#define SUITESPARSE_SUB_VERSION 2
#define SUITESPARSE_SUBSUB_VERSION 0
#define SUITESPARSE_VERSION \
SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

4620
colamd/ccolamd.c
View File

File diff suppressed because it is too large Load Diff

365
colamd/ccolamd.h
View File

@ -1,365 +0,0 @@
/* ========================================================================== */
/* === CCOLAMD/ccolamd.h ==================================================== */
/* ========================================================================== */
/* ----------------------------------------------------------------------------
* CCOLAMD Copyright (C), Univ. of Florida. Authors: Timothy A. Davis,
* Sivasankaran Rajamanickam, and Stefan Larimore
* See License.txt for the Version 2.1 of the GNU Lesser General Public License
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/*
* You must include this file (ccolamd.h) in any routine that uses ccolamd,
* csymamd, or the related macros and definitions.
*/
#ifndef CCOLAMD_H
#define CCOLAMD_H
/* make it easy for C++ programs to include CCOLAMD */
#ifdef __cplusplus
extern "C" {
#endif
/* for size_t definition: */
#include <stdlib.h>
/* ========================================================================== */
/* === CCOLAMD version ====================================================== */
/* ========================================================================== */
/* All versions of CCOLAMD will include the following definitions.
* As an example, to test if the version you are using is 1.3 or later:
*
* if (CCOLAMD_VERSION >= CCOLAMD_VERSION_CODE (1,3)) ...
*
* This also works during compile-time:
*
* #if CCOLAMD_VERSION >= CCOLAMD_VERSION_CODE (1,3)
* printf ("This is version 1.3 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*/
#define CCOLAMD_DATE "Nov 30, 2009"
#define CCOLAMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define CCOLAMD_MAIN_VERSION 2
#define CCOLAMD_SUB_VERSION 7
#define CCOLAMD_SUBSUB_VERSION 2
#define CCOLAMD_VERSION \
CCOLAMD_VERSION_CODE(CCOLAMD_MAIN_VERSION,CCOLAMD_SUB_VERSION)
/* ========================================================================== */
/* === Knob and statistics definitions ====================================== */
/* ========================================================================== */
/* size of the knobs [ ] array. Only knobs [0..3] are currently used. */
#define CCOLAMD_KNOBS 20
/* number of output statistics. Only stats [0..10] are currently used. */
#define CCOLAMD_STATS 20
/* knobs [0] and stats [0]: dense row knob and output statistic. */
#define CCOLAMD_DENSE_ROW 0
/* knobs [1] and stats [1]: dense column knob and output statistic. */
#define CCOLAMD_DENSE_COL 1
/* knobs [2]: aggressive absorption option */
#define CCOLAMD_AGGRESSIVE 2
/* knobs [3]: LU or Cholesky factorization option */
#define CCOLAMD_LU 3
/* stats [2]: memory defragmentation count output statistic */
#define CCOLAMD_DEFRAG_COUNT 2
/* stats [3]: ccolamd status: zero OK, > 0 warning or notice, < 0 error */
#define CCOLAMD_STATUS 3
/* stats [4..6]: error info, or info on jumbled columns */
#define CCOLAMD_INFO1 4
#define CCOLAMD_INFO2 5
#define CCOLAMD_INFO3 6
/* stats [7]: number of originally empty rows */
#define CCOLAMD_EMPTY_ROW 7
/* stats [8]: number of originally empty cols */
#define CCOLAMD_EMPTY_COL 8
/* stats [9]: number of rows with entries only in dense cols */
#define CCOLAMD_NEWLY_EMPTY_ROW 9
/* stats [10]: number of cols with entries only in dense rows */
#define CCOLAMD_NEWLY_EMPTY_COL 10
/* error codes returned in stats [3]: */
#define CCOLAMD_OK (0)
#define CCOLAMD_OK_BUT_JUMBLED (1)
#define CCOLAMD_ERROR_A_not_present (-1)
#define CCOLAMD_ERROR_p_not_present (-2)
#define CCOLAMD_ERROR_nrow_negative (-3)
#define CCOLAMD_ERROR_ncol_negative (-4)
#define CCOLAMD_ERROR_nnz_negative (-5)
#define CCOLAMD_ERROR_p0_nonzero (-6)
#define CCOLAMD_ERROR_A_too_small (-7)
#define CCOLAMD_ERROR_col_length_negative (-8)
#define CCOLAMD_ERROR_row_index_out_of_bounds (-9)
#define CCOLAMD_ERROR_out_of_memory (-10)
#define CCOLAMD_ERROR_invalid_cmember (-11)
#define CCOLAMD_ERROR_internal_error (-999)
/* ========================================================================== */
/* === Prototypes of user-callable routines ================================= */
/* ========================================================================== */
/* define UF_long */
#include "UFconfig.h"
size_t ccolamd_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
int nnz, /* nonzeros in A */
int n_row, /* number of rows in A */
int n_col /* number of columns in A */
) ;
size_t ccolamd_l_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
UF_long nnz, /* nonzeros in A */
UF_long n_row, /* number of rows in A */
UF_long n_col /* number of columns in A */
) ;
void ccolamd_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [CCOLAMD_KNOBS] /* parameter settings for ccolamd */
) ;
void ccolamd_l_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [CCOLAMD_KNOBS] /* parameter settings for ccolamd */
) ;
int ccolamd /* returns (1) if successful, (0) otherwise*/
( /* A and p arguments are modified on output */
int n_row, /* number of rows in A */
int n_col, /* number of columns in A */
int Alen, /* size of the array A */
int A [ ], /* row indices of A, of size Alen */
int p [ ], /* column pointers of A, of size n_col+1 */
double knobs [CCOLAMD_KNOBS],/* parameter settings for ccolamd */
int stats [CCOLAMD_STATS], /* ccolamd output statistics and error codes */
int cmember [ ] /* Constraint set of A, of size n_col */
) ;
UF_long ccolamd_l /* same as ccolamd, but with UF_long integers */
(
UF_long n_row,
UF_long n_col,
UF_long Alen,
UF_long A [ ],
UF_long p [ ],
double knobs [CCOLAMD_KNOBS],
UF_long stats [CCOLAMD_STATS],
UF_long cmember [ ]
) ;
int csymamd /* return (1) if OK, (0) otherwise */
(
int n, /* number of rows and columns of A */
int A [ ], /* row indices of A */
int p [ ], /* column pointers of A */
int perm [ ], /* output permutation, size n_col+1 */
double knobs [CCOLAMD_KNOBS],/* parameters (uses defaults if NULL) */
int stats [CCOLAMD_STATS], /* output statistics and error codes */
void * (*allocate) (size_t, size_t), /* pointer to calloc (ANSI C) or */
/* mxCalloc (for MATLAB mexFunction) */
void (*release) (void *), /* pointer to free (ANSI C) or */
/* mxFree (for MATLAB mexFunction) */
int cmember [ ], /* Constraint set of A */
int stype /* 0: use both parts, >0: upper, <0: lower */
) ;
UF_long csymamd_l /* same as csymamd, but with UF_long integers */
(
UF_long n,
UF_long A [ ],
UF_long p [ ],
UF_long perm [ ],
double knobs [CCOLAMD_KNOBS],
UF_long stats [CCOLAMD_STATS],
void * (*allocate) (size_t, size_t),
void (*release) (void *),
UF_long cmember [ ],
UF_long stype
) ;
void ccolamd_report
(
int stats [CCOLAMD_STATS]
) ;
void ccolamd_l_report
(
UF_long stats [CCOLAMD_STATS]
) ;
void csymamd_report
(
int stats [CCOLAMD_STATS]
) ;
void csymamd_l_report
(
UF_long stats [CCOLAMD_STATS]
) ;
/* ========================================================================== */
/* === Prototypes of "expert" routines ====================================== */
/* ========================================================================== */
/* These routines are meant to be used internally, or in a future version of
* UMFPACK. They appear here so that UMFPACK can use them, but they should not
* be called directly by the user.
*/
int ccolamd2
( /* A and p arguments are modified on output */
int n_row, /* number of rows in A */
int n_col, /* number of columns in A */
int Alen, /* size of the array A */
int A [ ], /* row indices of A, of size Alen */
int p [ ], /* column pointers of A, of size n_col+1 */
double knobs [CCOLAMD_KNOBS],/* parameter settings for ccolamd */
int stats [CCOLAMD_STATS], /* ccolamd output statistics and error codes */
/* each Front_ array is of size n_col+1: */
int Front_npivcol [ ], /* # pivot cols in each front */
int Front_nrows [ ], /* # of rows in each front (incl. pivot rows) */
int Front_ncols [ ], /* # of cols in each front (incl. pivot cols) */
int Front_parent [ ], /* parent of each front */
int Front_cols [ ], /* link list of pivot columns for each front */
int *p_nfr, /* total number of frontal matrices */
int InFront [ ], /* InFront [row] = f if row in front f */
int cmember [ ] /* Constraint set of A */
) ;
UF_long ccolamd2_l /* same as ccolamd2, but with UF_long integers */
(
UF_long n_row,
UF_long n_col,
UF_long Alen,
UF_long A [ ],
UF_long p [ ],
double knobs [CCOLAMD_KNOBS],
UF_long stats [CCOLAMD_STATS],
UF_long Front_npivcol [ ],
UF_long Front_nrows [ ],
UF_long Front_ncols [ ],
UF_long Front_parent [ ],
UF_long Front_cols [ ],
UF_long *p_nfr,
UF_long InFront [ ],
UF_long cmember [ ]
) ;
void ccolamd_apply_order
(
int Front [ ],
const int Order [ ],
int Temp [ ],
int nn,
int nfr
) ;
void ccolamd_l_apply_order
(
UF_long Front [ ],
const UF_long Order [ ],
UF_long Temp [ ],
UF_long nn,
UF_long nfr
) ;
void ccolamd_fsize
(
int nn,
int MaxFsize [ ],
int Fnrows [ ],
int Fncols [ ],
int Parent [ ],
int Npiv [ ]
) ;
void ccolamd_l_fsize
(
UF_long nn,
UF_long MaxFsize [ ],
UF_long Fnrows [ ],
UF_long Fncols [ ],
UF_long Parent [ ],
UF_long Npiv [ ]
) ;
void ccolamd_postorder
(
int nn,
int Parent [ ],
int Npiv [ ],
int Fsize [ ],
int Order [ ],
int Child [ ],
int Sibling [ ],
int Stack [ ],
int Front_cols [ ],
int cmember [ ]
) ;
void ccolamd_l_postorder
(
UF_long nn,
UF_long Parent [ ],
UF_long Npiv [ ],
UF_long Fsize [ ],
UF_long Order [ ],
UF_long Child [ ],
UF_long Sibling [ ],
UF_long Stack [ ],
UF_long Front_cols [ ],
UF_long cmember [ ]
) ;
int ccolamd_post_tree
(
int root,
int k,
int Child [ ],
const int Sibling [ ],
int Order [ ],
int Stack [ ]
) ;
UF_long ccolamd_l_post_tree
(
UF_long root,
UF_long k,
UF_long Child [ ],
const UF_long Sibling [ ],
UF_long Order [ ],
UF_long Stack [ ]
) ;
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN int (*ccolamd_printf) (const char *, ...) ;
#ifdef __cplusplus
}
#endif
#endif

25
colamd/ccolamd_global.c
View File

@ -1,25 +0,0 @@
/* ========================================================================== */
/* === ccolamd_global.c ===================================================== */
/* ========================================================================== */
/* ----------------------------------------------------------------------------
* CCOLAMD Copyright (C), Univ. of Florida. Authors: Timothy A. Davis,
* Sivasankaran Rajamanickam, and Stefan Larimore
* See License.txt for the Version 2.1 of the GNU Lesser General Public License
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Global variables for CCOLAMD */
#ifndef NPRINT
#ifdef MATLAB_MEX_FILE
#include "mex.h"
int (*ccolamd_printf) (const char *, ...) = mexPrintf ;
#else
#include <stdio.h>
int (*ccolamd_printf) (const char *, ...) = printf ;
#endif
#else
int (*ccolamd_printf) (const char *, ...) = ((void *) 0) ;
#endif

3611
colamd/colamd.c
View File

File diff suppressed because it is too large Load Diff

255
colamd/colamd.h
View File

@ -1,255 +0,0 @@
/* ========================================================================== */
/* === colamd/symamd prototypes and definitions ============================= */
/* ========================================================================== */
/* COLAMD / SYMAMD include file
You must include this file (colamd.h) in any routine that uses colamd,
symamd, or the related macros and definitions.
Authors:
The authors of the code itself are Stefan I. Larimore and Timothy A.
Davis (davis at cise.ufl.edu), University of Florida. The algorithm was
developed in collaboration with John Gilbert, Xerox PARC, and Esmond
Ng, Oak Ridge National Laboratory.
Acknowledgements:
This work was supported by the National Science Foundation, under
grants DMS-9504974 and DMS-9803599.
Notice:
Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use, copy, modify, and/or distribute
this program, provided that the Copyright, this License, and the
Availability of the original version is retained on all copies and made
accessible to the end-user of any code or package that includes COLAMD
or any modified version of COLAMD.
Availability:
The colamd/symamd library is available at
http://www.cise.ufl.edu/research/sparse/colamd/
This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h
file. It is required by the colamd.c, colamdmex.c, and symamdmex.c
files, and by any C code that calls the routines whose prototypes are
listed below, or that uses the colamd/symamd definitions listed below.
*/
#ifndef COLAMD_H
#define COLAMD_H
/* make it easy for C++ programs to include COLAMD */
#ifdef __cplusplus
extern "C" {
#endif
/* ========================================================================== */
/* === Include files ======================================================== */
/* ========================================================================== */
#include <stdlib.h>
/* ========================================================================== */
/* === COLAMD version ======================================================= */
/* ========================================================================== */
/* COLAMD Version 2.4 and later will include the following definitions.
* As an example, to test if the version you are using is 2.4 or later:
*
* #ifdef COLAMD_VERSION
* if (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4)) ...
* #endif
*
* This also works during compile-time:
*
* #if defined(COLAMD_VERSION) && (COLAMD_VERSION >= COLAMD_VERSION_CODE (2,4))
* printf ("This is version 2.4 or later\n") ;
* #else
* printf ("This is an early version\n") ;
* #endif
*
* Versions 2.3 and earlier of COLAMD do not include a #define'd version number.
*/
#define COLAMD_DATE "Nov 1, 2007"
#define COLAMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define COLAMD_MAIN_VERSION 2
#define COLAMD_SUB_VERSION 7
#define COLAMD_SUBSUB_VERSION 1
#define COLAMD_VERSION \
COLAMD_VERSION_CODE(COLAMD_MAIN_VERSION,COLAMD_SUB_VERSION)
/* ========================================================================== */
/* === Knob and statistics definitions ====================================== */
/* ========================================================================== */
/* size of the knobs [ ] array. Only knobs [0..1] are currently used. */
#define COLAMD_KNOBS 20
/* number of output statistics. Only stats [0..6] are currently used. */
#define COLAMD_STATS 20
/* knobs [0] and stats [0]: dense row knob and output statistic. */
#define COLAMD_DENSE_ROW 0
/* knobs [1] and stats [1]: dense column knob and output statistic. */
#define COLAMD_DENSE_COL 1
/* knobs [2]: aggressive absorption */
#define COLAMD_AGGRESSIVE 2
/* stats [2]: memory defragmentation count output statistic */
#define COLAMD_DEFRAG_COUNT 2
/* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */
#define COLAMD_STATUS 3
/* stats [4..6]: error info, or info on jumbled columns */
#define COLAMD_INFO1 4
#define COLAMD_INFO2 5
#define COLAMD_INFO3 6
/* error codes returned in stats [3]: */
#define COLAMD_OK (0)
#define COLAMD_OK_BUT_JUMBLED (1)
#define COLAMD_ERROR_A_not_present (-1)
#define COLAMD_ERROR_p_not_present (-2)
#define COLAMD_ERROR_nrow_negative (-3)
#define COLAMD_ERROR_ncol_negative (-4)
#define COLAMD_ERROR_nnz_negative (-5)
#define COLAMD_ERROR_p0_nonzero (-6)
#define COLAMD_ERROR_A_too_small (-7)
#define COLAMD_ERROR_col_length_negative (-8)
#define COLAMD_ERROR_row_index_out_of_bounds (-9)
#define COLAMD_ERROR_out_of_memory (-10)
#define COLAMD_ERROR_internal_error (-999)
/* ========================================================================== */
/* === Prototypes of user-callable routines ================================= */
/* ========================================================================== */
/* define UF_long */
#include "UFconfig.h"
size_t colamd_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
int nnz, /* nonzeros in A */
int n_row, /* number of rows in A */
int n_col /* number of columns in A */
) ;
size_t colamd_l_recommended /* returns recommended value of Alen, */
/* or 0 if input arguments are erroneous */
(
UF_long nnz, /* nonzeros in A */
UF_long n_row, /* number of rows in A */
UF_long n_col /* number of columns in A */
) ;
void colamd_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [COLAMD_KNOBS] /* parameter settings for colamd */
) ;
void colamd_l_set_defaults /* sets default parameters */
( /* knobs argument is modified on output */
double knobs [COLAMD_KNOBS] /* parameter settings for colamd */
) ;
int colamd /* returns (1) if successful, (0) otherwise*/
( /* A and p arguments are modified on output */
int n_row, /* number of rows in A */
int n_col, /* number of columns in A */
int Alen, /* size of the array A */
int A [], /* row indices of A, of size Alen */
int p [], /* column pointers of A, of size n_col+1 */
double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
int stats [COLAMD_STATS] /* colamd output statistics and error codes */
) ;
UF_long colamd_l /* returns (1) if successful, (0) otherwise*/
( /* A and p arguments are modified on output */
UF_long n_row, /* number of rows in A */
UF_long n_col, /* number of columns in A */
UF_long Alen, /* size of the array A */
UF_long A [], /* row indices of A, of size Alen */
UF_long p [], /* column pointers of A, of size n_col+1 */
double knobs [COLAMD_KNOBS],/* parameter settings for colamd */
UF_long stats [COLAMD_STATS]/* colamd output statistics and error codes */
) ;
int symamd /* return (1) if OK, (0) otherwise */
(
int n, /* number of rows and columns of A */
int A [], /* row indices of A */
int p [], /* column pointers of A */
int perm [], /* output permutation, size n_col+1 */
double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
int stats [COLAMD_STATS], /* output statistics and error codes */
void * (*allocate) (size_t, size_t),
/* pointer to calloc (ANSI C) or */
/* mxCalloc (for MATLAB mexFunction) */
void (*release) (void *)
/* pointer to free (ANSI C) or */
/* mxFree (for MATLAB mexFunction) */
) ;
UF_long symamd_l /* return (1) if OK, (0) otherwise */
(
UF_long n, /* number of rows and columns of A */
UF_long A [], /* row indices of A */
UF_long p [], /* column pointers of A */
UF_long perm [], /* output permutation, size n_col+1 */
double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
UF_long stats [COLAMD_STATS], /* output statistics and error codes */
void * (*allocate) (size_t, size_t),
/* pointer to calloc (ANSI C) or */
/* mxCalloc (for MATLAB mexFunction) */
void (*release) (void *)
/* pointer to free (ANSI C) or */
/* mxFree (for MATLAB mexFunction) */
) ;
void colamd_report
(
int stats [COLAMD_STATS]
) ;
void colamd_l_report
(
UF_long stats [COLAMD_STATS]
) ;
void symamd_report
(
int stats [COLAMD_STATS]
) ;
void symamd_l_report
(
UF_long stats [COLAMD_STATS]
) ;
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN int (*colamd_printf) (const char *, ...) ;
#ifdef __cplusplus
}
#endif
#endif /* COLAMD_H */

24
colamd/colamd_global.c
View File

@ -1,24 +0,0 @@
/* ========================================================================== */
/* === colamd_global.c ====================================================== */
/* ========================================================================== */
/* ----------------------------------------------------------------------------
* COLAMD, Copyright (C) 2007, Timothy A. Davis.
* See License.txt for the Version 2.1 of the GNU Lesser General Public License
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Global variables for COLAMD */
#ifndef NPRINT
#ifdef MATLAB_MEX_FILE
#include "mex.h"
int (*colamd_printf) (const char *, ...) = mexPrintf ;
#else
#include <stdio.h>
int (*colamd_printf) (const char *, ...) = printf ;
#endif
#else
int (*colamd_printf) (const char *, ...) = ((void *) 0) ;
#endif

38
colamd/main.cpp
View File

@ -1,38 +0,0 @@
#include "iostream"
#include "colamd.h"
using namespace std;
#define ALEN 100
void use_colamd()
{
int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
int p [ ] = {0, 3, 5, 9, 11} ;
int stats [COLAMD_STATS] ;
colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
for(int i = 0; i < 5; i++)
printf("%d ", p[i]);
printf("\n");
for(int i = 0; i < COLAMD_STATS; i++)
printf("%d ", stats[i]);
printf("\n");
}
int main()
{
/*
A:
[0 x x 0 x
x x 0 0 x
x 0 0 x x
0 0 x 0 0
x x x 0 x
]
*/
//int A [ALEN] = {0, 3, 2, 3, 1, 2, 0, 1, 3, 4, 3} ;
//int p [ ] = {0, 2, 4, 6, 10, 11} ;
use_colamd();
return 0;
}

24
configure.ac
View File

@ -4,12 +4,10 @@
AC_PREREQ(2.59)
AC_INIT(gtsam, 0.0.0, dellaert@cc.gatech.edu)
AM_INIT_AUTOMAKE(gtsam, 0.0.0)
AC_OUTPUT(Makefile CppUnitLite/Makefile colamd/Makefile spqr_mini/Makefile base/Makefile inference/Makefile linear/Makefile geometry/Makefile nonlinear/Makefile slam/Makefile tests/Makefile wrap/Makefile examples/Makefile)
AC_OUTPUT(Makefile CppUnitLite/Makefile base/Makefile inference/Makefile linear/Makefile geometry/Makefile nonlinear/Makefile slam/Makefile tests/Makefile wrap/Makefile examples/Makefile)
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADER([config.h])
AC_CONFIG_SRCDIR([CppUnitLite/Test.cpp])
AC_CONFIG_SRCDIR([colamd/colamd.c])
AC_CONFIG_SRCDIR([spqr_mini/spqr_front.cpp])
AC_CONFIG_SRCDIR([base/DSFVector.cpp])
AC_CONFIG_SRCDIR([geometry/Cal3_S2.cpp])
AC_CONFIG_SRCDIR([inference/SymbolicFactorGraph.cpp])
@ -152,4 +150,24 @@ AC_ARG_WITH([boost],
])
AC_SUBST([boost])
# ask for sparse library include directory
AC_ARG_WITH([sparse-inc],
[AS_HELP_STRING([--with-sparse-inc],
[specify the sparse library include directory (mandatory)])],
[SparseInc=$withval],
[AC_MSG_FAILURE(
[--with-sparse-inc has to be specified])
])
AC_SUBST([SparseInc])
# ask for sparse library lib directory
AC_ARG_WITH([sparse-lib],
[AS_HELP_STRING([--with-sparse-lib],
[specify the sparse library lib directory (mandatory)])],
[SparseLib=$withval],
[AC_MSG_FAILURE(
[--with-sparse-lib has to be specified])
])
AC_SUBST([SparseLib])
AC_OUTPUT

2
examples/Makefile.am
View File

@ -20,7 +20,7 @@ noinst_PROGRAMS += PlanarSLAMSelfContained # Solves SLAM example from tutorial
# rules to build local programs
#----------------------------------------------------------------------------------------------------
AM_LDFLAGS = $(BOOST_LDFLAGS)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
LDADD = ../libgtsam.la
AM_DEFAULT_SOURCE_EXT = .cpp

2
inference/FactorGraph-inl.h
View File

@ -20,7 +20,7 @@
#include <boost/tuple/tuple.hpp>
#include <boost/format.hpp>
#include <boost/graph/prim_minimum_spanning_tree.hpp>
#include <gtsam/colamd/ccolamd.h>
#include <ccolamd.h>
#include <gtsam/inference/FactorGraph.h>
#include <gtsam/inference/graph-inl.h>
#include <gtsam/base/DSF.h>

7
inference/Makefile.am
View File

@ -63,16 +63,16 @@ inferencedir = $(pkgincludedir)/inference
inference_HEADERS = $(headers)
noinst_LTLIBRARIES = libinference.la
libinference_la_SOURCES = $(sources)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
AM_CXXFLAGS =
#----------------------------------------------------------------------------------------------------
# rules to build local programs
#----------------------------------------------------------------------------------------------------
TESTS = $(check_PROGRAMS)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -L$(SparseLib) -lcolamd -lccolamd
LDADD = libinference.la ../base/libbase.la
LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la
LDADD += ../CppUnitLite/libCppUnitLite.a
AM_DEFAULT_SOURCE_EXT = .cpp
if USE_ACCELERATE_MACOS
@ -89,6 +89,5 @@ endif
if USE_LAPACK
AM_CXXFLAGS += -DGT_USE_LAPACK
LDADD += ../spqr_mini/libspqr_mini.la
endif

2
inference/SymbolMap.h
View File

@ -11,7 +11,7 @@
#define GTSAM_SYMBOL_BINARY
#define GTSAM_SYMBOL_SPECIAL
#include <gtsam/inference/Key.h>
#include <gtsam/nonlinear/Key.h>
#include <map>
#include <boost/unordered_map.hpp>

3
inference/inference-inl.h
View File

@ -10,7 +10,8 @@
#include <gtsam/inference/inference.h>
#include <gtsam/inference/FactorGraph-inl.h>
#include <gtsam/inference/BayesNet-inl.h>
#include <gtsam/colamd/ccolamd.h>
#include <ccolamd.h>
#include <boost/foreach.hpp>
#include <boost/format.hpp>

7
linear/Makefile.am
View File

@ -49,16 +49,16 @@ lineardir = $(pkgincludedir)/linear
linear_HEADERS = $(headers)
noinst_LTLIBRARIES = liblinear.la
liblinear_la_SOURCES = $(sources)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
AM_CXXFLAGS =
#----------------------------------------------------------------------------------------------------
# rules to build local programs
#----------------------------------------------------------------------------------------------------
TESTS = $(check_PROGRAMS)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -L$(SparseLib) -lcolamd -lccolamd
LDADD = liblinear.la ../inference/libinference.la ../base/libbase.la
LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la
LDADD += ../CppUnitLite/libCppUnitLite.a
AM_DEFAULT_SOURCE_EXT = .cpp
if USE_ACCELERATE_MACOS
@ -75,6 +75,5 @@ endif
if USE_LAPACK
AM_CXXFLAGS += -DGT_USE_LAPACK
LDADD += ../spqr_mini/libspqr_mini.la
endif

7
nonlinear/Makefile.am
View File

@ -39,16 +39,16 @@ nonlineardir = $(pkgincludedir)/nonlinear
nonlinear_HEADERS = $(headers)
noinst_LTLIBRARIES = libnonlinear.la
libnonlinear_la_SOURCES = $(sources)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
AM_CXXFLAGS =
#----------------------------------------------------------------------------------------------------
# rules to build local programs
#----------------------------------------------------------------------------------------------------
TESTS = $(check_PROGRAMS)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -L$(SparseLib) -lcolamd -lccolamd
LDADD = libnonlinear.la ../linear/liblinear.la ../inference/libinference.la ../base/libbase.la
LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la
LDADD += ../CppUnitLite/libCppUnitLite.a
AM_DEFAULT_SOURCE_EXT = .cpp
if USE_ACCELERATE_MACOS
@ -65,6 +65,5 @@ endif
if USE_LAPACK
AM_CXXFLAGS += -DGT_USE_LAPACK
LDADD += ../spqr_mini/libspqr_mini.la
endif

10
slam/Makefile.am
View File

@ -65,25 +65,21 @@ slamdir = $(pkgincludedir)/slam
slam_HEADERS = $(headers)
noinst_LTLIBRARIES = libslam.la
libslam_la_SOURCES = $(sources)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
#----------------------------------------------------------------------------------------------------
# rules to build local programs
#----------------------------------------------------------------------------------------------------
TESTS = $(check_PROGRAMS)
AM_DEFAULT_SOURCE_EXT = .cpp
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization)
AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -L$(SparseLib) -lcolamd -lccolamd
LDADD = libslam.la ../geometry/libgeometry.la ../nonlinear/libnonlinear.la ../linear/liblinear.la ../inference/libinference.la ../base/libbase.la
LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la
LDADD += ../CppUnitLite/libCppUnitLite.a
if USE_ACCELERATE_MACOS
AM_LDFLAGS += -Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate
endif
if USE_LAPACK
LDADD += ../spqr_mini/libspqr_mini.la
endif
# rule to run an executable
%.run: % $(LDADD)
./$^

25
spqr_mini/Makefile.am
View File

@ -1,25 +0,0 @@
# Only install if LAPACK enabled
if USE_LAPACK
sources = cholmod_error.c cholmod_common.c cholmod_memory.c spqr_front.cpp spqr_larftb.cpp
headers = UFconfig.h cholmod_common.h cholmod_internal.h cholmod_blas.h cholmod_core.h
headers += SuiteSparseQR_definitions.h SuiteSparseQR_subset.hpp spqr_subset.hpp spqr_larftb.h spqr_front.h
# Create a libtool library that is not installed
# It will be packaged in the toplevel libgtsam.la as specfied in ../Makefile.am
# The headers are not installed either
noinst_LTLIBRARIES = libspqr_mini.la
libspqr_mini_la_SOURCES = $(sources)
noinst_HEADERS = $(headers)
AM_CPPFLAGS = -DDLONG # Compiles cholmod in double/long mode
# On Mac, we compile using the BLAS/LAPACK headers in the Accelerate framework
if USE_ACCELERATE_MACOS
AM_CPPFLAGS += -I/System/Library/Frameworks/Accelerate.framework/Headers
endif
endif

69
spqr_mini/SuiteSparseQR_definitions.h
View File

@ -1,69 +0,0 @@
/* ========================================================================== */
/* === SuiteSparseQR_definitions.h ========================================== */
/* ========================================================================== */
/* Core definitions for both C and C++ programs. */
#ifndef SUITESPARSEQR_DEFINITIONS_H
#define SUITESPARSEQR_DEFINITIONS_H
/* ordering options */
#define SPQR_ORDERING_FIXED 0
#define SPQR_ORDERING_NATURAL 1
#define SPQR_ORDERING_COLAMD 2
#define SPQR_ORDERING_GIVEN 3 /* only used for C/C++ interface */
#define SPQR_ORDERING_CHOLMOD 4 /* CHOLMOD best-effort (COLAMD, METIS,...)*/
#define SPQR_ORDERING_AMD 5 /* AMD(A'*A) */
#define SPQR_ORDERING_METIS 6 /* metis(A'*A) */
#define SPQR_ORDERING_DEFAULT 7 /* SuiteSparseQR default ordering */
#define SPQR_ORDERING_BEST 8 /* try COLAMD, AMD, and METIS; pick best */
#define SPQR_ORDERING_BESTAMD 9 /* try COLAMD and AMD; pick best */
/* Let [m n] = size of the matrix after pruning singletons. The default
* ordering strategy is to use COLAMD if m <= 2*n. Otherwise, AMD(A'A) is
* tried. If there is a high fill-in with AMD then try METIS(A'A) and take
* the best of AMD and METIS. METIS is not tried if it isn't installed. */
/* tol options */
#define SPQR_DEFAULT_TOL (-2) /* if tol <= -2, the default tol is used */
#define SPQR_NO_TOL (-1) /* if -2 < tol < 0, then no tol is used */
/* for qmult, method can be 0,1,2,3: */
#define SPQR_QTX 0
#define SPQR_QX 1
#define SPQR_XQT 2
#define SPQR_XQ 3
/* system can be 0,1,2,3: Given Q*R=A*E from SuiteSparseQR_factorize: */
#define SPQR_RX_EQUALS_B 0 /* solve R*X=B or X = R\B */
#define SPQR_RETX_EQUALS_B 1 /* solve R*E'*X=B or X = E*(R\B) */
#define SPQR_RTX_EQUALS_B 2 /* solve R'*X=B or X = R'\B */
#define SPQR_RTX_EQUALS_ETB 3 /* solve R'*X=E'*B or X = R'\(E'*B) */
/* ========================================================================== */
/* === SuiteSparseQR version ================================================ */
/* ========================================================================== */
/*
All versions of SuiteSparseQR will include the following definitions.
As an example, to test if the version you are using is 1.3 or later:
if (SPQR_VERSION >= SPQR_VER_CODE (1,3)) ...
This also works during compile-time:
#if SPQR_VERSION >= SPQR_VER_CODE (1,3)
printf ("This is version 1.3 or later\n") ;
#else
printf ("This is version is earlier than 1.3\n") ;
#endif
*/
#define SPQR_DATE "Nov 30, 2009"
#define SPQR_VER_CODE(main,sub) ((main) * 1000 + (sub))
#define SPQR_MAIN_VERSION 1
#define SPQR_SUB_VERSION 2
#define SPQR_SUBSUB_VERSION 0
#define SPQR_VERSION SPQR_VER_CODE(SPQR_MAIN_VERSION,SPQR_SUB_VERSION)
#endif

604
spqr_mini/SuiteSparseQR_subset.hpp
View File

@ -1,604 +0,0 @@
// =============================================================================
// === SuiteSparseQR.hpp =======================================================
// =============================================================================
// User include file for C++ programs.
#ifndef SUITESPARSEQR_H
#define SUITESPARSEQR_H
// -----------------------------------------------------------------------------
// include files
// -----------------------------------------------------------------------------
#include "cholmod.h"
#include "UFconfig.h"
#include "SuiteSparseQR_definitions.h"
// =============================================================================
// === spqr_symbolic ===========================================================
// =============================================================================
// The contents of this object do not change during numeric factorization. The
// Symbolic object depends only on the pattern of the input matrix, and not its
// values. These contents also do not change with column pivoting for rank
// detection. This makes parallelism easier to manage, since all threads can
// have access to this object without synchronization.
//
// The total size of the Symbolic object is (10 + 2*m + anz + 2*n + 5*nf + rnz)
// Int's, where the user's input A matrix is m-by-n with anz nonzeros, nf <=
// MIN(m,n) is the number of frontal matrices, and rnz <= nnz(R) is the number
// of column indices used to represent the supernodal form of R (one Int per
// non-pivotal column index in the leading row of each block of R).
struct spqr_symbolic
{
// -------------------------------------------------------------------------
// row-form of the input matrix and its permutations
// -------------------------------------------------------------------------
// During symbolic analysis, the nonzero pattern of S = A(P,Q) is
// constructed, where A is the user's input matrix. Its numerical values
// are also constructed, but they do not become part of the Symbolic
// object. The matrix S is stored in row-oriented form. The rows of S are
// sorted according to their leftmost column index (via PLinv). Column
// indices in each row of S are in strictly ascending order, even though
// the input matrix A need not be sorted.
UF_long m, n, anz ; // S is m-by-n with anz entries
UF_long *Sp ; // size m+1, row pointers of S
UF_long *Sj ; // size anz = Sp [n], column indices of S
UF_long *Qfill ; // size n, fill-reducing column permutation.
// Qfill [k] = j if column k of A is column j of S.
UF_long *PLinv ; // size m, inverse row permutation that places S=A(P,Q)
// in increasing order of leftmost column index.
// PLinv [i] = k if row i of A is row k of S.
UF_long *Sleft ; // size n+2. The list of rows of S whose leftmost
// column index is j is given by
// Sleft [j] ... Sleft [j+1]-1. This can be empty (that is, Sleft
// [j] can equal Sleft [j+1]). Sleft [n] is the number of
// non-empty rows of S, and Sleft [n+1] == m. That is, Sleft [n]
// ... Sleft [n+1]-1 gives the empty rows of S, if any.
// -------------------------------------------------------------------------
// frontal matrices: pattern and tree
// -------------------------------------------------------------------------
// Each frontal matrix is fm-by-fn, with fnpiv pivot columns. The fn
// column indices are given by a set of size fnpiv pivot columns, defined
// by Super, followed by the pattern Rj [ Rp[f] ... Rp[f+1]-1 ].
// The row indices of the front are not kept. If the Householder vectors
// are not kept, the row indices are not needed. If the Householder
// vectors are kept, the row indices are computed dynamically during
// numerical factorization.
UF_long nf ; // number of frontal matrices; nf <= MIN (m,n)
UF_long maxfn ; // max # of columns in any front
// parent, child, and childp define the row merge tree or etree (A'A)
UF_long *Parent ; // size nf+1
UF_long *Child ; // size nf+1
UF_long *Childp ; // size nf+2
// The parent of a front f is Parent [f], or EMPTY if f=nf.
// A list of children of f can be obtained in the list
// Child [Childp [f] ... Childp [f+1]-1].
// Node nf in the tree is a placeholder; it does not represent a frontal
// matrix. All roots of the frontal "tree" (may be a forest) have the
// placeholder node nf as their parent. Thus, the tree of nodes 0:nf is
// truly a tree, with just one parent (node nf).
UF_long *Super ; // size nf+1. Super [f] gives the first pivot column
// in the front F. This refers to a column of S. The
// number of expected pivot columns in F is thus
// Super [f+1] - Super [f].
UF_long *Rp ; // size nf+1
UF_long *Rj ; // size rjsize; compressed supernodal form of R
UF_long *Post ; // size nf+1, post ordering of frontal tree. f=Post[k]
// gives the kth node in the postordered tree
UF_long rjsize ; // size of Rj
UF_long do_rank_detection ; // TRUE: allow for tol >= 0. FALSE: ignore tol
// the rest depends on whether or not rank-detection is allowed:
UF_long maxstack ; // max stack size (sequential case)
UF_long hisize ; // size of Hii
UF_long keepH ; // TRUE if H is present
UF_long *Hip ; // size nf+1. If H is kept, the row indices of frontal
// matrix f are in Hii [Hip [f] ... Hip [f] + Hm [f]],
// where Hii and Hm are stored in the numeric object.
// There is one block row of R per frontal matrix.
// The fn column indices of R are given by Rj [Rp [f] ... Rp [f+1]-1],
// where the first fp column indices are Super [f] ... Super [f+1]-1.
// The remaining column indices in Rj [...] are non-pivotal, and are
// in the range Super [f+1] to n. The number of rows of R is at
// most fp, but can be less if dead columns appear in the matrix.
// The number of columns in the contribution block C is always
// cn = fn - fp, where fn = Rp [f+1] - Rp [f].
UF_long ntasks ; // number of tasks in task graph
UF_long ns ; // number of stacks
// -------------------------------------------------------------------------
// the rest of the QR symbolic object is present only if ntasks > 1
// -------------------------------------------------------------------------
// Task tree (nodes 0:ntasks), including placeholder node
UF_long *TaskChildp ; // size ntasks+2
UF_long *TaskChild ; // size ntasks+1
UF_long *TaskStack ; // size ntasks+1
// list of fronts for each task
UF_long *TaskFront ; // size nf+1
UF_long *TaskFrontp ; // size ntasks+2
UF_long *On_stack ; // size nf+1, front f is on stack On_stack [f]
// size of each stack
UF_long *Stack_maxstack ; // size ns+2
} ;
// =============================================================================
// === spqr_numeric ============================================================
// =============================================================================
// The Numeric object contains the numerical values of the triangular/
// trapezoidal factor R, and optionally the Householder vectors H if they
// are kept.
template <typename Entry> struct spqr_numeric
{
// -------------------------------------------------------------------------
// Numeric R factor
// -------------------------------------------------------------------------
Entry **Rblock ; // size nf. R [f] is an (Entry *) pointer to the
// R block for front F. It is an upper trapezoidal
// of size Rm(f)-by-Rn(f), but only the upper
// triangular part is stored in column-packed format.
Entry **Stacks ; // size ns; an array of stacks holding the R and H
// factors and the current frontal matrix F at the head.
// This is followed by empty space, then the C blocks of
// prior frontal matrices at the bottom. When the
// factorization is complete, only the R and H part at
// the head of each stack is left.
UF_long *Stack_size ; // size ns; Stack_size [s] is the size of Stacks [s]
UF_long hisize ; // size of Hii
UF_long n ; // A is m-by-n
UF_long m ;
UF_long nf ; // number of frontal matrices
UF_long ntasks ; // number of tasks in task graph actually used
UF_long ns ; // number of stacks actually used
UF_long maxstack ; // size of sequential stack, if used
// -------------------------------------------------------------------------
// for rank detection and m < n case
// -------------------------------------------------------------------------
char *Rdead ; // size n, Rdead [k] = 1 if k is a dead pivot column,
// Rdead [k] = 0 otherwise. If no columns are dead,
// this is NULL. If m < n, then at least m-n columns
// will be dead.
UF_long rank ; // number of live pivot columns
UF_long rank1 ; // number of live pivot columns in first ntol columns
// of A
UF_long maxfrank ; // max number of rows in any R block
double norm_E_fro ; // 2-norm of w, the vector of dead column 2-norms
// -------------------------------------------------------------------------
// for keeping Householder vectors
// -------------------------------------------------------------------------
// The factorization is R = (H_s * ... * H_2 * H_1) * P_H
// where P_H is the permutation HPinv, and H_1, ... H_s are the Householder
// vectors (s = rjsize).
UF_long keepH ; // TRUE if H is present
UF_long rjsize ; // size of Hstair and HTau
UF_long *HStair ; // size rjsize. The list Hstair [Rp [f] ... Rp [f+1]-1]
// gives the staircase for front F
Entry *HTau ; // size rjsize. The list HTau [Rp [f] ... Rp [f+1]-1]
// gives the Householder coefficients for front F
UF_long *Hii ; // size hisize, row indices of H.
UF_long *HPinv ; // size m. HPinv [i] = k if row i of A and H is row k
// of R. This permutation includes QRsym->PLinv, and
// the permutation constructed via pivotal row ordering
// during factorization.
UF_long *Hm ; // size nf, Hm [f] = # of rows in front F
UF_long *Hr ; // size nf, Hr [f] = # of rows in R block of front F
UF_long maxfm ; // max (Hm [0:nf-1]), computed only if H kept
} ;
// =============================================================================
// === SuiteSparseQR_factorization =============================================
// =============================================================================
// A combined symbolic+numeric QR factorization of A or [A B],
// with singletons
template <typename Entry> struct SuiteSparseQR_factorization
{
// QR factorization of A or [A Binput] after singletons have been removed
double tol ; // tol used
spqr_symbolic *QRsym ;
spqr_numeric <Entry> *QRnum ;
// singletons, in compressed-row form; R is n1rows-by-n
UF_long *R1p ; // size n1rows+1
UF_long *R1j ;
Entry *R1x ;
UF_long r1nz ; // nnz (R1)
// combined singleton and fill-reducing permutation
UF_long *Q1fill ;
UF_long *P1inv ;
UF_long *HP1inv ; // NULL if n1cols == 0, in which case QRnum->HPinv
// serves in its place.
// Rmap and RmapInv are NULL if QR->rank == A->ncol
UF_long *Rmap ; // size n. Rmap [j] = k if column j of R is the kth
// live column and where k < QR->rank; otherwise, if
// j is a dead column, then k >= QR->rank.
UF_long *RmapInv ;
UF_long n1rows ; // number of singleton rows of [A B]
UF_long n1cols ; // number of singleton columns of [A B]
UF_long narows ; // number of rows of A
UF_long nacols ; // number of columns of A
UF_long bncols ; // number of columns of B
UF_long rank ; // rank estimate of A (n1rows + QRnum->rank1), ranges
// from 0 to min(m,n)
int allow_tol ; // if TRUE, do rank detection
} ;
// =============================================================================
// === Simple user-callable SuiteSparseQR functions ============================
// =============================================================================
// SuiteSparseQR Sparse QR factorization and solve
// SuiteSparseQR_qmult Q*X, Q'*X, X*Q, or X*Q' for X full or sparse
// returns rank(A) estimate, or EMPTY on failure
template <typename Entry> UF_long SuiteSparseQR
(
// inputs, not modified
int ordering, // all, except 3:given treated as 0:fixed
double tol, // only accept singletons above tol
UF_long econ, // number of rows of C and R to return; a value
// less than the rank r of A is treated as r, and
// a value greater than m is treated as m.
int getCTX, // if 0: return Z = C of size econ-by-bncols
// if 1: return Z = C' of size bncols-by-econ
// if 2: return Z = X of size econ-by-bncols
cholmod_sparse *A, // m-by-n sparse matrix
// B is either sparse or dense. If Bsparse is non-NULL, B is sparse and
// Bdense is ignored. If Bsparse is NULL and Bdense is non-NULL, then B is
// dense. B is not present if both are NULL.
cholmod_sparse *Bsparse,
cholmod_dense *Bdense,
// output arrays, neither allocated nor defined on input.
// Z is the matrix C, C', or X
cholmod_sparse **Zsparse,
cholmod_dense **Zdense,
cholmod_sparse **R, // the R factor
UF_long **E, // size n; fill-reducing ordering of A.
cholmod_sparse **H, // the Householder vectors (m-by-nh)
UF_long **HPinv, // size m; row permutation for H
cholmod_dense **HTau, // size nh, Householder coefficients
// workspace and parameters
cholmod_common *cc
) ;
// X = A\dense(B)
template <typename Entry> cholmod_dense *SuiteSparseQR
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
cholmod_sparse *A, // m-by-n sparse matrix
cholmod_dense *B, // m-by-nrhs
cholmod_common *cc // workspace and parameters
) ;
// X = A\dense(B) using default ordering and tolerance
template <typename Entry> cholmod_dense *SuiteSparseQR
(
cholmod_sparse *A, // m-by-n sparse matrix
cholmod_dense *B, // m-by-nrhs
cholmod_common *cc // workspace and parameters
) ;
// X = A\sparse(B)
template <typename Entry> cholmod_sparse *SuiteSparseQR
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
cholmod_sparse *A, // m-by-n sparse matrix
cholmod_sparse *B, // m-by-nrhs
cholmod_common *cc // workspace and parameters
) ;
// [Q,R,E] = qr(A), returning Q as a sparse matrix
template <typename Entry> UF_long SuiteSparseQR // returns rank(A) estimate
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
UF_long econ,
cholmod_sparse *A, // m-by-n sparse matrix
// outputs
cholmod_sparse **Q, // m-by-e sparse matrix where e=max(econ,rank(A))
cholmod_sparse **R, // e-by-n sparse matrix
UF_long **E, // permutation of 0:n-1, NULL if identity
cholmod_common *cc // workspace and parameters
) ;
// [Q,R,E] = qr(A), discarding Q
template <typename Entry> UF_long SuiteSparseQR // returns rank(A) estimate
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
UF_long econ,
cholmod_sparse *A, // m-by-n sparse matrix
// outputs
cholmod_sparse **R, // e-by-n sparse matrix
UF_long **E, // permutation of 0:n-1, NULL if identity
cholmod_common *cc // workspace and parameters
) ;
// [C,R,E] = qr(A,B), where C and B are dense
template <typename Entry> UF_long SuiteSparseQR
(
// inputs, not modified
int ordering, // all, except 3:given treated as 0:fixed
double tol, // only accept singletons above tol
UF_long econ, // number of rows of C and R to return
cholmod_sparse *A, // m-by-n sparse matrix
cholmod_dense *B, // m-by-nrhs dense matrix
// outputs
cholmod_dense **C, // C = Q'*B, an e-by-nrhs dense matrix
cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
UF_long **E, // permutation of 0:n-1, NULL if identity
cholmod_common *cc // workspace and parameters
) ;
// [C,R,E] = qr(A,B), where C and B are sparse
template <typename Entry> UF_long SuiteSparseQR
(
// inputs, not modified
int ordering, // all, except 3:given treated as 0:fixed
double tol, // only accept singletons above tol
UF_long econ, // number of rows of C and R to return
cholmod_sparse *A, // m-by-n sparse matrix
cholmod_sparse *B, // m-by-nrhs sparse matrix
// outputs
cholmod_sparse **C, // C = Q'*B, an e-by-nrhs sparse matrix
cholmod_sparse **R, // e-by-n sparse matrix where e=max(econ,rank(A))
UF_long **E, // permutation of 0:n-1, NULL if identity
cholmod_common *cc // workspace and parameters
) ;
// [Q,R,E] = qr(A) where Q is returned in Householder form
template <typename Entry> UF_long SuiteSparseQR
(
// inputs, not modified
int ordering, // all, except 3:given treated as 0:fixed
double tol, // only accept singletons above tol
UF_long econ, // number of rows of C and R to return
cholmod_sparse *A, // m-by-n sparse matrix
// outputs
cholmod_sparse **R, // the R factor
UF_long **E, // permutation of 0:n-1, NULL if identity
cholmod_sparse **H, // the Householder vectors (m-by-nh)
UF_long **HPinv, // size m; row permutation for H
cholmod_dense **HTau, // size nh, Householder coefficients
cholmod_common *cc // workspace and parameters
) ;
// =============================================================================
// === SuiteSparseQR_qmult =====================================================
// =============================================================================
// This function takes as input the matrix Q in Householder form, as returned
// by SuiteSparseQR (... H, HPinv, HTau, cc) above.
// returns Y of size m-by-n (NULL on failure)
template <typename Entry> cholmod_dense *SuiteSparseQR_qmult
(
// inputs, no modified
int method, // 0,1,2,3
cholmod_sparse *H, // either m-by-nh or n-by-nh
cholmod_dense *HTau, // size 1-by-nh
UF_long *HPinv, // size mh
cholmod_dense *Xdense, // size m-by-n
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> cholmod_sparse *SuiteSparseQR_qmult
(
// inputs, no modified
int method, // 0,1,2,3
cholmod_sparse *H, // either m-by-nh or n-by-nh
cholmod_dense *HTau, // size 1-by-nh
UF_long *HPinv, // size mh
cholmod_sparse *X,
// workspace and parameters
cholmod_common *cc
) ;
// =============================================================================
// === Expert user-callable SuiteSparseQR functions ============================
// =============================================================================
#ifndef NEXPERT
// These functions are "expert" routines, allowing reuse of the QR
// factorization for different right-hand-sides. They also allow the user to
// find the minimum 2-norm solution to an undertermined system of equations.
template <typename Entry>
SuiteSparseQR_factorization <Entry> *SuiteSparseQR_factorize
(
// inputs, not modified:
int ordering, // all, except 3:given treated as 0:fixed
double tol, // treat columns with 2-norm <= tol as zero
cholmod_sparse *A, // sparse matrix to factorize
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> cholmod_dense *SuiteSparseQR_solve // returns X
(
// inputs, not modified:
int system, // which system to solve
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
cholmod_dense *B, // right-hand-side, m-by-nrhs or n-by-nrhs
// workspace and parameters
cholmod_common *cc
) ;
template <typename Entry> cholmod_sparse *SuiteSparseQR_solve // returns X
(
// inputs, not modified:
int system, // which system to solve (0,1,2,3)
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
cholmod_sparse *Bsparse, // right-hand-side, m-by-nrhs or n-by-nrhs
// workspace and parameters
cholmod_common *cc
) ;
// returns Y of size m-by-n, or NULL on failure
template <typename Entry> cholmod_dense *SuiteSparseQR_qmult
(
// inputs, not modified
int method, // 0,1,2,3 (same as SuiteSparseQR_qmult)
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
cholmod_dense *Xdense, // size m-by-n with leading dimension ldx
// workspace and parameters
cholmod_common *cc
) ;
// returns Y of size m-by-n, or NULL on failure
template <typename Entry> cholmod_sparse *SuiteSparseQR_qmult
(
// inputs, not modified
int method, // 0,1,2,3
SuiteSparseQR_factorization <Entry> *QR, // of an m-by-n sparse matrix A
cholmod_sparse *Xsparse, // size m-by-n
// workspace and parameters
cholmod_common *cc
) ;
// free the QR object
template <typename Entry> int SuiteSparseQR_free
(
SuiteSparseQR_factorization <Entry> **QR, // of an m-by-n sparse matrix A
cholmod_common *cc
) ;
// find the min 2-norm solution to a sparse linear system
template <typename Entry> cholmod_dense *SuiteSparseQR_min2norm
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
cholmod_sparse *A,
cholmod_dense *B,
cholmod_common *cc
) ;
template <typename Entry> cholmod_sparse *SuiteSparseQR_min2norm
(
int ordering, // all, except 3:given treated as 0:fixed
double tol,
cholmod_sparse *A,
cholmod_sparse *B,
cholmod_common *cc
) ;
// symbolic QR factorization; no singletons exploited
template <typename Entry>
SuiteSparseQR_factorization <Entry> *SuiteSparseQR_symbolic
(
// inputs:
int ordering, // all, except 3:given treated as 0:fixed
int allow_tol, // if FALSE, tol is ignored by the numeric
// factorization, and no rank detection is performed
cholmod_sparse *A, // sparse matrix to factorize (A->x ignored)
cholmod_common *cc // workspace and parameters
) ;
// numeric QR factorization;
template <typename Entry> int SuiteSparseQR_numeric
(
// inputs:
double tol, // treat columns with 2-norm <= tol as zero
cholmod_sparse *A, // sparse matrix to factorize
// input/output
SuiteSparseQR_factorization <Entry> *QR,
cholmod_common *cc // workspace and parameters
) ;
#endif
// =============================================================================
// === high-resolution timing ==================================================
// =============================================================================
#ifdef TIMING
extern "C" {
#include <time.h>
#include <sys/time.h>
double spqr_time ( ) ; // returns current time in seconds
}
#endif
#endif

151
spqr_mini/UFconfig.h
View File

@ -1,151 +0,0 @@
/* ========================================================================== */
/* === UFconfig.h =========================================================== */
/* ========================================================================== */
/* Configuration file for SuiteSparse: a Suite of Sparse matrix packages
* (AMD, COLAMD, CCOLAMD, CAMD, CHOLMOD, UMFPACK, CXSparse, and others).
*
* UFconfig.h provides the definition of the long integer. On most systems,
* a C program can be compiled in LP64 mode, in which long's and pointers are
* both 64-bits, and int's are 32-bits. Windows 64, however, uses the LLP64
* model, in which int's and long's are 32-bits, and long long's and pointers
* are 64-bits.
*
* SuiteSparse packages that include long integer versions are
* intended for the LP64 mode. However, as a workaround for Windows 64
* (and perhaps other systems), the long integer can be redefined.
*
* If _WIN64 is defined, then the __int64 type is used instead of long.
*
* The long integer can also be defined at compile time. For example, this
* could be added to UFconfig.mk:
*
* CFLAGS = -O -D'UF_long=long long' -D'UF_long_max=9223372036854775801' \
* -D'UF_long_idd="lld"'
*
* This file defines UF_long as either long (on all but _WIN64) or
* __int64 on Windows 64. The intent is that a UF_long is always a 64-bit
* integer in a 64-bit code. ptrdiff_t might be a better choice than long;
* it is always the same size as a pointer.
*
* This file also defines the SUITESPARSE_VERSION and related definitions.
*
* Copyright (c) 2007, University of Florida. No licensing restrictions
* apply to this file or to the UFconfig directory. Author: Timothy A. Davis.
*/
#ifndef _UFCONFIG_H
#define _UFCONFIG_H
#ifdef __cplusplus
extern "C" {
#endif
#include <limits.h>
#include <stdlib.h>
/* ========================================================================== */
/* === UF_long ============================================================== */
/* ========================================================================== */
#ifndef UF_long
#ifdef _WIN64
#define UF_long __int64
#define UF_long_max _I64_MAX
#define UF_long_idd "I64d"
#else
#define UF_long long
#define UF_long_max LONG_MAX
#define UF_long_idd "ld"
#endif
#define UF_long_id "%" UF_long_idd
#endif
/* ========================================================================== */
/* === UFconfig parameters and functions ==================================== */
/* ========================================================================== */
/* SuiteSparse-wide parameters will be placed in this struct. So far, they
are only used by RBio. */
typedef struct UFconfig_struct
{
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 */
} UFconfig ;
void *UFmalloc /* pointer to allocated block of memory */
(
size_t nitems, /* number of items to malloc (>=1 is enforced) */
size_t size_of_item, /* sizeof each item */
int *ok, /* TRUE if successful, FALSE otherwise */
UFconfig *config /* SuiteSparse-wide configuration */
) ;
void *UFfree /* always returns NULL */
(
void *p, /* block to free */
UFconfig *config /* SuiteSparse-wide configuration */
) ;
/* ========================================================================== */
/* === SuiteSparse version ================================================== */
/* ========================================================================== */
/* SuiteSparse is not a package itself, but a collection of packages, some of
* which must be used together (UMFPACK requires AMD, CHOLMOD requires AMD,
* COLAMD, CAMD, and CCOLAMD, etc). A version number is provided here for the
* collection itself. The versions of packages within each version of
* SuiteSparse are meant to work together. Combining one packge from one
* version of SuiteSparse, with another package from another version of
* SuiteSparse, may or may not work.
*
* SuiteSparse Version 3.5.0 contains the following packages:
*
* AMD version 2.2.1
* BTF version 1.1.1
* CAMD version 2.2.1
* CCOLAMD version 2.7.2
* CHOLMOD version 1.7.2
* COLAMD version 2.7.2
* CSparse version 2.2.4
* CXSparse version 2.2.4
* KLU version 1.1.1
* LDL version 2.0.2
* RBio version 2.0.0
* SuiteSparseQR version 1.2.0
* UFcollection version 1.3.0
* UFconfig version number is the same as SuiteSparse
* UMFPACK version 5.5.0
* LINFACTOR version 1.1.0
* MESHND version 1.1.1
* SSMULT version 2.0.2
* MATLAB_Tools no specific version number
*
* Other package dependencies:
* BLAS required by CHOLMOD and UMFPACK
* LAPACK required by CHOLMOD
* METIS 4.0.1 required by CHOLMOD (optional) and KLU (optional)
*/
#define SUITESPARSE_DATE "Nov 30, 2009"
#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub))
#define SUITESPARSE_MAIN_VERSION 3
#define SUITESPARSE_SUB_VERSION 5
#define SUITESPARSE_SUBSUB_VERSION 1
#define SUITESPARSE_VERSION \
SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION)
#ifdef __cplusplus
}
#endif
#endif

456
spqr_mini/cholmod_blas.h
View File

@ -1,456 +0,0 @@
/* ========================================================================== */
/* === Include/cholmod_blas.h =============================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Include/cholmod_blas.h.
* Copyright (C) 2005-2006, Univ. of Florida. Author: Timothy A. Davis
* CHOLMOD/Include/cholmod_blas.h is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* This does not need to be included in the user's program. */
#ifndef CHOLMOD_BLAS_H
#define CHOLMOD_BLAS_H
/* ========================================================================== */
/* === Architecture ========================================================= */
/* ========================================================================== */
#if defined (__sun) || defined (MSOL2) || defined (ARCH_SOL2)
#define CHOLMOD_SOL2
#define CHOLMOD_ARCHITECTURE "Sun Solaris"
#elif defined (__sgi) || defined (MSGI) || defined (ARCH_SGI)
#define CHOLMOD_SGI
#define CHOLMOD_ARCHITECTURE "SGI Irix"
#elif defined (__linux) || defined (MGLNX86) || defined (ARCH_GLNX86)
#define CHOLMOD_LINUX
#define CHOLMOD_ARCHITECTURE "Linux"
#elif defined (__APPLE__)
#define CHOLMOD_MAC
#define CHOLMOD_ARCHITECTURE "Mac"
#elif defined (_AIX) || defined (MIBM_RS) || defined (ARCH_IBM_RS)
#define CHOLMOD_AIX
#define CHOLMOD_ARCHITECTURE "IBM AIX"
/* recent reports from IBM AIX seem to indicate that this is not needed: */
/* #define BLAS_NO_UNDERSCORE */
#elif defined (__alpha) || defined (MALPHA) || defined (ARCH_ALPHA)
#define CHOLMOD_ALPHA
#define CHOLMOD_ARCHITECTURE "Compaq Alpha"
#elif defined (_WIN32) || defined (WIN32) || defined (_WIN64) || defined (WIN64)
#if defined (__MINGW32__) || defined (__MINGW32__)
#define CHOLMOD_MINGW
#elif defined (__CYGWIN32__) || defined (__CYGWIN32__)
#define CHOLMOD_CYGWIN
#else
#define CHOLMOD_WINDOWS
#define BLAS_NO_UNDERSCORE
#endif
#define CHOLMOD_ARCHITECTURE "Microsoft Windows"
#elif defined (__hppa) || defined (__hpux) || defined (MHPUX) || defined (ARCH_HPUX)
#define CHOLMOD_HP
#define CHOLMOD_ARCHITECTURE "HP Unix"
#define BLAS_NO_UNDERSCORE
#elif defined (__hp700) || defined (MHP700) || defined (ARCH_HP700)
#define CHOLMOD_HP
#define CHOLMOD_ARCHITECTURE "HP 700 Unix"
#define BLAS_NO_UNDERSCORE
#else
/* If the architecture is unknown, and you call the BLAS, you may need to */
/* define BLAS_BY_VALUE, BLAS_NO_UNDERSCORE, and/or BLAS_CHAR_ARG yourself. */
#define CHOLMOD_ARCHITECTURE "unknown"
#endif
/* ========================================================================== */
/* === BLAS and LAPACK names ================================================ */
/* ========================================================================== */
/* Prototypes for the various versions of the BLAS. */
/* Determine if the 64-bit Sun Performance BLAS is to be used */
#if defined(CHOLMOD_SOL2) && !defined(NSUNPERF) && defined(BLAS64)
#define SUN64
#endif
#ifdef SUN64
#define BLAS_DTRSV dtrsv_64_
#define BLAS_DGEMV dgemv_64_
#define BLAS_DTRSM dtrsm_64_
#define BLAS_DGEMM dgemm_64_
#define BLAS_DSYRK dsyrk_64_
#define BLAS_DGER dger_64_
#define BLAS_DSCAL dscal_64_
#define LAPACK_DPOTRF dpotrf_64_
#define BLAS_ZTRSV ztrsv_64_
#define BLAS_ZGEMV zgemv_64_
#define BLAS_ZTRSM ztrsm_64_
#define BLAS_ZGEMM zgemm_64_
#define BLAS_ZHERK zherk_64_
#define BLAS_ZGER zgeru_64_
#define BLAS_ZSCAL zscal_64_
#define LAPACK_ZPOTRF zpotrf_64_
#elif defined (BLAS_NO_UNDERSCORE)
#define BLAS_DTRSV dtrsv
#define BLAS_DGEMV dgemv
#define BLAS_DTRSM dtrsm
#define BLAS_DGEMM dgemm
#define BLAS_DSYRK dsyrk
#define BLAS_DGER dger
#define BLAS_DSCAL dscal
#define LAPACK_DPOTRF dpotrf
#define BLAS_ZTRSV ztrsv
#define BLAS_ZGEMV zgemv
#define BLAS_ZTRSM ztrsm
#define BLAS_ZGEMM zgemm
#define BLAS_ZHERK zherk
#define BLAS_ZGER zgeru
#define BLAS_ZSCAL zscal
#define LAPACK_ZPOTRF zpotrf
#else
#define BLAS_DTRSV dtrsv_
#define BLAS_DGEMV dgemv_
#define BLAS_DTRSM dtrsm_
#define BLAS_DGEMM dgemm_
#define BLAS_DSYRK dsyrk_
#define BLAS_DGER dger_
#define BLAS_DSCAL dscal_
#define LAPACK_DPOTRF dpotrf_
#define BLAS_ZTRSV ztrsv_
#define BLAS_ZGEMV zgemv_
#define BLAS_ZTRSM ztrsm_
#define BLAS_ZGEMM zgemm_
#define BLAS_ZHERK zherk_
#define BLAS_ZGER zgeru_
#define BLAS_ZSCAL zscal_
#define LAPACK_ZPOTRF zpotrf_
#endif
/* ========================================================================== */
/* === BLAS and LAPACK integer arguments ==================================== */
/* ========================================================================== */
/* Compile CHOLMOD, UMFPACK, and SPQR with -DBLAS64 if you have a BLAS that
* uses 64-bit integers */
#if defined (LONGBLAS) || defined (BLAS64)
#define BLAS_INT UF_long
#else
#define BLAS_INT int
#endif
/* If the BLAS integer is smaller than the basic CHOLMOD integer, then we need
* to check for integer overflow when converting from Int to BLAS_INT. If
* any integer overflows, the externally-defined BLAS_OK variable is
* set to FALSE. BLAS_OK should be set to TRUE before calling any
* BLAS_* macro.
*/
#define CHECK_BLAS_INT (sizeof (BLAS_INT) < sizeof (Int))
#define EQ(K,k) (((BLAS_INT) K) == ((Int) k))
/* ========================================================================== */
/* === BLAS and LAPACK prototypes and macros ================================ */
/* ========================================================================== */
void BLAS_DGEMV (char *trans, BLAS_INT *m, BLAS_INT *n, double *alpha,
double *A, BLAS_INT *lda, double *X, BLAS_INT *incx, double *beta,
double *Y, BLAS_INT *incy) ;
#define BLAS_dgemv(trans,m,n,alpha,A,lda,X,incx,beta,Y,incy) \
{ \
BLAS_INT M = m, N = n, LDA = lda, INCX = incx, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (INCX,incx) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DGEMV (trans, &M, &N, alpha, A, &LDA, X, &INCX, beta, Y, &INCY) ; \
} \
}
void BLAS_ZGEMV (char *trans, BLAS_INT *m, BLAS_INT *n, double *alpha,
double *A, BLAS_INT *lda, double *X, BLAS_INT *incx, double *beta,
double *Y, BLAS_INT *incy) ;
#define BLAS_zgemv(trans,m,n,alpha,A,lda,X,incx,beta,Y,incy) \
{ \
BLAS_INT M = m, N = n, LDA = lda, INCX = incx, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (INCX,incx) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZGEMV (trans, &M, &N, alpha, A, &LDA, X, &INCX, beta, Y, &INCY) ; \
} \
}
void BLAS_DTRSV (char *uplo, char *trans, char *diag, BLAS_INT *n, double *A,
BLAS_INT *lda, double *X, BLAS_INT *incx) ;
#define BLAS_dtrsv(uplo,trans,diag,n,A,lda,X,incx) \
{ \
BLAS_INT N = n, LDA = lda, INCX = incx ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (LDA,lda) && EQ (INCX,incx))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DTRSV (uplo, trans, diag, &N, A, &LDA, X, &INCX) ; \
} \
}
void BLAS_ZTRSV (char *uplo, char *trans, char *diag, BLAS_INT *n, double *A,
BLAS_INT *lda, double *X, BLAS_INT *incx) ;
#define BLAS_ztrsv(uplo,trans,diag,n,A,lda,X,incx) \
{ \
BLAS_INT N = n, LDA = lda, INCX = incx ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (LDA,lda) && EQ (INCX,incx))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZTRSV (uplo, trans, diag, &N, A, &LDA, X, &INCX) ; \
} \
}
void BLAS_DTRSM (char *side, char *uplo, char *transa, char *diag, BLAS_INT *m,
BLAS_INT *n, double *alpha, double *A, BLAS_INT *lda, double *B,
BLAS_INT *ldb) ;
#define BLAS_dtrsm(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb) \
{ \
BLAS_INT M = m, N = n, LDA = lda, LDB = ldb ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (LDB,ldb))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DTRSM (side, uplo, transa, diag, &M, &N, alpha, A, &LDA, B, &LDB);\
} \
}
void BLAS_ZTRSM (char *side, char *uplo, char *transa, char *diag, BLAS_INT *m,
BLAS_INT *n, double *alpha, double *A, BLAS_INT *lda, double *B,
BLAS_INT *ldb) ;
#define BLAS_ztrsm(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb) \
{ \
BLAS_INT M = m, N = n, LDA = lda, LDB = ldb ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (LDB,ldb))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZTRSM (side, uplo, transa, diag, &M, &N, alpha, A, &LDA, B, &LDB);\
} \
}
void BLAS_DGEMM (char *transa, char *transb, BLAS_INT *m, BLAS_INT *n,
BLAS_INT *k, double *alpha, double *A, BLAS_INT *lda, double *B,
BLAS_INT *ldb, double *beta, double *C, BLAS_INT *ldc) ;
#define BLAS_dgemm(transa,transb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc) \
{ \
BLAS_INT M = m, N = n, K = k, LDA = lda, LDB = ldb, LDC = ldc ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (K,k) && \
EQ (LDA,lda) && EQ (LDB,ldb) && EQ (LDC,ldc))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DGEMM (transa, transb, &M, &N, &K, alpha, A, &LDA, B, &LDB, beta, \
C, &LDC) ; \
} \
}
void BLAS_ZGEMM (char *transa, char *transb, BLAS_INT *m, BLAS_INT *n,
BLAS_INT *k, double *alpha, double *A, BLAS_INT *lda, double *B,
BLAS_INT *ldb, double *beta, double *C, BLAS_INT *ldc) ;
#define BLAS_zgemm(transa,transb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc) \
{ \
BLAS_INT M = m, N = n, K = k, LDA = lda, LDB = ldb, LDC = ldc ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (K,k) && \
EQ (LDA,lda) && EQ (LDB,ldb) && EQ (LDC,ldc))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZGEMM (transa, transb, &M, &N, &K, alpha, A, &LDA, B, &LDB, beta, \
C, &LDC) ; \
} \
}
void BLAS_DSYRK (char *uplo, char *trans, BLAS_INT *n, BLAS_INT *k,
double *alpha, double *A, BLAS_INT *lda, double *beta, double *C,
BLAS_INT *ldc) ;
#define BLAS_dsyrk(uplo,trans,n,k,alpha,A,lda,beta,C,ldc) \
{ \
BLAS_INT N = n, K = k, LDA = lda, LDC = ldc ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (K,k) && EQ (LDA,lda) && \
EQ (LDC,ldc))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DSYRK (uplo, trans, &N, &K, alpha, A, &LDA, beta, C, &LDC) ; \
} \
} \
void BLAS_ZHERK (char *uplo, char *trans, BLAS_INT *n, BLAS_INT *k,
double *alpha, double *A, BLAS_INT *lda, double *beta, double *C,
BLAS_INT *ldc) ;
#define BLAS_zherk(uplo,trans,n,k,alpha,A,lda,beta,C,ldc) \
{ \
BLAS_INT N = n, K = k, LDA = lda, LDC = ldc ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (K,k) && EQ (LDA,lda) && \
EQ (LDC,ldc))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZHERK (uplo, trans, &N, &K, alpha, A, &LDA, beta, C, &LDC) ; \
} \
} \
void LAPACK_DPOTRF (char *uplo, BLAS_INT *n, double *A, BLAS_INT *lda,
BLAS_INT *info) ;
#define LAPACK_dpotrf(uplo,n,A,lda,info) \
{ \
BLAS_INT N = n, LDA = lda, INFO = 1 ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (LDA,lda))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
LAPACK_DPOTRF (uplo, &N, A, &LDA, &INFO) ; \
} \
info = INFO ; \
}
void LAPACK_ZPOTRF (char *uplo, BLAS_INT *n, double *A, BLAS_INT *lda,
BLAS_INT *info) ;
#define LAPACK_zpotrf(uplo,n,A,lda,info) \
{ \
BLAS_INT N = n, LDA = lda, INFO = 1 ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (LDA,lda))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
LAPACK_ZPOTRF (uplo, &N, A, &LDA, &INFO) ; \
} \
info = INFO ; \
}
/* ========================================================================== */
void BLAS_DSCAL (BLAS_INT *n, double *alpha, double *Y, BLAS_INT *incy) ;
#define BLAS_dscal(n,alpha,Y,incy) \
{ \
BLAS_INT N = n, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DSCAL (&N, alpha, Y, &INCY) ; \
} \
}
void BLAS_ZSCAL (BLAS_INT *n, double *alpha, double *Y, BLAS_INT *incy) ;
#define BLAS_zscal(n,alpha,Y,incy) \
{ \
BLAS_INT N = n, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (N,n) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZSCAL (&N, alpha, Y, &INCY) ; \
} \
}
void BLAS_DGER (BLAS_INT *m, BLAS_INT *n, double *alpha,
double *X, BLAS_INT *incx, double *Y, BLAS_INT *incy,
double *A, BLAS_INT *lda) ;
#define BLAS_dger(m,n,alpha,X,incx,Y,incy,A,lda) \
{ \
BLAS_INT M = m, N = n, LDA = lda, INCX = incx, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (INCX,incx) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_DGER (&M, &N, alpha, X, &INCX, Y, &INCY, A, &LDA) ; \
} \
}
void BLAS_ZGER (BLAS_INT *m, BLAS_INT *n, double *alpha,
double *X, BLAS_INT *incx, double *Y, BLAS_INT *incy,
double *A, BLAS_INT *lda) ;
#define BLAS_zgeru(m,n,alpha,X,incx,Y,incy,A,lda) \
{ \
BLAS_INT M = m, N = n, LDA = lda, INCX = incx, INCY = incy ; \
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDA,lda) && \
EQ (INCX,incx) && EQ (INCY,incy))) \
{ \
BLAS_OK = FALSE ; \
} \
if (!CHECK_BLAS_INT || BLAS_OK) \
{ \
BLAS_ZGER (&M, &N, alpha, X, &INCX, Y, &INCY, A, &LDA) ; \
} \
}
#endif

672
spqr_mini/cholmod_common.c
View File

@ -1,672 +0,0 @@
/* ========================================================================== */
/* === Core/cholmod_common ================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Core Module. Copyright (C) 2005-2006,
* Univ. of Florida. Author: Timothy A. Davis
* The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Core utility routines for the cholmod_common object:
*
* Primary routines:
* -----------------
* cholmod_start the first call to CHOLMOD
* cholmod_finish the last call to CHOLMOD
*
* Secondary routines:
* -------------------
* cholmod_defaults restore (most) default control parameters
* cholmod_allocate_work allocate (or reallocate) workspace in Common
* cholmod_free_work free workspace in Common
* cholmod_clear_flag clear Common->Flag in workspace
* cholmod_maxrank column dimension of Common->Xwork workspace
*
* The Common object is unique. It cannot be allocated or deallocated by
* CHOLMOD, since it contains the definition of the memory management routines
* used (pointers to malloc, free, realloc, and calloc, or their equivalent).
* The Common object contains workspace that is used between calls to
* CHOLMOD routines. This workspace allocated by CHOLMOD as needed, by
* cholmod_allocate_work and cholmod_free_work.
*/
#include "cholmod_internal.h"
#include "cholmod_core.h"
/* ========================================================================== */
/* === cholmod_start ======================================================== */
/* ========================================================================== */
/* Initialize Common default parameters and statistics. Sets workspace
* pointers to NULL.
*
* This routine must be called just once, prior to calling any other CHOLMOD
* routine. Do not call this routine after any other CHOLMOD routine (except
* cholmod_finish, to start a new CHOLMOD session), or a memory leak will
* occur.
*
* workspace: none
*/
int CHOLMOD(start)
(
cholmod_common *Common
)
{
int k ;
if (Common == NULL)
{
return (FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* user error handling routine */
/* ---------------------------------------------------------------------- */
Common->error_handler = NULL ;
/* ---------------------------------------------------------------------- */
/* integer and numerical types */
/* ---------------------------------------------------------------------- */
Common->itype = ITYPE ;
Common->dtype = DTYPE ;
/* ---------------------------------------------------------------------- */
/* default control parameters */
/* ---------------------------------------------------------------------- */
cholmod_l_defaults(Common) ;
Common->try_catch = FALSE ;
/* ---------------------------------------------------------------------- */
/* memory management routines */
/* ---------------------------------------------------------------------- */
/* The user can replace cholmod's memory management routines by redefining
* these function pointers. */
#ifndef NMALLOC
/* stand-alone ANSI C program */
Common->malloc_memory = malloc ;
Common->free_memory = free ;
Common->realloc_memory = realloc ;
Common->calloc_memory = calloc ;
#else
/* no memory manager defined at compile-time; MUST define one at run-time */
Common->malloc_memory = NULL ;
Common->free_memory = NULL ;
Common->realloc_memory = NULL ;
Common->calloc_memory = NULL ;
#endif
/* ---------------------------------------------------------------------- */
/* complex arithmetic routines */
/* ---------------------------------------------------------------------- */
// Common->complex_divide = CHOLMOD(divcomplex) ;
// Common->hypotenuse = CHOLMOD(hypot) ;
/* ---------------------------------------------------------------------- */
/* print routine */
/* ---------------------------------------------------------------------- */
#ifndef NPRINT
/* stand-alone ANSI C program */
Common->print_function = printf ;
#else
/* printing disabled */
Common->print_function = NULL ;
#endif
/* ---------------------------------------------------------------------- */
/* workspace */
/* ---------------------------------------------------------------------- */
/* This code assumes the workspace held in Common is not initialized. If
* it is, then a memory leak will occur because the pointers are
* overwritten with NULL. */
Common->nrow = 0 ;
Common->mark = EMPTY ;
Common->xworksize = 0 ;
Common->iworksize = 0 ;
Common->Flag = NULL ;
Common->Head = NULL ;
Common->Iwork = NULL ;
Common->Xwork = NULL ;
Common->no_workspace_reallocate = FALSE ;
/* ---------------------------------------------------------------------- */
/* statistics */
/* ---------------------------------------------------------------------- */
/* fl and lnz are computed in cholmod_analyze and cholmod_rowcolcounts */
Common->fl = EMPTY ;
Common->lnz = EMPTY ;
/* modfl is computed in cholmod_updown, cholmod_rowadd, and cholmod_rowdel*/
Common->modfl = EMPTY ;
/* all routines use status as their error-report code */
Common->status = CHOLMOD_OK ;
Common->malloc_count = 0 ; /* # calls to malloc minus # calls to free */
Common->memory_usage = 0 ; /* peak memory usage (in bytes) */
Common->memory_inuse = 0 ; /* current memory in use (in bytes) */
Common->nrealloc_col = 0 ;
Common->nrealloc_factor = 0 ;
Common->ndbounds_hit = 0 ;
Common->rowfacfl = 0 ;
Common->aatfl = EMPTY ;
/* Common->called_nd is TRUE if cholmod_analyze called or NESDIS */
Common->called_nd = FALSE ;
Common->blas_ok = TRUE ; /* false if BLAS int overflow occurs */
/* ---------------------------------------------------------------------- */
/* default SuiteSparseQR knobs and statististics */
/* ---------------------------------------------------------------------- */
for (k = 0 ; k < 4 ; k++) Common->SPQR_xstat [k] = 0 ;
for (k = 0 ; k < 10 ; k++) Common->SPQR_istat [k] = 0 ;
Common->SPQR_grain = 1 ; /* no Intel TBB multitasking, by default */
Common->SPQR_small = 1e6 ; /* target min task size for TBB */
Common->SPQR_shrink = 1 ; /* controls SPQR shrink realloc */
Common->SPQR_nthreads = 0 ; /* 0: let TBB decide how many threads to use */
DEBUG_INIT ("cholmod start", Common) ;
return (TRUE) ;
}
/* ========================================================================== */
/* === cholmod_defaults ===================================================== */
/* ========================================================================== */
/* Set Common default parameters, except for the function pointers.
*
* workspace: none
*/
int CHOLMOD(defaults)
(
cholmod_common *Common
)
{
Int i ;
RETURN_IF_NULL_COMMON (FALSE) ;
/* ---------------------------------------------------------------------- */
/* default control parameters */
/* ---------------------------------------------------------------------- */
Common->dbound = 0.0 ;
Common->grow0 = 1.2 ;
Common->grow1 = 1.2 ;
Common->grow2 = 5 ;
Common->maxrank = 8 ;
Common->final_asis = TRUE ;
Common->final_super = TRUE ;
Common->final_ll = FALSE ;
Common->final_pack = TRUE ;
Common->final_monotonic = TRUE ;
Common->final_resymbol = FALSE ;
/* use simplicial factorization if flop/nnz(L) < 40, supernodal otherwise */
Common->supernodal = CHOLMOD_AUTO ;
Common->supernodal_switch = 40 ;
Common->nrelax [0] = 4 ;
Common->nrelax [1] = 16 ;
Common->nrelax [2] = 48 ;
Common->zrelax [0] = 0.8 ;
Common->zrelax [1] = 0.1 ;
Common->zrelax [2] = 0.05 ;
Common->prefer_zomplex = FALSE ;
Common->prefer_upper = TRUE ;
Common->prefer_binary = FALSE ;
Common->quick_return_if_not_posdef = FALSE ;
/* METIS workarounds */
Common->metis_memory = 0.0 ; /* > 0 for memory guard (2 is reasonable) */
Common->metis_nswitch = 3000 ;
Common->metis_dswitch = 0.66 ;
Common->print = 3 ;
Common->precise = FALSE ;
/* ---------------------------------------------------------------------- */
/* default ordering methods */
/* ---------------------------------------------------------------------- */
/* Note that if the Partition module is not installed, the CHOLMOD_METIS
* and CHOLMOD_NESDIS methods will not be available. cholmod_analyze will
* report the CHOLMOD_NOT_INSTALLED error, and safely skip over them.
*/
#if (CHOLMOD_MAXMETHODS < 9)
#error "CHOLMOD_MAXMETHODS must be 9 or more (defined in cholmod_core.h)."
#endif
/* default strategy: try given, AMD, and then METIS if AMD reports high
* fill-in. NESDIS can be used instead, if Common->default_nesdis is TRUE.
*/
Common->nmethods = 0 ; /* use default strategy */
Common->default_nesdis = FALSE ; /* use METIS in default strategy */
Common->current = 0 ; /* current method being tried */
Common->selected = 0 ; /* the best method selected */
/* first, fill each method with default parameters */
for (i = 0 ; i <= CHOLMOD_MAXMETHODS ; i++)
{
/* CHOLMOD's default method is AMD for A or AA' */
Common->method [i].ordering = CHOLMOD_AMD ;
/* CHOLMOD nested dissection and minimum degree parameter */
Common->method [i].prune_dense = 10.0 ; /* dense row/col control */
/* min degree parameters (AMD, COLAMD, SYMAMD, CAMD, CCOLAMD, CSYMAMD)*/
Common->method [i].prune_dense2 = -1 ; /* COLAMD dense row control */
Common->method [i].aggressive = TRUE ; /* aggressive absorption */
Common->method [i].order_for_lu = FALSE ;/* order for Cholesky not LU */
/* CHOLMOD's nested dissection (METIS + constrained AMD) */
Common->method [i].nd_small = 200 ; /* small graphs aren't cut */
Common->method [i].nd_compress = TRUE ; /* compress graph & subgraphs */
Common->method [i].nd_camd = 1 ; /* use CAMD */
Common->method [i].nd_components = FALSE ; /* lump connected comp. */
Common->method [i].nd_oksep = 1.0 ; /* sep ok if < oksep*n */
/* statistics for each method are not yet computed */
Common->method [i].fl = EMPTY ;
Common->method [i].lnz = EMPTY ;
}
Common->postorder = TRUE ; /* follow ordering with weighted postorder */
/* Next, define some methods. The first five use default parameters. */
Common->method [0].ordering = CHOLMOD_GIVEN ; /* skip if UserPerm NULL */
Common->method [1].ordering = CHOLMOD_AMD ;
Common->method [2].ordering = CHOLMOD_METIS ;
Common->method [3].ordering = CHOLMOD_NESDIS ;
Common->method [4].ordering = CHOLMOD_NATURAL ;
/* CHOLMOD's nested dissection with large leaves of separator tree */
Common->method [5].ordering = CHOLMOD_NESDIS ;
Common->method [5].nd_small = 20000 ;
/* CHOLMOD's nested dissection with tiny leaves, and no AMD ordering */
Common->method [6].ordering = CHOLMOD_NESDIS ;
Common->method [6].nd_small = 4 ;
Common->method [6].nd_camd = 0 ; /* no CSYMAMD or CAMD */
/* CHOLMOD's nested dissection with no dense node removal */
Common->method [7].ordering = CHOLMOD_NESDIS ;
Common->method [7].prune_dense = -1. ;
/* COLAMD for A*A', AMD for A */
Common->method [8].ordering = CHOLMOD_COLAMD ;
return (TRUE) ;
}
/* ========================================================================== */
/* === cholmod_finish ======================================================= */
/* ========================================================================== */
/* The last call to CHOLMOD must be cholmod_finish. You may call this routine
* more than once, and can safely call any other CHOLMOD routine after calling
* it (including cholmod_start).
*
* The statistics and parameter settings in Common are preserved. The
* workspace in Common is freed. This routine is just another name for
* cholmod_free_work.
*/
int CHOLMOD(finish)
(
cholmod_common *Common
)
{
return (CHOLMOD(free_work) (Common)) ;
}
/* ========================================================================== */
/* === cholmod_allocate_work ================================================ */
/* ========================================================================== */
/* Allocate and initialize workspace for CHOLMOD routines, or increase the size
* of already-allocated workspace. If enough workspace is already allocated,
* then nothing happens.
*
* workspace: Flag (nrow), Head (nrow+1), Iwork (iworksize), Xwork (xworksize)
*/
int CHOLMOD(allocate_work)
(
/* ---- input ---- */
size_t nrow, /* # of rows in the matrix A */
size_t iworksize, /* size of Iwork */
size_t xworksize, /* size of Xwork */
/* --------------- */
cholmod_common *Common
)
{
double *W ;
Int *Head ;
Int i ;
size_t nrow1 ;
int ok = TRUE ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (FALSE) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* Allocate Flag (nrow) and Head (nrow+1) */
/* ---------------------------------------------------------------------- */
nrow = MAX (1, nrow) ;
/* nrow1 = nrow + 1 */
nrow1 = CHOLMOD(add_size_t) (nrow, 1, &ok) ;
if (!ok)
{
/* nrow+1 causes size_t overflow ; problem is too large */
Common->status = CHOLMOD_TOO_LARGE ;
CHOLMOD(free_work) (Common) ;
return (FALSE) ;
}
if (nrow > Common->nrow)
{
if (Common->no_workspace_reallocate)
{
/* CHOLMOD is not allowed to change the workspace here */
Common->status = CHOLMOD_INVALID ;
return (FALSE) ;
}
/* free the old workspace (if any) and allocate new space */
Common->Flag = CHOLMOD(free) (Common->nrow, sizeof (Int), Common->Flag,
Common) ;
Common->Head = CHOLMOD(free) (Common->nrow+1,sizeof (Int), Common->Head,
Common) ;
Common->Flag = CHOLMOD(malloc) (nrow, sizeof (Int), Common) ;
Common->Head = CHOLMOD(malloc) (nrow1, sizeof (Int), Common) ;
/* record the new size of Flag and Head */
Common->nrow = nrow ;
if (Common->status < CHOLMOD_OK)
{
CHOLMOD(free_work) (Common) ;
return (FALSE) ;
}
/* initialize Flag and Head */
Common->mark = EMPTY ;
CHOLMOD(clear_flag) (Common) ;
Head = Common->Head ;
for (i = 0 ; i <= (Int) (nrow) ; i++)
{
Head [i] = EMPTY ;
}
}
/* ---------------------------------------------------------------------- */
/* Allocate Iwork (iworksize) */
/* ---------------------------------------------------------------------- */
iworksize = MAX (1, iworksize) ;
if (iworksize > Common->iworksize)
{
if (Common->no_workspace_reallocate)
{
/* CHOLMOD is not allowed to change the workspace here */
Common->status = CHOLMOD_INVALID ;
return (FALSE) ;
}
/* free the old workspace (if any) and allocate new space.
* integer overflow safely detected in cholmod_malloc */
CHOLMOD(free) (Common->iworksize, sizeof (Int), Common->Iwork, Common) ;
Common->Iwork = CHOLMOD(malloc) (iworksize, sizeof (Int), Common) ;
/* record the new size of Iwork */
Common->iworksize = iworksize ;
if (Common->status < CHOLMOD_OK)
{
CHOLMOD(free_work) (Common) ;
return (FALSE) ;
}
/* note that Iwork does not need to be initialized */
}
/* ---------------------------------------------------------------------- */
/* Allocate Xwork (xworksize) and set it to ((double) 0.) */
/* ---------------------------------------------------------------------- */
/* make sure xworksize is >= 1 */
xworksize = MAX (1, xworksize) ;
if (xworksize > Common->xworksize)
{
if (Common->no_workspace_reallocate)
{
/* CHOLMOD is not allowed to change the workspace here */
Common->status = CHOLMOD_INVALID ;
return (FALSE) ;
}
/* free the old workspace (if any) and allocate new space */
CHOLMOD(free) (Common->xworksize, sizeof (double), Common->Xwork,
Common) ;
Common->Xwork = CHOLMOD(malloc) (xworksize, sizeof (double), Common) ;
/* record the new size of Xwork */
Common->xworksize = xworksize ;
if (Common->status < CHOLMOD_OK)
{
CHOLMOD(free_work) (Common) ;
return (FALSE) ;
}
/* initialize Xwork */
W = Common->Xwork ;
for (i = 0 ; i < (Int) xworksize ; i++)
{
W [i] = 0. ;
}
}
return (TRUE) ;
}
/* ========================================================================== */
/* === cholmod_free_work ==================================================== */
/* ========================================================================== */
/* Deallocate the CHOLMOD workspace.
*
* workspace: deallocates all workspace in Common
*/
int CHOLMOD(free_work)
(
cholmod_common *Common
)
{
RETURN_IF_NULL_COMMON (FALSE) ;
Common->Flag = CHOLMOD(free) (Common->nrow, sizeof (Int),
Common->Flag, Common) ;
Common->Head = CHOLMOD(free) (Common->nrow+1, sizeof (Int),
Common->Head, Common) ;
Common->Iwork = CHOLMOD(free) (Common->iworksize, sizeof (Int),
Common->Iwork, Common) ;
Common->Xwork = CHOLMOD(free) (Common->xworksize, sizeof (double),
Common->Xwork, Common) ;
Common->nrow = 0 ;
Common->iworksize = 0 ;
Common->xworksize = 0 ;
return (TRUE) ;
}
/* ========================================================================== */
/* === cholmod_clear_flag =================================================== */
/* ========================================================================== */
/* Increment mark to ensure Flag [0..nrow-1] < mark. If integer overflow
* occurs, or mark was initially negative, reset the entire array. This is
* not an error condition, but an intended function of the Flag workspace.
*
* workspace: Flag (nrow). Does not modify Flag if nrow is zero.
*/
UF_long cholmod_l_clear_flag
(
cholmod_common *Common
)
{
Int i, nrow, *Flag ;
RETURN_IF_NULL_COMMON (-1) ;
Common->mark++ ;
if (Common->mark <= 0)
{
nrow = Common->nrow ;
Flag = Common->Flag ;
PRINT2 (("reset Flag: nrow "ID"\n", nrow)) ;
PRINT2 (("reset Flag: mark %ld\n", Common->mark)) ;
for (i = 0 ; i < nrow ; i++)
{
Flag [i] = EMPTY ;
}
Common->mark = 0 ;
}
return (Common->mark) ;
}
/* ========================================================================== */
/* ==== cholmod_maxrank ===================================================== */
/* ========================================================================== */
/* Find a valid value of Common->maxrank. Returns 0 if error, or 2, 4, or 8
* if successful. */
size_t cholmod_l_maxrank /* returns validated value of Common->maxrank */
(
/* ---- input ---- */
size_t n, /* A and L will have n rows */
/* --------------- */
cholmod_common *Common
)
{
size_t maxrank ;
RETURN_IF_NULL_COMMON (0) ;
maxrank = Common->maxrank ;
if (n > 0)
{
/* Ensure maxrank*n*sizeof(double) does not result in integer overflow.
* If n is so large that 2*n*sizeof(double) results in integer overflow
* (n = 268,435,455 if an Int is 32 bits), then maxrank will be 0 or 1,
* but maxrank will be set to 2 below. 2*n will not result in integer
* overflow, and CHOLMOD will run out of memory or safely detect integer
* overflow elsewhere.
*/
maxrank = MIN (maxrank, Size_max / (n * sizeof (double))) ;
}
if (maxrank <= 2)
{
maxrank = 2 ;
}
else if (maxrank <= 4)
{
maxrank = 4 ;
}
else
{
maxrank = 8 ;
}
return (maxrank) ;
}
/* ========================================================================== */
/* === cholmod_dbound ======================================================= */
/* ========================================================================== */
/* Ensure the absolute value of a diagonal entry, D (j,j), is greater than
* Common->dbound. This routine is not meant for the user to call. It is used
* by the various LDL' factorization and update/downdate routines. The
* default value of Common->dbound is zero, and in that case this routine is not
* called at all. No change is made if D (j,j) is NaN. CHOLMOD does not call
* this routine if Common->dbound is NaN.
*/
double cholmod_l_dbound /* returns modified diagonal entry of D */
(
/* ---- input ---- */
double dj, /* diagonal entry of D, for LDL' factorization */
/* --------------- */
cholmod_common *Common
)
{
double dbound ;
RETURN_IF_NULL_COMMON (0) ;
if (!IS_NAN (dj))
{
dbound = Common->dbound ;
if (dj < 0)
{
if (dj > -dbound)
{
dj = -dbound ;
Common->ndbounds_hit++ ;
if (Common->status == CHOLMOD_OK)
{
ERROR (CHOLMOD_DSMALL, "diagonal below threshold") ;
}
}
}
else
{
if (dj < dbound)
{
dj = dbound ;
Common->ndbounds_hit++ ;
if (Common->status == CHOLMOD_OK)
{
ERROR (CHOLMOD_DSMALL, "diagonal below threshold") ;
}
}
}
}
return (dj) ;
}

0
spqr_mini/cholmod_common.h
View File

2290
spqr_mini/cholmod_core.h
View File

File diff suppressed because it is too large Load Diff

80
spqr_mini/cholmod_error.c
View File

@ -1,80 +0,0 @@
/* ========================================================================== */
/* === Core/cholmod_error =================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Core Module. Copyright (C) 2005-2006,
* Univ. of Florida. Author: Timothy A. Davis
* The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* CHOLMOD error-handling routine. */
#include "cholmod_internal.h"
#include "cholmod_core.h"
/* ========================================================================== */
/* ==== cholmod_error ======================================================= */
/* ========================================================================== */
/* An error has occurred. Set the status, optionally print an error message,
* and call the user error-handling routine (if it exists). If
* Common->try_catch is TRUE, then CHOLMOD is inside a try/catch block.
* The status is set, but no message is printed and the user error handler
* is not called. This is not (yet) an error, since CHOLMOD may recover.
*
* In the current version, this try/catch mechanism is used internally only in
* cholmod_analyze, which tries multiple ordering methods and picks the best
* one. If one or more ordering method fails, it keeps going. Only one
* ordering needs to succeed for cholmod_analyze to succeed.
*/
int CHOLMOD(error)
(
/* ---- input ---- */
int status, /* error status */
const char *file, /* name of source code file where error occured */
int line, /* line number in source code file where error occured*/
const char *message, /* error message */
/* --------------- */
cholmod_common *Common
)
{
RETURN_IF_NULL_COMMON (FALSE) ;
Common->status = status ;
if (!(Common->try_catch))
{
#ifndef NPRINT
/* print a warning or error message */
if (Common->print_function != NULL)
{
if (status > 0 && Common->print > 1)
{
(Common->print_function) ("CHOLMOD warning: %s\n", message) ;
fflush (stdout) ;
fflush (stderr) ;
}
else if (Common->print > 0)
{
(Common->print_function) ("CHOLMOD error: %s\n", message) ;
fflush (stdout) ;
fflush (stderr) ;
}
}
#endif
/* call the user error handler, if it exists */
if (Common->error_handler != NULL)
{
Common->error_handler (status, file, line, message) ;
}
}
return (TRUE) ;
}

400
spqr_mini/cholmod_internal.h
View File

@ -1,400 +0,0 @@
/* ========================================================================== */
/* === Include/cholmod_internal.h =========================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Include/cholmod_internal.h.
* Copyright (C) 2005-2006, Univ. of Florida. Author: Timothy A. Davis
* CHOLMOD/Include/cholmod_internal.h is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* CHOLMOD internal include file.
*
* This file contains internal definitions for CHOLMOD, not meant to be included
* in user code. They define macros that are not prefixed with CHOLMOD_. This
* file can safely #include'd in user code if you want to make use of the
* macros defined here, and don't mind the possible name conflicts with your
* code, however.
*
* Required by all CHOLMOD routines. Not required by any user routine that
* uses CHOLMOMD. Unless debugging is enabled, this file does not require any
* CHOLMOD module (not even the Core module).
*
* If debugging is enabled, all CHOLMOD modules require the Check module.
* Enabling debugging requires that this file be editted. Debugging cannot be
* enabled with a compiler flag. This is because CHOLMOD is exceedingly slow
* when debugging is enabled. Debugging is meant for development of CHOLMOD
* itself, not by users of CHOLMOD.
*/
#ifndef CHOLMOD_INTERNAL_H
#define CHOLMOD_INTERNAL_H
/* ========================================================================== */
/* === large file I/O ======================================================= */
/* ========================================================================== */
/* Definitions for large file I/O must come before any other #includes. If
* this causes problems (may not be portable to all platforms), then compile
* CHOLMOD with -DNLARGEFILE. You must do this for MATLAB 6.5 and earlier,
* for example. */
//#include "cholmod_io64.h"
/* ========================================================================== */
/* === debugging and basic includes ========================================= */
/* ========================================================================== */
/* turn off debugging */
#ifndef SPQR_NDEBUG
#define SPQR_NDEBUG
#endif
/* Uncomment this line to enable debugging. CHOLMOD will be very slow.
#undef SPQR_NDEBUG
*/
#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif
#if !defined(NPRINT) || !defined(SPQR_NDEBUG)
#include <stdio.h>
#endif
#include <stddef.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <stdlib.h>
/* ========================================================================== */
/* === basic definitions ==================================================== */
/* ========================================================================== */
/* Some non-conforming compilers insist on defining TRUE and FALSE. */
#undef TRUE
#undef FALSE
#define TRUE 1
#define FALSE 0
#define BOOLEAN(x) ((x) ? TRUE : FALSE)
/* NULL should already be defined, but ensure it is here. */
#ifndef NULL
#define NULL ((void *) 0)
#endif
/* FLIP is a "negation about -1", and is used to mark an integer i that is
* normally non-negative. FLIP (EMPTY) is EMPTY. FLIP of a number > EMPTY
* is negative, and FLIP of a number < EMTPY is positive. FLIP (FLIP (i)) = i
* for all integers i. UNFLIP (i) is >= EMPTY. */
#define EMPTY (-1)
#define FLIP(i) (-(i)-2)
#define UNFLIP(i) (((i) < EMPTY) ? FLIP (i) : (i))
/* MAX and MIN are not safe to use for NaN's */
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define MAX3(a,b,c) (((a) > (b)) ? (MAX (a,c)) : (MAX (b,c)))
#define MAX4(a,b,c,d) (((a) > (b)) ? (MAX3 (a,c,d)) : (MAX3 (b,c,d)))
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define IMPLIES(p,q) (!(p) || (q))
/* find the sign: -1 if x < 0, 1 if x > 0, zero otherwise.
* Not safe for NaN's */
#define SIGN(x) (((x) < 0) ? (-1) : (((x) > 0) ? 1 : 0))
/* round up an integer x to a multiple of s */
#define ROUNDUP(x,s) ((s) * (((x) + ((s) - 1)) / (s)))
#define ERROR(status,msg) \
CHOLMOD(error) (status, __FILE__, __LINE__, msg, Common)
/* Check a pointer and return if null. Set status to invalid, unless the
* status is already "out of memory" */
#define RETURN_IF_NULL(A,result) \
{ \
if ((A) == NULL) \
{ \
if (Common->status != CHOLMOD_OUT_OF_MEMORY) \
{ \
ERROR (CHOLMOD_INVALID, "argument missing") ; \
} \
return (result) ; \
} \
}
/* Return if Common is NULL or invalid */
#define RETURN_IF_NULL_COMMON(result) \
{ \
if (Common == NULL) \
{ \
return (result) ; \
} \
if (Common->itype != ITYPE || Common->dtype != DTYPE) \
{ \
Common->status = CHOLMOD_INVALID ; \
return (result) ; \
} \
}
#define IS_NAN(x) CHOLMOD_IS_NAN(x)
#define IS_ZERO(x) CHOLMOD_IS_ZERO(x)
#define IS_NONZERO(x) CHOLMOD_IS_NONZERO(x)
#define IS_LT_ZERO(x) CHOLMOD_IS_LT_ZERO(x)
#define IS_GT_ZERO(x) CHOLMOD_IS_GT_ZERO(x)
#define IS_LE_ZERO(x) CHOLMOD_IS_LE_ZERO(x)
/* 1e308 is a huge number that doesn't take many characters to print in a
* file, in CHOLMOD/Check/cholmod_read and _write. Numbers larger than this
* are interpretted as Inf, since sscanf doesn't read in Inf's properly.
* This assumes IEEE double precision arithmetic. DBL_MAX would be a little
* better, except that it takes too many digits to print in a file. */
#define HUGE_DOUBLE 1e308
/* ========================================================================== */
/* === int/UF_long and double/float definitions ============================= */
/* ========================================================================== */
/* CHOLMOD is designed for 3 types of integer variables:
*
* (1) all integers are int
* (2) most integers are int, some are UF_long
* (3) all integers are UF_long
*
* and two kinds of floating-point values:
*
* (1) double
* (2) float
*
* the complex types (ANSI-compatible complex, and MATLAB-compatable zomplex)
* are based on the double or float type, and are not selected here. They
* are typically selected via template routines.
*
* This gives 6 different modes in which CHOLMOD can be compiled (only the
* first two are currently supported):
*
* DINT double, int prefix: cholmod_
* DLONG double, UF_long prefix: cholmod_l_
* DMIX double, mixed int/UF_long prefix: cholmod_m_
* SINT float, int prefix: cholmod_si_
* SLONG float, UF_long prefix: cholmod_sl_
* SMIX float, mixed int/log prefix: cholmod_sm_
*
* These are selected with compile time flags (-DDLONG, for example). If no
* flag is selected, the default is DINT.
*
* All six versions use the same include files. The user-visible include files
* are completely independent of which int/UF_long/double/float version is being
* used. The integer / real types in all data structures (sparse, triplet,
* dense, common, and triplet) are defined at run-time, not compile-time, so
* there is only one "cholmod_sparse" data type. Void pointers are used inside
* that data structure to point to arrays of the proper type. Each data
* structure has an itype and dtype field which determines the kind of basic
* types used. These are defined in Include/cholmod_core.h.
*
* FUTURE WORK: support all six types (float, and mixed int/UF_long)
*
* UF_long is normally defined as long. However, for WIN64 it is __int64.
* It can also be redefined for other platforms, by modifying UFconfig.h.
*/
#include "UFconfig.h"
/* -------------------------------------------------------------------------- */
/* Size_max: the largest value of size_t */
/* -------------------------------------------------------------------------- */
#define Size_max ((size_t) (-1))
/* routines for doing arithmetic on size_t, and checking for overflow */
size_t cholmod_add_size_t (size_t a, size_t b, int *ok) ;
size_t cholmod_mult_size_t (size_t a, size_t k, int *ok) ;
size_t cholmod_l_add_size_t (size_t a, size_t b, int *ok) ;
size_t cholmod_l_mult_size_t (size_t a, size_t k, int *ok) ;
/* -------------------------------------------------------------------------- */
/* double (also complex double), UF_long */
/* -------------------------------------------------------------------------- */
#ifdef DLONG
#define Real double
#define Int UF_long
#define Int_max UF_long_max
#define CHOLMOD(name) cholmod_l_ ## name
#define LONG
#define DOUBLE
#define ITYPE CHOLMOD_LONG
#define DTYPE CHOLMOD_DOUBLE
#define ID UF_long_id
/* -------------------------------------------------------------------------- */
/* double, int/UF_long */
/* -------------------------------------------------------------------------- */
#elif defined (DMIX)
#error "mixed int/UF_long not yet supported"
/* -------------------------------------------------------------------------- */
/* single, int */
/* -------------------------------------------------------------------------- */
#elif defined (SINT)
#error "single-precision not yet supported"
/* -------------------------------------------------------------------------- */
/* single, UF_long */
/* -------------------------------------------------------------------------- */
#elif defined (SLONG)
#error "single-precision not yet supported"
/* -------------------------------------------------------------------------- */
/* single, int/UF_long */
/* -------------------------------------------------------------------------- */
#elif defined (SMIX)
#error "single-precision not yet supported"
/* -------------------------------------------------------------------------- */
/* double (also complex double), int: this is the default */
/* -------------------------------------------------------------------------- */
#else
#ifndef DINT
#define DINT
#endif
#define INT
#define DOUBLE
#define Real double
#define Int int
#define Int_max INT_MAX
#define CHOLMOD(name) cholmod_ ## name
#define ITYPE CHOLMOD_INT
#define DTYPE CHOLMOD_DOUBLE
#define ID "%d"
#endif
/* ========================================================================== */
/* === real/complex arithmetic ============================================== */
/* ========================================================================== */
//#include "cholmod_complexity.h"
/* ========================================================================== */
/* === Architecture and BLAS ================================================ */
/* ========================================================================== */
#define BLAS_OK Common->blas_ok
#include "cholmod_blas.h"
/* ========================================================================== */
/* === debugging definitions ================================================ */
/* ========================================================================== */
#ifndef SPQR_NDEBUG
#include <assert.h>
#include "cholmod.h"
/* The cholmod_dump routines are in the Check module. No CHOLMOD routine
* calls the cholmod_check_* or cholmod_print_* routines in the Check module,
* since they use Common workspace that may already be in use. Instead, they
* use the cholmod_dump_* routines defined there, which allocate their own
* workspace if they need it. */
#ifndef EXTERN
#define EXTERN extern
#endif
/* double, int */
EXTERN int cholmod_dump ;
EXTERN int cholmod_dump_malloc ;
UF_long cholmod_dump_sparse (cholmod_sparse *, const char *, cholmod_common *);
int cholmod_dump_factor (cholmod_factor *, const char *, cholmod_common *) ;
int cholmod_dump_triplet (cholmod_triplet *, const char *, cholmod_common *) ;
int cholmod_dump_dense (cholmod_dense *, const char *, cholmod_common *) ;
int cholmod_dump_subset (int *, size_t, size_t, const char *,
cholmod_common *) ;
int cholmod_dump_perm (int *, size_t, size_t, const char *, cholmod_common *) ;
int cholmod_dump_parent (int *, size_t, const char *, cholmod_common *) ;
void cholmod_dump_init (const char *, cholmod_common *) ;
int cholmod_dump_mem (const char *, UF_long, cholmod_common *) ;
void cholmod_dump_real (const char *, Real *, UF_long, UF_long, int, int,
cholmod_common *) ;
void cholmod_dump_super (UF_long, int *, int *, int *, int *, double *, int,
cholmod_common *) ;
int cholmod_dump_partition (UF_long, int *, int *, int *, int *, UF_long,
cholmod_common *) ;
int cholmod_dump_work(int, int, UF_long, cholmod_common *) ;
/* double, UF_long */
EXTERN int cholmod_l_dump ;
EXTERN int cholmod_l_dump_malloc ;
UF_long cholmod_l_dump_sparse (cholmod_sparse *, const char *,
cholmod_common *) ;
int cholmod_l_dump_factor (cholmod_factor *, const char *, cholmod_common *) ;
int cholmod_l_dump_triplet (cholmod_triplet *, const char *, cholmod_common *);
int cholmod_l_dump_dense (cholmod_dense *, const char *, cholmod_common *) ;
int cholmod_l_dump_subset (UF_long *, size_t, size_t, const char *,
cholmod_common *) ;
int cholmod_l_dump_perm (UF_long *, size_t, size_t, const char *,
cholmod_common *) ;
int cholmod_l_dump_parent (UF_long *, size_t, const char *, cholmod_common *) ;
void cholmod_l_dump_init (const char *, cholmod_common *) ;
int cholmod_l_dump_mem (const char *, UF_long, cholmod_common *) ;
void cholmod_l_dump_real (const char *, Real *, UF_long, UF_long, int, int,
cholmod_common *) ;
void cholmod_l_dump_super (UF_long, UF_long *, UF_long *, UF_long *, UF_long *,
double *, int, cholmod_common *) ;
int cholmod_l_dump_partition (UF_long, UF_long *, UF_long *, UF_long *,
UF_long *, UF_long, cholmod_common *) ;
int cholmod_l_dump_work(int, int, UF_long, cholmod_common *) ;
#define DEBUG_INIT(s,Common) { CHOLMOD(dump_init)(s, Common) ; }
#define ASSERT(expression) (assert (expression))
#define PRK(k,params) \
{ \
if (CHOLMOD(dump) >= (k) && Common->print_function != NULL) \
{ \
(Common->print_function) params ; \
} \
}
#define PRINT0(params) PRK (0, params)
#define PRINT1(params) PRK (1, params)
#define PRINT2(params) PRK (2, params)
#define PRINT3(params) PRK (3, params)
#define PRINTM(params) \
{ \
if (CHOLMOD(dump_malloc) > 0) \
{ \
printf params ; \
} \
}
#define DEBUG(statement) statement
#else
/* Debugging disabled (the normal case) */
#define PRK(k,params)
#define DEBUG_INIT(s,Common)
#define PRINT0(params)
#define PRINT1(params)
#define PRINT2(params)
#define PRINT3(params)
#define PRINTM(params)
#define ASSERT(expression)
#define DEBUG(statement)
#endif
#endif

563
spqr_mini/cholmod_memory.c
View File

@ -1,563 +0,0 @@
/* ========================================================================== */
/* === Core/cholmod_memory ================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Core Module. Copyright (C) 2005-2006,
* Univ. of Florida. Author: Timothy A. Davis
* The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Core memory management routines:
*
* Primary routines:
* -----------------
* cholmod_malloc malloc wrapper
* cholmod_free free wrapper
*
* Secondary routines:
* -------------------
* cholmod_calloc calloc wrapper
* cholmod_realloc realloc wrapper
* cholmod_realloc_multiple realloc wrapper for multiple objects
*
* The user may make use of these, just like malloc and free. You can even
* malloc an object and safely free it with cholmod_free, and visa versa
* (except that the memory usage statistics will be corrupted). These routines
* do differ from malloc and free. If cholmod_free is given a NULL pointer,
* for example, it does nothing (unlike the ANSI free). cholmod_realloc does
* not return NULL if given a non-NULL pointer and a nonzero size, even if it
* fails (it sets an error code in Common->status instead).
*
* CHOLMOD keeps track of the amount of memory it has allocated, and so the
* cholmod_free routine includes as a parameter the size of the object being
* freed. This is only used for memory usage statistics, which are very useful
* in finding memory leaks in your program. If you, the user of CHOLMOD, pass
* the wrong size, the only consequence is that the memory usage statistics
* will be invalid. This will causes assertions to fail if CHOLMOD is
* compiled with debugging enabled, but otherwise it will cause no errors.
*
* The cholmod_free_* routines for each CHOLMOD object keep track of the size
* of the blocks they free, so they do not require you to pass their sizes
* as a parameter.
*
* If a block of size zero is requested, these routines allocate a block of
* size one instead.
*/
#include "cholmod_internal.h"
#include "cholmod_core.h"
/* ========================================================================== */
/* === cholmod_add_size_t =================================================== */
/* ========================================================================== */
/* Safely compute a+b, and check for integer overflow. If overflow occurs,
* return 0 and set OK to FALSE. Also return 0 if OK is FALSE on input. */
size_t CHOLMOD(add_size_t) (size_t a, size_t b, int *ok)
{
size_t s = a + b ;
(*ok) = (*ok) && (s >= a) ;
return ((*ok) ? s : 0) ;
}
/* ========================================================================== */
/* === cholmod_mult_size_t ================================================== */
/* ========================================================================== */
/* Safely compute a*k, where k should be small, and check for integer overflow.
* If overflow occurs, return 0 and set OK to FALSE. Also return 0 if OK is
* FALSE on input. */
size_t CHOLMOD(mult_size_t) (size_t a, size_t k, int *ok)
{
size_t p = 0, s ;
while (*ok)
{
if (k % 2)
{
p = p + a ;
(*ok) = (*ok) && (p >= a) ;
}
k = k / 2 ;
if (!k) return (p) ;
s = a + a ;
(*ok) = (*ok) && (s >= a) ;
a = s ;
}
return (0) ;
}
/* ========================================================================== */
/* === cholmod_malloc ======================================================= */
/* ========================================================================== */
/* Wrapper around malloc routine. Allocates space of size MAX(1,n)*size, where
* size is normally a sizeof (...).
*
* This routine, cholmod_calloc, and cholmod_realloc do not set Common->status
* to CHOLMOD_OK on success, so that a sequence of cholmod_malloc's, _calloc's,
* or _realloc's can be used. If any of them fails, the Common->status will
* hold the most recent error status.
*
* Usage, for a pointer to int:
*
* p = cholmod_malloc (n, sizeof (int), Common)
*
* Uses a pointer to the malloc routine (or its equivalent) defined in Common.
*/
void *CHOLMOD(malloc) /* returns pointer to the newly malloc'd block */
(
/* ---- input ---- */
size_t n, /* number of items */
size_t size, /* size of each item */
/* --------------- */
cholmod_common *Common
)
{
void *p ;
size_t s ;
int ok = TRUE ;
RETURN_IF_NULL_COMMON (NULL) ;
if (size == 0)
{
ERROR (CHOLMOD_INVALID, "sizeof(item) must be > 0") ;
p = NULL ;
}
else if (n >= (Size_max / size) || n >= Int_max)
{
/* object is too big to allocate without causing integer overflow */
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
p = NULL ;
}
else
{
/* call malloc, or its equivalent */
s = CHOLMOD(mult_size_t) (MAX (1,n), size, &ok) ;
p = ok ? ((Common->malloc_memory) (s)) : NULL ;
if (p == NULL)
{
/* failure: out of memory */
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
}
else
{
/* success: increment the count of objects allocated */
Common->malloc_count++ ;
Common->memory_inuse += (n * size) ;
Common->memory_usage =
MAX (Common->memory_usage, Common->memory_inuse) ;
PRINTM (("cholmod_malloc %p %d cnt: %d inuse %d\n",
p, n*size, Common->malloc_count, Common->memory_inuse)) ;
}
}
return (p) ;
}
/* ========================================================================== */
/* === cholmod_free ========================================================= */
/* ========================================================================== */
/* Wrapper around free routine. Returns NULL, which can be assigned to the
* pointer being freed, as in:
*
* p = cholmod_free (n, sizeof (int), p, Common) ;
*
* In CHOLMOD, the syntax:
*
* cholmod_free (n, sizeof (int), p, Common) ;
*
* is used if p is a local pointer and the routine is returning shortly.
* Uses a pointer to the free routine (or its equivalent) defined in Common.
* Nothing is freed if the pointer is NULL.
*/
void *CHOLMOD(free) /* always returns NULL */
(
/* ---- input ---- */
size_t n, /* number of items */
size_t size, /* size of each item */
/* ---- in/out --- */
void *p, /* block of memory to free */
/* --------------- */
cholmod_common *Common
)
{
RETURN_IF_NULL_COMMON (NULL) ;
if (p != NULL)
{
/* only free the object if the pointer is not NULL */
/* call free, or its equivalent */
(Common->free_memory) (p) ;
Common->malloc_count-- ;
Common->memory_inuse -= (n * size) ;
PRINTM (("cholmod_free %p %d cnt: %d inuse %d\n",
p, n*size, Common->malloc_count, Common->memory_inuse)) ;
/* This assertion will fail if the user calls cholmod_malloc and
* cholmod_free with mismatched memory sizes. It shouldn't fail
* otherwise. */
ASSERT (IMPLIES (Common->malloc_count == 0, Common->memory_inuse == 0));
}
/* return NULL, and the caller should assign this to p. This avoids
* freeing the same pointer twice. */
return (NULL) ;
}
/* ========================================================================== */
/* === cholmod_calloc ======================================================= */
/* ========================================================================== */
/* Wrapper around calloc routine.
*
* Uses a pointer to the calloc routine (or its equivalent) defined in Common.
* This routine is identical to malloc, except that it zeros the newly allocated
* block to zero.
*/
void *CHOLMOD(calloc) /* returns pointer to the newly calloc'd block */
(
/* ---- input ---- */
size_t n, /* number of items */
size_t size, /* size of each item */
/* --------------- */
cholmod_common *Common
)
{
void *p ;
RETURN_IF_NULL_COMMON (NULL) ;
if (size == 0)
{
ERROR (CHOLMOD_INVALID, "sizeof(item) must be > 0") ;
p = NULL ;
}
else if (n >= (Size_max / size) || n >= Int_max)
{
/* object is too big to allocate without causing integer overflow */
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
p = NULL ;
}
else
{
/* call calloc, or its equivalent */
p = (Common->calloc_memory) (MAX (1,n), size) ;
if (p == NULL)
{
/* failure: out of memory */
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
}
else
{
/* success: increment the count of objects allocated */
Common->malloc_count++ ;
Common->memory_inuse += (n * size) ;
Common->memory_usage =
MAX (Common->memory_usage, Common->memory_inuse) ;
PRINTM (("cholmod_malloc %p %d cnt: %d inuse %d\n",
p, n*size, Common->malloc_count, Common->memory_inuse)) ;
}
}
return (p) ;
}
/* ========================================================================== */
/* === cholmod_realloc ====================================================== */
/* ========================================================================== */
/* Wrapper around realloc routine. Given a pointer p to a block of size
* (*n)*size memory, it changes the size of the block pointed to by p to be
* MAX(1,nnew)*size in size. It may return a pointer different than p. This
* should be used as (for a pointer to int):
*
* p = cholmod_realloc (nnew, sizeof (int), p, *n, Common) ;
*
* If p is NULL, this is the same as p = cholmod_malloc (...).
* A size of nnew=0 is treated as nnew=1.
*
* If the realloc fails, p is returned unchanged and Common->status is set
* to CHOLMOD_OUT_OF_MEMORY. If successful, Common->status is not modified,
* and p is returned (possibly changed) and pointing to a large block of memory.
*
* Uses a pointer to the realloc routine (or its equivalent) defined in Common.
*/
void *CHOLMOD(realloc) /* returns pointer to reallocated block */
(
/* ---- input ---- */
size_t nnew, /* requested # of items in reallocated block */
size_t size, /* size of each item */
/* ---- in/out --- */
void *p, /* block of memory to realloc */
size_t *n, /* current size on input, nnew on output if successful*/
/* --------------- */
cholmod_common *Common
)
{
size_t nold = (*n) ;
void *pnew ;
size_t s ;
int ok = TRUE ;
RETURN_IF_NULL_COMMON (NULL) ;
if (size == 0)
{
ERROR (CHOLMOD_INVALID, "sizeof(item) must be > 0") ;
p = NULL ;
}
else if (p == NULL)
{
/* A fresh object is being allocated. */
PRINT1 (("realloc fresh: %d %d\n", nnew, size)) ;
p = CHOLMOD(malloc) (nnew, size, Common) ;
*n = (p == NULL) ? 0 : nnew ;
}
else if (nold == nnew)
{
/* Nothing to do. Do not change p or n. */
PRINT1 (("realloc nothing: %d %d\n", nnew, size)) ;
}
else if (nnew >= (Size_max / size) || nnew >= Int_max)
{
/* failure: nnew is too big. Do not change p or n. */
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
}
else
{
/* The object exists, and is changing to some other nonzero size. */
/* call realloc, or its equivalent */
PRINT1 (("realloc : %d to %d, %d\n", nold, nnew, size)) ;
pnew = NULL ;
s = CHOLMOD(mult_size_t) (MAX (1,nnew), size, &ok) ;
pnew = ok ? ((Common->realloc_memory) (p, s)) : NULL ;
if (pnew == NULL)
{
/* Do not change p, since it still points to allocated memory */
if (nnew <= nold)
{
/* The attempt to reduce the size of the block from n to
* nnew has failed. The current block is not modified, so
* pretend to succeed, but do not change p. Do change
* CHOLMOD's notion of the size of the block, however. */
*n = nnew ;
PRINTM (("nnew <= nold failed, pretend to succeed\n")) ;
PRINTM (("cholmod_free %p %d cnt: %d inuse %d\n"
"cholmod_malloc %p %d cnt: %d inuse %d\n",
p, nold*size, Common->malloc_count-1,
Common->memory_inuse - nold*size,
p, nnew*size, Common->malloc_count,
Common->memory_inuse + (nnew-nold)*size)) ;
Common->memory_inuse += ((nnew-nold) * size) ;
}
else
{
/* Increasing the size of the block has failed.
* Do not change n. */
ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory") ;
}
}
else
{
/* success: return revised p and change the size of the block */
PRINTM (("cholmod_free %p %d cnt: %d inuse %d\n"
"cholmod_malloc %p %d cnt: %d inuse %d\n",
p, nold*size, Common->malloc_count-1,
Common->memory_inuse - nold*size,
pnew, nnew*size, Common->malloc_count,
Common->memory_inuse + (nnew-nold)*size)) ;
p = pnew ;
*n = nnew ;
Common->memory_inuse += ((nnew-nold) * size) ;
}
Common->memory_usage = MAX (Common->memory_usage, Common->memory_inuse);
}
return (p) ;
}
/* ========================================================================== */
/* === cholmod_realloc_multiple ============================================= */
/* ========================================================================== */
/* reallocate multiple blocks of memory, all of the same size (up to two integer
* and two real blocks). Either reallocations all succeed, or all are returned
* in the original size (they are freed if the original size is zero). The nnew
* blocks are of size 1 or more.
*/
int CHOLMOD(realloc_multiple)
(
/* ---- input ---- */
size_t nnew, /* requested # of items in reallocated blocks */
int nint, /* number of int/UF_long blocks */
int xtype, /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
/* ---- in/out --- */
void **I, /* int or UF_long block */
void **J, /* int or UF_long block */
void **X, /* complex or double block */
void **Z, /* zomplex case only: double block */
size_t *nold_p, /* current size of the I,J,X,Z blocks on input,
* nnew on output if successful */
/* --------------- */
cholmod_common *Common
)
{
double *xx, *zz ;
size_t i, j, x, z, nold ;
RETURN_IF_NULL_COMMON (FALSE) ;
if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
{
ERROR (CHOLMOD_INVALID, "invalid xtype") ;
return (FALSE) ;
}
nold = *nold_p ;
if (nint < 1 && xtype == CHOLMOD_PATTERN)
{
/* nothing to do */
return (TRUE) ;
}
i = nold ;
j = nold ;
x = nold ;
z = nold ;
if (nint > 0)
{
*I = CHOLMOD(realloc) (nnew, sizeof (Int), *I, &i, Common) ;
}
if (nint > 1)
{
*J = CHOLMOD(realloc) (nnew, sizeof (Int), *J, &j, Common) ;
}
switch (xtype)
{
case CHOLMOD_REAL:
*X = CHOLMOD(realloc) (nnew, sizeof (double), *X, &x, Common) ;
break ;
case CHOLMOD_COMPLEX:
*X = CHOLMOD(realloc) (nnew, 2*sizeof (double), *X, &x, Common) ;
break ;
case CHOLMOD_ZOMPLEX:
*X = CHOLMOD(realloc) (nnew, sizeof (double), *X, &x, Common) ;
*Z = CHOLMOD(realloc) (nnew, sizeof (double), *Z, &z, Common) ;
break ;
}
if (Common->status < CHOLMOD_OK)
{
/* one or more realloc's failed. Resize all back down to nold. */
if (nold == 0)
{
if (nint > 0)
{
*I = CHOLMOD(free) (i, sizeof (Int), *I, Common) ;
}
if (nint > 1)
{
*J = CHOLMOD(free) (j, sizeof (Int), *J, Common) ;
}
switch (xtype)
{
case CHOLMOD_REAL:
*X = CHOLMOD(free) (x, sizeof (double), *X, Common) ;
break ;
case CHOLMOD_COMPLEX:
*X = CHOLMOD(free) (x, 2*sizeof (double), *X, Common) ;
break ;
case CHOLMOD_ZOMPLEX:
*X = CHOLMOD(free) (x, sizeof (double), *X, Common) ;
*Z = CHOLMOD(free) (x, sizeof (double), *Z, Common) ;
break ;
}
}
else
{
if (nint > 0)
{
*I = CHOLMOD(realloc) (nold, sizeof (Int), *I, &i, Common) ;
}
if (nint > 1)
{
*J = CHOLMOD(realloc) (nold, sizeof (Int), *J, &j, Common) ;
}
switch (xtype)
{
case CHOLMOD_REAL:
*X = CHOLMOD(realloc) (nold, sizeof (double), *X, &x,
Common) ;
break ;
case CHOLMOD_COMPLEX:
*X = CHOLMOD(realloc) (nold, 2*sizeof (double), *X, &x,
Common) ;
break ;
case CHOLMOD_ZOMPLEX:
*X = CHOLMOD(realloc) (nold, sizeof (double), *X, &x,
Common) ;
*Z = CHOLMOD(realloc) (nold, sizeof (double), *Z, &z,
Common) ;
break ;
}
}
return (FALSE) ;
}
if (nold == 0)
{
/* New space was allocated. Clear the first entry so that valgrind
* doesn't complain about its access in change_complexity
* (Core/cholmod_complex.c). */
xx = *X ;
zz = *Z ;
switch (xtype)
{
case CHOLMOD_REAL:
xx [0] = 0 ;
break ;
case CHOLMOD_COMPLEX:
xx [0] = 0 ;
xx [1] = 0 ;
break ;
case CHOLMOD_ZOMPLEX:
xx [0] = 0 ;
zz [0] = 0 ;
break ;
}
}
/* all realloc's succeeded, change size to reflect realloc'ed size. */
*nold_p = nnew ;
return (TRUE) ;
}

654
spqr_mini/spqr_front.cpp
View File

@ -1,654 +0,0 @@
// =============================================================================
// === spqr_front ==============================================================
// =============================================================================
/* Given an m-by-n frontal matrix, use Householder reflections to reduce it
to upper trapezoidal form. Columns 0:npiv-1 are checked against tol.
0 # x x x x x x
1 # x x x x x x
2 # x x x x x x
3 # x x x x x x
4 - # x x x x x <- Stair [0] = 4
5 # x x x x x
6 # x x x x x
7 # x x x x x
8 - # x x x x <- Stair [1] = 8
9 # x x x x
10 # x x x x
11 # x x x x
12 - # x x x <- Stair [2] = 12
13 # x x x
14 - # x x <- Stair [3] = 14
15 - # x <- Stair [4] = 15
16 - # <- Stair [5] = 16
- <- Stair [6] = 17
Suppose npiv = 3, and no columns fall below tol:
0 R r r r r r r
1 h R r r r r r
2 h h R r r r r
3 h h h C c c c
4 - h h h C c c <- Stair [0] = 4
5 h h h h C c
6 h h h h h C
7 h h h h h h
8 - h h h h h <- Stair [1] = 8
9 h h h h h
10 h h h h h
11 h h h h h
12 - h h h h <- Stair [2] = 12
13 h h h h
14 - h h h <- Stair [3] = 14
15 - h h <- Stair [4] = 15
16 - h <- Stair [5] = 16
- <- Stair [6] = 17
where r is an entry in the R block, c is an entry in the C (contribution)
block, and h is an entry in the H block (Householder vectors). Diagonal
entries are capitalized.
Suppose the 2nd column has a norm <= tol; the result is a "squeezed" R:
0 R r r r r r r <- Stair [1] = 0 to denote a dead pivot column
1 h - R r r r r
2 h - h C c c c
3 h - h h C c c
4 - - h h h C c <- Stair [0] = 4
5 - h h h h C
6 - h h h h h
7 - h h h h h
8 - h h h h h
9 h h h h h
10 h h h h h
11 h h h h h
12 - h h h h <- Stair [2] = 12
13 h h h h
14 - h h h <- Stair [3] = 14
15 - h h <- Stair [4] = 15
16 - h <- Stair [5] = 16
- <- Stair [6] = 17
where "diagonal" entries are capitalized. The 2nd H vector is missing
(it is H2 = I, identity). The 2nd column of R is not deleted. The
contribution block is always triangular, but the first npiv columns of
the R can become "staggered". Columns npiv to n-1 in the R block are
always the same length.
If columns are found "dead", the staircase may be updated. Stair[k] is
set to zero if k is dead. Also, Stair[k] is increased, if necessary, to
ensure that R and C reside within the staircase. For example:
0 0 0
0 0 x
with npiv = 2 has a Stair = [ 0 1 2 ] on output, to reflect the C block:
- C c
- - C
A tol of zero means that any nonzero norm (however small) is accepted;
only exact zero columns are flagged as dead. A negative tol means that
the norms are ignored; a column is never flagged dead. The default tol
is set elsewhere as 20 * (m+1) * eps * max column 2-norm of A.
LAPACK's dlarf* routines are used to construct and apply the Householder
reflections. The panel size (block size) is provided as an input
parameter, which defines the number of Householder vectors in a panel.
However, when the front is small (or when the remaining part
of a front is small) the block size is increased to include the entire
front. "Small" is defined, below, as fronts with fewer than 5000 entries.
NOTE: this function does not check its inputs. If the caller runs out of
memory and passes NULL pointers, this function will segfault.
*/
#include "spqr_subset.hpp"
#include "iostream"
#define SMALL 5000
#define MINCHUNK 4
#define MINCHUNK_RATIO 4
// =============================================================================
// === spqr_private_house ======================================================
// =============================================================================
// Construct a Householder reflection H = I - tau * v * v' such that H*x is
// reduced to zero except for the first element. Returns X [0] = the first
// entry of H*x, and X [1:n-1] = the Householder vector V [1:n-1], where
// V [0] = 1. If X [1:n-1] is zero, then the H=I (tau = 0) is returned,
// and V [1:n-1] is all zero. In MATLAB (1-based indexing), ignoring the
// rescaling done in dlarfg/zlarfg, this function computes the following:
/*
function [x, tau] = house (x)
n = length (x) ;
beta = norm (x) ;
if (x (1) > 0)
beta = -beta ;
end
tau = (beta - x (1)) / beta ;
x (2:n) = x (2:n) / (x (1) - beta) ;
x (1) = beta ;
*/
// Note that for the complex case, the reflection must be applied as H'*x,
// which requires that tau be conjugated in spqr_private_apply1.
//
// This function performs about 3*n+2 flops
inline double spqr_private_larfg (Int n, double *X, cholmod_common *cc)
{
double tau = 0 ;
BLAS_INT N = n, one = 1 ;
if (CHECK_BLAS_INT && !EQ (N,n))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_DLARFG (&N, X, X + 1, &one, &tau) ;
}
return (tau) ;
}
inline Complex spqr_private_larfg (Int n, Complex *X, cholmod_common *cc)
{
Complex tau = 0 ;
BLAS_INT N = n, one = 1 ;
if (CHECK_BLAS_INT && !EQ (N,n))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_ZLARFG (&N, X, X + 1, &one, &tau) ;
}
return (tau) ;
}
template <typename Entry> Entry spqr_private_house // returns tau
(
// inputs, not modified
Int n,
// input/output
Entry *X, // size n
cholmod_common *cc
)
{
return (spqr_private_larfg (n, X, cc)) ;
}
// =============================================================================
// === spqr_private_apply1 =====================================================
// =============================================================================
// Apply a single Householder reflection; C = C - tau * v * v' * C. The
// algorithm used by dlarf is:
/*
function C = apply1 (C, v, tau)
w = C'*v ;
C = C - tau*v*w' ;
*/
// For the complex case, we need to apply H', which requires that tau be
// conjugated.
//
// If applied to a single column, this function performs 2*n-1 flops to
// compute w, and 2*n+1 to apply it to C, for a total of 4*n flops.
inline void spqr_private_larf (Int m, Int n, double *V, double tau,
double *C, Int ldc, double *W, cholmod_common *cc)
{
BLAS_INT M = m, N = n, LDC = ldc, one = 1 ;
char left = 'L' ;
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDC,ldc)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_DLARF (&left, &M, &N, V, &one, &tau, C, &LDC, W) ;
}
}
inline void spqr_private_larf (Int m, Int n, Complex *V, Complex tau,
Complex *C, Int ldc, Complex *W, cholmod_common *cc)
{
BLAS_INT M = m, N = n, LDC = ldc, one = 1 ;
char left = 'L' ;
Complex conj_tau = spqr_conj (tau) ;
if (CHECK_BLAS_INT && !(EQ (M,m) && EQ (N,n) && EQ (LDC,ldc)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_ZLARF (&left, &M, &N, V, &one, &conj_tau, C, &LDC, W) ;
}
}
template <typename Entry> void spqr_private_apply1
(
// inputs, not modified
Int m, // C is m-by-n
Int n,
Int ldc, // leading dimension of C
Entry *V, // size m, Householder vector V
Entry tau, // Householder coefficient
// input/output
Entry *C, // size m-by-n
// workspace, not defined on input or output
Entry *W, // size n
cholmod_common *cc
)
{
Entry vsave ;
if (m <= 0 || n <= 0)
{
return ; // nothing to do
}
vsave = V [0] ; // temporarily restore unit diagonal of V
V [0] = 1 ;
spqr_private_larf (m, n, V, tau, C, ldc, W, cc) ;
V [0] = vsave ; // restore V [0]
}
// =============================================================================
// === spqr_front ==============================================================
// =============================================================================
// Factorize a front F into a sequence of Householder vectors H, and an upper
// trapezoidal matrix R. R may be a squeezed upper trapezoidal matrix if any
// of the leading npiv columns are linearly dependent. Returns the row index
// rank that indicates the first entry in C, which is F (rank,npiv), or 0
// on error.
template <typename Entry> Int spqr_front
(
// input, not modified
Int m, // F is m-by-n with leading dimension m
Int n,
Int npiv, // number of pivot columns
double tol, // a column is flagged as dead if its norm is <= tol
Int ntol, // apply tol only to first ntol pivot columns
Int fchunk, // block size for compact WY Householder reflections,
// treated as 1 if fchunk <= 1
// input/output
Entry *F, // frontal matrix F of size m-by-n
Int *Stair, // size n, entries F (Stair[k]:m-1, k) are all zero,
// for each k = 0:n-1, and remain zero on output.
char *Rdead, // size npiv; all zero on input. If k is dead,
// Rdead [k] is set to 1
// output, not defined on input
Entry *Tau, // size n, Householder coefficients
// workspace, undefined on input and output
Entry *W, // size b*n, where b = min (fchunk,n,m)
// input/output
double *wscale,
double *wssq,
cholmod_common *cc // for cc->hypotenuse function
)
{
Entry tau ;
double wk ;
Entry *V ;
Int k, t, g, g1, nv, k1, k2, i, t0, vzeros, mleft, nleft, vsize, minchunk,
rank ;
// NOTE: inputs are not checked for NULL (except if debugging enabled)
ASSERT (F != NULL) ;
ASSERT (Stair != NULL) ;
ASSERT (Rdead != NULL) ;
ASSERT (Tau != NULL) ;
ASSERT (W != NULL) ;
ASSERT (m >= 0 && n >= 0) ;
npiv = MAX (0, npiv) ; // npiv must be between 0 and n
npiv = MIN (n, npiv) ;
g1 = 0 ; // row index of first queued-up Householder
k1 = 0 ; // pending Householders are in F (g1:t, k1:k2-1)
k2 = 0 ;
V = F ; // Householder vectors start here
g = 0 ; // number of good Householders found
nv = 0 ; // number of Householder reflections queued up
vzeros = 0 ; // number of explicit zeros in queued-up H's
t = 0 ; // staircase of current column
fchunk = MAX (fchunk, 1) ;
minchunk = MAX (MINCHUNK, fchunk/MINCHUNK_RATIO) ;
rank = MIN (m,npiv) ; // F (rank,npiv) is the first entry in C. rank
// is also the number of rows in the R block,
// and the number of good pivot columns found.
ntol = MIN (ntol, npiv) ; // Note ntol can be negative, which means to
// not use tol at all.
PR (("Front %ld by %ld with %ld pivots\n", m, n, npiv)) ;
for (k = 0 ; k < n ; k++)
{
// ---------------------------------------------------------------------
// reduce the kth column of F to eliminate all but "diagonal" F (g,k)
// ---------------------------------------------------------------------
// get the staircase for the kth column, and operate on the "diagonal"
// F (g,k); eliminate F (g+1:t-1, k) to zero
t0 = t ; // t0 = staircase of column k-1
t = Stair [k] ; // t = staircase of this column k
PR (("k %ld g %ld m %ld n %ld npiv %ld\n", k, g, m, n, npiv)) ;
if (g >= m)
{
// F (g,k) is outside the matrix, so we're done. If this happens
// when k < npiv, then there is no contribution block.
PR (("hit the wall, npiv: %ld k %ld rank %ld\n", npiv, k, rank)) ;
for ( ; k < npiv ; k++)
{
Rdead [k] = 1 ;
Stair [k] = 0 ; // remaining pivot columns all dead
Tau [k] = 0 ;
}
for ( ; k < n ; k++)
{
Stair [k] = m ; // non-pivotal columns
Tau [k] = 0 ;
}
ASSERT (nv == 0) ; // there can be no pending updates
return (rank) ;
}
// if t < g+1, then this column is all zero; fix staircase so that R is
// always inside the staircase.
t = MAX (g+1,t) ;
Stair [k] = t ;
// ---------------------------------------------------------------------
// If t just grew a lot since the last t, apply H now to all of F
// ---------------------------------------------------------------------
// vzeros is the number of zero entries in V after including the next
// Householder vector. If it would exceed 50% of the size of V,
// apply the pending Householder reflections now, but only if
// enough vectors have accumulated.
vzeros += nv * (t - t0) ;
if (nv >= minchunk)
{
vsize = (nv*(nv+1))/2 + nv*(t-g1-nv) ;
if (vzeros > MAX (16, vsize/2))
{
// apply pending block of Householder reflections
PR (("(1) apply k1 %ld k2 %ld\n", k1, k2)) ;
spqr_larftb (
0, // method 0: Left, Transpose
t0-g1, n-k2, nv, m, m,
V, // F (g1:t-1, k1:k1+nv-1)
&Tau [k1], // Tau (k1:k-1)
&F [INDEX (g1,k2,m)], // F (g1:t-1, k2:n-1)
W, cc) ; // size nv*nv + nv*(n-k2)
nv = 0 ; // clear queued-up Householder reflections
vzeros = 0 ;
}
}
// ---------------------------------------------------------------------
// find a Householder reflection that reduces column k
// ---------------------------------------------------------------------
tau = spqr_private_house (t-g, &F [INDEX (g,k,m)], cc) ;
// ---------------------------------------------------------------------
// check to see if the kth column is OK
// ---------------------------------------------------------------------
if (k < ntol && (wk = spqr_abs (F [INDEX (g,k,m)], cc)) <= tol)
{
// -----------------------------------------------------------------
// norm (F (g:t-1, k)) is too tiny; the kth pivot column is dead
// -----------------------------------------------------------------
// keep track of the 2-norm of w, the dead column 2-norms
if (wk != 0)
{
// see also LAPACK's dnrm2 function
if ((*wscale) == 0)
{
// this is the nonzero first entry in w
(*wssq) = 1 ;
}
if ((*wscale) < wk)
{
double rr = (*wscale) / wk ;
(*wssq) = 1 + (*wssq) * rr * rr ;
(*wscale) = wk ;
}
else
{
double rr = wk / (*wscale) ;
(*wssq) += rr * rr ;
}
}
// zero out F (g:m-1,k) and flag it as dead
for (i = g ; i < m ; i++)
{
// This is not strictly necessary. On output, entries outside
// the staircase are ignored.
F [INDEX (i,k,m)] = 0 ;
}
Stair [k] = 0 ;
Tau [k] = 0 ;
Rdead [k] = 1 ;
if (nv > 0)
{
// apply pending block of Householder reflections
PR (("(2) apply k1 %ld k2 %ld\n", k1, k2)) ;
spqr_larftb (
0, // method 0: Left, Transpose
t0-g1, n-k2, nv, m, m,
V, // F (g1:t-1, k1:k1+nv-1)
&Tau [k1], // Tau (k1:k-1)
&F [INDEX (g1,k2,m)], // F (g1:t-1, k2:n-1)
W, cc) ; // size nv*nv + nv*(n-k2)
nv = 0 ; // clear queued-up Householder reflections
vzeros = 0 ;
}
}
else
{
// -----------------------------------------------------------------
// one more good pivot column found
// -----------------------------------------------------------------
Tau [k] = tau ; // save the Householder coefficient
if (nv == 0)
{
// start the queue of pending Householder updates, and define
// the current panel as k1:k2-1
g1 = g ; // first row of V
k1 = k ; // first column of V
k2 = MIN (n, k+fchunk) ; // k2-1 is last col in panel
V = &F [INDEX (g1,k1,m)] ; // pending V starts here
// check for switch to unblocked code
mleft = m-g1 ; // number of rows left
nleft = n-k1 ; // number of columns left
if (mleft * (nleft-(fchunk+4)) < SMALL || mleft <= fchunk/2
|| fchunk <= 1)
{
// remaining matrix is small; switch to unblocked code by
// including the rest of the matrix in the panel. Always
// use unblocked code if fchunk <= 1.
k2 = n ;
}
}
nv++ ; // one more pending update; V is F (g1:t-1, k1:k1+nv-1)
// -----------------------------------------------------------------
// keep track of "pure" flops, for performance testing only
// -----------------------------------------------------------------
// The Householder vector is of length t-g, including the unit
// diagonal, and takes 3*(t-g) flops to compute. It will be
// applied as a block, but compute the "pure" flops by assuming
// that this single Householder vector is computed and then applied
// just by itself to the rest of the frontal matrix (columns
// k+1:n-1, or n-k-1 columns). Applying the Householder reflection
// to just one column takes 4*(t-g) flops. This computation only
// works if TBB is disabled, merely because it uses a global
// variable to keep track of the flop count. If TBB is used, this
// computation may result in a race condition; it is disabled in
// that case.
FLOP_COUNT ((t-g) * (3 + 4 * (n-k-1))) ;
// -----------------------------------------------------------------
// apply the kth Householder reflection to the current panel
// -----------------------------------------------------------------
// F (g:t-1, k+1:k2-1) -= v * tau * v' * F (g:t-1, k+1:k2-1), where
// v is stored in F (g:t-1,k). This applies just one reflection
// to the current panel.
PR (("apply 1: k %ld\n", k)) ;
spqr_private_apply1 (t-g, k2-k-1, m, &F [INDEX (g,k,m)], tau,
&F [INDEX (g,k+1,m)], W, cc) ;
g++ ; // one more pivot found
// -----------------------------------------------------------------
// apply the Householder reflections if end of panel reached
// -----------------------------------------------------------------
if (k == k2-1 || g == m) // or if last pivot is found
{
// apply pending block of Householder reflections
PR (("(3) apply k1 %ld k2 %ld\n", k1, k2)) ;
spqr_larftb (
0, // method 0: Left, Transpose
t-g1, n-k2, nv, m, m,
V, // F (g1:t-1, k1:k1+nv-1)
&Tau [k1], // Tau (k1:k2-1)
&F [INDEX (g1,k2,m)], // F (g1:t-1, k2:n-1)
W, cc) ; // size nv*nv + nv*(n-k2)
nv = 0 ; // clear queued-up Householder reflections
vzeros = 0 ;
}
}
// ---------------------------------------------------------------------
// determine the rank of the pivot columns
// ---------------------------------------------------------------------
if (k == npiv-1)
{
// the rank is the number of good columns found in the first
// npiv columns. It is also the number of rows in the R block.
// F (rank,npiv) is first entry in the C block.
rank = g ;
PR (("rank of Front pivcols: %ld\n", rank)) ;
}
}
if (CHECK_BLAS_INT && !cc->blas_ok)
{
// This cannot occur if the BLAS_INT and the Int are the same integer.
// In that case, CHECK_BLAS_INT is FALSE at compile-time, and the
// compiler will then remove this as dead code.
ERROR (CHOLMOD_INVALID, "problem too large for the BLAS") ;
return (0) ;
}
return (rank) ;
}
// =============================================================================
template Int spqr_front <double>
(
// input, not modified
Int m, // F is m-by-n with leading dimension m
Int n,
Int npiv, // number of pivot columns
double tol, // a column is flagged as dead if its norm is <= tol
Int ntol, // apply tol only to first ntol pivot columns
Int fchunk, // block size for compact WY Householder reflections,
// treated as 1 if fchunk <= 1 (in which case the
// unblocked code is used).
// input/output
double *F, // frontal matrix F of size m-by-n
Int *Stair, // size n, entries F (Stair[k]:m-1, k) are all zero,
// and remain zero on output.
char *Rdead, // size npiv; all zero on input. If k is dead,
// Rdead [k] is set to 1
// output, not defined on input
double *Tau, // size n, Householder coefficients
// workspace, undefined on input and output
double *W, // size b*n, where b = min (fchunk,n,m)
// input/output
double *wscale,
double *wssq,
cholmod_common *cc // for cc->hypotenuse function
) ;
// =============================================================================
template Int spqr_front <Complex>
(
// input, not modified
Int m, // F is m-by-n with leading dimension m
Int n,
Int npiv, // number of pivot columns
double tol, // a column is flagged as dead if its norm is <= tol
Int ntol, // apply tol only to first ntol pivot columns
Int fchunk, // block size for compact WY Householder reflections,
// treated as 1 if fchunk <= 1 (in which case the
// unblocked code is used).
// input/output
Complex *F, // frontal matrix F of size m-by-n
Int *Stair, // size n, entries F (Stair[k]:m-1, k) are all zero,
// and remain zero on output.
char *Rdead, // size npiv; all zero on input. If k is dead,
// Rdead [k] is set to 1
// output, not defined on input
Complex *Tau, // size n, Householder coefficients
// workspace, undefined on input and output
Complex *W, // size b*n, where b = min (fchunk,n,m)
// input/output
double *wscale,
double *wssq,
cholmod_common *cc // for cc->hypotenuse function
) ;

0
spqr_mini/spqr_front.h
View File

236
spqr_mini/spqr_larftb.cpp
View File

@ -1,236 +0,0 @@
// =============================================================================
// === spqr_larftb =============================================================
// =============================================================================
// Apply a set of Householder reflections to a matrix. Given the vectors
// V and coefficients Tau, construct the matrix T and then apply the updates.
// In MATLAB (1-based indexing), this function computes the following:
/*
function C = larftb (C, V, Tau, method)
[v k] = size (V) ;
[m n] = size (C) ;
% construct T for the compact WY representation
V = tril (V,-1) + eye (v,k) ;
T = zeros (k,k) ;
T (1,1) = Tau (1) ;
for j = 2:k
tau = Tau (j) ;
z = -tau * V (:, 1:j-1)' * V (:,j) ;
T (1:j-1,j) = T (1:j-1,1:j-1) * z ;
T (j,j) = tau ;
end
% apply the updates
if (method == 0)
C = C - V * T' * V' * C ; % method 0: Left, Transpose
elseif (method == 1)
C = C - V * T * V' * C ; % method 1: Left, No Transpose
elseif (method == 2)
C = C - C * V * T' * V' ; % method 2: Right, Transpose
elseif (method == 3)
C = C - C * V * T * V' ; % method 3: Right, No Transpose
end
*/
#include "spqr_subset.hpp"
inline void spqr_private_larft (char direct, char storev, Int n, Int k,
double *V, Int ldv, double *Tau, double *T, Int ldt, cholmod_common *cc)
{
BLAS_INT N = n, K = k, LDV = ldv, LDT = ldt ;
if (CHECK_BLAS_INT &&
!(EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) && EQ (LDT,ldt)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_DLARFT (&direct, &storev, &N, &K, V, &LDV, Tau, T, &LDT) ;
}
}
inline void spqr_private_larft (char direct, char storev, Int n, Int k,
Complex *V, Int ldv, Complex *Tau, Complex *T, Int ldt, cholmod_common *cc)
{
BLAS_INT N = n, K = k, LDV = ldv, LDT = ldt ;
if (CHECK_BLAS_INT &&
!(EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) && EQ (LDT,ldt)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_ZLARFT (&direct, &storev, &N, &K, V, &LDV, Tau, T, &LDT) ;
}
}
inline void spqr_private_larfb (char side, char trans, char direct, char storev,
Int m, Int n, Int k, double *V, Int ldv, double *T, Int ldt, double *C,
Int ldc, double *Work, Int ldwork, cholmod_common *cc)
{
BLAS_INT M = m, N = n, K = k, LDV = ldv, LDT = ldt, LDC = ldc,
LDWORK = ldwork ;
if (CHECK_BLAS_INT &&
!(EQ (M,m) && EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) &&
EQ (LDT,ldt) && EQ (LDV,ldv) && EQ (LDWORK,ldwork)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_DLARFB (&side, &trans, &direct, &storev, &M, &N, &K, V, &LDV,
T, &LDT, C, &LDC, Work, &LDWORK) ;
}
}
inline void spqr_private_larfb (char side, char trans, char direct, char storev,
Int m, Int n, Int k, Complex *V, Int ldv, Complex *T, Int ldt, Complex *C,
Int ldc, Complex *Work, Int ldwork, cholmod_common *cc)
{
char tr = (trans == 'T') ? 'C' : 'N' ; // change T to C
BLAS_INT M = m, N = n, K = k, LDV = ldv, LDT = ldt, LDC = ldc,
LDWORK = ldwork ;
if (CHECK_BLAS_INT &&
!(EQ (M,m) && EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) &&
EQ (LDT,ldt) && EQ (LDV,ldv) && EQ (LDWORK,ldwork)))
{
cc->blas_ok = FALSE ;
}
if (!CHECK_BLAS_INT || cc->blas_ok)
{
LAPACK_ZLARFB (&side, &tr, &direct, &storev, &M, &N, &K, V, &LDV,
T, &LDT, C, &LDC, Work, &LDWORK) ;
}
}
// =============================================================================
template <typename Entry> void spqr_larftb
(
// inputs, not modified (V is modified and then restored on output)
int method, // 0,1,2,3
Int m, // C is m-by-n
Int n,
Int k, // V is v-by-k
// for methods 0 and 1, v = m,
// for methods 2 and 3, v = n
Int ldc, // leading dimension of C
Int ldv, // leading dimension of V
Entry *V, // V is v-by-k, unit lower triangular (diag not stored)
Entry *Tau, // size k, the k Householder coefficients
// input/output
Entry *C, // C is m-by-n, with leading dimension ldc
// workspace, not defined on input or output
Entry *W, // for methods 0,1: size k*k + n*k
// for methods 2,3: size k*k + m*k
cholmod_common *cc
)
{
Entry *T, *Work ;
// -------------------------------------------------------------------------
// check inputs and split up workspace
// -------------------------------------------------------------------------
if (m <= 0 || n <= 0 || k <= 0)
{
return ; // nothing to do
}
T = W ; // triangular k-by-k matrix for block reflector
Work = W + k*k ; // workspace of size n*k or m*k for larfb
// -------------------------------------------------------------------------
// construct and apply the k-by-k upper triangular matrix T
// -------------------------------------------------------------------------
// larft and larfb are always used "Forward" and "Columnwise"
if (method == SPQR_QTX)
{
ASSERT (m >= k) ;
spqr_private_larft ('F', 'C', m, k, V, ldv, Tau, T, k, cc) ;
// Left, Transpose, Forward, Columwise:
spqr_private_larfb ('L', 'T', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
Work, n, cc) ;
}
else if (method == SPQR_QX)
{
ASSERT (m >= k) ;
spqr_private_larft ('F', 'C', m, k, V, ldv, Tau, T, k, cc) ;
// Left, No Transpose, Forward, Columwise:
spqr_private_larfb ('L', 'N', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
Work, n, cc) ;
}
else if (method == SPQR_XQT)
{
ASSERT (n >= k) ;
spqr_private_larft ('F', 'C', n, k, V, ldv, Tau, T, k, cc) ;
// Right, Transpose, Forward, Columwise:
spqr_private_larfb ('R', 'T', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
Work, m, cc) ;
}
else if (method == SPQR_XQ)
{
ASSERT (n >= k) ;
spqr_private_larft ('F', 'C', n, k, V, ldv, Tau, T, k, cc) ;
// Right, No Transpose, Forward, Columwise:
spqr_private_larfb ('R', 'N', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
Work, m, cc) ;
}
}
// =============================================================================
template void spqr_larftb <double>
(
// inputs, not modified (V is modified and then restored on output)
int method, // 0,1,2,3
Int m, // C is m-by-n
Int n,
Int k, // V is v-by-k
// for methods 0 and 1, v = m,
// for methods 2 and 3, v = n
Int ldc, // leading dimension of C
Int ldv, // leading dimension of V
double *V, // V is v-by-k, unit lower triangular (diag not stored)
double *Tau, // size k, the k Householder coefficients
// input/output
double *C, // C is m-by-n, with leading dimension ldc
// workspace, not defined on input or output
double *W, // for methods 0,1: size k*k + n*k
// for methods 2,3: size k*k + m*k
cholmod_common *cc
) ;
// =============================================================================
template void spqr_larftb <Complex>
(
// inputs, not modified (V is modified and then restored on output)
int method, // 0,1,2,3
Int m, // C is m-by-n
Int n,
Int k, // V is v-by-k
// for methods 0 and 1, v = m,
// for methods 2 and 3, v = n
Int ldc, // leading dimension of C
Int ldv, // leading dimension of V
Complex *V, // V is v-by-k, unit lower triangular (diag not stored)
Complex *Tau, // size k, the k Householder coefficients
// input/output
Complex *C, // C is m-by-n, with leading dimension ldc
// workspace, not defined on input or output
Complex *W, // for methods 0,1: size k*k + n*k
// for methods 2,3: size k*k + m*k
cholmod_common *cc
) ;

0
spqr_mini/spqr_larftb.h
View File

397
spqr_mini/spqr_subset.hpp
View File

@ -1,397 +0,0 @@
// =============================================================================
// === spqr.hpp ================================================================
// =============================================================================
// Internal definitions and non-user-callable routines. This should not be
// included in the user's code.
#ifndef SPQR_INTERNAL_H
#define SPQR_INTERNAL_H
// -----------------------------------------------------------------------------
// include files
// -----------------------------------------------------------------------------
#include "UFconfig.h"
extern "C" {
#include "cholmod_core.h"
#include "cholmod_blas.h"
}
#include "SuiteSparseQR_definitions.h"
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <cstring>
#include <complex>
typedef std::complex<double> Complex ;
// -----------------------------------------------------------------------------
// debugging and printing control
// -----------------------------------------------------------------------------
// force debugging off
#ifndef NDEBUG
#define NDEBUG
#endif
// force printing off
#ifndef NPRINT
#define NPRINT
#endif
// uncomment the following line to turn on debugging (SPQR will be slow!)
/*
#undef NDEBUG
*/
// uncomment the following line to turn on printing (LOTS of output!)
/*
#undef NPRINT
*/
// uncomment the following line to turn on expensive debugging (very slow!)
/*
#define DEBUG_EXPENSIVE
*/
// -----------------------------------------------------------------------------
// Int is defined at UF_long, from UFconfig.h
// -----------------------------------------------------------------------------
#define Int UF_long
#define Int_max UF_long_max
// -----------------------------------------------------------------------------
// basic macros
// -----------------------------------------------------------------------------
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define EMPTY (-1)
#define TRUE 1
#define FALSE 0
#define IMPLIES(p,q) (!(p) || (q))
// NULL should already be defined, but ensure it is here.
#ifndef NULL
#define NULL ((void *) 0)
#endif
// column-major indexing; A[i,j] is A (INDEX (i,j,lda))
#define INDEX(i,j,lda) ((i) + ((j)*(lda)))
// FLIP is a "negation about -1", and is used to mark an integer i that is
// normally non-negative. FLIP (EMPTY) is EMPTY. FLIP of a number > EMPTY
// is negative, and FLIP of a number < EMTPY is positive. FLIP (FLIP (i)) = i
// for all integers i. UNFLIP (i) is >= EMPTY.
#define EMPTY (-1)
#define FLIP(i) (-(i)-2)
#define UNFLIP(i) (((i) < EMPTY) ? FLIP (i) : (i))
// -----------------------------------------------------------------------------
// additional include files
// -----------------------------------------------------------------------------
#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif
#define ITYPE CHOLMOD_LONG
#define DTYPE CHOLMOD_DOUBLE
#define ID UF_long_id
// -----------------------------------------------------------------------------
#define ERROR(status,msg) \
printf ("CHOLMOD error: %s\n",msg) // Kai: disable cholmod_l_error to prevent from including tons of files
// Check a pointer and return if null. Set status to invalid, unless the
// status is already "out of memory"
#define RETURN_IF_NULL(A,result) \
{ \
if ((A) == NULL) \
{ \
if (cc->status != CHOLMOD_OUT_OF_MEMORY) \
{ \
ERROR (CHOLMOD_INVALID, NULL) ; \
} \
return (result) ; \
} \
}
// Return if Common is NULL or invalid
#define RETURN_IF_NULL_COMMON(result) \
{ \
if (cc == NULL) \
{ \
return (result) ; \
} \
if (cc->itype != ITYPE || cc->dtype != DTYPE) \
{ \
cc->status = CHOLMOD_INVALID ; \
return (result) ; \
} \
}
#define RETURN_IF_XTYPE_INVALID(A,result) \
{ \
if (A->xtype != xtype) \
{ \
ERROR (CHOLMOD_INVALID, "invalid xtype") ; \
return (result) ; \
} \
}
// -----------------------------------------------------------------------------
// debugging and printing macros
// -----------------------------------------------------------------------------
#ifndef NDEBUG
#ifdef MATLAB_MEX_FILE
// #define ASSERT(e) mxAssert (e, "error: ")
extern char spqr_mx_debug_string [200] ;
char *spqr_mx_id (int line) ;
#define ASSERT(e) \
((e) ? (void) 0 : \
mexErrMsgIdAndTxt (spqr_mx_id (__LINE__), \
"assert: (" #e ") file:" __FILE__ ))
#else
#include <assert.h>
#define ASSERT(e) assert (e)
#endif
#define DEBUG(e) e
#ifdef DEBUG_EXPENSIVE
#define DEBUG2(e) e
#define ASSERT2(e) ASSERT(e)
#else
#define DEBUG2(e)
#define ASSERT2(e)
#endif
#else
#define ASSERT(e)
#define ASSERT2(e)
#define DEBUG(e)
#define DEBUG2(e)
#endif
#ifndef NPRINT
#ifdef MATLAB_MEX_FILE
#define PR(e) mexPrintf e
#else
#define PR(e) printf e
#endif
#define PRVAL(e) spqrDebug_print (e)
#else
#define PR(e)
#define PRVAL(e)
#endif
// -----------------------------------------------------------------------------
// For counting flops; disabled if TBB is used or timing not enabled
// -----------------------------------------------------------------------------
#if defined(TIMING)
#define FLOP_COUNT(f) { if (cc->SPQR_grain <= 1) cc->other1 [0] += (f) ; }
#else
#define FLOP_COUNT(f)
#endif
template <typename Entry> void spqr_larftb
(
// inputs, not modified (V is modified and then restored on output)
int method, // 0,1,2,3
Int m, // C is m-by-n
Int n,
Int k, // V is v-by-k
// for methods 0 and 1, v = m,
// for methods 2 and 3, v = n
Int ldc, // leading dimension of C
Int ldv, // leading dimension of V
Entry *V, // V is v-by-k, unit lower triangular (diag not stored)
Entry *Tau, // size k, the k Householder coefficients
// input/output
Entry *C, // C is m-by-n, with leading dimension ldc
// workspace, not defined on input or output
Entry *W, // for methods 0,1: size k*k + n*k
// for methods 2,3: size k*k + m*k
cholmod_common *cc
) ;
// returns rank of F, or 0 on error
template <typename Entry> Int spqr_front
(
// input, not modified
Int m, // F is m-by-n with leading dimension m
Int n,
Int npiv, // number of pivot columns
double tol, // a column is flagged as dead if its norm is <= tol
Int ntol, // apply tol only to first ntol pivot columns
Int fchunk, // block size for compact WY Householder reflections,
// treated as 1 if fchunk <= 1
// input/output
Entry *F, // frontal matrix F of size m-by-n
Int *Stair, // size n, entries F (Stair[k]:m-1, k) are all zero,
// and remain zero on output.
char *Rdead, // size npiv; all zero on input. If k is dead,
// Rdead [k] is set to 1
// output, not defined on input
Entry *Tau, // size n, Householder coefficients
// workspace, undefined on input and output
Entry *W, // size b*(n+b), where b = min (fchunk,n,m)
// input/output
double *wscale,
double *wssq,
cholmod_common *cc // for cc->hypotenuse function
) ;
// =============================================================================
// === spqr_conj ===============================================================
// =============================================================================
inline double spqr_conj (double x)
{
return (x) ;
}
inline Complex spqr_conj (Complex x)
{
return (std::conj (x)) ;
}
// =============================================================================
// === spqr_abs ================================================================
// =============================================================================
inline double spqr_abs (double x, cholmod_common *cc) // cc is unused
{
return (fabs (x)) ;
}
inline double spqr_abs (Complex x, cholmod_common *cc)
{
return (cc->hypotenuse (x.real ( ), x.imag ( ))) ;
}
// =============================================================================
// === BLAS interface ==========================================================
// =============================================================================
// To compile SuiteSparseQR with 64-bit BLAS, use -DBLAS64. See also
// CHOLMOD/Include/cholmod_blas.h
extern "C" {
#include "cholmod_blas.h"
}
#ifdef SUN64
#define BLAS_DNRM2 dnrm2_64_
#define LAPACK_DLARF dlarf_64_
#define LAPACK_DLARFG dlarfg_64_
#define LAPACK_DLARFT dlarft_64_
#define LAPACK_DLARFB dlarfb_64_
#define BLAS_DZNRM2 dznrm2_64_
#define LAPACK_ZLARF zlarf_64_
#define LAPACK_ZLARFG zlarfg_64_
#define LAPACK_ZLARFT zlarft_64_
#define LAPACK_ZLARFB zlarfb_64_
#elif defined (BLAS_NO_UNDERSCORE)
#define BLAS_DNRM2 dnrm2
#define LAPACK_DLARF dlarf
#define LAPACK_DLARFG dlarfg
#define LAPACK_DLARFT dlarft
#define LAPACK_DLARFB dlarfb
#define BLAS_DZNRM2 dznrm2
#define LAPACK_ZLARF zlarf
#define LAPACK_ZLARFG zlarfg
#define LAPACK_ZLARFT zlarft
#define LAPACK_ZLARFB zlarfb
#else
#define BLAS_DNRM2 dnrm2_
#define LAPACK_DLARF dlarf_
#define LAPACK_DLARFG dlarfg_
#define LAPACK_DLARFT dlarft_
#define LAPACK_DLARFB dlarfb_
#define BLAS_DZNRM2 dznrm2_
#define LAPACK_ZLARF zlarf_
#define LAPACK_ZLARFG zlarfg_
#define LAPACK_ZLARFT zlarft_
#define LAPACK_ZLARFB zlarfb_
#endif
// =============================================================================
// === BLAS and LAPACK prototypes ==============================================
// =============================================================================
extern "C"
{
void LAPACK_DLARFT (char *direct, char *storev, BLAS_INT *n, BLAS_INT *k,
double *V, BLAS_INT *ldv, double *Tau, double *T, BLAS_INT *ldt) ;
void LAPACK_ZLARFT (char *direct, char *storev, BLAS_INT *n, BLAS_INT *k,
Complex *V, BLAS_INT *ldv, Complex *Tau, Complex *T, BLAS_INT *ldt) ;
void LAPACK_DLARFB (char *side, char *trans, char *direct, char *storev,
BLAS_INT *m, BLAS_INT *n, BLAS_INT *k, double *V, BLAS_INT *ldv,
double *T, BLAS_INT *ldt, double *C, BLAS_INT *ldc, double *Work,
BLAS_INT *ldwork) ;
void LAPACK_ZLARFB (char *side, char *trans, char *direct, char *storev,
BLAS_INT *m, BLAS_INT *n, BLAS_INT *k, Complex *V, BLAS_INT *ldv,
Complex *T, BLAS_INT *ldt, Complex *C, BLAS_INT *ldc, Complex *Work,
BLAS_INT *ldwork) ;
double BLAS_DNRM2 (BLAS_INT *n, double *X, BLAS_INT *incx) ;
double BLAS_DZNRM2 (BLAS_INT *n, Complex *X, BLAS_INT *incx) ;
void LAPACK_DLARFG (BLAS_INT *n, double *alpha, double *X, BLAS_INT *incx,
double *tau) ;
void LAPACK_ZLARFG (BLAS_INT *n, Complex *alpha, Complex *X, BLAS_INT *incx,
Complex *tau) ;
void LAPACK_DLARF (char *side, BLAS_INT *m, BLAS_INT *n, double *V,
BLAS_INT *incv, double *tau, double *C, BLAS_INT *ldc, double *Work) ;
void LAPACK_ZLARF (char *side, BLAS_INT *m, BLAS_INT *n, Complex *V,
BLAS_INT *incv, Complex *tau, Complex *C, BLAS_INT *ldc, Complex *Work) ;
}
#endif

2
tests/Makefile.am
View File

@ -36,7 +36,7 @@ noinst_PROGRAMS = timeGaussianFactorGraph timeLinearOnDataset
# rules to build unit tests
#----------------------------------------------------------------------------------------------------
TESTS = $(check_PROGRAMS)
AM_CPPFLAGS = -I$(boost) -I$(top_srcdir)/..
AM_CPPFLAGS = -I$(boost) -I$(SparseInc) -I$(top_srcdir)/..
AM_LDFLAGS = $(BOOST_LDFLAGS)
if USE_ACCELERATE_MACOS