diff options
Diffstat (limited to 'gsl-1.9/multiroots')
| -rw-r--r-- | gsl-1.9/multiroots/ChangeLog | 114 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/Makefile.am | 19 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/Makefile.in | 552 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/broyden.c | 456 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/convergence.c | 90 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/dnewton.c | 186 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/dogleg.c | 413 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/enorm.c | 33 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/fdfsolver.c | 176 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/fdjac.c | 96 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/fsolver.c | 163 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/gnewton.c | 226 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/gsl_multiroots.h | 177 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/hybrid.c | 661 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/hybridj.c | 608 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/newton.c | 151 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/test.c | 260 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/test_funcs.c | 744 | ||||
| -rw-r--r-- | gsl-1.9/multiroots/test_funcs.h | 76 |
19 files changed, 5201 insertions, 0 deletions
diff --git a/gsl-1.9/multiroots/ChangeLog b/gsl-1.9/multiroots/ChangeLog new file mode 100644 index 0000000..c84ad1a --- /dev/null +++ b/gsl-1.9/multiroots/ChangeLog @@ -0,0 +1,114 @@ +2007-01-26 Brian Gough <bjg@network-theory.co.uk> + + * fsolver.c (gsl_multiroot_fsolver_set): made vector argument x + const + + * fdfsolver.c (gsl_multiroot_fdfsolver_set): made vector argument + x const + +Tue Nov 12 22:26:40 2002 Brian Gough <bjg@network-theory.co.uk> + + * newton.c (newton_alloc): return error code, not null + + * hybridj.c (hybridj_alloc): return error code, not null + + * hybrid.c (hybrid_alloc): return error code, not null + + * gnewton.c (gnewton_alloc): return error code, not null + + * dnewton.c (dnewton_alloc): return error code, not null + + * broyden.c (broyden_alloc): return error code, not null + +Wed May 1 21:40:55 2002 Brian Gough <bjg@network-theory.co.uk> + + * fdfsolver.c (gsl_multiroot_fdfsolver_dx): new function to return + dx + (gsl_multiroot_fdfsolver_f): new function to return f + + * fsolver.c (gsl_multiroot_fsolver_dx): new function to return dx + (gsl_multiroot_fsolver_f): new function to return f + +Sat Jan 26 17:11:34 2002 Brian Gough <bjg@network-theory.co.uk> + + * hybrid.c (set): fix broken 'if' statement + +Thu Jan 10 19:26:55 2002 Brian Gough <bjg@network-theory.co.uk> + + * hybrid.c dnewton.c: return status values for user-function + (Theis Peter Hansen) + +Tue Jun 19 23:40:18 2001 Brian Gough <bjg@network-theory.co.uk> + + * hybrid.c: removed workspace for linalg calls, no longer needed + + * hybridj.c: removed workspace for linalg calls, no longer needed + +Wed Jun 6 13:31:18 2001 Brian Gough <bjg@network-theory.co.uk> + + * hybridj.c: updated to use new QR calling convention (now passes + workspace) + + * hybrid.c: updated to use new QR calling convention (now passes + workspace) + +Mon Apr 23 12:55:39 2001 Brian Gough <bjg@network-theory.co.uk> + + * Makefile.am (test_LDADD): added cblas lib + +Mon Apr 16 20:18:08 2001 Brian Gough <bjg@network-theory.co.uk> + + * dnewton.c (dnewton_iterate): removed unnecessary status variable + +Sun Feb 18 11:26:45 2001 Brian Gough <bjg@network-theory.co.uk> + + * fdfsolver.c fsolver.c: changed so that the solver _alloc + function no longer calls _set, the user must do that separately. + +Thu Nov 30 21:48:39 2000 Brian Gough <bjg@network-theory.co.uk> + + * newton.c (newton_iterate): return GSL_EBADFUNC if error in + function evaluation + + * hybridj.c (iterate): return GSL_EBADFUNC if error in + function evaluation + + * hybrid.c (iterate): return GSL_EBADFUNC if error in function + evaluation + + * gnewton.c (gnewton_iterate): return GSL_EBADFUNC if error in + function evaluation + + * fdjac.c (gsl_multiroot_fdjacobian): return GSL_EBADFUNC if error + in function evaluation + + * dnewton.c (dnewton_iterate): return GSL_EBADFUNC if error in + function evaluation + + * broyden.c (broyden_iterate): return GSL_EBADFUNC if error in + function evaluation + +Sun Aug 27 13:43:13 2000 Brian Gough <bjg@network-theory.co.uk> + + * hybridj.c hybrid.c dogleg.c: begin comments with a capital + letter to improve readability + +Sat Aug 26 16:12:19 2000 Brian Gough <bjg@network-theory.co.uk> + + * hybridj.c hybrid.c: renamed rdiag to tau, since it plays that + role here and is not the diagonal of R (see qr.c documentation for + more details) + +Wed Feb 23 15:36:39 2000 Brian Gough <bjg@network-theory.co.uk> + + * changed gsl_vector_copy to gsl_vector_cpy + +Fri Feb 18 18:45:02 2000 Brian Gough <bjg@network-theory.co.uk> + + * fixed various .c files to use permutation + +Wed Feb 16 21:13:24 2000 Brian Gough <bjg@network-theory.co.uk> + + * fixed Makefiles that include gsl_linalg.h to add + -I$(srcdir)/../permutation to their include path + diff --git a/gsl-1.9/multiroots/Makefile.am b/gsl-1.9/multiroots/Makefile.am new file mode 100644 index 0000000..fbc772d --- /dev/null +++ b/gsl-1.9/multiroots/Makefile.am @@ -0,0 +1,19 @@ +# -*-makefile-*- + +noinst_LTLIBRARIES = libgslmultiroots.la + +pkginclude_HEADERS = gsl_multiroots.h + +noinst_HEADERS = enorm.c dogleg.c + +INCLUDES= -I$(top_builddir) + +libgslmultiroots_la_SOURCES = fdjac.c fsolver.c fdfsolver.c convergence.c newton.c gnewton.c dnewton.c broyden.c hybrid.c hybridj.c + +check_PROGRAMS = test + +TESTS = $(check_PROGRAMS) + +test_SOURCES = test.c test_funcs.c test_funcs.h +test_LDADD = libgslmultiroots.la ../linalg/libgsllinalg.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../permutation/libgslpermutation.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + diff --git a/gsl-1.9/multiroots/Makefile.in b/gsl-1.9/multiroots/Makefile.in new file mode 100644 index 0000000..22e900e --- /dev/null +++ b/gsl-1.9/multiroots/Makefile.in @@ -0,0 +1,552 @@ +# Makefile.in generated by automake 1.9.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005 Free Software Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + +# -*-makefile-*- + + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +top_builddir = .. +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +INSTALL = @INSTALL@ +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +check_PROGRAMS = test$(EXEEXT) +subdir = multiroots +DIST_COMMON = $(noinst_HEADERS) $(pkginclude_HEADERS) \ + $(srcdir)/Makefile.am $(srcdir)/Makefile.in ChangeLog +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(SHELL) $(top_srcdir)/mkinstalldirs +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +LTLIBRARIES = $(noinst_LTLIBRARIES) +libgslmultiroots_la_LIBADD = +am_libgslmultiroots_la_OBJECTS = fdjac.lo fsolver.lo fdfsolver.lo \ + convergence.lo newton.lo gnewton.lo dnewton.lo broyden.lo \ + hybrid.lo hybridj.lo +libgslmultiroots_la_OBJECTS = $(am_libgslmultiroots_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) test_funcs.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslmultiroots.la ../linalg/libgsllinalg.la \ + ../blas/libgslblas.la ../cblas/libgslcblas.la \ + ../permutation/libgslpermutation.la ../matrix/libgslmatrix.la \ + ../vector/libgslvector.la ../block/libgslblock.la \ + ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la \ + ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la \ + ../utils/libutils.la +DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir) +depcomp = +am__depfiles_maybe = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +SOURCES = $(libgslmultiroots_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslmultiroots_la_SOURCES) $(test_SOURCES) +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = `echo $$p | sed -e 's|^.*/||'`; +am__installdirs = "$(DESTDIR)$(pkgincludedir)" +pkgincludeHEADERS_INSTALL = $(INSTALL_HEADER) +HEADERS = $(noinst_HEADERS) $(pkginclude_HEADERS) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +ECHO = @ECHO@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GSL_LT_CBLAS_VERSION = @GSL_LT_CBLAS_VERSION@ +GSL_LT_VERSION = @GSL_LT_VERSION@ +HAVE_AIX_IEEE_INTERFACE = @HAVE_AIX_IEEE_INTERFACE@ +HAVE_DARWIN86_IEEE_INTERFACE = @HAVE_DARWIN86_IEEE_INTERFACE@ +HAVE_DARWIN_IEEE_INTERFACE = @HAVE_DARWIN_IEEE_INTERFACE@ +HAVE_EXTENDED_PRECISION_REGISTERS = @HAVE_EXTENDED_PRECISION_REGISTERS@ +HAVE_FREEBSD_IEEE_INTERFACE = @HAVE_FREEBSD_IEEE_INTERFACE@ +HAVE_GNUM68K_IEEE_INTERFACE = @HAVE_GNUM68K_IEEE_INTERFACE@ +HAVE_GNUPPC_IEEE_INTERFACE = @HAVE_GNUPPC_IEEE_INTERFACE@ +HAVE_GNUSPARC_IEEE_INTERFACE = @HAVE_GNUSPARC_IEEE_INTERFACE@ +HAVE_GNUX86_IEEE_INTERFACE = @HAVE_GNUX86_IEEE_INTERFACE@ +HAVE_HPUX11_IEEE_INTERFACE = @HAVE_HPUX11_IEEE_INTERFACE@ +HAVE_HPUX_IEEE_INTERFACE = @HAVE_HPUX_IEEE_INTERFACE@ +HAVE_IEEE_COMPARISONS = @HAVE_IEEE_COMPARISONS@ +HAVE_IEEE_DENORMALS = @HAVE_IEEE_DENORMALS@ +HAVE_INLINE = @HAVE_INLINE@ +HAVE_IRIX_IEEE_INTERFACE = @HAVE_IRIX_IEEE_INTERFACE@ +HAVE_NETBSD_IEEE_INTERFACE = @HAVE_NETBSD_IEEE_INTERFACE@ +HAVE_OPENBSD_IEEE_INTERFACE = @HAVE_OPENBSD_IEEE_INTERFACE@ +HAVE_OS2EMX_IEEE_INTERFACE = @HAVE_OS2EMX_IEEE_INTERFACE@ +HAVE_PRINTF_LONGDOUBLE = @HAVE_PRINTF_LONGDOUBLE@ +HAVE_SOLARIS_IEEE_INTERFACE = @HAVE_SOLARIS_IEEE_INTERFACE@ +HAVE_SUNOS4_IEEE_INTERFACE = @HAVE_SUNOS4_IEEE_INTERFACE@ +HAVE_TRU64_IEEE_INTERFACE = @HAVE_TRU64_IEEE_INTERFACE@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAINT = @MAINT@ +MAINTAINER_MODE_FALSE = @MAINTAINER_MODE_FALSE@ +MAINTAINER_MODE_TRUE = @MAINTAINER_MODE_TRUE@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +RELEASED = @RELEASED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_STRIP = @ac_ct_STRIP@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +datadir = @datadir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +noinst_LTLIBRARIES = libgslmultiroots.la +pkginclude_HEADERS = gsl_multiroots.h +noinst_HEADERS = enorm.c dogleg.c +INCLUDES = -I$(top_builddir) +libgslmultiroots_la_SOURCES = fdjac.c fsolver.c fdfsolver.c convergence.c newton.c gnewton.c dnewton.c broyden.c hybrid.c hybridj.c +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c test_funcs.c test_funcs.h +test_LDADD = libgslmultiroots.la ../linalg/libgsllinalg.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../permutation/libgslpermutation.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \ + && exit 0; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps multiroots/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps multiroots/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +clean-noinstLTLIBRARIES: + -test -z "$(noinst_LTLIBRARIES)" || rm -f $(noinst_LTLIBRARIES) + @list='$(noinst_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +libgslmultiroots.la: $(libgslmultiroots_la_OBJECTS) $(libgslmultiroots_la_DEPENDENCIES) + $(LINK) $(libgslmultiroots_la_LDFLAGS) $(libgslmultiroots_la_OBJECTS) $(libgslmultiroots_la_LIBADD) $(LIBS) + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; for p in $$list; do \ + f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f $$p $$f"; \ + rm -f $$p $$f ; \ + done +test$(EXEEXT): $(test_OBJECTS) $(test_DEPENDENCIES) + @rm -f test$(EXEEXT) + $(LINK) $(test_LDFLAGS) $(test_OBJECTS) $(test_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +distclean-libtool: + -rm -f libtool +uninstall-info-am: +install-pkgincludeHEADERS: $(pkginclude_HEADERS) + @$(NORMAL_INSTALL) + test -z "$(pkgincludedir)" || $(mkdir_p) "$(DESTDIR)$(pkgincludedir)" + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \ + f=$(am__strip_dir) \ + echo " $(pkgincludeHEADERS_INSTALL) '$$d$$p' '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + $(pkgincludeHEADERS_INSTALL) "$$d$$p" "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +uninstall-pkgincludeHEADERS: + @$(NORMAL_UNINSTALL) + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + f=$(am__strip_dir) \ + echo " rm -f '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + rm -f "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +check-TESTS: $(TESTS) + @failed=0; all=0; xfail=0; xpass=0; skip=0; \ + srcdir=$(srcdir); export srcdir; \ + list='$(TESTS)'; \ + if test -n "$$list"; then \ + for tst in $$list; do \ + if test -f ./$$tst; then dir=./; \ + elif test -f $$tst; then dir=; \ + else dir="$(srcdir)/"; fi; \ + if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xpass=`expr $$xpass + 1`; \ + failed=`expr $$failed + 1`; \ + echo "XPASS: $$tst"; \ + ;; \ + *) \ + echo "PASS: $$tst"; \ + ;; \ + esac; \ + elif test $$? -ne 77; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xfail=`expr $$xfail + 1`; \ + echo "XFAIL: $$tst"; \ + ;; \ + *) \ + failed=`expr $$failed + 1`; \ + echo "FAIL: $$tst"; \ + ;; \ + esac; \ + else \ + skip=`expr $$skip + 1`; \ + echo "SKIP: $$tst"; \ + fi; \ + done; \ + if test "$$failed" -eq 0; then \ + if test "$$xfail" -eq 0; then \ + banner="All $$all tests passed"; \ + else \ + banner="All $$all tests behaved as expected ($$xfail expected failures)"; \ + fi; \ + else \ + if test "$$xpass" -eq 0; then \ + banner="$$failed of $$all tests failed"; \ + else \ + banner="$$failed of $$all tests did not behave as expected ($$xpass unexpected passes)"; \ + fi; \ + fi; \ + dashes="$$banner"; \ + skipped=""; \ + if test "$$skip" -ne 0; then \ + skipped="($$skip tests were not run)"; \ + test `echo "$$skipped" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$skipped"; \ + fi; \ + report=""; \ + if test "$$failed" -ne 0 && test -n "$(PACKAGE_BUGREPORT)"; then \ + report="Please report to $(PACKAGE_BUGREPORT)"; \ + test `echo "$$report" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$report"; \ + fi; \ + dashes=`echo "$$dashes" | sed s/./=/g`; \ + echo "$$dashes"; \ + echo "$$banner"; \ + test -z "$$skipped" || echo "$$skipped"; \ + test -z "$$report" || echo "$$report"; \ + echo "$$dashes"; \ + test "$$failed" -eq 0; \ + else :; fi + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-TESTS +check: check-am +all-am: Makefile $(LTLIBRARIES) $(HEADERS) +installdirs: + for dir in "$(DESTDIR)$(pkgincludedir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-libtool distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: install-pkgincludeHEADERS + +install-exec-am: + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-info-am uninstall-pkgincludeHEADERS + +.PHONY: CTAGS GTAGS all all-am check check-TESTS check-am clean \ + clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES ctags distclean distclean-compile \ + distclean-generic distclean-libtool distclean-tags distdir dvi \ + dvi-am html html-am info info-am install install-am \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man \ + install-pkgincludeHEADERS install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-info-am \ + uninstall-pkgincludeHEADERS + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gsl-1.9/multiroots/broyden.c b/gsl-1.9/multiroots/broyden.c new file mode 100644 index 0000000..dbec753 --- /dev/null +++ b/gsl-1.9/multiroots/broyden.c @@ -0,0 +1,456 @@ +/* multiroots/broyden.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +#include "enorm.c" + +/* Broyden's method. It is not an efficient or modern algorithm but + gives an example of a rank-1 update. + + C.G. Broyden, "A Class of Methods for Solving Nonlinear + Simultaneous Equations", Mathematics of Computation, vol 19 (1965), + p 577-593 + + */ + +typedef struct + { + gsl_matrix *H; + gsl_matrix *lu; + gsl_permutation *permutation; + gsl_vector *v; + gsl_vector *w; + gsl_vector *y; + gsl_vector *p; + gsl_vector *fnew; + gsl_vector *x_trial; + double phi; + } +broyden_state_t; + +static int broyden_alloc (void *vstate, size_t n); +static int broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static int broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static void broyden_free (void *vstate); + + +static int +broyden_alloc (void *vstate, size_t n) +{ + broyden_state_t *state = (broyden_state_t *) vstate; + gsl_vector *v, *w, *y, *fnew, *x_trial, *p; + gsl_permutation *perm; + gsl_matrix *m, *H; + + m = gsl_matrix_calloc (n, n); + + if (m == 0) + { + GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM); + } + + state->lu = m; + + perm = gsl_permutation_calloc (n); + + if (perm == 0) + { + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM); + } + + state->permutation = perm; + + H = gsl_matrix_calloc (n, n); + + if (H == 0) + { + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM); + } + + state->H = H; + + v = gsl_vector_calloc (n); + + if (v == 0) + { + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM); + } + + state->v = v; + + w = gsl_vector_calloc (n); + + if (w == 0) + { + gsl_vector_free (v); + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM); + } + + state->w = w; + + y = gsl_vector_calloc (n); + + if (y == 0) + { + gsl_vector_free (w); + gsl_vector_free (v); + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM); + } + + state->y = y; + + fnew = gsl_vector_calloc (n); + + if (fnew == 0) + { + gsl_vector_free (y); + gsl_vector_free (w); + gsl_vector_free (v); + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for fnew", GSL_ENOMEM); + } + + state->fnew = fnew; + + x_trial = gsl_vector_calloc (n); + + if (x_trial == 0) + { + gsl_vector_free (fnew); + gsl_vector_free (y); + gsl_vector_free (w); + gsl_vector_free (v); + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM); + } + + state->x_trial = x_trial; + + p = gsl_vector_calloc (n); + + if (p == 0) + { + gsl_vector_free (x_trial); + gsl_vector_free (fnew); + gsl_vector_free (y); + gsl_vector_free (w); + gsl_vector_free (v); + gsl_matrix_free (H); + gsl_permutation_free (perm); + gsl_matrix_free (m); + + GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM); + } + + state->p = p; + + return GSL_SUCCESS; +} + +static int +broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx) +{ + broyden_state_t *state = (broyden_state_t *) vstate; + size_t i, j, n = function->n; + int signum = 0; + + GSL_MULTIROOT_FN_EVAL (function, x, f); + + gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->lu); + gsl_linalg_LU_decomp (state->lu, state->permutation, &signum); + gsl_linalg_LU_invert (state->lu, state->permutation, state->H); + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + gsl_matrix_set(state->H,i,j,-gsl_matrix_get(state->H,i,j)); + + for (i = 0; i < n; i++) + { + gsl_vector_set (dx, i, 0.0); + } + + state->phi = enorm (f); + + return GSL_SUCCESS; +} + +static int +broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx) +{ + broyden_state_t *state = (broyden_state_t *) vstate; + + double phi0, phi1, t, lambda; + + gsl_matrix *H = state->H; + gsl_vector *p = state->p; + gsl_vector *y = state->y; + gsl_vector *v = state->v; + gsl_vector *w = state->w; + gsl_vector *fnew = state->fnew; + gsl_vector *x_trial = state->x_trial; + gsl_matrix *lu = state->lu; + gsl_permutation *perm = state->permutation; + + size_t i, j, iter; + + size_t n = function->n; + + /* p = H f */ + + for (i = 0; i < n; i++) + { + double sum = 0; + + for (j = 0; j < n; j++) + { + sum += gsl_matrix_get (H, i, j) * gsl_vector_get (f, j); + } + gsl_vector_set (p, i, sum); + } + + t = 1; + iter = 0; + + phi0 = state->phi; + +new_step: + + for (i = 0; i < n; i++) + { + double pi = gsl_vector_get (p, i); + double xi = gsl_vector_get (x, i); + gsl_vector_set (x_trial, i, xi + t * pi); + } + + { + int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + phi1 = enorm (fnew); + + iter++ ; + + if (phi1 > phi0 && iter < 10 && t > 0.1) + { + /* full step goes uphill, take a reduced step instead */ + + double theta = phi1 / phi0; + t *= (sqrt (1.0 + 6.0 * theta) - 1.0) / (3.0 * theta); + goto new_step; + } + + if (phi1 > phi0) + { + /* need to recompute Jacobian */ + int signum = 0; + + gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, lu); + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + gsl_matrix_set(lu,i,j,-gsl_matrix_get(lu,i,j)); + + gsl_linalg_LU_decomp (lu, perm, &signum); + gsl_linalg_LU_invert (lu, perm, H); + + gsl_linalg_LU_solve (lu, perm, f, p); + + t = 1; + + for (i = 0; i < n; i++) + { + double pi = gsl_vector_get (p, i); + double xi = gsl_vector_get (x, i); + gsl_vector_set (x_trial, i, xi + t * pi); + } + + { + int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + phi1 = enorm (fnew); + } + + /* y = f' - f */ + + for (i = 0; i < n; i++) + { + double yi = gsl_vector_get (fnew, i) - gsl_vector_get (f, i); + gsl_vector_set (y, i, yi); + } + + /* v = H y */ + + for (i = 0; i < n; i++) + { + double sum = 0; + + for (j = 0; j < n; j++) + { + sum += gsl_matrix_get (H, i, j) * gsl_vector_get (y, j); + } + + gsl_vector_set (v, i, sum); + } + + /* lambda = p . v */ + + lambda = 0; + + for (i = 0; i < n; i++) + { + lambda += gsl_vector_get (p, i) * gsl_vector_get (v, i); + } + + if (lambda == 0) + { + GSL_ERROR ("approximation to Jacobian has collapsed", GSL_EZERODIV) ; + } + + /* v' = v + t * p */ + + for (i = 0; i < n; i++) + { + double vi = gsl_vector_get (v, i) + t * gsl_vector_get (p, i); + gsl_vector_set (v, i, vi); + } + + /* w^T = p^T H */ + + for (i = 0; i < n; i++) + { + double sum = 0; + + for (j = 0; j < n; j++) + { + sum += gsl_matrix_get (H, j, i) * gsl_vector_get (p, j); + } + + gsl_vector_set (w, i, sum); + } + + /* Hij -> Hij - (vi wj / lambda) */ + + for (i = 0; i < n; i++) + { + double vi = gsl_vector_get (v, i); + + for (j = 0; j < n; j++) + { + double wj = gsl_vector_get (w, j); + double Hij = gsl_matrix_get (H, i, j) - vi * wj / lambda; + gsl_matrix_set (H, i, j, Hij); + } + } + + /* copy fnew into f */ + + gsl_vector_memcpy (f, fnew); + + /* copy x_trial into x */ + + gsl_vector_memcpy (x, x_trial); + + for (i = 0; i < n; i++) + { + double pi = gsl_vector_get (p, i); + gsl_vector_set (dx, i, t * pi); + } + + state->phi = phi1; + + return GSL_SUCCESS; +} + + +static void +broyden_free (void *vstate) +{ + broyden_state_t *state = (broyden_state_t *) vstate; + + gsl_matrix_free (state->H); + + gsl_matrix_free (state->lu); + gsl_permutation_free (state->permutation); + + gsl_vector_free (state->v); + gsl_vector_free (state->w); + gsl_vector_free (state->y); + gsl_vector_free (state->p); + + gsl_vector_free (state->fnew); + gsl_vector_free (state->x_trial); + +} + + +static const gsl_multiroot_fsolver_type broyden_type = +{"broyden", /* name */ + sizeof (broyden_state_t), + &broyden_alloc, + &broyden_set, + &broyden_iterate, + &broyden_free}; + +const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_broyden = &broyden_type; diff --git a/gsl-1.9/multiroots/convergence.c b/gsl-1.9/multiroots/convergence.c new file mode 100644 index 0000000..1934662 --- /dev/null +++ b/gsl-1.9/multiroots/convergence.c @@ -0,0 +1,90 @@ +/* multiroots/convergence.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> + +int +gsl_multiroot_test_delta (const gsl_vector * dx, const gsl_vector * x, + double epsabs, double epsrel) +{ + size_t i; + int ok = 1; + const size_t n = x->size ; + + if (epsrel < 0.0) + { + GSL_ERROR ("relative tolerance is negative", GSL_EBADTOL); + } + + for (i = 0 ; i < n ; i++) + { + double xi = gsl_vector_get(x,i); + double dxi = gsl_vector_get(dx,i); + double tolerance = epsabs + epsrel * fabs(xi) ; + + if (fabs(dxi) < tolerance) + { + ok = 1; + } + else + { + ok = 0; + break; + } + } + + if (ok) + return GSL_SUCCESS ; + + return GSL_CONTINUE; +} + +int +gsl_multiroot_test_residual (const gsl_vector * f, double epsabs) +{ + size_t i; + + double residual = 0; + + const size_t n = f->size; + + if (epsabs < 0.0) + { + GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL); + } + + for (i = 0 ; i < n ; i++) + { + double fi = gsl_vector_get(f, i); + + residual += fabs(fi); + } + + + if (residual < epsabs) + { + return GSL_SUCCESS; + } + + return GSL_CONTINUE ; +} + diff --git a/gsl-1.9/multiroots/dnewton.c b/gsl-1.9/multiroots/dnewton.c new file mode 100644 index 0000000..5a43a74 --- /dev/null +++ b/gsl-1.9/multiroots/dnewton.c @@ -0,0 +1,186 @@ +/* multiroots/dnewton.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +/* Newton method using a finite difference approximation to the jacobian. + The derivatives are estimated using a step size of + + h_i = sqrt(DBL_EPSILON) * x_i + + */ + +typedef struct + { + gsl_matrix * J; + gsl_matrix * lu; + gsl_permutation * permutation; + } +dnewton_state_t; + +static int dnewton_alloc (void * vstate, size_t n); +static int dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static int dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static void dnewton_free (void * vstate); + +static int +dnewton_alloc (void * vstate, size_t n) +{ + dnewton_state_t * state = (dnewton_state_t *) vstate; + gsl_permutation * p; + gsl_matrix * m, * J; + + m = gsl_matrix_calloc (n,n); + + if (m == 0) + { + GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM); + } + + state->lu = m ; + + p = gsl_permutation_calloc (n); + + if (p == 0) + { + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM); + } + + state->permutation = p ; + + J = gsl_matrix_calloc (n,n); + + if (J == 0) + { + gsl_permutation_free(p); + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM); + } + + state->J = J; + + return GSL_SUCCESS; +} + +static int +dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx) +{ + dnewton_state_t * state = (dnewton_state_t *) vstate; + size_t i, n = function->n ; + int status; + + status = GSL_MULTIROOT_FN_EVAL (function, x, f); + if (status) + return status; + + status = gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, + state->J); + if (status) + return status; + + for (i = 0; i < n; i++) + { + gsl_vector_set (dx, i, 0.0); + } + + return GSL_SUCCESS; +} + +static int +dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx) +{ + dnewton_state_t * state = (dnewton_state_t *) vstate; + + int signum ; + + size_t i; + + size_t n = function->n ; + + gsl_matrix_memcpy (state->lu, state->J); + + { + int status = gsl_linalg_LU_decomp (state->lu, state->permutation, &signum); + if (status) + return status; + } + + { + int status = gsl_linalg_LU_solve (state->lu, state->permutation, f, dx); + if (status) + return status; + } + + for (i = 0; i < n; i++) + { + double e = gsl_vector_get (dx, i); + double y = gsl_vector_get (x, i); + gsl_vector_set (dx, i, -e); + gsl_vector_set (x, i, y - e); + } + + { + int status = GSL_MULTIROOT_FN_EVAL (function, x, f); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->J); + + return GSL_SUCCESS; +} + + +static void +dnewton_free (void * vstate) +{ + dnewton_state_t * state = (dnewton_state_t *) vstate; + + gsl_matrix_free(state->J); + gsl_matrix_free(state->lu); + gsl_permutation_free(state->permutation); +} + + +static const gsl_multiroot_fsolver_type dnewton_type = +{"dnewton", /* name */ + sizeof (dnewton_state_t), + &dnewton_alloc, + &dnewton_set, + &dnewton_iterate, + &dnewton_free}; + +const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_dnewton = &dnewton_type; diff --git a/gsl-1.9/multiroots/dogleg.c b/gsl-1.9/multiroots/dogleg.c new file mode 100644 index 0000000..4a120ef --- /dev/null +++ b/gsl-1.9/multiroots/dogleg.c @@ -0,0 +1,413 @@ +/* multiroots/dogleg.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include "enorm.c" + +static void compute_diag (const gsl_matrix * J, gsl_vector * diag); +static void update_diag (const gsl_matrix * J, gsl_vector * diag); +static double compute_delta (gsl_vector * diag, gsl_vector * x); +static void compute_df (const gsl_vector * f_trial, const gsl_vector * f, gsl_vector * df); +static void compute_wv (const gsl_vector * qtdf, const gsl_vector *rdx, const gsl_vector *dx, const gsl_vector *diag, double pnorm, gsl_vector * w, gsl_vector * v); + +static double scaled_enorm (const gsl_vector * d, const gsl_vector * f); + +static double scaled_enorm (const gsl_vector * d, const gsl_vector * f) { + double e2 = 0 ; + size_t i, n = f->size ; + for (i = 0; i < n ; i++) { + double fi= gsl_vector_get(f, i); + double di= gsl_vector_get(d, i); + double u = di * fi; + e2 += u * u ; + } + return sqrt(e2); +} + +static double enorm_sum (const gsl_vector * a, const gsl_vector * b); + +static double enorm_sum (const gsl_vector * a, const gsl_vector * b) { + double e2 = 0 ; + size_t i, n = a->size ; + for (i = 0; i < n ; i++) { + double ai= gsl_vector_get(a, i); + double bi= gsl_vector_get(b, i); + double u = ai + bi; + e2 += u * u ; + } + return sqrt(e2); +} + +static void +compute_wv (const gsl_vector * qtdf, const gsl_vector *rdx, const gsl_vector *dx, const gsl_vector *diag, double pnorm, gsl_vector * w, gsl_vector * v) +{ + size_t i, n = qtdf->size; + + for (i = 0; i < n; i++) + { + double qtdfi = gsl_vector_get (qtdf, i); + double rdxi = gsl_vector_get (rdx, i); + double dxi = gsl_vector_get (dx, i); + double diagi = gsl_vector_get (diag, i); + + gsl_vector_set (w, i, (qtdfi - rdxi) / pnorm); + gsl_vector_set (v, i, diagi * diagi * dxi / pnorm); + } +} + + +static void +compute_df (const gsl_vector * f_trial, const gsl_vector * f, gsl_vector * df) +{ + size_t i, n = f->size; + + for (i = 0; i < n; i++) + { + double dfi = gsl_vector_get (f_trial, i) - gsl_vector_get (f, i); + gsl_vector_set (df, i, dfi); + } +} + +static void +compute_diag (const gsl_matrix * J, gsl_vector * diag) +{ + size_t i, j, n = diag->size; + + for (j = 0; j < n; j++) + { + double sum = 0; + for (i = 0; i < n; i++) + { + double Jij = gsl_matrix_get (J, i, j); + sum += Jij * Jij; + } + if (sum == 0) + sum = 1.0; + + gsl_vector_set (diag, j, sqrt (sum)); + } +} + +static void +update_diag (const gsl_matrix * J, gsl_vector * diag) +{ + size_t i, j, n = diag->size; + + for (j = 0; j < n; j++) + { + double cnorm, diagj, sum = 0; + for (i = 0; i < n; i++) + { + double Jij = gsl_matrix_get (J, i, j); + sum += Jij * Jij; + } + if (sum == 0) + sum = 1.0; + + cnorm = sqrt (sum); + diagj = gsl_vector_get (diag, j); + + if (cnorm > diagj) + gsl_vector_set (diag, j, cnorm); + } +} + +static double +compute_delta (gsl_vector * diag, gsl_vector * x) +{ + double Dx = scaled_enorm (diag, x); + double factor = 100; + + return (Dx > 0) ? factor * Dx : factor; +} + +static double +compute_actual_reduction (double fnorm, double fnorm1) +{ + double actred; + + if (fnorm1 < fnorm) + { + double u = fnorm1 / fnorm; + actred = 1 - u * u; + } + else + { + actred = -1; + } + + return actred; +} + +static double +compute_predicted_reduction (double fnorm, double fnorm1) +{ + double prered; + + if (fnorm1 < fnorm) + { + double u = fnorm1 / fnorm; + prered = 1 - u * u; + } + else + { + prered = 0; + } + + return prered; +} + +static void +compute_qtf (const gsl_matrix * q, const gsl_vector * f, gsl_vector * qtf) +{ + size_t i, j, N = f->size ; + + for (j = 0; j < N; j++) + { + double sum = 0; + for (i = 0; i < N; i++) + sum += gsl_matrix_get (q, i, j) * gsl_vector_get (f, i); + + gsl_vector_set (qtf, j, sum); + } +} + +static void +compute_rdx (const gsl_matrix * r, const gsl_vector * dx, gsl_vector * rdx) +{ + size_t i, j, N = dx->size ; + + for (i = 0; i < N; i++) + { + double sum = 0; + + for (j = i; j < N; j++) + { + sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, j); + } + + gsl_vector_set (rdx, i, sum); + } +} + + +static void +compute_trial_step (gsl_vector *x, gsl_vector * dx, gsl_vector * x_trial) +{ + size_t i, N = x->size; + + for (i = 0; i < N; i++) + { + double pi = gsl_vector_get (dx, i); + double xi = gsl_vector_get (x, i); + gsl_vector_set (x_trial, i, xi + pi); + } +} + +static int +newton_direction (const gsl_matrix * r, const gsl_vector * qtf, gsl_vector * p) +{ + const size_t N = r->size2; + size_t i; + int status; + + status = gsl_linalg_R_solve (r, qtf, p); + +#ifdef DEBUG + printf("rsolve status = %d\n", status); +#endif + + for (i = 0; i < N; i++) + { + double pi = gsl_vector_get (p, i); + gsl_vector_set (p, i, -pi); + } + + return status; +} + +static void +gradient_direction (const gsl_matrix * r, const gsl_vector * qtf, + const gsl_vector * diag, gsl_vector * g) +{ + const size_t M = r->size1; + const size_t N = r->size2; + + size_t i, j; + + for (j = 0; j < M; j++) + { + double sum = 0; + double dj; + + for (i = 0; i < N; i++) + { + sum += gsl_matrix_get (r, i, j) * gsl_vector_get (qtf, i); + } + + dj = gsl_vector_get (diag, j); + gsl_vector_set (g, j, -sum / dj); + } +} + +static void +minimum_step (double gnorm, const gsl_vector * diag, gsl_vector * g) +{ + const size_t N = g->size; + size_t i; + + for (i = 0; i < N; i++) + { + double gi = gsl_vector_get (g, i); + double di = gsl_vector_get (diag, i); + gsl_vector_set (g, i, (gi / gnorm) / di); + } +} + +static void +compute_Rg (const gsl_matrix * r, const gsl_vector * gradient, gsl_vector * Rg) +{ + const size_t N = r->size2; + size_t i, j; + + for (i = 0; i < N; i++) + { + double sum = 0; + + for (j = i; j < N; j++) + { + double gj = gsl_vector_get (gradient, j); + double rij = gsl_matrix_get (r, i, j); + sum += rij * gj; + } + + gsl_vector_set (Rg, i, sum); + } +} + +static void +scaled_addition (double alpha, gsl_vector * newton, double beta, gsl_vector * gradient, gsl_vector * p) +{ + const size_t N = p->size; + size_t i; + + for (i = 0; i < N; i++) + { + double ni = gsl_vector_get (newton, i); + double gi = gsl_vector_get (gradient, i); + gsl_vector_set (p, i, alpha * ni + beta * gi); + } +} + +static int +dogleg (const gsl_matrix * r, const gsl_vector * qtf, + const gsl_vector * diag, double delta, + gsl_vector * newton, gsl_vector * gradient, gsl_vector * p) +{ + double qnorm, gnorm, sgnorm, bnorm, temp; + + newton_direction (r, qtf, newton); + +#ifdef DEBUG + printf("newton = "); gsl_vector_fprintf(stdout, newton, "%g"); printf("\n"); +#endif + + qnorm = scaled_enorm (diag, newton); + + if (qnorm <= delta) + { + gsl_vector_memcpy (p, newton); +#ifdef DEBUG + printf("took newton (qnorm = %g <= delta = %g)\n", qnorm, delta); +#endif + return GSL_SUCCESS; + } + + gradient_direction (r, qtf, diag, gradient); + +#ifdef DEBUG + printf("grad = "); gsl_vector_fprintf(stdout, gradient, "%g"); printf("\n"); +#endif + + gnorm = enorm (gradient); + + if (gnorm == 0) + { + double alpha = delta / qnorm; + double beta = 0; + scaled_addition (alpha, newton, beta, gradient, p); +#ifdef DEBUG + printf("took scaled newton because gnorm = 0\n"); +#endif + return GSL_SUCCESS; + } + + minimum_step (gnorm, diag, gradient); + + compute_Rg (r, gradient, p); /* Use p as temporary space to compute Rg */ + +#ifdef DEBUG + printf("mingrad = "); gsl_vector_fprintf(stdout, gradient, "%g"); printf("\n"); + printf("Rg = "); gsl_vector_fprintf(stdout, p, "%g"); printf("\n"); +#endif + + temp = enorm (p); + sgnorm = (gnorm / temp) / temp; + + if (sgnorm > delta) + { + double alpha = 0; + double beta = delta; + scaled_addition (alpha, newton, beta, gradient, p); +#ifdef DEBUG + printf("took gradient\n"); +#endif + return GSL_SUCCESS; + } + + bnorm = enorm (qtf); + + { + double bg = bnorm / gnorm; + double bq = bnorm / qnorm; + double dq = delta / qnorm; + double dq2 = dq * dq; + double sd = sgnorm / delta; + double sd2 = sd * sd; + + double t1 = bg * bq * sd; + double u = t1 - dq; + double t2 = t1 - dq * sd2 + sqrt (u * u + (1-dq2) * (1 - sd2)); + + double alpha = dq * (1 - sd2) / t2; + double beta = (1 - alpha) * sgnorm; + +#ifdef DEBUG + printf("bnorm = %g\n", bnorm); + printf("gnorm = %g\n", gnorm); + printf("qnorm = %g\n", qnorm); + printf("delta = %g\n", delta); + printf("alpha = %g beta = %g\n", alpha, beta); + printf("took scaled combination of newton and gradient\n"); +#endif + + scaled_addition (alpha, newton, beta, gradient, p); + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multiroots/enorm.c b/gsl-1.9/multiroots/enorm.c new file mode 100644 index 0000000..4ad40a1 --- /dev/null +++ b/gsl-1.9/multiroots/enorm.c @@ -0,0 +1,33 @@ +/* multiroots/enorm.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +static double enorm (const gsl_vector * f); + +static double enorm (const gsl_vector * f) { + double e2 = 0 ; + size_t i, n = f->size ; + for (i = 0; i < n ; i++) { + double fi= gsl_vector_get(f, i); + e2 += fi * fi ; + } + return sqrt(e2); +} + + + diff --git a/gsl-1.9/multiroots/fdfsolver.c b/gsl-1.9/multiroots/fdfsolver.c new file mode 100644 index 0000000..0c34bba --- /dev/null +++ b/gsl-1.9/multiroots/fdfsolver.c @@ -0,0 +1,176 @@ +/* multiroots/fdfsolver.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> + +gsl_multiroot_fdfsolver * +gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, + size_t n) +{ + int status; + + gsl_multiroot_fdfsolver * s; + + s = (gsl_multiroot_fdfsolver *) malloc (sizeof (gsl_multiroot_fdfsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for multiroot solver struct", + GSL_ENOMEM, 0); + } + + s->x = gsl_vector_calloc (n); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->f = gsl_vector_calloc (n); + + if (s->f == 0) + { + gsl_vector_free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0); + } + + s->J = gsl_matrix_calloc (n,n); + + if (s->J == 0) + { + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0); + } + + s->dx = gsl_vector_calloc (n); + + if (s->dx == 0) + { + gsl_matrix_free (s->J); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0); + } + + s->state = malloc (T->size); + + if (s->state == 0) + { + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for multiroot solver state", + GSL_ENOMEM, 0); + } + + s->type = T ; + + status = (s->type->alloc)(s->state, n); + + if (status != GSL_SUCCESS) + { + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to set solver", status, 0); + } + + s->fdf = NULL; + + return s; +} + +int +gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s, + gsl_multiroot_function_fdf * f, + const gsl_vector * x) +{ + if (s->x->size != f->n) + { + GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN); + } + + if (x->size != f->n) + { + GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN); + } + + s->fdf = f; + gsl_vector_memcpy(s->x,x); + + return (s->type->set) (s->state, s->fdf, s->x, s->f, s->J, s->dx); +} + +int +gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * s) +{ + return (s->type->iterate) (s->state, s->fdf, s->x, s->f, s->J, s->dx); +} + +void +gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * s) +{ + (s->type->free) (s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); +} + +const char * +gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * s) +{ + return s->type->name; +} + +gsl_vector * +gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * s) +{ + return s->x; +} + +gsl_vector * +gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * s) +{ + return s->dx; +} + +gsl_vector * +gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * s) +{ + return s->f; +} diff --git a/gsl-1.9/multiroots/fdjac.c b/gsl-1.9/multiroots/fdjac.c new file mode 100644 index 0000000..6628589 --- /dev/null +++ b/gsl-1.9/multiroots/fdjac.c @@ -0,0 +1,96 @@ +/* multiroots/fdjac.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <gsl/gsl_multiroots.h> + +int +gsl_multiroot_fdjacobian (gsl_multiroot_function * F, + const gsl_vector * x, const gsl_vector * f, + double epsrel, gsl_matrix * jacobian) +{ + const size_t n = x->size; + const size_t m = f->size; + const size_t n1 = jacobian->size1; + const size_t n2 = jacobian->size2; + + if (m != n1 || n != n2) + { + GSL_ERROR ("function and jacobian are not conformant", GSL_EBADLEN); + } + + { + size_t i,j; + gsl_vector *x1, *f1; + + x1 = gsl_vector_alloc (n); + + if (x1 == 0) + { + GSL_ERROR ("failed to allocate space for x1 workspace", GSL_ENOMEM); + } + + f1 = gsl_vector_alloc (m); + + if (f1 == 0) + { + gsl_vector_free (x1); + + GSL_ERROR ("failed to allocate space for f1 workspace", GSL_ENOMEM); + } + + gsl_vector_memcpy (x1, x); /* copy x into x1 */ + + for (j = 0; j < n; j++) + { + double xj = gsl_vector_get (x, j); + double dx = epsrel * fabs (xj); + + if (dx == 0) + { + dx = epsrel; + } + + gsl_vector_set (x1, j, xj + dx); + + { + int status = GSL_MULTIROOT_FN_EVAL (F, x1, f1); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + gsl_vector_set (x1, j, xj); + + for (i = 0; i < m; i++) + { + double g1 = gsl_vector_get (f1, i); + double g0 = gsl_vector_get (f, i); + gsl_matrix_set (jacobian, i, j, (g1 - g0) / dx); + } + } + + gsl_vector_free (x1); + gsl_vector_free (f1); + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multiroots/fsolver.c b/gsl-1.9/multiroots/fsolver.c new file mode 100644 index 0000000..bb334a4 --- /dev/null +++ b/gsl-1.9/multiroots/fsolver.c @@ -0,0 +1,163 @@ +/* multiroots/fsolver.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> + +gsl_multiroot_fsolver * +gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, + size_t n) +{ + int status; + + gsl_multiroot_fsolver * s; + + s = (gsl_multiroot_fsolver *) malloc (sizeof (gsl_multiroot_fsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for multiroot solver struct", + GSL_ENOMEM, 0); + } + + s->x = gsl_vector_calloc (n); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->f = gsl_vector_calloc (n); + + if (s->f == 0) + { + gsl_vector_free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0); + } + + s->dx = gsl_vector_calloc (n); + + if (s->dx == 0) + { + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0); + } + + s->state = malloc (T->size); + + if (s->state == 0) + { + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for multiroot solver state", + GSL_ENOMEM, 0); + } + + s->type = T ; + + status = (s->type->alloc)(s->state, n); + + if (status != GSL_SUCCESS) + { + (s->type->free)(s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to set solver", status, 0); + } + + s->function = NULL; + + return s; +} + +int +gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s, + gsl_multiroot_function * f, + const gsl_vector * x) +{ + if (s->x->size != f->n) + { + GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN); + } + + if (x->size != f->n) + { + GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN); + } + + s->function = f; + gsl_vector_memcpy(s->x,x); + + return (s->type->set) (s->state, s->function, s->x, s->f, s->dx); +} + +int +gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * s) +{ + return (s->type->iterate) (s->state, s->function, s->x, s->f, s->dx); +} + +void +gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s) +{ + (s->type->free) (s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); +} + +const char * +gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * s) +{ + return s->type->name; +} + +gsl_vector * +gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * s) +{ + return s->x; +} + +gsl_vector * +gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * s) +{ + return s->dx; +} + +gsl_vector * +gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * s) +{ + return s->f; +} diff --git a/gsl-1.9/multiroots/gnewton.c b/gsl-1.9/multiroots/gnewton.c new file mode 100644 index 0000000..13ade76 --- /dev/null +++ b/gsl-1.9/multiroots/gnewton.c @@ -0,0 +1,226 @@ +/* multiroots/gnewton.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +#include "enorm.c" + +/* Simple globally convergent Newton method (rejects uphill steps) */ + +typedef struct + { + double phi; + gsl_vector * x_trial; + gsl_vector * d; + gsl_matrix * lu; + gsl_permutation * permutation; + } +gnewton_state_t; + +static int gnewton_alloc (void * vstate, size_t n); +static int gnewton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static void gnewton_free (void * vstate); + +static int +gnewton_alloc (void * vstate, size_t n) +{ + gnewton_state_t * state = (gnewton_state_t *) vstate; + gsl_vector * d, * x_trial ; + gsl_permutation * p; + gsl_matrix * m; + + m = gsl_matrix_calloc (n,n); + + if (m == 0) + { + GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM); + } + + state->lu = m ; + + p = gsl_permutation_calloc (n); + + if (p == 0) + { + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM); + } + + state->permutation = p ; + + d = gsl_vector_calloc (n); + + if (d == 0) + { + gsl_permutation_free(p); + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM); + } + + state->d = d; + + x_trial = gsl_vector_calloc (n); + + if (x_trial == 0) + { + gsl_vector_free(d); + gsl_permutation_free(p); + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM); + } + + state->x_trial = x_trial; + + return GSL_SUCCESS; +} + + +static int +gnewton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + gnewton_state_t * state = (gnewton_state_t *) vstate; + size_t i, n = FDF->n ; + + GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J); + + for (i = 0; i < n; i++) + { + gsl_vector_set (dx, i, 0.0); + } + + state->phi = enorm(f); + + return GSL_SUCCESS; +} + +static int +gnewton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + gnewton_state_t * state = (gnewton_state_t *) vstate; + + int signum ; + double t, phi0, phi1; + + size_t i; + + size_t n = fdf->n ; + + gsl_matrix_memcpy (state->lu, J); + + gsl_linalg_LU_decomp (state->lu, state->permutation, &signum); + + gsl_linalg_LU_solve (state->lu, state->permutation, f, state->d); + + t = 1; + + phi0 = state->phi; + +new_step: + + for (i = 0; i < n; i++) + { + double di = gsl_vector_get (state->d, i); + double xi = gsl_vector_get (x, i); + gsl_vector_set (state->x_trial, i, xi - t*di); + } + + { + int status = GSL_MULTIROOT_FN_EVAL_F (fdf, state->x_trial, f); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + phi1 = enorm (f); + + if (phi1 > phi0 && t > GSL_DBL_EPSILON) + { + /* full step goes uphill, take a reduced step instead */ + + double theta = phi1 / phi0; + double u = (sqrt(1.0 + 6.0 * theta) - 1.0) / (3.0 * theta); + + t *= u ; + + goto new_step; + } + + /* copy x_trial into x */ + + gsl_vector_memcpy (x, state->x_trial); + + for (i = 0; i < n; i++) + { + double di = gsl_vector_get (state->d, i); + gsl_vector_set (dx, i, -t*di); + } + + { + int status = GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + state->phi = phi1; + + return GSL_SUCCESS; +} + + +static void +gnewton_free (void * vstate) +{ + gnewton_state_t * state = (gnewton_state_t *) vstate; + + gsl_vector_free(state->d); + gsl_vector_free(state->x_trial); + gsl_matrix_free(state->lu); + gsl_permutation_free(state->permutation); +} + + +static const gsl_multiroot_fdfsolver_type gnewton_type = +{"gnewton", /* name */ + sizeof (gnewton_state_t), + &gnewton_alloc, + &gnewton_set, + &gnewton_iterate, + &gnewton_free}; + +const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_gnewton = &gnewton_type; diff --git a/gsl-1.9/multiroots/gsl_multiroots.h b/gsl-1.9/multiroots/gsl_multiroots.h new file mode 100644 index 0000000..6ea5681 --- /dev/null +++ b/gsl-1.9/multiroots/gsl_multiroots.h @@ -0,0 +1,177 @@ +/* multiroots/gsl_multiroots.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_MULTIROOTS_H__ +#define __GSL_MULTIROOTS_H__ + +#include <stdlib.h> +#include <gsl/gsl_types.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +/* Definition of vector-valued functions with parameters based on gsl_vector */ + +struct gsl_multiroot_function_struct +{ + int (* f) (const gsl_vector * x, void * params, gsl_vector * f); + size_t n; + void * params; +}; + +typedef struct gsl_multiroot_function_struct gsl_multiroot_function ; + +#define GSL_MULTIROOT_FN_EVAL(F,x,y) (*((F)->f))(x,(F)->params,(y)) + +int gsl_multiroot_fdjacobian (gsl_multiroot_function * F, + const gsl_vector * x, const gsl_vector * f, + double epsrel, gsl_matrix * jacobian); + + +typedef struct + { + const char *name; + size_t size; + int (*alloc) (void *state, size_t n); + int (*set) (void *state, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); + int (*iterate) (void *state, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); + void (*free) (void *state); + } +gsl_multiroot_fsolver_type; + +typedef struct + { + const gsl_multiroot_fsolver_type * type; + gsl_multiroot_function * function ; + gsl_vector * x ; + gsl_vector * f ; + gsl_vector * dx ; + void *state; + } +gsl_multiroot_fsolver; + +gsl_multiroot_fsolver * +gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, + size_t n); + +void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * s); + +int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * s, + gsl_multiroot_function * f, + const gsl_vector * x); + +int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * s); + +const char * gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * s); +gsl_vector * gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * s); +gsl_vector * gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * s); +gsl_vector * gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * s); + +/* Definition of vector-valued functions and gradient with parameters + based on gsl_vector */ + +struct gsl_multiroot_function_fdf_struct +{ + int (* f) (const gsl_vector * x, void * params, gsl_vector * f); + int (* df) (const gsl_vector * x, void * params, gsl_matrix * df); + int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix *df); + size_t n; + void * params; +}; + +typedef struct gsl_multiroot_function_fdf_struct gsl_multiroot_function_fdf ; + +#define GSL_MULTIROOT_FN_EVAL_F(F,x,y) ((*((F)->f))(x,(F)->params,(y))) +#define GSL_MULTIROOT_FN_EVAL_DF(F,x,dy) ((*((F)->df))(x,(F)->params,(dy))) +#define GSL_MULTIROOT_FN_EVAL_F_DF(F,x,y,dy) ((*((F)->fdf))(x,(F)->params,(y),(dy))) + +typedef struct + { + const char *name; + size_t size; + int (*alloc) (void *state, size_t n); + int (*set) (void *state, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); + int (*iterate) (void *state, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); + void (*free) (void *state); + } +gsl_multiroot_fdfsolver_type; + +typedef struct + { + const gsl_multiroot_fdfsolver_type * type; + gsl_multiroot_function_fdf * fdf ; + gsl_vector * x; + gsl_vector * f; + gsl_matrix * J; + gsl_vector * dx; + void *state; + } +gsl_multiroot_fdfsolver; + +gsl_multiroot_fdfsolver * +gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, + size_t n); + +int +gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * s, + gsl_multiroot_function_fdf * fdf, + const gsl_vector * x); + +int +gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * s); + +void +gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * s); + +const char * gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * s); +gsl_vector * gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * s); +gsl_vector * gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * s); +gsl_vector * gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * s); + +int gsl_multiroot_test_delta (const gsl_vector * dx, const gsl_vector * x, + double epsabs, double epsrel); + +int gsl_multiroot_test_residual (const gsl_vector * f, double epsabs); + +GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_dnewton; +GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_broyden; +GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_hybrid; +GSL_VAR const gsl_multiroot_fsolver_type * gsl_multiroot_fsolver_hybrids; + +GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_newton; +GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_gnewton; +GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_hybridj; +GSL_VAR const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_hybridsj; + + +__END_DECLS + +#endif /* __GSL_MULTIROOTS_H__ */ diff --git a/gsl-1.9/multiroots/hybrid.c b/gsl-1.9/multiroots/hybrid.c new file mode 100644 index 0000000..ddb65e5 --- /dev/null +++ b/gsl-1.9/multiroots/hybrid.c @@ -0,0 +1,661 @@ +/* multiroots/hybrid.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +#include "dogleg.c" + +typedef struct +{ + size_t iter; + size_t ncfail; + size_t ncsuc; + size_t nslow1; + size_t nslow2; + double fnorm; + double delta; + gsl_matrix *J; + gsl_matrix *q; + gsl_matrix *r; + gsl_vector *tau; + gsl_vector *diag; + gsl_vector *qtf; + gsl_vector *newton; + gsl_vector *gradient; + gsl_vector *x_trial; + gsl_vector *f_trial; + gsl_vector *df; + gsl_vector *qtdf; + gsl_vector *rdx; + gsl_vector *w; + gsl_vector *v; +} +hybrid_state_t; + +static int hybrid_alloc (void *vstate, size_t n); +static int hybrid_set (void *vstate, gsl_multiroot_function * func, + gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static int hybrids_set (void *vstate, gsl_multiroot_function * func, + gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static int set (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx, int scale); +static int hybrid_iterate (void *vstate, gsl_multiroot_function * func, + gsl_vector * x, gsl_vector * f, gsl_vector * dx); +static void hybrid_free (void *vstate); +static int iterate (void *vstate, gsl_multiroot_function * func, + gsl_vector * x, gsl_vector * f, gsl_vector * dx, + int scale); + +static int +hybrid_alloc (void *vstate, size_t n) +{ + hybrid_state_t *state = (hybrid_state_t *) vstate; + gsl_matrix *J, *q, *r; + gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial, + *df, *qtdf, *rdx, *w, *v; + + J = gsl_matrix_calloc (n, n); + + if (J == 0) + { + GSL_ERROR ("failed to allocate space for J", GSL_ENOMEM); + } + + state->J = J; + + q = gsl_matrix_calloc (n, n); + + if (q == 0) + { + gsl_matrix_free (J); + + GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM); + } + + state->q = q; + + r = gsl_matrix_calloc (n, n); + + if (r == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + + GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM); + } + + state->r = r; + + tau = gsl_vector_calloc (n); + + if (tau == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + + GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM); + } + + state->tau = tau; + + diag = gsl_vector_calloc (n); + + if (diag == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + + GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM); + } + + state->diag = diag; + + qtf = gsl_vector_calloc (n); + + if (qtf == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + + GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM); + } + + state->qtf = qtf; + + newton = gsl_vector_calloc (n); + + if (newton == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + + GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM); + } + + state->newton = newton; + + gradient = gsl_vector_calloc (n); + + if (gradient == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + + GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM); + } + + state->gradient = gradient; + + x_trial = gsl_vector_calloc (n); + + if (x_trial == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + + GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM); + } + + state->x_trial = x_trial; + + f_trial = gsl_vector_calloc (n); + + if (f_trial == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + + GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM); + } + + state->f_trial = f_trial; + + df = gsl_vector_calloc (n); + + if (df == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + + GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM); + } + + state->df = df; + + qtdf = gsl_vector_calloc (n); + + if (qtdf == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + + GSL_ERROR ("failed to allocate space for qtdf", GSL_ENOMEM); + } + + state->qtdf = qtdf; + + + rdx = gsl_vector_calloc (n); + + if (rdx == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + + GSL_ERROR ("failed to allocate space for rdx", GSL_ENOMEM); + } + + state->rdx = rdx; + + w = gsl_vector_calloc (n); + + if (w == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + gsl_vector_free (rdx); + + GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM); + } + + state->w = w; + + v = gsl_vector_calloc (n); + + if (v == 0) + { + gsl_matrix_free (J); + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + gsl_vector_free (rdx); + gsl_vector_free (w); + + GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM); + } + + state->v = v; + + return GSL_SUCCESS; +} + +static int +hybrid_set (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx) +{ + int status = set (vstate, func, x, f, dx, 0); + return status; +} + +static int +hybrids_set (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx) +{ + int status = set (vstate, func, x, f, dx, 1); + return status; +} + +static int +set (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx, int scale) +{ + hybrid_state_t *state = (hybrid_state_t *) vstate; + + gsl_matrix *J = state->J; + gsl_matrix *q = state->q; + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + + int status; + + status = GSL_MULTIROOT_FN_EVAL (func, x, f); + + if (status) + { + return status; + } + + status = gsl_multiroot_fdjacobian (func, x, f, GSL_SQRT_DBL_EPSILON, J); + + if (status) + { + return status; + } + + state->iter = 1; + state->fnorm = enorm (f); + state->ncfail = 0; + state->ncsuc = 0; + state->nslow1 = 0; + state->nslow2 = 0; + + gsl_vector_set_all (dx, 0.0); + + /* Store column norms in diag */ + + if (scale) + compute_diag (J, diag); + else + gsl_vector_set_all (diag, 1.0); + + /* Set delta to factor |D x| or to factor if |D x| is zero */ + + state->delta = compute_delta (diag, x); + + /* Factorize J into QR decomposition */ + + status = gsl_linalg_QR_decomp (J, tau); + + if (status) + { + return status; + } + + status = gsl_linalg_QR_unpack (J, tau, q, r); + + return status; +} + +static int +hybrid_iterate (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx) +{ + int status = iterate (vstate, func, x, f, dx, 0); + return status; +} + +static int +hybrids_iterate (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx) +{ + int status = iterate (vstate, func, x, f, dx, 1); + return status; +} + +static int +iterate (void *vstate, gsl_multiroot_function * func, gsl_vector * x, + gsl_vector * f, gsl_vector * dx, int scale) +{ + hybrid_state_t *state = (hybrid_state_t *) vstate; + + const double fnorm = state->fnorm; + + gsl_matrix *J = state->J; + gsl_matrix *q = state->q; + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + gsl_vector *qtf = state->qtf; + gsl_vector *x_trial = state->x_trial; + gsl_vector *f_trial = state->f_trial; + gsl_vector *df = state->df; + gsl_vector *qtdf = state->qtdf; + gsl_vector *rdx = state->rdx; + gsl_vector *w = state->w; + gsl_vector *v = state->v; + + double prered, actred; + double pnorm, fnorm1, fnorm1p; + double ratio; + double p1 = 0.1, p5 = 0.5, p001 = 0.001, p0001 = 0.0001; + + /* Compute qtf = Q^T f */ + + compute_qtf (q, f, qtf); + + /* Compute dogleg step */ + + dogleg (r, qtf, diag, state->delta, state->newton, state->gradient, dx); + + /* Take a trial step */ + + compute_trial_step (x, dx, state->x_trial); + + pnorm = scaled_enorm (diag, dx); + + if (state->iter == 1) + { + if (pnorm < state->delta) + { + state->delta = pnorm; + } + } + + /* Evaluate function at x + p */ + + { + int status = GSL_MULTIROOT_FN_EVAL (func, x_trial, f_trial); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + /* Set df = f_trial - f */ + + compute_df (f_trial, f, df); + + /* Compute the scaled actual reduction */ + + fnorm1 = enorm (f_trial); + + actred = compute_actual_reduction (fnorm, fnorm1); + + /* Compute rdx = R dx */ + + compute_rdx (r, dx, rdx); + + /* Compute the scaled predicted reduction phi1p = |Q^T f + R dx| */ + + fnorm1p = enorm_sum (qtf, rdx); + + prered = compute_predicted_reduction (fnorm, fnorm1p); + + /* Compute the ratio of the actual to predicted reduction */ + + if (prered > 0) + { + ratio = actred / prered; + } + else + { + ratio = 0; + } + + /* Update the step bound */ + + if (ratio < p1) + { + state->ncsuc = 0; + state->ncfail++; + state->delta *= p5; + } + else + { + state->ncfail = 0; + state->ncsuc++; + + if (ratio >= p5 || state->ncsuc > 1) + state->delta = GSL_MAX (state->delta, pnorm / p5); + if (fabs (ratio - 1) <= p1) + state->delta = pnorm / p5; + } + + /* Test for successful iteration */ + + if (ratio >= p0001) + { + gsl_vector_memcpy (x, x_trial); + gsl_vector_memcpy (f, f_trial); + state->fnorm = fnorm1; + state->iter++; + } + + /* Determine the progress of the iteration */ + + state->nslow1++; + if (actred >= p001) + state->nslow1 = 0; + + if (actred >= p1) + state->nslow2 = 0; + + if (state->ncfail == 2) + { + gsl_multiroot_fdjacobian (func, x, f, GSL_SQRT_DBL_EPSILON, J); + + state->nslow2++; + + if (state->iter == 1) + { + if (scale) + compute_diag (J, diag); + state->delta = compute_delta (diag, x); + } + else + { + if (scale) + update_diag (J, diag); + } + + /* Factorize J into QR decomposition */ + + gsl_linalg_QR_decomp (J, tau); + gsl_linalg_QR_unpack (J, tau, q, r); + + return GSL_SUCCESS; + } + + /* Compute qtdf = Q^T df, w = (Q^T df - R dx)/|dx|, v = D^2 dx/|dx| */ + + compute_qtf (q, df, qtdf); + + compute_wv (qtdf, rdx, dx, diag, pnorm, w, v); + + /* Rank-1 update of the jacobian Q'R' = Q(R + w v^T) */ + + gsl_linalg_QR_update (q, r, w, v); + + /* No progress as measured by jacobian evaluations */ + + if (state->nslow2 == 5) + { + return GSL_ENOPROGJ; + } + + /* No progress as measured by function evaluations */ + + if (state->nslow1 == 10) + { + return GSL_ENOPROG; + } + + return GSL_SUCCESS; +} + + +static void +hybrid_free (void *vstate) +{ + hybrid_state_t *state = (hybrid_state_t *) vstate; + + gsl_vector_free (state->v); + gsl_vector_free (state->w); + gsl_vector_free (state->rdx); + gsl_vector_free (state->qtdf); + gsl_vector_free (state->df); + gsl_vector_free (state->f_trial); + gsl_vector_free (state->x_trial); + gsl_vector_free (state->gradient); + gsl_vector_free (state->newton); + gsl_vector_free (state->qtf); + gsl_vector_free (state->diag); + gsl_vector_free (state->tau); + gsl_matrix_free (state->r); + gsl_matrix_free (state->q); + gsl_matrix_free (state->J); +} + +static const gsl_multiroot_fsolver_type hybrid_type = { + "hybrid", /* name */ + sizeof (hybrid_state_t), + &hybrid_alloc, + &hybrid_set, + &hybrid_iterate, + &hybrid_free +}; + +static const gsl_multiroot_fsolver_type hybrids_type = { + "hybrids", /* name */ + sizeof (hybrid_state_t), + &hybrid_alloc, + &hybrids_set, + &hybrids_iterate, + &hybrid_free +}; + +const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrid = &hybrid_type; +const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrids = + &hybrids_type; diff --git a/gsl-1.9/multiroots/hybridj.c b/gsl-1.9/multiroots/hybridj.c new file mode 100644 index 0000000..a44a164 --- /dev/null +++ b/gsl-1.9/multiroots/hybridj.c @@ -0,0 +1,608 @@ +/* multiroots/hybridj.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +#include "dogleg.c" + +typedef struct + { + size_t iter; + size_t ncfail; + size_t ncsuc; + size_t nslow1; + size_t nslow2; + double fnorm; + double delta; + gsl_matrix *q; + gsl_matrix *r; + gsl_vector *tau; + gsl_vector *diag; + gsl_vector *qtf; + gsl_vector *newton; + gsl_vector *gradient; + gsl_vector *x_trial; + gsl_vector *f_trial; + gsl_vector *df; + gsl_vector *qtdf; + gsl_vector *rdx; + gsl_vector *w; + gsl_vector *v; + } +hybridj_state_t; + +static int hybridj_alloc (void *vstate, size_t n); +static int hybridj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int hybridsj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale); +static int hybridj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static void hybridj_free (void *vstate); +static int iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale); + +static int +hybridj_alloc (void *vstate, size_t n) +{ + hybridj_state_t *state = (hybridj_state_t *) vstate; + gsl_matrix *q, *r; + gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial, + *df, *qtdf, *rdx, *w, *v; + + q = gsl_matrix_calloc (n, n); + + if (q == 0) + { + GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM); + } + + state->q = q; + + r = gsl_matrix_calloc (n, n); + + if (r == 0) + { + gsl_matrix_free (q); + + GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM); + } + + state->r = r; + + tau = gsl_vector_calloc (n); + + if (tau == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + + GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM); + } + + state->tau = tau; + + diag = gsl_vector_calloc (n); + + if (diag == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + + GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM); + } + + state->diag = diag; + + qtf = gsl_vector_calloc (n); + + if (qtf == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + + GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM); + } + + state->qtf = qtf; + + newton = gsl_vector_calloc (n); + + if (newton == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + + GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM); + } + + state->newton = newton; + + gradient = gsl_vector_calloc (n); + + if (gradient == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + + GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM); + } + + state->gradient = gradient; + + x_trial = gsl_vector_calloc (n); + + if (x_trial == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + + GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM); + } + + state->x_trial = x_trial; + + f_trial = gsl_vector_calloc (n); + + if (f_trial == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + + GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM); + } + + state->f_trial = f_trial; + + df = gsl_vector_calloc (n); + + if (df == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + + GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM); + } + + state->df = df; + + qtdf = gsl_vector_calloc (n); + + if (qtdf == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + + GSL_ERROR ("failed to allocate space for qtdf", GSL_ENOMEM); + } + + state->qtdf = qtdf; + + + rdx = gsl_vector_calloc (n); + + if (rdx == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + + GSL_ERROR ("failed to allocate space for rdx", GSL_ENOMEM); + } + + state->rdx = rdx; + + w = gsl_vector_calloc (n); + + if (w == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + gsl_vector_free (rdx); + + GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM); + } + + state->w = w; + + v = gsl_vector_calloc (n); + + if (v == 0) + { + gsl_matrix_free (q); + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (qtdf); + gsl_vector_free (rdx); + gsl_vector_free (w); + + GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM); + } + + state->v = v; + + return GSL_SUCCESS; +} + +static int +hybridj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = set (vstate, fdf, x, f, J, dx, 0); + return status ; +} + +static int +hybridsj_set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = set (vstate, fdf, x, f, J, dx, 1); + return status ; +} + +static int +set (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale) +{ + hybridj_state_t *state = (hybridj_state_t *) vstate; + + gsl_matrix *q = state->q; + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + + GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J); + + state->iter = 1; + state->fnorm = enorm (f); + state->ncfail = 0; + state->ncsuc = 0; + state->nslow1 = 0; + state->nslow2 = 0; + + gsl_vector_set_all (dx, 0.0); + + /* Store column norms in diag */ + + if (scale) + compute_diag (J, diag); + else + gsl_vector_set_all (diag, 1.0); + + /* Set delta to factor |D x| or to factor if |D x| is zero */ + + state->delta = compute_delta (diag, x); + + /* Factorize J into QR decomposition */ + + gsl_linalg_QR_decomp (J, tau); + gsl_linalg_QR_unpack (J, tau, q, r); + + return GSL_SUCCESS; +} + +static int +hybridj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = iterate (vstate, fdf, x, f, J, dx, 0); + return status; +} + +static int +hybridsj_iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = iterate (vstate, fdf, x, f, J, dx, 1); + return status; +} + +static int +iterate (void *vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale) +{ + hybridj_state_t *state = (hybridj_state_t *) vstate; + + const double fnorm = state->fnorm; + + gsl_matrix *q = state->q; + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + gsl_vector *qtf = state->qtf; + gsl_vector *x_trial = state->x_trial; + gsl_vector *f_trial = state->f_trial; + gsl_vector *df = state->df; + gsl_vector *qtdf = state->qtdf; + gsl_vector *rdx = state->rdx; + gsl_vector *w = state->w; + gsl_vector *v = state->v; + + double prered, actred; + double pnorm, fnorm1, fnorm1p; + double ratio; + double p1 = 0.1, p5 = 0.5, p001 = 0.001, p0001 = 0.0001; + + /* Compute qtf = Q^T f */ + + compute_qtf (q, f, qtf); + + /* Compute dogleg step */ + + dogleg (r, qtf, diag, state->delta, state->newton, state->gradient, dx); + + /* Take a trial step */ + + compute_trial_step (x, dx, state->x_trial); + + pnorm = scaled_enorm (diag, dx); + + if (state->iter == 1) + { + if (pnorm < state->delta) + { + state->delta = pnorm; + } + } + + /* Evaluate function at x + p */ + + { + int status = GSL_MULTIROOT_FN_EVAL_F (fdf, x_trial, f_trial); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + /* Set df = f_trial - f */ + + compute_df (f_trial, f, df); + + /* Compute the scaled actual reduction */ + + fnorm1 = enorm (f_trial); + + actred = compute_actual_reduction (fnorm, fnorm1); + + /* Compute rdx = R dx */ + + compute_rdx (r, dx, rdx); + + /* Compute the scaled predicted reduction phi1p = |Q^T f + R dx| */ + + fnorm1p = enorm_sum (qtf, rdx); + + prered = compute_predicted_reduction (fnorm, fnorm1p); + + /* Compute the ratio of the actual to predicted reduction */ + + if (prered > 0) + { + ratio = actred / prered; + } + else + { + ratio = 0; + } + + /* Update the step bound */ + + if (ratio < p1) + { + state->ncsuc = 0; + state->ncfail++; + state->delta *= p5; + } + else + { + state->ncfail = 0; + state->ncsuc++; + + if (ratio >= p5 || state->ncsuc > 1) + state->delta = GSL_MAX (state->delta, pnorm / p5); + if (fabs (ratio - 1) <= p1) + state->delta = pnorm / p5; + } + + /* Test for successful iteration */ + + if (ratio >= p0001) + { + gsl_vector_memcpy (x, x_trial); + gsl_vector_memcpy (f, f_trial); + state->fnorm = fnorm1; + state->iter++; + } + + /* Determine the progress of the iteration */ + + state->nslow1++; + if (actred >= p001) + state->nslow1 = 0; + + if (actred >= p1) + state->nslow2 = 0; + + if (state->ncfail == 2) + { + { + int status = GSL_MULTIROOT_FN_EVAL_DF (fdf, x, J); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + state->nslow2++; + + if (state->iter == 1) + { + if (scale) + compute_diag (J, diag); + state->delta = compute_delta (diag, x); + } + else + { + if (scale) + update_diag (J, diag); + } + + /* Factorize J into QR decomposition */ + + gsl_linalg_QR_decomp (J, tau); + gsl_linalg_QR_unpack (J, tau, q, r); + return GSL_SUCCESS; + } + + /* Compute qtdf = Q^T df, w = (Q^T df - R dx)/|dx|, v = D^2 dx/|dx| */ + + compute_qtf (q, df, qtdf); + + compute_wv (qtdf, rdx, dx, diag, pnorm, w, v); + + /* Rank-1 update of the jacobian Q'R' = Q(R + w v^T) */ + + gsl_linalg_QR_update (q, r, w, v); + + /* No progress as measured by jacobian evaluations */ + + if (state->nslow2 == 5) + { + return GSL_ENOPROGJ; + } + + /* No progress as measured by function evaluations */ + + if (state->nslow1 == 10) + { + return GSL_ENOPROG; + } + + return GSL_SUCCESS; +} + + +static void +hybridj_free (void *vstate) +{ + hybridj_state_t *state = (hybridj_state_t *) vstate; + + gsl_vector_free (state->v); + gsl_vector_free (state->w); + gsl_vector_free (state->rdx); + gsl_vector_free (state->qtdf); + gsl_vector_free (state->df); + gsl_vector_free (state->f_trial); + gsl_vector_free (state->x_trial); + gsl_vector_free (state->gradient); + gsl_vector_free (state->newton); + gsl_vector_free (state->qtf); + gsl_vector_free (state->diag); + gsl_vector_free (state->tau); + gsl_matrix_free (state->r); + gsl_matrix_free (state->q); +} + +static const gsl_multiroot_fdfsolver_type hybridj_type = +{ + "hybridj", /* name */ + sizeof (hybridj_state_t), + &hybridj_alloc, + &hybridj_set, + &hybridj_iterate, + &hybridj_free +}; + +static const gsl_multiroot_fdfsolver_type hybridsj_type = +{ + "hybridsj", /* name */ + sizeof (hybridj_state_t), + &hybridj_alloc, + &hybridsj_set, + &hybridsj_iterate, + &hybridj_free +}; + +const gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridj = &hybridj_type; +const gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridsj = &hybridsj_type; diff --git a/gsl-1.9/multiroots/newton.c b/gsl-1.9/multiroots/newton.c new file mode 100644 index 0000000..88509d0 --- /dev/null +++ b/gsl-1.9/multiroots/newton.c @@ -0,0 +1,151 @@ +/* multiroots/newton.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multiroots.h> +#include <gsl/gsl_linalg.h> + +typedef struct + { + gsl_matrix * lu; + gsl_permutation * permutation; + } +newton_state_t; + +static int newton_alloc (void * vstate, size_t n); +static int newton_set (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static void newton_free (void * vstate); + +static int +newton_alloc (void * vstate, size_t n) +{ + newton_state_t * state = (newton_state_t *) vstate; + gsl_permutation * p; + gsl_matrix * m; + + m = gsl_matrix_calloc (n,n); + + if (m == 0) + { + GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM); + } + + state->lu = m ; + + p = gsl_permutation_calloc (n); + + if (p == 0) + { + gsl_matrix_free(m); + + GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM); + } + + state->permutation = p ; + + return GSL_SUCCESS; +} + +static int +newton_set (void * vstate, gsl_multiroot_function_fdf * FDF, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + newton_state_t * state = (newton_state_t *) vstate; + + size_t i, n = FDF->n ; + + state = 0 ; /* avoid warnings about unused parameters */ + + GSL_MULTIROOT_FN_EVAL_F_DF (FDF, x, f, J); + + for (i = 0; i < n; i++) + { + gsl_vector_set (dx, i, 0.0); + } + + return GSL_SUCCESS; +} + +static int +newton_iterate (void * vstate, gsl_multiroot_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + newton_state_t * state = (newton_state_t *) vstate; + + int signum; + + size_t i; + + size_t n = fdf->n ; + + gsl_matrix_memcpy (state->lu, J); + + gsl_linalg_LU_decomp (state->lu, state->permutation, &signum); + + gsl_linalg_LU_solve (state->lu, state->permutation, f, dx); + + for (i = 0; i < n; i++) + { + double e = gsl_vector_get (dx, i); + double y = gsl_vector_get (x, i); + gsl_vector_set (dx, i, -e); + gsl_vector_set (x, i, y - e); + } + + { + int status = GSL_MULTIROOT_FN_EVAL_F_DF (fdf, x, f, J); + + if (status != GSL_SUCCESS) + { + return GSL_EBADFUNC; + } + } + + return GSL_SUCCESS; +} + + +static void +newton_free (void * vstate) +{ + newton_state_t * state = (newton_state_t *) vstate; + + gsl_matrix_free(state->lu); + + gsl_permutation_free(state->permutation); +} + + +static const gsl_multiroot_fdfsolver_type newton_type = +{"newton", /* name */ + sizeof (newton_state_t), + &newton_alloc, + &newton_set, + &newton_iterate, + &newton_free}; + +const gsl_multiroot_fdfsolver_type * gsl_multiroot_fdfsolver_newton = &newton_type; diff --git a/gsl-1.9/multiroots/test.c b/gsl-1.9/multiroots/test.c new file mode 100644 index 0000000..029f2b7 --- /dev/null +++ b/gsl-1.9/multiroots/test.c @@ -0,0 +1,260 @@ +/* multiroots/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <stdio.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_multiroots.h> + +#include <gsl/gsl_ieee_utils.h> + +#include "test_funcs.h" +int test_fdf (const char * desc, gsl_multiroot_function_fdf * function, initpt_function initpt, double factor, const gsl_multiroot_fdfsolver_type * T); +int test_f (const char * desc, gsl_multiroot_function_fdf * fdf, initpt_function initpt, double factor, const gsl_multiroot_fsolver_type * T); + + +int +main (void) +{ + const gsl_multiroot_fsolver_type * fsolvers[5] ; + const gsl_multiroot_fsolver_type ** T1 ; + + const gsl_multiroot_fdfsolver_type * fdfsolvers[5] ; + const gsl_multiroot_fdfsolver_type ** T2 ; + + double f; + + fsolvers[0] = gsl_multiroot_fsolver_dnewton; + fsolvers[1] = gsl_multiroot_fsolver_broyden; + fsolvers[2] = gsl_multiroot_fsolver_hybrid; + fsolvers[3] = gsl_multiroot_fsolver_hybrids; + fsolvers[4] = 0; + + fdfsolvers[0] = gsl_multiroot_fdfsolver_newton; + fdfsolvers[1] = gsl_multiroot_fdfsolver_gnewton; + fdfsolvers[2] = gsl_multiroot_fdfsolver_hybridj; + fdfsolvers[3] = gsl_multiroot_fdfsolver_hybridsj; + fdfsolvers[4] = 0; + + gsl_ieee_env_setup(); + + + f = 1.0 ; + + T1 = fsolvers ; + + while (*T1 != 0) + { + test_f ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T1); + test_f ("Roth", &roth, roth_initpt, f, *T1); + test_f ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T1); + test_f ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T1); + test_f ("Powell singular", &powellsing, powellsing_initpt, f, *T1); + test_f ("Wood", &wood, wood_initpt, f, *T1); + test_f ("Helical", &helical, helical_initpt, f, *T1); + test_f ("Discrete BVP", &dbv, dbv_initpt, f, *T1); + test_f ("Trig", &trig, trig_initpt, f, *T1); + T1++; + } + + T2 = fdfsolvers ; + + while (*T2 != 0) + { + test_fdf ("Rosenbrock", &rosenbrock, rosenbrock_initpt, f, *T2); + test_fdf ("Roth", &roth, roth_initpt, f, *T2); + test_fdf ("Powell badly scaled", &powellscal, powellscal_initpt, f, *T2); + test_fdf ("Brown badly scaled", &brownscal, brownscal_initpt, f, *T2); + test_fdf ("Powell singular", &powellsing, powellsing_initpt, f, *T2); + test_fdf ("Wood", &wood, wood_initpt, f, *T2); + test_fdf ("Helical", &helical, helical_initpt, f, *T2); + test_fdf ("Discrete BVP", &dbv, dbv_initpt, f, *T2); + test_fdf ("Trig", &trig, trig_initpt, f, *T2); + T2++; + } + + exit (gsl_test_summary ()); +} + +void scale (gsl_vector * x, double factor); + +void +scale (gsl_vector * x, double factor) +{ + size_t i, n = x->size; + + if (gsl_vector_isnull(x)) + { + for (i = 0; i < n; i++) + { + gsl_vector_set (x, i, factor); + } + } + else + { + for (i = 0; i < n; i++) + { + double xi = gsl_vector_get(x, i); + gsl_vector_set(x, i, factor * xi); + } + } +} + +int +test_fdf (const char * desc, gsl_multiroot_function_fdf * function, + initpt_function initpt, double factor, + const gsl_multiroot_fdfsolver_type * T) +{ + int status; + double residual = 0; + size_t i, n = function->n, iter = 0; + + gsl_vector *x = gsl_vector_alloc (n); + gsl_matrix *J = gsl_matrix_alloc (n, n); + + gsl_multiroot_fdfsolver *s; + + (*initpt) (x); + + if (factor != 1.0) scale(x, factor); + + s = gsl_multiroot_fdfsolver_alloc (T, n); + gsl_multiroot_fdfsolver_set (s, function, x); + + do + { + iter++; + status = gsl_multiroot_fdfsolver_iterate (s); + + if (status) + break ; + + status = gsl_multiroot_test_residual (s->f, 0.0000001); + } + while (status == GSL_CONTINUE && iter < 1000); + +#ifdef DEBUG + printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); + printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); +#endif + + +#ifdef TEST_JACOBIAN + { + double r,sum; size_t j; + + gsl_multiroot_function f1 ; + f1.f = function->f ; + f1.n = function->n ; + f1.params = function->params ; + + gsl_multiroot_fdjacobian (&f1, s->x, s->f, GSL_SQRT_DBL_EPSILON, J); + + /* compare J and s->J */ + + r=0;sum=0; + for (i = 0; i < n; i++) + for (j = 0; j< n ; j++) + { + double u = gsl_matrix_get(J,i,j); + double su = gsl_matrix_get(s->J, i, j); + r = fabs(u - su)/(1e-6 + 1e-6 * fabs(u)); sum+=r; + if (fabs(u - su) > 1e-6 + 1e-6 * fabs(u)) + printf("broken jacobian %g\n", r); + } + printf("avg r = %g\n", sum/(n*n)); + } +#endif + + for (i = 0; i < n ; i++) + { + residual += fabs(gsl_vector_get(s->f, i)); + } + + gsl_multiroot_fdfsolver_free (s); + gsl_matrix_free(J); + gsl_vector_free(x); + + gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual); + + return status; +} + + +int +test_f (const char * desc, gsl_multiroot_function_fdf * fdf, + initpt_function initpt, double factor, + const gsl_multiroot_fsolver_type * T) +{ + int status; + size_t i, n = fdf->n, iter = 0; + double residual = 0; + + gsl_vector *x; + + gsl_multiroot_fsolver *s; + gsl_multiroot_function function; + + function.f = fdf->f; + function.params = fdf->params; + function.n = n ; + + x = gsl_vector_alloc (n); + + (*initpt) (x); + + if (factor != 1.0) scale(x, factor); + + s = gsl_multiroot_fsolver_alloc (T, n); + gsl_multiroot_fsolver_set (s, &function, x); + +/* printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); */ +/* printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); */ + + do + { + iter++; + status = gsl_multiroot_fsolver_iterate (s); + + if (status) + break ; + + status = gsl_multiroot_test_residual (s->f, 0.0000001); + } + while (status == GSL_CONTINUE && iter < 1000); + +#ifdef DEBUG + printf("x "); gsl_vector_fprintf (stdout, s->x, "%g"); printf("\n"); + printf("f "); gsl_vector_fprintf (stdout, s->f, "%g"); printf("\n"); +#endif + + for (i = 0; i < n ; i++) + { + residual += fabs(gsl_vector_get(s->f, i)); + } + + gsl_multiroot_fsolver_free (s); + gsl_vector_free(x); + + gsl_test(status, "%s on %s (%g), %u iterations, residual = %.2g", T->name, desc, factor, iter, residual); + + return status; +} diff --git a/gsl-1.9/multiroots/test_funcs.c b/gsl-1.9/multiroots/test_funcs.c new file mode 100644 index 0000000..2edb2bb --- /dev/null +++ b/gsl-1.9/multiroots/test_funcs.c @@ -0,0 +1,744 @@ +/* multiroots/test_funcs.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_multiroots.h> + +#include "test_funcs.h" + +/* For information on testing see the following paper, + + J.J More, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained + Optimization Software", ACM Transactions on Mathematical Software, + Vol 7, No 1, (1981) p 17-41 + + */ + +/* Rosenbrock Function */ + +gsl_multiroot_function_fdf rosenbrock = +{&rosenbrock_f, + &rosenbrock_df, + &rosenbrock_fdf, + 2, 0}; + +void +rosenbrock_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, -1.2); + gsl_vector_set (x, 1, 1.0); +} + +int +rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double y0 = 1 - x0; + double y1 = 10 * (x1 - x0 * x0); + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + + double df00 = -1; + double df01 = 0; + double df10 = -20 * x0; + double df11 = 10; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +rosenbrock_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + rosenbrock_f (x, params, f); + rosenbrock_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Freudenstein and Roth function */ + +gsl_multiroot_function_fdf roth = +{&roth_f, + &roth_df, + &roth_fdf, + 2, 0}; + +void +roth_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, 4.5); /* changed from the value in the paper */ + gsl_vector_set (x, 1, 3.5); /* otherwise the problem is too hard */ +} + +int +roth_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double y0 = -13.0 + x0 + ((5.0 - x1)*x1 - 2.0)*x1; + double y1 = -29.0 + x0 + ((x1 + 1.0)*x1 - 14.0)*x1; + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +roth_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x1 = gsl_vector_get (x, 1); + + double df00 = 1; + double df01 = -3 * x1 * x1 + 10 * x1 - 2; + double df10 = 1; + double df11 = 3 * x1 * x1 + 2 * x1 - 14; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +roth_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + roth_f (x, params, f); + roth_df (x, params, df); + + return GSL_SUCCESS; +} + + + +/* Powell badly scaled function */ + +gsl_multiroot_function_fdf powellscal = +{&powellscal_f, + &powellscal_df, + &powellscal_fdf, + 2, 0}; + +void +powellscal_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, 0.0); + gsl_vector_set (x, 1, 1.0); +} + +int +powellscal_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double y0 = 10000.0 * x0 * x1 - 1.0; + double y1 = exp (-x0) + exp (-x1) - 1.0001; + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double df00 = 10000.0 * x1, df01 = 10000.0 * x0; + double df10 = -exp (-x0), df11 = -exp (-x1); + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +powellscal_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + powellscal_f (x, params, f); + powellscal_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Brown badly scaled function */ + +gsl_multiroot_function_fdf brownscal = +{&brownscal_f, + &brownscal_df, + &brownscal_fdf, + 2, 0}; + +void +brownscal_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, 1.0); + gsl_vector_set (x, 1, 1.0); +} + +int +brownscal_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double y0 = x0 - 1e6; + double y1 = x0 * x1 - 2; + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double df00 = 1.0, df01 = 0.0; + double df10 = x1, df11 = x0; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +brownscal_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + brownscal_f (x, params, f); + brownscal_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Powell Singular Function */ + +gsl_multiroot_function_fdf powellsing = +{&powellsing_f, + &powellsing_df, + &powellsing_fdf, + 4, 0}; + +void +powellsing_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, 3.0); + gsl_vector_set (x, 1, -1.0); + gsl_vector_set (x, 2, 0.0); + gsl_vector_set (x, 3, 1.0); +} + +int +powellsing_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + + double y0 = x0 + 10 * x1; + double y1 = sqrt (5.0) * (x2 - x3); + double y2 = pow (x1 - 2 * x2, 2.0); + double y3 = sqrt (10.0) * pow (x0 - x3, 2.0); + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + gsl_vector_set (f, 2, y2); + gsl_vector_set (f, 3, y3); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + + double df00 = 1, df01 = 10, df02 = 0, df03 = 0; + double df10 = 0, df11 = 0, df12 = sqrt (5.0), df13 = -df12; + double df20 = 0, df21 = 2 * (x1 - 2 * x2), df22 = -2 * df21, df23 = 0; + double df30 = 2 * sqrt (10.0) * (x0 - x3), df31 = 0, df32 = 0, df33 = -df30; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 0, 2, df02); + gsl_matrix_set (df, 0, 3, df03); + + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + gsl_matrix_set (df, 1, 2, df12); + gsl_matrix_set (df, 1, 3, df13); + + gsl_matrix_set (df, 2, 0, df20); + gsl_matrix_set (df, 2, 1, df21); + gsl_matrix_set (df, 2, 2, df22); + gsl_matrix_set (df, 2, 3, df23); + + gsl_matrix_set (df, 3, 0, df30); + gsl_matrix_set (df, 3, 1, df31); + gsl_matrix_set (df, 3, 2, df32); + gsl_matrix_set (df, 3, 3, df33); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +powellsing_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + powellsing_f (x, params, f); + powellsing_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Wood function */ + +gsl_multiroot_function_fdf wood = +{&wood_f, + &wood_df, + &wood_fdf, + 4, 0}; + +void +wood_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, -3.0); + gsl_vector_set (x, 1, -1.0); + gsl_vector_set (x, 2, -3.0); + gsl_vector_set (x, 3, -1.0); +} + +int +wood_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + + double t1 = x1 - x0 * x0; + double t2 = x3 - x2 * x2; + + double y0 = -200.0 * x0 * t1 - (1 - x0); + double y1 = 200.0 * t1 + 20.2 * (x1 - 1) + 19.8 * (x3 - 1); + double y2 = -180.0 * x2 * t2 - (1 - x2); + double y3 = 180.0 * t2 + 20.2 * (x3 - 1) + 19.8 * (x1 - 1); + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + gsl_vector_set (f, 2, y2); + gsl_vector_set (f, 3, y3); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +wood_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + + double t1 = x1 - 3 * x0 * x0; + double t2 = x3 - 3 * x2 * x2; + + double df00 = -200.0 * t1 + 1, df01 = -200.0 * x0, df02 = 0, df03 = 0; + double df10 = -400.0*x0, df11 = 200.0 + 20.2, df12 = 0, df13 = 19.8; + double df20 = 0, df21 = 0, df22 = -180.0 * t2 + 1, df23 = -180.0 * x2; + double df30 = 0, df31 = 19.8, df32 = -2 * 180 * x2, df33 = 180.0 + 20.2; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 0, 2, df02); + gsl_matrix_set (df, 0, 3, df03); + + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + gsl_matrix_set (df, 1, 2, df12); + gsl_matrix_set (df, 1, 3, df13); + + gsl_matrix_set (df, 2, 0, df20); + gsl_matrix_set (df, 2, 1, df21); + gsl_matrix_set (df, 2, 2, df22); + gsl_matrix_set (df, 2, 3, df23); + + gsl_matrix_set (df, 3, 0, df30); + gsl_matrix_set (df, 3, 1, df31); + gsl_matrix_set (df, 3, 2, df32); + gsl_matrix_set (df, 3, 3, df33); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +wood_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + wood_f (x, params, f); + wood_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Helical Valley Function */ + +gsl_multiroot_function_fdf helical = +{&helical_f, + &helical_df, + &helical_fdf, + 3, 0}; + +void +helical_initpt (gsl_vector * x) +{ + gsl_vector_set (x, 0, -1.0); + gsl_vector_set (x, 1, 0.0); + gsl_vector_set (x, 2, 0.0); +} + +int +helical_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + + double t1, t2; + double y0, y1, y2; + + if (x0 > 0) + { + t1 = atan(x1/x0) / (2.0 * M_PI); + } + else if (x0 < 0) + { + t1 = 0.5 + atan(x1/x0) / (2.0 * M_PI); + } + else + { + t1 = 0.25 * (x1 > 0 ? +1 : -1); + } + + t2 = sqrt(x0*x0 + x1*x1) ; + + y0 = 10 * (x2 - 10 * t1); + y1 = 10 * (t2 - 1); + y2 = x2 ; + + gsl_vector_set (f, 0, y0); + gsl_vector_set (f, 1, y1); + gsl_vector_set (f, 2, y2); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +helical_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + + double t = x0 * x0 + x1 * x1 ; + double t1 = 2 * M_PI * t ; + double t2 = sqrt(t) ; + + double df00 = 100*x1/t1, df01 = -100.0 * x0/t1, df02 = 10.0; + double df10 = 10*x0/t2, df11 = 10*x1/t2, df12 = 0; + double df20 = 0, df21 = 0, df22 = 1.0; + + gsl_matrix_set (df, 0, 0, df00); + gsl_matrix_set (df, 0, 1, df01); + gsl_matrix_set (df, 0, 2, df02); + + gsl_matrix_set (df, 1, 0, df10); + gsl_matrix_set (df, 1, 1, df11); + gsl_matrix_set (df, 1, 2, df12); + + gsl_matrix_set (df, 2, 0, df20); + gsl_matrix_set (df, 2, 1, df21); + gsl_matrix_set (df, 2, 2, df22); + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +helical_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + helical_f (x, params, f); + helical_df (x, params, df); + + return GSL_SUCCESS; +} + + +/* Discrete Boundary Value Function */ + +#define N 10 + +gsl_multiroot_function_fdf dbv = +{&dbv_f, + &dbv_df, + &dbv_fdf, + N, 0}; + +void +dbv_initpt (gsl_vector * x) +{ + size_t i; + double h = 1.0 / (N + 1.0); + + for (i = 0; i < N; i++) + { + double t = (i + 1) * h; + double z = t * (t - 1); + gsl_vector_set (x, i, z); + } +} + +int +dbv_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + size_t i; + + double h = 1.0 / (N + 1.0); + + for (i = 0; i < N; i++) + { + double z, ti = (i + 1) * h; + double xi = 0, xim1 = 0, xip1 = 0; + + xi = gsl_vector_get (x, i); + + if (i > 0) + xim1 = gsl_vector_get (x, i - 1); + + if (i < N - 1) + xip1 = gsl_vector_get (x, i + 1); + + z = 2 * xi - xim1 - xip1 + h * h * pow(xi + ti + 1, 3.0) / 2.0; + + gsl_vector_set (f, i, z); + + } + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +dbv_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + size_t i, j; + + double h = 1.0 / (N + 1.0); + + for (i = 0; i < N; i++) + for (j = 0; j < N; j++) + gsl_matrix_set (df, i, j, 0.0); + + for (i = 0; i < N; i++) + { + double dz_dxi, ti = (i + 1) * h; + + double xi = gsl_vector_get (x, i); + + dz_dxi = 2.0 + (3.0 / 2.0) * h * h * pow(xi + ti + 1, 2.0) ; + + gsl_matrix_set (df, i, i, dz_dxi); + + if (i > 0) + gsl_matrix_set (df, i, i-1, -1.0); + + if (i < N - 1) + gsl_matrix_set (df, i, i+1, -1.0); + + } + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +dbv_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + dbv_f (x, params, f); + dbv_df (x, params, df); + + return GSL_SUCCESS; +} + +/* Trigonometric Function */ + +gsl_multiroot_function_fdf trig = +{&trig_f, + &trig_df, + &trig_fdf, + N, 0}; + +void +trig_initpt (gsl_vector * x) +{ + size_t i; + + for (i = 0; i < N; i++) /* choose an initial point which converges */ + { + gsl_vector_set (x, i, 0.05); + } +} + +int +trig_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + size_t i; + double sum = 0; + + for (i = 0; i < N; i++) + { + sum += cos(gsl_vector_get(x,i)); + } + + for (i = 0; i < N; i++) + { + double xi = gsl_vector_get (x,i); + double z = N - sum + (i + 1) * (1 - cos(xi)) - sin(xi); + + gsl_vector_set (f, i, z); + } + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +trig_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + size_t i, j; + + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + double dz; + double xi = gsl_vector_get(x, i); + double xj = gsl_vector_get(x, j); + + if (j == i) + dz = sin(xi) + (i + 1) * sin(xi) - cos(xi); + else + dz = sin(xj); + + gsl_matrix_set(df, i, j, dz); + } + } + + params = 0; /* avoid warning about unused parameters */ + + return GSL_SUCCESS; +} + +int +trig_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + trig_f (x, params, f); + trig_df (x, params, df); + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multiroots/test_funcs.h b/gsl-1.9/multiroots/test_funcs.h new file mode 100644 index 0000000..ba7ee48 --- /dev/null +++ b/gsl-1.9/multiroots/test_funcs.h @@ -0,0 +1,76 @@ +/* multiroots/test_funcs.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + + +typedef void (*initpt_function) (gsl_vector * x); + +extern gsl_multiroot_function_fdf rosenbrock; +void rosenbrock_initpt (gsl_vector * x); +int rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f); +int rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df); +int rosenbrock_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf roth; +void roth_initpt (gsl_vector * x); +int roth_f (const gsl_vector * x, void *params, gsl_vector * f); +int roth_df (const gsl_vector * x, void *params, gsl_matrix * df); +int roth_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf brownscal; +void brownscal_initpt (gsl_vector * x); +int brownscal_f (const gsl_vector * x, void *params, gsl_vector * f); +int brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df); +int brownscal_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf powellscal; +void powellscal_initpt (gsl_vector * x); +int powellscal_f (const gsl_vector * x, void *params, gsl_vector * f); +int powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df); +int powellscal_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf powellsing; +void powellsing_initpt (gsl_vector * x); +int powellsing_f (const gsl_vector * x, void *params, gsl_vector * f); +int powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df); +int powellsing_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf wood; +void wood_initpt (gsl_vector * x); +int wood_f (const gsl_vector * x, void *params, gsl_vector * f); +int wood_df (const gsl_vector * x, void *params, gsl_matrix * df); +int wood_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf helical; +void helical_initpt (gsl_vector * x); +int helical_f (const gsl_vector * x, void *params, gsl_vector * f); +int helical_df (const gsl_vector * x, void *params, gsl_matrix * df); +int helical_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf dbv; +void dbv_initpt (gsl_vector * x); +int dbv_f (const gsl_vector * x, void *params, gsl_vector * f); +int dbv_df (const gsl_vector * x, void *params, gsl_matrix * df); +int dbv_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + +extern gsl_multiroot_function_fdf trig; +void trig_initpt (gsl_vector * x); +int trig_f (const gsl_vector * x, void *params, gsl_vector * f); +int trig_df (const gsl_vector * x, void *params, gsl_matrix * df); +int trig_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * df); + |