--- /dev/null
+c===========================================================================
+c
+c This file is part of TISEAN
+c
+c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+c
+c TISEAN is free software; you can redistribute it and/or modify
+c it under the terms of the GNU General Public License as published by
+c the Free Software Foundation; either version 2 of the License, or
+c (at your option) any later version.
+c
+c TISEAN is distributed in the hope that it will be useful,
+c but WITHOUT ANY WARRANTY; without even the implied warranty of
+c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+c GNU General Public License for more details.
+c
+c You should have received a copy of the GNU General Public License
+c along with TISEAN; if not, write to the Free Software
+c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+c
+c===========================================================================
+c constrained randomization
+c author T. Schreiber (1999)
+c===========================================================================
+ parameter(nx=100000,mx=20)
+ double precision time
+ dimension x(nx,mx), y(nx,mx), xx(nx,mx), icol(mx)
+ character*72 file, fout, comment
+ common nmax,cost,temp,cmin,rate,x
+ external rand
+ data wr/0.9/, nsur/1/
+ data iverb/15/
+
+ call whatido("constrained randomization",iverb)
+ call what_cost()
+ call what_cool()
+ call what_permute()
+ rr=rand(ican("I",0)/real(2**22))
+ nsur=min(999,ican("n",nsur))
+ nmax=ican("l",nx)
+ nexcl=ican("x",0)
+ mcmax=ican("m",0)
+ call columns(mc,mx,icol)
+ if(mcmax.eq.0) mcmax=max(1,mc)
+ wr=fcan("u",wr)
+ call opts_cost(mcmax)
+ call opts_cool()
+ call opts_permute()
+ isout=igetout(fout,iverb)
+
+ call nthstring(1,file)
+ call xreadfile(nmax,mcmax,nx,xx,nexcl,icol,file,iverb)
+ call cost_transform(nmax,mcmax,nx,xx)
+ if(file.eq."-") file="stdin"
+ if(isout.eq.1) call addsuff(fout,file,"_rnd")
+ if(nsur.gt.1) call suffix(fout,"_000")
+
+ do 10 isur=1,nsur
+ do 20 n=1,nmax
+ do 20 m=1,mcmax
+ 20 x(n,m)=xx(n,m)
+ rate=1
+ cost=r1mach(2)
+ cmin=cost
+ if(nsur.gt.1) write(fout(index(fout," ")-3:72),'(i3.3)') isur
+ temp=cool_init()
+ call cost_init()
+ call permute_init()
+ call cost_full(iv_vcost(iverb))
+ cmin=cost
+ time=0
+ 1 time=time+1.
+ call permute(n1,n2)
+ cmax=cost-temp*log(rand(0.0)) ! maximal acceptable cost
+ call cost_update(n1,n2,cmax,iaccept,iv_vmatch(iverb))
+ tnew=cool(iaccept,iend,iv_cool(iverb))
+ if(tnew.ne.temp.or.cost.lt.cmin*wr) then
+ cc=cost
+ call cost_full(iv_vcost(iverb))
+ if(iv_match(iverb).eq.1) write(istderr(),*)
+ . "cost function mismatch: ", abs((cc-cost)/cost)
+ endif
+ temp=tnew
+ if(cost.lt.cmin*wr) then
+ if(iv_cost(iverb).eq.1) write(istderr(),*)
+ . "after ",real(time)," steps at T=",temp," cost: ",cost
+ cmin=cost
+ call cost_inverse(nmax,mcmax,nx,x,y)
+ write(comment,'(8h# cost: ,g15.5)') cost
+ call xwritecfile(nmax,mcmax,nx,y,fout,iverb,comment)
+ endif
+ if(iend.ne.1) goto 1
+ write(comment,'(8h# cost: ,g15.5)') cost
+ call writecfile(nmax,mcmax,nx,y,fout,iverb,comment)
+ 10 continue
+ end
+
+ subroutine usage()
+c usage message
+
+ call whatineed("[-n# -u# -I# -o outfile -l# -x# -c# -V# -h]"//
+ . " [cost opt.] [cooling opt.] [permutation opt.] file")
+ call popt("n","number of surrogates (1)")
+ call popt("u","improvement factor before write (0.9)")
+ call popt("I","seed for random numbers (0)")
+ call popt("l","number of values to be read (all)")
+ call popt("x","number of values to be skipped (0)")
+ call popt("c","column to be read (1 or file,#)")
+ call pout("file_rnd(_nnn)")
+ call pall()
+ call ptext("Verbosity levels (add what you want):")
+ call ptext(" 1 = input/output" )
+ call ptext(" 2 = current cost if improved")
+ call ptext(" 4 = cost mismatch")
+ call ptext(" 8 = temperature etc. at cooling")
+ call ptext(" 16 = verbose cost if improved")
+ call ptext(" 32 = verbose cost mismatch")
+ write(istderr(),'()')
+ call usage_cost()
+ write(istderr(),'()')
+ call usage_cool()
+ write(istderr(),'()')
+ call usage_permute()
+ write(istderr(),'()')
+ stop
+ end
+