source/maxlik/src_FPy/dfp.f

   1       subroutine dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua,
   2      1               rest,iest,lint,n,*)
   3 c
   4 c this subroutine performs minimization of a certain objective function,
   5 c f, by the davidon-fletcher-powell method. diagonal and upper-diagonal
   6 c elements of the matrix v are stored in the vector v. linesearch is
   7 c included in the subroutine dirmin. the variables are converted into
   8 c their working form (vector x), and vice-versa by the subroutine convrt
   9 c calculation terminates normally if both tests on delta(x) and delta(g)
  10 c are satisfied. maximum number of cycles is maxit.
  11 c
  12 c ************** written by a.liwo, institute of chemistry
  13 c **************                    university of gdansk
  14 c **************                    ul.sobieskiego 18
  15 c **************                    80952 gdansk, poland
  16 c
  17 c ************** first r-32 version (fortran iv):   august 1986
  18 c ************** ibm pc adaptation  (fortran 77):   may 1988
  19 c ************** unix version       (fortran 77):   november 1990
  20 c
  21       implicit real*8 (a-h,o-z)
  22       include 'DIMENSIONS.ZSCOPT'
  23       parameter (maxdim=max_paropt)
  24       parameter (maxvdim=maxdim*(maxdim+1)/2)
  25       logical fail,init,qua,rest,lint
  26       common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
  27      1               vy(maxdim),g(maxdim)
  28       common /dirit/ nprint
  29
  30       print *,"DFP n",n
  31       call convrt(n,x,.true.)
  32       nvdim=n*(n+1)/2
  33       isave=nsave
  34       kkk=1
  35       if (.not.rest) goto 20
  36
  37       print *,'reading restart information from the restart file'
  38       rewind nfile
  39       nstep=nvdim/100
  40       nstep1=nstep
  41       if (nstep*100.lt.nvdim) nstep1=nstep+1
  42       read(nfile) kkk,f,dftau0
  43       do 222 i=1,n
  44   222 read (nfile) g(i),y(i),x(i),s(i)
  45       do 111 i=1,nstep1
  46       ns=(i-1)*100+1
  47       ne=i*100
  48       if (ne.gt.nvdim) ne=nvdim
  49   111 read (nfile) (v(j),j=ns,ne)
  50       endfile nfile
  51       print *,'restart information has been read in'
  52       if (kkk.gt.1) kk1=kkk
  53       goto 21
  54    20 f=fun(n,x,nfun)
  55       call grad(n,iest,g,x)
  56       do 201 i=1,n
  57   201 y(i)=g(i)
  58    21 n1=n+1
  59       nfun=0
  60       do 2 it=1,maxit
  61       if (kkk.gt.1) goto 17
  62       kk1=2
  63       ind=0
  64       do 3 i=1,n
  65       i1=i+1
  66       ind=ind+1
  67       v(ind)=1.
  68       if (i1.gt.n) goto 3
  69       do 4 j=i1,n
  70       ind=ind+1
  71     4 v(ind)=0
  72     3 continue
  73       do 5 i=1,n
  74     5 d(i)=-g(i)
  75       dftau0=dirder(n,g,d)
  76       call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua)
  77       if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n)
  78       if (fail) return 1
  79       if (.not.init) goto 1
  80       if (nprint.gt.-1) write (iout,1100) it,nfun
  81       if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
  82  1100 format (//1x,' * * * * * * * * no point remarkably lower in ene',
  83      1     'rgy can be found.'/1x,'calculation terminated in cycle',i3,
  84      2     ' after',i6,' function & gradient evaluations')
  85       call convrt(n,x,.false.)
  86 *     nprint=2
  87       call csave(lint)
  88       return
  89     1 if (n.eq.1) return
  90       xnorm=0.0d0
  91       gnorm=0.0d0
  92       xmax=0.0d0
  93       gmax=0.0d0
  94       do 163 i=1,n
  95       dxi=dabs(s(i))
  96       gi=dabs(y(i))
  97       if (dxi.gt.xmax) xmax=dxi
  98       if (gi.gt.gmax) gmax=gi
  99       xnorm=xnorm+dxi*dxi
 100   163 gnorm=gnorm+gi*gi
 101       xnorm=dsqrt(xnorm/n)
 102       gnorm=dsqrt(gnorm/n)
 103       kk=1
 104       if(nprint.gt.0)write(iout,1000)it,kk,xmax,xnorm,gmax,gnorm,f,nfun
 105       if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n)
 106  1070 format (/5x,' * * * * gradients'/100(1x,i5,5x,1pe12.4/))
 107       if (nfun.lt.isave) goto 23
 108       call convrt (n,x,.false.)
 109       call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
 110       call csave(lint)
 111       isave=isave+nsave
 112    23 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 17
 113       if (nprint.gt.-1) write (iout,1010) it,nfun
 114       if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
 115       return
 116    17 do 6 kk=kk1,n1
 117       do 7 k=1,n
 118       gk=y(k)
 119       y(k)=gk-g(k)
 120     7 g(k)=gk
 121       sy=0.
 122       do 8 k=1,n
 123     8 sy=sy+s(k)*y(k)
 124       if (dabs(sy).gt.tolrec) goto 162
 125       if (nprint.gt.-1) write (iout,1030)
 126  1030 format (//' * * * delta(x) is perpendicular to delta(g). ',
 127      1       'the program will try gradient direction * * *')
 128       goto 2
 129   162 do 9 i=1,n
 130       vyi=0.
 131       do 10 j=1,i
 132       ind=(2*n-j)*(j-1)/2+i
 133    10 vyi=vyi+v(ind)*y(j)
 134       i1=i+1
 135       if (i1.gt.n) goto 9
 136       do 11 j=i1,n
 137       ind=(2*n-i)*(i-1)/2+j
 138    11 vyi=vyi+v(ind)*y(j)
 139     9 vy(i)=vyi
 140       yvy=0.
 141       do 12 i=1,n
 142    12 yvy=yvy+y(i)*vy(i)
 143       ind=0
 144       do 13 i=1,n
 145       si=s(i)/sy
 146       vyi=vy(i)/yvy
 147       do 13 j=i,n
 148       ind=ind+1
 149       v(ind)=v(ind)+si*s(j)-vyi*vy(j)
 150    13 continue
 151       if (nprint.gt.4) call matout(nvdim,n)
 152       do 18 i=1,n
 153       if (v((2*n-i)*(i-1)/2+i).gt.0.) goto 18
 154       if (nprint.gt.-1) write (iout,1080) i
 155       goto 2
 156    18 continue
 157  1080 format (/1x,i3,'-th diagonal element of the matrix v is negative.'
 158      1      ,' restart with gradient direction follows.')
 159       do 14 i=1,n
 160       di=0.
 161       do 15 j=1,i
 162       ind=(2*n-j)*(j-1)/2+i
 163    15 di=di-v(ind)*g(j)
 164       i1=i+1
 165       if (i1.gt.n) goto 14
 166       do 16 j=i1,n
 167       ind=(2*n-i)*(i-1)/2+j
 168    16 di=di-v(ind)*g(j)
 169    14 d(i)=di
 170       dftau0=dirder(n,g,d)
 171       do 202 i=1,n
 172   202 y(i)=g(i)
 173       call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua)
 174       if (fail) return 1
 175       if (.not.init) goto 19
 176       if (nprint.gt.-1) write (iout,1110)
 177  1110 format (/1x,'energy from dirmin differs but slightly from the'
 178      1 ,' initial value. restart with gradient direction follows.')
 179       goto 2
 180    19 xnorm=0.0d0
 181       gnorm=0.0d0
 182       xmax=0.0d0
 183       gmax=0.0d0
 184       do 161 i=1,n
 185       dxi=dabs(s(i))
 186       gi=dabs(y(i))
 187       if (xmax.lt.dxi) xmax=dxi
 188       if (gmax.lt.gi) gmax=gi
 189       xnorm=xnorm+dxi*dxi
 190   161 gnorm=gnorm+gi*gi
 191       xnorm=dsqrt(xnorm/n)
 192       gnorm=dsqrt(gnorm/n)
 193       if (nprint.gt.0) write (iout,1000) it,kk,xmax,xnorm,gmax,gnorm,f,
 194      1                                   nfun
 195       if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n)
 196  1040 format (/1x,' * * * * parameters:'/100(1x,i5,5x,1pe14.6/))
 197       if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n)
 198  1000 format (/1x,'iteration:',i4,' direction:',i4,'  max. element of ',
 199      1 'delta(x): ',1pe14.6,' rms of delta(x): ',1pe14.6/' maximum ',
 200      2 'element of gradient: ',1pe14.6,' rms of gradient: ',1pe14.6,
 201      3 ' energy: ',1pe14.6,' fun. eval:',i6)
 202 c
 203       if (nfun.lt.isave) goto 22
 204       call convrt(n,x,.false.)
 205       call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
 206       call csave(lint)
 207       isave=isave+nsave
 208    22 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 6
 209       if (nprint.gt.-1) write (iout,1010) it,nfun
 210       if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
 211       if (nprint.gt.-1) write (iout,1111)
 212  1111 format (//11x,'matrix v')
 213       if (nprint.gt.-1) call matout(nvdim,n)
 214  1010 format (//1x,'* * * * * * convergence has been achieved in ',
 215      1 'cycle',i4,' after',i6,' function & gradient evaluations')
 216       call convrt(n,x,.false.)
 217 *     nprint=2
 218       call csave(lint)
 219       return
 220     6 continue
 221     2 kkk=1
 222 c if this point is reached, convergence hasn't been achieved.
 223       if (nprint.gt.-1) write (iout,1020) maxit,nfun
 224       if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
 225       if (nprint.gt.-1) write (iout,1111)
 226       if (nprint.gt.-1) call matout(nvdim,n)
 227       call convrt(n,x,.false.)
 228 *     nprint=2
 229       call csave(lint)
 230  1020 format (//1x,' * * * * * * * convergence has not been achieved',
 231      1 ' * * * * * * * '/1x,'number of cycles done:',i3,' number of ',
 232      2 'function & gradient evaluations:',i6)
 233       return
 234       end
 235 *******************************************************************************
 236       subroutine grad(n,iest,g,x)
 237       implicit real*8 (a-h,o-z)
 238       include 'DIMENSIONS.ZSCOPT'
 239       dimension g(max_paropt),x(max_paropt)
 240       integer uiparm,nf,iest
 241       external fdum
 242       double precision obg(max_paropt)
 243       call targetgrad(n,x,nf,g,uiparm,obg,fdum)
 244       return
 245       end
 246 *******************************************************************************
 247       subroutine dirmin(iout,iest,nfun,n,f,dftau0,d,g,x,s,fail,init,qua)
 248       implicit real*8 (a-h,o-z)
 249       include 'DIMENSIONS.ZSCOPT'
 250       parameter (maxdim=max_paropt)
 251       logical fail,init,qua,weiter
 252       common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov,
 253      1               xmin
 254       dimension x(n),s(n),d(n),g(n),z(maxdim),h(maxdim)
 255 *     write (2,*) 'dirmin has been called'
 256       fail=.false.
 257       init=.true.
 258       weiter=.false.
 259       amx=dabs(d(1))
 260       do 22 i=2,n
 261       pom=dabs(d(i))
 262       if (pom.gt.amx) amx=pom
 263    22 continue
 264       tau=dmin1(1.0d0,amov/amx)
 265       do 1 i=1,n
 266       s(i)=x(i)
 267       h(i)=g(i)
 268     1 z(i)=x(i)
 269 *     write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n)
 270       if (nprint.gt.1) write (iout,1020) dftau0
 271       do 2 ilin=1,maxlin
 272 *     write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n)
 273       if (weiter) goto 10
 274       do 3 id=1,maxd
 275       do 4 iexp=1,maxexp
 276       do 5 i=1,n
 277     5 x(i)=s(i)+tau*d(i)
 278       f1=fun(n,x,nfun)
 279       call grad(n,iest,g,x)
 280       dftau=dirder(n,g,d)
 281       if (nprint.gt.1) write (iout,1000) ilin,tau,iexp,f1,f,dftau
 282  1000 format (1x,'dirmin iter: ',i4,' tau=',1pe12.4,
 283      1 '  contr. step:',i4/' actual f value: ',1pe14.6,' old value: ',
 284      2 1pe14.6,' derivative: ',1pe14.6)
 285       if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n)
 286       if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n)
 287       if (dabs(f1-f).gt.tolqu.and.amx*tau.ge.xmin) goto 21
 288       if (f1.lt.f) goto 14
 289       do 24 i=1,n
 290       x(i)=s(i)
 291       g(i)=h(i)
 292    24 s(i)=x(i)-z(i)
 293       return
 294    21 if (f1.lt.f.or.dftau.gt.0.) goto 6
 295     4 tau=tau/expand
 296       if (nprint.gt.-1) write(iout,1010) f
 297       do 25 i=1,n
 298       g(i)=h(i)
 299    25 x(i)=s(i)
 300       if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
 301       fail=.true.
 302       return
 303  1010 format (//' * * * * * dirmin failed. least function value is: ',
 304      1        1pe14.6,' * * * * *')
 305  1080 format (/' * * * * parameters:'/100(1x,i5,5x,1pe14.6/))
 306  1070 format (/' * * * * gradients:'/100(1x,i5,5x,1pe14.6/))
 307  1020 format (/1x,'initial value of directional derivative: ',1pe14.6)
 308     6 if (dftau.gt.0.) goto 10
 309       init=.false.
 310       if (dabs(dftau/amx).le.dermin) goto 14
 311       dftau0=dftau
 312       do 9 i=1,n
 313       h(i)=g(i)
 314     9 s(i)=x(i)
 315       tau=tau*expand
 316       f=f1
 317     3 continue
 318       if (nprint.gt.-1) write (iout,1030)
 319  1030 format (//1x,' * * * * * error in dirmin. positive directional ',
 320      1  'derivative cannot be found * * * * *')
 321       if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
 322       fail=.true.
 323       return
 324    10 term=f1-f-dftau0*tau
 325       if (qua.or.dabs(dftau0/amx).lt.0.01.and.dabs(dftau/amx).lt.0.01)
 326      1 goto 30
 327       difd=tau*(dftau-dftau0)
 328       tt=tau*tau
 329       ttt=tt*tau
 330       aa=(tau*difd-2.*term)/ttt
 331       if (aa.eq.0.) goto 30
 332       bb=(3.*term-difd)/tt
 333       delta=4.*bb*bb-12.*aa*dftau0
 334       if (delta.lt.0.) goto 30
 335       delta=dsqrt(delta)
 336       tau1=(-2.*bb-delta)/(6.*aa)
 337       tau2=(-2.*bb+delta)/(6.*aa)
 338       if (nprint.gt.1) write (iout,1113) ilin,tau1,tau2
 339  1113 format (/1x,'dirmin iter.:',i3,' cubic interpolation',
 340      1 '  tau1=',1pe12.4,'  tau2= ',1pe12.4)
 341       if (tau1.gt.0.0d0.and.tau1.le.tau) goto 16
 342       if (tau2.le.0.0d0.or.tau2.gt.tau) goto 30
 343       tau1=tau2
 344       goto 16
 345    30 tau1=-dftau0*tau*tau/(2.*term)
 346       if (nprint.gt.1) write (iout,1040) ilin,tau1
 347       if (tau1.lt.tau.and.tau1.gt.0.) goto 16
 348       if (dabs(amx*tau1).lt.xmin) goto 145
 349       weiter=.false.
 350       tau=tau/expand
 351       if (f1.gt.f) goto 2
 352       do 17 i=1,n
 353       s(i)=x(i)
 354       h(i)=g(i)
 355    17 d(i)=-d(i)
 356       f=f1
 357       dftau0=-dftau
 358       goto 2
 359  1040 format (/1x,'dirmin iteration:',i4,' tau after quadr. interpo',
 360      1   'lation: ',1pe12.4)
 361    16 do 11 i=1,n
 362    11 x(i)=s(i)+tau1*d(i)
 363       if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n)
 364       if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n)
 365       f2=fun(n,x,nfun)
 366       call grad(n,iest,g,x)
 367       dftau1=dirder(n,g,d)
 368       dftn=dabs(dftau1/amx)
 369       if (nprint.gt.1) write (iout,1112) f2,dftau1
 370  1112 format (1x,'function value: ',1pe14.6,' derivative: ',1pe14.6)
 371       if (f2.lt.f.or.dftau1.gt.0.) goto 20
 372       if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).gt.xmin.and.dftn.ge.
 373      1 dermin) goto 27
 374       do 28 i=1,n
 375       x(i)=s(i)
 376       g(i)=h(i)
 377    28 s(i)=x(i)-z(i)
 378       goto 18
 379    27 weiter=.false.
 380       tau=tau1/expand
 381       goto 2
 382    20 weiter=.true.
 383       if (f2.lt.f) init=.false.
 384       if (f2.lt.f) goto 26
 385       if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).ge.xmin.and.dftn.ge.
 386      1 dermin) goto 31
 387       do 32 i=1,n
 388       x(i)=s(i)
 389       s(i)=x(i)-z(i)
 390    32 g(i)=h(i)
 391       goto 18
 392    31 f1=f2
 393       dftau=dftau1
 394       tau=tau1
 395       goto 2
 396    26 taudf=tau1*dftau1
 397       if (nprint.gt.1) write (iout,1111) taudf
 398  1111 format (/1x,'difference: ',1pe10.2)
 399       if (dabs(taudf).lt.tolqu.or.dftn.lt.dermin) goto 144
 400       if (dftau1.lt.0.) goto 13
 401       do 12 i=1,n
 402       d(i)=-d(i)
 403       h(i)=g(i)
 404    12 s(i)=x(i)
 405       dftau=-dftau0
 406       dftau0=-dftau1
 407       f1=f
 408       f=f2
 409       tau=tau1
 410       goto 2
 411    13 continue
 412       do 19 i=1,n
 413       h(i)=g(i)
 414    19 s(i)=x(i)
 415       dftau0=dftau1
 416       f=f2
 417       tau=tau-tau1
 418     2 continue
 419       if (nprint.gt.-1) write (iout,1050)
 420  1050 format (//1x,' * * * * nonconvergence in dirmin, calculation ',
 421      1 'terminated. * * * *')
 422       fail=.true.
 423       if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
 424       return
 425   144 f=f2
 426       goto 145
 427    14 f=f1
 428   145 do 15 i=1,n
 429    15 s(i)=x(i)-z(i)
 430    18 if (nprint.gt.1) write (iout,1060) ilin
 431  1060 format (/1x,'dirmin converged in',i3,' iterations')
 432       return
 433       end
 434 *******************************************************************************
 435       double precision function dirder(n,g,d)
 436       implicit real*8 (a-h,o-z)
 437       dimension d(n),g(n)
 438 c calculates the directional derivative. d is the direction vector,
 439 c g - gradient vector
 440       yy=0.0d0
 441       do 1 i=1,n
 442       yy=yy+d(i)*g(i)
 443     1 continue
 444       dirder=yy
 445       return
 446       end
 447 *******************************************************************************
 448       double precision function fun(n,x,nfun)
 449       implicit real*8 (a-h,o-z)
 450       include 'DIMENSIONS.ZSCOPT'
 451       parameter (maxdim=max_paropt)
 452       dimension x(max_paropt)
 453       external fdum
 454       integer uiparm
 455       double precision f,obf
 456       print *,"FUN: n=",n
 457       nfun=nfun+1
 458       call targetfun(n,x,nfun,f,uiparm,obf,fdum)
 459       fun=f
 460       return
 461       end
 462 *******************************************************************************
 463       subroutine convrt(n,x,dir)
 464       implicit real*8 (a-h,o-z)
 465       dimension x(n)
 466       return
 467       end
 468 *******************************************************************************
 469       subroutine dfpopt(inp,iout,ifirst,n,xopt)
 470       implicit real*8 (a-h,o-z)
 471 c
 472 c  this subroutine reads in the iteration parameters, calls
 473 c  the minimization procedures (davejr or dfp or both of them one by
 474 c  one) and writes the final coordinates to file ares.xyz.
 475 c
 476       include 'DIMENSIONS.ZSCOPT'
 477       parameter (maxdim=max_paropt)
 478       parameter (maxvdim=maxdim*(maxdim+1)/2)
 479       common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
 480      1               vy(maxdim),g(maxdim)
 481       logical conv,qua,rest,lint
 482       common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov,
 483      1               xmin
 484 c
 485       integer option,odfp,onr
 486       dimension xopt(maxdim)
 487       data odfp /4h dfp/,onr /4h  nr/, ncu /4h cub/,nrest /4hrest/
 488       print *,'ce2 has been called'
 489       print *,"DFPOPT: n",n
 490 c
 491 c      if (ifirst.eq.0) goto 1
 492 c      read (inp,1001) option,icub,irest
 493 c      rest=.false.
 494 c      if (irest.eq.nrest) rest=.true.
 495 c      qua=.true.
 496 c      if (icub.eq.ncu) qua=.false.
 497 c      read (inp,1010) maxit,itolx,itolg,maxlin,maxexp,maxd,itolq,itolr,
 498 c     1                iderm,nprint,ixmin,nsave,nfile,iest
 499 c      read (inp,1030) expand,amov
 500 c      if (nsave.eq.0) nsave=10000
 501 c      nfile=4
 502 c      if (itolx.eq.0) itolx=3
 503 c      if (itolg.eq.0) itolg=1
 504 c      if (itolr.eq.0) itolr=6
 505 c      if (iderm.eq.0) iderm=7
 506 c      if (itolq.eq.0) itolq=4
 507 c      if (ixmin.eq.0) ixmin=5
 508 c      if (maxit.eq.0) maxit=5
 509 c      if (maxlin.eq.0) maxlin=10
 510 c      if (maxd.eq.0) maxd=10
 511 c      if (maxexp.eq.0) maxexp=10
 512 c      if (expand.eq.0.) expand=2.0d0
 513 c      if (amov.eq.0.) amov=.25d0
 514 c      tolx=10.0d0**(-itolx)
 515 c      tolg=10.0d0**(-itolg)
 516 c      xmin=10.0d0**(-ixmin)
 517 c      tolqu=10.0d0**(-itolq)
 518 c      dermin=10.0d0**(-iderm)
 519 c      tolrec=10.0d0**(-itolr)
 520 c      if (option.eq.odfp) goto 4
 521       do i=1,n
 522         x(i)=xopt(i)
 523       enddo
 524       nprint=2
 525       nsave=100000
 526       nfile=89
 527       maxit=100
 528       maxlin=10
 529       maxd=10
 530       maxexp=10
 531       expand=2.0d0
 532       amov=.25d0
 533       tolx=1.0d-3
 534       tolg=1.0d-1
 535       xmin=1.0d-5
 536       tolqu=1.0d-3
 537       dermin=1.0d-7
 538       tolrec=1.0d-6
 539       iest=0
 540       qua=.false.
 541       rest=.false.
 542
 543     4 write (iout,2000) maxit,tolx,tolg,maxlin,maxexp,maxd,tolqu
 544       write (iout,2002) expand,amov,xmin,dermin,tolrec
 545       if (iest.gt.0) write (iout,2060)
 546       write (iout,2040) nsave,nfile
 547       if (rest) write (iout,2050) nfile
 548       if (qua) write (iout,2020)
 549       if (.not.qua) write (iout,2030)
 550     1 call dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua,rest,
 551      1         iest,lint,n,*2)
 552       do 300 i=1,n
 553   300   xopt(i)=x(i)
 554       return
 555     2 write (iout,2003)
 556       do 301 i=1,n
 557   301   xopt(i)=x(i)
 558  1001 format (3a4)
 559  1030 format (8f10.5)
 560  1010 format (14i5)
 561  2000 format (//'welcome to dfp routine'//'the following ',
 562      1 'options are used:'//
 563      2 'maximum number of cycles:',i4/'tolerance on coordinates:',
 564      3 1pe10.2/'tolerance on gradient:',1pe10.2/
 565      4 'maximum number of iterations in linesearch:',i4/
 566      5 'maximum number of unsuccesful function evaluations:',i4/
 567      6 'maximum number of steps in positive derivative search:',i4/
 568      7 'tolerance in linesearch:',1pe10.2)
 569  2002 format ('expansion coefficient in linesearch:',1pe10.2/
 570      9 'maximum movement of atoms:',1pe10.2/
 571      a 'minimum change of coordinates considered significant:',1pe10.2/
 572      c 'minimum value of directional derivative considered different ',
 573      d 'from zero:',1pe10.2/
 574      e 'lower bound of the scalar product (delta(x),delta(g)):',
 575      f 1pe10.2)
 576  2003 format (/'not quite succesful... perhaps you will have a '
 577      1         'better time'/)
 578 c
 579  2060 format (/'estimated partials used.')
 580  2020 format (/'quadratic interpolation will be used in line',
 581      1 ' minimization')
 582  2030 format (/'cubic or quadratic interpolation will be used ',
 583      1 'in line minimization')
 584  2040 format (/'coordinates, gradients, and matrix v will be saved',
 585      1 ' appr. every',i4,' energy and gradient evaluations on file',i3)
 586  2050 format (//'restart to be done using information saved ',
 587      1 'on file',i3/)
 588       return
 589       end
 590 *******************************************************************************
 591       subroutine csave(lint)
 592       logical lint
 593       return
 594       end
 595 *******************************************************************************
 596       subroutine psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
 597       implicit real*8 (a-h,o-z)
 598       include 'DIMENSIONS.ZSCOPT'
 599       parameter (maxdim=max_paropt)
 600       parameter (maxvdim=maxdim*(maxdim+1)/2)
 601       logical fail,init,qua,rest,lint
 602       common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
 603      1               vy(maxdim),g(maxdim)
 604       common /dirit/ nprint
 605       common /restart/ restfile
 606       character*16 restfile
 607       common inp,iout,ipn,iprint
 608       rewind nfile
 609       nstep=nvdim/100
 610       nstep1=nstep
 611       if (nstep*100.lt.nvdim) nstep1=nstep+1
 612       write (nfile) kk,f,dftau0
 613       do 222 i=1,n
 614   222 write (nfile) g(i),y(i),x(i),s(i)
 615       do 111 i=1,nstep1
 616       ns=(i-1)*100+1
 617       ne=i*100
 618       if (ne.gt.nvdim) ne=nvdim
 619   111 write (nfile) (v(j),j=ns,ne)
 620       endfile nfile
 621       write (*,1010) it,nfun,f,restfile
 622  1010 format (/1x,'cycle: ',i2,' energy evaluations: ',i5,' value:',
 623      1 f20.10,'.'/' information for restart saved on file ',a16)
 624       return
 625       end
 626 *******************************************************************************
 627       subroutine matout(ndim,n)
 628       implicit real*8 (a-h,o-z)
 629       include 'DIMENSIONS.ZSCOPT'
 630       parameter (maxdim=max_paropt)
 631       parameter (maxvdim=maxdim*(maxdim+1)/2)
 632       logical fail,init,qua,rest,lint
 633       common/bigmat/ a(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
 634      1               vy(maxdim),g(maxdim)
 635       character*6 line6 / '------' /
 636       character*12 line12 / '------------' /
 637       dimension b(8)
 638       data iout /2/
 639       do 1 i=1,n,6
 640       nlim=min0(i+5,n)
 641       write (iout,1000) (k,k=i,nlim)
 642       write (iout,1020) line6,(line12,k=i,nlim)
 643  1000 format (/7x,6(i7,5x))
 644  1020 format (a6,6a12)
 645       do 2 j=i,n
 646       jlim=min0(j,nlim)
 647       jlim1=jlim-i+1
 648       do 3 k=i,jlim
 649     3 b(k-i+1)=a(n*(k-1)-k*(k-1)/2+j)
 650       write (iout,1010) j,(b(k),k=1,jlim1)
 651     2 continue
 652     1 continue
 653  1010 format (i3,3x,6(1pe11.4,1x))
 654       return
 655       end
 656
 657