http://server.ccl.net/cca/software/SOURCES/FORTRAN/fit-any-function/fit/simplex2.shtml
CCL simplex2.f
```      dimension rea(40)
dimension parvert(maxpar*(maxpar+1)),refl(maxpar),ssd(maxpar+1),
+           expa(maxpar),cont(maxpar),var((maxvar+1)*maxpoints),
+           par(maxpar),x(maxvar+1),center(maxpar)
c
sdconv  = .false.
parconv = .true.
npar    = maxpar
nvar    = maxvar
npoints = 0
c
c     initialize the random number generator
c
seed    = 72767
xxx     = rand(seed)
c
c     testing to see if limits are exceeded.
c
if ((npar.gt.maxpar).or.(nvar.gt.(maxvar+1)).or.(npoints.gt.
+        maxpoints)) then
write (*,1001) npar,maxpar,nvar,maxvar,npoints,maxpoints
stop
end if
c
c     reading the initial values for the parameters :
c
do i=1,npar
parvert(i)=rea(i)
end do
c
c     reading all the values for the variables :
c
ymean=0.
ivar=0
if (nrea.ge.(nvar+1)) then
npoints=npoints+1
do j=1,(nvar+1)
ivar=ivar+1
var(ivar)=rea(j)
if (j.eq.(nvar+1)) ymean=ymean+var(nvar+1)
end do
endif
goto 10
20 continue
ymean=ymean/npoints
c
c     generating the initial simplex :
c
c             first two vertex at +10% of first par and +/- 10% of all
c                            others.
c
lastpar=npar
do i1=1,npar
i=i1*npar
do j=1,npar
parvert(i+j)=parvert(j)*(0.8+0.4*rand(0))
enddo
enddo
lastpar=npar*(npar+1)
c
c     now that we have the initial vertex ready calculate the SSD for
c     each. first loop over all vertex :
c
ncycles=-1
ibest=1
iworst=1
100 ncycles=ncycles+1
do i=1,(npar+1)
ssd(i)=0.
ifirstpar=(i-1)*npar
c
c             transfer all parameters from PARVERT to PAR :
c
do j=1,npar
par(j)=parvert(ifirstpar+j)
end do
c
c             now call stddev for each vertex
c
call stddev(par,var,ssd(i),x)
end do
c
c        finding the best and worst vertex :
c
do i=1,npar+1
if (ssd(i).lt.ssd(ibest)) then
ibest=i
else
if (ssd(i).gt.ssd(iworst)) iworst=i
end if
end do
ifirstbest=(ibest-1)*npar
ifirstworst=(iworst-1)*npar
c
c     at this point IWORST contains the index of the worst vertex,
c     and IBEST that of the best.
c     test for convergence in any cycle other than the first.
c
if (ncycles.ne.0) then
c
c        first condition is in SSD of best, RMS.LE.0.01% of YMEAN
c
cond=abs(ymean/10000.)
rms=sqrt(ssd(ibest)/npoints)
if (rms.le.cond) sdconv=.true.
c
c        also stop if best and worst vertex are within 0.001% of each other
c
parconv=.true.
do i=1,npar
i1=i*npar
do j=1,npar
if (parconv) then
test=parvert(i1+j)-parvert(j)
test=abs(100000.*test/parvert(j))
if (test.gt.1.0) parconv=.false.
endif
end do
end do
c
if ((parconv).or.(sdconv)) goto 999
endif
c
c     reflect the worst vertex and obtain its SSD :
c
ifirstworst=(iworst-1)*npar
ifirstbest=(ibest-1)*npar
do j=1,npar
center(j)=0.
do i=1,(npar+1)
ifirstpar=(i-1)*npar
if (i.ne.iworst) center(j)=center(j)+
+                           parvert(ifirstpar+j)
end do
center(j)=parvert(ifirstworst+j)-(center(j)/(npar))
refl(j)=parvert(ifirstworst+j)-(2.0*center(j))
end do
call stddev(refl,var,sdr,x)
c
c     compare the reflexed vertex with the best :
c
if (sdr.le.ssd(ibest)) then
c
c            test expanded
c
do i=1,npar
expa(i)=refl(i)-center(i)
end do
call stddev(expa,var,sde,x)
if (sde.le.sdr) then
c
c                here if expanded better than reflected better than best
c
ssd(iworst)=sde
do i=1,npar
parvert(ifirstworst+i)=expa(i)
end do
goto 100
else
c
c                here if expanded not better than reflected but reflected
c                better than best
c
ssd(iworst)=sdr
do i=1,npar
parvert(ifirstworst+i)=refl(i)
end do
goto 100
end if
endif
c
c     if reflected not better than best :
c
if (sdr.lt.ssd(iworst)) then
c
c          if better than worst :
c
ssd(iworst)=sdr
do i=1,npar
parvert(ifirstworst+i)=refl(i)
end do
goto 100
else
c
c          if worse than worst try contracted
c
do i=1,npar
cont(i)=parvert(ifirstworst+i)-(center(i)/2.)
end do
call stddev(cont,var,sdc,x)
if (sdc.le.ssd(iworst)) then
c
c                 here if contracted better than worst
c
ssd(iworst)=sdc
do i=1,npar
parvert(ifirstworst+i)=cont(i)
end do
goto 100
else
c
c                 no improvement, try shrinking
c
do i=1,(npar+1)
if (i.ne.ibest) then
ifp=(i-1)*npar
do j=1,npar
temp=parvert(ifp+j)
best=parvert(ifirstbest+j)
parvert(ifp+j)=temp-(0.5*(temp-best))
end do
end if
end do
goto 100
end if
end if
999 write (*,1003) ncycles
if (sdconv) write (*,1004)
if (parconv) write (*,1005)
write (*,'(1h )')
do j=1,npar
write (*,1006) j, parvert(ifirstbest+j)
end do
write (*,1007) rms,ssd(ibest)/npoints
c
c     to build the format for table :
c
nvp3=nvar+3
if (nvp3.gt.9) write (mf,'(i2)') nvp3
if (nvp3.le.9) write (mf,'(i1)') nvp3
wf='(1x,'//mf//'(d9.3,1x))'
c
c     to build the heading :
c
do i=1,132
hdng(i:i)=' '
end do
ihdng=2
do i=1,nvar
hdng(ihdng+3:ihdng+4)='x('
write (mf,'(i2)') i
hdng(ihdng+5:ihdng+6)=mf(1:2)
hdng(ihdng+7:ihdng+7)=')'
ihdng=ihdng+10
end do
hdng(ihdng+5:ihdng+5)='y'
ihdng=ihdng+10
hdng(ihdng+3:ihdng+7)='ycalc'
ihdng=ihdng+10
hdng(ihdng+4:ihdng+6)='dif'
ihdng=ihdng+6
write (*,'(a)') hdng(1:ihdng)
write (*,'(1h )')
c
do i=1,npar
par(i)=parvert(ifirstbest+i)
end do
do i=1,npoints
ifirstpoint=(i-1)*(nvar+1)
do j=1,(nvar+1)
x(j)=var(ifirstpoint+j)
end do
yc=ycalc(par,x)
dy=x(nvar+1)-yc
write (*,fmt=wf) (x(j),j=1,(nvar+1)),yc,dy
end do
end if
stop
1001 format (/,5x,'npar=',i3,' maxpar=',i3,';  nvar=',i3,' maxvar=',
+        i3,';   npoints=',i3,' maxpoints=',i3)
1002 format(a256)
1003 format (/,5x,'Convergence reached after',I4,' cycles',/)
1004 format (10x,'RMS less than 0.1% of mean Y')
1005 format (10x,'All the Vertex are within 0.01%')
1006 format (5x,'a(',i2,') = ',D13.5)
1007 format (/,10x,'RMS = ',D13.5,'  Std Dev = ',D13.5,///)
end
c---------------------------------------------------------------------
subroutine stddev(par,var,sd,x)
c---------------------------------------------------------------------
c
c     It will loop over all points and calculate the summation of the
c     standard deviation for a given set of parameters.
c
implicit double precision(a-h,o-z)
common /nums/npar,nvar,npoints
dimension par(2),var(2),x(2)
sd=0.
c
c     looping over all points :
c
do i=1,npoints
ifirstv=(i-1)*(nvar+1)
c
c        transfering the variables from VAR to X :
c
do j=1,(nvar+1)
x(j)=var(ifirstv+j)
end do
c
c        now calling the YCALC function to evaluate Y :
c
yc=ycalc(par,x)
dy=yc-x(nvar+1)
sd=sd+(dy*dy)
end do
return
end
c---------------------------------------------------------------------
double precision function ycalc(a,x)
c---------------------------------------------------------------------
c
c     Evaluates YCALC for a given set of parameters and a set of
c     variables.
c
implicit double precision(a-h,o-z)
dimension a(*),x(*)
c
c     the definition of Y should be inserted here :
c
```
 Modified: Fri Dec 11 17:00:00 1992 GMT Page accessed 3491 times since Sat Apr 17 22:01:38 1999 GMT