# HG changeset patch
# User Jaroslav Hajek
# Date 1222628975 -7200
# Node ID 614c795265e28c559e6a0cd0b54178a43825f24d
# Parent 7376fd3b25cc09dd7bc1ddeb2665115f05f00e63
initial implementation of fsolve
* * *
remove old fsolve code
diff --git a/configure.in b/configure.in
--- a/configure.in
+++ b/configure.in
@@ -1999,7 +1999,7 @@
libcruft/Makerules libcruft/amos/Makefile libcruft/blas/Makefile
libcruft/daspk/Makefile libcruft/dasrt/Makefile
libcruft/dassl/Makefile libcruft/fftpack/Makefile
- libcruft/lapack/Makefile libcruft/minpack/Makefile
+ libcruft/lapack/Makefile
libcruft/misc/Makefile libcruft/odepack/Makefile
libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile
libcruft/qrupdate/Makefile
diff --git a/libcruft/Makefile.in b/libcruft/Makefile.in
--- a/libcruft/Makefile.in
+++ b/libcruft/Makefile.in
@@ -42,7 +42,7 @@
# dirname is substituted by configure and may be the empty string.
CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \
- @FFT_DIR@ @LAPACK_DIR@ lapack-xtra minpack \
+ @FFT_DIR@ @LAPACK_DIR@ lapack-xtra \
misc odepack ordered-qz quadpack qrupdate ranlib \
slatec-err slatec-fn villad
diff --git a/libcruft/minpack/Makefile.in b/libcruft/minpack/Makefile.in
deleted file mode 100644
--- a/libcruft/minpack/Makefile.in
+++ /dev/null
@@ -1,34 +0,0 @@
-# Makefile for octave's libcruft/minpack directory
-#
-# Copyright (C) 1993, 1994, 1995, 2007 John W. Eaton
-#
-# This file is part of Octave.
-#
-# Octave 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 3 of the License, or (at
-# your option) any later version.
-#
-# Octave 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 Octave; see the file COPYING. If not, see
-# .
-
-TOPDIR = ../..
-
-srcdir = @srcdir@
-top_srcdir = @top_srcdir@
-VPATH = @srcdir@
-
-EXTERNAL_DISTFILES = $(DISTFILES)
-
-FSRC = dogleg.f dpmpar.f enorm.f fdjac1.f hybrd.f hybrd1.f \
- hybrj.f hybrj1.f qform.f qrfac.f r1mpyq.f r1updt.f
-
-include $(TOPDIR)/Makeconf
-
-include ../Makerules
diff --git a/libcruft/minpack/dogleg.f b/libcruft/minpack/dogleg.f
deleted file mode 100644
--- a/libcruft/minpack/dogleg.f
+++ /dev/null
@@ -1,177 +0,0 @@
- SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
- INTEGER N,LR
- DOUBLE PRECISION DELTA
- DOUBLE PRECISION R(LR),DIAG(N),QTB(N),X(N),WA1(N),WA2(N)
-C **********
-C
-C SUBROUTINE DOGLEG
-C
-C GIVEN AN M BY N MATRIX A, AN N BY N NONSINGULAR DIAGONAL
-C MATRIX D, AN M-VECTOR B, AND A POSITIVE NUMBER DELTA, THE
-C PROBLEM IS TO DETERMINE THE CONVEX COMBINATION X OF THE
-C GAUSS-NEWTON AND SCALED GRADIENT DIRECTIONS THAT MINIMIZES
-C (A*X - B) IN THE LEAST SQUARES SENSE, SUBJECT TO THE
-C RESTRICTION THAT THE EUCLIDEAN NORM OF D*X BE AT MOST DELTA.
-C
-C THIS SUBROUTINE COMPLETES THE SOLUTION OF THE PROBLEM
-C IF IT IS PROVIDED WITH THE NECESSARY INFORMATION FROM THE
-C QR FACTORIZATION OF A. THAT IS, IF A = Q*R, WHERE Q HAS
-C ORTHOGONAL COLUMNS AND R IS AN UPPER TRIANGULAR MATRIX,
-C THEN DOGLEG EXPECTS THE FULL UPPER TRIANGLE OF R AND
-C THE FIRST N COMPONENTS OF (Q TRANSPOSE)*B.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
-C
-C WHERE
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R.
-C
-C R IS AN INPUT ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER
-C TRIANGULAR MATRIX R STORED BY ROWS.
-C
-C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(N+1))/2.
-C
-C DIAG IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE
-C DIAGONAL ELEMENTS OF THE MATRIX D.
-C
-C QTB IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST
-C N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B.
-C
-C DELTA IS A POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER
-C BOUND ON THE EUCLIDEAN NORM OF D*X.
-C
-C X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED
-C CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION AND THE
-C SCALED GRADIENT DIRECTION.
-C
-C WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N.
-C
-C SUBPROGRAMS CALLED
-C
-C MINPACK-SUPPLIED ... DPMPAR,ENORM
-C
-C FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,DSQRT
-C
-C MINPACK. VERSION OF JULY 1978.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,J,JJ,JP1,K,L
- DOUBLE PRECISION ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM,
- * TEMP,ZERO
- DOUBLE PRECISION DPMPAR,ENORM
- DATA ONE,ZERO /1.0D0,0.0D0/
-C
-C EPSMCH IS THE MACHINE PRECISION.
-C
- EPSMCH = DPMPAR(1)
-C
-C FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.
-C
- JJ = (N*(N + 1))/2 + 1
- DO 50 K = 1, N
- J = N - K + 1
- JP1 = J + 1
- JJ = JJ - K
- L = JJ + 1
- SUM = ZERO
- IF (N .LT. JP1) GO TO 20
- DO 10 I = JP1, N
- SUM = SUM + R(L)*X(I)
- L = L + 1
- 10 CONTINUE
- 20 CONTINUE
- TEMP = R(JJ)
- IF (TEMP .NE. ZERO) GO TO 40
- L = J
- DO 30 I = 1, J
- TEMP = DMAX1(TEMP,DABS(R(L)))
- L = L + N - I
- 30 CONTINUE
- TEMP = EPSMCH*TEMP
- IF (TEMP .EQ. ZERO) TEMP = EPSMCH
- 40 CONTINUE
- X(J) = (QTB(J) - SUM)/TEMP
- 50 CONTINUE
-C
-C TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.
-C
- DO 60 J = 1, N
- WA1(J) = ZERO
- WA2(J) = DIAG(J)*X(J)
- 60 CONTINUE
- QNORM = ENORM(N,WA2)
- IF (QNORM .LE. DELTA) GO TO 140
-C
-C THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.
-C NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.
-C
- L = 1
- DO 80 J = 1, N
- TEMP = QTB(J)
- DO 70 I = J, N
- WA1(I) = WA1(I) + R(L)*TEMP
- L = L + 1
- 70 CONTINUE
- WA1(J) = WA1(J)/DIAG(J)
- 80 CONTINUE
-C
-C CALCULATE THE NORM OF THE SCALED GRADIENT DIRECTION,
-C NORMALIZE, AND RESCALE THE GRADIENT.
-C
- GNORM = ENORM(N,WA1)
- SGNORM = ZERO
- ALPHA = DELTA/QNORM
- IF (GNORM .EQ. ZERO) GO TO 120
- DO 90 J = 1, N
- WA1(J) = (WA1(J)/GNORM)/DIAG(J)
- 90 CONTINUE
-C
-C CALCULATE THE POINT ALONG THE SCALED GRADIENT
-C AT WHICH THE QUADRATIC IS MINIMIZED.
-C
- L = 1
- DO 110 J = 1, N
- SUM = ZERO
- DO 100 I = J, N
- SUM = SUM + R(L)*WA1(I)
- L = L + 1
- 100 CONTINUE
- WA2(J) = SUM
- 110 CONTINUE
- TEMP = ENORM(N,WA2)
- SGNORM = (GNORM/TEMP)/TEMP
-C
-C TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.
-C
- ALPHA = ZERO
- IF (SGNORM .GE. DELTA) GO TO 120
-C
-C THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.
-C FINALLY, CALCULATE THE POINT ALONG THE DOGLEG
-C AT WHICH THE QUADRATIC IS MINIMIZED.
-C
- BNORM = ENORM(N,QTB)
- TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA)
- TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2
- * + DSQRT((TEMP-(DELTA/QNORM))**2
- * +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2))
- ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP
- 120 CONTINUE
-C
-C FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON
-C DIRECTION AND THE SCALED GRADIENT DIRECTION.
-C
- TEMP = (ONE - ALPHA)*DMIN1(SGNORM,DELTA)
- DO 130 J = 1, N
- X(J) = TEMP*WA1(J) + ALPHA*X(J)
- 130 CONTINUE
- 140 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE DOGLEG.
-C
- END
diff --git a/libcruft/minpack/dpmpar.f b/libcruft/minpack/dpmpar.f
deleted file mode 100644
--- a/libcruft/minpack/dpmpar.f
+++ /dev/null
@@ -1,52 +0,0 @@
- DOUBLE PRECISION FUNCTION DPMPAR(I)
- INTEGER I
-C **********
-C
-C FUNCTION DPMPAR
-C
-C THIS FUNCTION PROVIDES DOUBLE PRECISION MACHINE PARAMETERS
-C WHEN THE APPROPRIATE SET OF DATA STATEMENTS IS ACTIVATED (BY
-C REMOVING THE C FROM COLUMN 1) AND ALL OTHER DATA STATEMENTS ARE
-C RENDERED INACTIVE. MOST OF THE PARAMETER VALUES WERE OBTAINED
-C FROM THE CORRESPONDING BELL LABORATORIES PORT LIBRARY FUNCTION.
-C
-C THE FUNCTION STATEMENT IS
-C
-C DOUBLE PRECISION FUNCTION DPMPAR(I)
-C
-C WHERE
-C
-C I IS AN INTEGER INPUT VARIABLE SET TO 1, 2, OR 3 WHICH
-C SELECTS THE DESIRED MACHINE PARAMETER. IF THE MACHINE HAS
-C T BASE B DIGITS AND ITS SMALLEST AND LARGEST EXPONENTS ARE
-C EMIN AND EMAX, RESPECTIVELY, THEN THESE PARAMETERS ARE
-C
-C DPMPAR(1) = B**(1 - T), THE MACHINE PRECISION,
-C
-C DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
-C
-C DPMPAR(3) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
-C
-C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C Modified Mon Aug 28 14:46:17 CDT 1989 by John W. Eaton
-C (address@hidden) to use D1MACH
-C
-C **********
-C
- DOUBLE PRECISION D1MACH
-C
- IF ( I .EQ. 1 ) THEN
- DPMPAR = D1MACH(4)
- ELSEIF ( I . EQ. 2 ) THEN
- DPMPAR = D1MACH(1)
- ELSEIF ( I .EQ. 3 ) THEN
- DPMPAR = D1MACH(2)
- ENDIF
-C
- RETURN
-C
-C LAST CARD OF FUNCTION DPMPAR.
-C
- END
diff --git a/libcruft/minpack/enorm.f b/libcruft/minpack/enorm.f
deleted file mode 100644
--- a/libcruft/minpack/enorm.f
+++ /dev/null
@@ -1,108 +0,0 @@
- DOUBLE PRECISION FUNCTION ENORM(N,X)
- INTEGER N
- DOUBLE PRECISION X(N)
-C **********
-C
-C FUNCTION ENORM
-C
-C GIVEN AN N-VECTOR X, THIS FUNCTION CALCULATES THE
-C EUCLIDEAN NORM OF X.
-C
-C THE EUCLIDEAN NORM IS COMPUTED BY ACCUMULATING THE SUM OF
-C SQUARES IN THREE DIFFERENT SUMS. THE SUMS OF SQUARES FOR THE
-C SMALL AND LARGE COMPONENTS ARE SCALED SO THAT NO OVERFLOWS
-C OCCUR. NON-DESTRUCTIVE UNDERFLOWS ARE PERMITTED. UNDERFLOWS
-C AND OVERFLOWS DO NOT OCCUR IN THE COMPUTATION OF THE UNSCALED
-C SUM OF SQUARES FOR THE INTERMEDIATE COMPONENTS.
-C THE DEFINITIONS OF SMALL, INTERMEDIATE AND LARGE COMPONENTS
-C DEPEND ON TWO CONSTANTS, RDWARF AND RGIANT. THE MAIN
-C RESTRICTIONS ON THESE CONSTANTS ARE THAT RDWARF**2 NOT
-C UNDERFLOW AND RGIANT**2 NOT OVERFLOW. THE CONSTANTS
-C GIVEN HERE ARE SUITABLE FOR EVERY KNOWN COMPUTER.
-C
-C THE FUNCTION STATEMENT IS
-C
-C DOUBLE PRECISION FUNCTION ENORM(N,X)
-C
-C WHERE
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE.
-C
-C X IS AN INPUT ARRAY OF LENGTH N.
-C
-C SUBPROGRAMS CALLED
-C
-C FORTRAN-SUPPLIED ... DABS,DSQRT
-C
-C MINPACK. VERSION OF OCTOBER 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I
- DOUBLE PRECISION AGIANT,FLOATN,ONE,RDWARF,RGIANT,S1,S2,S3,XABS,
- * X1MAX,X3MAX,ZERO
- DATA ONE,ZERO,RDWARF,RGIANT /1.0D0,0.0D0,3.834D-20,1.304D19/
- S1 = ZERO
- S2 = ZERO
- S3 = ZERO
- X1MAX = ZERO
- X3MAX = ZERO
- FLOATN = N
- AGIANT = RGIANT/FLOATN
- DO 90 I = 1, N
- XABS = DABS(X(I))
- IF (XABS .GT. RDWARF .AND. XABS .LT. AGIANT) GO TO 70
- IF (XABS .LE. RDWARF) GO TO 30
-C
-C SUM FOR LARGE COMPONENTS.
-C
- IF (XABS .LE. X1MAX) GO TO 10
- S1 = ONE + S1*(X1MAX/XABS)**2
- X1MAX = XABS
- GO TO 20
- 10 CONTINUE
- S1 = S1 + (XABS/X1MAX)**2
- 20 CONTINUE
- GO TO 60
- 30 CONTINUE
-C
-C SUM FOR SMALL COMPONENTS.
-C
- IF (XABS .LE. X3MAX) GO TO 40
- S3 = ONE + S3*(X3MAX/XABS)**2
- X3MAX = XABS
- GO TO 50
- 40 CONTINUE
- IF (XABS .NE. ZERO) S3 = S3 + (XABS/X3MAX)**2
- 50 CONTINUE
- 60 CONTINUE
- GO TO 80
- 70 CONTINUE
-C
-C SUM FOR INTERMEDIATE COMPONENTS.
-C
- S2 = S2 + XABS**2
- 80 CONTINUE
- 90 CONTINUE
-C
-C CALCULATION OF NORM.
-C
- IF (S1 .EQ. ZERO) GO TO 100
- ENORM = X1MAX*DSQRT(S1+(S2/X1MAX)/X1MAX)
- GO TO 130
- 100 CONTINUE
- IF (S2 .EQ. ZERO) GO TO 110
- IF (S2 .GE. X3MAX)
- * ENORM = DSQRT(S2*(ONE+(X3MAX/S2)*(X3MAX*S3)))
- IF (S2 .LT. X3MAX)
- * ENORM = DSQRT(X3MAX*((S2/X3MAX)+(X3MAX*S3)))
- GO TO 120
- 110 CONTINUE
- ENORM = X3MAX*DSQRT(S3)
- 120 CONTINUE
- 130 CONTINUE
- RETURN
-C
-C LAST CARD OF FUNCTION ENORM.
-C
- END
diff --git a/libcruft/minpack/fdjac1.f b/libcruft/minpack/fdjac1.f
deleted file mode 100644
--- a/libcruft/minpack/fdjac1.f
+++ /dev/null
@@ -1,151 +0,0 @@
- SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,
- * WA1,WA2)
- INTEGER N,LDFJAC,IFLAG,ML,MU
- DOUBLE PRECISION EPSFCN
- DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N),WA1(N),WA2(N)
- EXTERNAL FCN
-C **********
-C
-C SUBROUTINE FDJAC1
-C
-C THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION
-C TO THE N BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED
-C PROBLEM OF N FUNCTIONS IN N VARIABLES. IF THE JACOBIAN HAS
-C A BANDED FORM, THEN FUNCTION EVALUATIONS ARE SAVED BY ONLY
-C APPROXIMATING THE NONZERO TERMS.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,
-C WA1,WA2)
-C
-C WHERE
-C
-C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
-C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED
-C IN AN EXTERNAL STATEMENT IN THE USER CALLING
-C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
-C
-C SUBROUTINE FCN(N,X,FVEC,IFLAG)
-C INTEGER N,IFLAG
-C DOUBLE PRECISION X(N),FVEC(N)
-C ----------
-C CALCULATE THE FUNCTIONS AT X AND
-C RETURN THIS VECTOR IN FVEC.
-C ----------
-C RETURN
-C END
-C
-C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
-C THE USER WANTS TO TERMINATE EXECUTION OF FDJAC1.
-C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF FUNCTIONS AND VARIABLES.
-C
-C X IS AN INPUT ARRAY OF LENGTH N.
-C
-C FVEC IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE
-C FUNCTIONS EVALUATED AT X.
-C
-C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
-C APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X.
-C
-C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
-C
-C IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE
-C THE EXECUTION OF FDJAC1. SEE DESCRIPTION OF FCN.
-C
-C ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES
-C THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE
-C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET
-C ML TO AT LEAST N - 1.
-C
-C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE
-C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS
-C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE
-C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS
-C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE
-C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE
-C PRECISION.
-C
-C MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES
-C THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE
-C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET
-C MU TO AT LEAST N - 1.
-C
-C WA1 AND WA2 ARE WORK ARRAYS OF LENGTH N. IF ML + MU + 1 IS AT
-C LEAST N, THEN THE JACOBIAN IS CONSIDERED DENSE, AND WA2 IS
-C NOT REFERENCED.
-C
-C SUBPROGRAMS CALLED
-C
-C MINPACK-SUPPLIED ... DPMPAR
-C
-C FORTRAN-SUPPLIED ... DABS,DMAX1,DSQRT
-C
-C MINPACK. VERSION OF JUNE 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,J,K,MSUM
- DOUBLE PRECISION EPS,EPSMCH,H,TEMP,ZERO
- DOUBLE PRECISION DPMPAR
- DATA ZERO /0.0D0/
-C
-C EPSMCH IS THE MACHINE PRECISION.
-C
- EPSMCH = DPMPAR(1)
-C
- EPS = DSQRT(DMAX1(EPSFCN,EPSMCH))
- MSUM = ML + MU + 1
- IF (MSUM .LT. N) GO TO 40
-C
-C COMPUTATION OF DENSE APPROXIMATE JACOBIAN.
-C
- DO 20 J = 1, N
- TEMP = X(J)
- H = EPS*DABS(TEMP)
- IF (H .EQ. ZERO) H = EPS
- X(J) = TEMP + H
- CALL FCN(N,X,WA1,IFLAG)
- IF (IFLAG .LT. 0) GO TO 30
- X(J) = TEMP
- DO 10 I = 1, N
- FJAC(I,J) = (WA1(I) - FVEC(I))/H
- 10 CONTINUE
- 20 CONTINUE
- 30 CONTINUE
- GO TO 110
- 40 CONTINUE
-C
-C COMPUTATION OF BANDED APPROXIMATE JACOBIAN.
-C
- DO 90 K = 1, MSUM
- DO 60 J = K, N, MSUM
- WA2(J) = X(J)
- H = EPS*DABS(WA2(J))
- IF (H .EQ. ZERO) H = EPS
- X(J) = WA2(J) + H
- 60 CONTINUE
- CALL FCN(N,X,WA1,IFLAG)
- IF (IFLAG .LT. 0) GO TO 100
- DO 80 J = K, N, MSUM
- X(J) = WA2(J)
- H = EPS*DABS(WA2(J))
- IF (H .EQ. ZERO) H = EPS
- DO 70 I = 1, N
- FJAC(I,J) = ZERO
- IF (I .GE. J - MU .AND. I .LE. J + ML)
- * FJAC(I,J) = (WA1(I) - FVEC(I))/H
- 70 CONTINUE
- 80 CONTINUE
- 90 CONTINUE
- 100 CONTINUE
- 110 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE FDJAC1.
-C
- END
diff --git a/libcruft/minpack/hybrd.f b/libcruft/minpack/hybrd.f
deleted file mode 100644
--- a/libcruft/minpack/hybrd.f
+++ /dev/null
@@ -1,459 +0,0 @@
- SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,DIAG,
- * MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC,R,LR,
- * QTF,WA1,WA2,WA3,WA4)
- INTEGER N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR
- DOUBLE PRECISION XTOL,EPSFCN,FACTOR
- DOUBLE PRECISION X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),
- * QTF(N),WA1(N),WA2(N),WA3(N),WA4(N)
- EXTERNAL FCN
-C **********
-C
-C SUBROUTINE HYBRD
-C
-C THE PURPOSE OF HYBRD IS TO FIND A ZERO OF A SYSTEM OF
-C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
-C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A
-C SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS
-C THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,
-C DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,
-C LDFJAC,R,LR,QTF,WA1,WA2,WA3,WA4)
-C
-C WHERE
-C
-C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
-C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED
-C IN AN EXTERNAL STATEMENT IN THE USER CALLING
-C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
-C
-C SUBROUTINE FCN(N,X,FVEC,IFLAG)
-C INTEGER N,IFLAG
-C DOUBLE PRECISION X(N),FVEC(N)
-C ----------
-C CALCULATE THE FUNCTIONS AT X AND
-C RETURN THIS VECTOR IN FVEC.
-C ---------
-C RETURN
-C END
-C
-C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
-C THE USER WANTS TO TERMINATE EXECUTION OF HYBRD.
-C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF FUNCTIONS AND VARIABLES.
-C
-C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
-C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
-C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
-C
-C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE FUNCTIONS EVALUATED AT THE OUTPUT X.
-C
-C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION
-C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE
-C ITERATES IS AT MOST XTOL.
-C
-C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
-C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV
-C BY THE END OF AN ITERATION.
-C
-C ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES
-C THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE
-C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET
-C ML TO AT LEAST N - 1.
-C
-C MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES
-C THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE
-C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET
-C MU TO AT LEAST N - 1.
-C
-C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE
-C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS
-C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE
-C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS
-C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE
-C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE
-C PRECISION.
-C
-C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE
-C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG
-C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS IMPLICIT
-C (MULTIPLICATIVE) SCALE FACTORS FOR THE VARIABLES.
-C
-C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE
-C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,
-C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER
-C VALUES OF MODE ARE EQUIVALENT TO MODE = 1.
-C
-C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE
-C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF
-C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE
-C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE
-C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.
-C
-C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED
-C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,
-C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST
-C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND
-C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE
-C FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS
-C OF FCN WITH IFLAG = 0 ARE MADE.
-C
-C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
-C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
-C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
-C INFO IS SET AS FOLLOWS.
-C
-C INFO = 0 IMPROPER INPUT PARAMETERS.
-C
-C INFO = 1 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES
-C IS AT MOST TOL.
-C
-C INFO = 2 NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED
-C MAXFEV.
-C
-C INFO = 3 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
-C THE APPROXIMATE SOLUTION X IS POSSIBLE.
-C
-C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS, AS
-C MEASURED BY THE IMPROVEMENT FROM THE LAST
-C FIVE JACOBIAN EVALUATIONS.
-C
-C INFO = 5 ITERATION IS NOT MAKING GOOD PROGRESS, AS
-C MEASURED BY THE IMPROVEMENT FROM THE LAST
-C TEN ITERATIONS.
-C
-C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
-C CALLS TO FCN.
-C
-C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
-C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION
-C OF THE FINAL APPROXIMATE JACOBIAN.
-C
-C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
-C
-C R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE
-C UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION
-C OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE.
-C
-C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(N+1))/2.
-C
-C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE VECTOR (Q TRANSPOSE)*FVEC.
-C
-C WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N.
-C
-C SUBPROGRAMS CALLED
-C
-C USER-SUPPLIED ...... FCN
-C
-C MINPACK-SUPPLIED ... DOGLEG,DPMPAR,ENORM,FDJAC1,
-C QFORM,QRFAC,R1MPYQ,R1UPDT
-C
-C FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,MIN0,MOD
-C
-C MINPACK. VERSION OF SEPTEMBER 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,IFLAG,ITER,J,JM1,L,MSUM,NCFAIL,NCSUC,NSLOW1,NSLOW2
- INTEGER IWA(1)
- LOGICAL JEVAL,SING
- DOUBLE PRECISION ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,
- * PRERED,P1,P5,P001,P0001,RATIO,SUM,TEMP,XNORM,
- * ZERO
- DOUBLE PRECISION DPMPAR,ENORM
- DATA ONE,P1,P5,P001,P0001,ZERO
- * /1.0D0,1.0D-1,5.0D-1,1.0D-3,1.0D-4,0.0D0/
-C
-C EPSMCH IS THE MACHINE PRECISION.
-C
- EPSMCH = DPMPAR(1)
-C
- INFO = 0
- IFLAG = 0
- NFEV = 0
-C
-C CHECK THE INPUT PARAMETERS FOR ERRORS.
-C
- IF (N .LE. 0 .OR. XTOL .LT. ZERO .OR. MAXFEV .LE. 0
- * .OR. ML .LT. 0 .OR. MU .LT. 0 .OR. FACTOR .LE. ZERO
- * .OR. LDFJAC .LT. N .OR. LR .LT. (N*(N + 1))/2) GO TO 300
- IF (MODE .NE. 2) GO TO 20
- DO 10 J = 1, N
- IF (DIAG(J) .LE. ZERO) GO TO 300
- 10 CONTINUE
- 20 CONTINUE
-C
-C EVALUATE THE FUNCTION AT THE STARTING POINT
-C AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- CALL FCN(N,X,FVEC,IFLAG)
- NFEV = 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM = ENORM(N,FVEC)
-C
-C DETERMINE THE NUMBER OF CALLS TO FCN NEEDED TO COMPUTE
-C THE JACOBIAN MATRIX.
-C
- MSUM = MIN0(ML+MU+1,N)
-C
-C INITIALIZE ITERATION COUNTER AND MONITORS.
-C
- ITER = 1
- NCSUC = 0
- NCFAIL = 0
- NSLOW1 = 0
- NSLOW2 = 0
-C
-C BEGINNING OF THE OUTER LOOP.
-C
- 30 CONTINUE
- JEVAL = .TRUE.
-C
-C CALCULATE THE JACOBIAN MATRIX.
-C
- IFLAG = 2
- CALL FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1,
- * WA2)
- NFEV = NFEV + MSUM
- IF (IFLAG .LT. 0) GO TO 300
-C
-C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
-C
- CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)
-C
-C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
-C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
-C
- IF (ITER .NE. 1) GO TO 70
- IF (MODE .EQ. 2) GO TO 50
- DO 40 J = 1, N
- DIAG(J) = WA2(J)
- IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
- 40 CONTINUE
- 50 CONTINUE
-C
-C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
-C AND INITIALIZE THE STEP BOUND DELTA.
-C
- DO 60 J = 1, N
- WA3(J) = DIAG(J)*X(J)
- 60 CONTINUE
- XNORM = ENORM(N,WA3)
- DELTA = FACTOR*XNORM
- IF (DELTA .EQ. ZERO) DELTA = FACTOR
- 70 CONTINUE
-C
-C FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.
-C
- DO 80 I = 1, N
- QTF(I) = FVEC(I)
- 80 CONTINUE
- DO 120 J = 1, N
- IF (FJAC(J,J) .EQ. ZERO) GO TO 110
- SUM = ZERO
- DO 90 I = J, N
- SUM = SUM + FJAC(I,J)*QTF(I)
- 90 CONTINUE
- TEMP = -SUM/FJAC(J,J)
- DO 100 I = J, N
- QTF(I) = QTF(I) + FJAC(I,J)*TEMP
- 100 CONTINUE
- 110 CONTINUE
- 120 CONTINUE
-C
-C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.
-C
- SING = .FALSE.
- DO 150 J = 1, N
- L = J
- JM1 = J - 1
- IF (JM1 .LT. 1) GO TO 140
- DO 130 I = 1, JM1
- R(L) = FJAC(I,J)
- L = L + N - I
- 130 CONTINUE
- 140 CONTINUE
- R(L) = WA1(J)
- IF (WA1(J) .EQ. ZERO) SING = .TRUE.
- 150 CONTINUE
-C
-C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.
-C
- CALL QFORM(N,N,FJAC,LDFJAC,WA1)
-C
-C RESCALE IF NECESSARY.
-C
- IF (MODE .EQ. 2) GO TO 170
- DO 160 J = 1, N
- DIAG(J) = DMAX1(DIAG(J),WA2(J))
- 160 CONTINUE
- 170 CONTINUE
-C
-C BEGINNING OF THE INNER LOOP.
-C
- 180 CONTINUE
-C
-C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
-C
- IF (NPRINT .LE. 0) GO TO 190
- IFLAG = 0
- IF (MOD(ITER-1,NPRINT) .EQ. 0) CALL FCN(N,X,FVEC,IFLAG)
- IF (IFLAG .LT. 0) GO TO 300
- 190 CONTINUE
-C
-C DETERMINE THE DIRECTION P.
-C
- CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)
-C
-C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
-C
- DO 200 J = 1, N
- WA1(J) = -WA1(J)
- WA2(J) = X(J) + WA1(J)
- WA3(J) = DIAG(J)*WA1(J)
- 200 CONTINUE
- PNORM = ENORM(N,WA3)
-C
-C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
-C
- IF (ITER .EQ. 1) DELTA = DMIN1(DELTA,PNORM)
-C
-C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- CALL FCN(N,WA2,WA4,IFLAG)
- NFEV = NFEV + 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM1 = ENORM(N,WA4)
-C
-C COMPUTE THE SCALED ACTUAL REDUCTION.
-C
- ACTRED = -ONE
- IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
-C
-C COMPUTE THE SCALED PREDICTED REDUCTION.
-C
- L = 1
- DO 220 I = 1, N
- SUM = ZERO
- DO 210 J = I, N
- SUM = SUM + R(L)*WA1(J)
- L = L + 1
- 210 CONTINUE
- WA3(I) = QTF(I) + SUM
- 220 CONTINUE
- TEMP = ENORM(N,WA3)
- PRERED = ZERO
- IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2
-C
-C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
-C REDUCTION.
-C
- RATIO = ZERO
- IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED
-C
-C UPDATE THE STEP BOUND.
-C
- IF (RATIO .GE. P1) GO TO 230
- NCSUC = 0
- NCFAIL = NCFAIL + 1
- DELTA = P5*DELTA
- GO TO 240
- 230 CONTINUE
- NCFAIL = 0
- NCSUC = NCSUC + 1
- IF (RATIO .GE. P5 .OR. NCSUC .GT. 1)
- * DELTA = DMAX1(DELTA,PNORM/P5)
- IF (DABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5
- 240 CONTINUE
-C
-C TEST FOR SUCCESSFUL ITERATION.
-C
- IF (RATIO .LT. P0001) GO TO 260
-C
-C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
-C
- DO 250 J = 1, N
- X(J) = WA2(J)
- WA2(J) = DIAG(J)*X(J)
- FVEC(J) = WA4(J)
- 250 CONTINUE
- XNORM = ENORM(N,WA2)
- FNORM = FNORM1
- ITER = ITER + 1
- 260 CONTINUE
-C
-C DETERMINE THE PROGRESS OF THE ITERATION.
-C
- NSLOW1 = NSLOW1 + 1
- IF (ACTRED .GE. P001) NSLOW1 = 0
- IF (JEVAL) NSLOW2 = NSLOW2 + 1
- IF (ACTRED .GE. P1) NSLOW2 = 0
-C
-C TEST FOR CONVERGENCE.
-C
- IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1
- IF (INFO .NE. 0) GO TO 300
-C
-C TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
-C
- IF (NFEV .GE. MAXFEV) INFO = 2
- IF (P1*DMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3
- IF (NSLOW2 .EQ. 5) INFO = 4
- IF (NSLOW1 .EQ. 10) INFO = 5
- IF (INFO .NE. 0) GO TO 300
-C
-C CRITERION FOR RECALCULATING JACOBIAN APPROXIMATION
-C BY FORWARD DIFFERENCES.
-C
- IF (NCFAIL .EQ. 2) GO TO 290
-C
-C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN
-C AND UPDATE QTF IF NECESSARY.
-C
- DO 280 J = 1, N
- SUM = ZERO
- DO 270 I = 1, N
- SUM = SUM + FJAC(I,J)*WA4(I)
- 270 CONTINUE
- WA2(J) = (SUM - WA3(J))/PNORM
- WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM)
- IF (RATIO .GE. P0001) QTF(J) = SUM
- 280 CONTINUE
-C
-C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.
-C
- CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)
- CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)
- CALL R1MPYQ(1,N,QTF,1,WA2,WA3)
-C
-C END OF THE INNER LOOP.
-C
- JEVAL = .FALSE.
- GO TO 180
- 290 CONTINUE
-C
-C END OF THE OUTER LOOP.
-C
- GO TO 30
- 300 CONTINUE
-C
-C TERMINATION, EITHER NORMAL OR USER IMPOSED.
-C
- IF (IFLAG .LT. 0) INFO = IFLAG
- IFLAG = 0
- IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,IFLAG)
- RETURN
-C
-C LAST CARD OF SUBROUTINE HYBRD.
-C
- END
diff --git a/libcruft/minpack/hybrd1.f b/libcruft/minpack/hybrd1.f
deleted file mode 100644
--- a/libcruft/minpack/hybrd1.f
+++ /dev/null
@@ -1,123 +0,0 @@
- SUBROUTINE HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA)
- INTEGER N,INFO,LWA
- DOUBLE PRECISION TOL
- DOUBLE PRECISION X(N),FVEC(N),WA(LWA)
- EXTERNAL FCN
-C **********
-C
-C SUBROUTINE HYBRD1
-C
-C THE PURPOSE OF HYBRD1 IS TO FIND A ZERO OF A SYSTEM OF
-C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
-C OF THE POWELL HYBRID METHOD. THIS IS DONE BY USING THE
-C MORE GENERAL NONLINEAR EQUATION SOLVER HYBRD. THE USER
-C MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE FUNCTIONS.
-C THE JACOBIAN IS THEN CALCULATED BY A FORWARD-DIFFERENCE
-C APPROXIMATION.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA)
-C
-C WHERE
-C
-C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
-C CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED
-C IN AN EXTERNAL STATEMENT IN THE USER CALLING
-C PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
-C
-C SUBROUTINE FCN(N,X,FVEC,IFLAG)
-C INTEGER N,IFLAG
-C DOUBLE PRECISION X(N),FVEC(N)
-C ----------
-C CALCULATE THE FUNCTIONS AT X AND
-C RETURN THIS VECTOR IN FVEC.
-C ---------
-C RETURN
-C END
-C
-C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
-C THE USER WANTS TO TERMINATE EXECUTION OF HYBRD1.
-C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF FUNCTIONS AND VARIABLES.
-C
-C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
-C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
-C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
-C
-C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE FUNCTIONS EVALUATED AT THE OUTPUT X.
-C
-C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
-C WHEN THE ALGORITHM ESTIMATES THAT THE RELATIVE ERROR
-C BETWEEN X AND THE SOLUTION IS AT MOST TOL.
-C
-C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
-C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
-C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
-C INFO IS SET AS FOLLOWS.
-C
-C INFO = 0 IMPROPER INPUT PARAMETERS.
-C
-C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR
-C BETWEEN X AND THE SOLUTION IS AT MOST TOL.
-C
-C INFO = 2 NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED
-C 200*(N+1).
-C
-C INFO = 3 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
-C THE APPROXIMATE SOLUTION X IS POSSIBLE.
-C
-C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS.
-C
-C WA IS A WORK ARRAY OF LENGTH LWA.
-C
-C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(3*N+13))/2.
-C
-C SUBPROGRAMS CALLED
-C
-C USER-SUPPLIED ...... FCN
-C
-C MINPACK-SUPPLIED ... HYBRD
-C
-C MINPACK. VERSION OF JULY 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER INDEX,J,LR,MAXFEV,ML,MODE,MU,NFEV,NPRINT
- DOUBLE PRECISION EPSFCN,FACTOR,ONE,XTOL,ZERO
- DATA FACTOR,ONE,ZERO /1.0D2,1.0D0,0.0D0/
- INFO = 0
-C
-C CHECK THE INPUT PARAMETERS FOR ERRORS.
-C
- IF (N .LE. 0 .OR. TOL .LT. ZERO .OR. LWA .LT. (N*(3*N + 13))/2)
- * GO TO 20
-C
-C CALL HYBRD.
-C
- MAXFEV = 200*(N + 1)
- XTOL = TOL
- ML = N - 1
- MU = N - 1
- EPSFCN = ZERO
- MODE = 2
- DO 10 J = 1, N
- WA(J) = ONE
- 10 CONTINUE
- NPRINT = 0
- LR = (N*(N + 1))/2
- INDEX = 6*N + LR
- CALL HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,WA(1),MODE,
- * FACTOR,NPRINT,INFO,NFEV,WA(INDEX+1),N,WA(6*N+1),LR,
- * WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1))
- IF (INFO .EQ. 5) INFO = 4
- 20 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE HYBRD1.
-C
- END
diff --git a/libcruft/minpack/hybrj.f b/libcruft/minpack/hybrj.f
deleted file mode 100644
--- a/libcruft/minpack/hybrj.f
+++ /dev/null
@@ -1,441 +0,0 @@
- SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,MODE,
- * FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,WA1,WA2,
- * WA3,WA4)
- INTEGER N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR
- DOUBLE PRECISION XTOL,FACTOR
- DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N),DIAG(N),R(LR),
- * QTF(N),WA1(N),WA2(N),WA3(N),WA4(N)
- EXTERNAL FCN
-C **********
-C
-C SUBROUTINE HYBRJ
-C
-C THE PURPOSE OF HYBRJ IS TO FIND A ZERO OF A SYSTEM OF
-C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
-C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A
-C SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,
-C MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,
-C WA1,WA2,WA3,WA4)
-C
-C WHERE
-C
-C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
-C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST
-C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER
-C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
-C
-C SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
-C INTEGER N,LDFJAC,IFLAG
-C DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N)
-C ----------
-C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND
-C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC.
-C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND
-C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC.
-C ---------
-C RETURN
-C END
-C
-C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
-C THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ.
-C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF FUNCTIONS AND VARIABLES.
-C
-C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
-C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
-C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
-C
-C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE FUNCTIONS EVALUATED AT THE OUTPUT X.
-C
-C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
-C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION
-C OF THE FINAL APPROXIMATE JACOBIAN.
-C
-C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
-C
-C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION
-C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE
-C ITERATES IS AT MOST XTOL.
-C
-C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
-C OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1
-C HAS REACHED MAXFEV.
-C
-C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE
-C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG
-C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS IMPLICIT
-C (MULTIPLICATIVE) SCALE FACTORS FOR THE VARIABLES.
-C
-C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE
-C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,
-C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER
-C VALUES OF MODE ARE EQUIVALENT TO MODE = 1.
-C
-C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE
-C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF
-C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE
-C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE
-C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.
-C
-C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED
-C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,
-C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST
-C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND
-C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE
-C FOR PRINTING. FVEC AND FJAC SHOULD NOT BE ALTERED.
-C IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS OF FCN
-C WITH IFLAG = 0 ARE MADE.
-C
-C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
-C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
-C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
-C INFO IS SET AS FOLLOWS.
-C
-C INFO = 0 IMPROPER INPUT PARAMETERS.
-C
-C INFO = 1 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES
-C IS AT MOST TOL.
-C
-C INFO = 2 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS
-C REACHED MAXFEV.
-C
-C INFO = 3 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
-C THE APPROXIMATE SOLUTION X IS POSSIBLE.
-C
-C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS, AS
-C MEASURED BY THE IMPROVEMENT FROM THE LAST
-C FIVE JACOBIAN EVALUATIONS.
-C
-C INFO = 5 ITERATION IS NOT MAKING GOOD PROGRESS, AS
-C MEASURED BY THE IMPROVEMENT FROM THE LAST
-C TEN ITERATIONS.
-C
-C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
-C CALLS TO FCN WITH IFLAG = 1.
-C
-C NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
-C CALLS TO FCN WITH IFLAG = 2.
-C
-C R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE
-C UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION
-C OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE.
-C
-C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(N+1))/2.
-C
-C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE VECTOR (Q TRANSPOSE)*FVEC.
-C
-C WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N.
-C
-C SUBPROGRAMS CALLED
-C
-C USER-SUPPLIED ...... FCN
-C
-C MINPACK-SUPPLIED ... DOGLEG,DPMPAR,ENORM,
-C QFORM,QRFAC,R1MPYQ,R1UPDT
-C
-C FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,MOD
-C
-C MINPACK. VERSION OF SEPTEMBER 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,IFLAG,ITER,J,JM1,L,NCFAIL,NCSUC,NSLOW1,NSLOW2
- INTEGER IWA(1)
- LOGICAL JEVAL,SING
- DOUBLE PRECISION ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,
- * PRERED,P1,P5,P001,P0001,RATIO,SUM,TEMP,XNORM,
- * ZERO
- DOUBLE PRECISION DPMPAR,ENORM
- DATA ONE,P1,P5,P001,P0001,ZERO
- * /1.0D0,1.0D-1,5.0D-1,1.0D-3,1.0D-4,0.0D0/
-C
-C EPSMCH IS THE MACHINE PRECISION.
-C
- EPSMCH = DPMPAR(1)
-C
- INFO = 0
- IFLAG = 0
- NFEV = 0
- NJEV = 0
-C
-C CHECK THE INPUT PARAMETERS FOR ERRORS.
-C
- IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. XTOL .LT. ZERO
- * .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO
- * .OR. LR .LT. (N*(N + 1))/2) GO TO 300
- IF (MODE .NE. 2) GO TO 20
- DO 10 J = 1, N
- IF (DIAG(J) .LE. ZERO) GO TO 300
- 10 CONTINUE
- 20 CONTINUE
-C
-C EVALUATE THE FUNCTION AT THE STARTING POINT
-C AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
- NFEV = 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM = ENORM(N,FVEC)
-C
-C INITIALIZE ITERATION COUNTER AND MONITORS.
-C
- ITER = 1
- NCSUC = 0
- NCFAIL = 0
- NSLOW1 = 0
- NSLOW2 = 0
-C
-C BEGINNING OF THE OUTER LOOP.
-C
- 30 CONTINUE
- JEVAL = .TRUE.
-C
-C CALCULATE THE JACOBIAN MATRIX.
-C
- IFLAG = 2
- CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
- NJEV = NJEV + 1
- IF (IFLAG .LT. 0) GO TO 300
-C
-C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
-C
- CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)
-C
-C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
-C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
-C
- IF (ITER .NE. 1) GO TO 70
- IF (MODE .EQ. 2) GO TO 50
- DO 40 J = 1, N
- DIAG(J) = WA2(J)
- IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
- 40 CONTINUE
- 50 CONTINUE
-C
-C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
-C AND INITIALIZE THE STEP BOUND DELTA.
-C
- DO 60 J = 1, N
- WA3(J) = DIAG(J)*X(J)
- 60 CONTINUE
- XNORM = ENORM(N,WA3)
- DELTA = FACTOR*XNORM
- IF (DELTA .EQ. ZERO) DELTA = FACTOR
- 70 CONTINUE
-C
-C FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.
-C
- DO 80 I = 1, N
- QTF(I) = FVEC(I)
- 80 CONTINUE
- DO 120 J = 1, N
- IF (FJAC(J,J) .EQ. ZERO) GO TO 110
- SUM = ZERO
- DO 90 I = J, N
- SUM = SUM + FJAC(I,J)*QTF(I)
- 90 CONTINUE
- TEMP = -SUM/FJAC(J,J)
- DO 100 I = J, N
- QTF(I) = QTF(I) + FJAC(I,J)*TEMP
- 100 CONTINUE
- 110 CONTINUE
- 120 CONTINUE
-C
-C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.
-C
- SING = .FALSE.
- DO 150 J = 1, N
- L = J
- JM1 = J - 1
- IF (JM1 .LT. 1) GO TO 140
- DO 130 I = 1, JM1
- R(L) = FJAC(I,J)
- L = L + N - I
- 130 CONTINUE
- 140 CONTINUE
- R(L) = WA1(J)
- IF (WA1(J) .EQ. ZERO) SING = .TRUE.
- 150 CONTINUE
-C
-C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.
-C
- CALL QFORM(N,N,FJAC,LDFJAC,WA1)
-C
-C RESCALE IF NECESSARY.
-C
- IF (MODE .EQ. 2) GO TO 170
- DO 160 J = 1, N
- DIAG(J) = DMAX1(DIAG(J),WA2(J))
- 160 CONTINUE
- 170 CONTINUE
-C
-C BEGINNING OF THE INNER LOOP.
-C
- 180 CONTINUE
-C
-C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
-C
- IF (NPRINT .LE. 0) GO TO 190
- IFLAG = 0
- IF (MOD(ITER-1,NPRINT) .EQ. 0)
- * CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
- IF (IFLAG .LT. 0) GO TO 300
- 190 CONTINUE
-C
-C DETERMINE THE DIRECTION P.
-C
- CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)
-C
-C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
-C
- DO 200 J = 1, N
- WA1(J) = -WA1(J)
- WA2(J) = X(J) + WA1(J)
- WA3(J) = DIAG(J)*WA1(J)
- 200 CONTINUE
- PNORM = ENORM(N,WA3)
-C
-C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
-C
- IF (ITER .EQ. 1) DELTA = DMIN1(DELTA,PNORM)
-C
-C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- CALL FCN(N,WA2,WA4,FJAC,LDFJAC,IFLAG)
- NFEV = NFEV + 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM1 = ENORM(N,WA4)
-C
-C COMPUTE THE SCALED ACTUAL REDUCTION.
-C
- ACTRED = -ONE
- IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
-C
-C COMPUTE THE SCALED PREDICTED REDUCTION.
-C
- L = 1
- DO 220 I = 1, N
- SUM = ZERO
- DO 210 J = I, N
- SUM = SUM + R(L)*WA1(J)
- L = L + 1
- 210 CONTINUE
- WA3(I) = QTF(I) + SUM
- 220 CONTINUE
- TEMP = ENORM(N,WA3)
- PRERED = ZERO
- IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2
-C
-C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
-C REDUCTION.
-C
- RATIO = ZERO
- IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED
-C
-C UPDATE THE STEP BOUND.
-C
- IF (RATIO .GE. P1) GO TO 230
- NCSUC = 0
- NCFAIL = NCFAIL + 1
- DELTA = P5*DELTA
- GO TO 240
- 230 CONTINUE
- NCFAIL = 0
- NCSUC = NCSUC + 1
- IF (RATIO .GE. P5 .OR. NCSUC .GT. 1)
- * DELTA = DMAX1(DELTA,PNORM/P5)
- IF (DABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5
- 240 CONTINUE
-C
-C TEST FOR SUCCESSFUL ITERATION.
-C
- IF (RATIO .LT. P0001) GO TO 260
-C
-C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
-C
- DO 250 J = 1, N
- X(J) = WA2(J)
- WA2(J) = DIAG(J)*X(J)
- FVEC(J) = WA4(J)
- 250 CONTINUE
- XNORM = ENORM(N,WA2)
- FNORM = FNORM1
- ITER = ITER + 1
- 260 CONTINUE
-C
-C DETERMINE THE PROGRESS OF THE ITERATION.
-C
- NSLOW1 = NSLOW1 + 1
- IF (ACTRED .GE. P001) NSLOW1 = 0
- IF (JEVAL) NSLOW2 = NSLOW2 + 1
- IF (ACTRED .GE. P1) NSLOW2 = 0
-C
-C TEST FOR CONVERGENCE.
-C
- IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1
- IF (INFO .NE. 0) GO TO 300
-C
-C TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
-C
- IF (NFEV .GE. MAXFEV) INFO = 2
- IF (P1*DMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3
- IF (NSLOW2 .EQ. 5) INFO = 4
- IF (NSLOW1 .EQ. 10) INFO = 5
- IF (INFO .NE. 0) GO TO 300
-C
-C CRITERION FOR RECALCULATING JACOBIAN.
-C
- IF (NCFAIL .EQ. 2) GO TO 290
-C
-C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN
-C AND UPDATE QTF IF NECESSARY.
-C
- DO 280 J = 1, N
- SUM = ZERO
- DO 270 I = 1, N
- SUM = SUM + FJAC(I,J)*WA4(I)
- 270 CONTINUE
- WA2(J) = (SUM - WA3(J))/PNORM
- WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM)
- IF (RATIO .GE. P0001) QTF(J) = SUM
- 280 CONTINUE
-C
-C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.
-C
- CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)
- CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)
- CALL R1MPYQ(1,N,QTF,1,WA2,WA3)
-C
-C END OF THE INNER LOOP.
-C
- JEVAL = .FALSE.
- GO TO 180
- 290 CONTINUE
-C
-C END OF THE OUTER LOOP.
-C
- GO TO 30
- 300 CONTINUE
-C
-C TERMINATION, EITHER NORMAL OR USER IMPOSED.
-C
- IF (IFLAG .LT. 0) INFO = IFLAG
- IFLAG = 0
- IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
- RETURN
-C
-C LAST CARD OF SUBROUTINE HYBRJ.
-C
- END
diff --git a/libcruft/minpack/hybrj1.f b/libcruft/minpack/hybrj1.f
deleted file mode 100644
--- a/libcruft/minpack/hybrj1.f
+++ /dev/null
@@ -1,127 +0,0 @@
- SUBROUTINE HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA)
- INTEGER N,LDFJAC,INFO,LWA
- DOUBLE PRECISION TOL
- DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N),WA(LWA)
- EXTERNAL FCN
-C **********
-C
-C SUBROUTINE HYBRJ1
-C
-C THE PURPOSE OF HYBRJ1 IS TO FIND A ZERO OF A SYSTEM OF
-C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
-C OF THE POWELL HYBRID METHOD. THIS IS DONE BY USING THE
-C MORE GENERAL NONLINEAR EQUATION SOLVER HYBRJ. THE USER
-C MUST PROVIDE A SUBROUTINE WHICH CALCULATES THE FUNCTIONS
-C AND THE JACOBIAN.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA)
-C
-C WHERE
-C
-C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
-C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST
-C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER
-C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
-C
-C SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
-C INTEGER N,LDFJAC,IFLAG
-C DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N)
-C ----------
-C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND
-C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC.
-C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND
-C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC.
-C ---------
-C RETURN
-C END
-C
-C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
-C THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ1.
-C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF FUNCTIONS AND VARIABLES.
-C
-C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
-C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
-C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
-C
-C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
-C THE FUNCTIONS EVALUATED AT THE OUTPUT X.
-C
-C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
-C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION
-C OF THE FINAL APPROXIMATE JACOBIAN.
-C
-C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
-C
-C TOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
-C WHEN THE ALGORITHM ESTIMATES THAT THE RELATIVE ERROR
-C BETWEEN X AND THE SOLUTION IS AT MOST TOL.
-C
-C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
-C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
-C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
-C INFO IS SET AS FOLLOWS.
-C
-C INFO = 0 IMPROPER INPUT PARAMETERS.
-C
-C INFO = 1 ALGORITHM ESTIMATES THAT THE RELATIVE ERROR
-C BETWEEN X AND THE SOLUTION IS AT MOST TOL.
-C
-C INFO = 2 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS
-C REACHED 100*(N+1).
-C
-C INFO = 3 TOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
-C THE APPROXIMATE SOLUTION X IS POSSIBLE.
-C
-C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS.
-C
-C WA IS A WORK ARRAY OF LENGTH LWA.
-C
-C LWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(N+13))/2.
-C
-C SUBPROGRAMS CALLED
-C
-C USER-SUPPLIED ...... FCN
-C
-C MINPACK-SUPPLIED ... HYBRJ
-C
-C MINPACK. VERSION OF JULY 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER J,LR,MAXFEV,MODE,NFEV,NJEV,NPRINT
- DOUBLE PRECISION FACTOR,ONE,XTOL,ZERO
- DATA FACTOR,ONE,ZERO /1.0D2,1.0D0,0.0D0/
- INFO = 0
-C
-C CHECK THE INPUT PARAMETERS FOR ERRORS.
-C
- IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. TOL .LT. ZERO
- * .OR. LWA .LT. (N*(N + 13))/2) GO TO 20
-C
-C CALL HYBRJ.
-C
- MAXFEV = 100*(N + 1)
- XTOL = TOL
- MODE = 2
- DO 10 J = 1, N
- WA(J) = ONE
- 10 CONTINUE
- NPRINT = 0
- LR = (N*(N + 1))/2
- CALL HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,WA(1),MODE,
- * FACTOR,NPRINT,INFO,NFEV,NJEV,WA(6*N+1),LR,WA(N+1),
- * WA(2*N+1),WA(3*N+1),WA(4*N+1),WA(5*N+1))
- IF (INFO .EQ. 5) INFO = 4
- 20 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE HYBRJ1.
-C
- END
diff --git a/libcruft/minpack/qform.f b/libcruft/minpack/qform.f
deleted file mode 100644
--- a/libcruft/minpack/qform.f
+++ /dev/null
@@ -1,95 +0,0 @@
- SUBROUTINE QFORM(M,N,Q,LDQ,WA)
- INTEGER M,N,LDQ
- DOUBLE PRECISION Q(LDQ,M),WA(M)
-C **********
-C
-C SUBROUTINE QFORM
-C
-C THIS SUBROUTINE PROCEEDS FROM THE COMPUTED QR FACTORIZATION OF
-C AN M BY N MATRIX A TO ACCUMULATE THE M BY M ORTHOGONAL MATRIX
-C Q FROM ITS FACTORED FORM.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE QFORM(M,N,Q,LDQ,WA)
-C
-C WHERE
-C
-C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF ROWS OF A AND THE ORDER OF Q.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF COLUMNS OF A.
-C
-C Q IS AN M BY M ARRAY. ON INPUT THE FULL LOWER TRAPEZOID IN
-C THE FIRST MIN(M,N) COLUMNS OF Q CONTAINS THE FACTORED FORM.
-C ON OUTPUT Q HAS BEEN ACCUMULATED INTO A SQUARE MATRIX.
-C
-C LDQ IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY Q.
-C
-C WA IS A WORK ARRAY OF LENGTH M.
-C
-C SUBPROGRAMS CALLED
-C
-C FORTRAN-SUPPLIED ... MIN0
-C
-C MINPACK. VERSION OF JANUARY 1979.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,J,JM1,K,L,MINMN,NP1
- DOUBLE PRECISION ONE,SUM,TEMP,ZERO
- DATA ONE,ZERO /1.0D0,0.0D0/
-C
-C ZERO OUT UPPER TRIANGLE OF Q IN THE FIRST MIN(M,N) COLUMNS.
-C
- MINMN = MIN0(M,N)
- IF (MINMN .LT. 2) GO TO 30
- DO 20 J = 2, MINMN
- JM1 = J - 1
- DO 10 I = 1, JM1
- Q(I,J) = ZERO
- 10 CONTINUE
- 20 CONTINUE
- 30 CONTINUE
-C
-C INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX.
-C
- NP1 = N + 1
- IF (M .LT. NP1) GO TO 60
- DO 50 J = NP1, M
- DO 40 I = 1, M
- Q(I,J) = ZERO
- 40 CONTINUE
- Q(J,J) = ONE
- 50 CONTINUE
- 60 CONTINUE
-C
-C ACCUMULATE Q FROM ITS FACTORED FORM.
-C
- DO 120 L = 1, MINMN
- K = MINMN - L + 1
- DO 70 I = K, M
- WA(I) = Q(I,K)
- Q(I,K) = ZERO
- 70 CONTINUE
- Q(K,K) = ONE
- IF (WA(K) .EQ. ZERO) GO TO 110
- DO 100 J = K, M
- SUM = ZERO
- DO 80 I = K, M
- SUM = SUM + Q(I,J)*WA(I)
- 80 CONTINUE
- TEMP = SUM/WA(K)
- DO 90 I = K, M
- Q(I,J) = Q(I,J) - TEMP*WA(I)
- 90 CONTINUE
- 100 CONTINUE
- 110 CONTINUE
- 120 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE QFORM.
-C
- END
diff --git a/libcruft/minpack/qrfac.f b/libcruft/minpack/qrfac.f
deleted file mode 100644
--- a/libcruft/minpack/qrfac.f
+++ /dev/null
@@ -1,164 +0,0 @@
- SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
- INTEGER M,N,LDA,LIPVT
- INTEGER IPVT(LIPVT)
- LOGICAL PIVOT
- DOUBLE PRECISION A(LDA,N),SIGMA(N),ACNORM(N),WA(N)
-C **********
-C
-C SUBROUTINE QRFAC
-C
-C THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS WITH COLUMN
-C PIVOTING (OPTIONAL) TO COMPUTE A QR FACTORIZATION OF THE
-C M BY N MATRIX A. THAT IS, QRFAC DETERMINES AN ORTHOGONAL
-C MATRIX Q, A PERMUTATION MATRIX P, AND AN UPPER TRAPEZOIDAL
-C MATRIX R WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE,
-C SUCH THAT A*P = Q*R. THE HOUSEHOLDER TRANSFORMATION FOR
-C COLUMN K, K = 1,2,...,MIN(M,N), IS OF THE FORM
-C
-C T
-C I - (1/U(K))*U*U
-C
-C WHERE U HAS ZEROS IN THE FIRST K-1 POSITIONS. THE FORM OF
-C THIS TRANSFORMATION AND THE METHOD OF PIVOTING FIRST
-C APPEARED IN THE CORRESPONDING LINPACK SUBROUTINE.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
-C
-C WHERE
-C
-C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF ROWS OF A.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF COLUMNS OF A.
-C
-C A IS AN M BY N ARRAY. ON INPUT A CONTAINS THE MATRIX FOR
-C WHICH THE QR FACTORIZATION IS TO BE COMPUTED. ON OUTPUT
-C THE STRICT UPPER TRAPEZOIDAL PART OF A CONTAINS THE STRICT
-C UPPER TRAPEZOIDAL PART OF R, AND THE LOWER TRAPEZOIDAL
-C PART OF A CONTAINS A FACTORED FORM OF Q (THE NON-TRIVIAL
-C ELEMENTS OF THE U VECTORS DESCRIBED ABOVE).
-C
-C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.
-C
-C PIVOT IS A LOGICAL INPUT VARIABLE. IF PIVOT IS SET TRUE,
-C THEN COLUMN PIVOTING IS ENFORCED. IF PIVOT IS SET FALSE,
-C THEN NO COLUMN PIVOTING IS DONE.
-C
-C IPVT IS AN INTEGER OUTPUT ARRAY OF LENGTH LIPVT. IPVT
-C DEFINES THE PERMUTATION MATRIX P SUCH THAT A*P = Q*R.
-C COLUMN J OF P IS COLUMN IPVT(J) OF THE IDENTITY MATRIX.
-C IF PIVOT IS FALSE, IPVT IS NOT REFERENCED.
-C
-C LIPVT IS A POSITIVE INTEGER INPUT VARIABLE. IF PIVOT IS FALSE,
-C THEN LIPVT MAY BE AS SMALL AS 1. IF PIVOT IS TRUE, THEN
-C LIPVT MUST BE AT LEAST N.
-C
-C SIGMA IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE
-C DIAGONAL ELEMENTS OF R.
-C
-C ACNORM IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE
-C NORMS OF THE CORRESPONDING COLUMNS OF THE INPUT MATRIX A.
-C IF THIS INFORMATION IS NOT NEEDED, THEN ACNORM CAN COINCIDE
-C WITH SIGMA.
-C
-C WA IS A WORK ARRAY OF LENGTH N. IF PIVOT IS FALSE, THEN WA
-C CAN COINCIDE WITH SIGMA.
-C
-C SUBPROGRAMS CALLED
-C
-C MINPACK-SUPPLIED ... DPMPAR,ENORM
-C
-C FORTRAN-SUPPLIED ... DMAX1,DSQRT,MIN0
-C
-C MINPACK. VERSION OF DECEMBER 1978.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,J,JP1,K,KMAX,MINMN
- DOUBLE PRECISION AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
- DOUBLE PRECISION DPMPAR,ENORM
- DATA ONE,P05,ZERO /1.0D0,5.0D-2,0.0D0/
-C
-C EPSMCH IS THE MACHINE PRECISION.
-C
- EPSMCH = DPMPAR(1)
-C
-C COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
-C
- DO 10 J = 1, N
- ACNORM(J) = ENORM(M,A(1,J))
- SIGMA(J) = ACNORM(J)
- WA(J) = SIGMA(J)
- IF (PIVOT) IPVT(J) = J
- 10 CONTINUE
-C
-C REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
-C
- MINMN = MIN0(M,N)
- DO 110 J = 1, MINMN
- IF (.NOT.PIVOT) GO TO 40
-C
-C BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
-C
- KMAX = J
- DO 20 K = J, N
- IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
- 20 CONTINUE
- IF (KMAX .EQ. J) GO TO 40
- DO 30 I = 1, M
- TEMP = A(I,J)
- A(I,J) = A(I,KMAX)
- A(I,KMAX) = TEMP
- 30 CONTINUE
- SIGMA(KMAX) = SIGMA(J)
- WA(KMAX) = WA(J)
- K = IPVT(J)
- IPVT(J) = IPVT(KMAX)
- IPVT(KMAX) = K
- 40 CONTINUE
-C
-C COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
-C J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
-C
- AJNORM = ENORM(M-J+1,A(J,J))
- IF (AJNORM .EQ. ZERO) GO TO 100
- IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
- DO 50 I = J, M
- A(I,J) = A(I,J)/AJNORM
- 50 CONTINUE
- A(J,J) = A(J,J) + ONE
-C
-C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
-C AND UPDATE THE NORMS.
-C
- JP1 = J + 1
- IF (N .LT. JP1) GO TO 100
- DO 90 K = JP1, N
- SUM = ZERO
- DO 60 I = J, M
- SUM = SUM + A(I,J)*A(I,K)
- 60 CONTINUE
- TEMP = SUM/A(J,J)
- DO 70 I = J, M
- A(I,K) = A(I,K) - TEMP*A(I,J)
- 70 CONTINUE
- IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
- TEMP = A(J,K)/SIGMA(K)
- SIGMA(K) = SIGMA(K)*DSQRT(DMAX1(ZERO,ONE-TEMP**2))
- IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
- SIGMA(K) = ENORM(M-J,A(JP1,K))
- WA(K) = SIGMA(K)
- 80 CONTINUE
- 90 CONTINUE
- 100 CONTINUE
- SIGMA(J) = -AJNORM
- 110 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE QRFAC.
-C
- END
diff --git a/libcruft/minpack/r1mpyq.f b/libcruft/minpack/r1mpyq.f
deleted file mode 100644
--- a/libcruft/minpack/r1mpyq.f
+++ /dev/null
@@ -1,92 +0,0 @@
- SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
- INTEGER M,N,LDA
- DOUBLE PRECISION A(LDA,N),V(N),W(N)
-C **********
-C
-C SUBROUTINE R1MPYQ
-C
-C GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE
-C Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS
-C
-C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
-C
-C AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH
-C ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY.
-C Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE
-C GV, GW ROTATIONS IS SUPPLIED.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
-C
-C WHERE
-C
-C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF ROWS OF A.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF COLUMNS OF A.
-C
-C A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX
-C TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q
-C DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A.
-C
-C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
-C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.
-C
-C V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE
-C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I)
-C DESCRIBED ABOVE.
-C
-C W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE
-C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I)
-C DESCRIBED ABOVE.
-C
-C SUBROUTINES CALLED
-C
-C FORTRAN-SUPPLIED ... DABS,DSQRT
-C
-C MINPACK. VERSION OF DECEMBER 1978.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
-C
-C **********
- INTEGER I,J,NMJ,NM1
- DOUBLE PRECISION COS,ONE,SIN,TEMP
- DATA ONE /1.0D0/
-C
-C APPLY THE FIRST SET OF GIVENS ROTATIONS TO A.
-C
- NM1 = N - 1
- IF (NM1 .LT. 1) GO TO 50
- DO 20 NMJ = 1, NM1
- J = N - NMJ
- IF (DABS(V(J)) .GT. ONE) COS = ONE/V(J)
- IF (DABS(V(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
- IF (DABS(V(J)) .LE. ONE) SIN = V(J)
- IF (DABS(V(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
- DO 10 I = 1, M
- TEMP = COS*A(I,J) - SIN*A(I,N)
- A(I,N) = SIN*A(I,J) + COS*A(I,N)
- A(I,J) = TEMP
- 10 CONTINUE
- 20 CONTINUE
-C
-C APPLY THE SECOND SET OF GIVENS ROTATIONS TO A.
-C
- DO 40 J = 1, NM1
- IF (DABS(W(J)) .GT. ONE) COS = ONE/W(J)
- IF (DABS(W(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
- IF (DABS(W(J)) .LE. ONE) SIN = W(J)
- IF (DABS(W(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
- DO 30 I = 1, M
- TEMP = COS*A(I,J) + SIN*A(I,N)
- A(I,N) = -SIN*A(I,J) + COS*A(I,N)
- A(I,J) = TEMP
- 30 CONTINUE
- 40 CONTINUE
- 50 CONTINUE
- RETURN
-C
-C LAST CARD OF SUBROUTINE R1MPYQ.
-C
- END
diff --git a/libcruft/minpack/r1updt.f b/libcruft/minpack/r1updt.f
deleted file mode 100644
--- a/libcruft/minpack/r1updt.f
+++ /dev/null
@@ -1,207 +0,0 @@
- SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)
- INTEGER M,N,LS
- LOGICAL SING
- DOUBLE PRECISION S(LS),U(M),V(N),W(M)
-C **********
-C
-C SUBROUTINE R1UPDT
-C
-C GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U,
-C AND AN N-VECTOR V, THE PROBLEM IS TO DETERMINE AN
-C ORTHOGONAL MATRIX Q SUCH THAT
-C
-C T
-C (S + U*V )*Q
-C
-C IS AGAIN LOWER TRAPEZOIDAL.
-C
-C THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1)
-C TRANSFORMATIONS
-C
-C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
-C
-C WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE
-C WHICH ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES,
-C RESPECTIVELY. Q ITSELF IS NOT ACCUMULATED, RATHER THE
-C INFORMATION TO RECOVER THE GV, GW ROTATIONS IS RETURNED.
-C
-C THE SUBROUTINE STATEMENT IS
-C
-C SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)
-C
-C WHERE
-C
-C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF ROWS OF S.
-C
-C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
-C OF COLUMNS OF S. N MUST NOT EXCEED M.
-C
-C S IS AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER
-C TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS
-C THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE.
-C
-C LS IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
-C (N*(2*M-N+1))/2.
-C
-C U IS AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE
-C VECTOR U.
-C
-C V IS AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR
-C V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO
-C RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE.
-C
-C W IS AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION
-C NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED
-C ABOVE.
-C
-C SING IS A LOGICAL OUTPUT VARIABLE. SING IS SET TRUE IF ANY
-C OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE
-C V IS SET FALSE.
-C
-C SUBPROGRAMS CALLED
-C
-C MINPACK-SUPPLIED ... DPMPAR
-C
-C FORTRAN-SUPPLIED ... DABS,DSQRT
-C
-C MINPACK. VERSION OF DECEMBER 1978.
-C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE,
-C JOHN L. NAZARETH
-C
-C **********
- INTEGER I,J,JJ,L,NMJ,NM1
- DOUBLE PRECISION COS,COTAN,GIANT,ONE,P5,P25,SIN,TAN,TAU,TEMP,
- * ZERO
- DOUBLE PRECISION DPMPAR
- DATA ONE,P5,P25,ZERO /1.0D0,5.0D-1,2.5D-1,0.0D0/
-C
-C GIANT IS THE LARGEST MAGNITUDE.
-C
- GIANT = DPMPAR(3)
-C
-C INITIALIZE THE DIAGONAL ELEMENT POINTER.
-C
- JJ = (N*(2*M - N + 1))/2 - (M - N)
-C
-C MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W.
-C
- L = JJ
- DO 10 I = N, M
- W(I) = S(L)
- L = L + 1
- 10 CONTINUE
-C
-C ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR
-C IN SUCH A WAY THAT A SPIKE IS INTRODUCED INTO W.
-C
- NM1 = N - 1
- IF (NM1 .LT. 1) GO TO 70
- DO 60 NMJ = 1, NM1
- J = N - NMJ
- JJ = JJ - (M - J + 1)
- W(J) = ZERO
- IF (V(J) .EQ. ZERO) GO TO 50
-C
-C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
-C J-TH ELEMENT OF V.
-C
- IF (DABS(V(N)) .GE. DABS(V(J))) GO TO 20
- COTAN = V(N)/V(J)
- SIN = P5/DSQRT(P25+P25*COTAN**2)
- COS = SIN*COTAN
- TAU = ONE
- IF (DABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
- GO TO 30
- 20 CONTINUE
- TAN = V(J)/V(N)
- COS = P5/DSQRT(P25+P25*TAN**2)
- SIN = COS*TAN
- TAU = SIN
- 30 CONTINUE
-C
-C APPLY THE TRANSFORMATION TO V AND STORE THE INFORMATION
-C NECESSARY TO RECOVER THE GIVENS ROTATION.
-C
- V(N) = SIN*V(J) + COS*V(N)
- V(J) = TAU
-C
-C APPLY THE TRANSFORMATION TO S AND EXTEND THE SPIKE IN W.
-C
- L = JJ
- DO 40 I = J, M
- TEMP = COS*S(L) - SIN*W(I)
- W(I) = SIN*S(L) + COS*W(I)
- S(L) = TEMP
- L = L + 1
- 40 CONTINUE
- 50 CONTINUE
- 60 CONTINUE
- 70 CONTINUE
-C
-C ADD THE SPIKE FROM THE RANK 1 UPDATE TO W.
-C
- DO 80 I = 1, M
- W(I) = W(I) + V(N)*U(I)
- 80 CONTINUE
-C
-C ELIMINATE THE SPIKE.
-C
- SING = .FALSE.
- IF (NM1 .LT. 1) GO TO 140
- DO 130 J = 1, NM1
- IF (W(J) .EQ. ZERO) GO TO 120
-C
-C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
-C J-TH ELEMENT OF THE SPIKE.
-C
- IF (DABS(S(JJ)) .GE. DABS(W(J))) GO TO 90
- COTAN = S(JJ)/W(J)
- SIN = P5/DSQRT(P25+P25*COTAN**2)
- COS = SIN*COTAN
- TAU = ONE
- IF (DABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
- GO TO 100
- 90 CONTINUE
- TAN = W(J)/S(JJ)
- COS = P5/DSQRT(P25+P25*TAN**2)
- SIN = COS*TAN
- TAU = SIN
- 100 CONTINUE
-C
-C APPLY THE TRANSFORMATION TO S AND REDUCE THE SPIKE IN W.
-C
- L = JJ
- DO 110 I = J, M
- TEMP = COS*S(L) + SIN*W(I)
- W(I) = -SIN*S(L) + COS*W(I)
- S(L) = TEMP
- L = L + 1
- 110 CONTINUE
-C
-C STORE THE INFORMATION NECESSARY TO RECOVER THE
-C GIVENS ROTATION.
-C
- W(J) = TAU
- 120 CONTINUE
-C
-C TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S.
-C
- IF (S(JJ) .EQ. ZERO) SING = .TRUE.
- JJ = JJ + (M - J + 1)
- 130 CONTINUE
- 140 CONTINUE
-C
-C MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S.
-C
- L = JJ
- DO 150 I = N, M
- S(L) = W(I)
- L = L + 1
- 150 CONTINUE
- IF (S(JJ) .EQ. ZERO) SING = .TRUE.
- RETURN
-C
-C LAST CARD OF SUBROUTINE R1UPDT.
-C
- END
diff --git a/liboctave/Makefile.in b/liboctave/Makefile.in
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -73,14 +73,14 @@
SPARSE_MX_OP_INC := $(shell $(AWK) -f $(srcdir)/sparse-mk-ops.awk prefix=smx list_h_files=1 $(srcdir)/sparse-mx-ops)
-OPT_BASE := $(addsuffix -opts, DASPK DASRT DASSL LSODE NLEqn Quad)
+OPT_BASE := $(addsuffix -opts, DASPK DASRT DASSL LSODE Quad)
OPT_IN := $(addsuffix .in, $(OPT_BASE))
OPT_INC := $(addsuffix .h, $(OPT_BASE))
INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DAERT.h \
DAERTFunc.h DASPK.h DASRT.h DASSL.h FEGrid.h \
- LinConst.h LP.h LSODE.h NLConst.h NLEqn.h \
- NLFunc.h NLP.h ODE.h ODEFunc.h ODES.h ODESFunc.h \
+ LinConst.h LP.h LSODE.h \
+ ODE.h ODEFunc.h ODES.h ODESFunc.h \
Objective.h QP.h Quad.h Range.h base-dae.h \
base-de.h base-min.h byte-swap.h cmd-edit.h cmd-hist.h \
data-conv.h dir-ops.h file-ops.h file-stat.h functor.h getopt.h \
@@ -143,7 +143,7 @@
SPARSE_MX_OP_SRC := $(shell $(AWK) -f $(srcdir)/sparse-mk-ops.awk prefix=smx list_cc_files=1 $(srcdir)/sparse-mx-ops)
LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DASPK.cc DASRT.cc \
- DASSL.cc FEGrid.cc LinConst.cc LSODE.cc NLEqn.cc ODES.cc \
+ DASSL.cc FEGrid.cc LinConst.cc LSODE.cc ODES.cc \
Quad.cc Range.cc data-conv.cc dir-ops.cc \
file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \
lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \
diff --git a/liboctave/NLConst.h b/liboctave/NLConst.h
deleted file mode 100644
--- a/liboctave/NLConst.h
+++ /dev/null
@@ -1,72 +0,0 @@
-/*
-
-Copyright (C) 1993, 1994, 1995, 1996, 1997, 2002, 2004, 2005, 2007
- John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#if !defined (octave_NLConst_h)
-#define octave_NLConst_h 1
-
-class ColumnVector;
-
-#include "Bounds.h"
-#include "NLFunc.h"
-
-class
-NLConst : public Bounds, public NLFunc
-{
-public:
-
- NLConst (void)
- : Bounds (), NLFunc () { }
-
- NLConst (octave_idx_type n)
- : Bounds (n), NLFunc () { }
-
- NLConst (const ColumnVector& lb, const NLFunc f, const ColumnVector& ub)
- : Bounds (lb, ub), NLFunc (f) { }
-
- NLConst (const NLConst& a)
- : Bounds (a.lb, a.ub), NLFunc (a.fun, a.jac) { }
-
- NLConst& operator = (const NLConst& a)
- {
- if (this != &a)
- {
- Bounds::operator = (a);
- NLFunc::operator = (a);
- }
- return *this;
- }
-
- ~NLConst (void) { }
-
-private:
-
- void error (const char *msg);
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/liboctave/NLEqn-opts.in b/liboctave/NLEqn-opts.in
deleted file mode 100644
--- a/liboctave/NLEqn-opts.in
+++ /dev/null
@@ -1,43 +0,0 @@
-# Copyright (C) 2002 John W. Eaton
-#
-# This file is part of Octave.
-#
-# Octave 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 3 of the License, or (at
-# your option) any later version.
-#
-# Octave 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 Octave; see the file COPYING. If not, see
-# .
-
-CLASS = "NLEqn"
-
-FCN_NAME = "fsolve"
-
-INCLUDE = "dColVector.h"
-INCLUDE = "NLFunc.h"
-
-
-DOC_STRING
-When called with two arguments, this function allows you set options
-parameters for the function @code{fsolve}. Given one argument,
address@hidden returns the value of the corresponding option. If
-no arguments are supplied, the names of all the available options and
-their current values are displayed.
-END_DOC_STRING
-
-OPTION
- NAME = "tolerance"
- DOC_ITEM
-Nonnegative relative tolerance.
- END_DOC_ITEM
- TYPE = "double"
- INIT_VALUE = "::sqrt (DBL_EPSILON)"
- SET_EXPR = "(val > 0.0) ? val : ::sqrt (DBL_EPSILON)"
-END_OPTION
diff --git a/liboctave/NLEqn.cc b/liboctave/NLEqn.cc
deleted file mode 100644
--- a/liboctave/NLEqn.cc
+++ /dev/null
@@ -1,249 +0,0 @@
-/*
-
-Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004,
- 2005, 2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include
-#endif
-
-#include "NLEqn.h"
-#include "dMatrix.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "quit.h"
-
-typedef octave_idx_type (*hybrd1_fcn_ptr) (octave_idx_type*, double*, double*, octave_idx_type*);
-
-typedef octave_idx_type (*hybrj1_fcn_ptr) (octave_idx_type*, double*, double*, double*, octave_idx_type*, octave_idx_type*);
-
-extern "C"
-{
- F77_RET_T
- F77_FUNC (hybrd1, HYBRD1) (hybrd1_fcn_ptr, const octave_idx_type&, double*,
- double*, const double&, octave_idx_type&, double*,
- const octave_idx_type&);
-
-
- F77_RET_T
- F77_FUNC (hybrj1, HYBRJ1) (hybrj1_fcn_ptr, const octave_idx_type&, double*,
- double*, double*, const octave_idx_type&, const
- double&, octave_idx_type&, double*, const octave_idx_type&);
-}
-
-static NLFunc::nonlinear_fcn user_fun;
-static NLFunc::jacobian_fcn user_jac;
-
-// error handling
-
-void
-NLEqn::error (const char* msg)
-{
- (*current_liboctave_error_handler) ("fatal NLEqn error: %s", msg);
-}
-
-// Other operations
-
-octave_idx_type
-hybrd1_fcn (octave_idx_type *n, double *x, double *fvec, octave_idx_type *iflag)
-{
- BEGIN_INTERRUPT_WITH_EXCEPTIONS;
-
- octave_idx_type nn = *n;
- ColumnVector tmp_f (nn);
- ColumnVector tmp_x (nn);
-
- for (octave_idx_type i = 0; i < nn; i++)
- tmp_x.elem (i) = x[i];
-
- tmp_f = (*user_fun) (tmp_x);
-
- if (tmp_f.length () == 0)
- *iflag = -1;
- else
- {
- for (octave_idx_type i = 0; i < nn; i++)
- fvec[i] = tmp_f.elem (i);
- }
-
- END_INTERRUPT_WITH_EXCEPTIONS;
-
- return 0;
-}
-
-octave_idx_type
-hybrj1_fcn (octave_idx_type *n, double *x, double *fvec, double *fjac,
- octave_idx_type *ldfjac, octave_idx_type *iflag)
-{
- BEGIN_INTERRUPT_WITH_EXCEPTIONS;
-
- octave_idx_type nn = *n;
- ColumnVector tmp_x (nn);
-
- for (octave_idx_type i = 0; i < nn; i++)
- tmp_x.elem (i) = x[i];
-
- octave_idx_type flag = *iflag;
- if (flag == 1)
- {
- ColumnVector tmp_f (nn);
-
- tmp_f = (*user_fun) (tmp_x);
-
- if (tmp_f.length () == 0)
- *iflag = -1;
- else
- {
- for (octave_idx_type i = 0; i < nn; i++)
- fvec[i] = tmp_f.elem (i);
- }
- }
- else
- {
- Matrix tmp_fj (nn, nn);
-
- tmp_fj = (*user_jac) (tmp_x);
-
- if (tmp_fj.rows () == 0 || tmp_fj.columns () == 0)
- *iflag = -1;
- else
- {
- octave_idx_type ld = *ldfjac;
- for (octave_idx_type j = 0; j < nn; j++)
- for (octave_idx_type i = 0; i < nn; i++)
- fjac[j*ld+i] = tmp_fj.elem (i, j);
- }
- }
-
- END_INTERRUPT_WITH_EXCEPTIONS;
-
- return 0;
-}
-
-ColumnVector
-NLEqn::solve (octave_idx_type& info)
-{
- ColumnVector retval;
-
- octave_idx_type n = x.capacity ();
-
- if (n == 0)
- {
- error ("equation set not initialized");
- return retval;
- }
-
- double tol = tolerance ();
-
- retval = x;
- double *px = retval.fortran_vec ();
-
- user_fun = fun;
- user_jac = jac;
-
- if (jac)
- {
- Array fvec (n);
- double *pfvec = fvec.fortran_vec ();
-
- octave_idx_type lwa = (n*(n+13))/2;
- Array wa (lwa);
- double *pwa = wa.fortran_vec ();
-
- Array fjac (n*n);
- double *pfjac = fjac.fortran_vec ();
-
- F77_XFCN (hybrj1, HYBRJ1, (hybrj1_fcn, n, px, pfvec, pfjac, n,
- tol, info, pwa, lwa));
-
- solution_status = info;
-
- fval = ColumnVector (fvec);
- }
- else
- {
- Array fvec (n);
- double *pfvec = fvec.fortran_vec ();
-
- octave_idx_type lwa = (n*(3*n+13))/2;
- Array wa (lwa);
- double *pwa = wa.fortran_vec ();
-
- F77_XFCN (hybrd1, HYBRD1, (hybrd1_fcn, n, px, pfvec, tol, info,
- pwa, lwa));
-
- solution_status = info;
-
- fval = ColumnVector (fvec);
- }
-
- return retval;
-}
-
-std::string
-NLEqn::error_message (void) const
-{
- std::string retval;
-
- std::string prefix;
-
- octave_idx_type info = solution_status;
- if (info < 0)
- info = -info;
-
- switch (info)
- {
- case 0:
- retval = "improper input parameters";
- break;
-
- case 1:
- retval = "solution converged within specified tolerance";
- break;
-
- case 2:
- retval = "number of function calls exceeded limit";
- break;
-
- case 3:
- retval = "no further improvement possible (tolerance may be too small)";
- break;
-
- case 4:
- retval = "iteration is not making good progress";
- break;
-
- default:
- retval = "unknown error state";
- break;
- }
-
- if (solution_status < 0)
- retval = std::string ("user requested termination: ") + retval;
-
- return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/liboctave/NLEqn.h b/liboctave/NLEqn.h
deleted file mode 100644
--- a/liboctave/NLEqn.h
+++ /dev/null
@@ -1,117 +0,0 @@
-/*
-
-Copyright (C) 1993, 1994, 1995, 1996, 1997, 2002, 2004, 2005, 2006,
- 2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#if !defined (octave_NLEqn_h)
-#define octave_NLEqn_h 1
-
-#include
-
-#include "NLEqn-opts.h"
-#include "lo-ieee.h"
-#include "lo-math.h"
-
-class
-OCTAVE_API
-NLEqn : public NLFunc, public NLEqn_options
-{
-public:
-
- NLEqn (void)
- : NLFunc (), NLEqn_options (), x (), fval (),
- solution_status (0) { }
-
- NLEqn (const ColumnVector& xx, const NLFunc f)
- : NLFunc (f), NLEqn_options (), x (xx), fval (x.numel (), octave_NaN),
- solution_status (0) { }
-
- NLEqn (const NLEqn& a)
- : NLFunc (a.fun, a.jac), NLEqn_options (), x (a.x), fval (a.fval),
- solution_status (a.solution_status) { }
-
- NLEqn& operator = (const NLEqn& a)
- {
- if (this != &a)
- {
- NLFunc::operator = (a);
- NLEqn_options::operator = (a);
-
- x = a.x;
- fval = a.fval;
- solution_status = a.solution_status;
- }
- return *this;
- }
-
- ~NLEqn (void) { }
-
- void set_states (const ColumnVector& xx) { x = xx; }
-
- ColumnVector states (void) const { return x; }
-
- octave_idx_type size (void) const { return x.capacity (); }
-
- ColumnVector solve (void)
- {
- octave_idx_type info;
- return solve (info);
- }
-
- ColumnVector solve (const ColumnVector& xvec)
- {
- set_states (xvec);
- octave_idx_type info;
- return solve (info);
- }
-
- ColumnVector solve (const ColumnVector& xvec, octave_idx_type& info)
- {
- set_states (xvec);
- return solve (info);
- }
-
- ColumnVector solve (octave_idx_type& info);
-
- octave_idx_type solution_state (void) const { return solution_status; }
-
- bool solution_ok (void) const { return solution_status == 1; }
-
- ColumnVector function_value (void) const { return fval; }
-
- std::string error_message (void) const;
-
-private:
-
- ColumnVector x;
- ColumnVector fval;
- octave_idx_type solution_status;
-
- void error (const char* msg);
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/liboctave/NLFunc.h b/liboctave/NLFunc.h
deleted file mode 100644
--- a/liboctave/NLFunc.h
+++ /dev/null
@@ -1,89 +0,0 @@
-/*
-
-Copyright (C) 1993, 1994, 1995, 1996, 1997, 2005, 2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#if !defined (octave_NLFunc_h)
-#define octave_NLFunc_h 1
-
-class ColumnVector;
-class Matrix;
-
-class
-NLFunc
-{
-public:
-
- typedef ColumnVector (*nonlinear_fcn) (const ColumnVector&);
- typedef Matrix (*jacobian_fcn) (const ColumnVector&);
-
- NLFunc (void)
- : fun (0), jac (0) { }
-
- NLFunc (const nonlinear_fcn f)
- : fun (f), jac (0) { }
-
- NLFunc (const nonlinear_fcn f, const jacobian_fcn j)
- : fun (f), jac (j) { }
-
- NLFunc (const NLFunc& a)
- : fun (a.fun), jac (a.jac) { }
-
- NLFunc& operator = (const NLFunc& a)
- {
- if (this != &a)
- {
- fun = a.fun;
- jac = a.jac;
- }
- return *this;
- }
-
- ~NLFunc (void) { }
-
- nonlinear_fcn function (void) const { return fun; }
-
- NLFunc& set_function (const nonlinear_fcn f)
- {
- fun = f;
- return *this;
- }
-
- jacobian_fcn jacobian_function (void) const { return jac; }
-
- NLFunc& set_jacobian_function (const jacobian_fcn j)
- {
- jac = j;
- return *this;
- }
-
-protected:
-
- nonlinear_fcn fun;
- jacobian_fcn jac;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/liboctave/NLP.h b/liboctave/NLP.h
deleted file mode 100644
--- a/liboctave/NLP.h
+++ /dev/null
@@ -1,111 +0,0 @@
-/*
-
-Copyright (C) 1993, 1994, 1995, 1996, 1997, 2005, 2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#if !defined (octave_NLP_h)
-#define octave_NLP_h 1
-
-#include "dColVector.h"
-#include "Objective.h"
-#include "Bounds.h"
-#include "LinConst.h"
-#include "NLConst.h"
-#include "base-min.h"
-
-class
-NLP : public base_minimizer
-{
-public:
-
- NLP (void)
- : base_minimizer (), phi (), bnds (), lc (), nlc () { }
-
- NLP (const ColumnVector& x, const Objective& obj)
- : base_minimizer (x), phi (obj), bnds (), lc (), nlc () { }
-
- NLP (const ColumnVector& x, const Objective& obj, const Bounds& b)
- : base_minimizer (x), phi (obj), bnds (b), lc (), nlc () { }
-
- NLP (const ColumnVector& x, const Objective& obj, const Bounds& b,
- const LinConst& l)
- : base_minimizer (x), phi (obj), bnds (b), lc (l), nlc () { }
-
- NLP (const ColumnVector& x, const Objective& obj, const Bounds& b,
- const LinConst& l, const NLConst& nl)
- : base_minimizer (x), phi (obj), bnds (b), lc (l), nlc (nl) { }
-
- NLP (const ColumnVector& x, const Objective& obj, const LinConst& l)
- : base_minimizer (x), phi (obj), bnds (), lc (l), nlc () { }
-
- NLP (const ColumnVector& x, const Objective& obj, const LinConst& l,
- const NLConst& nl)
- : base_minimizer (x), phi (obj), bnds (), lc (l), nlc (nl) { }
-
- NLP (const ColumnVector& x, const Objective& obj, const NLConst& nl)
- : base_minimizer (x), phi (obj), bnds (), lc (), nlc (nl) { }
-
- NLP (const ColumnVector& x, const Objective& obj, const Bounds& b,
- const NLConst& nl)
- : base_minimizer (x), phi (obj), bnds (b), lc (), nlc (nl) { }
-
- NLP (const NLP& a)
- : base_minimizer (a), phi (a.phi), bnds (a.bnds), lc (a.lc), nlc (a.nlc)
- { }
-
- NLP& operator = (const NLP& a)
- {
- if (this != &a)
- {
- base_minimizer::operator = (a);
-
- phi = a.phi;
- bnds = a.bnds;
- lc = a.lc;
- nlc = a.nlc;
- }
- return *this;
- }
-
- virtual ~NLP (void) { }
-
- Objective objective (void) const { return phi; }
-
- Bounds bounds (void) const { return bnds; }
-
- LinConst linear_constraints (void) const { return lc; }
-
- NLConst nonlinear_constraints (void) const { return nlc; }
-
-protected:
-
- Objective phi;
- Bounds bnds;
- LinConst lc;
- NLConst nlc;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -0,0 +1,7 @@
+2008-09-28 Jaroslav Hajek
+
+ * optimization/__fdjac__.m: New function file.
+ * optimization/__dogleg__.m: New function file.
+ * optimization/fsolve.m: New function file.
+ * optimization/Makefile.in: Include the new sources.
+
diff --git a/scripts/optimization/Makefile.in b/scripts/optimization/Makefile.in
--- a/scripts/optimization/Makefile.in
+++ b/scripts/optimization/Makefile.in
@@ -35,6 +35,9 @@
SOURCES = \
__fsolve_defopts__.m \
fzero.m \
+ __fdjac__.m \
+ __dogleg__.m \
+ fsolve.m \
glpk.m \
glpkmex.m \
lsqnonneg.m \
diff --git a/scripts/optimization/__dogleg__.m b/scripts/optimization/__dogleg__.m
new file mode 100644
--- /dev/null
+++ b/scripts/optimization/__dogleg__.m
@@ -0,0 +1,60 @@
+## Copyright (C) 2008 Jaroslav Hajek
+##
+## This file is part of Octave.
+##
+## Octave 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 3 of the License, or (at
+## your option) any later version.
+##
+## Octave 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 Octave; see the file COPYING. If not, see
+## .
+
+## -*- texinfo -*-
+## @deftypefn{Function File} address@hidden = __dogleg__ (@var{r}, @var{b}, @var{x}, @var{d}, @var{delta})
+## Solves the double dogleg trust-region problem:
+## Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
+## x being a convex combination of the gauss-newton and scaled gradient.
+## @end deftypefn
+
+## TODO: error checks
+## TODO: handle singularity, or leave it up to mldivide?
+
+function x = __dogleg__ (r, b, d, delta)
+ # get Gauss-Newton direction
+ x = r \ b;
+ xn = norm (d .* x);
+ if (xn > delta)
+ # GN is too big, get scaled gradient
+ s = (r' * b) ./ d;
+ sn = norm (s);
+ if (sn > 0)
+ # normalize and rescale
+ s = (s / sn) ./ d;
+ # get the line minimizer in s direction
+ tn = norm (r*s);
+ snm = (sn / tn) / tn;
+ if (snm < delta)
+ # get the dogleg path minimizer
+ bn = norm (b);
+ dxn = delta/xn; snmd = snm/delta;
+ t = (bn/sn) * (bn/xn) * snmd;
+ t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
+ alpha = dxn*(1-snmd^2) / t;
+ else
+ alpha = 0;
+ endif
+ else
+ snm = 0;
+ endif
+ # form the appropriate convex combination
+ x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
+ endif
+endfunction
+
diff --git a/scripts/optimization/__fdjac__.m b/scripts/optimization/__fdjac__.m
new file mode 100644
--- /dev/null
+++ b/scripts/optimization/__fdjac__.m
@@ -0,0 +1,38 @@
+## Copyright (C) 2008 Jaroslav Hajek
+##
+## This file is part of Octave.
+##
+## Octave 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 3 of the License, or (at
+## your option) any later version.
+##
+## Octave 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 Octave; see the file COPYING. If not, see
+## .
+
+## -*- texinfo -*-
+## @deftypefn{Function File} address@hidden, @var{fjac}]} = __fdjac__ (@var{fcn}, @var{x}, @var{err})
+## @end deftypefn
+
+function [fvec, fjac] = __fdjac__ (fcn, x, err = 0)
+ err = sqrt (max (eps, err));
+ fvec = fcn (x);
+ fv = fvec(:);
+ h = max (abs (x(:))*err, (fv*eps)/realmax);
+ h(h == 0) = err;
+ fjac = zeros (length (fv), numel (x));
+ for i = 1:numel (x)
+ x1 = x;
+ x1(i) += h(i);
+ fjac(:,i) = (fcn (x1)(:) - fv) / h(i);
+ endfor
+endfunction
+
+
+
diff --git a/scripts/optimization/fsolve.m b/scripts/optimization/fsolve.m
new file mode 100644
--- /dev/null
+++ b/scripts/optimization/fsolve.m
@@ -0,0 +1,258 @@
+## Copyright (C) 2008 VZLU Prague, a.s.
+##
+## This file is part of Octave.
+##
+## Octave 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 3 of the License, or (at
+## your option) any later version.
+##
+## Octave 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 Octave; see the file COPYING. If not, see
+## .
+##
+## Author: Jaroslav Hajek
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {} fsolve(@var{fcn},@var{x0},@var{options})
+## @deftypefnx{Function File} address@hidden, @var{fvec}, @var{info}]} = fsolve (@var{fcn}, @dots{})
+## @end deftypefn
+
+function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options)
+
+ if (nargin < 3)
+ options = struct ();
+ endif
+
+ xsiz = size (x0);
+ n = numel (x0);
+
+ has_jac = strcmp (optimget (options, "Jacobian", "off"), "on");
+ maxiter = optimget (options, "MaxIter", Inf);
+ maxfev = optimget (options, "MaxFunEvals", Inf);
+ outfcn = optimget (options, "OutputFcn");
+ funvalchk = strcmp (optimget (options, "FunValCheck", "off"), "on");
+
+ if (funvalchk)
+ # replace fun with a guarded version
+ fun = @(x) guarded_eval (fun, x);
+ endif
+
+ # These defaults are rather stringent. I think that normally, user prefers
+ # accuracy to performance.
+
+ macheps = eps (class (x0));
+
+ tolx = optimget (options, "TolX", 1e1*macheps);
+ tolf = optimget (options, "TolFun",1e2*macheps);
+
+ factor = 100;
+ # TODO: TypicalX corresponds to user scaling (???)
+ autodg = true;
+
+ niter = 1; nfev = 0;
+
+ x = x0(:);
+ info = 0;
+
+ # outer loop
+ while (niter < maxiter && nfev < maxfev && ! info)
+
+ # calc func value and jacobian (possibly via FD)
+ # handle arbitrary shapes of x and f and remember them
+ if (has_jac)
+ [fvec, fjac] = fcn (reshape (x, xsiz));
+ nfev ++;
+ else
+ [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
+ nfev += 1 + length (x);
+ endif
+ fsiz = size (fvec);
+ fvec = fvec(:);
+ fn = norm (fvec);
+
+ # get QR factorization of the jacobian
+ [q, r] = qr (fjac);
+
+ # get column norms, use them as scaling factor
+ jcn = norm (fjac, 'columns').';
+ if (niter == 1)
+ if (autodg)
+ dg = jcn;
+ dg(dg == 0) = 1;
+ endif
+ xn = norm (dg .* x);
+ delta = factor * xn;
+ endif
+
+ # rescale if necessary
+ if (autodg)
+ dg = max (dg, jcn);
+ endif
+
+ nfail = 0;
+ nsuc = 0;
+ # inner loop
+ while (niter <= maxiter && nfev < maxfev && ! info)
+
+ qtf = q'*fvec;
+
+ # get TR model (dogleg) minimizer
+ p = - __dogleg__ (r, qtf, dg, delta);
+ pn = norm (dg .* p);
+
+ if (niter == 1)
+ delta = min (delta, pn);
+ endif
+
+ fvec1 = fcn (reshape (x + p, xsiz)) (:);
+ fn1 = norm (fvec1);
+
+ if (fn1 < fn)
+ # scaled actual reduction
+ actred = 1 - (fn1/fn)^2;
+ else
+ actred = -1;
+ endif
+
+ # scaled predicted reduction, and ratio
+ w = qtf + r*p;
+ t = norm (w);
+ if (t < fn)
+ prered = 1 - (t/fn)^2;
+ ratio = actred / prered;
+ else
+ prered = 0;
+ ratio = 0;
+ endif
+
+ # update delta
+ if (ratio < 0.1)
+ nsuc = 0;
+ nfail ++;
+ delta *= 0.5;
+ if (delta <= sqrt (macheps)*xn)
+ # trust region became uselessly small
+ info = -3;
+ break;
+ endif
+ else
+ nfail = 0;
+ nsuc ++;
+ if (abs (1-ratio) <= 0.1)
+ delta = 2*pn;
+ elseif (ratio >= 0.5 || nsuc > 1)
+ delta = max (delta, 2*pn);
+ endif
+ endif
+
+ if (ratio >= 1e-4)
+ # successful iteration
+ x += p;
+ xn = norm (dg .* x);
+ fvec = fvec1;
+ fn = fn1;
+ niter ++;
+ endif
+
+ # Tests for termination conditions. A mysterious place, anything can
+ # happen if you change something here...
+
+ # The rule of thumb (which I'm not sure M*b is quite following) is that
+ # for a tolerance that depends on scaling, only 0 makes sense as a
+ # default value. But 0 usually means uselessly long iterations,
+ # so we need scaling-independent tolerances wherever possible.
+
+ # XXX: why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector of
+ # perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by tolf ~
+ # eps we demand as much accuracy as we can expect.
+ if (fn <= tolf*n*xn)
+ info = 1;
+ # The following tests done only after successful step.
+ elseif (actred > 0)
+ # This one is classic. Note that we use scaled variables again, but
+ # compare to scaled step, so nothing bad.
+ if (pn <= tolx*xn)
+ info = 2;
+ # Again a classic one. It seems weird to use the same tolf for two
+ # different tests, but that's what M*b manual appears to say.
+ elseif (actred < tolf)
+ info = 3
+ endif
+ endif
+
+ # criterion for recalculating jacobian
+ if (nfail == 2)
+ break;
+ endif
+
+ # compute the scaled Broyden update
+ u = (fvec1 - q*w) / pn;
+ v = dg .* ((dg .* p) / pn);
+
+ # update the QR factorization
+ [q, r] = qrupdate (q, r, u, v);
+
+ endwhile
+ endwhile
+
+ # restore original shapes
+ x = reshape (x, xsiz);
+ fvec = reshape (fvec, fsiz);
+
+ output.iterations = niter;
+ output.funcCount = niter + 2;
+
+endfunction
+
+# an assistant function that evaluates a function handle and checks for bad
+# results.
+function fx = guarded_eval (fun, x)
+ fx = fun (x);
+ if (! all (isreal (fx)))
+ error ("fsolve:notreal", "fsolve: non-real value encountered");
+ elseif (any (isnan (fx)))
+ error ("fsolve:isnan", "fsolve: NaN value encountered");
+ endif
+endfunction
+
+%!function retval = f (p)
+%! x = p(1);
+%! y = p(2);
+%! z = p(3);
+%! retval = zeros (3, 1);
+%! retval(1) = sin(x) + y**2 + log(z) - 7;
+%! retval(2) = 3*x + 2**y -z**3 + 1;
+%! retval(3) = x + y + z - 5;
+%!test
+%! x_opt = [ 0.599054;
+%! 2.395931;
+%! 2.005014 ];
+%! tol = 1.0e-5;
+%! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
+%! assert (info > 0);
+%! assert (norm (x - x_opt, 1) < tol);
+%! assert (norm (fval) < tol);
+
+%!function retval = f (p)
+%! x = p(1);
+%! y = p(2);
+%! z = p(3);
+%! w = p(4);
+%! retval = zeros (4, 1);
+%! retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
+%! retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
+%! retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
+%! retval(4) = x^2 + 2*y^3 + z - w - 4;
+%!test
+%! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
+%! tol = 1.0e-5;
+%! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
+%! assert (info > 0);
+%! assert (norm (x - x_opt, 1) < tol);
+%! assert (norm (fval) < tol);
diff --git a/src/DLD-FUNCTIONS/fsolve.cc b/src/DLD-FUNCTIONS/fsolve.cc
deleted file mode 100644
--- a/src/DLD-FUNCTIONS/fsolve.cc
+++ /dev/null
@@ -1,646 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006,
- 2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave 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 3 of the License, or (at your
-option) any later version.
-
-Octave 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 Octave; see the file COPYING. If not, see
-.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include
-#endif
-
-#include
-#include
-#include
-
-#include
-#include
-#include
-
-#include "dNDArray.h"
-#include "NLEqn.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-map.h"
-#include "oct-obj.h"
-#include "ov-fcn.h"
-#include "ov-cell.h"
-#include "pager.h"
-#include "unwind-prot.h"
-#include "utils.h"
-#include "variables.h"
-
-#include "NLEqn-opts.cc"
-
-// Global pointer for user defined function required by hybrd1.
-static octave_function *fsolve_fcn;
-
-// Global pointer for optional user defined jacobian function.
-static octave_function *fsolve_jac;
-
-// Original dimensions of X0.
-static dim_vector x_dims;
-
-// Have we warned about imaginary values returned from user function?
-static bool warned_fcn_imaginary = false;
-static bool warned_jac_imaginary = false;
-
-// Is this a recursive call?
-static int call_depth = 0;
-
-octave_idx_type
-hybrd_info_to_fsolve_info (octave_idx_type info)
-{
- switch (info)
- {
- case -1:
- case 1:
- break;
-
- case 2:
- info = 0;
- break;
-
- case 3:
- case 4:
- case 5:
- info = -2;
- break;
-
- default:
- {
- std::ostringstream buf;
- buf << "fsolve: unexpected value of INFO from MINPACK (= "
- << info << ")";
- std::string msg = buf.str ();
- warning (msg.c_str ());
- }
- break;
- }
-
- return info;
-}
-
-ColumnVector
-fsolve_user_function (const ColumnVector& x)
-{
- ColumnVector retval;
-
- octave_idx_type n = x.length ();
-
- octave_value_list args;
- args.resize (1);
-
- if (n > 1)
- {
- NDArray m (ArrayN (x, x_dims));
- octave_value vars (m);
- args(0) = vars;
- }
- else
- {
- double d = x (0);
- octave_value vars (d);
- args(0) = vars;
- }
-
- if (fsolve_fcn)
- {
- octave_value_list tmp = fsolve_fcn->do_multi_index_op (1, args);
-
- if (tmp.length () > 0 && tmp(0).is_defined ())
- {
- if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
- {
- warning ("fsolve: ignoring imaginary part returned from user-supplied function");
- warned_fcn_imaginary = true;
- }
-
- retval = ColumnVector (tmp(0).vector_value (false, true));
-
- if (error_state || retval.length () <= 0)
- gripe_user_supplied_eval ("fsolve");
- else if (retval.length () != x.length ())
- error ("fsolve: unable to solve non-square systems");
- }
- else
- gripe_user_supplied_eval ("fsolve");
- }
-
- return retval;
-}
-
-Matrix
-fsolve_user_jacobian (const ColumnVector& x)
-{
- Matrix retval;
-
- octave_idx_type n = x.length ();
-
- octave_value_list args;
- args.resize (1);
-
- if (n > 1)
- {
- NDArray m (ArrayN (x, x_dims));
- octave_value vars (m);
- args(0) = vars;
- }
- else
- {
- double d = x(0);
- octave_value vars (d);
- args(0) = vars;
- }
-
- if (fsolve_fcn)
- {
- octave_value_list tmp = fsolve_jac->do_multi_index_op (1, args);
-
- if (tmp.length () > 0 && tmp(0).is_defined ())
- {
- if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
- {
- warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function");
- warned_fcn_imaginary = true;
- }
-
- retval = tmp(0).matrix_value ();
-
- if (error_state || retval.length () <= 0)
- gripe_user_supplied_eval ("fsolve");
- else if (! (retval.rows () == x.length ()
- && retval.columns () == x.length ()))
- error ("fsolve: invalid Jacobian matrix dimensions");
- }
- else
- gripe_user_supplied_eval ("fsolve");
- }
-
- return retval;
-}
-
-static std::set
-make_unimplemented_options (void)
-{
- static bool initialized = false;
-
- static std::set options;
-
- if (! initialized)
- {
- initialized = true;
-
- options.insert ("Largescale");
- options.insert ("DerivativeCheck");
- options.insert ("Diagnostics");
- options.insert ("DiffMaxChange");
- options.insert ("DiffMinChange");
- options.insert ("Display");
- options.insert ("FunValCheck");
- options.insert ("Jacobian");
- options.insert ("MaxFunEvals");
- options.insert ("MaxIter");
- options.insert ("OutputFcn");
- options.insert ("PlotFcns");
- options.insert ("TolFun");
- options.insert ("TolX");
- options.insert ("TypicalX");
- options.insert ("JacobMult");
- options.insert ("JacobPattern");
- options.insert ("MaxPCGIter");
- options.insert ("PrecondBandwidth");
- options.insert ("TolPCG");
- options.insert ("NonlEqnAlgorithm");
- options.insert ("LineSearchType");
- }
-
- return options;
-}
-
-static void
-override_options (NLEqn_options& opts, const Octave_map& option_map)
-{
- string_vector keys = option_map.keys ();
-
- for (octave_idx_type i = 0; i < keys.length (); i++)
- {
- std::string key = keys(i);
-
- // FIXME -- we should be using case-insensitive comparisons.
-
- if (key == "TolX")
- {
- if (option_map.contains (key))
- {
- Cell c = option_map.contents (key);
-
- if (c.numel () == 1)
- {
- octave_value val = c(0);
-
- if (! val.is_empty ())
- {
- double dval = val.double_value ();
-
- if (! error_state)
- opts.set_tolerance (dval);
- else
- gripe_wrong_type_arg ("fsolve", val);
- }
- }
- else
- error ("fsolve: invalid value for %s option", key.c_str ());
- }
- }
- else
- {
- static std::set unimplemented_options
- = make_unimplemented_options ();
-
- if (unimplemented_options.find (key) != unimplemented_options.end ())
- {
- if (option_map.contains (key))
- {
- Cell c = option_map.contents (key);
-
- if (c.numel () == 1)
- {
- octave_value val = c(0);
-
- if (! val.is_empty ())
- warning_with_id ("Octave:fsolve-unimplemented-option",
- "fsolve: option `%s' not implemented",
- key.c_str ());
- }
- }
- }
- }
- }
-}
-
-#define FSOLVE_ABORT() \
- do \
- { \
- unwind_protect::run_frame ("Ffsolve"); \
- return retval; \
- } \
- while (0)
-
-#define FSOLVE_ABORT1(msg) \
- do \
- { \
- ::error ("fsolve: " msg); \
- FSOLVE_ABORT (); \
- } \
- while (0)
-
-#define FSOLVE_ABORT2(fmt, arg) \
- do \
- { \
- ::error ("fsolve: " fmt, arg); \
- FSOLVE_ABORT (); \
- } \
- while (0)
-
-DEFUN_DLD (fsolve, args, nargout,
- "-*- texinfo -*-\n\
address@hidden {Loadable Function} address@hidden, @var{fval}, @var{info}] =} fsolve (@var{fcn}, @var{x0}, @var{options})\n\
-Given @var{fcn}, the name of a function of the form @code{f (@var{x})}\n\
-and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\
-equations such that @code{f(@var{x}) == 0}.\n\
-\n\
-On return, @var{fval} contains the value of the function @var{fcn}\n\
-evaluated at @var{x}, and @var{info} may be one of the following values:\n\
-\n\
address@hidden @asis\n\
-\n\
address@hidden 1\n\
-Algorithm converged with relative error between two consecutive iterates\n\
-less than or equal to the specified tolerance (see @code{fsolve_options}).\n\
address@hidden 0\n\
-Iteration limit exceeded.\n\
address@hidden -1\n\
-Error in user-supplied function.\n\
address@hidden -2\n\
-Algorithm failed to converge.\n\
address@hidden table\n\
-\n\
-If @var{fcn} is a two-element string array, or a two element cell array\n\
-containing either the function name or inline or function handle. The\n\
-first element names the function @math{f} described above, and the second\n\
-element names a function of the form @code{j (@var{x})} to compute the\n\
-Jacobian matrix with elements\n\
address@hidden
-$$ J = {\\partial f_i \\over \\partial x_j} $$\n\
address@hidden tex\n\
address@hidden
-\n\
address@hidden
- df_i\n\
-jac(i,j) = ----\n\
- dx_j\n\
address@hidden example\n\
address@hidden ifinfo\n\
-\n\
-You can use the function @code{fsolve_options} to set optional\n\
-parameters for @code{fsolve}.\n\
-\n\
-If the optional argument @var{options} is provided, it is expected to\n\
-be a structure of the form returned by @code{optimset}. Options\n\
-specified in this structure override any set globally by\n\
address@hidden, fsolve_options}.\n\
address@hidden deftypefn")
-{
- octave_value_list retval;
-
- warned_fcn_imaginary = false;
- warned_jac_imaginary = false;
-
- unwind_protect::begin_frame ("Ffsolve");
-
- unwind_protect_int (call_depth);
- call_depth++;
-
- if (call_depth > 1)
- FSOLVE_ABORT1 ("invalid recursive call");
-
- int nargin = args.length ();
-
- if ((nargin == 2 || nargin == 3) && nargout < 4)
- {
- std::string fcn_name, fname, jac_name, jname;
- fsolve_fcn = 0;
- fsolve_jac = 0;
-
- octave_value f_arg = args(0);
-
- if (f_arg.is_cell ())
- {
- Cell c = f_arg.cell_value ();
- if (c.length() == 1)
- f_arg = c(0);
- else if (c.length() == 2)
- {
- if (c(0).is_function_handle () || c(0).is_inline_function ())
- fsolve_fcn = c(0).function_value ();
- else
- {
- fcn_name = unique_symbol_name ("__fsolve_fcn__");
- fname = "function y = ";
- fname.append (fcn_name);
- fname.append (" (x) y = ");
- fsolve_fcn = extract_function
- (c(0), "fsolve", fcn_name, fname, "; endfunction");
- }
-
- if (fsolve_fcn)
- {
- if (c(1).is_function_handle () || c(1).is_inline_function ())
- fsolve_jac = c(1).function_value ();
- else
- {
- jac_name = unique_symbol_name ("__fsolve_jac__");
- jname = "function y = ";
- jname.append (jac_name);
- jname.append (" (x) jac = ");
- fsolve_jac = extract_function
- (c(1), "fsolve", jac_name, jname, "; endfunction");
-
- if (!fsolve_jac)
- {
- if (fcn_name.length())
- clear_function (fcn_name);
- fsolve_fcn = 0;
- }
- }
- }
- }
- else
- FSOLVE_ABORT1 ("incorrect number of elements in cell array");
- }
-
- if (! fsolve_fcn && ! f_arg.is_cell())
- {
- if (f_arg.is_function_handle () || f_arg.is_inline_function ())
- fsolve_fcn = f_arg.function_value ();
- else
- {
- switch (f_arg.rows ())
- {
- case 1:
- do
- {
- fcn_name = unique_symbol_name ("__fsolve_fcn__");
- fname = "function y = ";
- fname.append (fcn_name);
- fname.append (" (x) y = ");
- fsolve_fcn = extract_function
- (f_arg, "fsolve", fcn_name, fname, "; endfunction");
- }
- while (0);
- break;
-
- case 2:
- {
- string_vector tmp = f_arg.all_strings ();
-
- if (! error_state)
- {
- fcn_name = unique_symbol_name ("__fsolve_fcn__");
- fname = "function y = ";
- fname.append (fcn_name);
- fname.append (" (x) y = ");
- fsolve_fcn = extract_function
- (tmp(0), "fsolve", fcn_name, fname, "; endfunction");
-
- if (fsolve_fcn)
- {
- jac_name = unique_symbol_name ("__fsolve_jac__");
- jname = "function y = ";
- jname.append (jac_name);
- jname.append (" (x) jac = ");
- fsolve_jac = extract_function
- (tmp(1), "fsolve", jac_name, jname,
- "; endfunction");
-
- if (!fsolve_jac)
- {
- if (fcn_name.length())
- clear_function (fcn_name);
- fsolve_fcn = 0;
- }
- }
- }
- }
- }
- }
- }
-
- if (error_state || ! fsolve_fcn)
- FSOLVE_ABORT ();
-
- NDArray xa = args(1).array_value ();
- x_dims = xa.dims ();
- ColumnVector x (xa);
-
- if (error_state)
- FSOLVE_ABORT1 ("expecting vector as second argument");
-
- if (nargin > 3)
- warning ("fsolve: ignoring extra arguments");
-
- if (nargout > 3)
- warning ("fsolve: can't compute path output yet");
-
- NLFunc nleqn_fcn (fsolve_user_function);
- if (fsolve_jac)
- nleqn_fcn.set_jacobian_function (fsolve_user_jacobian);
-
- NLEqn nleqn (x, nleqn_fcn);
-
- NLEqn_options local_fsolve_opts (fsolve_opts);
-
- if (nargin > 2)
- {
- Octave_map option_map = args(2).map_value ();
-
- if (! error_state)
- override_options (local_fsolve_opts, option_map);
- else
- error ("fsolve: expecting optimset structure as third argument");
- }
-
- if (! error_state)
- {
- nleqn.set_options (local_fsolve_opts);
-
- octave_idx_type info;
- ColumnVector soln = nleqn.solve (info);
-
- if (fcn_name.length())
- clear_function (fcn_name);
- if (jac_name.length())
- clear_function (jac_name);
-
- if (! error_state)
- {
- retval(2) = static_cast (hybrd_info_to_fsolve_info (info));
- retval(1) = nleqn.function_value ();
- retval(0) = NDArray (ArrayN (soln.reshape (x_dims)));
-
- if (! nleqn.solution_ok () && nargout < 2)
- {
- std::string msg = nleqn.error_message ();
- error ("fsolve: %s", msg.c_str ());
- }
- }
- }
- }
- else
- print_usage ();
-
- unwind_protect::run_frame ("Ffsolve");
-
- return retval;
-}
-
-/*
-%!function retval = f (p)
-%! x = p(1);
-%! y = p(2);
-%! z = p(3);
-%! retval = zeros (3, 1);
-%! retval(1) = sin(x) + y**2 + log(z) - 7;
-%! retval(2) = 3*x + 2**y -z**3 + 1;
-%! retval(3) = x + y + z - 5;
-%!test
-%! x_opt = [ 0.599054;
-%! 2.395931;
-%! 2.005014 ];
-%! tol = 1.0e-5;
-%! [x, fval, info] = fsolve ("f", [ 0.5; 2.0; 2.5 ]);
-%! info_bad = (info <= 0);
-%! solution_bad = sum (abs (x - x_opt) > tol);
-%! value_bad = sum (abs (fval) > tol);
-%! if (info_bad)
-%! printf_assert ("info bad\n");
-%! else
-%! printf_assert ("info good\n");
-%! endif
-%! if (solution_bad)
-%! printf_assert ("solution bad\n");
-%! else
-%! printf_assert ("solution good\n");
-%! endif
-%! if (value_bad)
-%! printf_assert ("value bad\n");
-%! else
-%! printf_assert ("value good\n");
-%! endif
-%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
-
-%!function retval = f (p)
-%! x = p(1);
-%! y = p(2);
-%! z = p(3);
-%! w = p(4);
-%! retval = zeros (4, 1);
-%! retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
-%! retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
-%! retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
-%! retval(4) = x^2 + 2*y^3 + z - w - 4;
-%!test
-%! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
-%! tol = 1.0e-5;
-%! [x, fval, info] = fsolve ("f", [-1, 1, 2, -1]);
-%! info_bad = (info <= 0);
-%! solution_bad = sum (abs (x - x_opt) > tol);
-%! value_bad = sum (abs (fval) > tol);
-%! if (info_bad)
-%! printf_assert ("info bad\n");
-%! else
-%! printf_assert ("info good\n");
-%! endif
-%! if (solution_bad)
-%! printf_assert ("solution bad\n");
-%! else
-%! printf_assert ("solution good\n");
-%! endif
-%! if (value_bad)
-%! printf_assert ("value bad\n");
-%! else
-%! printf_assert ("value good\n");
-%! endif
-%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
-
-%!test
-%! fsolve_options ("tolerance", eps);
-%! assert(fsolve_options ("tolerance") == eps);
-
-%!error fsolve_options ("foo", 1, 2);
-*/
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff --git a/src/Makefile.in b/src/Makefile.in
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -57,7 +57,7 @@
endif
endif
-OPT_BASE := $(addsuffix -opts, DASPK DASRT DASSL LSODE NLEqn Quad)
+OPT_BASE := $(addsuffix -opts, DASPK DASRT DASSL LSODE Quad)
OPT_HANDLERS := $(addsuffix .cc, $(OPT_BASE))
OPT_IN := $(addprefix ../liboctave/, $(addsuffix .in, $(OPT_BASE)))
OPT_INC := $(addprefix ../liboctave/, $(addsuffix .h, $(OPT_BASE)))
@@ -66,7 +66,7 @@
chol.cc ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
dasrt.cc dassl.cc det.cc dispatch.cc dlmread.cc dmperm.cc eig.cc \
expm.cc fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc \
- fltk_backend.cc fsolve.cc \
+ fltk_backend.cc \
gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
givens.cc hess.cc hex2num.cc inv.cc kron.cc lookup.cc lsode.cc \
lu.cc luinc.cc matrix_type.cc max.cc md5sum.cc pinv.cc qr.cc \