plnsrch.f :::::::::::::: c The following code is a slightly modified version of Stephen E. Wright's c original code. subroutine lnsrch(type,w1,we1,xy1,xye1,w2,we2,xy2,xye2, # av,am,bv,bm,aev,aem,bev,bem,dm, # cm,cv,cem,cev, # zlo,zhi,zelo,zehi, # lambda,z,ze,fval,notify, # wdim,wedim,zdim,zedim,xydim,ntime, # wk,wklen) c c INPUT: integer zdim,wdim,zedim,wedim,xydim,ntime double precision w1(wdim,ntime),w2(wdim,ntime) double precision we1(wedim),we2(wedim) double precision xy1(xydim,ntime),xy2(xydim,ntime) double precision xye1(xydim),xye2(xydim) double precision av(wdim,ntime),am(wdim) double precision aev(wedim),aem(wedim) double precision bv(zdim,ntime),bm(zdim) double precision bev(zedim),bem(zedim) double precision cv(zdim),cm(zdim*xydim) double precision cev(zedim),cem(zedim*xydim) double precision dm(zdim*wdim) double precision zlo(zdim),zhi(zdim),zelo(zedim),zehi(zedim) character*(*) type c c OUTPUT: double precision lambda,z(zdim,ntime),ze(zedim),fval integer notify c c WORKSPACE PASSED THRU SUBROUTINE CALL: integer wklen double precision wk(wklen) c c LOCAL VARIABLES: double precision ct,atv,atm integer zlft,zrt,zmid,brk1,brk2,slope,zelft,zert,zemid, : ebrk1,ebrk2,eslope,btil,dtil,betil,detil c c c*********************************************************************** c c If type='min' then this routine solves the problem c c minimax sum [ a *w + (1/2)w *A w + b *z - (1/2)z *B z - z *Dw ] c t=1,N t t t t t t t t t t c c + a *w + (1/2)w *A w + b *z - (1/2)z *B z c e e e e e e e e e e c c - c e e c c subject to (w,w ) in [(w1,we1),(w2,we2)] c e c and zlo <= z <= zhi , zelo<= z <= zehi , c e c c T T c where = [ sum (C z + c)*x ] + (C z + c )*x c e e t=1,N t t-1 e e e N c c and x = Ex + Fw + b for t=1,...,N , c t t-1 t c c x = x = F w + b . c 0 e e e e c c (Note: w,w ,x,x are given, so E,E , and F are not needed.) c e e e c c c First LNPREP reformulates the problem as c _ _ _ c minimax c + a*lambda + (1/2)lambda*A*lambda c _ _ c + [ sum b*z - (1/2)z *Bz + lambda*(d*z ) ] c t=1,N t t t t c _ _ c + b *z - (1/2)z *B z + lambda*(d *z ) c e e e e e e e c c subject to lambda in [0,1] and zlo <= z <= zhi , zelo<= z <= zehi. c e c c c Next, SEARCH finds a saddle-pair (lambda,z). The solution of the original c problem is then ( lambda*(w1,we1) + (1-lambda)*(w2,we2) ; z) c with saddlevalue "fval". c c If type='max' then the routine solves the problem c c maximin sum [ a*w - (1/2)w *Aw + b*z + (1/2)z *Bz - z *Dw ] c t=1,N t t t t t t t t c c + a *w - (1/2)w *A w + b *z + (1/2)z *B z c e e e e e e e e e e c c - c e e c c subject to (w,w ) in [(w1,we1),(w2,we2)] c e c and zlo <= z <= zhi , zelo<= z <= zehi , c e c c where = [ sum (C z + c)*x ] + (C z + c )*x c e e t=1,N t t+1 e e e 1 c c T T c and x = E x + F w + b for t=1,...,N , c t t+1 t c T c x = x = F w + b . c N+1 e e e e c c by a similar reformulation. c c If type does not equal either 'max' or 'min' then 'min' is assumed. c - - - - - - - - - - - - - - - - - - - - - - - - c zlft = 1 zrt = zlft + zdim*ntime zmid = zrt + zdim*ntime brk1 = zmid + zdim*ntime brk2 = brk1 + zdim*ntime slope = brk2 + zdim*ntime btil = slope + zdim*ntime dtil = btil + zdim*ntime zelft = dtil + zdim*ntime zert = zelft + zedim zemid = zert + zedim ebrk1 = zemid + zedim ebrk2 = ebrk1 + zedim eslope = ebrk2 + zedim betil = eslope + zedim detil = betil + zedim c if ( (detil+zedim) .le. wklen ) then if ( type .ne. 'max' ) then print*, '* In the primal line search ' else print*, '* In the dual line search ' endif call lnprep(w1,we1,xy1,xye1,w2,we2,xy2,xye2,av,am,bv,bm, # aev,aem,bev,bem,dm,cm,cv,cem,cev, # zlo,zhi,zelo,zehi,ct,atv,atm, # wk(zlft),wk(zrt),wk(zmid),wk(brk1),wk(brk2), # wk(slope),wk(btil),wk(dtil), # wk(zelft),wk(zert),wk(zemid),wk(ebrk1), # wk(ebrk2),wk(eslope),wk(betil),wk(detil), # wdim,zdim,wedim,zedim,xydim,ntime,type) call search(ct,atv,atm,bem,bm,wk(btil),wk(dtil),wk(betil), # wk(detil),wk(zlft),wk(zrt),wk(zmid),wk(brk1), # wk(brk2),wk(slope),wk(zelft),wk(zert),wk(zemid), # wk(ebrk1),wk(ebrk2),wk(eslope), # lambda,z,ze,fval,zdim,zedim,ntime) notify = 0 if ( type .eq. 'max' ) fval = -fval else c ! THERE'S NOT ENOUGH WORKING SPACE ALLOCATED !!!! notify = 1 endif c c c END OF SUBROUTINE 'lnsrch' return end c c c======================================================================= c======================================================================= c subroutine lndata(dtil,btil,bm,zhi,zlo, # brk1,brk2,slope,zlft,zrt,zmid) double precision dtil,btil,bm double precision zhi,zlo,brk1,brk2,slope,zlft,zrt,zmid double precision eps c*********************************************************************** c Let EPS be a distinguishability constant. eps = 1.0D-10 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if ( dabs(dtil) .lt. eps ) then c !! In this case we consider dtil to be zero so that lambda c has no effect on the value of z, i.e. rho is constant. brk1 = -1.0d0 brk2 = -1.0d0 slope = 0.0d0 if ( btil .ge. bm*zhi ) then zlft = zhi else if ( btil .le. bm*zlo ) then zlft = zlo else zlft = btil/bm endif zmid = zlft zrt = zlft else if ( bm .gt. 0.0d0 ) then c !! In this case bm is positive so that rho is smooth c with two linear sections and a quadratic section in between. zmid = btil/bm slope = dtil/bm if (dtil .gt. 0.0d0 ) then zlft = zlo zrt = zhi brk1 = (bm*zlo - btil)/dtil brk2 = (bm*zhi - btil)/dtil else zlft = zhi zrt = zlo brk1 = (bm*zhi - btil)/dtil brk2 = (bm*zlo - btil)/dtil endif else c !! In this case bm is zero so that rho is nonsmooth c with two linear sections. zmid = zhi brk1 = -btil/dtil brk2 = -btil/dtil if ( dtil .gt. 0.0d0 ) then zlft = zlo zrt = zhi else zlft = zhi zrt = zlo endif endif c c END OF SUBROUTINE lndata return end :::::::::::::: pmain.f :::::::::::::: PROGRAM FGMDRV c c This is the driver routine for the DYNFGM program. It allocates the c storage space based on upper limits of the dimensions of the various c vectors, and then calls DYNFGM. c c======================================================================= c c Constant And Variable Declarations c PARAMETER ( mudim = 20, : muedim = 20, : mxydim = 20, : mvdim = 20, : mvedim = 20, : nmax = 16, : mxmode = 16 ) c c Description of the above constants: c mudim : maximum primal dimension c muedim : maximum primal endpoint dimension c mxydim : maximum trajectory dimension c mvdim : maximum dual dimension c mvedim : maximum dual endpoint dimension c nmax : maximum number of subdivisions of the unit interval (time) c mxmode : maximum number of elements generated in each iteration c PARAMETER ( maxnfe = 16 * (muedim+mvedim+2*mxydim : + nmax*(mudim+mvdim+2*mxydim)) : + 2*16*16, : mxvln = 8*(nmax*mudim+muedim), : mxuln = 8*(nmax*mvdim+mvedim), : mxinpt = 4*mxydim*mxydim + mudim + mvdim, : mxdata = mxydim*(mudim+muedim+mxydim+mvdim+mvedim+4) : + 4*(muedim+mvedim) + 2*(nmax+1)*(mudim+mvdim) : + 2*(muedim+mvedim+mudim+mvdim) + mudim*mvdim ) c c Description of the above constants: c maxnfe : space needed to store finite sets (controls, trajectories, c function values, and barycentric coordinates) c mxvln : workspace required by LNSRCH for dual linesearch c mxuln : workspace required by LNSRCH for primal linesearch c mxdata : space needed to store the data (u+,u-,P,p,Q,q,A,B,b, etc) c mxinpt : workspace required by computations within INPUT c c c Now maxwrk should be set to the largest of maxsad, mxvln, mxuln c and mxinpt. E.g. if nmax, mudim, muedim, xydim, mvdim, mvedim c are small ( nmax < 30, *dim < 6 ) maxsad may be the largest; if only c xydim is large then mxinpt may be the largest; if *dim are small and c nmax is large then may need to use either mxvln or mxuln. PARAMETER ( maxwrk = mxvln ) c c----------------------------------------------------------------------- c PARAMETER ( NZDLEN = maxnfe + maxwrk + mxdata ) DOUBLE PRECISION ZZD(NZDLEN) real tarray(2), tt c c----------------------------------------------------------------------- c CALL DYNFGM(ZZD,NZDLEN) c c----------------------------------------------------------------------- tt = dtime (tarray) print*, 'time:', tarray, tt c c End of program FGMDRV STOP END :::::::::::::: pobjvl.f :::::::::::::: SUBROUTINE OBJVL(U,V,PV,QV,PM,QM,DM,NU,NV,VALUE) c The objective value of the minimax problem is commputed for the optimal c solution found INTEGER NU,NV DOUBLE PRECISION U(NU),V(NV),PV(NU),QV(NV),PM(NU,NU),QM(NV,NV), : DM(NV,NU) DOUBLE PRECISION W1,W2,W3,S1,S2,VALUE W1=0.0D0 W3=0.0D0 S1=0.0D0 DO 300 I=1,NU DO 100 J=1,NU W1=W1+PM(I,J)*U(I)*U(J) 100 CONTINUE DO 200 J=1,NV W3=W3+V(J)*DM(J,I)*U(I) 200 CONTINUE S1=S1+PV(I)*U(I) 300 CONTINUE W2=0.0D0 S2=0.0D0 DO 500 I=1,NV DO 400 J=1,NV W2=W2+QM(I,J)*V(I)*V(J) 400 CONTINUE S2=S2+QV(I)*V(I) 500 CONTINUE VALUE = S1 + S2 + W1/2.0d0 - W2/2.0d0 - W3 RETURN c END OF subroutine 'OBJVAL' END :::::::::::::: poutput.f :::::::::::::: SUBROUTINE OUTPT0(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA,EPS,EPSNRM, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) INTEGER UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME DOUBLE PRECISION U(UDIM,TIME),V(VDIM,TIME) DOUBLE PRECISION Y(XYDIM,TIME),X(XYDIM,TIME) DOUBLE PRECISION UE(VEDIM),VE(VEDIM),XE(XYDIM),YE(XYDIM) DOUBLE PRECISION ALPHA,BETA,EPS,EPSNRM C INTEGER PROPT(3,5) COMMON /PRDATA/ PROPT C INTEGER OUT C----------------------------------------------------------------------- C 4001 FORMAT(2x,'OUTER ITERATION 0: gap = ',G12.6, $ ' relative gap = ',G12.6) C DO 191 JJ=1,3 OUT = PROPT(JJ,5) IF ( PROPT(JJ,1) .GT. 1 ) THEN c i.e. if using this file IF ( PROPT(JJ,2) .GE. 2 ) THEN c i.e. if we want readability then format the output WRITE(OUT,'(80("="))') IF ( PROPT(JJ,3) .GE. 2 ) THEN WRITE(OUT,4001) EPS, EPSNRM ELSE WRITE(OUT,'(2x,"OUTER LOOP ITERATION 0:")') ENDIF IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(/," u vector:")') WRITE(OUT,'(6(E12.2))') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(" ue vector:")') WRITE(OUT,'(6(E12.2))') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" x vector:")') WRITE(OUT,'(6(E12.2))') ((X(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(" xe vector:")') WRITE(OUT,'(6(E12.2))') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(34x,"primal objective value:",G20.14)') ALPHA IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(/," v vector:")') WRITE(OUT,'(6(E12.2))') ((V(I,K),I=1,VDIM),K=1,TIME) WRITE(OUT,'(" ve vector:")') WRITE(OUT,'(6(E12.2))') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" y vector:")') WRITE(OUT,'(6(E12.2))') ((Y(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(" ye vector:")') WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(36x,"dual objective value:",G20.14)') BETA ELSE IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT,*) EPS, EPSNRM IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((X(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2) WRITE(OUT,'(G26.19)') ALPHA IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((V(I,K),I=1,VDIM), K=1,TIME) WRITE(OUT,'(G26.19)') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((Y(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT,'(G26.19)') BETA ENDIF ENDIF 191 CONTINUE C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C End of subroutine 'OUTPT0' RETURN END SUBROUTINE OUTPT1(MU,DELTA) DOUBLE PRECISION DELTA INTEGER MU C INTEGER PROPT(3,5) COMMON /PRDATA/ PROPT C INTEGER OUT DO 191 JJ=1,3 IF ( PROPT(JJ,1) .GE. 2 ) THEN IF ( PROPT(JJ,2) .GE. 2 ) THEN OUT = PROPT(JJ,5) WRITE(OUT,'(80("="))') WRITE(OUT,'(2X,"OUTER LOOP ITERATION ",I3,":",11X, $ "Proximation parameter = ",D12.2)') MU,DELTA ENDIF ENDIF 191 CONTINUE RETURN C end of subroutine 'OUTPT1' END SUBROUTINE OUTPT2(INFORM,NU) INTEGER NU,INFORM C INTEGER PROPT(3,5) COMMON /PRDATA/ PROPT C INTEGER OUT C DO 191 J=1,3 OUT= PROPT(J,5) IF (PROPT(J,1) .GE. 2) THEN c >>>I.e., if using this file in all iterations then ... c IF (PROPT(J,2) .GE. 2) THEN c >>>I.e., if we want the file to be readable then format it. IF (INFORM .EQ. 0) THEN WRITE(OUT,1002) NU 1002 FORMAT(20X,'Inner algorithm converged in ',I2, $ ' iterations') WRITE(OUT,'(80("="))') ELSE IF (INFORM .EQ. 1) THEN WRITE(OUT,1003) NU 1003 FORMAT(20X,'Inner algorithm failed to converge in ', $ I2,' iterations') WRITE(OUT,'(80("="))') ENDIF ENDIF ENDIF 191 CONTINUE C RETURN C end of subroutine 'OUTPT2' END SUBROUTINE OUTPT3(U,V,X,Y,UE,VE,XE,YE, $ ALPHA,BETA,EPS,EPSNRM,RATIO, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) INTEGER UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME DOUBLE PRECISION U(UDIM,TIME),V(VDIM,TIME) DOUBLE PRECISION X(XYDIM,TIME),Y(XYDIM,TIME) DOUBLE PRECISION UE(UEDIM),VE(VEDIM) DOUBLE PRECISION XE(XYDIM),YE(XYDIM) DOUBLE PRECISION ALPHA,BETA,EPS,EPSNRM,RATIO C INTEGER PROPT(3,5) COMMON /PRDATA/ PROPT C INTEGER OUT DO 193 JJ=1,3 OUT = PROPT(JJ,5) IF (PROPT(JJ,1).GT.1) THEN c i.e. if using this file, and printing Best vectors then IF ( PROPT(JJ,2) .GE. 2 ) THEN c i.e. if we want readability then format the output WRITE(OUT,'(80("-"))') IF ( PROPT(JJ,3) .GE. 2 ) THEN WRITE(OUT,4000) EPS, EPSNRM, RATIO 4000 FORMAT(2x,'gap = ',G13.6,3x,'relative gap = ', $ G13.6,3x,'ratio = ',G13.6) ENDIF IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(/," u vector:")') WRITE(OUT,'(6(E12.2))') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(" ue vector:")') WRITE(OUT,'(6(E12.2))') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" x vector:")') WRITE(OUT,'(6(E12.2))') ((X(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(" xe vector:")') WRITE(OUT,'(6(E12.2))') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(34x,"primal objective value:",G20.14)') ALPHA IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(/," v vector:")') WRITE(OUT,'(6(E12.2))') ((V(I,K),I=1,VDIM),K=1,TIME) WRITE(OUT,'(" ve vector:")') WRITE(OUT,'(6(E12.2))') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" y vector:")') WRITE(OUT,'(6(E12.2))') ((Y(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(" ye vector:")') WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(36x,"dual objective value:",G20.14)') BETA ELSE IF ( PROPT(JJ,3) .GE. 2 ) THEN WRITE(OUT,*) EPS, EPSNRM, RATIO ENDIF IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(G26.19)') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((X(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2) WRITE(OUT,'(G26.19)') ALPHA IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(G26.19)') ((V(I,K),I=1,VDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((Y(J,K),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT,'(G26.19)') BETA ENDIF ENDIF 193 CONTINUE C C----------------------------------------------------------------------- C RETURN C End of subroutine 'OUTPT3' END SUBROUTINE OUTPT4(MU,INFORM,EPS,EPSNRM,ALPHA,BETA, : U,V,UE,VE,DVAL,UDIM,UEDIM,VDIM,VEDIM,N) C====================================================================== C VARIABLES PASSED IN: C ------------------- C INTEGER MU,INFORM INTEGER UDIM,UEDIM,VDIM,VEDIM,N DOUBLE PRECISION U(UDIM,N),V(VDIM,N) DOUBLE PRECISION UE(UEDIM), VE(VEDIM) DOUBLE PRECISION EPS,EPSNRM,ALPHA,BETA,DVAL C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Common Blocks: C 1. PRDATA contains "PRint DATA". The numbers PROPT(j,5) (for j=1,2,3) C give the unit numbers for the various output files. Typically, j=1 C corresponds to the standard output unit, i.e. PROPT(1,5)=6. The other C entries of PROPT should represent the various printing options chosen C by the user: it allows for three output units with four integer C descriptions for each. INTEGER PROPT(3,5) COMMON /PRDATA/ PROPt C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C local variables: C --------------- C INTEGER OUt C====================================================================== C DO 191 J=1,3 OUT= PROPT(J,5) C IF (PROPT(J,1) .GE. 2) THEN C >>>I.e., if using this file in all iterations then ... C IF (PROPT(J,2) .GE. 2) THEN C >>>I.e., if we want the file to be readable then format it. IF (INFORM .EQ. 0) THEN WRITE(OUT,1002) MU 1002 FORMAT(20X,'Outer algorithm converged in ',I3, $ ' iterations') WRITE(OUT,1001) DVAL 1001 FORMAT(20X,'Additive constant is ',G26.19) WRITE(OUT,'(80("="))') ELSE IF (INFORM .EQ. 1) THEN WRITE(OUT,1003) MU 1003 FORMAT(20X,'Outer algorithm failed to converge in ', $ I3,' iterations') WRITE(OUT,'(80("="))') ENDIF ELSE C >>>We are printing to a free format file, so print the additive C constant and signal the END-OF-DATA. WRITE(OUT,*) DVAL WRITE(OUT,'("ENDDATA")') ENDIF ELSE IF (PROPT(J,1) .EQ. 1) THEN C >>>Only the final solution is to be printed. WRITE(OUT,'(G26.19)') ((U(I,K),I=1,UDIM),K=1,N) WRITE(OUT,'(G26.19)') (UE(I),I=1,UEDIM) WRITE(OUT,'(G26.19)') ((V(I,K),I=1,VDIM),K=1,N) WRITE(OUT,'(G26.19)') (VE(I),I=1,VEDIM) WRITE(OUT,*) ( ALPHA + BETA )/2.0D0 WRITE(OUT,*) DVAL ENDIF 191 CONTINUE C C C C End of subroutine 'OUTPT4' RETURN END SUBROUTINE OUTLNS(U,V,X,Y,UE,VE,XE,YE, $ ALPHA,BETA,EPS,EPSNRM,RATIO,USTEP,VSTEP, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) INTEGER UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME DOUBLE PRECISION U(UDIM,TIME),V(VDIM,TIME) DOUBLE PRECISION X(XYDIM,TIME),Y(XYDIM,TIME) DOUBLE PRECISION UE(UEDIM),VE(VEDIM) DOUBLE PRECISION XE(XYDIM),YE(XYDIM) DOUBLE PRECISION ALPHA,BETA,EPS,EPSNRM,RATIO,USTEP,VSTEP C INTEGER PROPT(3,5) COMMON /PRDATA/ PROPT C INTEGER OUT C DO 193 JJ=1,3 OUT = PROPT(JJ,5) IF (PROPT(JJ,1).GT.1) THEN c i.e. if using this file, and printing Best vectors then IF ( PROPT(JJ,2) .GE. 2 ) THEN c i.e. if we want readability then format the output WRITE(OUT,'(80("-"))') IF ( PROPT(JJ,3) .GE. 2 ) THEN WRITE(OUT,4000) EPS, EPSNRM, RATIO, USTEP, VSTEP 4000 FORMAT(2x,'gap = ',G13.6,3x,'relative gap = ', $ G13.6,3x,'ratio = ',G13.6,/,2x,'steps: ', $ G11.4,',',1X,G11.4) ENDIF IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(/," u vector:")') WRITE(OUT,'(6(E12.2))') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(" ue vector:")') WRITE(OUT,'(6(E12.2))') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" x vector:")') WRITE(OUT,'(6(E12.2))') ((X(J,K),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(" xe vector:")') WRITE(OUT,'(6(E12.2))') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(34x,"primal objective value:",G20.14)') ALPHA IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(/," v vector:")') WRITE(OUT,'(6(E12.2))') ((V(I,K),I=1,VDIM),K=1,TIME) WRITE(OUT,'(" ve vector:")') WRITE(OUT,'(6(E12.2))') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" y vector:")') WRITE(OUT,'(6(E12.2))') ((Y(J,K),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(" ye vector:")') WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(36x,"dual objective value:",G20.14)') BETA ELSE IF ( PROPT(JJ,3) .GE. 2 ) THEN WRITE(OUT,*) EPS, EPSNRM, RATIO, USTEP, VSTEP ENDIF IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(G26.19)') ((U(J,K),J=1,UDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (UE(J), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((X(J,K),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(G26.19)') (XE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2) WRITE(OUT,'(G26.19)') ALPHA IF (MOD(PROPT(JJ,4),2).EQ.1) THEN WRITE(OUT,'(G26.19)') ((V(I,K),I=1,VDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (VE(J), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((Y(J,K),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(6(E12.2))') (YE(J), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT,'(G26.19)') BETA ENDIF ENDIF 193 CONTINUE C C----------------------------------------------------------------------- C RETURN C End of subroutine 'OUTLNS' END :::::::::::::: pprox.f :::::::::::::: SUBROUTINE PROX(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PMAT,PVEC,QMAT,QVEC, $ DMAT,PEMAT,PEVEC,QEMAT,QEVEC,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ PD,PV,QD,QV,PED,PEV,QED,QEV,ALPHA,BETA,DVAL, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ MAXUFE,MAXVFE,MUMAX,NUMAX,DEPS) C C C----------------------------------------------------------------------- C/////////////////////////////////////////////////////////////////////// C/////////////////////////////////////////////////////////////////////// C----------------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Variables passed thru subroutine call C INTEGER LENWRK INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME INTEGER MAXUFE,MAXVFE,NUMAX,MUMAX INTEGER BESTU,BESTV INTEGER MU DOUBLE PRECISION DEPS C C DOUBLE PRECISION U(UDIM,TIME,MAXUFE),UE(UEDIM,MAXUFE) DOUBLE PRECISION V(VDIM,TIME,MAXVFE),VE(VEDIM,MAXVFE) DOUBLE PRECISION X(XYDIM,TIME,MAXUFE),XE(XYDIM,MAXUFE) DOUBLE PRECISION Y(XYDIM,TIME,MAXVFE),YE(XYDIM,MAXVFE) C DOUBLE PRECISION UUPR(UDIM),ULOW(UDIM),UEUPR(UEDIM),UELOW(UEDIM) DOUBLE PRECISION VUPR(VDIM),VLOW(VDIM),VEUPR(VEDIM),VELOW(VEDIM) C DOUBLE PRECISION PMAT(UDIM),PVEC(UDIM,TIME) DOUBLE PRECISION QMAT(VDIM),QVEC(VDIM,TIME) DOUBLE PRECISION DMAT(VDIM,UDIM) DOUBLE PRECISION PEMAT(UEDIM),PEVEC(UEDIM) DOUBLE PRECISION QEMAT(VEDIM),QEVEC(VEDIM) C DOUBLE PRECISION AMAT(XYDIM,XYDIM) DOUBLE PRECISION BMAT(XYDIM,UDIM),BVEC(XYDIM) DOUBLE PRECISION CMAT(VDIM,XYDIM),CVEC(XYDIM) DOUBLE PRECISION BEMAT(XYDIM,UEDIM),BEVEC(XYDIM) DOUBLE PRECISION CEMAT(VEDIM,XYDIM),CEVEC(XYDIM) C DOUBLE PRECISION XI(MAXUFE),ETA(MAXVFE) C DOUBLE PRECISION PD(UDIM),PV(UDIM,TIME) DOUBLE PRECISION PED(UEDIM),PEV(UEDIM) DOUBLE PRECISION QD(VDIM),QV(VDIM,TIME) DOUBLE PRECISION QED(VEDIM),QEV(VEDIM) C DOUBLE PRECISION ALPHA(MAXUFE),BETA(MAXVFE),DVAL C DOUBLE PRECISION WORK(LENWRK) C c*********************************************************************** c c This subroutine takes as input the data for the following problem: c c minimax J(u,u ;v,v ) over ( U x U ) x ( V x V ). c e e e e c c Here U is a closed box in K-by-N dimensional Euclidean space and U is c e c a closed box in K -dimensional Euclidean space, so that the primal vectors c e c have the form (u,u ), where u is of the form u(i,t) with i=1,...,K and c e c t=1,...,N, and u has the form u (i) with i=1,...,K . c e e e c Similarly, V is a closed box of points of the form v(j,t) with j=1,...,L c and t=1,...,N, and V is defined similarly to U . c e e c c J has the form c N c J(u,u ;v,v ) = [ sum J (u(t),v(t)) ] + J (u ,v ) - < u,u ; v,v > c e e t=1 t e e e e e c c c 1 1 c where J (u(t),v(t)) = p u + q v + -u Pu - -v Qv - v Du , c t t t t t 2 t t 2 t t t t c c 1 1 c J (u ,v ) = p u + q v + - u P u - -v Q v c e e e e e e e 2 e e e 2 e e e c c N T T c < u,u ; v,v > = [ sum (C v +c)x ] + (C v +c )x c e e t=1 t t-1 e e e e c c and x = Ax + Bu + b for t=1,...,N c t t-1 t c c x = x = B u + b . c 0 e e e e c c Associated with this saddle-point problem are the primal and dual problems: c c min f(u,u ) and max g(v,v ) c over UxU e over VxV e c e e c c where f(u,u ) = max J(u,u ;v,v ) and g(v,v ) = min J(u,u ;v,v ) . c e over VxV e e e over UxU e e c e e c c It is assumed that P , P , Q , Q , are all diagonal matrices with nonnegative c e e c entries on the diagonal. Also associated with the problem are the maps c F(u,u ) = argmax { J(u,u ;v,v ):(v,v ) in VxV } c e e e e e c and G(v,v ) = argmin { J(u,u ;v,v ) : (u,u ) in UxU } . c e e e e e c In the documentation below, the term "alternative element" will always refer c to an element given by applying F or G to the current "best" guess. c c c This routine runs the outer loop associated with Proximal-Point Algorithm. c It calls the subroutine DFGM to solve a subproblem of the above form for c which the matrices P , P , Q and Q are positive definite. c e e c It is expected that the method employed within DFGM will be similar in c structure to the Finite-Generation Method. For this reason, the data c structures have been defined so as to keep track of "finite sets" of primal c and dual vectors. The structure of the algorithm is as follows: c c Given: an initial guess for the vectors u, u , v, v , and stopping c e e c criteria "deps" (double precision floating point number) c and a positive integer "mumax". c Step 0. (a) Initialization : compute the function values for the c initial guesses and compute the trajectories x and y c of the initial guesses. c (b) Set mu=1 and delta=1/10. c Decision Step: If the normalized duality gap is less than 'DEPS', c then stop: the pair (u,ue) and the pair (v,ve) are both c EPSILON-optimal. Otherwise, if mu>>>>>> INITIALIZATION <<<<<<<< C BESTU = 1 BESTV = 1 ALTU = 2 ALTV = 2 OLDU = MAXUFE OLDV = MAXVFE C C Call routine to calculate the trajectories X and Y for U , V C CALL TRAJ(U(1,1,BESTU),UE(1,BESTU),V(1,1,BESTV),VE(1,BESTV), : X(1,1,BESTU),XE(1,BESTU),Y(1,1,BESTV),YE(1,BESTV), : AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C C Compute the function value ALPHA for U = BESTU and BETA for V = BESTV. C CALL FGVAL(U(1,1,ALTU),V(1,1,ALTV),UE(1,ALTU),VE(1,ALTV), $ ALPHA(BESTU),BETA(BESTV), $ U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV), $ XE(1,BESTU),YE(1,BESTV), $ UUPR,ULOW,VLOW,VUPR,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C Compute the Duality Gap for BESTU,BESTV. C EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF OLDEPS = EPS C CALL OUTPT0(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA,EPS,EPSNRM, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) C C Number Of Iterations: MU = 0 C C Proximation parameter: cccc DELTA = 1.D0 DELTA = 1.D-1 C C//////////////////////////// C>>>>>>> Begin Loop <<<<<<<< C>>>>>>> ========== <<<<<<<< C//////////////////////////// C 1000 CONTINUE C IF (( EPSNRM .GT. DEPS ) .AND. ( MU .LT. MUMAX )) THEN C C i.e. While ( EPS > DEPS and MU < MUMAX ) C ===== --- ---- === -- ----- C Do the following: C == === ========= C MU = MU + 1 DO 40 K=1,TIME DO 10 J=1,UDIM 10 U(J,K,OLDU) = U(J,K,BESTU) DO 20 I=1,VDIM 20 V(I,K,OLDV) = V(I,K,BESTV) DO 30 J=1,XYDIM X(J,K,OLDU) = X(J,K,BESTU) 30 Y(J,K,OLDV) = Y(J,K,BESTV) 40 CONTINUE DO 50 J=1,UEDIM 50 UE(J,OLDU) = UE(J,BESTU) DO 60 I=1,VEDIM 60 VE(I,OLDV) = VE(I,BESTV) DO 70 J=1,XYDIM XE(J,OLDU) = XE(J,BESTU) YE(J,OLDV) = YE(J,BESTV) 70 CONTINUE C CALL OUTPT1(MU,DELTA) C DO 80 J=1,UDIM 80 PD(J) = DELTA+PMAT(J) DO 90 I=1,VDIM 90 QD(I) = DELTA+QMAT(I) DO 120 K=1,TIME DO 100 J=1,UDIM PV(J,K) = PVEC(J,K)-DELTA*U(J,K,BESTU) 100 CONTINUE DO 110 I=1,VDIM QV(I,K) = QVEC(I,K)+DELTA*V(I,K,BESTV) 110 CONTINUE 120 CONTINUE DO 130 J=1,UEDIM PED(J) = DELTA+PEMAT(J) PEV(J) = PEVEC(J)-DELTA*UE(J,BESTU) 130 CONTINUE DO 140 I=1,VEDIM QED(I) = DELTA+QEMAT(I) QEV(I) = QEVEC(I)+DELTA*VE(I,BESTV) 140 CONTINUE C C CALL DFGM(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PD,PV,QD,QV, $ DMAT,PED,PEV,QED,QEV,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ ALPHA,BETA,BESTU,BESTV,EPS,EPSNRM,INFORM,NU, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ 8,8,NUMAX,DEPS/1.d8) C CALL OUTPT2(INFORM,NU) C DDD = 1.d-10 C RESTART = MU -1 cccc DUM = 1.d4 DUM = 1.d2/2.d0**(RESTART - 4) DO 1170 K=1,TIME DO 1150 J=1,UDIM IF( U(J,K,OLDU)-U(J,K,BESTU) .GT. DDD ) : DUM = DMIN1((U(J,K,BESTU)-ULOW(J))/ : (U(J,K,OLDU)-U(J,K,BESTU)), DUM) 1150 IF( U(J,K,BESTU)-U(J,K,OLDU) .GT. DDD ) : DUM = DMIN1((UUPR(J)-U(J,K,BESTU))/ : (U(J,K,BESTU)-U(J,K,OLDU)), DUM) 1170 CONTINUE DO 1180 J=1,UEDIM IF( UE(J,OLDU)-UE(J,BESTU) .GT. DDD ) : DUM = DMIN1((UE(J,BESTU)-UELOW(J))/ : (UE(J,OLDU)-UE(J,BESTU)), DUM) 1180 IF( UE(J,BESTU)-UE(J,OLDU) .GT. DDD ) : DUM = DMIN1((UEUPR(J)-UE(J,BESTU))/ : (UE(J,BESTU)-UE(J,OLDU)), DUM) C DO 170 K=1,TIME DO 150 J=1,UDIM 150 U(J,K,1) = U(J,K,BESTU) + DUM * : (U(J,K,BESTU)-U(J,K,OLDU)) 170 CONTINUE DO 180 J=1,UEDIM 180 UE(J,1) = UE(J,BESTU) + DUM * : (UE(J,BESTU)-UE(J,OLDU)) C cccc DUM = 1.d4 DUM = 1.d2/2.d0**(RESTART - 4) DO 2170 K=1,TIME DO 2150 I=1,VDIM IF( V(I,K,OLDV)-V(I,K,BESTV) .GT. DDD ) : DUM = DMIN1((V(I,K,BESTV)-VLOW(I))/ : (V(I,K,OLDV)-V(I,K,BESTV)), DUM) 2150 IF( V(I,K,BESTV)-V(I,K,OLDV) .GT. DDD ) : DUM = DMIN1((VUPR(I)-V(I,K,BESTV))/ : (V(I,K,BESTV)-V(I,K,OLDV)), DUM) 2170 CONTINUE DO 2180 I=1,VEDIM IF( VE(I,OLDV)-VE(I,BESTV) .GT. DDD ) : DUM = DMIN1((VE(I,BESTV)-VELOW(I))/ : (VE(I,OLDV)-VE(I,BESTV)), DUM) 2180 IF( VE(I,BESTV)-VE(I,OLDV) .GT. DDD ) : DUM = DMIN1((VEUPR(I)-VE(I,BESTV))/ : (VE(I,BESTV)-VE(I,OLDV)), DUM) C DO 270 K=1,TIME DO 250 I=1,VDIM 250 V(I,K,1) = V(I,K,BESTV) + DUM * : (V(I,K,BESTV)-V(I,K,OLDV)) 270 CONTINUE DO 280 I=1,VEDIM 280 VE(I,1) = VE(I,BESTV) + DUM * : (VE(I,BESTV)-VE(I,OLDV)) C UDIR = 1 VDIR = 1 CALL TRAJ(U(1,1,UDIR ),UE(1,UDIR ),V(1,1,VDIR ),VE(1,VDIR ), : X(1,1,UDIR ),XE(1,UDIR ),Y(1,1,VDIR ),YE(1,VDIR ), : AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C EPSOLD = OLDEPS C CALL LINE(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ BESTU,BESTV,UDIR,VDIR,1,1,ALTU,ALTV, $ UUPR,ULOW,VUPR,VLOW,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT, $ PEMAT,PEVEC,QEMAT,QEVEC, $ AMAT,BMAT,BVEC,CMAT,CVEC, $ BEMAT,BEVEC,CEMAT,CEVEC, $ MAXUFE,MAXVFE, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME, $ WORK,LENWRK,QUIT) C BESTU = 1 BESTV = 1 C C Compute the Duality Gap for BESTU,BESTV. C EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF C OLDEPS = EPS C C Update DELTA if convergence was quick enough. C IF (INFORM .EQ. 0 .AND. $ NU .LE. MAX(TIME*UDIM+UEDIM,TIME*VDIM+VEDIM)) THEN DELTA = DELTA/10.0d0 ENDIF C C Continue iterating ! GO TO 1000 ENDIF C////////////////////////// C>>>>>>> End Loop <<<<<<<< C>>>>>>> ======== <<<<<<<< C////////////////////////// C C C IF ( EPSNRM .LE. DEPS ) THEN C Successful termination ! INFORM = 0 ELSE C Maximum number of iterations exceeded ! INFORM = 1 ENDIF C CALL OUTPT4(MU,INFORM,EPS,EPSNRM,ALPHA,BETA, : U,V,UE,VE,DVAL,UDIM,UEDIM,VDIM,VEDIM,TIME) C======================================================================= C RETURN C End of subroutine 'DFGM' END :::::::::::::: proxo.f :::::::::::::: SUBROUTINE PROX(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PMAT,PVEC,QMAT,QVEC, $ DMAT,PEMAT,PEVEC,QEMAT,QEVEC,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ PD,PV,QD,QV,PED,PEV,QED,QEV,ALPHA,BETA,DVAL, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ MAXUFE,MAXVFE,MUMAX,NUMAX,DEPS) C C C----------------------------------------------------------------------- C/////////////////////////////////////////////////////////////////////// C/////////////////////////////////////////////////////////////////////// C----------------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Variables passed thru subroutine call C INTEGER LENWRK INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME INTEGER MAXUFE,MAXVFE,NUMAX,MUMAX INTEGER BESTU,BESTV INTEGER MU DOUBLE PRECISION DEPS C C DOUBLE PRECISION U(UDIM,TIME,MAXUFE),UE(UEDIM,MAXUFE) DOUBLE PRECISION V(VDIM,TIME,MAXVFE),VE(VEDIM,MAXVFE) DOUBLE PRECISION X(XYDIM,TIME,MAXUFE),XE(XYDIM,MAXUFE) DOUBLE PRECISION Y(XYDIM,TIME,MAXVFE),YE(XYDIM,MAXVFE) C DOUBLE PRECISION UUPR(UDIM),ULOW(UDIM),UEUPR(UEDIM),UELOW(UEDIM) DOUBLE PRECISION VUPR(VDIM),VLOW(VDIM),VEUPR(VEDIM),VELOW(VEDIM) C DOUBLE PRECISION PMAT(UDIM),PVEC(UDIM,TIME) DOUBLE PRECISION QMAT(VDIM),QVEC(VDIM,TIME) DOUBLE PRECISION DMAT(VDIM,UDIM) DOUBLE PRECISION PEMAT(UEDIM),PEVEC(UEDIM) DOUBLE PRECISION QEMAT(VEDIM),QEVEC(VEDIM) C DOUBLE PRECISION AMAT(XYDIM,XYDIM) DOUBLE PRECISION BMAT(XYDIM,UDIM),BVEC(XYDIM) DOUBLE PRECISION CMAT(VDIM,XYDIM),CVEC(XYDIM) DOUBLE PRECISION BEMAT(XYDIM,UEDIM),BEVEC(XYDIM) DOUBLE PRECISION CEMAT(VEDIM,XYDIM),CEVEC(XYDIM) C DOUBLE PRECISION XI(MAXUFE),ETA(MAXVFE) C DOUBLE PRECISION PD(UDIM),PV(UDIM,TIME) DOUBLE PRECISION PED(UEDIM),PEV(UEDIM) DOUBLE PRECISION QD(VDIM),QV(VDIM,TIME) DOUBLE PRECISION QED(VEDIM),QEV(VEDIM) C DOUBLE PRECISION ALPHA(MAXUFE),BETA(MAXVFE),DVAL C DOUBLE PRECISION WORK(LENWRK) C c*********************************************************************** c c This subroutine takes as input the data for the following problem: c c minimax J(u,u ;v,v ) over ( U x U ) x ( V x V ). c e e e e c c Here U is a closed box in K-by-N dimensional Euclidean space and U is c e c a closed box in K -dimensional Euclidean space, so that the primal vectors c e c have the form (u,u ), where u is of the form u(i,t) with i=1,...,K and c e c t=1,...,N, and u has the form u (i) with i=1,...,K . c e e e c Similarly, V is a closed box of points of the form v(j,t) with j=1,...,L c and t=1,...,N, and V is defined similarly to U . c e e c c J has the form c N c J(u,u ;v,v ) = [ sum J (u(t),v(t)) ] + J (u ,v ) - < u,u ; v,v > c e e t=1 t e e e e e c c c 1 1 c where J (u(t),v(t)) = p u + q v + -u Pu - -v Qv - v Du , c t t t t t 2 t t 2 t t t t c c 1 1 c J (u ,v ) = p u + q v + - u P u - -v Q v c e e e e e e e 2 e e e 2 e e e c c N T T c < u,u ; v,v > = [ sum (C v +c)x ] + (C v +c )x c e e t=1 t t-1 e e e e c c and x = Ax + Bu + b for t=1,...,N c t t-1 t c c x = x = B u + b . c 0 e e e e c c Associated with this saddle-point problem are the primal and dual problems: c c min f(u,u ) and max g(v,v ) c over UxU e over VxV e c e e c c where f(u,u ) = max J(u,u ;v,v ) and g(v,v ) = min J(u,u ;v,v ) . c e over VxV e e e over UxU e e c e e c c It is assumed that P , P , Q , Q , are all diagonal matrices with nonnegative c e e c entries on the diagonal. Also associated with the problem are the maps c F(u,u ) = argmax { J(u,u ;v,v ):(v,v ) in VxV } c e e e e e c and G(v,v ) = argmin { J(u,u ;v,v ) : (u,u ) in UxU } . c e e e e e c In the documentation below, the term "alternative element" will always refer c to an element given by applying F or G to the current "best" guess. c c c This routine runs the outer loop associated with Proximal-Point Algorithm. c It calls the subroutine DFGM to solve a subproblem of the above form for c which the matrices P , P , Q and Q are positive definite. c e e c It is expected that the method employed within DFGM will be similar in c structure to the Finite-Generation Method. For this reason, the data c structures have been defined so as to keep track of "finite sets" of primal c and dual vectors. The structure of the algorithm is as follows: c c Given: an initial guess for the vectors u, u , v, v , and stopping c e e c criteria "deps" (double precision floating point number) c and a positive integer "mumax". c Step 0. (a) Initialization : compute the function values for the c initial guesses and compute the trajectories x and y c of the initial guesses. c (b) Set mu=1 and delta=1/10. c Decision Step: If the normalized duality gap is less than 'DEPS', c then stop: the pair (u,ue) and the pair (v,ve) are both c EPSILON-optimal. Otherwise, if mu>>>>>> INITIALIZATION <<<<<<<< C BESTU = 1 BESTV = 1 ALTU = 2 ALTV = 2 OLDU = MAXUFE OLDV = MAXVFE C C Call routine to calculate the trajectories X and Y for U , V C CALL TRAJ(U(1,1,BESTU),UE(1,BESTU),V(1,1,BESTV),VE(1,BESTV), : X(1,1,BESTU),XE(1,BESTU),Y(1,1,BESTV),YE(1,BESTV), : AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C C Compute the function value ALPHA for U = BESTU and BETA for V = BESTV. C CALL FGVAL(U(1,1,ALTU),V(1,1,ALTV),UE(1,ALTU),VE(1,ALTV), $ ALPHA(BESTU),BETA(BESTV), $ U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV), $ XE(1,BESTU),YE(1,BESTV), $ UUPR,ULOW,VLOW,VUPR,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C Compute the Duality Gap for BESTU,BESTV. C EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF OLDEPS = EPS C CALL OUTPT0(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA,EPS,EPSNRM, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) C C Number Of Iterations: MU = 0 C C Proximation parameter: cccc DELTA = 1.D0 DELTA = 1.D-1 C C//////////////////////////// C>>>>>>> Begin Loop <<<<<<<< C>>>>>>> ========== <<<<<<<< C//////////////////////////// C 1000 CONTINUE C IF (( EPSNRM .GT. DEPS ) .AND. ( MU .LT. MUMAX )) THEN C C i.e. While ( EPS > DEPS and MU < MUMAX ) C ===== --- ---- === -- ----- C Do the following: C == === ========= C MU = MU + 1 C CALL OUTPT1(MU,DELTA) C DO 80 J=1,UDIM 80 PD(J) = DELTA+PMAT(J) DO 90 I=1,VDIM 90 QD(I) = DELTA+QMAT(I) DO 120 K=1,TIME DO 100 J=1,UDIM PV(J,K) = PVEC(J,K)-DELTA*U(J,K,BESTU) 100 CONTINUE DO 110 I=1,VDIM QV(I,K) = QVEC(I,K)+DELTA*V(I,K,BESTV) 110 CONTINUE 120 CONTINUE DO 130 J=1,UEDIM PED(J) = DELTA+PEMAT(J) PEV(J) = PEVEC(J)-DELTA*UE(J,BESTU) 130 CONTINUE DO 140 I=1,VEDIM QED(I) = DELTA+QEMAT(I) QEV(I) = QEVEC(I)+DELTA*VE(I,BESTV) 140 CONTINUE C C CALL DFGM(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PD,PV,QD,QV, $ DMAT,PED,PEV,QED,QEV,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ ALPHA,BETA,BESTU,BESTV,EPS,EPSNRM,INFORM,NU, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ 8,8,NUMAX,DEPS/1.d8) C CALL OUTPT2(INFORM,NU) C IF ( BESTU .GT. 1 ) THEN DO 170 K=1,TIME DO 150 J=1,UDIM 150 U(J,K,1) = U(J,K,BESTU) DO 160 J=1,XYDIM 160 X(J,K,1) = X(J,K,BESTU) 170 CONTINUE DO 180 J=1,UEDIM 180 UE(J,1) = UE(J,BESTU) DO 190 J=1,XYDIM 190 XE(J,1) = XE(J,BESTU) BESTU = 1 ENDIF C IF ( BESTV .GT. 1 ) THEN DO 220 K=1,TIME DO 200 I=1,VDIM 200 V(I,K,1) = V(I,K,BESTV) DO 210 I=1,XYDIM 210 Y(I,K,1) = Y(I,K,BESTV) 220 CONTINUE DO 230 I=1,VEDIM 230 VE(I,1) = VE(I,BESTV) DO 240 I=1,XYDIM 240 YE(I,1) = YE(I,BESTV) BESTV = 1 ENDIF C UDIR = BESTU VDIR = BESTV C EPSOLD = OLDEPS C C C - - - - - - - - - - - - - - - - - - - - - - - - - CALL FGVAL(U(1,1,ALTU),V(1,1,ALTV),UE(1,ALTU),VE(1,ALTV), $ ALPHA(BESTU),BETA(BESTV), $ U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV), $ XE(1,BESTU),YE(1,BESTV), $ UUPR,ULOW,VLOW,VUPR,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C Compute the Duality Gap for BESTU,BESTV. C RATIO = (ALPHA(BESTU) - BETA(BESTV))/EPSOLD EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF C USTEP = 1. VSTEP = 1. C CALL OUTLNS(U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV),XE(1,BESTU),YE(1,BESTV), $ ALPHA(BESTU),BETA(BESTV),EPS,EPSNRM,RATIO,USTEP,VSTEP, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) C C----------------------------------------------------------------------- C C OLDEPS = EPS C C Update DELTA if convergence was quick enough. C cccc IF (INFORM .EQ. 0 .AND. cccc $ NU .LE. MAX(TIME*UDIM+UEDIM,TIME*VDIM+VEDIM)) THEN cccc DELTA = DELTA/10.0d0 cccc ENDIF C C Continue iterating ! GO TO 1000 ENDIF C////////////////////////// C>>>>>>> End Loop <<<<<<<< C>>>>>>> ======== <<<<<<<< C////////////////////////// C C C IF ( EPSNRM .LE. DEPS ) THEN C Successful termination ! INFORM = 0 ELSE C Maximum number of iterations exceeded ! INFORM = 1 ENDIF C CALL OUTPT4(MU,INFORM,EPS,EPSNRM,ALPHA,BETA, : U,V,UE,VE,DVAL,UDIM,UEDIM,VDIM,VEDIM,TIME) C======================================================================= C RETURN C End of subroutine 'DFGM' END :::::::::::::: prxbz.f :::::::::::::: SUBROUTINE PROX(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PMAT,PVEC,QMAT,QVEC, $ DMAT,PEMAT,PEVEC,QEMAT,QEVEC,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ PD,PV,QD,QV,PED,PEV,QED,QEV,ALPHA,BETA,DVAL, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ MAXUFE,MAXVFE,MUMAX,NUMAX,DEPS) C C C----------------------------------------------------------------------- C/////////////////////////////////////////////////////////////////////// C/////////////////////////////////////////////////////////////////////// C----------------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Variables passed thru subroutine call C INTEGER LENWRK INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME INTEGER MAXUFE,MAXVFE,NUMAX,MUMAX INTEGER BESTU,BESTV INTEGER MU DOUBLE PRECISION DEPS C C DOUBLE PRECISION U(UDIM,TIME,MAXUFE),UE(UEDIM,MAXUFE) DOUBLE PRECISION V(VDIM,TIME,MAXVFE),VE(VEDIM,MAXVFE) DOUBLE PRECISION X(XYDIM,TIME,MAXUFE),XE(XYDIM,MAXUFE) DOUBLE PRECISION Y(XYDIM,TIME,MAXVFE),YE(XYDIM,MAXVFE) C DOUBLE PRECISION UUPR(UDIM),ULOW(UDIM),UEUPR(UEDIM),UELOW(UEDIM) DOUBLE PRECISION VUPR(VDIM),VLOW(VDIM),VEUPR(VEDIM),VELOW(VEDIM) C DOUBLE PRECISION PMAT(UDIM),PVEC(UDIM,TIME) DOUBLE PRECISION QMAT(VDIM),QVEC(VDIM,TIME) DOUBLE PRECISION DMAT(VDIM,UDIM) DOUBLE PRECISION PEMAT(UEDIM),PEVEC(UEDIM) DOUBLE PRECISION QEMAT(VEDIM),QEVEC(VEDIM) C DOUBLE PRECISION AMAT(XYDIM,XYDIM) DOUBLE PRECISION BMAT(XYDIM,UDIM),BVEC(XYDIM) DOUBLE PRECISION CMAT(VDIM,XYDIM),CVEC(XYDIM) DOUBLE PRECISION BEMAT(XYDIM,UEDIM),BEVEC(XYDIM) DOUBLE PRECISION CEMAT(VEDIM,XYDIM),CEVEC(XYDIM) C DOUBLE PRECISION XI(MAXUFE),ETA(MAXVFE) C DOUBLE PRECISION PD(UDIM),PV(UDIM,TIME) DOUBLE PRECISION PED(UEDIM),PEV(UEDIM) DOUBLE PRECISION QD(VDIM),QV(VDIM,TIME) DOUBLE PRECISION QED(VEDIM),QEV(VEDIM) C DOUBLE PRECISION ALPHA(MAXUFE),BETA(MAXVFE),DVAL C DOUBLE PRECISION WORK(LENWRK) C c*********************************************************************** c c This subroutine takes as input the data for the following problem: c c minimax J(u,u ;v,v ) over ( U x U ) x ( V x V ). c e e e e c c Here U is a closed box in K-by-N dimensional Euclidean space and U is c e c a closed box in K -dimensional Euclidean space, so that the primal vectors c e c have the form (u,u ), where u is of the form u(i,t) with i=1,...,K and c e c t=1,...,N, and u has the form u (i) with i=1,...,K . c e e e c Similarly, V is a closed box of points of the form v(j,t) with j=1,...,L c and t=1,...,N, and V is defined similarly to U . c e e c c J has the form c N c J(u,u ;v,v ) = [ sum J (u(t),v(t)) ] + J (u ,v ) - < u,u ; v,v > c e e t=1 t e e e e e c c c 1 1 c where J (u(t),v(t)) = p u + q v + -u Pu - -v Qv - v Du , c t t t t t 2 t t 2 t t t t c c 1 1 c J (u ,v ) = p u + q v + - u P u - -v Q v c e e e e e e e 2 e e e 2 e e e c c N T T c < u,u ; v,v > = [ sum (C v +c)x ] + (C v +c )x c e e t=1 t t-1 e e e e c c and x = Ax + Bu + b for t=1,...,N c t t-1 t c c x = x = B u + b . c 0 e e e e c c Associated with this saddle-point problem are the primal and dual problems: c c min f(u,u ) and max g(v,v ) c over UxU e over VxV e c e e c c where f(u,u ) = max J(u,u ;v,v ) and g(v,v ) = min J(u,u ;v,v ) . c e over VxV e e e over UxU e e c e e c c It is assumed that P , P , Q , Q , are all diagonal matrices with nonnegative c e e c entries on the diagonal. Also associated with the problem are the maps c F(u,u ) = argmax { J(u,u ;v,v ):(v,v ) in VxV } c e e e e e c and G(v,v ) = argmin { J(u,u ;v,v ) : (u,u ) in UxU } . c e e e e e c In the documentation below, the term "alternative element" will always refer c to an element given by applying F or G to the current "best" guess. c c c This routine runs the outer loop associated with Proximal-Point Algorithm. c It calls the subroutine DFGM to solve a subproblem of the above form for c which the matrices P , P , Q and Q are positive definite. c e e c It is expected that the method employed within DFGM will be similar in c structure to the Finite-Generation Method. For this reason, the data c structures have been defined so as to keep track of "finite sets" of primal c and dual vectors. The structure of the algorithm is as follows: c c Given: an initial guess for the vectors u, u , v, v , and stopping c e e c criteria "deps" (double precision floating point number) c and a positive integer "mumax". c Step 0. (a) Initialization : compute the function values for the c initial guesses and compute the trajectories x and y c of the initial guesses. c (b) Set mu=1 and delta=1/10. c Decision Step: If the normalized duality gap is less than 'DEPS', c then stop: the pair (u,ue) and the pair (v,ve) are both c EPSILON-optimal. Otherwise, if mu>>>>>> INITIALIZATION <<<<<<<< C TEPS = 1.D-10 TUDIM = UDIM*TIME + UEDIM TVDIM = VDIM*TIME + VEDIM BESTU = 1 BESTV = 1 ALTU = 2 ALTV = 2 ZBESTU = 9 ZALTU = 10 ZGRDU = 11 ZCGU = 12 ZBESTV = 9 ZALTV = 10 ZGRDV = 11 ZCGV = 12 C C Call routine to calculate the trajectories X and Y for U , V C CALL TRAJ(U(1,1,BESTU),UE(1,BESTU),V(1,1,BESTV),VE(1,BESTV), : X(1,1,BESTU),XE(1,BESTU),Y(1,1,BESTV),YE(1,BESTV), : AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C C Compute the function value ALPHA for U = BESTU and BETA for V = BESTV. C CALL FGVAL(U(1,1,ALTU),V(1,1,ALTV),UE(1,ALTU),VE(1,ALTV), $ ALPHA(BESTU),BETA(BESTV), $ U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV), $ XE(1,BESTU),YE(1,BESTV), $ UUPR,ULOW,VLOW,VUPR,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C Compute the Duality Gap for BESTU,BESTV. C EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF OLDEPS = EPS C CALL OUTPT0(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA,EPS,EPSNRM, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) C C Number Of Iterations: MU = 0 cccc RESTART = 0 RESTART = -1 C C Proximation parameter: cccc DELTA = 1.D0 DELTA = 1.D-1 C C//////////////////////////// C>>>>>>> Begin Loop <<<<<<<< C>>>>>>> ========== <<<<<<<< C//////////////////////////// C 1000 CONTINUE C IF (( EPSNRM .GT. DEPS ) .AND. ( MU .LT. MUMAX )) THEN C C i.e. While ( EPS > DEPS and MU < MUMAX ) C ===== --- ---- === -- ----- C Do the following: C == === ========= DO 40 K=1,TIME DO 10 J=1,UDIM 10 U(J,K,ZBESTU) = U(J,K,BESTU) DO 20 I=1,VDIM 20 V(I,K,ZBESTV) = V(I,K,BESTV) DO 30 J=1,XYDIM X(J,K,ZBESTU) = X(J,K,BESTU) 30 Y(J,K,ZBESTV) = Y(J,K,BESTV) 40 CONTINUE DO 50 J=1,UEDIM 50 UE(J,ZBESTU) = UE(J,BESTU) DO 60 I=1,VEDIM 60 VE(I,ZBESTV) = VE(I,BESTV) DO 70 J=1,XYDIM XE(J,ZBESTU) = XE(J,BESTU) YE(J,ZBESTV) = YE(J,BESTV) 70 CONTINUE C MU = MU + 1 RESTART = RESTART + 1 C CALL OUTPT1(MU,DELTA) C DO 80 J=1,UDIM 80 PD(J) = DELTA+PMAT(J) DO 90 I=1,VDIM 90 QD(I) = DELTA+QMAT(I) DO 120 K=1,TIME DO 100 J=1,UDIM PV(J,K) = PVEC(J,K)-DELTA*U(J,K,BESTU) 100 CONTINUE DO 110 I=1,VDIM QV(I,K) = QVEC(I,K)+DELTA*V(I,K,BESTV) 110 CONTINUE 120 CONTINUE DO 130 J=1,UEDIM PED(J) = DELTA+PEMAT(J) PEV(J) = PEVEC(J)-DELTA*UE(J,BESTU) 130 CONTINUE DO 140 I=1,VEDIM QED(I) = DELTA+QEMAT(I) QEV(I) = QEVEC(I)+DELTA*VE(I,BESTV) 140 CONTINUE C C IF ( BESTU .GT. 1 ) THEN DO 170 K=1,TIME DO 150 J=1,UDIM 150 U(J,K,1) = U(J,K,BESTU) DO 160 J=1,XYDIM 160 X(J,K,1) = X(J,K,BESTU) 170 CONTINUE DO 180 J=1,UEDIM 180 UE(J,1) = UE(J,BESTU) DO 190 J=1,XYDIM 190 XE(J,1) = XE(J,BESTU) BESTU = 1 ENDIF C IF ( BESTV .GT. 1 ) THEN DO 220 K=1,TIME DO 200 I=1,VDIM 200 V(I,K,1) = V(I,K,BESTV) DO 210 I=1,XYDIM 210 Y(I,K,1) = Y(I,K,BESTV) 220 CONTINUE DO 230 I=1,VEDIM 230 VE(I,1) = VE(I,BESTV) DO 240 I=1,XYDIM 240 YE(I,1) = YE(I,BESTV) BESTV = 1 ENDIF C C CALL DFGM(U,V,X,Y,UE,VE,XE,YE,UUPR,ULOW,UEUPR,UELOW, $ VUPR,VLOW,VEUPR,VELOW,PD,PV,QD,QV, $ DMAT,PED,PEV,QED,QEV,AMAT,BMAT,BVEC, $ CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC,XI,ETA, $ ALPHA,BETA,BESTU,BESTV,EPS,EPSNRM,INFORM,NU, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ 8,8,NUMAX,DEPS/1.d8) C C CALL OUTPT2(INFORM,NU) C C DO 940 K=1,TIME DO 910 J=1,UDIM 910 U(J,K,ZGRDU) = U(J,K,BESTU) DO 920 I=1,VDIM 920 V(I,K,ZGRDV) = V(I,K,BESTV) DO 930 J=1,XYDIM X(J,K,ZGRDU) = X(J,K,BESTU) 930 Y(J,K,ZGRDV) = Y(J,K,BESTV) 940 CONTINUE DO 950 J=1,UEDIM 950 UE(J,ZGRDU) = UE(J,BESTU) DO 960 I=1,VEDIM 960 VE(I,ZGRDV) = VE(I,BESTV) DO 970 J=1,XYDIM XE(J,ZGRDU) = XE(J,BESTU) YE(J,ZGRDV) = YE(J,BESTV) 970 CONTINUE C C Compute the pointers for the U and V finite sets. cccc IF ( RESTART .NE. 1 ) THEN IF ( RESTART .GT. 1 ) THEN T1 = 0.D0 T2 = 0.D0 T3 = 0.D0 T4 = 0.D0 DO 340 K=1,TIME DO 310 J=1,UDIM TX = U(J,K,ZBESTU) - U(J,K,ZOLDU) TY = U(J,K,ZOGRDU) - U(J,K,ZGRDU) + TX IF ( U(J,K,ZGRDU).LT.UUPR(J)-TEPS .AND. $ U(J,K,ZGRDU).GT.ULOW(J)+TEPS ) THEN TG = U(J,K,ZGRDU) - U(J,K,ZBESTU) T1 = T1 + TX*TG T2 = T2 + (TX - TY)*TG ENDIF T3 = T3 + TX*TY T4 = T4 + TY*(TX - TY) 310 CONTINUE DO 320 I=1,VDIM TX = V(I,K,ZBESTV) - V(I,K,ZOLDV) TY = V(I,K,ZGRDV) - V(I,K,ZOGRDV) - TX IF ( V(I,K,ZGRDV).LT.VUPR(I)-TEPS .AND. $ V(I,K,ZGRDV).GT.VLOW(I)+TEPS ) THEN TG = V(I,K,ZGRDV) - V(I,K,ZBESTV) T1 = T1 - TX*TG T2 = T2 - (TX + TY)*TG ENDIF T3 = T3 + TX*TY T4 = T4 + TY*(TX + TY) 320 CONTINUE 340 CONTINUE DO 350 J=1,UEDIM TX = UE(J,ZBESTU) - UE(J,ZOLDU) TY = UE(J,ZOGRDU) - UE(J,ZGRDU) + TX IF ( UE(J,ZGRDU).LT.UEUPR(J)-TEPS .AND. $ UE(J,ZGRDU).GT.UELOW(J)+TEPS ) THEN TG = UE(J,ZGRDU) - UE(J,ZBESTU) T1 = T1 + TX*TG T2 = T2 + (TX - TY)*TG ENDIF T3 = T3 + TX*TY T4 = T4 + TY*(TX - TY) 350 CONTINUE DO 360 I=1,VEDIM TX = VE(I,ZBESTV) - VE(I,ZOLDV) TY = VE(I,ZGRDV) - VE(I,ZOGRDV) - TX IF ( VE(I,ZGRDV).LT.VEUPR(I)-TEPS .AND. $ VE(I,ZGRDV).GT.VELOW(I)+TEPS ) THEN TG = VE(I,ZGRDV) - VE(I,ZBESTV) T1 = T1 - TX*TG T2 = T2 - (TX + TY)*TG ENDIF T3 = T3 + TX*TY T4 = T4 + TY*(TX + TY) 360 CONTINUE C cccc Preventing something divided by zero! IF(DABS(T3).LT.TEPS) GOTO 4444 C C Compute ZCGu C DO 4 J = 1, UDIM DO 4 K = 1, TIME IF ( U(J,K,ZGRDU).LT.UUPR(J)-TEPS .AND. $ U(J,K,ZGRDU).GT.ULOW(J)+TEPS ) THEN TX = U(J,K,ZBESTU) - U(J,K,ZOLDU) TY = U(J,K,ZOGRDU) - U(J,K,ZGRDU) + TX U(J,K,ZCGU) = U(J,K,ZGRDU) - (TX*T4*T1)/(T3*T3) + $ ((TX - TY)*T1 + TX*T2)/T3 ccc print*, 'inprimal','j=',j,'k=',k ELSE U(J,K,ZCGU) = U(J,K,ZGRDU) ENDIF 4 CONTINUE DO 3 J = 1, UEDIM IF ( UE(J,ZGRDU).LT.UEUPR(J)-TEPS .AND. $ UE(J,ZGRDU).GT.UELOW(J)+TEPS ) THEN TX = UE(J,ZBESTU) - UE(J,ZOLDU) TY = UE(J,ZOGRDU) - UE(J,ZGRDU) + TX UE(J,ZCGU) = UE(J,ZGRDU) - (TX*T4*T1)/(T3*T3) + $ ((TX - TY)*T1 + TX*T2)/T3 ccc print*, 'inprimal','j=',j ELSE UE(J,ZCGU) = UE(J,ZGRDU) ENDIF 3 CONTINUE C C C C Compute ZCGv C DO 24 I = 1, VDIM DO 24 K = 1, TIME IF ( V(I,K,ZGRDV).LT.VUPR(I)-TEPS .AND. $ V(I,K,ZGRDV).GT.VLOW(I)+TEPS ) THEN TX = V(I,K,ZBESTV) - V(I,K,ZOLDV) TY = V(I,K,ZGRDV) - V(I,K,ZOGRDV) - TX V(I,K,ZCGV) = V(I,K,ZGRDV) - (TX*T4*T1)/(T3*T3) + $ ((TX - TY)*T1 + TX*T2)/T3 ccc print*, 'indual','i=',i,'k=',k ELSE V(I,K,ZCGV) = V(I,K,ZGRDV) ENDIF 24 CONTINUE DO 23 I = 1, VEDIM IF ( VE(I,ZGRDV).LT.VEUPR(I)-TEPS .AND. $ VE(I,ZGRDV).GT.VELOW(I)+TEPS ) THEN TX = VE(I,ZBESTV) - VE(I,ZOLDV) TY = VE(I,ZGRDV) - VE(I,ZOGRDV) - TX VE(I,ZCGV) = VE(I,ZGRDV) - (TX*T4*T1)/(T3*T3) + $ ((TX - TY)*T1 + TX*T2)/T3 ccc print*, 'indual','i=',i ELSE VE(I,ZCGV) = VE(I,ZGRDV) ENDIF 23 CONTINUE C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C cccc DUM = 1.d0/2.d0**(RESTART - 4) DUM = 1.d2/2.d0**(RESTART - 4) DO 1170 K=1,TIME DO 1150 J=1,UDIM IF( U(J,K,ZGRDU)-U(J,K,ZCGU) .GT. TEPS) : DUM = DMIN1((U(J,K,ZGRDU)-ULOW(J))/ : (U(J,K,ZGRDU)-U(J,K,ZCGU)), DUM) 1150 IF( U(J,K,ZCGU)-U(J,K,ZGRDU) .GT. TEPS) : DUM = DMIN1((UUPR(J)-U(J,K,ZGRDU))/ : (U(J,K,ZCGU)-U(J,K,ZGRDU)), DUM) 1170 CONTINUE DO 1180 J=1,UEDIM IF( UE(J,ZGRDU)-UE(J,ZCGU) .GT. TEPS) : DUM = DMIN1((UE(J,ZGRDU)-UELOW(J))/ : (UE(J,ZGRDU)-UE(J,ZCGU)), DUM) 1180 IF( UE(J,ZCGU)-UE(J,ZGRDU) .GT. TEPS) : DUM = DMIN1((UEUPR(J)-UE(J,ZGRDU))/ : (UE(J,ZCGU)-UE(J,ZGRDU)), DUM) C ccc print*, 'the primal dum =',dum DO 3170 K=1,TIME DO 3150 J=1,UDIM 3150 U(J,K,1) = U(J,K,ZGRDU) + DUM * : (U(J,K,ZCGU)-U(J,K,ZGRDU)) 3170 CONTINUE DO 3180 J=1,UEDIM 3180 UE(J,1) = UE(J,ZGRDU) + DUM * : (UE(J,ZCGU)-UE(J,ZGRDU)) C cccc DUM = 1.d0/2.d0**(RESTART - 4) DUM = 1.d2/2.d0**(RESTART - 4) DO 2170 K=1,TIME DO 2150 I=1,VDIM IF( V(I,K,ZGRDV)-V(I,K,ZCGV) .GT. TEPS) : DUM = DMIN1((V(I,K,ZGRDV)-VLOW(I))/ : (V(I,K,ZGRDV)-V(I,K,ZCGV)), DUM) 2150 IF( V(I,K,ZCGV)-V(I,K,ZGRDV) .GT. TEPS) : DUM = DMIN1((VUPR(I)-V(I,K,ZGRDV))/ : (V(I,K,ZCGV)-V(I,K,ZGRDV)), DUM) 2170 CONTINUE DO 2180 I=1,VEDIM IF( VE(I,ZGRDV)-VE(I,ZCGV) .GT. TEPS) : DUM = DMIN1((VE(I,ZGRDV)-VELOW(I))/ : (VE(I,ZGRDV)-VE(I,ZCGV)), DUM) 2180 IF( VE(I,ZCGV)-VE(I,ZGRDV) .GT. TEPS) : DUM = DMIN1((VEUPR(I)-VE(I,ZGRDV))/ : (VE(I,ZCGV)-VE(I,ZGRDV)), DUM) C ccc print*, 'the dual dum =',dum DO 3270 K=1,TIME DO 3250 I=1,VDIM 3250 V(I,K,1) = V(I,K,ZGRDV) + DUM * : (V(I,K,ZCGV)-V(I,K,ZGRDV)) 3270 CONTINUE DO 3280 I=1,VEDIM 3280 VE(I,1) = VE(I,ZGRDV) + DUM * : (VE(I,ZCGV)-VE(I,ZGRDV)) C CALL TRAJ(U(1,1,1),UE(1,1),V(1,1,1),VE(1,1), : X(1,1,1),XE(1,1),Y(1,1,1),YE(1,1), : AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C EPSOLD = OLDEPS C BESTU = 1 BESTV = 1 ALTU = 2 ALTV = 2 C CALL LINE(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ ZGRDU,ZGRDV,1,1,BESTU,BESTV,ALTU,ALTV, $ UUPR,ULOW,VUPR,VLOW,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT, $ PEMAT,PEVEC,QEMAT,QEVEC, $ AMAT,BMAT,BVEC,CMAT,CVEC, $ BEMAT,BEVEC,CEMAT,CEVEC, $ MAXUFE,MAXVFE, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME, $ WORK,LENWRK,QUIT) C C ENDIF 4444 CONTINUE IF(RESTART.LE.1 .OR. DABS(T3).LT.TEPS) THEN C EPSOLD = OLDEPS C C C Compute the function value ALPHA for U = BESTU and BETA for V = BESTV. C CALL FGVAL(U(1,1,ALTU),V(1,1,ALTV),UE(1,ALTU),VE(1,ALTV), $ ALPHA(BESTU),BETA(BESTV), $ U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV), $ XE(1,BESTU),YE(1,BESTV), $ UUPR,ULOW,VLOW,VUPR,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C C Compute the Duality Gap for BESTU,BESTV. C RATIO = (ALPHA(BESTU) - BETA(BESTV))/OLDEPS EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF C USTEP = 0. VSTEP = 0. C CALL OUTLNS(U(1,1,BESTU),V(1,1,BESTV),X(1,1,BESTU),Y(1,1,BESTV), $ UE(1,BESTU),VE(1,BESTV),XE(1,BESTU),YE(1,BESTV), $ ALPHA(BESTU),BETA(BESTV),EPS,EPSNRM,RATIO,USTEP,VSTEP, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME) C C ENDIF C Compute the Duality Gap for BESTU,BESTV. C EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF C OLDEPS = EPS C C Update DELTA if convergence was quick enough. C cccc IF (INFORM .EQ. 0 .AND. cccc $ NU .LE. MAX(TIME*UDIM+UEDIM,TIME*VDIM+VEDIM)) THEN cccc DELTA = DELTA/10.0d0 cccccccc RESTART = 0 cccc RESTART = -1 cccc ENDIF C ccccccc IF(RESTART.EQ.1 .AND. RATIO.LT..2D0) RESTART = 0 C C ZOLDU = ZBESTU ZOALTU = ZALTU ZOGRDU = ZGRDU ZOCGU = ZCGU IF ( ZOLDU .EQ. 9 ) THEN ZBESTU = 13 ZALTU = 14 ZGRDU = 15 ZCGU = 16 ELSE ZBESTU = 9 ZALTU = 10 ZGRDU = 11 ZCGU = 12 ENDIF C ZOLDV = ZBESTV ZOALTV = ZALTV ZOGRDV = ZGRDV ZOCGV = ZCGV IF ( ZOLDV .EQ. 9 ) THEN ZBESTV = 13 ZALTV = 14 ZGRDV = 15 ZCGV = 16 ELSE ZBESTV = 9 ZALTV = 10 ZGRDV = 11 ZCGV = 12 ENDIF C C Continue iterating ! GO TO 1000 ENDIF C////////////////////////// C>>>>>>> End Loop <<<<<<<< C>>>>>>> ======== <<<<<<<< C////////////////////////// C C C IF ( EPSNRM .LE. DEPS ) THEN C Successful termination ! INFORM = 0 ELSE C Maximum number of iterations exceeded ! INFORM = 1 ENDIF C CALL OUTPT4(MU,INFORM,EPS,EPSNRM,ALPHA,BETA, : U,V,UE,VE,DVAL,UDIM,UEDIM,VDIM,VEDIM,TIME) C======================================================================= C RETURN C End of subroutine 'DFGM' END :::::::::::::: psearch.f :::::::::::::: c The following code is a modified version of Stephen E. Wright's c original code. subroutine search(ct,atv,atm,bem,bm,btil,dtil,betil,detil, # zlft,zrt,zmid,brk1,brk2,slope, # zelft,zert,zemid,ebrk1,ebrk2,eslope, # lambda,z,ze,fval,zdim,zedim,ntime) c c INPUT: integer zdim,zedim,ntime double precision ct,atv,atm double precision bm(zdim),bem(zedim) double precision btil(zdim,ntime),dtil(zdim,ntime) double precision betil(zedim),detil(zedim) double precision zlft(zdim,ntime),zrt(zdim,ntime) double precision zmid(zdim,ntime) double precision brk1(zdim,ntime),brk2(zdim,ntime) double precision slope(zdim,ntime) double precision zelft(zedim),zert(zedim),zemid(zedim) double precision ebrk1(zedim),ebrk2(zedim),eslope(zedim) c c OUTPUT: double precision lambda,z(zdim,ntime),ze(zedim),fval c c dummy variables and counters: double precision a,b,fplus,fminus,tt,tfplus,tfminus double precision sum1,eps,ztemp integer n,i,k,ik,jk c c*********************************************************************** c This routine solves the problem c minimize f(lambda) = c c ct + atv*lambda + (1/2)lambda*lambda*atm c + sum / (btil(t)+lambda*dtil(t))*z(lambda;t) \ c t=1,N \ - (1/2)z(lambda;t)*bm *z(lambda;t) / c + { (betil+lambda*detil)*ze(lambda) - (1/2)ze(lambda)*bem*ze(lambda) } c c subject to 0.0<= lambda <=1.0 , where c / zlft(i,t) if lambda <= brk1(i,t) c z(lambda;i,t) = < zrt(i,t) if lambda > brk2(i,t) c \ zmid(i,t) + lambda*slope(i,t) otherwise. c and c / zelft(i) if lambda <= ebrk1(i) c ze(lambda;i) = < zert(i) if lambda > ebrk2(i) c \ zemid(i) + lambda*eslope(i) otherwise. c c It is assumed that the data were prepared by 'LNPREP' so that, for example, c z and ze are continuous if brk1 < brk2. This also guarantees certain c relationships among the data that enable us to evaluate the derivatives c with fewer arithmetical operations. c c First the right-hand derivative at 0 and the left-hand derivative at 1 are c computed. If the former is positive we exit with lambda=0.0; if the latter c is negative then we take lambda=1.0. Then check the values of the c right-hand and left-hand derivatives at the break points to reduce the c interval of uncertainty from [0,1] to [a,b] where there is no break point c inbetween. And finally, use an interpolation to get the optimal c lambda. c c The routine also returns the function value 'fval' as well as the vectors c z(lambda), ze(lambda). c c The following code was written by Ciyou Zhu in July 1989 c at University of Washington. c======================================================================= c EPS is a distinguishability constant. eps = 1.0d-11 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c ====== c SEARCH c ====== c a = 0.0d0 b = 1.0d0 c c COMPUTE the RH-derivative at a, LH-derivative at b fplus = atv fminus = atv + atm do 11 j=1,ntime do 10 i=1,zdim if ( 0.0d0 .lt. brk1(i,j) ) then fplus = fplus + dtil(i,j)*zlft(i,j) else if ( 0.0d0 .ge. brk2(i,j) ) then fplus = fplus + dtil(i,j)*zrt(i,j) else fplus = fplus + dtil(i,j)*zmid(i,j) endif if ( 1.0d0 .le. brk1(i,j) ) then fminus = fminus + dtil(i,j)*zlft(i,j) else if ( 1.0d0 .gt. brk2(i,j) ) then fminus = fminus + dtil(i,j)*zrt(i,j) else fminus = fminus + dtil(i,j)*(zmid(i,j) + slope(i,j)) endif 10 continue 11 continue do 12 i=1,zedim if ( 0.0d0 .lt. ebrk1(i) ) then fplus = fplus + detil(i)*zelft(i) else if ( 0.0d0 .ge. ebrk2(i) ) then fplus = fplus + detil(i)*zert(i) else fplus = fplus + detil(i)*zemid(i) endif if ( 1.0d0 .le. ebrk1(i) ) then fminus = fminus + detil(i)*zelft(i) else if ( 1.0d0 .gt. ebrk2(i) ) then fminus = fminus + detil(i)*zert(i) else fminus = fminus + detil(i)*(zemid(i) + eslope(i)) endif 12 continue c c if ( fplus .gt. -eps ) then c !! Stop with lambda=0.0d0; compute LH-derivative at lambda print*,' lambda=0 is optimal !' lambda = 0.0d0 fminus = atv do 150 j=1,ntime do 150 i=1,zdim if ( 0.0d0 .le. brk1(i,j)) then fminus = fminus + zlft(i,j)*dtil(i,j) else if ( 0.0d0 .gt. brk2(i,j) ) then fminus = fminus + zrt(i,j)*dtil(i,j) else fminus= fminus + zmid(i,j)*dtil(i,j) endif 150 continue do 151 i=1,zedim if ( 0.0d0 .le. ebrk1(i)) then fminus = fminus + zelft(i)*detil(i) else if ( 0.0d0 .gt. ebrk2(i) ) then fminus = fminus + zert(i)*detil(i) else fminus= fminus + zemid(i)*detil(i) endif 151 continue c else if ( fminus .lt. eps ) then c !! Stop with lambda=1.0d0; compute RH-derivative at lambda print*,' lambda=1 is optimal !' lambda = 1.0d0 fplus = atv + atm do 110 j=1,ntime do 110 i=1,zdim if ( 1.0d0 .lt. brk1(i,j) ) then fplus = fplus + zlft(i,j)*dtil(i,j) else if ( 1.0d0 .ge. brk2(i,j) ) then fplus = fplus + zrt(i,j)*dtil(i,j) else fplus = fplus + (zmid(i,j)+slope(i,j))*dtil(i,j) endif 110 continue do 111 i=1,zedim if ( 1.0d0 .lt. ebrk1(i) ) then fplus = fplus + zelft(i)*detil(i) else if ( 1.0d0 .ge. ebrk2(i) ) then fplus = fplus + zert(i)*detil(i) else fplus = fplus + (zemid(i)+eslope(i))*detil(i) endif 111 continue else c !! Search inside the interval [0,1] for fplus=0.0d0 c c Check the break points n = 0 do 500 k=1,2 do 500 jk=1,ntime do 500 ik=1,zdim if ( k .eq. 1 ) tt = brk1(ik,jk) if ( k .eq. 2 ) tt = brk2(ik,jk) if ( a .lt. tt .and. tt .lt. b ) then c !!COMPUTE the right-hand and left-hand derivatives at tt n = n + 1 tfplus = atv + tt*atm tfminus = atv + tt*atm do 280 j=1,ntime do 280 i=1,zdim if ( tt.lt. brk1(i,j) ) then tfplus = tfplus + zlft(i,j)*dtil(i,j) else if ( tt.ge. brk2(i,j) ) then tfplus = tfplus + zrt(i,j)*dtil(i,j) else tfplus = tfplus + (zmid(i,j)+slope(i,j)*tt)*dtil(i,j) endif if ( tt.le. brk1(i,j) ) then tfminus = tfminus + zlft(i,j)*dtil(i,j) else if ( tt.gt. brk2(i,j) ) then tfminus = tfminus + zrt(i,j)*dtil(i,j) else tfminus = tfminus+(zmid(i,j)+slope(i,j)*tt)*dtil(i,j) endif 280 continue do 281 i=1,zedim if ( tt.lt. ebrk1(i) ) then tfplus = tfplus + zelft(i)*detil(i) else if ( tt.ge. ebrk2(i) ) then tfplus = tfplus + zert(i)*detil(i) else tfplus = tfplus + (zemid(i)+eslope(i)*tt)*detil(i) endif if ( tt.le. ebrk1(i) ) then tfminus = tfminus + zelft(i)*detil(i) else if ( tt.gt. ebrk2(i) ) then tfminus = tfminus + zert(i)*detil(i) else tfminus = tfminus + (zemid(i)+eslope(i)*tt)*detil(i) endif 281 continue c c Reduce the interval if ( tfplus .lt. -eps ) then a = tt fplus = tfplus else if (tfminus .gt. eps ) then b = tt fminus = tfminus else lambda = tt fplus = tfplus fminus = tfminus print*,' one of the break points is optimal !' goto 82 endif endif 500 continue do 501 k=1,2 do 501 ik=1,zedim if ( k .eq. 1 ) tt = ebrk1(ik) if ( k .eq. 2 ) tt = ebrk2(ik) if ( a .lt. tt .and. tt .lt. b ) then c !!COMPUTE the right-hand and left-hand derivatives at tt n = n + 1 tfplus = atv + tt*atm tfminus = atv + tt*atm do 480 j=1,ntime do 480 i=1,zdim if ( tt.lt. brk1(i,j) ) then tfplus = tfplus + zlft(i,j)*dtil(i,j) else if ( tt.ge. brk2(i,j) ) then tfplus = tfplus + zrt(i,j)*dtil(i,j) else tfplus = tfplus + (zmid(i,j)+slope(i,j)*tt)*dtil(i,j) endif if ( tt.le. brk1(i,j) ) then tfminus = tfminus + zlft(i,j)*dtil(i,j) else if ( tt.gt. brk2(i,j) ) then tfminus = tfminus + zrt(i,j)*dtil(i,j) else tfminus = tfminus+(zmid(i,j)+slope(i,j)*tt)*dtil(i,j) endif 480 continue do 481 i=1,zedim if ( tt.lt. ebrk1(i) ) then tfplus = tfplus + zelft(i)*detil(i) else if ( tt.ge. ebrk2(i) ) then tfplus = tfplus + zert(i)*detil(i) else tfplus = tfplus + (zemid(i)+eslope(i)*tt)*detil(i) endif if ( tt.le. ebrk1(i) ) then tfminus = tfminus + zelft(i)*detil(i) else if ( tt.gt. ebrk2(i) ) then tfminus = tfminus + zert(i)*detil(i) else tfminus = tfminus + (zemid(i)+eslope(i)*tt)*detil(i) endif 481 continue c c Reduce the interval if ( tfplus .lt. -eps ) then a = tt fplus = tfplus else if (tfminus .gt. eps ) then b = tt fminus = tfminus else lambda = tt fplus = tfplus fminus = tfminus print*,' one of the break points is optimal !' goto 82 endif endif 501 continue lambda = ( a * fminus - b * fplus ) / ( fminus - fplus ) c c COMPUTE the right-hand and left-hand derivatives at lambda fplus = atv + lambda*atm fminus = atv + lambda*atm do 80 j=1,ntime do 80 i=1,zdim if ( lambda .lt. brk1(i,j) ) then fplus = fplus + zlft(i,j)*dtil(i,j) else if ( lambda .ge. brk2(i,j) ) then fplus = fplus + zrt(i,j)*dtil(i,j) else fplus = fplus + (zmid(i,j)+slope(i,j)*lambda)*dtil(i,j) endif if ( lambda .le. brk1(i,j) ) then fminus = fminus + zlft(i,j)*dtil(i,j) else if ( lambda .gt. brk2(i,j) ) then fminus = fminus + zrt(i,j)*dtil(i,j) else fminus = fminus+(zmid(i,j)+slope(i,j)*lambda)*dtil(i,j) endif 80 continue do 81 i=1,zedim if ( lambda .lt. ebrk1(i) ) then fplus = fplus + zelft(i)*detil(i) else if ( lambda .ge. ebrk2(i) ) then fplus = fplus + zert(i)*detil(i) else fplus = fplus + (zemid(i)+eslope(i)*lambda)*detil(i) endif if ( lambda .le. ebrk1(i) ) then fminus = fminus + zelft(i)*detil(i) else if ( lambda .gt. ebrk2(i) ) then fminus = fminus + zert(i)*detil(i) else fminus = fminus + (zemid(i)+eslope(i)*lambda)*detil(i) endif 81 continue print*,' ',n,' break points are checked !' endif 82 continue c c End of Search c ===-==-====== c c c Calculate z(lambda) and the optimal value fval. fval = ct + atv*lambda + lambda*lambda*atm/2.0d0 do 300 j=1,ntime do 300 i=1,zdim if ( lambda .lt. brk1(i,j) ) then ztemp = zlft(i,j) else if ( lambda .gt. brk2(i,j) ) then ztemp = zrt(i,j) else if ( (fminus .lt. -eps) .and. (lambda .eq. brk1(i,j)) # .and. ( brk1(i,j) .eq. brk2(i,j)) ) then tt = fminus/(fminus-fplus) ztemp = (1.0d0-tt)*zlft(i,j) + tt*zrt(i,j) else ztemp = zmid(i,j) + slope(i,j)*lambda endif z(i,j) = ztemp sum1 = btil(i,j) + dtil(i,j)*lambda - ztemp*bm(i)/2.0d0 fval = fval + sum1*ztemp 300 continue do 301 i=1,zedim if ( lambda .lt. ebrk1(i) ) then ztemp = zelft(i) else if ( lambda .gt. ebrk2(i) ) then ztemp = zert(i) else if ( (fminus .lt. -eps) .and. (lambda .eq. ebrk1(i)) # .and. ( ebrk1(i) .eq. ebrk2(i)) ) then tt = fminus/(fminus-fplus) ztemp = (1.0d0-tt)*zelft(i) + tt*zert(i) else ztemp = zemid(i) + eslope(i)*lambda endif ze(i) = ztemp sum1 = betil(i) + detil(i)*lambda - ztemp*bem(i)/2.0d0 fval = fval + sum1*ztemp 301 continue c c c c END OF SUBROUTINE search return end :::::::::::::: puvsets.f :::::::::::::: c The following code is a modified version of Stephen E. Wright's c original code. SUBROUTINE UVSETS(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA,XI,ETA, $ UUPR,ULOW,VUPR,VLOW,UEUPR,UELOW,VEUPR,VELOW, $ PMAT,PVEC,QMAT,QVEC,DMAT,PEMAT,PEVEC,QEMAT,QEVEC, $ AMAT,BMAT,BVEC,CMAT,CVEC,BEMAT,BEVEC,CEMAT,CEVEC, $ BESTU,BESTV,ALTU,ALTV,UDIR,VDIR,OLDU,OLDV, $ GRDU,GRDV,CGU,CGV,UTIME,VTIME, $ OALTU,OALTV,OGRDU,OGRDV,OCGU,OCGV,OUDIR,OVDIR, $ NUFE,NVFE,NU,MAXUFE,MAXVFE, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME,WORK,LENWRK,QUIT) C*********************************************************************** C C This routine is required by the DYNFGM program to do the following: C a. Build the finite-element sets (within each iteration) and C compute the corresponding trajectories. C b. Return (i) the function values for various elements of the C U and V arrays; C (ii) the total number of elements currently in the C U and V sets; C (iii) the positions in the U and V arrays in which C the direction vectors (to be calculated in the C routine 'DIRECT') should be put; C (iv) the positions in the U and V arrays of the (old) C "best" elements of the previous iteration; C (v) the positions in the U and V arrays in which C the NEW "best" elements and their "alternative" C elements (to be calculated in the routine 'LINE') C should be placed. C c. Initialize the arrays XI,ETA of "barycentric coordinates". C d. Print any output which is deemed necessary at this point in the C algorithm. C*********************************************************************** C DECLARATIONS FOR VARIABLES PASSED THRU THE SUBROUTINE CALL: C C Variables whose values may be changed within UVSETS: C U,V,UE,VE -- finite sets of decision vectors; C X,Y,XE,YE -- the corresponding trajectories; C XI,ETA -- vectors of coefficients giving a convex combination of C the elements in the finite sets; C ALPHA,BETA -- arrays giving the function values associated with C the elements of the finite sets; C NUFE,NVFE -- Integers indicating the size of the current set of C U or V elements (which are supposedly stored in the C first NUFE or NVFE positions); C QUIT -- Flag whose value is .FALSE. when UVSETS is called. If it C is decided that the algorithm should be immediately C terminated (due to some criterion specified within C UVSETS by the user), then UVSETS should change QUIT to C be .TRUE. and print out any appropriate information. INTEGER UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME,MAXUFE,MAXVFE DOUBLE PRECISION U(UDIM,TIME,MAXUFE),V(VDIM,TIME,MAXVFE) DOUBLE PRECISION X(XYDIM,TIME,MAXUFE),Y(XYDIM,TIME,MAXVFE) DOUBLE PRECISION UE(UEDIM,MAXUFE),VE(VEDIM,MAXVFE) DOUBLE PRECISION XE(XYDIM,MAXUFE),YE(XYDIM,MAXVFE) DOUBLE PRECISION XI(MAXUFE),ETA(MAXVFE) DOUBLE PRECISION ALPHA(MAXUFE),BETA(MAXVFE) INTEGER NUFE,NVFE LOGICAL QUIT C C Pointers (indicating positions in the finite sets) C 1. UDIR,VDIR. On entering UVSETS, UDIR indicates the position of C the direction U-vector found in the call to DIRECT during C the previous iteration (unless NU=1, in which case it is 0). C Upon exiting, UDIR should indicate the desired position in C which the new direction vector should be placed during the C next call to the routine DIRECT; VDIR is used similarly. C 2. BESTU,BESTV. On entering UVSETS, BESTU indicates the position of C the "best" U-vector found in the call to LINE during the C previous iteration (unless NU=1, in which case it indicates C the initial "guess"). The vector in this position will be C used in the current iteration, so it must be saved in the C U array.(The position to which it is moved must be returned C by UVSETS in the pointer OLDU.) Upon exiting, BESTU should C indicate the position in which the new "best" vector will C be placed during the next call to the routine LINE. The C description for BESTV is similar. C 3. ALTU,ALTV. On entering UVSETS, ALTU indicates the position of C the "alternative" U-vector which gives the minimum defining C the function value for the previous "best" V-vector. Upon C exiting, ALTU should indicate the position in which the new C "alternate" U-vector will be placed during the next call to C the routine LINE. Similarly for ALTV. C 4. OLDU,OLDV. When UVSETS is entered, the current "best" U-vector C (as calculated in the previous call to LINE) is indicated by C the pointer BESTU. This vector must be saved in the U array C because it will be used in the next call to LINE. The position C in which it may be found after exiting from UVSETS is passed in C the variable OLDU. Similarly for OLDV. INTEGER BESTU,BESTV INTEGER ALTU,ALTV INTEGER UDIR,VDIR INTEGER OLDU,OLDV INTEGER GRDU,GRDV,CGU,CGV,UTIME,VTIME INTEGER OALTU,OALTV,OGRDU,OGRDV,OCGU,OCGV,OUDIR,OVDIR C C Data (values not to be changed in UVSETS) DOUBLE PRECISION UUPR(UDIM),ULOW(UDIM) DOUBLE PRECISION VUPR(VDIM),VLOW(VDIM) DOUBLE PRECISION UEUPR(UEDIM),UELOW(UEDIM) DOUBLE PRECISION VEUPR(VEDIM),VELOW(VEDIM) DOUBLE PRECISION PMAT(UDIM),PVEC(UDIM,TIME) DOUBLE PRECISION QMAT(VDIM),QVEC(VDIM,TIME) DOUBLE PRECISION DMAT(VDIM,UDIM) DOUBLE PRECISION PEMAT(UDIM),PEVEC(UDIM),QEMAT(UDIM),QEVEC(UDIM) DOUBLE PRECISION AMAT(XYDIM,XYDIM) DOUBLE PRECISION BMAT(XYDIM,UDIM),BVEC(XYDIM) DOUBLE PRECISION CMAT(VDIM,XYDIM),CVEC(XYDIM) DOUBLE PRECISION BEMAT(XYDIM,UEDIM),BEVEC(XYDIM) DOUBLE PRECISION CEMAT(VEDIM,XYDIM),CEVEC(XYDIM) C C C Dimensions for the above arrays C C Workspace for use in UVSETS and other routines called during the C algorithm (NOTE: values placed in this array cannot be saved between C successive calls to UVSETS.) INTEGER LENWRK DOUBLE PRECISION WORK(LENWRK) C C*********************************************************************** C*********************************************************************** C C ALL CODE FROM THIS POINT ON IS USER-SUPPLIED !! THE ABOVE GIVES THE C STANDARD DESCRIPTION (OR TEMPLATE) OF WHAT IS EXPECTED OF THIS C SUBROUTINE BY THE 'DYNFGM' PROGRAM AND THE NATURE OF THE DATA PASSED C IN AND OUT. THIS ROUTINE MAY COMMUNICATE WITH OTHER USER-SUPPLIED C ROUTINES THROUGH THE USE OF N A M E D COMMON BLOCKS. C --------- C*********************************************************************** C*********************************************************************** C C THE FOLLOWING CODE WAS WRITTEN BY Ciyou Zhu IN JULY 1989 C --- --------- ---- --- ------- -- ================= -- =========== C AT UNIVERSITY OF WASHINGTON C -- ======================== C C Common Blocks: C 1. PRDATA contains "PRint DATA". The numbers PROPT(j,5) (for j=1,2,3) C give the unit numbers for the various output files. Typically, j=1 C corresponds to the standard output unit, i.e. PROPT(1,5)=6. The other C entries of PROPT should represent the various printing options chosen C by the user: it allows for three output units with four integer C descriptions for each. COMMON /PRDATA/ PROPT INTEGER PROPT(3,5) C 2. SETSIZ contains information concerning the nature of the finite sets. C MODE gives the number of elements to be generated recursively from C the elements U(ALTU) and V(ALTV). If KEEPU is true, then all the C U elements of the past are kept. KEEPV has a similar function. COMMON /SETSIZ/ MODE,KEEPU,KEEPV INTEGER MODE LOGICAL KEEPU,KEEPV C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Local variables: INTEGER OUT INTEGER JT,NT,TUDIM,TVDIM INTEGER NBESTU,NBESTV,NALTU,NALTV,NGRDU,NGRDV DOUBLE PRECISION TX,TY,TG,T1,T2,T3,TEPS C C TEPS is the constant used to prevent overflowing when C calculating the CG points. C TEPS = 1.D-10 C TUDIM = UDIM*TIME + UEDIM TVDIM = VDIM*TIME + VEDIM C C----------------------------------------------------------------------- C DO 191 JJ=1,3 IF ( PROPT(JJ,1) .GE. 2 ) THEN IF ( PROPT(JJ,2) .GE. 2 ) THEN OUT = PROPT(JJ,5) WRITE(OUT,'(80("="))') WRITE(OUT,'(4x,"INNER ITERATION ",I2,":")') NU ENDIF ENDIF 191 CONTINUE C C Compute the Gradient Points: C CALL FGVAL(U(1,1,GRDU),V(1,1,GRDV), : UE(1,GRDU),VE(1,GRDV), : ALPHA(ALTU),BETA(ALTV), : U(1,1,ALTU),V(1,1,ALTV), : X(1,1,ALTU),Y(1,1,ALTV), : UE(1,ALTU),VE(1,ALTV), : XE(1,ALTU),YE(1,ALTV), : UUPR,ULOW,VLOW,VUPR, : UEUPR,UELOW,VEUPR,VELOW, : PMAT,PVEC,QMAT,QVEC,DMAT, : PEMAT,PEVEC,QEMAT,QEVEC, : BMAT,BVEC,CMAT,CVEC, : BEMAT,BEVEC,CEMAT,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C - - - - - - - - - - - - - - - - - - - - - - - C C Compute the Trajectories GRDx , GRDy , for GRDu , GRDv C C CALL TRAJ(U(1,1,GRDU),UE(1,GRDU), : V(1,1,GRDV),VE(1,GRDV), : X(1,1,GRDU),XE(1,GRDU), : Y(1,1,GRDV),YE(1,GRDV), : AMAT,BMAT,CMAT,BVEC,CVEC, : BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C C Check whether to restart by using the dual information C IF (BETA(ALTV).GT.BETA(BESTV).OR.ALPHA(ALTU).LT.ALPHA(BESTU)) : THEN IF (BETA(ALTV).GT.BETA(BESTV)) THEN NBESTV = ALTV NALTU = GRDU NGRDV = BESTV VTIME = NU - 1 PRINT *, ' Dual restarted !' ELSE NBESTV = BESTV NALTU = ALTU NGRDV = GRDV ENDIF IF (ALPHA(ALTU).LT.ALPHA(BESTU)) THEN NBESTU = ALTU NALTV = GRDV NGRDU = BESTU UTIME = NU - 1 PRINT *, ' Primal restarted !' ELSE NBESTU = BESTU NALTV = ALTV NGRDU = GRDU ENDIF C BESTU = NBESTU BESTV = NBESTV ALTU = NALTU ALTV = NALTV GRDU = NGRDU GRDV = NGRDV C CALL FGVAL(U(1,1,GRDU),V(1,1,GRDV), : UE(1,GRDU),VE(1,GRDV), : ALPHA(ALTU),BETA(ALTV), : U(1,1,ALTU),V(1,1,ALTV), : X(1,1,ALTU),Y(1,1,ALTV), : UE(1,ALTU),VE(1,ALTV), : XE(1,ALTU),YE(1,ALTV), : UUPR,ULOW,VLOW,VUPR, : UEUPR,UELOW,VEUPR,VELOW, : PMAT,PVEC,QMAT,QVEC,DMAT, : PEMAT,PEVEC,QEMAT,QEVEC, : BMAT,BVEC,CMAT,CVEC, : BEMAT,BEVEC,CEMAT,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C - - - - - - - - - - - - - - - - - - - - - - - C C Compute the Trajectories GRDx , GRDy , for GRDu , GRDv C C CALL TRAJ(U(1,1,GRDU),UE(1,GRDU), : V(1,1,GRDV),VE(1,GRDV), : X(1,1,GRDU),XE(1,GRDU), : Y(1,1,GRDV),YE(1,GRDV), : AMAT,BMAT,CMAT,BVEC,CVEC, : BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C ENDIF C C Compute the size of the U and V finite sets. NUFE = 8 NVFE = 8 C C C Compute the pointers for the U and V finite sets. C OUDIR = UDIR 111 CONTINUE IF ( MOD( NU - UTIME , TUDIM ) .EQ. 1 ) THEN UDIR = GRDU ELSE UDIR = CGU C C Compute CGu C TX = 0. TY = 0. DO 6 JT = 1, UDIM DO 6 NT = 1, TIME-1 TG = PMAT(JT)*(U(JT,NT,BESTU)-U(JT,NT,OLDU)) DO 7 I = 1, XYDIM 7 TG = TG - BMAT(I,JT)*(Y(I,NT+1,BESTV)-Y(I,NT+1,OLDV)) DO 9 I = 1, VDIM 9 TG = TG - DMAT(I,JT)*(V(I,NT,ALTV)-V(I,NT,OALTV)) TX = TX + TG*(U(JT,NT,BESTU)-U(JT,NT,GRDU)) TY = TY + TG*(U(JT,NT,OUDIR)-U(JT,NT,BESTU)) 6 CONTINUE DO 16 JT = 1, UDIM TG = PMAT(JT)*(U(JT,TIME,BESTU)-U(JT,TIME,OLDU)) DO 17 I = 1, XYDIM 17 TG = TG - BMAT(I,JT)*(YE(I,BESTV)-YE(I,OLDV)) DO 19 I = 1, VDIM 19 TG = TG - DMAT(I,JT)*(V(I,TIME,ALTV)-V(I,TIME,OALTV)) TX = TX + TG*(U(JT,TIME,BESTU)-U(JT,TIME,GRDU)) TY = TY + TG*(U(JT,TIME,OUDIR)-U(JT,TIME,BESTU)) 16 CONTINUE DO 1 JT = 1, UEDIM TG = PEMAT(JT)*(UE(JT,BESTU)-UE(JT,OLDU)) DO 8 I = 1, XYDIM 8 TG = TG - BEMAT(I,JT)*(Y(I,1,BESTV)-Y(I,1,OLDV)) TX = TX + TG*(UE(JT,BESTU)-UE(JT,GRDU)) TY = TY + TG*(UE(JT,OUDIR)-UE(JT,BESTU)) 1 CONTINUE IF ( TX .LT. 0. ) TX = 0. IF ( TY .LT. TEPS ) THEN T1 = 0.D0 ELSE T1 = TX/TY ENDIF DO 4 JT = 1, UDIM DO 4 NT = 1, TIME U(JT,NT,CGU)=(U(JT,NT,GRDU)+T1*U(JT,NT,OUDIR))/(1.D0+T1) 4 CONTINUE DO 3 JT = 1, UEDIM UE(JT,CGU) = (UE(JT,GRDU)+T1*UE(JT,OUDIR))/(1.D0+T1) 3 CONTINUE C ENDIF OLDU = BESTU OALTU = ALTU OGRDU = GRDU OCGU = CGU IF ( OLDU .LT. 4 ) THEN BESTU = 5 ALTU = 6 GRDU = 7 CGU = 8 ELSE BESTU = 1 ALTU = 2 GRDU = 3 CGU = 4 ENDIF C OVDIR = VDIR 222 CONTINUE IF ( MOD( NU - VTIME , TVDIM ) .EQ. 1 ) THEN VDIR = GRDV ELSE VDIR = CGV C C Compute CGv C TX = 0. TY = 0. DO 26 JT = 1, VDIM DO 26 NT = 2, TIME TG = QMAT(JT)*(V(JT,NT,BESTV)-V(JT,NT,OLDV)) DO 27 I = 1, XYDIM 27 TG = TG + CMAT(JT,I)*(X(I,NT-1,BESTU)-X(I,NT-1,OLDU)) DO 29 I = 1, UDIM 29 TG = TG + DMAT(JT,I)*(U(I,NT,ALTU)-U(I,NT,OALTU)) TX = TX + TG*(V(JT,NT,BESTV)-V(JT,NT,GRDV)) TY = TY + TG*(V(JT,NT,OUDIR)-V(JT,NT,BESTV)) 26 CONTINUE DO 36 JT = 1, VDIM TG = QMAT(JT)*(V(JT,1,BESTV)-V(JT,1,OLDV)) DO 37 I = 1, XYDIM 37 TG = TG + CMAT(JT,I)*(XE(I,BESTU)-XE(I,OLDU)) DO 39 I = 1, UDIM 39 TG = TG + DMAT(JT,I)*(U(I,1,ALTU)-U(I,1,OALTU)) TX = TX + TG*(V(JT,1,BESTV)-V(JT,1,GRDV)) TY = TY + TG*(V(JT,1,OUDIR)-V(JT,1,BESTV)) 36 CONTINUE DO 21 JT = 1, VEDIM TG = QEMAT(JT)*(VE(JT,BESTV)-VE(JT,OLDV)) DO 28 I = 1, XYDIM 28 TG = TG + CEMAT(JT,I)*(X(I,TIME,BESTV)-X(I,TIME,OLDV)) TX = TX + TG*(VE(JT,BESTV)-VE(JT,GRDV)) TY = TY + TG*(VE(JT,OUDIR)-VE(JT,BESTV)) 21 CONTINUE IF ( TX .LT. 0. ) TX = 0. IF ( TY .LT. TEPS ) THEN T1 = 0.D0 ELSE T1 = TX/TY ENDIF DO 24 JT = 1, VDIM DO 24 NT = 1, TIME V(JT,NT,CGV)=(V(JT,NT,GRDV)+T1*V(JT,NT,OVDIR))/(1.D0+T1) 24 CONTINUE DO 23 JT = 1, VEDIM VE(JT,CGV) = (VE(JT,GRDV)+T1*VE(JT,OVDIR))/(1.D0+T1) 23 CONTINUE C ENDIF OLDV = BESTV OALTV = ALTV OGRDV = GRDV OCGV = CGV IF ( OLDV .LT. 4 ) THEN BESTV = 5 ALTV = 6 GRDV = 7 CGV = 8 ELSE BESTV = 1 ALTV = 2 GRDV = 3 CGV = 4 ENDIF C - - - - - - - - - - - - - - - - - - - - - - - C C Compute the Trajectories OCGx , OCGy , for OCGu , OCGv C C IF (MOD(NU-UTIME,TUDIM).NE.1 .OR. MOD(NU-VTIME,TVDIM).NE.1) THEN CALL TRAJ(U(1,1,OCGU),UE(1,OCGU), : V(1,1,OCGV),VE(1,OCGV), : X(1,1,OCGU),XE(1,OCGU), : Y(1,1,OCGV),YE(1,OCGV), : AMAT,BMAT,CMAT,BVEC,CVEC, : BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) ENDIF C C C----------------------------------------------------------------------- C RETURN C End Of Subroutine 'UVSETS' END :::::::::::::: pxy.f :::::::::::::: c********************* c* subroutine TRAJ * c********************* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c The subroutine TRAJ transforms the control vectors 'u' and 'v' into c the state vectors 'x' and 'y'. The transformations are: c c x(t) = A x(t-1) + B u(t) + b c xe = Be ue + be c where t = 1, ntime ; xe: = x(0) c c y(t) = A' y(t+1) + C' v(t) + c c ye = Ce' ve + ce c where t = 1, ntime ; ye: = y(ntime+1) . c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE TRAJ(U,UE,V,VE,X,XE,Y,YE,AM,BM,CM,BV,CV, : BME,CME,BVE,CVE,NDU,NDUE,NXY,NDV,NDVE,NTIME) c c======================================================================= c IMPLICIT DOUBLE PRECISION (A-H,O-Z) c INTEGER NDU, NDV, NDUE, NDVE, NTIME, NXY c DOUBLE PRECISION U(NDU,NTIME), V(NDV,NTIME), UE(NDUE), VE(NDVE), : X(NXY,NTIME), Y(NXY,NTIME), XE(NXY), YE(NXY) c DOUBLE PRECISION AM(NXY,NXY), BM(NXY,NDU), CM(NDV,NXY), : BME(NXY,NDUE), CME(NDVE,NXY), : BV(NXY), CV(NXY), BVE(NXY), CVE(NXY) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c DOUBLE PRECISION SUM c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c c Compute the vectors X(t) from U(t): c DO 50 I=1,NXY SUM=0.0D0 DO 60 J=1,NDUE 60 SUM=SUM+BME(I,J)*UE(J) XE(I)=SUM+BVE(I) 50 CONTINUE c c DO 100 I=1,NXY SUM=0.0D0 DO 110 J=1,NXY 110 SUM=SUM+AM(I,J)*XE(J) DO 120 J=1,NDU 120 SUM=SUM+BM(I,J)*U(J,1) X(I,1)=SUM+BV(I) 100 CONTINUE c c DO 200 K=2,NTIME c DO 210 I=1,NXY SUM=0.0D0 DO 220 J=1,NXY 220 SUM=SUM+AM(I,J)*X(J,K-1) DO 230 J=1,NDU 230 SUM=SUM+BM(I,J)*U(J,K) X(I,K)=SUM+BV(I) 210 CONTINUE c 200 CONTINUE c c---------------------------------------------- c Compute the vectors Y(t) from V(t): c DO 250 I=1,NXY SUM=0.0D0 DO 260 J=1,NDVE 260 SUM=SUM+CME(J,I)*VE(J) YE(I)=SUM+CVE(I) 250 CONTINUE c c DO 300 I=1,NXY SUM=0.0D0 DO 310 J=1,NXY 310 SUM=SUM+AM(J,I)*YE(J) DO 320 J=1,NDV 320 SUM=SUM+CM(J,I)*V(J,NTIME) Y(I,NTIME)=SUM+CV(I) 300 CONTINUE c c DO 400 K=NTIME-1,1,-1 c DO 410 I=1,NXY SUM=0.0D0 DO 420 J=1,NXY 420 SUM=SUM+AM(J,I)*Y(J,K+1) DO 430 J=1,NDV 430 SUM=SUM+CM(J,I)*V(J,K) Y(I,K)=SUM+CV(I) 410 CONTINUE c 400 CONTINUE c c======================================================================= c c RETURN c END OF SUBROUTINE 'TRAJ' END :::::::::::::: search.f :::::::::::::: subroutine search(ct,atv,atm,bem,bm,btil,dtil,betil,detil, # zlft,zrt,zmid,brk1,brk2,slope, # zelft,zert,zemid,ebrk1,ebrk2,eslope, # lambda,z,ze,fval,zdim,zedim,ntime) c c INPUT: integer zdim,zedim,ntime double precision ct,atv,atm double precision bm(zdim),bem(zedim) double precision btil(zdim,ntime),dtil(zdim,ntime) double precision betil(zedim),detil(zedim) double precision zlft(zdim,ntime),zrt(zdim,ntime) double precision zmid(zdim,ntime) double precision brk1(zdim,ntime),brk2(zdim,ntime) double precision slope(zdim,ntime) double precision zelft(zedim),zert(zedim),zemid(zedim) double precision ebrk1(zedim),ebrk2(zedim),eslope(zedim) c c OUTPUT: double precision lambda,z(zdim,ntime),ze(zedim),fval c c dummy variables and counters: double precision lower,upper,a,b,fplus,fminus double precision fibo0,fibo1,fibo2,sum1,eps,ztemp integer i,n,k c c*********************************************************************** c This routine solves the problem c minimize f(lambda) = c c ct + atv*lambda + (1/2)lambda*lambda*atm c + sum / (btil(t)+lambda*dtil(t))*z(lambda;t) \ c t=1,N \ - (1/2)z(lambda;t)*bm*z(lambda;t) / c + { (betil+lambda*detil)*ze(lambda) - (1/2)ze(lambda)*bem*ze(lambda) } c c subject to 0.0<= lambda <=1.0 , where c / zlft(i,t) if lambda <= brk1(i,t) c z(lambda;i,t) = < zrt(i,t) if lambda > brk2(i,t) c \ zmid(i,t) + lambda*slope(i,t) otherwise. c and c / zelft(i) if lambda <= ebrk1(i) c ze(lambda;i) = < zert(i) if lambda > ebrk2(i) c \ zemid(i) + lambda*eslope(i) otherwise. c c It is assumed that the data were prepared by 'LNPREP' so that, for example, c z and ze are continuous. This also guarantees certain relationships among c the data that enable us to evaluate the derivatives with fewer arithmetical c operations. c c First the right-hand derivative at 0 and the left-hand derivative at 1 are c computed. If the former is positive we exit with lambda=0.0; if the latter c is negative then we take lambda=1.0. Then a "Bisection" search is used to c reduce the interval of uncertainty from [0,1] to [a,b] where b-a = 1.0d-12. c We check whether the right c hand derivative at the left test-point is positive, since the derivatives c require less computation. This also allows for greater accuracy since the c derivative tends toward zero (or else the test-points tend to a breakpoint). c c The routine also returns the function value 'fval' as well as the vectors c z(lambda), ze(lambda). c c======================================================================= c EPS is a distinguishability constant. eps = 1.0d-11 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c ================= c BISECTION SEARCH c ================= c a = 0.0d0 b = 1.0d0 c c COMPUTE the RH-derivative at a, LH-derivative at b fplus = atv fminus = atv + atm do 11 j=1,ntime do 10 i=1,zdim if ( 0.0d0 .lt. brk1(i,j) ) then fplus = fplus + dtil(i,j)*zlft(i,j) else if ( 0.0d0 .ge. brk2(i,j) ) then fplus = fplus + dtil(i,j)*zrt(i,j) else fplus = fplus + dtil(i,j)*zmid(i,j) endif if ( 1.0d0 .le. brk1(i,j) ) then fminus = fminus + dtil(i,j)*zlft(i,j) else if ( 1.0d0 .gt. brk2(i,j) ) then fminus = fminus + dtil(i,j)*zrt(i,j) else fminus = fminus + dtil(i,j)*(zmid(i,j) + slope(i,j)) endif 10 continue 11 continue do 12 i=1,zedim if ( 0.0d0 .lt. ebrk1(i) ) then fplus = fplus + detil(i)*zelft(i) else if ( 0.0d0 .ge. ebrk2(i) ) then fplus = fplus + detil(i)*zert(i) else fplus = fplus + detil(i)*zemid(i) endif if ( 1.0d0 .le. ebrk1(i) ) then fminus = fminus + detil(i)*zelft(i) else if ( 1.0d0 .gt. ebrk2(i) ) then fminus = fminus + detil(i)*zert(i) else fminus = fminus + detil(i)*(zemid(i) + eslope(i)) endif 12 continue c c if ( fplus .gt. -eps ) then c !! Stop with lambda=0.0d0; compute LH-derivative at lambda lambda = 0.0d0 fminus = atv do 150 j=1,ntime do 150 i=1,zdim if ( 0.0d0 .le. brk1(i,j)) then fminus = fminus + zlft(i,j)*dtil(i,j) else if ( 0.0d0 .gt. brk2(i,j) ) then fminus = fminus + zrt(i,j)*dtil(i,j) else fminus= fminus + zmid(i,j)*dtil(i,j) endif 150 continue do 151 i=1,zedim if ( 0.0d0 .le. ebrk1(i)) then fminus = fminus + zelft(i)*detil(i) else if ( 0.0d0 .gt. ebrk2(i) ) then fminus = fminus + zert(i)*detil(i) else fminus= fminus + zemid(i)*detil(i) endif 151 continue c else if ( fminus .lt. eps ) then c !! Stop with lambda=1.0d0; compute RH-derivative at lambda lambda = 1.0d0 fplus = atv + atm do 110 j=1,ntime do 110 i=1,zdim if ( 1.0d0 .lt. brk1(i,j) ) then fplus = fplus + zlft(i,j)*dtil(i,j) else if ( 1.0d0 .ge. brk2(i,j) ) then fplus = fplus + zrt(i,j)*dtil(i,j) else fplus = fplus + (zmid(i,j)+slope(i,j))*dtil(i,j) endif 110 continue do 111 i=1,zedim if ( 1.0d0 .lt. ebrk1(i) ) then fplus = fplus + zelft(i)*detil(i) else if ( 1.0d0 .ge. ebrk2(i) ) then fplus = fplus + zert(i)*detil(i) else fplus = fplus + (zemid(i)+eslope(i))*detil(i) endif 111 continue else c !! Search inside the interval [0,1] for fplus=0.0d0 c n =40 c lower = (a + b)/2.0d0 c COMPUTE the function values at lower c k = 1 60 fplus = atv + lower*atm do 51 j=1,ntime do 50 i=1,zdim if ( lower .lt. brk1(i,j) ) then fplus = fplus + dtil(i,j)*zlft(i,j) else if ( lower .ge. brk2(i,j) ) then fplus = fplus + dtil(i,j)*zrt(i,j) else fplus = fplus + dtil(i,j)*(zmid(i,j)+slope(i,j)*lower) endif 50 continue 51 continue do 52 i=1,zedim if ( lower .lt. ebrk1(i) ) then fplus = fplus + detil(i)*zelft(i) else if ( lower .ge. ebrk2(i) ) then fplus = fplus + detil(i)*zert(i) else fplus = fplus + detil(i)*(zemid(i)+eslope(i)*lower) endif 52 continue if ( dabs(fplus) .le. eps ) then c !! LOWER is optimal lambda = lower else if ( fplus .lt. -eps ) then c !! So we move A up, and find a new LOWER. a = lower lower = (a + b)/2.0d0 else if ( fplus .gt. eps ) then c !! So we move B down, and find a new LOWER. b = lower lower = (a + b)/2.0d0 endif k = k + 1 if ( k .lt. n) then go to 60 else lambda = (a+b)/2.0d0 endif endif c COMPUTE the right-hand derivative at lambda fplus = atv + lambda*atm fminus = atv + lambda*atm do 80 j=1,ntime do 80 i=1,zdim if ( lambda .lt. brk1(i,j) ) then fplus = fplus + zlft(i,j)*dtil(i,j) else if ( lambda .ge. brk2(i,j) ) then fplus = fplus + zrt(i,j)*dtil(i,j) else fplus = fplus + (zmid(i,j)+slope(i,j)*lambda)*dtil(i,j) endif if ( lambda .le. brk1(i,j) ) then fminus = fminus + zlft(i,j)*dtil(i,j) else if ( lambda .gt. brk2(i,j) ) then fminus = fminus + zrt(i,j)*dtil(i,j) else fminus = fminus+(zmid(i,j)+slope(i,j)*lambda)*dtil(i,j) endif 80 continue do 81 i=1,zedim if ( lambda .lt. ebrk1(i) ) then fplus = fplus + zelft(i)*detil(i) else if ( lambda .ge. ebrk2(i) ) then fplus = fplus + zert(i)*detil(i) else fplus = fplus + (zemid(i)+eslope(i)*lambda)*detil(i) endif if ( lambda .le. ebrk1(i) ) then fminus = fminus + zelft(i)*detil(i) else if ( lambda .gt. ebrk2(i) ) then fminus = fminus + zert(i)*detil(i) else fminus = fminus + (zemid(i)+eslope(i)*lambda)*detil(i) endif 81 continue endif c c End of Bisection Search c ===-==-=========-====== c c c Calculate z(lambda) and the optimal value fval. fval = ct + atv*lambda + lambda*lambda*atm/2.0d0 do 300 j=1,ntime do 300 i=1,zdim if ( lambda .lt. brk1(i,j) ) then ztemp = zlft(i,j) else if ( lambda .gt. brk2(i,j) ) then ztemp = zrt(i,j) else if ( (fminus .lt. -eps) .and. (lambda .eq. brk1(i,j)) # .and. ( brk1(i,j) .eq. brk2(i,j)) ) then tt = fminus/(fminus-fplus) ztemp = (1.0d0-tt)*zlft(i,j) + tt*zrt(i,j) else ztemp = zmid(i,j) + slope(i,j)*lambda endif z(i,j) = ztemp sum1 = btil(i,j) + dtil(i,j)*lambda - ztemp*bm(i)/2.0d0 fval = fval + sum1*ztemp 300 continue do 301 i=1,zedim if ( lambda .lt. ebrk1(i) ) then ztemp = zelft(i) else if ( lambda .gt. ebrk2(i) ) then ztemp = zert(i) else if ( (fminus .lt. -eps) .and. (lambda .eq. ebrk1(i)) # .and. ( ebrk1(i) .eq. ebrk2(i)) ) then tt = fminus/(fminus-fplus) ztemp = (1.0d0-tt)*zelft(i) + tt*zert(i) else ztemp = zemid(i) + eslope(i)*lambda endif ze(i) = ztemp sum1 = betil(i) + detil(i)*lambda - ztemp*bem(i)/2.0d0 fval = fval + sum1*ztemp 301 continue c c c c END OF SUBROUTINE search return end :::::::::::::: tempx :::::::::::::: rm data.* nice -19 a.out7pz0.o1 nice -19 pbdynpz0.b1 nice -19 podynpz1.o1 nice -19 pbdynpz1.b1 nice -19 podynpz2.o1 nice -19 pbdynpz2.b1 nice -19 podynpz3.o1 nice -19 pbdynpz3.b1 nice -19 podynpz4.o1 nice -19 pbdynpz4.b1 nice -19 podynpz5.o1 nice -19 pbdynpz5.b1 nice -19 podynpz6.o1 nice -19 pbdynpz6.b1 nice -19 podynpz7.o1 nice -19 pbdynpz7.b1 nice -19 podynpz8.o1 nice -19 pbdynpz8.b1 nice -19 podynpz9.o1 nice -19 pbdynpz9.b1 nice -19 pbdz3pz0.e1 nice -19 pbdz3pz1.e1 nice -19 pbdz3pz2.e1 nice -19 pbdz3pz3.e1 nice -19 pbdz3pz4.e1 nice -19 pbdz3pz5.e1 nice -19 pbdz3pz6.e1 nice -19 pbdz3pz7.e1 nice -19 pbdz3pz8.e1 nice -19 pbdz3pz9.e1 nice -19 grep Outer pz*.*1 > piter1.7.1 nice -19 grep time: pz*.*1 > ptime1.7.1 rm pz*.*1 nice -19 podynpz0.o2 nice -19 pbdynpz0.b2 nice -19 podynpz1.o2 nice -19 pbdynpz1.b2 nice -19 podynpz2.o2 nice -19 pbdynpz2.b2 nice -19 podynpz3.o2 nice -19 pbdynpz3.b2 nice -19 podynpz4.o2 nice -19 pbdynpz4.b2 nice -19 podynpz5.o2 nice -19 pbdynpz5.b2 nice -19 podynpz6.o2 nice -19 pbdynpz6.b2 nice -19 podynpz7.o2 nice -19 pbdynpz7.b2 nice -19 podynpz8.o2 nice -19 pbdynpz8.b2 nice -19 podynpz9.o2 nice -19 pbdynpz9.b2 nice -19 pbdz3pz0.e2 nice -19 pbdz3pz1.e2 nice -19 pbdz3pz2.e2 nice -19 pbdz3pz3.e2 nice -19 pbdz3pz4.e2 nice -19 pbdz3pz5.e2 nice -19 pbdz3pz6.e2 nice -19 pbdz3pz7.e2 nice -19 pbdz3pz8.e2 nice -19 pbdz3pz9.e2 nice -19 grep Outer pz*.*2 > piter1.7.2 nice -19 grep time: pz*.*2 > ptime1.7.2 rm pz*.*2 rm data.* nice -19 a.out8pz0.o1 nice -19 pbdynpz0.b1 nice -19 podynpz1.o1 nice -19 pbdynpz1.b1 nice -19 podynpz2.o1 nice -19 pbdynpz2.b1 nice -19 podynpz3.o1 nice -19 pbdynpz3.b1 nice -19 podynpz4.o1 nice -19 pbdynpz4.b1 nice -19 podynpz5.o1 nice -19 pbdynpz5.b1 nice -19 podynpz6.o1 nice -19 pbdynpz6.b1 nice -19 podynpz7.o1 nice -19 pbdynpz7.b1 nice -19 podynpz8.o1 nice -19 pbdynpz8.b1 nice -19 podynpz9.o1 nice -19 pbdynpz9.b1 nice -19 pbdz3pz0.e1 nice -19 pbdz3pz1.e1 nice -19 pbdz3pz2.e1 nice -19 pbdz3pz3.e1 nice -19 pbdz3pz4.e1 nice -19 pbdz3pz5.e1 nice -19 pbdz3pz6.e1 nice -19 pbdz3pz7.e1 nice -19 pbdz3pz8.e1 nice -19 pbdz3pz9.e1 nice -19 grep Outer pz*.*1 > piter1.8.1 nice -19 grep time: pz*.*1 > ptime1.8.1 rm pz*.*1 nice -19 podynpz0.o2 nice -19 pbdynpz0.b2 nice -19 podynpz1.o2 nice -19 pbdynpz1.b2 nice -19 podynpz2.o2 nice -19 pbdynpz2.b2 nice -19 podynpz3.o2 nice -19 pbdynpz3.b2 nice -19 podynpz4.o2 nice -19 pbdynpz4.b2 nice -19 podynpz5.o2 nice -19 pbdynpz5.b2 nice -19 podynpz6.o2 nice -19 pbdynpz6.b2 nice -19 podynpz7.o2 nice -19 pbdynpz7.b2 nice -19 podynpz8.o2 nice -19 pbdynpz8.b2 nice -19 podynpz9.o2 nice -19 pbdynpz9.b2 nice -19 pbdz3pz0.e2 nice -19 pbdz3pz1.e2 nice -19 pbdz3pz2.e2 nice -19 pbdz3pz3.e2 nice -19 pbdz3pz4.e2 nice -19 pbdz3pz5.e2 nice -19 pbdz3pz6.e2 nice -19 pbdz3pz7.e2 nice -19 pbdz3pz8.e2 nice -19 pbdz3pz9.e2 nice -19 grep Outer pz*.*2 > piter1.8.2 nice -19 grep time: pz*.*2 > ptime1.8.2 rm pz*.*2 rm data.* nice -19 a.out7pz0.o1 nice -19 pbdynpz0.b1 nice -19 podynpz1.o1 nice -19 pbdynpz1.b1 nice -19 podynpz2.o1 nice -19 pbdynpz2.b1 nice -19 podynpz3.o1 nice -19 pbdynpz3.b1 nice -19 podynpz4.o1 nice -19 pbdynpz4.b1 nice -19 podynpz5.o1 nice -19 pbdynpz5.b1 nice -19 podynpz6.o1 nice -19 pbdynpz6.b1 nice -19 podynpz7.o1 nice -19 pbdynpz7.b1 nice -19 podynpz8.o1 nice -19 pbdynpz8.b1 nice -19 podynpz9.o1 nice -19 pbdynpz9.b1 nice -19 pbdz3pz0.e1 nice -19 pbdz3pz1.e1 nice -19 pbdz3pz2.e1 nice -19 pbdz3pz3.e1 nice -19 pbdz3pz4.e1 nice -19 pbdz3pz5.e1 nice -19 pbdz3pz6.e1 nice -19 pbdz3pz7.e1 nice -19 pbdz3pz8.e1 nice -19 pbdz3pz9.e1 nice -19 grep Outer pz*.*1 > piter2.7.1 nice -19 grep time: pz*.*1 > ptime2.7.1 rm pz*.*1 nice -19 podynpz0.o2 nice -19 pbdynpz0.b2 nice -19 podynpz1.o2 nice -19 pbdynpz1.b2 nice -19 podynpz2.o2 nice -19 pbdynpz2.b2 nice -19 podynpz3.o2 nice -19 pbdynpz3.b2 nice -19 podynpz4.o2 nice -19 pbdynpz4.b2 nice -19 podynpz5.o2 nice -19 pbdynpz5.b2 nice -19 podynpz6.o2 nice -19 pbdynpz6.b2 nice -19 podynpz7.o2 nice -19 pbdynpz7.b2 nice -19 podynpz8.o2 nice -19 pbdynpz8.b2 nice -19 podynpz9.o2 nice -19 pbdynpz9.b2 nice -19 pbdz3pz0.e2 nice -19 pbdz3pz1.e2 nice -19 pbdz3pz2.e2 nice -19 pbdz3pz3.e2 nice -19 pbdz3pz4.e2 nice -19 pbdz3pz5.e2 nice -19 pbdz3pz6.e2 nice -19 pbdz3pz7.e2 nice -19 pbdz3pz8.e2 nice -19 pbdz3pz9.e2 nice -19 grep Outer pz*.*2 > piter2.7.2 nice -19 grep time: pz*.*2 > ptime2.7.2 rm pz*.*2 rm data.* nice -19 a.out8pz0.o1 nice -19 pbdynpz0.b1 nice -19 podynpz1.o1 nice -19 pbdynpz1.b1 nice -19 podynpz2.o1 nice -19 pbdynpz2.b1 nice -19 podynpz3.o1 nice -19 pbdynpz3.b1 nice -19 podynpz4.o1 nice -19 pbdynpz4.b1 nice -19 podynpz5.o1 nice -19 pbdynpz5.b1 nice -19 podynpz6.o1 nice -19 pbdynpz6.b1 nice -19 podynpz7.o1 nice -19 pbdynpz7.b1 nice -19 podynpz8.o1 nice -19 pbdynpz8.b1 nice -19 podynpz9.o1 nice -19 pbdynpz9.b1 nice -19 pbdz3pz0.e1 nice -19 pbdz3pz1.e1 nice -19 pbdz3pz2.e1 nice -19 pbdz3pz3.e1 nice -19 pbdz3pz4.e1 nice -19 pbdz3pz5.e1 nice -19 pbdz3pz6.e1 nice -19 pbdz3pz7.e1 nice -19 pbdz3pz8.e1 nice -19 pbdz3pz9.e1 nice -19 grep Outer pz*.*1 > piter2.8.1 nice -19 grep time: pz*.*1 > ptime2.8.1 rm pz*.*1 nice -19 podynpz0.o2 nice -19 pbdynpz0.b2 nice -19 podynpz1.o2 nice -19 pbdynpz1.b2 nice -19 podynpz2.o2 nice -19 pbdynpz2.b2 nice -19 podynpz3.o2 nice -19 pbdynpz3.b2 nice -19 podynpz4.o2 nice -19 pbdynpz4.b2 nice -19 podynpz5.o2 nice -19 pbdynpz5.b2 nice -19 podynpz6.o2 nice -19 pbdynpz6.b2 nice -19 podynpz7.o2 nice -19 pbdynpz7.b2 nice -19 podynpz8.o2 nice -19 pbdynpz8.b2 nice -19 podynpz9.o2 nice -19 pbdynpz9.b2 nice -19 pbdz3pz0.e2 nice -19 pbdz3pz1.e2 nice -19 pbdz3pz2.e2 nice -19 pbdz3pz3.e2 nice -19 pbdz3pz4.e2 nice -19 pbdz3pz5.e2 nice -19 pbdz3pz6.e2 nice -19 pbdz3pz7.e2 nice -19 pbdz3pz8.e2 nice -19 pbdz3pz9.e2 nice -19 grep Outer pz*.*2 > piter2.8.2 nice -19 grep time: pz*.*2 > ptime2.8.2 rm pz*.*2 rm data.* #################################################### End 92007-2.for