      SUBROUTINE SIMPLX(F,N,XOUT,FOUT,NEV,ACC,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C  ROUTINE SEARCHES FOR A MINIMUM IN N-DIM SPACE OF PARAMETERS
C
C  --INPUT--
C   F   - FUNCTION, IN A CALLING PROGRAM F(X), WHERE X IS AN ARRAY OF N
C         PARAMETERS, MUST APPEAR IN AN EXTERNAL
C   N   - NUMBER OF PARAMETERS
C   XOUT- GUESS FOR A MINIMUM
C   NEV - LIMIT ON FUNCTION EVALUATIONS
C   ACC - REQUIRED ACCURACY, .GT.0. ABSOLUTE, .LT.0. RELATIVE
C   WK  - WORK ARRAY, MUST BE DIMENSIONED WK(N*N+4*N+1), AT AN ENTRY
C         WK(1) IS A SCALE PARAMETER, IF WK(1).LT.0. THEN WK(I)
C         I=1,...,N  CONTAINS SCALE PARAMETERS, ATTENTION - CANNOT BE TO
C         SMALL
C
C  --OUTPUT--
C   XOUT- LOCATION OF MINIMUM
C   FOUT- VALUE OF MINIMUM
C   NEV - ACTUAL NUMBER OF EVALUTIONS
C   ACC - ESTIMATED ACCURACY
C
      EXTERNAL F
      DIMENSION XOUT(*),WK(*)
C
      CALL SIMPLE(F,N,N+1,XOUT,FOUT,NEV,ACC,WK(N+2),WK,WK(N*N+2*N+2),
     S  WK(N*N+3*N+2))
C
      RETURN
      END


      SUBROUTINE SIMPLE(F,N,M,XOUT,FOUT,NEV,ACC,X,FX,XR,XRR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F
      LOGICAL L
      DIMENSION X(N,M)
      DIMENSION XOUT(N),FX(M),XR(N),XRR(N)
C
      L=ACC.LT.0.D0
      PN=1.D0/N
      N1=N+1
      PN1=1.D0/N1
C
C  INITIAL POINTS
      IF(FX(1).EQ.0.D0)FX(1)=1.D0
C
      IF(FX(1).GT.0.D0)GOTO 1
C
      FX(1)=-FX(1)
      GOTO 3
C
 1    DO 2 I=1,N
 2    FX(I)=FX(1)
C
 3    DO 4 I=1,N
      FR=XOUT(I) !-PN1*FX(I)
      DO 4 J=1,N1
 4    X(I,J)=FR
C
      DO 5 J=1,N
 5    X(J,J+1)=X(J,J+1)+FX(J)
C
C
C  VALUES AT INITIAL POINTS
      DO 10 J=1,N1
 10   FX(J)=F(X(1,J))
      IEV=N1+1
C
C  THE SEARCH SETS IN
 20   JL=1
      JH=N1
      DO 30 J=1,N1
      IF(FX(J).GT.FX(JH))JH=J
      IF(FX(J).LT.FX(JL))JL=J
 30   CONTINUE
C
C  TESTING FOR EXIT
      FRR=0.D0
      DO 32 J=1,N1
 32   FRR=FRR+(FX(J)-FX(JL))**2
      FRR=SQRT(PN*FRR)
      FR=ACC
      IF(L)FR=ABS(FX(JL)*ACC)
      IF(FRR.LE.FR.OR.IEV.GE.NEV)GOTO 190
C
C  EVALUATING A MEAN
 38   DO 50 I=1,N
      FR=0.D0
      DO 40 J=1,N1
 40   FR=FR+X(I,J)
 50   XOUT(I)=(FR-X(I,JH))*PN
C
C  REFLECTED POINT
      DO 60 I=1,N
 60   XR(I)=2.D0*XOUT(I)-X(I,JH)
      FR=F(XR)
      IEV=IEV+1
C
C  TESTING FOR EXPANSION
      IF(FR.LT.FX(JL))GOTO 130
C
C  TESTING FOR REFLECTED POINT
      DO 70 J=1,N1
      IF(FR.LE.FX(J).AND.J.NE.JH)GOTO 170
 70   CONTINUE
C
C  TESTING FOR A DIRECT REDUCED ATTEMPT
      IF(FR.GT.FX(JH))GOTO 90
C
C  INSERTION OF THE REFLECTED POINT
      FX(JH)=FR
      DO 80 I=1,N
 80   X(I,JH)=XR(I)
C
C  REDUCED ATTEMPT
 90   DO 100 I=1,N
 100  XRR(I)=.5D0*(X(I,JH)+XOUT(I))
      FRR=F(XRR)
      IEV=IEV+1
C
C  TESTING FOR THE EXPANDED POINT
      IF(FRR.LE.FX(JH))GOTO 150
C
C  REDUCING SIMPLEX
      DO 120 J=1,N1
      IF(J.EQ.JL)GOTO 120
      DO 110 I=1,N
 110  X(I,J)=.5D0*(X(I,J)+X(I,JL))
      FX(J)=F(X(1,J))
 120  CONTINUE
      IEV=IEV+N
      GOTO 20
C
C  EXPANDING IN THE DIRECTION
 130  DO 140 I=1,N
 140  XRR(I)=4.D0*XOUT(I)-3.D0*X(I,JH)
      FRR=F(XRR)
      IEV=IEV+1
C
C  TESTING FOR THE REFLECTED POINT
      IF(FRR.GE.FX(JL))GOTO 170
C
C  ACCEPTING EXPANDED POINT
 150  FX(JH)=FRR
      DO 160 I=1,N
 160  X(I,JH)=XRR(I)
      GOTO 20
C
C  ACCEPTING REFLECTED POINT
 170  FX(JH)=FR
      DO 180 I=1,N
 180  X(I,JH)=XR(I)
      GOTO 20
C
C  EXIT
 190  NEV=IEV
C
      DO 210 I=1,N
      FR=0.D0
      DO 200 J=1,N1
 200  FR=FR+X(I,J)
 210  XOUT(I)=FR*PN1
      FOUT=F(XOUT)
C
      IF(FOUT.LE.FX(JL))GOTO 230
C
      FOUT=FX(JL)
      DO 220 I=1,N
 220  XOUT(I)=X(I,JL)
C
 230  ACC=FRR
      IF(L.AND.FOUT.NE.0.)ACC=ACC/ABS(FOUT)
C
      RETURN
      END





