diff --git a/Makefile.am b/Makefile.am index 100e6f32a..918f20854 100644 --- a/Makefile.am +++ b/Makefile.am @@ -9,10 +9,10 @@ ACLOCAL_AMFLAGS = -I m4 AUTOMAKE_OPTIONS = foreign # All the sub-directories that need to be built -SUBDIRS = CppUnitLite colamd ldl spqr_mini base geometry inference linear nonlinear slam . tests wrap +SUBDIRS = CppUnitLite colamd spqr_mini base geometry inference linear nonlinear slam . tests wrap # And the corresponding libraries produced -SUBLIBS = colamd/libcolamd.la ldl/libldl.la \ +SUBLIBS = colamd/libcolamd.la \ base/libbase.la geometry/libgeometry.la inference/libinference.la \ linear/liblinear.la nonlinear/libnonlinear.la slam/libslam.la diff --git a/base/Makefile.am b/base/Makefile.am index 51109d268..c99daabd2 100644 --- a/base/Makefile.am +++ b/base/Makefile.am @@ -50,6 +50,10 @@ if USE_LAPACK AM_CPPFLAGS += -DGT_USE_LAPACK endif +if USE_LDL +AM_CPPFLAGS += -DGT_USE_LDL +endif + # On Mac, we compile using the BLAS/LAPACK headers in the Accelerate framework if USE_ACCELERATE_MACOS AM_CPPFLAGS += -F Accelerate @@ -61,7 +65,10 @@ endif TESTS = $(check_PROGRAMS) AM_DEFAULT_SOURCE_EXT = .cpp AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -LDADD = libbase.la ../ldl/libldl.la ../CppUnitLite/libCppUnitLite.a +LDADD = libbase.la ../CppUnitLite/libCppUnitLite.a +if USE_LDL +LDADD += libldl.la +endif if USE_BLAS_LINUX AM_LDFLAGS += -lcblas -latlas diff --git a/base/Matrix.cpp b/base/Matrix.cpp index 1e634b4f8..06502a6ae 100644 --- a/base/Matrix.cpp +++ b/base/Matrix.cpp @@ -34,7 +34,9 @@ extern "C" { #include #include -#include +#ifdef GT_USE_LDL +#include +#endif #include "Matrix.h" #include "Vector.h" @@ -972,6 +974,7 @@ Matrix RtR(const Matrix &A) } /* ************************************************************************* */ +#ifdef GT_USE_LDL Vector solve_ldl(const Matrix& M, const Vector& rhs) { int N = M.size1(); // size of the matrix @@ -1052,6 +1055,7 @@ Vector solve_ldl(const Matrix& M, const Vector& rhs) { return result; } +#endif /* * This is not ultra efficient, but not terrible, either. diff --git a/base/Matrix.h b/base/Matrix.h index ca71d5592..a2a517ed2 100644 --- a/base/Matrix.h +++ b/base/Matrix.h @@ -376,7 +376,9 @@ Matrix RtR(const Matrix& A); Matrix cholesky_inverse(const Matrix &A); /** Solve Ax=b with S.P.D. matrix using Davis' LDL code */ +#ifdef GT_USE_LDL Vector solve_ldl(const Matrix& A, const Vector& b); +#endif /** * Numerical exponential map, naive approach, not industrial strength !!! diff --git a/base/tests/testMatrix.cpp b/base/tests/testMatrix.cpp index b8dfea6f8..3314fbe7a 100644 --- a/base/tests/testMatrix.cpp +++ b/base/tests/testMatrix.cpp @@ -7,7 +7,6 @@ #include #include -//#include #include #include #include @@ -966,6 +965,7 @@ TEST( matrix, transposeMultiplyAdd ) } /* ************************************************************************* */ +#ifdef GT_USE_LDL TEST( matrix, LDL_factorization ) { // run demo inside Matrix.cpp code @@ -994,6 +994,7 @@ TEST( matrix, LDL_factorization ) { CHECK(assert_equal(expected, actual)); } +#endif /* ************************************************************************* */ TEST( matrix, linear_dependent ) diff --git a/configure.ac b/configure.ac index 22b54b079..aaf4b50e9 100644 --- a/configure.ac +++ b/configure.ac @@ -4,12 +4,11 @@ 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 ldl/Makefile spqr_mini/Makefile base/Makefile inference/Makefile linear/Makefile geometry/Makefile nonlinear/Makefile slam/Makefile tests/Makefile wrap/Makefile) +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) 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([ldl/ldl.cpp]) AC_CONFIG_SRCDIR([spqr_mini/spqr_front.cpp]) AC_CONFIG_SRCDIR([base/DSFVector.cpp]) AC_CONFIG_SRCDIR([geometry/Cal3_S2.cpp]) @@ -60,7 +59,7 @@ AC_ARG_ENABLE([blas], no) blas=false ;; *) AC_MSG_ERROR([bad value ${enableval} for --enable-blas]) ;; esac],[blas=false]) - +ak AM_CONDITIONAL([USE_BLAS], test x$blas = xtrue) AM_CONDITIONAL([USE_BLAS_MACOS], [test x$blas = xtrue && test x$ISMAC = xtrue]) AM_CONDITIONAL([USE_BLAS_LINUX], [test x$blas = xtrue && test x$ISMAC = xfalse]) @@ -92,6 +91,17 @@ AC_ARG_ENABLE([spqr], AM_CONDITIONAL([USE_SPQR], [test x$spqr = xtrue]) +# enable using LDL library from SuiteSparse +AC_ARG_ENABLE([ldl], + [ --enable-ldl Enable LDL library support], + [case "${enableval}" in + yes) ldl=true ;; + no) ldl=false ;; + *) AC_MSG_ERROR([bad value ${enableval} for --enable-debug]) ;; + esac],[ldl=false]) + +AM_CONDITIONAL([USE_LDL], [test x$ldl = xtrue]) + # enable profiling AC_ARG_ENABLE([profiling], [ --enable-profiling Enable profiling], diff --git a/geometry/Makefile.am b/geometry/Makefile.am index f6e66cf9d..b08d89810 100644 --- a/geometry/Makefile.am +++ b/geometry/Makefile.am @@ -39,8 +39,11 @@ AM_CPPFLAGS = -I$(boost) -I$(top_srcdir) -I$(top_srcdir)/base #---------------------------------------------------------------------------------------------------- TESTS = $(check_PROGRAMS) AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) -LDADD = libgeometry.la ../base/libbase.la ../ldl/libldl.la ../CppUnitLite/libCppUnitLite.a +LDADD = libgeometry.la ../base/libbase.la ../CppUnitLite/libCppUnitLite.a AM_DEFAULT_SOURCE_EXT = .cpp +if USE_LDL +LDADD += libldl.la +endif # rule to run an executable %.run: % $(LDADD) diff --git a/inference/Makefile.am b/inference/Makefile.am index a38dde6ca..7a1e7ac11 100644 --- a/inference/Makefile.am +++ b/inference/Makefile.am @@ -62,8 +62,11 @@ AM_CXXFLAGS = TESTS = $(check_PROGRAMS) AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) LDADD = libinference.la ../base/libbase.la -LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la ../ldl/libldl.la +LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la AM_DEFAULT_SOURCE_EXT = .cpp +if USE_LDL +LDADD += libldl.la +endif # rule to run an executable %.run: % $(LDADD) diff --git a/ldl/Makefile.am b/ldl/Makefile.am deleted file mode 100644 index c2bbfca7b..000000000 --- a/ldl/Makefile.am +++ /dev/null @@ -1,8 +0,0 @@ -#---------------------------------------------------------------------------------------------------- -# ldl -# replaced Makefile with automake for easy linking -#---------------------------------------------------------------------------------------------------- - -noinst_LTLIBRARIES = libldl.la -libldl_la_SOURCES = ldl.cpp -noinst_HEADERS = ldl.h UFconfig.h diff --git a/ldl/README.txt b/ldl/README.txt deleted file mode 100644 index 7be8dd1f0..000000000 --- a/ldl/README.txt +++ /dev/null @@ -1,136 +0,0 @@ -LDL Version 2.0: a sparse LDL' factorization and solve package. - Written in C, with both a C and MATLAB mexFunction interface. - -These routines are not terrifically fast (they do not use dense matrix kernels), -but the code is very short and concise. The purpose is to illustrate the -algorithms in a very concise and readable manner, primarily for educational -purposes. Although the code is very concise, this package is slightly faster -than the built-in sparse Cholesky factorization in MATLAB 6.5 (chol), when -using the same input permutation. - -Requires UFconfig, in the ../UFconfig directory relative to this directory. - -Quick start (Unix, or Windows with Cygwin): - - To compile, test, and install LDL, you may wish to first obtain a copy of - AMD v2.0 from http://www.cise.ufl.edu/research/sparse, and place it in the - ../AMD directory, relative to this directory. Next, type "make", which - will compile the LDL library and three demo main programs (one of which - requires AMD). It will also compile the LDL MATLAB mexFunction (if you - have MATLAB). Typing "make clean" will remove non-essential files. - AMD v2.0 or later is required. Its use is optional. - -Quick start (for MATLAB users); - - To compile, test, and install the LDL mexFunctions (ldlsparse and - ldlsymbol), start MATLAB in this directory and type ldl_install. - This works on any system supported by MATLAB. - --------------------------------------------------------------------------------- - -LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved. - -LDL License: - - Your use or distribution of LDL or any modified version of - LDL implies that you agree to this License. - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 - USA - - Permission is hereby granted to use or copy this program under the - terms of the GNU LGPL, provided that the Copyright, this License, - and the Availability of the original version is retained on all copies. - User documentation of any code that uses this code or any modified - version of this code must cite the Copyright, this License, the - Availability note, and "Used by permission." Permission to modify - the code and to distribute modified code is granted, provided the - Copyright, this License, and the Availability note are retained, - and a notice that the code was modified is included. - -Availability: - - http://www.cise.ufl.edu/research/sparse/ldl - -Acknowledgements: - - This work was supported by the National Science Foundation, under - grant CCR-0203270. - - Portions of this work were done while on sabbatical at Stanford University - and Lawrence Berkeley National Laboratory (with funding from the SciDAC - program). I would like to thank Gene Golub, Esmond Ng, and Horst Simon - for making this sabbatical possible. I would like to thank Pete Stewart - for his comments on a draft of this software and paper. - --------------------------------------------------------------------------------- -Files and directories in this distribution: --------------------------------------------------------------------------------- - - Documentation, and compiling: - - README.txt this file - Makefile for compiling LDL - ChangeLog changes since V1.0 (Dec 31, 2003) - License license - lesser.txt the GNU LGPL license - - ldl_userguide.pdf user guide in PDF - ldl_userguide.ps user guide in postscript - ldl_userguide.tex user guide in Latex - ldl.bib bibliography for user guide - - The LDL library itself: - - ldl.c the C-callable routines - ldl.h include file for any code that calls LDL - - A simple C main program that demonstrates how to use LDL: - - ldlsimple.c a stand-alone C program, uses the basic features of LDL - ldlsimple.out output of ldlsimple - - ldllsimple.c long integer version of ldlsimple.c - - Demo C program, for testing LDL and providing an example of its use - - ldlmain.c a stand-alone C main program that uses and tests LDL - Matrix a directory containing matrices used by ldlmain.c - ldlmain.out output of ldlmain - ldlamd.out output of ldlamd (ldlmain.c compiled with AMD) - ldllamd.out output of ldllamd (ldlmain.c compiled with AMD, long) - - MATLAB-related, not required for use in a regular C program - - Contents.m a list of the MATLAB-callable routines - ldl.m MATLAB help file for the LDL mexFunction - ldldemo.m MATLAB demo of how to use the LDL mexFunction - ldldemo.out diary output of ldldemo - ldltest.m to test the LDL mexFunction - ldltest.out diary output of ldltest - ldlmex.c the LDL mexFunction for MATLAB - ldlrow.m the numerical algorithm that LDL is based on - ldlmain2.m compiles and runs ldlmain.c as a MATLAB mexFunction - ldlmain2.out output of ldlmain2.m - ldlsymbolmex.c symbolic factorization using LDL (see SYMBFACT, ETREE) - ldlsymbol.m help file for the LDLSYMBOL mexFunction - - ldl_install.m compile, install, and test LDL functions - ldl_make.m compile LDL (ldlsparse and ldlsymbol) - - ldlsparse.m help for ldlsparse - -See ldl.c for a description of how to use the code from a C program. Type -"help ldl" in MATLAB for information on how to use LDL in a MATLAB program. diff --git a/ldl/UFconfig.h b/ldl/UFconfig.h deleted file mode 100644 index 7b5e79e54..000000000 --- a/ldl/UFconfig.h +++ /dev/null @@ -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 - -/* ========================================================================== */ -/* === 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.4.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.1 - * CSparse version 2.2.3 - * CXSparse version 2.2.3 - * KLU version 1.1.0 - * BTF version 1.1.0 - * LDL version 2.0.1 - * UFconfig version number is the same as SuiteSparse - * UMFPACK version 5.4.0 - * RBio version 1.1.2 - * UFcollection version 1.2.0 - * LINFACTOR version 1.1.0 - * MESHND version 1.1.1 - * SSMULT version 2.0.0 - * MATLAB_Tools no specific version number - * SuiteSparseQR version 1.1.2 - * - * 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 "May 20, 2009" -#define SUITESPARSE_VER_CODE(main,sub) ((main) * 1000 + (sub)) -#define SUITESPARSE_MAIN_VERSION 3 -#define SUITESPARSE_SUB_VERSION 4 -#define SUITESPARSE_SUBSUB_VERSION 0 -#define SUITESPARSE_VERSION \ - SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION,SUITESPARSE_SUB_VERSION) - -#ifdef __cplusplus -} -#endif -#endif diff --git a/ldl/UFconfig.mk b/ldl/UFconfig.mk deleted file mode 100644 index 3ca3e3940..000000000 --- a/ldl/UFconfig.mk +++ /dev/null @@ -1,354 +0,0 @@ -#=============================================================================== -# UFconfig.mk: common configuration file for the SuiteSparse -#=============================================================================== - -# This file contains all configuration settings for all packages authored or -# co-authored by Tim Davis at the University of Florida: -# -# Package Version Description -# ------- ------- ----------- -# AMD 1.2 or later approximate minimum degree ordering -# COLAMD 2.4 or later column approximate minimum degree ordering -# CCOLAMD 1.0 or later constrained column approximate minimum degree ordering -# CAMD any constrained approximate minimum degree ordering -# UMFPACK 4.5 or later sparse LU factorization, with the BLAS -# CHOLMOD any sparse Cholesky factorization, update/downdate -# KLU 0.8 or later sparse LU factorization, BLAS-free -# BTF 0.8 or later permutation to block triangular form -# LDL 1.2 or later concise sparse LDL' -# LPDASA any linear program solve (dual active set algorithm) -# CXSparse any extended version of CSparse (int/long, real/complex) -# SuiteSparseQR any sparse QR factorization -# -# The UFconfig directory and the above packages should all appear in a single -# directory, in order for the Makefile's within each package to find this file. -# -# To enable an option of the form "# OPTION = ...", edit this file and -# delete the "#" in the first column of the option you wish to use. - -#------------------------------------------------------------------------------ -# Generic configuration -#------------------------------------------------------------------------------ - -# C compiler and compiler flags: These will normally not give you optimal -# performance. You should select the optimization parameters that are best -# for your system. On Linux, use "CFLAGS = -O3 -fexceptions" for example. -CC ?= cc -# CFLAGS = -O (for example; see below for details) - -# C++ compiler (also uses CFLAGS) -CPLUSPLUS = g++ - -# ranlib, and ar, for generating libraries -RANLIB = ranlib -AR = ar cr - -# delete and rename a file -RM = rm -f -MV = mv -f - -# Fortran compiler (not normally required) -F77 = f77 -F77FLAGS = -O -F77LIB = - -# C and Fortran libraries -LIB = -lm - -# For compiling MATLAB mexFunctions (MATLAB 7.5 or later) -MEX = mex -O -largeArrayDims -lmwlapack -lmwblas - -# For compiling MATLAB mexFunctions (MATLAB 7.3 and 7.4) -# MEX = mex -O -largeArrayDims -lmwlapack - -# For MATLAB 7.2 or earlier, you must use one of these options: -# MEX = mex -O -lmwlapack -# MEX = mex -O - -# Which version of MAKE you are using (default is "make") -# MAKE = make -# MAKE = gmake - -#------------------------------------------------------------------------------ -# BLAS and LAPACK configuration: -#------------------------------------------------------------------------------ - -# UMFPACK and CHOLMOD both require the BLAS. CHOLMOD also requires LAPACK. -# See Kazushige Goto's BLAS at http://www.cs.utexas.edu/users/flame/goto/ or -# http://www.tacc.utexas.edu/~kgoto/ for the best BLAS to use with CHOLMOD. -# LAPACK is at http://www.netlib.org/lapack/ . You can use the standard -# Fortran LAPACK along with Goto's BLAS to obtain very good performance. -# CHOLMOD gets a peak numeric factorization rate of 3.6 Gflops on a 3.2 GHz -# Pentium 4 (512K cache, 4GB main memory) with the Goto BLAS, and 6 Gflops -# on a 2.5Ghz dual-core AMD Opteron. - -# These settings will probably not work, since there is no fixed convention for -# naming the BLAS and LAPACK library (*.a or *.so) files. - -# Using the Goto BLAS: -# BLAS = -lgoto -lgfortran -lgfortranbegin -lg2c - -# This is probably slow ... it might connect to the Standard Reference BLAS: -BLAS = -lblas -lgfortran -lgfortranbegin -lg2c -LAPACK = -llapack - -# Using non-optimized versions: -# BLAS = -lblas_plain -lgfortran -lgfortranbegin -lg2c -# LAPACK = -llapack_plain - -# The BLAS might not contain xerbla, an error-handling routine for LAPACK and -# the BLAS. Also, the standard xerbla requires the Fortran I/O library, and -# stops the application program if an error occurs. A C version of xerbla -# distributed with this software (UFconfig/xerbla/libcerbla.a) includes a -# Fortran-callable xerbla routine that prints nothing and does not stop the -# application program. This is optional. -# XERBLA = ../../UFconfig/xerbla/libcerbla.a - -# If you wish to use the XERBLA in LAPACK and/or the BLAS instead, -# use this option: -XERBLA = - -# If you wish to use the Fortran UFconfig/xerbla/xerbla.f instead, use this: -# XERBLA = ../../UFconfig/xerbla/libxerbla.a - -#------------------------------------------------------------------------------ -# METIS, optionally used by CHOLMOD -#------------------------------------------------------------------------------ - -# If you do not have METIS, or do not wish to use it in CHOLMOD, you must -# compile CHOLMOD with the -DNPARTITION flag. You must also use the -# "METIS =" option, below. - -# The path is relative to where it is used, in CHOLMOD/Lib, CHOLMOD/MATLAB, etc. -# You may wish to use an absolute path. METIS is optional. Compile -# CHOLMOD with -DNPARTITION if you do not wish to use METIS. -METIS_PATH = ../../metis-4.0 -METIS = ../../metis-4.0/libmetis.a - -# If you use CHOLMOD_CONFIG = -DNPARTITION then you must use the following -# options: -# METIS_PATH = -# METIS = - -#------------------------------------------------------------------------------ -# UMFPACK configuration: -#------------------------------------------------------------------------------ - -# Configuration flags for UMFPACK. See UMFPACK/Source/umf_config.h for details. -# -# -DNBLAS do not use the BLAS. UMFPACK will be very slow. -# -D'LONGBLAS=long' or -DLONGBLAS='long long' defines the integers used by -# LAPACK and the BLAS (defaults to 'int') -# -DNSUNPERF do not use the Sun Perf. Library (default is use it on Solaris) -# -DNPOSIX do not use POSIX routines sysconf and times. -# -DGETRUSAGE use getrusage -# -DNO_TIMER do not use any timing routines -# -DNRECIPROCAL do not multiply by the reciprocal -# -DNO_DIVIDE_BY_ZERO do not divide by zero - -UMFPACK_CONFIG = - -#------------------------------------------------------------------------------ -# CHOLMOD configuration -#------------------------------------------------------------------------------ - -# CHOLMOD Library Modules, which appear in libcholmod.a: -# Core requires: none -# Check requires: Core -# Cholesky requires: Core, AMD, COLAMD. optional: Partition, Supernodal -# MatrixOps requires: Core -# Modify requires: Core -# Partition requires: Core, CCOLAMD, METIS. optional: Cholesky -# Supernodal requires: Core, BLAS, LAPACK -# -# CHOLMOD test/demo Modules (all are GNU GPL, do not appear in libcholmod.a): -# Tcov requires: Core, Check, Cholesky, MatrixOps, Modify, Supernodal -# optional: Partition -# Valgrind same as Tcov -# Demo requires: Core, Check, Cholesky, MatrixOps, Supernodal -# optional: Partition -# -# Configuration flags: -# -DNCHECK do not include the Check module. License GNU LGPL -# -DNCHOLESKY do not include the Cholesky module. License GNU LGPL -# -DNPARTITION do not include the Partition module. License GNU LGPL -# also do not include METIS. -# -DNGPL do not include any GNU GPL Modules in the CHOLMOD library: -# -DNMATRIXOPS do not include the MatrixOps module. License GNU GPL -# -DNMODIFY do not include the Modify module. License GNU GPL -# -DNSUPERNODAL do not include the Supernodal module. License GNU GPL -# -# -DNPRINT do not print anything. -# -D'LONGBLAS=long' or -DLONGBLAS='long long' defines the integers used by -# LAPACK and the BLAS (defaults to 'int') -# -DNSUNPERF for Solaris only. If defined, do not use the Sun -# Performance Library - -CHOLMOD_CONFIG = - -#------------------------------------------------------------------------------ -# SuiteSparseQR configuration: -#------------------------------------------------------------------------------ - -# The SuiteSparseQR library can be compiled with the following options: -# -# -DNPARTITION do not include the CHOLMOD partition module -# -DNEXPERT do not include the functions in SuiteSparseQR_expert.cpp -# -DTIMING enable timing and flop counts -# -DHAVE_TBB enable the use of Intel's Threading Building Blocks (TBB) - -# default, without timing, without TBB: -SPQR_CONFIG = -# with timing and TBB: -# SPQR_CONFIG = -DTIMING -DHAVE_TBB -# with timing -# SPQR_CONFIG = -DTIMING - -# with TBB, you must select this: -# TBB = -ltbb -# without TBB: -TBB = - -# with timing, you must include the timing library: -# RTLIB = -lrt -# without timing -RTLIB = - -#------------------------------------------------------------------------------ -# Linux -#------------------------------------------------------------------------------ - -# Using default compilers: -# CC = gcc -CFLAGS = -O3 -fexceptions - -# alternatives: -# CFLAGS = -g -fexceptions \ - -Wall -W -Wshadow -Wmissing-prototypes -Wstrict-prototypes \ - -Wredundant-decls -Wnested-externs -Wdisabled-optimization -ansi -# CFLAGS = -O3 -fexceptions \ - -Wall -W -Werror -Wshadow -Wmissing-prototypes -Wstrict-prototypes \ - -Wredundant-decls -Wnested-externs -Wdisabled-optimization -ansi -# CFLAGS = -O3 -fexceptions -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -# CFLAGS = -O3 -# CFLAGS = -O3 -g -fexceptions -# CFLAGS = -g -fexceptions \ - -Wall -W -Wshadow \ - -Wredundant-decls -Wdisabled-optimization -ansi - -# consider: -# -fforce-addr -fmove-all-movables -freduce-all-givs -ftsp-ordering -# -frename-registers -ffast-math -funroll-loops - -# Using the Goto BLAS: -# BLAS = -lgoto -lfrtbegin -lg2c $(XERBLA) -lpthread - -# Using Intel's icc and ifort compilers: -# (does not work for mexFunctions unless you add a mexopts.sh file) -# F77 = ifort -# CC = icc -# CFLAGS = -O3 -xN -vec_report=0 -# CFLAGS = -g -# old (broken): CFLAGS = -ansi -O3 -ip -tpp7 -xW -vec_report0 - -# 64bit: -# F77FLAGS = -O -m64 -# CFLAGS = -O3 -fexceptions -m64 -# BLAS = -lgoto64 -lfrtbegin -lg2c -lpthread $(XERBLA) -# LAPACK = -llapack64 - - -# SUSE Linux 10.1, AMD Opteron, with GOTO Blas -# F77 = gfortran -# BLAS = -lgoto_opteron64 -lgfortran - -# SUSE Linux 10.1, Intel Pentium, with GOTO Blas -# F77 = gfortran -# BLAS = -lgoto -lgfortran - -#------------------------------------------------------------------------------ -# Solaris -#------------------------------------------------------------------------------ - -# 32-bit -# CFLAGS = -KPIC -dalign -xc99=%none -Xc -xlibmieee -xO5 -xlibmil -m32 - -# 64-bit -# CFLAGS = -fast -KPIC -xc99=%none -xlibmieee -xlibmil -m64 -Xc - -# FFLAGS = -fast -KPIC -dalign -xlibmil -m64 - -# The Sun Performance Library includes both LAPACK and the BLAS: -# BLAS = -xlic_lib=sunperf -# LAPACK = - - -#------------------------------------------------------------------------------ -# Compaq Alpha -#------------------------------------------------------------------------------ - -# 64-bit mode only -# CFLAGS = -O2 -std1 -# BLAS = -ldxml -# LAPACK = - -#------------------------------------------------------------------------------ -# Macintosh -#------------------------------------------------------------------------------ - -# CC = gcc -# CFLAGS = -O3 -fno-common -no-cpp-precomp -fexceptions -# LIB = -lstdc++ -# BLAS = -framework Accelerate -# LAPACK = -framework Accelerate - -#------------------------------------------------------------------------------ -# IBM RS 6000 -#------------------------------------------------------------------------------ - -# BLAS = -lessl -# LAPACK = - -# 32-bit mode: -# CFLAGS = -O4 -qipa -qmaxmem=16384 -qproto -# F77FLAGS = -O4 -qipa -qmaxmem=16384 - -# 64-bit mode: -# CFLAGS = -O4 -qipa -qmaxmem=16384 -q64 -qproto -# F77FLAGS = -O4 -qipa -qmaxmem=16384 -q64 -# AR = ar -X64 - -#------------------------------------------------------------------------------ -# SGI IRIX -#------------------------------------------------------------------------------ - -# BLAS = -lscsl -# LAPACK = - -# 32-bit mode -# CFLAGS = -O - -# 64-bit mode (32 bit int's and 64-bit long's): -# CFLAGS = -64 -# F77FLAGS = -64 - -# SGI doesn't have ranlib -# RANLIB = echo - -#------------------------------------------------------------------------------ -# AMD Opteron (64 bit) -#------------------------------------------------------------------------------ - -# BLAS = -lgoto_opteron64 -lg2c -# LAPACK = -llapack_opteron64 - -# SUSE Linux 10.1, AMD Opteron -# F77 = gfortran -# BLAS = -lgoto_opteron64 -lgfortran -# LAPACK = -llapack_opteron64 - -#------------------------------------------------------------------------------ -# remove object files and profile output -#------------------------------------------------------------------------------ - -CLEAN = *.o *.obj *.ln *.bb *.bbg *.da *.tcov *.gcov gmon.out *.bak *.d *.gcda *.gcno diff --git a/ldl/ldl.cpp b/ldl/ldl.cpp deleted file mode 100644 index a9b35c846..000000000 --- a/ldl/ldl.cpp +++ /dev/null @@ -1,507 +0,0 @@ -/* ========================================================================== */ -/* === ldl.c: sparse LDL' factorization and solve package =================== */ -/* ========================================================================== */ - -/* LDL: a simple set of routines for sparse LDL' factorization. These routines - * are not terrifically fast (they do not use dense matrix kernels), but the - * code is very short. The purpose is to illustrate the algorithms in a very - * concise manner, primarily for educational purposes. Although the code is - * very concise, this package is slightly faster than the built-in sparse - * Cholesky factorization in MATLAB 7.0 (chol), when using the same input - * permutation. - * - * The routines compute the LDL' factorization of a real sparse symmetric - * matrix A (or PAP' if a permutation P is supplied), and solve upper - * and lower triangular systems with the resulting L and D factors. If A is - * positive definite then the factorization will be accurate. A can be - * indefinite (with negative values on the diagonal D), but in this case no - * guarantee of accuracy is provided, since no numeric pivoting is performed. - * - * The n-by-n sparse matrix A is in compressed-column form. The nonzero values - * in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with corresponding row - * indices in Ai [Ap [j] ... Ap [j+1]-1]. Ap [0] = 0 is required, and thus - * nz = Ap [n] is the number of nonzeros in A. Ap is an int array of size n+1. - * The int array Ai and the double array Ax are of size nz. This data structure - * is identical to the one used by MATLAB, except for the following - * generalizations. The row indices in each column of A need not be in any - * particular order, although they must be in the range 0 to n-1. Duplicate - * entries can be present; any duplicates are summed. That is, if row index i - * appears twice in a column j, then the value of A (i,j) is the sum of the two - * entries. The data structure used here for the input matrix A is more - * flexible than MATLAB's, which requires sorted columns with no duplicate - * entries. - * - * Only the diagonal and upper triangular part of A (or PAP' if a permutation - * P is provided) is accessed. The lower triangular parts of the matrix A or - * PAP' can be present, but they are ignored. - * - * The optional input permutation is provided as an array P of length n. If - * P [k] = j, the row and column j of A is the kth row and column of PAP'. - * If P is present then the factorization is LDL' = PAP' or L*D*L' = A(P,P) in - * 0-based MATLAB notation. If P is not present (a null pointer) then no - * permutation is performed, and the factorization is LDL' = A. - * - * The lower triangular matrix L is stored in the same compressed-column - * form (an int Lp array of size n+1, an int Li array of size Lp [n], and a - * double array Lx of the same size as Li). It has a unit diagonal, which is - * not stored. The row indices in each column of L are always returned in - * ascending order, with no duplicate entries. This format is compatible with - * MATLAB, except that it would be more convenient for MATLAB to include the - * unit diagonal of L. Doing so here would add additional complexity to the - * code, and is thus omitted in the interest of keeping this code short and - * readable. - * - * The elimination tree is held in the Parent [0..n-1] array. It is normally - * not required by the user, but it is required by ldl_numeric. The diagonal - * matrix D is held as an array D [0..n-1] of size n. - * - * -------------------- - * C-callable routines: - * -------------------- - * - * ldl_symbolic: Given the pattern of A, computes the Lp and Parent arrays - * required by ldl_numeric. Takes time proportional to the number of - * nonzeros in L. Computes the inverse Pinv of P if P is provided. - * Also returns Lnz, the count of nonzeros in each column of L below - * the diagonal (this is not required by ldl_numeric). - * ldl_numeric: Given the pattern and numerical values of A, the Lp array, - * the Parent array, and P and Pinv if applicable, computes the - * pattern and numerical values of L and D. - * ldl_lsolve: Solves Lx=b for a dense vector b. - * ldl_dsolve: Solves Dx=b for a dense vector b. - * ldl_ltsolve: Solves L'x=b for a dense vector b. - * ldl_perm: Computes x=Pb for a dense vector b. - * ldl_permt: Computes x=P'b for a dense vector b. - * ldl_valid_perm: checks the validity of a permutation vector - * ldl_valid_matrix: checks the validity of the sparse matrix A - * - * ---------------------------- - * Limitations of this package: - * ---------------------------- - * - * In the interest of keeping this code simple and readable, ldl_symbolic and - * ldl_numeric assume their inputs are valid. You can check your own inputs - * prior to calling these routines with the ldl_valid_perm and ldl_valid_matrix - * routines. Except for the two ldl_valid_* routines, no routine checks to see - * if the array arguments are present (non-NULL). Like all C routines, no - * routine can determine if the arrays are long enough and don't overlap. - * - * The ldl_numeric does check the numerical factorization, however. It returns - * n if the factorization is successful. If D (k,k) is zero, then k is - * returned, and L is only partially computed. - * - * No pivoting to control fill-in is performed, which is often critical for - * obtaining good performance. I recommend that you compute the permutation P - * using AMD or SYMAMD (approximate minimum degree ordering routines), or an - * appropriate graph-partitioning based ordering. See the ldldemo.m routine for - * an example in MATLAB, and the ldlmain.c stand-alone C program for examples of - * how to find P. Routines for manipulating compressed-column matrices are - * available in UMFPACK. AMD, SYMAMD, UMFPACK, and this LDL package are all - * available at http://www.cise.ufl.edu/research/sparse. - * - * ------------------------- - * Possible simplifications: - * ------------------------- - * - * These routines could be made even simpler with a few additional assumptions. - * If no input permutation were performed, the caller would have to permute the - * matrix first, but the computation of Pinv, and the use of P and Pinv could be - * removed. If only the diagonal and upper triangular part of A or PAP' are - * present, then the tests in the "if (i < k)" statement in ldl_symbolic and - * "if (i <= k)" in ldl_numeric, are always true, and could be removed (i can - * equal k in ldl_symbolic, but then the body of the if statement would - * correctly do no work since Flag [k] == k). If we could assume that no - * duplicate entries are present, then the statement Y [i] += Ax [p] could be - * replaced with Y [i] = Ax [p] in ldl_numeric. - * - * -------------------------- - * Description of the method: - * -------------------------- - * - * LDL computes the symbolic factorization by finding the pattern of L one row - * at a time. It does this based on the following theory. Consider a sparse - * system Lx=b, where L, x, and b, are all sparse, and where L comes from a - * Cholesky (or LDL') factorization. The elimination tree (etree) of L is - * defined as follows. The parent of node j is the smallest k > j such that - * L (k,j) is nonzero. Node j has no parent if column j of L is completely zero - * below the diagonal (j is a root of the etree in this case). The nonzero - * pattern of x is the union of the paths from each node i to the root, for - * each nonzero b (i). To compute the numerical solution to Lx=b, we can - * traverse the columns of L corresponding to nonzero values of x. This - * traversal does not need to be done in the order 0 to n-1. It can be done in - * any "topological" order, such that x (i) is computed before x (j) if i is a - * descendant of j in the elimination tree. - * - * The row-form of the LDL' factorization is shown in the MATLAB function - * ldlrow.m in this LDL package. Note that row k of L is found via a sparse - * triangular solve of L (1:k-1, 1:k-1) \ A (1:k-1, k), to use 1-based MATLAB - * notation. Thus, we can start with the nonzero pattern of the kth column of - * A (above the diagonal), follow the paths up to the root of the etree of the - * (k-1)-by-(k-1) leading submatrix of L, and obtain the pattern of the kth row - * of L. Note that we only need the leading (k-1)-by-(k-1) submatrix of L to - * do this. The elimination tree can be constructed as we go. - * - * The symbolic factorization does the same thing, except that it discards the - * pattern of L as it is computed. It simply counts the number of nonzeros in - * each column of L and then constructs the Lp index array when it's done. The - * symbolic factorization does not need to do this in topological order. - * Compare ldl_symbolic with the first part of ldl_numeric, and note that the - * while (len > 0) loop is not present in ldl_symbolic. - * - * LDL Version 1.3, Copyright (c) 2006 by Timothy A Davis, - * University of Florida. All Rights Reserved. Developed while on sabbatical - * at Stanford University and Lawrence Berkeley National Laboratory. Refer to - * the README file for the License. Available at - * http://www.cise.ufl.edu/research/sparse. - */ - -#include "ldl.h" - -/* ========================================================================== */ -/* === ldl_symbolic ========================================================= */ -/* ========================================================================== */ - -/* The input to this routine is a sparse matrix A, stored in column form, and - * an optional permutation P. The output is the elimination tree - * and the number of nonzeros in each column of L. Parent [i] = k if k is the - * parent of i in the tree. The Parent array is required by ldl_numeric. - * Lnz [k] gives the number of nonzeros in the kth column of L, excluding the - * diagonal. - * - * One workspace vector (Flag) of size n is required. - * - * If P is NULL, then it is ignored. The factorization will be LDL' = A. - * Pinv is not computed. In this case, neither P nor Pinv are required by - * ldl_numeric. - * - * If P is not NULL, then it is assumed to be a valid permutation. If - * row and column j of A is the kth pivot, the P [k] = j. The factorization - * will be LDL' = PAP', or A (p,p) in MATLAB notation. The inverse permutation - * Pinv is computed, where Pinv [j] = k if P [k] = j. In this case, both P - * and Pinv are required as inputs to ldl_numeric. - * - * The floating-point operation count of the subsequent call to ldl_numeric - * is not returned, but could be computed after ldl_symbolic is done. It is - * the sum of (Lnz [k]) * (Lnz [k] + 2) for k = 0 to n-1. - */ - -void LDL_symbolic -( - LDL_int n, /* A and L are n-by-n, where n >= 0 */ - LDL_int Ap [ ], /* input of size n+1, not modified */ - LDL_int Ai [ ], /* input of size nz=Ap[n], not modified */ - LDL_int Lp [ ], /* output of size n+1, not defined on input */ - LDL_int Parent [ ], /* output of size n, not defined on input */ - LDL_int Lnz [ ], /* output of size n, not defined on input */ - LDL_int Flag [ ], /* workspace of size n, not defn. on input or output */ - LDL_int P [ ], /* optional input of size n */ - LDL_int Pinv [ ] /* optional output of size n (used if P is not NULL) */ -) -{ - LDL_int i, k, p, kk, p2 ; - if (P) - { - /* If P is present then compute Pinv, the inverse of P */ - for (k = 0 ; k < n ; k++) - { - Pinv [P [k]] = k ; - } - } - for (k = 0 ; k < n ; k++) - { - /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ - Parent [k] = -1 ; /* parent of k is not yet known */ - Flag [k] = k ; /* mark node k as visited */ - Lnz [k] = 0 ; /* count of nonzeros in column k of L */ - kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */ - p2 = Ap [kk+1] ; - for (p = Ap [kk] ; p < p2 ; p++) - { - /* A (i,k) is nonzero (original or permuted A) */ - i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; - if (i < k) - { - /* follow path from i to root of etree, stop at flagged node */ - for ( ; Flag [i] != k ; i = Parent [i]) - { - /* find parent of i if not yet determined */ - if (Parent [i] == -1) Parent [i] = k ; - Lnz [i]++ ; /* L (k,i) is nonzero */ - Flag [i] = k ; /* mark i as visited */ - } - } - } - } - /* construct Lp index array from Lnz column counts */ - Lp [0] = 0 ; - for (k = 0 ; k < n ; k++) - { - Lp [k+1] = Lp [k] + Lnz [k] ; - } -} - - -/* ========================================================================== */ -/* === ldl_numeric ========================================================== */ -/* ========================================================================== */ - -/* Given a sparse matrix A (the arguments n, Ap, Ai, and Ax) and its symbolic - * analysis (Lp and Parent, and optionally P and Pinv), compute the numeric LDL' - * factorization of A or PAP'. The outputs of this routine are arguments Li, - * Lx, and D. It also requires three size-n workspaces (Y, Pattern, and Flag). - */ - -LDL_int LDL_numeric /* returns n if successful, k if D (k,k) is zero */ -( - LDL_int n, /* A and L are n-by-n, where n >= 0 */ - LDL_int Ap [ ], /* input of size n+1, not modified */ - LDL_int Ai [ ], /* input of size nz=Ap[n], not modified */ - double Ax [ ], /* input of size nz=Ap[n], not modified */ - LDL_int Lp [ ], /* input of size n+1, not modified */ - LDL_int Parent [ ], /* input of size n, not modified */ - LDL_int Lnz [ ], /* output of size n, not defn. on input */ - LDL_int Li [ ], /* output of size lnz=Lp[n], not defined on input */ - double Lx [ ], /* output of size lnz=Lp[n], not defined on input */ - double D [ ], /* output of size n, not defined on input */ - double Y [ ], /* workspace of size n, not defn. on input or output */ - LDL_int Pattern [ ],/* workspace of size n, not defn. on input or output */ - LDL_int Flag [ ], /* workspace of size n, not defn. on input or output */ - LDL_int P [ ], /* optional input of size n */ - LDL_int Pinv [ ] /* optional input of size n */ -) -{ - double yi, l_ki ; - LDL_int i, k, p, kk, p2, len, top ; - for (k = 0 ; k < n ; k++) - { - /* compute nonzero Pattern of kth row of L, in topological order */ - Y [k] = 0.0 ; /* Y(0:k) is now all zero */ - top = n ; /* stack for pattern is empty */ - Flag [k] = k ; /* mark node k as visited */ - Lnz [k] = 0 ; /* count of nonzeros in column k of L */ - kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */ - p2 = Ap [kk+1] ; - for (p = Ap [kk] ; p < p2 ; p++) - { - i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; /* get A(i,k) */ - if (i <= k) - { - Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */ - for (len = 0 ; Flag [i] != k ; i = Parent [i]) - { - Pattern [len++] = i ; /* L(k,i) is nonzero */ - Flag [i] = k ; /* mark i as visited */ - } - while (len > 0) Pattern [--top] = Pattern [--len] ; - } - } - /* compute numerical values kth row of L (a sparse triangular solve) */ - D [k] = Y [k] ; /* get D(k,k) and clear Y(k) */ - Y [k] = 0.0 ; - for ( ; top < n ; top++) - { - i = Pattern [top] ; /* Pattern [top:n-1] is pattern of L(:,k) */ - yi = Y [i] ; /* get and clear Y(i) */ - Y [i] = 0.0 ; - p2 = Lp [i] + Lnz [i] ; - for (p = Lp [i] ; p < p2 ; p++) - { - Y [Li [p]] -= Lx [p] * yi ; - } - l_ki = yi / D [i] ; /* the nonzero entry L(k,i) */ - D [k] -= l_ki * yi ; - Li [p] = k ; /* store L(k,i) in column form of L */ - Lx [p] = l_ki ; - Lnz [i]++ ; /* increment count of nonzeros in col i */ - } - if (D [k] == 0.0) return (k) ; /* failure, D(k,k) is zero */ - } - return (n) ; /* success, diagonal of D is all nonzero */ -} - - -/* ========================================================================== */ -/* === ldl_lsolve: solve Lx=b ============================================== */ -/* ========================================================================== */ - -void LDL_lsolve -( - LDL_int n, /* L is n-by-n, where n >= 0 */ - double X [ ], /* size n. right-hand-side on input, soln. on output */ - LDL_int Lp [ ], /* input of size n+1, not modified */ - LDL_int Li [ ], /* input of size lnz=Lp[n], not modified */ - double Lx [ ] /* input of size lnz=Lp[n], not modified */ -) -{ - LDL_int j, p, p2 ; - for (j = 0 ; j < n ; j++) - { - p2 = Lp [j+1] ; - for (p = Lp [j] ; p < p2 ; p++) - { - X [Li [p]] -= Lx [p] * X [j] ; - } - } -} - - -/* ========================================================================== */ -/* === ldl_dsolve: solve Dx=b ============================================== */ -/* ========================================================================== */ - -void LDL_dsolve -( - LDL_int n, /* D is n-by-n, where n >= 0 */ - double X [ ], /* size n. right-hand-side on input, soln. on output */ - double D [ ] /* input of size n, not modified */ -) -{ - LDL_int j ; - for (j = 0 ; j < n ; j++) - { - X [j] /= D [j] ; - } -} - - -/* ========================================================================== */ -/* === ldl_ltsolve: solve L'x=b ============================================ */ -/* ========================================================================== */ - -void LDL_ltsolve -( - LDL_int n, /* L is n-by-n, where n >= 0 */ - double X [ ], /* size n. right-hand-side on input, soln. on output */ - LDL_int Lp [ ], /* input of size n+1, not modified */ - LDL_int Li [ ], /* input of size lnz=Lp[n], not modified */ - double Lx [ ] /* input of size lnz=Lp[n], not modified */ -) -{ - int j, p, p2 ; - for (j = n-1 ; j >= 0 ; j--) - { - p2 = Lp [j+1] ; - for (p = Lp [j] ; p < p2 ; p++) - { - X [j] -= Lx [p] * X [Li [p]] ; - } - } -} - - -/* ========================================================================== */ -/* === ldl_perm: permute a vector, x=Pb ===================================== */ -/* ========================================================================== */ - -void LDL_perm -( - LDL_int n, /* size of X, B, and P */ - double X [ ], /* output of size n. */ - double B [ ], /* input of size n. */ - LDL_int P [ ] /* input permutation array of size n. */ -) -{ - LDL_int j ; - for (j = 0 ; j < n ; j++) - { - X [j] = B [P [j]] ; - } -} - - -/* ========================================================================== */ -/* === ldl_permt: permute a vector, x=P'b =================================== */ -/* ========================================================================== */ - -void LDL_permt -( - LDL_int n, /* size of X, B, and P */ - double X [ ], /* output of size n. */ - double B [ ], /* input of size n. */ - LDL_int P [ ] /* input permutation array of size n. */ -) -{ - LDL_int j ; - for (j = 0 ; j < n ; j++) - { - X [P [j]] = B [j] ; - } -} - - -/* ========================================================================== */ -/* === ldl_valid_perm: check if a permutation vector is valid =============== */ -/* ========================================================================== */ - -LDL_int LDL_valid_perm /* returns 1 if valid, 0 otherwise */ -( - LDL_int n, - LDL_int P [ ], /* input of size n, a permutation of 0:n-1 */ - LDL_int Flag [ ] /* workspace of size n */ -) -{ - LDL_int j, k ; - if (n < 0 || !Flag) - { - return (0) ; /* n must be >= 0, and Flag must be present */ - } - if (!P) - { - return (1) ; /* If NULL, P is assumed to be the identity perm. */ - } - for (j = 0 ; j < n ; j++) - { - Flag [j] = 0 ; /* clear the Flag array */ - } - for (k = 0 ; k < n ; k++) - { - j = P [k] ; - if (j < 0 || j >= n || Flag [j] != 0) - { - return (0) ; /* P is not valid */ - } - Flag [j] = 1 ; - } - return (1) ; /* P is valid */ -} - - -/* ========================================================================== */ -/* === ldl_valid_matrix: check if a sparse matrix is valid ================== */ -/* ========================================================================== */ - -/* This routine checks to see if a sparse matrix A is valid for input to - * ldl_symbolic and ldl_numeric. It returns 1 if the matrix is valid, 0 - * otherwise. A is in sparse column form. The numerical values in column j - * are stored in Ax [Ap [j] ... Ap [j+1]-1], with row indices in - * Ai [Ap [j] ... Ap [j+1]-1]. The Ax array is not checked. - */ - -LDL_int LDL_valid_matrix -( - LDL_int n, - LDL_int Ap [ ], - LDL_int Ai [ ] -) -{ - LDL_int j, p ; - if (n < 0 || !Ap || !Ai || Ap [0] != 0) - { - return (0) ; /* n must be >= 0, and Ap and Ai must be present */ - } - for (j = 0 ; j < n ; j++) - { - if (Ap [j] > Ap [j+1]) - { - return (0) ; /* Ap must be monotonically nondecreasing */ - } - } - for (p = 0 ; p < Ap [n] ; p++) - { - if (Ai [p] < 0 || Ai [p] >= n) - { - return (0) ; /* row indices must be in the range 0 to n-1 */ - } - } - return (1) ; /* matrix is valid */ -} diff --git a/ldl/ldl.h b/ldl/ldl.h deleted file mode 100644 index 2ef71e07f..000000000 --- a/ldl/ldl.h +++ /dev/null @@ -1,104 +0,0 @@ -/* ========================================================================== */ -/* === ldl.h: include file for the LDL package ============================= */ -/* ========================================================================== */ - -/* LDL Copyright (c) Timothy A Davis, - * University of Florida. All Rights Reserved. See README for the License. - */ - -#include "UFconfig.h" - -//#ifdef LDL_LONG -//#define LDL_int UF_long -//#define LDL_ID UF_long_id -// -//#define LDL_symbolic ldl_l_symbolic -//#define LDL_numeric ldl_l_numeric -//#define LDL_lsolve ldl_l_lsolve -//#define LDL_dsolve ldl_l_dsolve -//#define LDL_ltsolve ldl_l_ltsolve -//#define LDL_perm ldl_l_perm -//#define LDL_permt ldl_l_permt -//#define LDL_valid_perm ldl_l_valid_perm -//#define LDL_valid_matrix ldl_l_valid_matrix -// -//#else -#define LDL_int int -#define LDL_ID "%d" - -//#define LDL_symbolic ldl_symbolic -//#define LDL_numeric ldl_numeric -//#define LDL_lsolve ldl_lsolve -//#define LDL_dsolve ldl_dsolve -//#define LDL_ltsolve ldl_ltsolve -//#define LDL_perm ldl_perm -//#define LDL_permt ldl_permt -//#define LDL_valid_perm ldl_valid_perm -//#define LDL_valid_matrix ldl_valid_matrix - -//#endif - -/* ========================================================================== */ -/* === int version ========================================================== */ -/* ========================================================================== */ - -void LDL_symbolic (int n, int Ap [ ], int Ai [ ], int Lp [ ], - int Parent [ ], int Lnz [ ], int Flag [ ], int P [ ], int Pinv [ ]) ; - -int LDL_numeric (int n, int Ap [ ], int Ai [ ], double Ax [ ], - int Lp [ ], int Parent [ ], int Lnz [ ], int Li [ ], double Lx [ ], - double D [ ], double Y [ ], int Pattern [ ], int Flag [ ], - int P [ ], int Pinv [ ]) ; - -void LDL_lsolve (int n, double X [ ], int Lp [ ], int Li [ ], - double Lx [ ]) ; - -void LDL_dsolve (int n, double X [ ], double D [ ]) ; - -void LDL_ltsolve (int n, double X [ ], int Lp [ ], int Li [ ], - double Lx [ ]) ; - -void LDL_perm (int n, double X [ ], double B [ ], int P [ ]) ; -void LDL_permt (int n, double X [ ], double B [ ], int P [ ]) ; - -int LDL_valid_perm (int n, int P [ ], int Flag [ ]) ; -int LDL_valid_matrix ( int n, int Ap [ ], int Ai [ ]) ; - -/* ========================================================================== */ -/* === long version ========================================================= */ -/* ========================================================================== */ - -//void ldl_l_symbolic (UF_long n, UF_long Ap [ ], UF_long Ai [ ], UF_long Lp [ ], -// UF_long Parent [ ], UF_long Lnz [ ], UF_long Flag [ ], UF_long P [ ], -// UF_long Pinv [ ]) ; -// -//UF_long ldl_l_numeric (UF_long n, UF_long Ap [ ], UF_long Ai [ ], double Ax [ ], -// UF_long Lp [ ], UF_long Parent [ ], UF_long Lnz [ ], UF_long Li [ ], -// double Lx [ ], double D [ ], double Y [ ], UF_long Pattern [ ], -// UF_long Flag [ ], UF_long P [ ], UF_long Pinv [ ]) ; -// -//void ldl_l_lsolve (UF_long n, double X [ ], UF_long Lp [ ], UF_long Li [ ], -// double Lx [ ]) ; -// -//void ldl_l_dsolve (UF_long n, double X [ ], double D [ ]) ; -// -//void ldl_l_ltsolve (UF_long n, double X [ ], UF_long Lp [ ], UF_long Li [ ], -// double Lx [ ]) ; -// -//void ldl_l_perm (UF_long n, double X [ ], double B [ ], UF_long P [ ]) ; -//void ldl_l_permt (UF_long n, double X [ ], double B [ ], UF_long P [ ]) ; -// -//UF_long ldl_l_valid_perm (UF_long n, UF_long P [ ], UF_long Flag [ ]) ; -//UF_long ldl_l_valid_matrix ( UF_long n, UF_long Ap [ ], UF_long Ai [ ]) ; - -/* ========================================================================== */ -/* === LDL version ========================================================== */ -/* ========================================================================== */ - -#define LDL_DATE "Nov 1, 2007" -#define LDL_VERSION_CODE(main,sub) ((main) * 1000 + (sub)) -#define LDL_MAIN_VERSION 2 -#define LDL_SUB_VERSION 0 -#define LDL_SUBSUB_VERSION 1 -#define LDL_VERSION LDL_VERSION_CODE(LDL_MAIN_VERSION,LDL_SUB_VERSION) - diff --git a/ldl/ldlsimple.cpp b/ldl/ldlsimple.cpp deleted file mode 100644 index 87a355565..000000000 --- a/ldl/ldlsimple.cpp +++ /dev/null @@ -1,69 +0,0 @@ -/* ========================================================================== */ -/* === ldlsimple.c: a simple LDL main program ============================== */ -/* ========================================================================== */ - -/* LDLSIMPLE: this is a very simple main program that illustrates the basic - * usage of the LDL routines. The output of this program is in ldlsimple.out. - * This program factorizes the matrix - * - * A =[ ... - * 1.7 0 0 0 0 0 0 0 .13 0 - * 0 1. 0 0 .02 0 0 0 0 .01 - * 0 0 1.5 0 0 0 0 0 0 0 - * 0 0 0 1.1 0 0 0 0 0 0 - * 0 .02 0 0 2.6 0 .16 .09 .52 .53 - * 0 0 0 0 0 1.2 0 0 0 0 - * 0 0 0 0 .16 0 1.3 0 0 .56 - * 0 0 0 0 .09 0 0 1.6 .11 0 - * .13 0 0 0 .52 0 0 .11 1.4 0 - * 0 .01 0 0 .53 0 .56 0 0 3.1 ] ; - * - * and then solves a system Ax=b whose true solution is - * x = [.1 .2 .3 .4 .5 .6 .7 .8 .9 1]' ; - * - * Note that Li and Lx are statically allocated, with length 13. This is just - * enough to hold the factor L, but normally this size is not known until after - * ldl_symbolic has analyzed the matrix. The size of Li and Lx must be greater - * than or equal to lnz = Lp [N], which is 13 for this matrix L. - * - * LDL Version 1.3, Copyright (c) 2006 by Timothy A Davis, - * University of Florida. All Rights Reserved. See README for the License. - */ - -#include -#include "ldl.h" -#define N 10 /* A is 10-by-10 */ -#define ANZ 19 /* # of nonzeros on diagonal and upper triangular part of A */ -#define LNZ 13 /* # of nonzeros below the diagonal of L */ - -int main (void) -{ - /* only the upper triangular part of A is required */ - int Ap [N+1] = {0, 1, 2, 3, 4, 6, 7, 9, 11, 15, ANZ}, - Ai [ANZ] = {0, 1, 2, 3, 1,4, 5, 4,6, 4,7, 0,4,7,8, 1,4,6,9 } ; - double Ax [ANZ] = {1.7, 1., 1.5, 1.1, .02,2.6, 1.2, .16,1.3, .09,1.6, - .13,.52,.11,1.4, .01,.53,.56,3.1}, - b [N] = {.287, .22, .45, .44, 2.486, .72, 1.55, 1.424, 1.621, 3.759}; - double Lx [LNZ], D [N], Y [N] ; - int Li [LNZ], Lp [N+1], Parent [N], Lnz [N], Flag [N], Pattern [N], d, i ; - - /* factorize A into LDL' (P and Pinv not used) */ - LDL_symbolic (N, Ap, Ai, Lp, Parent, Lnz, Flag, NULL, NULL) ; - printf ("Nonzeros in L, excluding diagonal: %d\n", Lp [N]) ; - d = LDL_numeric (N, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Pattern, - Flag, NULL, NULL) ; - - if (d == N) - { - /* solve Ax=b, overwriting b with the solution x */ - LDL_lsolve (N, b, Lp, Li, Lx) ; - LDL_dsolve (N, b, D) ; - LDL_ltsolve (N, b, Lp, Li, Lx) ; - for (i = 0 ; i < N ; i++) printf ("x [%d] = %g\n", i, b [i]) ; - } - else - { - printf ("LDL_numeric failed, D (%d,%d) is zero\n", d, d) ; - } - return (0) ; -} diff --git a/linear/Makefile.am b/linear/Makefile.am index f3432190e..13afc1cb6 100644 --- a/linear/Makefile.am +++ b/linear/Makefile.am @@ -50,8 +50,11 @@ AM_CXXFLAGS = TESTS = $(check_PROGRAMS) AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) LDADD = liblinear.la ../inference/libinference.la ../base/libbase.la -LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la ../ldl/libldl.la +LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la AM_DEFAULT_SOURCE_EXT = .cpp +if USE_LDL +LDADD += libldl.la +endif # rule to run an executable %.run: % $(LDADD) diff --git a/matlab/example/alex01_simple1.m b/matlab/example/alex01_simple1.m index 809e2420f..8f3eb5706 100644 --- a/matlab/example/alex01_simple1.m +++ b/matlab/example/alex01_simple1.m @@ -1,4 +1,6 @@ % Script to perform SQP on a simple example from the SQP tutorial +% This makes use of LDL solving of a full quadratic programming +% problem % % Problem: % min(x) f(x) = (x2-2)^2 - x1^2 @@ -32,7 +34,7 @@ for i=1:maxIt dcx = [8*x1; 2*x2]; dL = dfx - lam * dcx; - % update the hessian (BFGS) + % update the hessian (BFGS) - note: this does not use the full lagrangian if (i>1) Bis = Bi*s; y = dfx - prev_dfx; diff --git a/nonlinear/ConstraintOptimizer.cpp b/nonlinear/ConstraintOptimizer.cpp index 85f007c0c..67c36c2b9 100644 --- a/nonlinear/ConstraintOptimizer.cpp +++ b/nonlinear/ConstraintOptimizer.cpp @@ -3,6 +3,8 @@ * @author Alex Cunningham */ +/** IMPORTANT NOTE: this file is only compiled when LDL is available */ + #include using namespace std; diff --git a/nonlinear/ConstraintOptimizer.h b/nonlinear/ConstraintOptimizer.h index 68612d333..b7051223c 100644 --- a/nonlinear/ConstraintOptimizer.h +++ b/nonlinear/ConstraintOptimizer.h @@ -4,6 +4,13 @@ * @author Alex Cunningham */ +/** + * IMPORTANT NOTE: This is an EXPERIMENTAL system that is not ready for use! + * DO NOT USE if you actually wanted to accomplish something + * + * REQUIRES --enable-ldl flag to be set for this class to work, not compiled otherwise + */ + #pragma once #include @@ -76,4 +83,3 @@ namespace gtsam { } // \namespace gtsam - diff --git a/nonlinear/Makefile.am b/nonlinear/Makefile.am index 37949408d..9d02f24f6 100644 --- a/nonlinear/Makefile.am +++ b/nonlinear/Makefile.am @@ -26,8 +26,10 @@ headers += NonlinearEquality.h check_PROGRAMS += tests/testNonlinearConstraint # SQP +if USE_LDL sources += ConstraintOptimizer.cpp check_PROGRAMS += tests/testConstraintOptimizer +endif #---------------------------------------------------------------------------------------------------- # Create a libtool library that is not installed @@ -47,8 +49,11 @@ AM_CXXFLAGS = TESTS = $(check_PROGRAMS) AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) LDADD = libnonlinear.la ../linear/liblinear.la ../inference/libinference.la ../base/libbase.la -LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la ../ldl/libldl.la +LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la AM_DEFAULT_SOURCE_EXT = .cpp +if USE_LDL +LDADD += libldl.la +endif # rule to run an executable %.run: % $(LDADD) diff --git a/nonlinear/tests/testConstraintOptimizer.cpp b/nonlinear/tests/testConstraintOptimizer.cpp index de29cf0b4..3eabab300 100644 --- a/nonlinear/tests/testConstraintOptimizer.cpp +++ b/nonlinear/tests/testConstraintOptimizer.cpp @@ -4,6 +4,8 @@ * @author Alex Cunningham */ +/** IMPORTANT NOTE: this file is only compiled when LDL is available */ + #include #include diff --git a/slam/Makefile.am b/slam/Makefile.am index 76fce1340..5e7c2f1c7 100644 --- a/slam/Makefile.am +++ b/slam/Makefile.am @@ -67,7 +67,10 @@ TESTS = $(check_PROGRAMS) AM_DEFAULT_SOURCE_EXT = .cpp AM_LDFLAGS = $(BOOST_LDFLAGS) $(boost_serialization) LDADD = libslam.la ../geometry/libgeometry.la ../nonlinear/libnonlinear.la ../linear/liblinear.la ../inference/libinference.la ../base/libbase.la -LDADD += ../CppUnitLite/libCppUnitLite.a ../ldl/libldl.la ../colamd/libcolamd.la +LDADD += ../CppUnitLite/libCppUnitLite.a ../colamd/libcolamd.la +if USE_LDL +LDADD += libldl.la +endif if USE_LAPACK LDADD += ../spqr_mini/libspqr_mini.la diff --git a/tests/Makefile.am b/tests/Makefile.am index 973a80b73..94e95d5c6 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -3,7 +3,7 @@ # More elaborate unit tests that test functionality with slam examples #---------------------------------------------------------------------------------------------------- -check_PROGRAMS = testBayesNetPreconditioner testConstraintOptimizer +check_PROGRAMS = testBayesNetPreconditioner check_PROGRAMS += testGaussianBayesNet testGaussianFactor testGaussianFactorGraph check_PROGRAMS += testGaussianISAM testGraph check_PROGRAMS += testInference testIterative testGaussianJunctionTree @@ -11,6 +11,10 @@ check_PROGRAMS += testNonlinearEquality testNonlinearFactor testNonlinearFactorG check_PROGRAMS += testNonlinearOptimizer testSQP testSubgraphPreconditioner check_PROGRAMS += testSymbolicBayesNet testSymbolicFactorGraph testTupleConfig +if USE_LDL +check_PROGRAMS = testConstraintOptimizer +endif + # Timing tests noinst_PROGRAMS = timeGaussianFactorGraph