From 25bd1c840f38d0fd2fed2c469ae75e9d4c083155 Mon Sep 17 00:00:00 2001 From: Alex Cunningham Date: Thu, 22 Apr 2010 22:17:08 +0000 Subject: [PATCH] Added Tim Davis' LDL library to use for solving quadratic programming problem. Currently, code compiles and executes some sample code in a test for Matrix. Also fixed some small issues with colamd. --- Makefile.am | 13 +- colamd/Makefile | 5 +- cpp/Makefile.am | 4 +- cpp/Matrix.cpp | 56 ++++ cpp/Matrix.h | 3 + cpp/testMatrix.cpp | 30 ++ ldl/Makefile | 39 +++ ldl/README.txt | 136 +++++++++ ldl/UFconfig.h | 118 ++++++++ ldl/UFconfig.mk | 354 ++++++++++++++++++++++ ldl/ldl.cpp | 507 ++++++++++++++++++++++++++++++++ ldl/ldl.h | 104 +++++++ ldl/ldlsimple.cpp | 69 +++++ matlab/example/alex01_simple1.m | 11 +- 14 files changed, 1437 insertions(+), 12 deletions(-) create mode 100644 ldl/Makefile create mode 100644 ldl/README.txt create mode 100644 ldl/UFconfig.h create mode 100644 ldl/UFconfig.mk create mode 100644 ldl/ldl.cpp create mode 100644 ldl/ldl.h create mode 100644 ldl/ldlsimple.cpp diff --git a/Makefile.am b/Makefile.am index dc5bce42c..ea66ff12a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,7 +3,7 @@ ACLOCAL_AMFLAGS = -I m4 # make automake install some standard but missing files AUTOMAKE_OPTIONS = foreign -SUBDIRS = colamd CppUnitLite cpp wrap +SUBDIRS = colamd CppUnitLite ldl cpp wrap # install the matlab toolbox install-exec-hook: @@ -15,7 +15,11 @@ install-exec-hook: install -d $(prefix)/include/colamd && \ cp -f colamd/colamd.h $(prefix)/include/colamd/ && \ cp -f colamd/UFconfig.h $(prefix)/include/colamd/ && \ - cp -f colamd/libcolamd.a $(prefix)/lib/ + cp -f colamd/libcolamd.a $(prefix)/lib/ && \ + install -d $(prefix)/include/ldl && \ + cp -f ldl/ldl.h $(prefix)/include/ldl/ && \ + cp -f ldl/UFconfig.h $(prefix)/include/ldl/ && \ + cp -f ldl/libldl.a $(prefix)/lib/ dist-hook: mkdir $(distdir)/CppUnitLite @@ -27,5 +31,10 @@ dist-hook: cp -p $(srcdir)/colamd/*.h $(distdir)/colamd cp -p $(srcdir)/colamd/*.c $(distdir)/colamd cp -p $(srcdir)/colamd/*.cpp $(distdir)/colamd + mkdir $(distdir)/ldl + cp -p $(srcdir)/ldl/Makefile $(distdir)/ldl + cp -p $(srcdir)/ldl/*.h $(distdir)/ldl + cp -p $(srcdir)/ldl/*.c $(distdir)/ldl + cp -p $(srcdir)/ldl/*.cpp $(distdir)/ldl mkdir $(distdir)/matlab cp -p $(srcdir)/matlab/*.m $(distdir)/matlab diff --git a/colamd/Makefile b/colamd/Makefile index 319373ac9..e91b5567e 100755 --- a/colamd/Makefile +++ b/colamd/Makefile @@ -6,7 +6,6 @@ all: libcolamd.a CC=gcc CXX=gcc -CPP=gcc CXXFLAGS += -O2 CXXFLAGS += -fPIC @@ -20,8 +19,8 @@ library = libcolamd.a # This should probably be fixed, but will work for 64 bit machines now $(library): $(objects) - $(CPP) $(CXXFLAGS) -c -o colamd.o colamd.c - $(CPP) $(CXXFLAGS) -c -o colamd_global.o colamd_global.c + $(CXX) $(CXXFLAGS) -c -o colamd.o colamd.c + $(CXX) $(CXXFLAGS) -c -o colamd_global.o colamd_global.c ar crsv $@ $(objects) ranlib $(library) diff --git a/cpp/Makefile.am b/cpp/Makefile.am index aa6d0d172..f8f47fdd7 100644 --- a/cpp/Makefile.am +++ b/cpp/Makefile.am @@ -283,7 +283,7 @@ AM_CXXFLAGS = -I$(boost) -fPIC lib_LTLIBRARIES = libgtsam.la libgtsam_la_SOURCES = $(sources) libgtsam_la_CPPFLAGS = $(AM_CXXFLAGS) -libgtsam_la_LDFLAGS = -version-info $(version) -L../colamd -lcolamd #-lboost_serialization-mt +libgtsam_la_LDFLAGS = -version-info $(version) ../ldl/libldl.a -L../colamd -lcolamd # -lboost_serialization-mt # enable debug if --enable-debug is set in configure if DEBUG @@ -294,7 +294,7 @@ endif include_HEADERS = $(headers) AM_CXXFLAGS += -I.. -AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serialization) #-L. -lgtsam +AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite ../ldl/libldl.a $(BOOST_LDFLAGS) $(boost_serialization) # adding cblas implementation - split into a default linux version using the # autotools script, and a mac version that is hardcoded diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp index 6261952c4..ef618671a 100644 --- a/cpp/Matrix.cpp +++ b/cpp/Matrix.cpp @@ -28,6 +28,8 @@ #include #include +#include + #include "Matrix.h" #include "Vector.h" #include "svdcmp.h" @@ -891,6 +893,60 @@ Matrix RtR(const Matrix &A) return trans(LLt(A)); } +/* ************************************************************************* */ +Vector solve_ldl(const Matrix& M, const Vector& rhs) { + // execute test in here + + // create a matrix (from ldlsimple.c example) + const int N = 10, // size of the matrix + ANZ = 19, // # of nonzeros on diagonal and upper triangular part of A + LNZ = 13; // # of nonzeros below the diagonal of L + Matrix A = Matrix_(N, N, + 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, .13, 0.0, + 0.0, 1., 0.0, 0.0, .02, 0.0, 0.0, 0.0, 0.0, .01, + 0.0, 0.0, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, .02, 0.0, 0.0, 2.6, 0.0, .16, .09, .52, .53, + 0.0, 0.0, 0.0, 0.0, 0.0, 1.2, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, .16, 0.0, 1.3, 0.0, 0.0, .56, + 0.0, 0.0, 0.0, 0.0, .09, 0.0, 0.0, 1.6, .11, 0.0, + .13, 0.0, 0.0, 0.0, .52, 0.0, 0.0, .11, 1.4, 0.0, + 0.0, .01, 0.0, 0.0, .53, 0.0, .56, 0.0, 0.0, 3.1); + + // convert to the right format + + /* 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); + + cout << "Nonzeros in L, excluding diagonal: " << Lp [N] << endl; + 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++) + cout << "x [" << i << "] = " << b[i] << endl; + } + + // copy solution out + Vector result = Vector_(N, b); + + return result; +} + /* * This is not ultra efficient, but not terrible, either. */ diff --git a/cpp/Matrix.h b/cpp/Matrix.h index 3bba01839..ad4aeeae5 100644 --- a/cpp/Matrix.h +++ b/cpp/Matrix.h @@ -359,6 +359,9 @@ Matrix RtR(const Matrix& A); /** Return the inverse of a S.P.D. matrix. Inversion is done via Cholesky decomposition. */ Matrix cholesky_inverse(const Matrix &A); +/** Solve Ax=b with S.P.D. matrix using Davis' LDL code */ +Vector solve_ldl(const Matrix& A, const Vector& b); + /** * Numerical exponential map, naive approach, not industrial strength !!! * @param A matrix to exponentiate diff --git a/cpp/testMatrix.cpp b/cpp/testMatrix.cpp index 7902183cd..e889bc9be 100644 --- a/cpp/testMatrix.cpp +++ b/cpp/testMatrix.cpp @@ -7,6 +7,7 @@ #include #include +//#include #include #include #include @@ -935,6 +936,35 @@ TEST( matrix, transposeMultiplyAdd ) CHECK(assert_equal(expected, x)); } +/* ************************************************************************* */ +TEST( matrix, LDL_factorization ) { + + // run demo inside Matrix.cpp code + + // create a matrix (from ldlsimple.c example) + Matrix A = Matrix_(10, 10, + 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, .13, 0.0, + 0.0, 1., 0.0, 0.0, .02, 0.0, 0.0, 0.0, 0.0, .01, + 0.0, 0.0, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, .02, 0.0, 0.0, 2.6, 0.0, .16, .09, .52, .53, + 0.0, 0.0, 0.0, 0.0, 0.0, 1.2, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, .16, 0.0, 1.3, 0.0, 0.0, .56, + 0.0, 0.0, 0.0, 0.0, .09, 0.0, 0.0, 1.6, .11, 0.0, + .13, 0.0, 0.0, 0.0, .52, 0.0, 0.0, .11, 1.4, 0.0, + 0.0, .01, 0.0, 0.0, .53, 0.0, .56, 0.0, 0.0, 3.1); + Vector b = Vector_(10, .287, .22, .45, .44, 2.486, .72, 1.55, 1.424, 1.621, 3.759); + + // perform LDL solving + Vector actual = solve_ldl(A, b); + + // check solution + Vector expected = Vector_(10, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1.0); + + CHECK(assert_equal(expected, actual)); + +} + /* ************************************************************************* */ int main() { TestResult tr; diff --git a/ldl/Makefile b/ldl/Makefile new file mode 100644 index 000000000..9f2f0384e --- /dev/null +++ b/ldl/Makefile @@ -0,0 +1,39 @@ +#------------------------------------------------------------------------------- +# Makefile for the LDL library +#------------------------------------------------------------------------------- + +default: all + +include ./UFconfig.mk + +CXXFLAGS += -O2 +CXXFLAGS += -fPIC -I./ + +CXX = g++ + +all: demo + +#------------------------------------------------------------------------------- +# the ldl library: +#------------------------------------------------------------------------------- + +libldl.a: ldl.cpp ldl.h + $(CXX) $(CXXFLAGS) -c ldl.cpp -o ldl.o + $(AR) libldl.a ldl.o + $(RANLIB) libldl.a + +demo: ldlsimple.cpp libldl.a + $(CXX) ldlsimple.cpp -o ldldemo -L./ -lldl + +clean: + $(RM) $(CLEAN) libldl.a ldldemo + +check: + echo 'no check for ldl' + +distdir: + +install: all + +distclean: clean + diff --git a/ldl/README.txt b/ldl/README.txt new file mode 100644 index 000000000..7be8dd1f0 --- /dev/null +++ b/ldl/README.txt @@ -0,0 +1,136 @@ +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 new file mode 100644 index 000000000..7b5e79e54 --- /dev/null +++ b/ldl/UFconfig.h @@ -0,0 +1,118 @@ +/* ========================================================================== */ +/* === 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 new file mode 100644 index 000000000..bda4c7b1d --- /dev/null +++ b/ldl/UFconfig.mk @@ -0,0 +1,354 @@ +#=============================================================================== +# 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 new file mode 100644 index 000000000..a9b35c846 --- /dev/null +++ b/ldl/ldl.cpp @@ -0,0 +1,507 @@ +/* ========================================================================== */ +/* === 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 new file mode 100644 index 000000000..2ef71e07f --- /dev/null +++ b/ldl/ldl.h @@ -0,0 +1,104 @@ +/* ========================================================================== */ +/* === 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 new file mode 100644 index 000000000..87a355565 --- /dev/null +++ b/ldl/ldlsimple.cpp @@ -0,0 +1,69 @@ +/* ========================================================================== */ +/* === 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/matlab/example/alex01_simple1.m b/matlab/example/alex01_simple1.m index fa36f1e88..cfa648369 100644 --- a/matlab/example/alex01_simple1.m +++ b/matlab/example/alex01_simple1.m @@ -30,32 +30,33 @@ for i=1:maxIt % evaluate derivatives in x dfx = [2*x1; 2*(x2-2)]; dcx = [8*x1; 2*x2]; + dL = dfx - lam * dcx; % update the hessian (BFGS) if (i>1) Bis = Bi*s; - y = dfx - prev_dfx; + y = dL - prev_dL; Bi = Bi + (y*y')/(y'*s) - (Bis*Bis')/(s'*Bis); end - prev_dfx = dfx; + prev_dL = dL; % evaluate hessians in x ddfx = diag([2, 2]); ddcx = diag([8, 2]); % construct and solve CQP subproblem - dL = dfx - lam * dcx; + Bgn0 = dfx * dfx'; Bgn1 = dfx * dfx' - lam * dcx * dcx'; % GN approx 1 Bgn2 = dL * dL'; % GN approx 2 Ba = ddfx - lam * ddcx; % analytic hessians - B = Bi; + B = ddfx; g = dfx; h = -cx; [delta lambda] = solveCQP(B, -dcx, -dcx', g, h); % update - s = 0.2*delta; + s = 0.5*delta; x = x + s lam = lambda