Mercurial > octave-nkf
changeset 8306:43795cf108d0
initial implementation of fsolve
remove old fsolve code
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sun, 28 Sep 2008 21:09:35 +0200 |
| parents | 368b504777a8 |
| children | ec969f3b8955 |
| files | configure.in libcruft/Makefile.in libcruft/minpack/Makefile.in libcruft/minpack/dogleg.f libcruft/minpack/dpmpar.f libcruft/minpack/enorm.f libcruft/minpack/fdjac1.f libcruft/minpack/hybrd.f libcruft/minpack/hybrd1.f libcruft/minpack/hybrj.f libcruft/minpack/hybrj1.f libcruft/minpack/qform.f libcruft/minpack/qrfac.f libcruft/minpack/r1mpyq.f libcruft/minpack/r1updt.f liboctave/Makefile.in liboctave/NLConst.h liboctave/NLEqn-opts.in liboctave/NLEqn.cc liboctave/NLEqn.h liboctave/NLFunc.h liboctave/NLP.h scripts/ChangeLog scripts/optimization/Makefile.in scripts/optimization/__dogleg__.m scripts/optimization/__fdjac__.m scripts/optimization/fsolve.m src/DLD-FUNCTIONS/fsolve.cc src/Makefile.in |
| diffstat | 29 files changed, 413 insertions(+), 3565 deletions(-) [+] |
line wrap: on
line diff
--- a/configure.in Fri Oct 31 08:06:45 2008 +0100 +++ b/configure.in Sun Sep 28 21:09:35 2008 +0200 @@ -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
--- a/libcruft/Makefile.in Fri Oct 31 08:06:45 2008 +0100 +++ b/libcruft/Makefile.in Sun Sep 28 21:09:35 2008 +0200 @@ -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
--- a/libcruft/minpack/Makefile.in Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -# <http://www.gnu.org/licenses/>. - -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
--- a/libcruft/minpack/dogleg.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/dpmpar.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 (chpf127@emx.utexas.edu) 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
--- a/libcruft/minpack/enorm.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/fdjac1.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/hybrd.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/hybrd1.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/hybrj.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/hybrj1.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/qform.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/qrfac.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/r1mpyq.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/libcruft/minpack/r1updt.f Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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
--- a/liboctave/Makefile.in Fri Oct 31 08:06:45 2008 +0100 +++ b/liboctave/Makefile.in Sun Sep 28 21:09:35 2008 +0200 @@ -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 \
--- a/liboctave/NLConst.h Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#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: *** -*/
--- a/liboctave/NLEqn-opts.in Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -# <http://www.gnu.org/licenses/>. - -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, -@code{fsolve_options} 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
--- a/liboctave/NLEqn.cc Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#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<double> fvec (n); - double *pfvec = fvec.fortran_vec (); - - octave_idx_type lwa = (n*(n+13))/2; - Array<double> wa (lwa); - double *pwa = wa.fortran_vec (); - - Array<double> 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<double> fvec (n); - double *pfvec = fvec.fortran_vec (); - - octave_idx_type lwa = (n*(3*n+13))/2; - Array<double> 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: *** -*/
--- a/liboctave/NLEqn.h Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#if !defined (octave_NLEqn_h) -#define octave_NLEqn_h 1 - -#include <cfloat> - -#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: *** -*/
--- a/liboctave/NLFunc.h Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#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: *** -*/
--- a/liboctave/NLP.h Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#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: *** -*/
--- a/scripts/ChangeLog Fri Oct 31 08:06:45 2008 +0100 +++ b/scripts/ChangeLog Sun Sep 28 21:09:35 2008 +0200 @@ -1,3 +1,10 @@ +2008-09-28 Jaroslav Hajek <highegg@gmail.com> + + * 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. + 2008-09-28 Jaroslav Hajek <highegg@gmail.com> * optimization/fzero.m: Replace tabs by spaces.
--- a/scripts/optimization/Makefile.in Fri Oct 31 08:06:45 2008 +0100 +++ b/scripts/optimization/Makefile.in Sun Sep 28 21:09:35 2008 +0200 @@ -35,6 +35,9 @@ SOURCES = \ __fsolve_defopts__.m \ fzero.m \ + __fdjac__.m \ + __dogleg__.m \ + fsolve.m \ glpk.m \ glpkmex.m \ lsqnonneg.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/__dogleg__.m Sun Sep 28 21:09:35 2008 +0200 @@ -0,0 +1,60 @@ +## Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com> +## +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {@var{x}} = __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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/__fdjac__.m Sun Sep 28 21:09:35 2008 +0200 @@ -0,0 +1,38 @@ +## Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com> +## +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{fval}, @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 + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/fsolve.m Sun Sep 28 21:09:35 2008 +0200 @@ -0,0 +1,297 @@ +## 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 +## <http://www.gnu.org/licenses/>. +## +## Author: Jaroslav Hajek <highegg@gmail.com> + +# -*- texinfo -*- +# @deftypefn{Function File} {} fsolve(@var{fcn}, @var{x0}, @var{options}) +# @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{}) +# Solves a system of nonlinear equations defined by the function @var{fcn}. +# @var{fcn} should accepts a vector (array) defining the unknown variables, +# and return a vector of left-hand sides of the equations. Right-hand sides +# are defined to be zeros. +# In other words, this function attempts to determine a vector @var{X} such +# that @code{@var{fcn}(@var{X})} gives (approximately) all zeros. +# @var{x0} determines a starting guess. The shape of @var{x0} is preserved +# in all calls to @var{fcn}, but otherwise it is treated as a column vector. +# @var{options} is a structure specifying additional options. Currently, fsolve +# recognizes these options: FunValCheck, OutputFcn, TolX, TolFun, MaxIter, +# MaxFunEvals and Jacobian. +# +# If Jacobian is 'on', it specifies that @var{fcn}, called with 2 output arguments, +# also returns the Jacobian matrix of right-hand sides at the requested point. +# TolX specifies the termination tolerance in the unknown variables, while +# TolFun is a tolerance for equations. Default is @code{1e1*eps} +# for TolX and @code{1e2*eps} for TolFun. +# For description of the other options, see @code{optimset}. +# +# On return, @var{fval} contains the value of the function @var{fcn} +# evaluated at @var{x}, and @var{info} may be one of the following values: +# +# @table @asis +# @item 1 +# Converged to a solution point. Relative residual error is less than specified +# by TolFun. +# @item 2 +# Last relative step size was less that TolX. +# @item 3 +# Last relative decrease in residual was less than TolF. +# @item 0 +# Iteration limit exceeded. +# @item -3 +# The trust region radius became excessively small. +# @end table +# +# Note: If you only have a single nonlinear equation of one variable, using +# @code{fzero} is usually a much better idea. +# @seealso{fzero,optimset} +# @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; + # FIXME: 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);
--- a/src/DLD-FUNCTIONS/fsolve.cc Fri Oct 31 08:06:45 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <algorithm> -#include <set> -#include <string> - -#include <iomanip> -#include <iostream> -#include <sstream> - -#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<double> (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<double> (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<std::string> -make_unimplemented_options (void) -{ - static bool initialized = false; - - static std::set<std::string> 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<std::string> 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\ -@deftypefn {Loadable Function} {[@var{x}, @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\ -@table @asis\n\ -\n\ -@item 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\ -@item 0\n\ -Iteration limit exceeded.\n\ -@item -1\n\ -Error in user-supplied function.\n\ -@item -2\n\ -Algorithm failed to converge.\n\ -@end 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\ -@tex\n\ -$$ J = {\\partial f_i \\over \\partial x_j} $$\n\ -@end tex\n\ -@ifinfo\n\ -\n\ -@example\n\ - df_i\n\ -jac(i,j) = ----\n\ - dx_j\n\ -@end example\n\ -@end 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\ -@code{optimset, fsolve_options}.\n\ -@end 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<double> (hybrd_info_to_fsolve_info (info)); - retval(1) = nleqn.function_value (); - retval(0) = NDArray (ArrayN<double> (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 <Invalid call to fsolve_options.*> fsolve_options ("foo", 1, 2); -*/ - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/Makefile.in Fri Oct 31 08:06:45 2008 +0100 +++ b/src/Makefile.in Sun Sep 28 21:09:35 2008 +0200 @@ -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 \