Mercurial > forge
changeset 11418:ed31a7de4a80 octave-forge
admin: closing down nonfree section for good
| author | carandraug |
|---|---|
| date | Mon, 28 Jan 2013 20:15:39 +0000 |
| parents | 1737676a0a47 |
| children | c8b7b87a88a4 |
| files | nonfree/CONTENTS nonfree/Makefile nonfree/spline-gsvspl/COPYING nonfree/spline-gsvspl/DESCRIPTION nonfree/spline-gsvspl/INDEX nonfree/spline-gsvspl/inst/csaps.m nonfree/spline-gsvspl/src/.svnignore nonfree/spline-gsvspl/src/Makeconf.in nonfree/spline-gsvspl/src/Makefile nonfree/spline-gsvspl/src/autogen.sh nonfree/spline-gsvspl/src/configure.base nonfree/spline-gsvspl/src/gcvspl.cc nonfree/spline-gsvspl/src/gcvsplf.f |
| diffstat | 13 files changed, 0 insertions(+), 2458 deletions(-) [+] |
line wrap: on
line diff
--- a/nonfree/CONTENTS Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -Packages which are not freely modifiable and redistributable -with modifications, or which depend on code which is not -free. This includes functions which only permit non-commercial -use and functions which must be kept together as a package. -Functions in all other directories must be freely redistributable -with modifications. Functions in non-free must be freely -redistributable for non-commercial use. Functions of unknown -license should not be included anywhere, since no license implies -default license implies no rights to redistribute.
--- a/nonfree/Makefile Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -sinclude ../Makeconf - - -ifeq ($(PKGDIR),) -PKGDIR = ../packages/nonfree -else -PKGDIR += /nonfree -endif - -# Determine which subdirectories are to be installed. Of those, determine -# which have their own Makefile. -SUBMAKEDIRS = $(dir $(wildcard */Makefile)) -NOINSTALLDIRS = $(dir $(wildcard */NOINSTALL)) -BUILDDIRS = $(filter-out $(NOINSTALLDIRS), $(SUBMAKEDIRS)) -INSTALLDIRS = $(filter-out $(BUILDDIRS), $(filter-out .svn/ $(NOINSTALLDIRS), $(dir $(wildcard */.)))) - -.PHONY: all package clean distclean $(SUBMAKEDIRS) - -ifdef OCTAVE_FORGE -all: - -package: checkpkgdir $(patsubst %,dopkg/%,$(BUILDDIRS)) $(patsubst %,dopkg2/%,$(INSTALLDIRS)) - -checkpkgdir: - @if [ ! -d $(PKGDIR) ]; then mkdir $(PKGDIR); fi - -pkg/%: - @if [ -e "$(opkg)/Makefile" ]; then \ - $(MAKE) dopkg/$(opkg); \ - else \ - $(MAKE) dopkg2/$(opkg); \ - fi - -dopkg/%: - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) pre-pkg - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) pkg/$(opkg) - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) post-pkg - -dopkg2/%: - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) pre-pkg -f../../Makeconf - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) pkg/$(opkg) -f../../Makeconf - @$(MAKE) PKGDIR=$(PKGDIR) -C $(opkg) post-pkg -f../../Makeconf -else -package: - @echo not yet configured - -all: - @echo not yet configured -endif - -clean: $(SUBMAKEDIRS) - -# Propogate make to the subdirectory if the goal is a valid target -# in the subdirectory Makefile. -$(SUBMAKEDIRS): - @echo Processing extra/$@ - @if test -z "$(MAKECMDGOALS)" ; then \ - cd $@ && $(MAKE) ; \ - elif grep "^$(MAKECMDGOALS) *[:]" $@Makefile >/dev/null; then \ - cd $@ && $(MAKE) $(MAKECMDGOALS) ; \ - fi
--- a/nonfree/spline-gsvspl/COPYING Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -MEMO: GCVSPL software package - -(C) COPYRIGHT 1985, 1986: H.J. Woltring - Philips Medical Systems Division, Eindhoven - University of Nijmegen (The Netherlands) - -DATE: 1986-05-12 - -NB: This software is copyrighted, and may be copied for excercise, -study and use without authorization from the copyright owner(s), in -compliance with paragraph 16b of the Dutch Copyright Act of 1912 -("Auteurswet 1912"). Within the constraints of this legislation, all -forms of academic and research-oriented excercise, study, and use are -allowed, including any necessary modifications. Copying and use as -object for commercial exploitation are not allowed without permission -of the copyright owners, including those upon whose work the package -is based. - -[see also: -http://isb.ri.ccf.org/biomch-l/archives/biomch-l-1994-06/00093.html -http://www.utc.edu/Human-Movement/3-d/herman.htm -Kai Habel]
--- a/nonfree/spline-gsvspl/DESCRIPTION Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -Name: spline-gcvspl -Version: 1.0.8 -Date: 2009-05-03 -Author: Joerg Specht -Maintainer: Joerg Specht -Title: Spline function based on GCVSPL package -Description: B-spline data smoothing using generalized cross-validation and mean squared prediction or explicit user smoothing -Categories: Interpolation -Depends: octave (>= 2.9.7) -License: Non Commercial Use Only -Url: http://octave.sf.net
--- a/nonfree/spline-gsvspl/INDEX Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -analysis >> Data analysis -Spline functions - csaps - gcvspl
--- a/nonfree/spline-gsvspl/inst/csaps.m Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,85 +0,0 @@ -## -*- texinfo -*- -## @deftypefn{Function File}{[@var{yi}, @var{p}] =} csaps(@var{x}, @var{y}, @var{p}, @var{xi}, @var{w}=[]) -## @deftypefnx{Function File}{[@var{pp}, @var{p}] =} csaps(@var{x}, @var{y}, @var{p}=-1, [], @var{w}=[]) -## -## Cubic spline approximation (smoothing)@* -## approximate [x,y] weighted w at xi -## -## @table @asis -## @item @var{p}<0 -## automatic smoothing -## @item @var{p}=0 -## maximum smoothing: straight line -## @item @var{p}=1 -## no smoothing: interpolation -## @end table -## -## @seealso{csapi, ppval, gcvspl} -## @end deftypefn - -## Author: Joerg Specht - -## This program is granted to the public domain. -## -## THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND -## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE -## FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS -## OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -## HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -## LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -## OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF -## SUCH DAMAGE. - -function [r,p]=csaps(x,y,p,xi,w) - if(nargin < 5) - w = []; - if(nargin < 4) - xi = []; - if(nargin < 3) - p = -1; - endif - endif - endif - - if(columns(x) > 1) - x = x.'; - y = y.'; - w = w.'; - endif - - [x,i] = sort(x); - y = y(i); - - if(p < 0) - md = 2; - p = 1; - else - md = 1; - ## csaps uses p=[0..1]; gcvspl uses p=[Inf..0] - p = 1/p - 1; - endif - if(length(xi) > 0) - if(columns(xi) > 1) - transposed = 1; - xi = xi.'; - else - transposed = 0; - endif - [yi,wk] = gcvspl(x,y,xi,w,[],2,md,p); - if(transposed) - r = yi.'; - else - r = yi; - endif - else - [D,wk] = gcvspl(x,y,x(1:length(x)-1),w,[],2,md,p,[3,2,1,0]); - ## gcvspl() produces derivates D, mkpp() needs coefficients P - pp = mkpp(x,[D(:,1)/6,D(:,2)/2,D(:,3),D(:,4)]); - r = pp; - endif - p = 1/(wk(4)+1); - -endfunction
--- a/nonfree/spline-gsvspl/src/.svnignore Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -autom4te.cache -configure -
--- a/nonfree/spline-gsvspl/src/Makeconf.in Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ - -## Makeconf is automatically generated from Makeconf.base and Makeconf.add -## in the various subdirectories. To regenerate, use ./autogen.sh to -## create a new ./Makeconf.in, then use ./configure to generate a new -## Makeconf. - -OCTAVE_FORGE = 1 - -SHELL = @SHELL@ - -canonical_host_type = @canonical_host_type@ -prefix = @prefix@ -exec_prefix = @exec_prefix@ -bindir = @bindir@ -mandir = @mandir@ -libdir = @libdir@ -datadir = @datadir@ -infodir = @infodir@ -includedir = @includedir@ -datarootdir = @datarootdir@ -INSTALL = @INSTALL@ -INSTALL_PROGRAM = @INSTALL_PROGRAM@ -INSTALL_SCRIPT = @INSTALL_SCRIPT@ -INSTALL_DATA = @INSTALL_DATA@ -INSTALLOCT=octinst.sh - -DESTDIR = - -RANLIB = @RANLIB@ -STRIP = @STRIP@ -LN_S = @LN_S@ - -AWK = @AWK@ - -# Most octave programs will be compiled with $(MKOCTFILE). Those which -# cannot use mkoctfile directly can request the flags that mkoctfile -# would use as follows: -# FLAG = $(shell $(MKOCTFILE) -p FLAG) -# The following flags are for compiling programs that are independent -# of Octave. How confusing. -CC = @CC@ -CFLAGS = @CFLAGS@ -CPPFLAGS = @CPPFLAGS@ -CPICFLAG = @CPICFLAG@ -CXX = @CXX@ -CXXFLAGS = @CXXFLAGS@ -CXXPICFLAG = @CXXPICFLAG@ -F77 = @F77@ -FFLAGS = @FFLAGS@ -FPICFLAG = @FPICFLAG@ - -OCTAVE = @OCTAVE@ -OCTAVE_VERSION = @OCTAVE_VERSION@ -MKOCTFILE = @MKOCTFILE@ -DHAVE_OCTAVE_$(ver) -v -SHLEXT = @SHLEXT@ - -ver = @ver@ -MPATH = @mpath@ -OPATH = @opath@ -XPATH = @xpath@ -ALTMPATH = @altmpath@ -ALTOPATH = @altopath@ - -%.o: %.c ; $(MKOCTFILE) -c $< -%.o: %.f ; $(MKOCTFILE) -c $< -%.o: %.cc ; $(MKOCTFILE) -c $< -%.oct: %.cc ; $(MKOCTFILE) $<
--- a/nonfree/spline-gsvspl/src/Makefile Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -sinclude ./Makeconf - -all: gcvspl.oct - -gcvspl.oct: gcvspl.cc gcvsplf.f - $(MKOCTFILE) -v gcvspl.cc gcvsplf.f - -clean: ; -$(RM) *.o octave-core core *.oct *~ -
--- a/nonfree/spline-gsvspl/src/autogen.sh Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -#! /bin/sh - -## Generate ./configure -rm -f configure.in -echo "dnl --- DO NOT EDIT --- Automatically generated by autogen.sh" > configure.in -cat configure.base >> configure.in -cat <<EOF >> configure.in - AC_OUTPUT(\$CONFIGURE_OUTPUTS) - dnl XXX FIXME XXX chmod is not in autoconf's list of portable functions - - echo " " - echo " \"\\\$prefix\" is \$prefix" - echo " \"\\\$exec_prefix\" is \$exec_prefix" - AC_MSG_RESULT([\$STATUS_MSG - -find . -name NOINSTALL -print # shows which toolboxes won't be installed -]) -EOF - -autoconf configure.in > configure.tmp -if [ diff configure.tmp configure > /dev/null 2>&1 ]; then - rm -f configure.tmp; -else - mv -f configure.tmp configure - chmod 0755 configure -fi -rm -f configure.in
--- a/nonfree/spline-gsvspl/src/configure.base Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,354 +0,0 @@ -dnl The configure script is generated by autogen.sh from configure.base -dnl and the various configure.add files in the source tree. Edit -dnl configure.base and reprocess rather than modifying ./configure. - -dnl autoconf 2.13 certainly doesn't work! What is the minimum requirement? -AC_PREREQ(2.2) - -AC_INIT(configure.base) - -PACKAGE=octave-forge -MAJOR_VERSION=0 -MINOR_VERSION=1 -PATCH_LEVEL=0 - -dnl Kill caching --- this ought to be the default -define([AC_CACHE_LOAD], )dnl -define([AC_CACHE_SAVE], )dnl - -dnl uncomment to put support files in another directory -dnl AC_CONFIG_AUX_DIR(admin) - -VERSION=$MAJOR_VERSION.$MINOR_VERSION.$PATCH_LEVEL -AC_SUBST(PACKAGE) -AC_SUBST(VERSION) - -dnl need to find admin files, so keep track of the top dir. -TOPDIR=`pwd` -AC_SUBST(TOPDIR) - -dnl if mkoctfile doesn't work, then we need the following: -dnl AC_PROG_CXX -dnl AC_PROG_F77 - -dnl Need C compiler regardless so define it in a way that -dnl makes autoconf happy and we can override whatever we -dnl need with mkoctfile -p. -dnl XXX FIXME XXX should use mkoctfile to get CC and CFLAGS -AC_PROG_CC - -dnl XXX FIXME XXX need tests for -p -c -s in mkoctfile. - -dnl ******************************************************************* -dnl Sort out mkoctfile version number and install paths - -dnl XXX FIXME XXX latest octave has octave-config so we don't -dnl need to discover things here. Doesn't have --exe-site-dir -dnl but defines --oct-site-dir and --m-site-dir - -dnl Check for mkoctfile -AC_CHECK_PROG(MKOCTFILE,mkoctfile,mkoctfile) -test -z "$MKOCTFILE" && AC_MSG_WARN([no mkoctfile found on path]) - -AC_SUBST(ver) -AC_SUBST(subver) -AC_SUBST(mpath) -AC_SUBST(opath) -AC_SUBST(xpath) -AC_SUBST(altpath) -AC_SUBST(altmpath) -AC_SUBST(altopath) - -AC_ARG_WITH(path, - [ --with-path install path prefix], - [ path=$withval ]) -AC_ARG_WITH(mpath, - [ --with-mpath override path for m-files], - [mpath=$withval]) -AC_ARG_WITH(opath, - [ --with-opath override path for oct-files], - [opath=$withval]) -AC_ARG_WITH(xpath, - [ --with-xpath override path for executables], - [xpath=$withval]) -AC_ARG_WITH(altpath, - [ --with-altpath alternative functions install path prefix], - [ altpath=$withval ]) -AC_ARG_WITH(altmpath, - [ --with-altmpath override path for alternative m-files], - [altmpath=$withval]) -AC_ARG_WITH(altopath, - [ --with-altopath override path for alternative oct-files], - [altopath=$withval]) - -if test -n "$path" ; then - test -z "$mpath" && mpath=$path - test -z "$opath" && opath=$path/oct - test -z "$xpath" && xpath=$path/bin - test -z "$altpath" && altpath=$path-alternatives -fi - -if test -n "$altpath" ; then - test -z "$altmpath" && altmpath=$altpath - test -z "$altopath" && altopath=$altpath/oct -fi - -dnl Don't query if path/ver are given in the configure environment -#if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$altmpath" || test -z "$altopath" || test -z "$ver" ; then -if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$ver" ; then - dnl Construct program to get mkoctfile version and local install paths - cat > conftest.cc <<EOF -#include <octave/config.h> -#include <octave/version.h> -#include <octave/defaults.h> - -#define INFOV "\nINFOV=" OCTAVE_VERSION "\n" - -#define INFOH "\nINFOH=" OCTAVE_CANONICAL_HOST_TYPE "\n" - -#ifdef OCTAVE_LOCALVERFCNFILEDIR -# define INFOM "\nINFOM=" OCTAVE_LOCALVERFCNFILEDIR "\n" -#else -# define INFOM "\nINFOM=" OCTAVE_LOCALFCNFILEPATH "\n" -#endif - -#ifdef OCTAVE_LOCALVEROCTFILEDIR -# define INFOO "\nINFOO=" OCTAVE_LOCALVEROCTFILEDIR "\n" -#else -# define INFOO "\nINFOO=" OCTAVE_LOCALOCTFILEPATH "\n" -#endif - -#ifdef OCTAVE_LOCALVERARCHLIBDIR -# define INFOX "\nINFOX=" OCTAVE_LOCALVERARCHLIBDIR "\n" -#else -# define INFOX "\nINFOX=" OCTAVE_LOCALARCHLIBDIR "\n" -#endif - -const char *infom = INFOM; -const char *infoo = INFOO; -const char *infox = INFOX; -const char *infoh = INFOH; -const char *infov = INFOV; -EOF - - dnl Compile program perhaps with a special version of mkoctfile - $MKOCTFILE conftest.cc || AC_MSG_ERROR(Could not run $MKOCTFILE) - - dnl Strip the config info from the compiled file - eval `strings conftest.o | grep "^INFO.=" | sed -e "s,//.*$,,"` - rm -rf conftest* - - dnl set the appropriate variables if they are not already set - ver=`echo $INFOV | sed -e "s/\.//" -e "s/\..*$//"` - subver=`echo $INFOV | sed -e "[s/^[^.]*[.][^.]*[.]//]"` - alt_mbase=`echo $INFOM | sed -e "[s,\/[^\/]*$,,]"` - alt_obase=`echo $INFOO | sed -e "[s,/site.*$,/site,]"` - test -z "$mpath" && mpath=$INFOM/octave-forge - test -z "$opath" && opath=$INFOO/octave-forge - test -z "$xpath" && xpath=$INFOX - test -z "$altmpath" && altmpath=$alt_mbase/octave-forge-alternatives/m - test -z "$altopath" && altopath=$alt_obase/octave-forge-alternatives/oct/$INFOH -fi - -dnl ******************************************************************* - -dnl XXX FIXME XXX Should we allow the user to override these? -dnl Do we even need them? The individual makefiles can call mkoctfile -p -dnl themselves, so the only reason to keep them is for configure, and -dnl for those things which are not built using mkoctfile (e.g., aurecord) -dnl but it is not clear we should be using octave compile flags for those. - -dnl C compiler and flags -AC_MSG_RESULT([retrieving compile and link flags from $MKOCTFILE]) -CC=`$MKOCTFILE -p CC` -CFLAGS=`$MKOCTFILE -p CFLAGS` -CPPFLAGS=`$MKOCTFILE -p CPPFLAGS` -CPICFLAG=`$MKOCTFILE -p CPICFLAG` -LDFLAGS=`$MKOCTFILE -p LDFLAGS` -LIBS=`$MKOCTFILE -p LIBS` -AC_SUBST(CC) -AC_SUBST(CFLAGS) -AC_SUBST(CPPFLAGS) -AC_SUBST(CPICFLAG) - -dnl Fortran compiler and flags -F77=`$MKOCTFILE -p F77` -FFLAGS=`$MKOCTFILE -p FFLAGS` -FPICFLAG=`$MKOCTFILE -p FPICFLAG` -AC_SUBST(F77) -AC_SUBST(FFLAGS) -AC_SUBST(FPICFLAG) - -dnl C++ compiler and flags -CXX=`$MKOCTFILE -p CXX` -CXXFLAGS=`$MKOCTFILE -p CXXFLAGS` -CXXPICFLAG=`$MKOCTFILE -p CXXPICFLAG` -AC_SUBST(CXX) -AC_SUBST(CXXFLAGS) -AC_SUBST(CXXPICFLAG) - -dnl ******************************************************************* - -dnl Check for features of your version of mkoctfile. -dnl All checks should be designed so that the default -dnl action if the tests are not performed is to do whatever -dnl is appropriate for the most recent version of Octave. - -dnl Define the following macro: -dnl OF_CHECK_LIB(lib,fn,true,false,helpers) -dnl This is just like AC_CHECK_LIB, but it doesn't update LIBS -AC_DEFUN(OF_CHECK_LIB, -[save_LIBS="$LIBS" -AC_CHECK_LIB($1,$2,$3,$4,$5) -LIBS="$save_LIBS" -]) - -dnl Define the following macro: -dnl TRY_MKOCTFILE(msg,program,action_if_true,action_if_false) -dnl -AC_DEFUN(TRY_MKOCTFILE, -[AC_MSG_CHECKING($1) -cat > conftest.cc << EOF -#include <octave/config.h> -$2 -EOF -ac_try="$MKOCTFILE -c conftest.cc" -if AC_TRY_EVAL(ac_try) ; then - AC_MSG_RESULT(yes) - $3 -else - AC_MSG_RESULT(no) - $4 -fi -]) - -dnl -dnl Check if F77_FUNC works with MKOCTFILE -dnl -TRY_MKOCTFILE([for F77_FUNC], -[int F77_FUNC (hello, HELLO) (const int &n);],, -[MKOCTFILE="$MKOCTFILE -DF77_FUNC=F77_FCN"]) - -dnl ********************************************************** - -dnl Evaluate an expression in octave -dnl -dnl OCTAVE_EVAL(expr,var) -> var=expr -dnl -AC_DEFUN(OCTAVE_EVAL, -[AC_MSG_CHECKING([for $1 in Octave]) -$2=`echo "disp($1)" | $OCTAVE -qf` -AC_MSG_RESULT($$2) -AC_SUBST($2) -]) - -dnl Check status of an octave variable -dnl -dnl OCTAVE_CHECK_EXIST(variable,action_if_true,action_if_false) -dnl -AC_DEFUN(OCTAVE_CHECK_EXIST, -[AC_MSG_CHECKING([for $1 in Octave]) -if test `echo 'disp(exist("$1"))' | $OCTAVE -qf`X != 0X ; then - AC_MSG_RESULT(yes) - $2 -else - AC_MSG_RESULT(no) - $3 -fi -]) - -dnl should check that $(OCTAVE) --version matches $(MKOCTFILE) --version -AC_CHECK_PROG(OCTAVE,octave,octave) -OCTAVE_EVAL(OCTAVE_VERSION,OCTAVE_VERSION) - -dnl grab canonical host type so we can write system specific install stuff -OCTAVE_EVAL(octave_config_info('canonical_host_type'),canonical_host_type) - -dnl grab SHLEXT from octave config -OCTAVE_EVAL(octave_config_info('SHLEXT'),SHLEXT) - -AC_PROG_LN_S - -AC_PROG_RANLIB - -dnl Use $(COPY_FLAGS) to set options for cp when installing .oct files. -COPY_FLAGS="-Rfp" -case "$canonical_host_type" in - *-*-linux*) - COPY_FLAGS="-fdp" - ;; -esac -AC_SUBST(COPY_FLAGS) - -dnl Use $(STRIP) in the makefile to strip executables. If not found, -dnl STRIP expands to ':', which in the makefile does nothing. -dnl Don't need this for .oct files since mkoctfile handles them directly -STRIP=${STRIP-strip} -AC_CHECK_PROG(STRIP,$STRIP,$STRIP,:) - -dnl Strip on windows, don't strip on Mac OS/X or IRIX -dnl For the rest, you can force strip using MKOCTFILE="mkoctfile -s" -dnl or avoid strip using STRIP=: before ./configure -case "$canonical_host_type" in - powerpc-apple-darwin*|*-sgi-*) - STRIP=: - ;; - *-cygwin-*|*-mingw-*) - MKOCTFILE="$MKOCTFILE -s" - ;; -esac - -AC_DEFUN(AC_CHECK_QHULL_VERSION, -[AC_MSG_CHECKING([for qh_qhull in -lqhull with qh_version]) -AC_CACHE_VAL(ac_cv_lib_qhull_version, -[changequote(, )dnl -cat > conftest.c <<EOF -#include <stdio.h> -char qh_version[] = "version"; -char qh_qhull(); -int -main(argc, argv) - int argc; - char *argv[]; -{ - qh_qhull(); - return 0; -} -EOF -changequote([, ])dnl -ac_try="${CC-cc} $CFLAGS $CPPFLAGS $LDFLAGS conftest.c -o conftest -lqhull $LIBS" -if AC_TRY_EVAL(ac_try) && test -s conftest ; then - ac_cv_lib_qhull_version=yes -else - ac_cv_lib_qhull_version=no -fi -rm -f conftest.c conftest.o conftest -])dnl -if test "$ac_cv_lib_qhull_version" = "yes"; then - AC_MSG_RESULT(yes) - ifelse([$1], , , [$1]) -else - AC_MSG_RESULT(no) - ifelse([$2], , , [$2]) -fi -]) - -CONFIGURE_OUTPUTS="Makeconf" -STATUS_MSG=" -octave commands will install into the following directories: - m-files: $mpath - oct-files: $opath - binaries: $xpath -alternatives: - m-files: $altmpath - oct-files: $altopath - -shell commands will install into the following directories: - binaries: $bindir - man pages: $mandir - libraries: $libdir - headers: $includedir - -octave-forge is configured with - octave: $OCTAVE (version $OCTAVE_VERSION) - mkoctfile: $MKOCTFILE for Octave $subver"
--- a/nonfree/spline-gsvspl/src/gcvspl.cc Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,361 +0,0 @@ -/* -** Author: Joerg Specht -** -** This program is granted to the public domain. -** -** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND -** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE -** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS -** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF -** SUCH DAMAGE. -*/ - -// use: mkoctfile gcvspl.cc gcvsplf.f -#define NDEBUG - -#include <octave/config.h> -#include <octave/defun-dld.h> -#include <octave/error.h> -#include <octave/help.h> -#include <octave/oct-obj.h> -#include <octave/pager.h> -#include <octave/symtab.h> -#include <octave/variables.h> - -extern "C" { - /* by f2c -A ...: */ - extern int gcvspl_(double *x, double *y, int *ny, - double *wx, double *wy, int *m, int *n, - int *k, int *md, double *val, double *c, - int *nc, double *wk, int *ier); - extern double splder_(int *ider, int *m, int *n, double *t, - double *x, double *c, int *l, double *q); -} - -#ifndef NDEBUG -#ifndef DEBUGLEVEL -#define DEBUGLEVEL 255 -#endif -static int do_debug(const char *fmt, ...) { - va_list args; - fprintf(stderr, "debug: "); - va_start(args, fmt); - vfprintf(stderr, fmt, args); - va_end(args); - fprintf(stderr, "\n"); - fflush(NULL); - return 1; -} -#define DEBUG(level, fmtetc) ((level) & DEBUGLEVEL) && do_debug fmtetc -#else -#define DEBUG(level, fmtetc) -#endif - -#define ASSERT(what) \ - if(!(what)) { \ - (*current_liboctave_error_handler)( \ - "gcvspl: assertion failed: "#what); \ - return retval; \ - } - -DEFUN_DLD(gcvspl, args, nargout, -"B-spline data smoothing subroutine.\n\n\ -[yf,wk]=gcvspl(x(n,1),y(n,k), xf(nf,1)=x, wx(n,1)=[],wy(1,k)=[], m=2,md=2,val=1, ider=[0])\n\ -\n\ -uses GCVSPL.FOR, 1986-05-12\n\ -from http://www.netlib.org/gcv/index.html\n\ -for B-spline data smoothing using generalized cross-validation\n\ - and mean squared prediction or explicit user smoothing\n\ -by H.J. Woltring, University of Nijmegen,\n\ - Philips Medical Systems, Eindhoven (The Netherlands)\n\ -\n\ -Purpose:\n\ - Natural B-spline data smoothing subroutine, using the Generali-\n\ - zed Cross-Validation and Mean-Squared Prediction Error Criteria\n\ - of Craven & Wahba (1979). Alternatively, the amount of smoothing\n\ - can be given explicitly, or it can be based on the effective\n\ - number of degrees of freedom in the smoothing process as defined\n\ - by Wahba (1980). The model assumes uncorrelated, additive noise\n\ - and essentially smooth, underlying functions. The noise may be\n\ - non-stationary, and the independent co-ordinates may be spaced\n\ - non-equidistantly. Multiple datasets, with common independent\n\ - variables and weight factors are accomodated.\n\ - A full description of the package is provided in:\n\ - H.J. Woltring (1986), A FORTRAN package for generalized,\n\ - cross-validatory spline smoothing and differentiation.\n\ - Advances in Engineering Software 8(2):104-113\n\ -\n\ -Meaning of parameters:\n\ - X(N,1) Independent variables: strictly increasing knot\n\ - sequence, with X(I-1).lt.X(I), I=2,...,N.\n\ - Y(N,K) Input data to be smoothed (or interpolated).\n\ - XF(NF,1)Points where the function should be approximated.\n\ - WX(N,1) Weight factor array; WX(I) corresponds with\n\ - the relative inverse variance of point Y(I,*).\n\ - If no relative weighting information is\n\ - available, the WX(I) should be set to ONE.\n\ - All WX(I).gt.ZERO, I=1,...,N.\n\ - WY(1,K) Weight factor array; WY(J) corresponds with\n\ - the relative inverse variance of point Y(*,J).\n\ - If no relative weighting information is\n\ - available, the WY(J) should be set to ONE.\n\ - All WY(J).gt.ZERO, J=1,...,K.\n\ - NB: The effective weight for point Y(I,J) is\n\ - equal to WX(I)*WY(J).\n\ - M Half order of the required B-splines (spline\n\ - degree 2*M-1), with M.gt.0. The values M =\n\ - 1,2,3,4 correspond to linear, cubic, quintic,\n\ - and heptic splines, respectively. N.ge.2*M.\n\ - MD Optimization mode switch:\n\ - |MD| = 1: Prior given value for p in VAL\n\ - (VAL.ge.ZERO). This is the fastest\n\ - use of GCVSPL, since no iteration\n\ - is performed in p.\n\ - |MD| = 2: Generalized cross validation.\n\ - |MD| = 3: True predicted mean-squared error,\n\ - with prior given variance in VAL.\n\ - |MD| = 4: Prior given number of degrees of\n\ - freedom in VAL (ZERO.le.VAL.le.N-M).\n\ - After return from MD.ne.1, the same number of\n\ - degrees of freedom can be obtained, for identical\n\ - weight factors and knot positions, by selecting\n\ - |MD|=1, and by copying the value of p from WK(4)\n\ - into VAL. In this way, no iterative optimization\n\ - is required when processing other data in Y.\n\ - VAL Mode value, as described above under MD.\n\ - IDER Derivative order required, with 0.le.IDER\n\ - and IDER.le.2*M. If IDER.eq.0, the function\n\ - value is returned; otherwise, the IDER-th\n\ - derivative of the spline is returned.\n\ -\n\ -Return values:\n\ - YF(NF,1)Approximated values at XF.\n\ - WK(IWK) On normal exit, the first 6 values of WK are\n\ - assigned as follows:\n\ - WK(1) = Generalized Cross Validation value\n\ - WK(2) = Mean Squared Residual.\n\ - WK(3) = Estimate of the number of degrees of\n\ - freedom of the residual sum of squares\n\ - per dataset, with 0.lt.WK(3).lt.N-M.\n\ - WK(4) = Smoothing parameter p, multiplicative\n\ - with the splines' derivative constraint.\n\ - WK(5) = Estimate of the true mean squared error\n\ - (different formula for |MD| = 3).\n\ - WK(6) = Gauss-Markov error variance.\n\ -\n\ - If WK(4) --> 0 , WK(3) --> 0 , and an inter-\n\ - polating spline is fitted to the data (p --> 0).\n\ - A very small value > 0 is used for p, in order\n\ - to avoid division by zero in the GCV function.\n\ -\n\ - If WK(4) --> inf, WK(3) --> N-M, and a least-\n\ - squares polynomial of order M (degree M-1) is\n\ - fitted to the data (p --> inf). For numerical\n\ - reasons, a very high value is used for p.\n\ -\n\ - Upon return, the contents of WK can be used for\n\ - covariance propagation in terms of the matrices\n\ - B and WE: see the source listings. The variance\n\ - estimate for dataset J follows as WK(6)/WY(J).\n\ -\n\ -Remarks:\n\ - (1) GCVSPL calculates a natural spline of order 2*M (degree\n\ - 2*M-1) which smoothes or interpolates a given set of data\n\ - points, using statistical considerations to determine the\n\ - amount of smoothing required (Craven & Wahba, 1979). If the\n\ - error variance is a priori known, it should be supplied to\n\ - the routine in VAL, for |MD|=3. The degree of smoothing is\n\ - then determined to minimize an unbiased estimate of the true\n\ - mean squared error. On the other hand, if the error variance\n\ - is not known, one may select |MD|=2. The routine then deter-\n\ - mines the degree of smoothing to minimize the generalized\n\ - cross validation function. This is asymptotically the same\n\ - as minimizing the true predicted mean squared error (Craven &\n\ - Wahba, 1979). If the estimates from |MD|=2 or 3 do not appear\n\ - suitable to the user (as apparent from the smoothness of the\n\ - M-th derivative or from the effective number of degrees of\n\ - freedom returned in WK(3) ), the user may select an other\n\ - value for the noise variance if |MD|=3, or a reasonably large\n\ - number of degrees of freedom if |MD|=4. If |MD|=1, the proce-\n\ - dure is non-iterative, and returns a spline for the given\n\ - value of the smoothing parameter p as entered in VAL.\n\ -\n\ - (2) The number of arithmetic operations and the amount of\n\ - storage required are both proportional to N, so very large\n\ - datasets may be accomodated. The data points do not have\n\ - to be equidistant in the independant variable X or uniformly\n\ - weighted in the dependant variable Y. However, the data\n\ - points in X must be strictly increasing. Multiple dataset\n\ - processing (K.gt.1) is numerically more efficient dan\n\ - separate processing of the individual datasets (K.eq.1).\n\ -\n\ - (3) If |MD|=3 (a priori known noise variance), any value of\n\ - N.ge.2*M is acceptable. However, it is advisable for N-2*M\n\ - be rather large (at least 20) if |MD|=2 (GCV).\n\ -\n\ - (4) For |MD| > 1, GCVSPL tries to iteratively minimize the\n\ - selected criterion function. This minimum is unique for |MD|\n\ - = 4, but not necessarily for |MD| = 2 or 3. Consequently, \n\ - local optima rather that the global optimum might be found,\n\ - and some actual findings suggest that local optima might\n\ - yield more meaningful results than the global optimum if N\n\ - is small. Therefore, the user has some control over the\n\ - search procedure. If MD > 1, the iterative search starts\n\ - from a value which yields a number of degrees of freedom\n\ - which is approximately equal to N/2, until the first (local)\n\ - minimum is found via a golden section search procedure\n\ - (Utreras, 1980). If MD < -1, the value for p contained in\n\ - WK(4) is used instead. Thus, if MD = 2 or 3 yield too noisy\n\ - an estimate, the user might try |MD| = 1 or 4, for suitably\n\ - selected values for p or for the number of degrees of\n\ - freedom, and then run GCVSPL with MD = -2 or -3. The con-\n\ - tents of N, M, K, X, WX, WY, and WK are assumed unchanged\n\ - if MD < 0.\n\ -\n\ - (5) GCVSPL calculates the spline coefficient array C(N,K);\n\ - this array can be used to calculate the spline function\n\ - value and any of its derivatives up to the degree 2*M-1\n\ - at any argument T within the knot range, using subrou-\n\ - tines SPLDER and SEARCH, and the knot array X(N). Since\n\ - the splines are constrained at their Mth derivative, only\n\ - the lower spline derivatives will tend to be reliable\n\ - estimates of the underlying, true signal derivatives.\n\ -\n\ - (6) GCVSPL combines elements of subroutine CRVO5 by Utre-\n\ - ras (1980), subroutine SMOOTH by Lyche et al. (1983), and\n\ - subroutine CUBGCV by Hutchinson (1985). The trace of the\n\ - influence matrix is assessed in a similar way as described\n\ - by Hutchinson & de Hoog (1985). The major difference is\n\ - that the present approach utilizes non-symmetrical B-spline\n\ - design matrices as described by Lyche et al. (1983); there-\n\ - fore, the original algorithm by Erisman & Tinney (1975) has\n\ - been used, rather than the symmetrical version adopted by\n\ - Hutchinson & de Hoog.\n\ -\n\ -References:\n\ - P. Craven & G. Wahba (1979), Smoothing noisy data with\n\ - spline functions. Numerische Mathematik 31, 377-403.\n\ -\n\ - A.M. Erisman & W.F. Tinney (1975), On computing certain\n\ - elements of the inverse of a sparse matrix. Communications\n\ - of the ACM 18(3), 177-179.\n\ -\n\ - M.F. Hutchinson & F.R. de Hoog (1985), Smoothing noisy data\n\ - with spline functions. Numerische Mathematik 47(1), 99-106.\n\ -\n\ - M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO Division of\n\ - Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601,\n\ - Australia.\n\ -\n\ - T. Lyche, L.L. Schumaker, & K. Sepehrnoori (1983), Fortran\n\ - subroutines for computing smoothing and interpolating natural\n\ - splines. Advances in Engineering Software 5(1), 2-5.\n\ -\n\ - F. Utreras (1980), Un paquete de programas para ajustar curvas\n\ - mediante funciones spline. Informe Tecnico MA-80-B-209, Depar-\n\ - tamento de Matematicas, Faculdad de Ciencias Fisicas y Matema-\n\ - ticas, Universidad de Chile, Santiago.\n\ -\n\ - Wahba, G. (1980). Numerical and statistical methods for mildly,\n\ - moderately and severely ill-posed problems with noisy data.\n\ - Technical report nr. 595 (February 1980). Department of Statis-\n\ - tics, University of Madison (WI), U.S.A.\ -") -{ - octave_value_list retval; - DEBUG(1, ("A")); - int nargs = args.length(); - ASSERT(nargs >= 2 && nargs <= 9); - Matrix x=args(0).matrix_value(); - Matrix y=args(1).matrix_value(); - Matrix x_target = nargs > 2 ? args(2).matrix_value() : x; - // int ny=y.rows(); --> ny=n - ASSERT(!error_state); - int n=x.rows(); - int k=y.columns(); - Matrix wx = (nargs > 3 && !args(3).is_zero_by_zero()) - ? args(3).matrix_value() : Matrix(n,1,1.0); - Matrix wy = (nargs > 4 && !args(4).is_zero_by_zero()) - ? args(4).matrix_value() : Matrix(1,k,1.0); - int m = nargs > 5 ? (int)rint(args(5).double_value()) : 2; - int md = nargs > 6 ? (int)rint(args(6).double_value()) : 2; - double val = nargs > 7 ? args(7).double_value() : 1.0; - Matrix mider = nargs > 8 ? args(8).matrix_value() : Matrix(1,1,0.0); - - // int nc=n; - OCTAVE_LOCAL_BUFFER(double,c,n*k); // Matrix c(n,k); - OCTAVE_LOCAL_BUFFER(double,wk,6*(n*m+1)+n);// work array, return status - int ier; - - DEBUG(1, ("B")); - ASSERT(!error_state); - ASSERT(x.columns() == 1); - ASSERT(y.rows() == n); - ASSERT(wx.rows() == n && wx.columns() == 1); - ASSERT(wy.rows() == 1 && wy.columns() == k); - ASSERT(m > 0); - ASSERT(n >= 2*m); - ASSERT(k >= 1); - ASSERT(md >= 1 && md <= 4); - ASSERT(val >= 0); - ASSERT(mider.rows() == 1); - int nider = mider.columns(); - OCTAVE_LOCAL_BUFFER(int,ider,nider); - for(int i=0; i<nider; i++) { - ider[i] = (int)rint(mider.xelem(0,i)); - ASSERT(0 <= ider[i] && ider[i] <= 2*m); - } - - DEBUG(2, - ("gcvspl_(x,y,ny=%d,wx,wy,m=%d,n=%d,k=%d,md=%d,val=%g,c,nc=%d,wk,ier)", - n, m, n, k, md, val, n)); - double *x_fortran = x.fortran_vec(); - gcvspl_(x_fortran, y.fortran_vec(), &n, - wx.fortran_vec(), wy.fortran_vec(), &m, &n, &k, - &md, &val, c, &n, wk, &ier); - - if(ier != 0) { - (*current_liboctave_error_handler) - (ier==1 ? "M<=0 || N<2*M" - : ier==2 ? "knots not sorted or negative weight" - : ier==3 ? "wrong mode parameter or value" - : "unknown error"); - return retval; - } - - DEBUG(1, ("D")); - int l=0; // index in x, just simplifies search - OCTAVE_LOCAL_BUFFER(double,q,2*m);// work array - int nxt=x_target.rows(); - ASSERT(x_target.columns() == 1); - Matrix y_target(nxt,k*nider); - DEBUG(1, ("E")); - double *xf=x_target.fortran_vec(); - double *yf=y_target.fortran_vec(); - for(int i=0; i<nxt; i++) { // next point - double xt=xf[i]; - for(int j=0; j<k; j++) // next curve - for(int o=0; o<nider; o++) // next derivate - yf[nxt*(nider*j+o)+i] = splder_(&ider[o], &m, &n, &xt, x_fortran, c+n*j, &l, q); - } - - DEBUG(1, ("F")); - DEBUG(2, ("nargout=%d", nargout)); - retval(0) = y_target; - if(nargout > 1) { - DEBUG(1, ("G")); - RowVector wkv(6); - double *wkvf=wkv.fortran_vec(); - for(int i=0; i<6; i++) - wkvf[i] = wk[i]; - retval(1) = wkv; - } - DEBUG(1, ("H")); - return retval; -}
--- a/nonfree/spline-gsvspl/src/gcvsplf.f Sun Jan 27 15:12:41 2013 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1445 +0,0 @@ -C from: http://www.netlib.org/gcv/index.html -C "for B-spline data smoothing using generalized cross-validation -C and mean squared prediction or explicit user smoothing" -C "by H.J. Woltring, University of Nijmegen, -C Philips Medical Systems, Eindhoven (The Netherlands)" -C -C*********************************************************************** -C -C GCVSPL.FOR, 1986-05-12 -C -C*********************************************************************** -C -C SUBROUTINE GCVSPL (REAL*8) -C -C Purpose: -C ******* -C -C Natural B-spline data smoothing subroutine, using the Generali- -C zed Cross-Validation and Mean-Squared Prediction Error Criteria -C of Craven & Wahba (1979). Alternatively, the amount of smoothing -C can be given explicitly, or it can be based on the effective -C number of degrees of freedom in the smoothing process as defined -C by Wahba (1980). The model assumes uncorrelated, additive noise -C and essentially smooth, underlying functions. The noise may be -C non-stationary, and the independent co-ordinates may be spaced -C non-equidistantly. Multiple datasets, with common independent -C variables and weight factors are accomodated. -C -C -C Calling convention: -C ****************** -C -C CALL GCVSPL ( X, Y, NY, WX, WY, M, N, K, MD, VAL, C, NC, WK, IER ) -C -C Meaning of parameters: -C ********************* -C -C X(N) ( I ) Independent variables: strictly increasing knot -C sequence, with X(I-1).lt.X(I), I=2,...,N. -C Y(NY,K) ( I ) Input data to be smoothed (or interpolated). -C NY ( I ) First dimension of array Y(NY,K), with NY.ge.N. -C WX(N) ( I ) Weight factor array; WX(I) corresponds with -C the relative inverse variance of point Y(I,*). -C If no relative weighting information is -C available, the WX(I) should be set to ONE. -C All WX(I).gt.ZERO, I=1,...,N. -C WY(K) ( I ) Weight factor array; WY(J) corresponds with -C the relative inverse variance of point Y(*,J). -C If no relative weighting information is -C available, the WY(J) should be set to ONE. -C All WY(J).gt.ZERO, J=1,...,K. -C NB: The effective weight for point Y(I,J) is -C equal to WX(I)*WY(J). -C M ( I ) Half order of the required B-splines (spline -C degree 2*M-1), with M.gt.0. The values M = -C 1,2,3,4 correspond to linear, cubic, quintic, -C and heptic splines, respectively. -C N ( I ) Number of observations per dataset, with N.ge.2*M. -C K ( I ) Number of datasets, with K.ge.1. -C MD ( I ) Optimization mode switch: -C |MD| = 1: Prior given value for p in VAL -C (VAL.ge.ZERO). This is the fastest -C use of GCVSPL, since no iteration -C is performed in p. -C |MD| = 2: Generalized cross validation. -C |MD| = 3: True predicted mean-squared error, -C with prior given variance in VAL. -C |MD| = 4: Prior given number of degrees of -C freedom in VAL (ZERO.le.VAL.le.N-M). -C MD < 0: It is assumed that the contents of -C X, W, M, N, and WK have not been -C modified since the previous invoca- -C tion of GCVSPL. If MD < -1, WK(4) -C is used as an initial estimate for -C the smoothing parameter p. -C Other values for |MD|, and inappropriate values -C for VAL will result in an error condition, or -C cause a default value for VAL to be selected. -C After return from MD.ne.1, the same number of -C degrees of freedom can be obtained, for identical -C weight factors and knot positions, by selecting -C |MD|=1, and by copying the value of p from WK(4) -C into VAL. In this way, no iterative optimization -C is required when processing other data in Y. -C VAL ( I ) Mode value, as described above under MD. -C C(NC,K) ( O ) Spline coefficients, to be used in conjunction -C with function SPLDER. NB: the dimensions of C -C in GCVSPL and in SPLDER are different -C only a single column of C(N,K) is needed, and the -C proper column C(1,J), with J=1...K should be used -C when calling SPLDER. -C NC ( I ) First dimension of array C(NC,K), NC.ge.N. -C WK(IWK) (I/W/O) Work vector, with length IWK.ge.6*(N*M+1)+N. -C On normal exit, the first 6 values of WK are -C assigned as follows: -C -C WK(1) = Generalized Cross Validation value -C WK(2) = Mean Squared Residual. -C WK(3) = Estimate of the number of degrees of -C freedom of the residual sum of squares -C per dataset, with 0.lt.WK(3).lt.N-M. -C WK(4) = Smoothing parameter p, multiplicative -C with the splines' derivative constraint. -C WK(5) = Estimate of the true mean squared error -C (different formula for |MD| = 3). -C WK(6) = Gauss-Markov error variance. -C -C If WK(4) --> 0 , WK(3) --> 0 , and an inter- -C polating spline is fitted to the data (p --> 0). -C A very small value > 0 is used for p, in order -C to avoid division by zero in the GCV function. -C -C If WK(4) --> inf, WK(3) --> N-M, and a least- -C squares polynomial of order M (degree M-1) is -C fitted to the data (p --> inf). For numerical -C reasons, a very high value is used for p. -C -C Upon return, the contents of WK can be used for -C covariance propagation in terms of the matrices -C B and WE: see the source listings. The variance -C estimate for dataset J follows as WK(6)/WY(J). -C -C IER ( O ) Error parameter: -C -C IER = 0: Normal exit -C IER = 1: M.le.0 .or. N.lt.2*M -C IER = 2: Knot sequence is not strictly -C increasing, or some weight -C factor is not positive. -C IER = 3: Wrong mode parameter or value. -C -C Remarks: -C ******* -C -C (1) GCVSPL calculates a natural spline of order 2*M (degree -C 2*M-1) which smoothes or interpolates a given set of data -C points, using statistical considerations to determine the -C amount of smoothing required (Craven & Wahba, 1979). If the -C error variance is a priori known, it should be supplied to -C the routine in VAL, for |MD|=3. The degree of smoothing is -C then determined to minimize an unbiased estimate of the true -C mean squared error. On the other hand, if the error variance -C is not known, one may select |MD|=2. The routine then deter- -C mines the degree of smoothing to minimize the generalized -C cross validation function. This is asymptotically the same -C as minimizing the true predicted mean squared error (Craven & -C Wahba, 1979). If the estimates from |MD|=2 or 3 do not appear -C suitable to the user (as apparent from the smoothness of the -C M-th derivative or from the effective number of degrees of -C freedom returned in WK(3) ), the user may select an other -C value for the noise variance if |MD|=3, or a reasonably large -C number of degrees of freedom if |MD|=4. If |MD|=1, the proce- -C dure is non-iterative, and returns a spline for the given -C value of the smoothing parameter p as entered in VAL. -C -C (2) The number of arithmetic operations and the amount of -C storage required are both proportional to N, so very large -C datasets may be accomodated. The data points do not have -C to be equidistant in the independant variable X or uniformly -C weighted in the dependant variable Y. However, the data -C points in X must be strictly increasing. Multiple dataset -C processing (K.gt.1) is numerically more efficient dan -C separate processing of the individual datasets (K.eq.1). -C -C (3) If |MD|=3 (a priori known noise variance), any value of -C N.ge.2*M is acceptable. However, it is advisable for N-2*M -C be rather large (at least 20) if |MD|=2 (GCV). -C -C (4) For |MD| > 1, GCVSPL tries to iteratively minimize the -C selected criterion function. This minimum is unique for |MD| -C = 4, but not necessarily for |MD| = 2 or 3. Consequently, -C local optima rather that the global optimum might be found, -C and some actual findings suggest that local optima might -C yield more meaningful results than the global optimum if N -C is small. Therefore, the user has some control over the -C search procedure. If MD > 1, the iterative search starts -C from a value which yields a number of degrees of freedom -C which is approximately equal to N/2, until the first (local) -C minimum is found via a golden section search procedure -C (Utreras, 1980). If MD < -1, the value for p contained in -C WK(4) is used instead. Thus, if MD = 2 or 3 yield too noisy -C an estimate, the user might try |MD| = 1 or 4, for suitably -C selected values for p or for the number of degrees of -C freedom, and then run GCVSPL with MD = -2 or -3. The con- -C tents of N, M, K, X, WX, WY, and WK are assumed unchanged -C if MD < 0. -C -C (5) GCVSPL calculates the spline coefficient array C(N,K); -C this array can be used to calculate the spline function -C value and any of its derivatives up to the degree 2*M-1 -C at any argument T within the knot range, using subrou- -C tines SPLDER and SEARCH, and the knot array X(N). Since -C the splines are constrained at their Mth derivative, only -C the lower spline derivatives will tend to be reliable -C estimates of the underlying, true signal derivatives. -C -C (6) GCVSPL combines elements of subroutine CRVO5 by Utre- -C ras (1980), subroutine SMOOTH by Lyche et al. (1983), and -C subroutine CUBGCV by Hutchinson (1985). The trace of the -C influence matrix is assessed in a similar way as described -C by Hutchinson & de Hoog (1985). The major difference is -C that the present approach utilizes non-symmetrical B-spline -C design matrices as described by Lyche et al. (1983); there- -C fore, the original algorithm by Erisman & Tinney (1975) has -C been used, rather than the symmetrical version adopted by -C Hutchinson & de Hoog. -C -C References: -C ********** -C -C P. Craven & G. Wahba (1979), Smoothing noisy data with -C spline functions. Numerische Mathematik 31, 377-403. -C -C A.M. Erisman & W.F. Tinney (1975), On computing certain -C elements of the inverse of a sparse matrix. Communications -C of the ACM 18(3), 177-179. -C -C M.F. Hutchinson & F.R. de Hoog (1985), Smoothing noisy data -C with spline functions. Numerische Mathematik 47(1), 99-106. -C -C M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO Division of -C Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601, -C Australia. -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori (1983), Fortran -C subroutines for computing smoothing and interpolating natural -C splines. Advances in Engineering Software 5(1), 2-5. -C -C F. Utreras (1980), Un paquete de programas para ajustar curvas -C mediante funciones spline. Informe Tecnico MA-80-B-209, Depar- -C tamento de Matematicas, Faculdad de Ciencias Fisicas y Matema- -C ticas, Universidad de Chile, Santiago. -C -C Wahba, G. (1980). Numerical and statistical methods for mildly, -C moderately and severely ill-posed problems with noisy data. -C Technical report nr. 595 (February 1980). Department of Statis- -C tics, University of Madison (WI), U.S.A. -C -C Subprograms required: -C ******************** -C -C BASIS, PREP, SPLC, BANDET, BANSOL, TRINV -C -C*********************************************************************** -C - SUBROUTINE GCVSPL ( X, Y, NY, WX, WY, M, N, K, MD, VAL, C, NC, - 1 WK, IER ) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( RATIO=2D0, TAU=1.618033983D0, IBWE=7, - 1 ZERO=0D0, HALF=5D-1 , ONE=1D0, TOL=1D-6, - 2 EPS=1D-15, EPSINV=ONE/EPS ) - DIMENSION X(N), Y(NY,K), WX(N), WY(K), C(NC,K), WK(N+6*(N*M+1)) - SAVE M2, NM1, EL - DATA M2, NM1, EL / 2*0, 0D0 / -C -C*** Parameter check and work array initialization -C - IER = 0 -C*** Check on mode parameter - IF ((IABS(MD).GT.4) .OR.( MD.EQ. 0 ) .OR. - 1 ((IABS(MD).EQ.1).AND.( VAL.LT.ZERO)).OR. - 2 ((IABS(MD).EQ.3).AND.( VAL.LT.ZERO)).OR. - 3 ((IABS(MD).EQ.4).AND.((VAL.LT.ZERO) .OR.(VAL.GT.N-M)))) THEN - IER = 3 - RETURN - ENDIF -C*** Check on M and N - IF (MD.GT.0) THEN - M2 = 2 * M - NM1 = N - 1 - ELSE - IF ((M2.NE.2*M).OR.(NM1.NE.N-1)) THEN - IER = 3 - RETURN - ENDIF - ENDIF - IF ((M.LE.0).OR.(N.LT.M2)) THEN - IER = 1 - RETURN - ENDIF -C*** Check on knot sequence and weights - IF (WX(1).LE.ZERO) IER = 2 - DO 10 I=2,N - IF ((WX(I).LE.ZERO).OR.(X(I-1).GE.X(I))) IER = 2 - IF (IER.NE.0) RETURN - 10 CONTINUE - DO 15 J=1,K - IF (WY(J).LE.ZERO) IER = 2 - IF (IER.NE.0) RETURN - 15 CONTINUE -C -C*** Work array parameters (address information for covariance -C*** propagation by means of the matrices STAT, B, and WE). NB: -C*** BWE cannot be used since it is modified by function TRINV. -C - NM2P1 = N*(M2+1) - NM2M1 = N*(M2-1) -C ISTAT = 1 -C IBWE = ISTAT + 6 - IB = IBWE + NM2P1 - IWE = IB + NM2M1 -C IWK = IWE + NM2P1 -C -C*** Compute the design matrices B and WE, the ratio -C*** of their L1-norms, and check for iterative mode. -C - IF (MD.GT.0) THEN - CALL BASIS ( M, N, X, WK(IB), R1, WK(IBWE) ) - CALL PREP ( M, N, X, WX, WK(IWE), EL ) - EL = EL / R1 - ENDIF - IF (IABS(MD).NE.1) GO TO 20 -C*** Prior given value for p - R1 = VAL - GO TO 100 -C -C*** Iterate to minimize the GCV function (|MD|=2), -C*** the MSE function (|MD|=3), or to obtain the prior -C*** given number of degrees of freedom (|MD|=4). -C - 20 IF (MD.LT.-1) THEN - R1 = WK(4) - ELSE - R1 = ONE / EL - ENDIF - R2 = R1 * RATIO - GF2 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R2,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - 40 GF1 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R1,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - IF (GF1.GT.GF2) GO TO 50 - IF (WK(4).LE.ZERO) GO TO 100 - R2 = R1 - GF2 = GF1 - R1 = R1 / RATIO - GO TO 40 - 50 R3 = R2 * RATIO - 60 GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - IF (GF3.GT.GF2) GO TO 70 - IF (WK(4).GE.EPSINV) GO TO 100 - R2 = R3 - GF2 = GF3 - R3 = R3 * RATIO - GO TO 60 - 70 R2 = R3 - GF2 = GF3 - ALPHA = (R2-R1) / TAU - R4 = R1 + ALPHA - R3 = R2 - ALPHA - GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - GF4 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R4,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - 80 IF (GF3.LE.GF4) THEN - R2 = R4 - GF2 = GF4 - ERR = (R2-R1) / (R1+R2) - IF ((ERR*ERR+ONE.EQ.ONE).OR.(ERR.LE.TOL)) GO TO 90 - R4 = R3 - GF4 = GF3 - ALPHA = ALPHA / TAU - R3 = R2 - ALPHA - GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - ELSE - R1 = R3 - GF1 = GF3 - ERR = (R2-R1) / (R1+R2) - IF ((ERR*ERR+ONE.EQ.ONE).OR.(ERR.LE.TOL)) GO TO 90 - R3 = R4 - GF3 = GF4 - ALPHA = ALPHA / TAU - R4 = R1 + ALPHA - GF4 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R4,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) - ENDIF - GO TO 80 - 90 R1 = HALF * (R1+R2) -C -C*** Calculate final spline coefficients -C - 100 GF1 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R1,EPS,C,NC, - 1 WK,WK(IB),WK(IWE),EL,WK(IBWE)) -C -C*** Ready -C - RETURN - END -C BASIS.FOR, 1985-06-03 -C -C*********************************************************************** -C -C SUBROUTINE BASIS (REAL*8) -C -C Purpose: -C ******* -C -C Subroutine to assess a B-spline tableau, stored in vectorized -C form. -C -C Calling convention: -C ****************** -C -C CALL BASIS ( M, N, X, B, BL, Q ) -C -C Meaning of parameters: -C ********************* -C -C M ( I ) Half order of the spline (degree 2*M-1), -C M > 0. -C N ( I ) Number of knots, N >= 2*M. -C X(N) ( I ) Knot sequence, X(I-1) < X(I), I=2,N. -C B(1-M:M-1,N) ( O ) Output tableau. Element B(J,I) of array -C B corresponds with element b(i,i+j) of -C the tableau matrix B. -C BL ( O ) L1-norm of B. -C Q(1-M:M) ( W ) Internal work array. -C -C Remark: -C ****** -C -C This subroutine is an adaptation of subroutine BASIS from the -C paper by Lyche et al. (1983). No checking is performed on the -C validity of M and N. If the knot sequence is not strictly in- -C creasing, division by zero may occur. -C -C Reference: -C ********* -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - SUBROUTINE BASIS ( M, N, X, B, BL, Q ) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( ZERO=0D0, ONE=1D0 ) - DIMENSION X(N), B(1-M:M-1,N), Q(1-M:M) -C - IF (M.EQ.1) THEN -C*** Linear spline - DO 3 I=1,N - B(0,I) = ONE - 3 CONTINUE - BL = ONE - RETURN - ENDIF -C -C*** General splines -C - MM1 = M - 1 - MP1 = M + 1 - M2 = 2 * M - DO 15 L=1,N -C*** 1st row - DO 5 J=-MM1,M - Q(J) = ZERO - 5 CONTINUE - Q(MM1) = ONE - IF ((L.NE.1).AND.(L.NE.N)) - 1 Q(MM1) = ONE / ( X(L+1) - X(L-1) ) -C*** Successive rows - ARG = X(L) - DO 13 I=3,M2 - IR = MP1 - I - V = Q(IR) - IF (L.LT.I) THEN -C*** Left-hand B-splines - DO 6 J=L+1,I - U = V - V = Q(IR+1) - Q(IR) = U + (X(J)-ARG)*V - IR = IR + 1 - 6 CONTINUE - ENDIF - J1 = MAX0(L-I+1,1) - J2 = MIN0(L-1,N-I) - IF (J1.LE.J2) THEN -C*** Ordinary B-splines - IF (I.LT.M2) THEN - DO 8 J=J1,J2 - Y = X(I+J) - U = V - V = Q(IR+1) - Q(IR) = U + (V-U)*(Y-ARG)/(Y-X(J)) - IR = IR + 1 - 8 CONTINUE - ELSE - DO 10 J=J1,J2 - U = V - V = Q(IR+1) - Q(IR) = (ARG-X(J))*U + (X(I+J)-ARG)*V - IR = IR + 1 - 10 CONTINUE - ENDIF - ENDIF - NMIP1 = N - I + 1 - IF (NMIP1.LT.L) THEN -C*** Right-hand B-splines - DO 12 J=NMIP1,L-1 - U = V - V = Q(IR+1) - Q(IR) = (ARG-X(J))*U + V - IR = IR + 1 - 12 CONTINUE - ENDIF - 13 CONTINUE - DO 14 J=-MM1,MM1 - B(J,L) = Q(J) - 14 CONTINUE - 15 CONTINUE -C -C*** Zero unused parts of B -C - DO 17 I=1,MM1 - DO 16 K=I,MM1 - B(-K, I) = ZERO - B( K,N+1-I) = ZERO - 16 CONTINUE - 17 CONTINUE -C -C*** Assess L1-norm of B -C - BL = 0D0 - DO 19 I=1,N - DO 18 K=-MM1,MM1 - BL = BL + ABS(B(K,I)) - 18 CONTINUE - 19 CONTINUE - BL = BL / N -C -C*** Ready -C - RETURN - END -C PREP.FOR, 1985-07-04 -C -C*********************************************************************** -C -C SUBROUTINE PREP (REAL*8) -C -C Purpose: -C ******* -C -C To compute the matrix WE of weighted divided difference coeffi- -C cients needed to set up a linear system of equations for sol- -C ving B-spline smoothing problems, and its L1-norm EL. The matrix -C WE is stored in vectorized form. -C -C Calling convention: -C ****************** -C -C CALL PREP ( M, N, X, W, WE, EL ) -C -C Meaning of parameters: -C ********************* -C -C M ( I ) Half order of the B-spline (degree -C 2*M-1), with M > 0. -C N ( I ) Number of knots, with N >= 2*M. -C X(N) ( I ) Strictly increasing knot array, with -C X(I-1) < X(I), I=2,N. -C W(N) ( I ) Weight matrix (diagonal), with -C W(I).gt.0.0, I=1,N. -C WE(-M:M,N) ( O ) Array containing the weighted divided -C difference terms in vectorized format. -C Element WE(J,I) of array E corresponds -C with element e(i,i+j) of the matrix -C W**-1 * E. -C EL ( O ) L1-norm of WE. -C -C Remark: -C ****** -C -C This subroutine is an adaptation of subroutine PREP from the paper -C by Lyche et al. (1983). No checking is performed on the validity -C of M and N. Division by zero may occur if the knot sequence is -C not strictly increasing. -C -C Reference: -C ********* -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - SUBROUTINE PREP ( M, N, X, W, WE, EL ) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( ZERO=0D0, ONE=1D0 ) - DIMENSION X(N), W(N), WE((2*M+1)*N) -C -C*** Calculate the factor F1 -C - M2 = 2 * M - MP1 = M + 1 - M2M1 = M2 - 1 - M2P1 = M2 + 1 - NM = N - M - F1 = -ONE - IF (M.NE.1) THEN - DO 5 I=2,M - F1 = -F1 * I - 5 CONTINUE - DO 6 I=MP1,M2M1 - F1 = F1 * I - 6 CONTINUE - END IF -C -C*** Columnwise evaluation of the unweighted design matrix E -C - I1 = 1 - I2 = M - JM = MP1 - DO 17 J=1,N - INC = M2P1 - IF (J.GT.NM) THEN - F1 = -F1 - F = F1 - ELSE - IF (J.LT.MP1) THEN - INC = 1 - F = F1 - ELSE - F = F1 * (X(J+M)-X(J-M)) - END IF - END IF - IF ( J.GT.MP1) I1 = I1 + 1 - IF (I2.LT. N) I2 = I2 + 1 - JJ = JM -C*** Loop for divided difference coefficients - FF = F - Y = X(I1) - I1P1 = I1 + 1 - DO 11 I=I1P1,I2 - FF = FF / (Y-X(I)) - 11 CONTINUE - WE(JJ) = FF - JJ = JJ + M2 - I2M1 = I2 - 1 - IF (I1P1.LE.I2M1) THEN - DO 14 L=I1P1,I2M1 - FF = F - Y = X(L) - DO 12 I=I1,L-1 - FF = FF / (Y-X(I)) - 12 CONTINUE - DO 13 I=L+1,I2 - FF = FF / (Y-X(I)) - 13 CONTINUE - WE(JJ) = FF - JJ = JJ + M2 - 14 CONTINUE - END IF - FF = F - Y = X(I2) - DO 16 I=I1,I2M1 - FF = FF / (Y-X(I)) - 16 CONTINUE - WE(JJ) = FF - JJ = JJ + M2 - JM = JM + INC - 17 CONTINUE -C -C*** Zero the upper left and lower right corners of E -C - KL = 1 - N2M = M2P1*N + 1 - DO 19 I=1,M - KU = KL + M - I - DO 18 K=KL,KU - WE( K) = ZERO - WE(N2M-K) = ZERO - 18 CONTINUE - KL = KL + M2P1 - 19 CONTINUE -C -C*** Weighted matrix WE = W**-1 * E and its L1-norm -C - 20 JJ = 0 - EL = 0D0 - DO 22 I=1,N - WI = W(I) - DO 21 J=1,M2P1 - JJ = JJ + 1 - WE(JJ) = WE(JJ) / WI - EL = EL + ABS(WE(JJ)) - 21 CONTINUE - 22 CONTINUE - EL = EL / N -C -C*** Ready -C - RETURN - END -C SPLC.FOR, 1985-12-12 -C -C Author: H.J. Woltring -C -C Organizations: University of Nijmegen, and -C Philips Medical Systems, Eindhoven -C (The Netherlands) -C -C*********************************************************************** -C -C FUNCTION SPLC (REAL*8) -C -C Purpose: -C ******* -C -C To assess the coefficients of a B-spline and various statistical -C parameters, for a given value of the regularization parameter p. -C -C Calling convention: -C ****************** -C -C FV = SPLC ( M, N, K, Y, NY, WX, WY, MODE, VAL, P, EPS, C, NC, -C 1 STAT, B, WE, EL, BWE) -C -C Meaning of parameters: -C ********************* -C -C SPLC ( O ) GCV function value if |MODE|.eq.2, -C MSE value if |MODE|.eq.3, and absolute -C difference with the prior given number of -C degrees of freedom if |MODE|.eq.4. -C M ( I ) Half order of the B-spline (degree 2*M-1), -C with M > 0. -C N ( I ) Number of observations, with N >= 2*M. -C K ( I ) Number of datasets, with K >= 1. -C Y(NY,K) ( I ) Observed measurements. -C NY ( I ) First dimension of Y(NY,K), with NY.ge.N. -C WX(N) ( I ) Weight factors, corresponding to the -C relative inverse variance of each measure- -C ment, with WX(I) > 0.0. -C WY(K) ( I ) Weight factors, corresponding to the -C relative inverse variance of each dataset, -C with WY(J) > 0.0. -C MODE ( I ) Mode switch, as described in GCVSPL. -C VAL ( I ) Prior variance if |MODE|.eq.3, and -C prior number of degrees of freedom if -C |MODE|.eq.4. For other values of MODE, -C VAL is not used. -C P ( I ) Smoothing parameter, with P >= 0.0. If -C P.eq.0.0, an interpolating spline is -C calculated. -C EPS ( I ) Relative rounding tolerance*10.0. EPS is -C the smallest positive number such that -C EPS/10.0 + 1.0 .ne. 1.0. -C C(NC,K) ( O ) Calculated spline coefficient arrays. NB: -C the dimensions of in GCVSPL and in SPLDER -C are different -C column of C(N,K) is needed, and the proper -C column C(1,J), with J=1...K, should be used -C when calling SPLDER. -C NC ( I ) First dimension of C(NC,K), with NC.ge.N. -C STAT(6) ( O ) Statistics array. See the description in -C subroutine GCVSPL. -C B (1-M:M-1,N) ( I ) B-spline tableau as evaluated by subroutine -C BASIS. -C WE( -M:M ,N) ( I ) Weighted B-spline tableau (W**-1 * E) as -C evaluated by subroutine PREP. -C EL ( I ) L1-norm of the matrix WE as evaluated by -C subroutine PREP. -C BWE(-M:M,N) ( O ) Central 2*M+1 bands of the inverted -C matrix ( B + p * W**-1 * E )**-1 -C -C Remarks: -C ******* -C -C This subroutine combines elements of subroutine SPLC0 from the -C paper by Lyche et al. (1983), and of subroutine SPFIT1 by -C Hutchinson (1985). -C -C References: -C ********** -C -C M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO division of -C Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601, -C Australia. -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - FUNCTION SPLC( M, N, K, Y, NY, WX, WY, MODE, VAL, P, EPS, C, NC, - 1 STAT, B, WE, EL, BWE) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( ZERO=0D0, ONE=1D0, TWO=2D0 ) - DIMENSION Y(NY,K), WX(N), WY(K), C(NC,K), STAT(6), - 1 B(1-M:M-1,N), WE(-M:M,N), BWE(-M:M,N) -C -C*** Check on p-value -C - DP = P - STAT(4) = P - PEL = P * EL -C*** Pseudo-interpolation if p is too small - IF (PEL.LT.EPS) THEN - DP = EPS / EL - STAT(4) = ZERO - ENDIF -C*** Pseudo least-squares polynomial if p is too large - IF (PEL*EPS.GT.ONE) THEN - DP = ONE / (EL*EPS) - STAT(4) = DP - ENDIF -C -C*** Calculate BWE = B + p * W**-1 * E -C - DO 40 I=1,N - KM = -MIN0(M,I-1) - KP = MIN0(M,N-I) - DO 30 L=KM,KP - IF (IABS(L).EQ.M) THEN - BWE(L,I) = DP * WE(L,I) - ELSE - BWE(L,I) = B(L,I) + DP * WE(L,I) - ENDIF - 30 CONTINUE - 40 CONTINUE -C -C*** Solve BWE * C = Y, and assess TRACE [ B * BWE**-1 ] -C - CALL BANDET ( BWE, M, N ) - CALL BANSOL ( BWE, Y, NY, C, NC, M, N, K ) - STAT(3) = TRINV ( WE, BWE, M, N ) * DP - TRN = STAT(3) / N -C -C*** Compute mean-squared weighted residual -C - ESN = ZERO - DO 70 J=1,K - DO 60 I=1,N - DT = -Y(I,J) - KM = -MIN0(M-1,I-1) - KP = MIN0(M-1,N-I) - DO 50 L=KM,KP - DT = DT + B(L,I)*C(I+L,J) - 50 CONTINUE - ESN = ESN + DT*DT*WX(I)*WY(J) - 60 CONTINUE - 70 CONTINUE - ESN = ESN / (N*K) -C -C*** Calculate statistics and function value -C - STAT(6) = ESN / TRN - STAT(1) = STAT(6) / TRN - STAT(2) = ESN -C STAT(3) = trace [p*B * BWE**-1] -C STAT(4) = P - IF (IABS(MODE).NE.3) THEN -C*** Unknown variance: GCV - STAT(5) = STAT(6) - ESN - IF (IABS(MODE).EQ.1) SPLC = ZERO - IF (IABS(MODE).EQ.2) SPLC = STAT(1) - IF (IABS(MODE).EQ.4) SPLC = DABS( STAT(3) - VAL ) - ELSE -C*** Known variance: estimated mean squared error - STAT(5) = ESN - VAL*(TWO*TRN - ONE) - SPLC = STAT(5) - ENDIF -C - RETURN - END -C BANDET.FOR, 1985-06-03 -C -C*********************************************************************** -C -C SUBROUTINE BANDET (REAL*8) -C -C Purpose: -C ******* -C -C This subroutine computes the LU decomposition of an N*N matrix -C E. It is assumed that E has M bands above and M bands below the -C diagonal. The decomposition is returned in E. It is assumed that -C E can be decomposed without pivoting. The matrix E is stored in -C vectorized form in the array E(-M:M,N), where element E(J,I) of -C the array E corresponds with element e(i,i+j) of the matrix E. -C -C Calling convention: -C ****************** -C -C CALL BANDET ( E, M, N ) -C -C Meaning of parameters: -C ********************* -C -C E(-M:M,N) (I/O) Matrix to be decomposed. -C M, N ( I ) Matrix dimensioning parameters, -C M >= 0, N >= 2*M. -C -C Remark: -C ****** -C -C No checking on the validity of the input data is performed. -C If (M.le.0), no action is taken. -C -C*********************************************************************** -C - SUBROUTINE BANDET ( E, M, N ) -C - IMPLICIT REAL*8 (A-H,O-Z) - DIMENSION E(-M:M,N) -C - IF (M.LE.0) RETURN - DO 40 I=1,N - DI = E(0,I) - MI = MIN0(M,I-1) - IF (MI.GE.1) THEN - DO 10 K=1,MI - DI = DI - E(-K,I)*E(K,I-K) - 10 CONTINUE - E(0,I) = DI - ENDIF - LM = MIN0(M,N-I) - IF (LM.GE.1) THEN - DO 30 L=1,LM - DL = E(-L,I+L) - KM = MIN0(M-L,I-1) - IF (KM.GE.1) THEN - DU = E(L,I) - DO 20 K=1,KM - DU = DU - E( -K, I)*E(L+K,I-K) - DL = DL - E(-L-K,L+I)*E( K,I-K) - 20 CONTINUE - E(L,I) = DU - ENDIF - E(-L,I+L) = DL / DI - 30 CONTINUE - ENDIF - 40 CONTINUE -C -C*** Ready -C - RETURN - END -C BANSOL.FOR, 1985-12-12 -C -C*********************************************************************** -C -C SUBROUTINE BANSOL (REAL*8) -C -C Purpose: -C ******* -C -C This subroutine solves systems of linear equations given an LU -C decomposition of the design matrix. Such a decomposition is pro- -C vided by subroutine BANDET, in vectorized form. It is assumed -C that the design matrix is not singular. -C -C Calling convention: -C ****************** -C -C CALL BANSOL ( E, Y, NY, C, NC, M, N, K ) -C -C Meaning of parameters: -C ********************* -C -C E(-M:M,N) ( I ) Input design matrix, in LU-decomposed, -C vectorized form. Element E(J,I) of the -C array E corresponds with element -C e(i,i+j) of the N*N design matrix E. -C Y(NY,K) ( I ) Right hand side vectors. -C C(NC,K) ( O ) Solution vectors. -C NY, NC, M, N, K ( I ) Dimensioning parameters, with M >= 0, -C N > 2*M, and K >= 1. -C -C Remark: -C ****** -C -C This subroutine is an adaptation of subroutine BANSOL from the -C paper by Lyche et al. (1983). No checking is performed on the -C validity of the input parameters and data. Division by zero may -C occur if the system is singular. -C -C Reference: -C ********* -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - SUBROUTINE BANSOL ( E, Y, NY, C, NC, M, N, K ) -C - IMPLICIT REAL*8 (A-H,O-Z) - DIMENSION E(-M:M,N), Y(NY,K), C(NC,K) -C -C*** Check on special cases: M=0, M=1, M>1 -C - NM1 = N - 1 - IF (M-1 .lt. 0) GO TO 10 - IF (M-1 .eq. 0) GO TO 40 - GO TO 80 -C -C*** M = 0: Diagonal system -C - 10 DO 30 I=1,N - DO 20 J=1,K - C(I,J) = Y(I,J) / E(0,I) - 20 CONTINUE - 30 CONTINUE - RETURN -C -C*** M = 1: Tridiagonal system -C - 40 DO 70 J=1,K - C(1,J) = Y(1,J) - DO 50 I=2,N - C(I,J) = Y(I,J) - E(-1,I)*C(I-1,J) - 50 CONTINUE - C(N,J) = C(N,J) / E(0,N) - DO 60 I=NM1,1,-1 - C(I,J) = (C(I,J) - E( 1,I)*C(I+1,J)) / E(0,I) - 60 CONTINUE - 70 CONTINUE - RETURN -C -C*** M > 1: General system -C - 80 DO 130 J=1,K - C(1,J) = Y(1,J) - DO 100 I=2,N - MI = MIN0(M,I-1) - D = Y(I,J) - DO 90 L=1,MI - D = D - E(-L,I)*C(I-L,J) - 90 CONTINUE - C(I,J) = D - 100 CONTINUE - C(N,J) = C(N,J) / E(0,N) - DO 120 I=NM1,1,-1 - MI = MIN0(M,N-I) - D = C(I,J) - DO 110 L=1,MI - D = D - E( L,I)*C(I+L,J) - 110 CONTINUE - C(I,J) = D / E(0,I) - 120 CONTINUE - 130 CONTINUE - RETURN -C - END -C TRINV.FOR, 1985-06-03 -C -C*********************************************************************** -C -C FUNCTION TRINV (REAL*8) -C -C Purpose: -C ******* -C -C To calculate TRACE [ B * E**-1 ], where B and E are N * N -C matrices with bandwidth 2*M+1, and where E is a regular matrix -C in LU-decomposed form. B and E are stored in vectorized form, -C compatible with subroutines BANDET and BANSOL. -C -C Calling convention: -C ****************** -C -C TRACE = TRINV ( B, E, M, N ) -C -C Meaning of parameters: -C ********************* -C -C B(-M:M,N) ( I ) Input array for matrix B. Element B(J,I) -C corresponds with element b(i,i+j) of the -C matrix B. -C E(-M:M,N) (I/O) Input array for matrix E. Element E(J,I) -C corresponds with element e(i,i+j) of the -C matrix E. This matrix is stored in LU- -C decomposed form, with L unit lower tri- -C angular, and U upper triangular. The unit -C diagonal of L is not stored. Upon return, -C the array E holds the central 2*M+1 bands -C of the inverse E**-1, in similar ordering. -C M, N ( I ) Array and matrix dimensioning parameters -C (M.gt.0, N.ge.2*M+1). -C TRINV ( O ) Output function value TRACE [ B * E**-1 ] -C -C Reference: -C ********* -C -C A.M. Erisman & W.F. Tinney, On computing certain elements of the -C inverse of a sparse matrix. Communications of the ACM 18(1975), -C nr. 3, pp. 177-179. -C -C*********************************************************************** -C - REAL*8 FUNCTION TRINV ( B, E, M, N ) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( ZERO=0D0, ONE=1D0 ) - DIMENSION B(-M:M,N), E(-M:M,N) -C -C*** Assess central 2*M+1 bands of E**-1 and store in array E -C - E(0,N) = ONE / E(0,N) - DO 40 I=N-1,1,-1 - MI = MIN0(M,N-I) - DD = ONE / E(0,I) -C*** Save Ith column of L and Ith row of U, and normalize U row - DO 10 K=1,MI - E( K,N) = E( K, I) * DD - E(-K,1) = E(-K,K+I) - 10 CONTINUE - DD = DD + DD -C*** Invert around Ith pivot - DO 30 J=MI,1,-1 - DU = ZERO - DL = ZERO - DO 20 K=1,MI - DU = DU - E( K,N)*E(J-K,I+K) - DL = DL - E(-K,1)*E(K-J,I+J) - 20 CONTINUE - E( J, I) = DU - E(-J,J+I) = DL - DD = DD - (E(J,N)*DL + E(-J,1)*DU) - 30 CONTINUE - E(0,I) = 5D-1 * DD - 40 CONTINUE -C -C*** Assess TRACE [ B * E**-1 ] and clear working storage -C - DD = ZERO - DO 60 I=1,N - MN = -MIN0(M,I-1) - MP = MIN0(M,N-I) - DO 50 K=MN,MP - DD = DD + B(K,I)*E(-K,K+I) - 50 CONTINUE - 60 CONTINUE - TRINV = DD - DO 70 K=1,M - E( K,N) = ZERO - E(-K,1) = ZERO - 70 CONTINUE -C -C*** Ready -C - RETURN - END -C SPLDER.FOR, 1985-06-11 -C -C*********************************************************************** -C -C FUNCTION SPLDER (REAL*8) -C -C Purpose: -C ******* -C -C To produce the value of the function (IDER.eq.0) or of the -C IDERth derivative (IDER.gt.0) of a 2M-th order B-spline at -C the point T. The spline is described in terms of the half -C order M, the knot sequence X(N), N.ge.2*M, and the spline -C coefficients C(N). -C -C Calling convention: -C ****************** -C -C SVIDER = SPLDER ( IDER, M, N, T, X, C, L, Q ) -C -C Meaning of parameters: -C ********************* -C -C SPLDER ( O ) Function or derivative value. -C IDER ( I ) Derivative order required, with 0.le.IDER -C and IDER.le.2*M. If IDER.eq.0, the function -C value is returned; otherwise, the IDER-th -C derivative of the spline is returned. -C M ( I ) Half order of the spline, with M.gt.0. -C N ( I ) Number of knots and spline coefficients, -C with N.ge.2*M. -C T ( I ) Argument at which the spline or its deri- -C vative is to be evaluated, with X(1).le.T -C and T.le.X(N). -C X(N) ( I ) Strictly increasing knot sequence array, -C X(I-1).lt.X(I), I=2,...,N. -C C(N) ( I ) Spline coefficients, as evaluated by -C subroutine GVCSPL. -C L (I/O) L contains an integer such that: -C X(L).le.T and T.lt.X(L+1) if T is within -C the range X(1).le.T and T.lt.X(N). If -C T.lt.X(1), L is set to 0, and if T.ge.X(N), -C L is set to N. The search for L is facili- -C tated if L has approximately the right -C value on entry. -C Q(2*M) ( W ) Internal work array. -C -C Remark: -C ****** -C -C This subroutine is an adaptation of subroutine SPLDER of -C the paper by Lyche et al. (1983). No checking is performed -C on the validity of the input parameters. -C -C Reference: -C ********* -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - REAL*8 FUNCTION SPLDER ( IDER, M, N, T, X, C, L, Q ) -C - IMPLICIT REAL*8 (A-H,O-Z) - PARAMETER ( ZERO=0D0, ONE=1D0 ) - DIMENSION X(N), C(N), Q(2*M) -C -C*** Derivatives of IDER.ge.2*M are alway zero -C - M2 = 2 * M - K = M2 - IDER - IF (K.LT.1) THEN - SPLDER = ZERO - RETURN - ENDIF -C -C*** Search for the interval value L -C - CALL SEARCH ( N, X, T, L ) -C -C*** Initialize parameters and the 1st row of the B-spline -C*** coefficients tableau -C - TT = T - MP1 = M + 1 - NPM = N + M - M2M1 = M2 - 1 - K1 = K - 1 - NK = N - K - LK = L - K - LK1 = LK + 1 - LM = L - M - JL = L + 1 - JU = L + M2 - II = N - M2 - ML = -L - DO 2 J=JL,JU - IF ((J.GE.MP1).AND.(J.LE.NPM)) THEN - Q(J+ML) = C(J-M) - ELSE - Q(J+ML) = ZERO - ENDIF - 2 CONTINUE -C -C*** The following loop computes differences of the B-spline -C*** coefficients. If the value of the spline is required, -C*** differencing is not necessary. -C - IF (IDER.GT.0) THEN - JL = JL - M2 - ML = ML + M2 - DO 6 I=1,IDER - JL = JL + 1 - II = II + 1 - J1 = MAX0(1,JL) - J2 = MIN0(L,II) - MI = M2 - I - J = J2 + 1 - IF (J1.LE.J2) THEN - DO 3 JIN=J1,J2 - J = J - 1 - JM = ML + J - Q(JM) = (Q(JM) - Q(JM-1)) / (X(J+MI) - X(J)) - 3 CONTINUE - ENDIF - IF (JL.GE.1) GO TO 6 - I1 = I + 1 - J = ML + 1 - IF (I1.LE.ML) THEN - DO 5 JIN=I1,ML - J = J - 1 - Q(J) = -Q(J-1) - 5 CONTINUE - ENDIF - 6 CONTINUE - DO 7 J=1,K - Q(J) = Q(J+IDER) - 7 CONTINUE - ENDIF -C -C*** Compute lower half of the evaluation tableau -C - IF (K1.GE.1) THEN - DO 14 I=1,K1 - NKI = NK + I - IR = K - JJ = L - KI = K - I - NKI1 = NKI + 1 -C*** Right-hand B-splines - IF (L.GE.NKI1) THEN - DO 9 J=NKI1,L - Q(IR) = Q(IR-1) + (TT-X(JJ))*Q(IR) - JJ = JJ - 1 - IR = IR - 1 - 9 CONTINUE - ENDIF -C*** Middle B-splines - LK1I = LK1 + I - J1 = MAX0(1,LK1I) - J2 = MIN0(L, NKI) - IF (J1.LE.J2) THEN - DO 11 J=J1,J2 - XJKI = X(JJ+KI) - Z = Q(IR) - Q(IR) = Z + (XJKI-TT)*(Q(IR-1)-Z)/(XJKI-X(JJ)) - IR = IR - 1 - JJ = JJ - 1 - 11 CONTINUE - ENDIF -C*** Left-hand B-splines - IF (LK1I.LE.0) THEN - JJ = KI - LK1I1 = 1 - LK1I - DO 13 J=1,LK1I1 - Q(IR) = Q(IR) + (X(JJ)-TT)*Q(IR-1) - JJ = JJ - 1 - IR = IR - 1 - 13 CONTINUE - ENDIF - 14 CONTINUE - ENDIF -C -C*** Compute the return value -C - Z = Q(K) -C*** Multiply with factorial if IDER.gt.0 - IF (IDER.GT.0) THEN - DO 16 J=K,M2M1 - Z = Z * J - 16 CONTINUE - ENDIF - SPLDER = Z -C -C*** Ready -C - RETURN - END -C SEARCH.FOR, 1985-06-03 -C -C*********************************************************************** -C -C SUBROUTINE SEARCH (REAL*8) -C -C Purpose: -C ******* -C -C Given a strictly increasing knot sequence X(1) < ... < X(N), -C where N >= 1, and a real number T, this subroutine finds the -C value L such that X(L) <= T < X(L+1). If T < X(1), L = 0; -C if X(N) <= T, L = N. -C -C Calling convention: -C ****************** -C -C CALL SEARCH ( N, X, T, L ) -C -C Meaning of parameters: -C ********************* -C -C N ( I ) Knot array dimensioning parameter. -C X(N) ( I ) Stricly increasing knot array. -C T ( I ) Input argument whose knot interval is to -C be found. -C L (I/O) Knot interval parameter. The search procedure -C is facilitated if L has approximately the -C right value on entry. -C -C Remark: -C ****** -C -C This subroutine is an adaptation of subroutine SEARCH from -C the paper by Lyche et al. (1983). No checking is performed -C on the input parameters and data; the algorithm may fail if -C the input sequence is not strictly increasing. -C -C Reference: -C ********* -C -C T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines -C for computing smoothing and interpolating natural splines. -C Advances in Engineering Software 5(1983)1, pp. 2-5. -C -C*********************************************************************** -C - SUBROUTINE SEARCH ( N, X, T, L ) -C - IMPLICIT REAL*8 (A-H,O-Z) - DIMENSION X(N) -C - IF (T.LT.X(1)) THEN -C*** Out of range to the left - L = 0 - RETURN - ENDIF - IF (T.GE.X(N)) THEN -C*** Out of range to the right - L = N - RETURN - ENDIF -C*** Validate input value of L - L = MAX0(L,1) - IF (L.GE.N) L = N-1 -C -C*** Often L will be in an interval adjoining the interval found -C*** in a previous call to search -C - IF (T.GE.X(L)) GO TO 5 - L = L - 1 - IF (T.GE.X(L)) RETURN -C -C*** Perform bisection -C - IL = 1 - 3 IU = L - 4 L = (IL+IU) / 2 - IF (IU-IL.LE.1) RETURN - IF (T.LT.X(L)) GO TO 3 - IL = L - GO TO 4 - 5 IF (T.LT.X(L+1)) RETURN - L = L + 1 - IF (T.LT.X(L+1)) RETURN - IL = L + 1 - IU = N - GO TO 4 -C - END