:::::::::::::: readme.p :::::::::::::: Disclaimer: The program is for EXPERIMENTATION PURPOSE ONLY. It has NOT been refined. Actually it was adapted from Stephen E. Wright's code of the finite-generation method. (See Reference 22 in the paper.) I changed the subroutines, but didn't find time to change all the comments correspondingly. So some of THE COMMENTS BECOME QUITE INACCURATE." --- Ciyou Zhu Following is a clue on how to generate the test results in the paper: 0) All the files of the program are packed in one by commands like more pc*.* >pc.sendout You should first restore each individual file with the CORRECT name in a working directory. The control files with the name pc*.* should contain exactly 18 lines. 1) Issue the following commands to compile: f77 pgdt7.f -o a.out7 f77 pgdt8.f -o a.out8 make -f mpbznw (This will generate MPPA executable file pbdznew.) make -f mponw (This will generate PPA executable file podynew.) make -f mpnw (This will generate RPPA executable file pdynew.) mv pbdznew pbdyn mv podynew podyn mv pdynew pbdz3 2) Make the file with the name tempx executable. MAKE SURE THAT IN THE WORKING DIRECTORY THERE ARE NO OTHER FILES WITH THE NAME data.* OR pz*.*. The next command will ERASE all such files. Then issue the command: tempx After tempx finished its job, the numbers of iteration in the paper will be sitting in the files with the name piter*.*.*. Note that the files with the name pz*.b* contain the results of MPPA; the files with the name pz*.o* contain the results of PPA; the files with the name pz*.e* contain the results of RPPA. :::::::::::::: pc0.1 :::::::::::::: name of problem Program #0 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.0 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc0.2 :::::::::::::: name of problem Program #10 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.0 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc1.1 :::::::::::::: name of problem Program #1 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.1 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc1.2 :::::::::::::: name of problem Program #11 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.1 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc2.1 :::::::::::::: name of problem Program #2 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.2 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc2.2 :::::::::::::: name of problem Program #12 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.2 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc3.1 :::::::::::::: name of problem Program #3 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of inner iterations 60 max # of outer iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.3 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc3.2 :::::::::::::: name of problem Program #13 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of inner iterations 60 max # of outer iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.3 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc4.1 :::::::::::::: name of problem Program #4 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.4 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc4.2 :::::::::::::: name of problem Program #14 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.4 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc5.1 :::::::::::::: name of problem Program #5 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.5 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc5.2 :::::::::::::: name of problem Program #15 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.5 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc6.1 :::::::::::::: name of problem Program #6 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.6 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc6.2 :::::::::::::: name of problem Program #16 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.6 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc7.1 :::::::::::::: name of problem Program #7 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.7 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc7.2 :::::::::::::: name of problem Program #17 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.7 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc8.1 :::::::::::::: name of problem Program #8 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.8 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc8.2 :::::::::::::: name of problem Program #18 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.8 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc9.1 :::::::::::::: name of problem Program #9 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 4 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.9 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: pc9.2 :::::::::::::: name of problem Program #19 continuous time? true primal control dimension 20 primal endpoint dimension 20 trajectory dimension 20 dual control dimension 20 dual endpoint dimension 20 n-# of time periods 16 max # of outer iterations 60 max # of inner iterations 30 duality gap CUTOFF 1.d-8 mode-# of elements gen'd 1 initial u given in input file false initial v given in input file false name of Input Data file data.9 name of Auxiliary Output file aux.out.t print options (2,2,3,0)(0,3,1,7)(0,3,2,0) ****************************** :::::::::::::: mpbznw :::::::::::::: FC = f77 SOR1 = prxbz.f pdfgm.f pfgval.f pmain.f pdynfgm.f pxy.f pinput.f poutput.f OBJ1 = prxbz.o pdfgm.o pfgval.o pmain.o pdynfgm.o pxy.o pinput.o poutput.o SOR2 = pinit.f puvsets.f pline.f OBJ2 = pinit.o puvsets.o pline.o SOR3 = pcoeff.f pobjvl.f plnprep.f plnsrch.f search.f OBJ3 = pcoeff.o pobjvl.o plnprep.o plnsrch.o search.o pbdznew : $(OBJ1) $(OBJ2) $(OBJ3) $(FC) $(OBJ1) $(OBJ2) $(OBJ3) -o pbdznew :::::::::::::: mpnw :::::::::::::: FC = f77 SOR1 = pprox.f pdfgm.f pfgval.f pmain.f pdynfgm.f pxy.f pinput.f poutput.f OBJ1 = pprox.o pdfgm.o pfgval.o pmain.o pdynfgm.o pxy.o pinput.o poutput.o SOR2 = pinit.f puvsets.f pline.f OBJ2 = pinit.o puvsets.o pline.o SOR3 = pcoeff.f pobjvl.f plnprep.f plnsrch.f search.f OBJ3 = pcoeff.o pobjvl.o plnprep.o plnsrch.o search.o pdynew : $(OBJ1) $(OBJ2) $(OBJ3) $(FC) $(OBJ1) $(OBJ2) $(OBJ3) -o pdynew :::::::::::::: mponw :::::::::::::: FC = f77 SOR1 = proxo.f pdfgm.f pfgval.f pmain.f pdynfgm.f pxy.f pinput.f poutput.f OBJ1 = proxo.o pdfgm.o pfgval.o pmain.o pdynfgm.o pxy.o pinput.o poutput.o SOR2 = pinit.f puvsets.f pline.f OBJ2 = pinit.o puvsets.o pline.o SOR3 = pcoeff.f pobjvl.f plnprep.f plnsrch.f search.f OBJ3 = pcoeff.o pobjvl.o plnprep.o plnsrch.o search.o podynew : $(OBJ1) $(OBJ2) $(OBJ3) $(FC) $(OBJ1) $(OBJ2) $(OBJ3) -o podynew :::::::::::::: pcoeff.f :::::::::::::: subroutine coeff(pbvec,pbmat,qbvec,qbmat,dbmat, $ u,x,v,y,ue,xe,ve,ye, $ pevec,pvec,pemat,pmat, $ qvec,qevec,qmat,qemat,dmat, $ cmat,cvec,cemat,cevec, $ n,nufe,nvfe, $ udim,xydim,vdim,uedim,vedim) c c/////////////////////////////////////////////////////////////////////////// c// NOTE : THIS ROUTINE ASSUMES THE MATRICES " P, P , Q, and Q " ARE // c// e e // c// ALL DIAGONAL AND ARE GIVEN HERE BY A 1-DIM'L ARRAY GIVING // c// THE DIAGONAL ELEMENTS. // c/////////////////////////////////////////////////////////////////////////// c c*********************************************************************** c c This subroutine takes as input the data for the following problem: c c minimax J(u,u ;v,v ) over (co U x U ) x (co V x V ). c e e e e c c Here U is a finite set of points in Euclidean space of the form c (u,u ), where u is of the form u(j,t) with j=1,...,k and t=1,...,N. c e c Similarly, V is a finite set of points of the form v(i,t) with i=1,...,l c and t=1,...,N. 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 c e e e e e c c c 1 1 c where J (u(t),v(t)) = pu + qv + -u Pu - -v Qv - v Du c c 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 c The output consists of the coefficients of the equivalent problem: c c minimax G(g,h) over g(1)+...+g(m) = 1 , g(j)>0 for j=1,...,m c h(1)+...+h(n) = 1 , h(i)>0 for i=1,...,n c c where _ 1 _ _ 1 _ _ c G(g,h) = pg + -gPg + hq - -hQh - gDh . c 2 2 c c The transformation is given by: c _ N c p(j) = p *u(j) + sum p*u(t;j) c e e t=1 c _ N c P(j,jj) = u(j)*P u(jj) + sum u(t;j)*Pu(t;jj) c e e e t=1 c _ N c q(i) = sum q*v(t;i) + q *v(i) c t=1 e e c _ N c Q(i,ii) = sum v(t;i)*Qv(t;ii) + v(i)*Q v(ii) c t=1 e e e c c _ N N T c D(i,j) = sum v(t;i)*Du(t;j) + sum x(t-1;j)*[C v(t;i) + c] c t=1 t=1 c T c + x(N;j)*[C v(N+1;i) + c ] c e e c c Here u(t;j) is the k-dimensional vector corresponding to the t-th c component of the i-th element of the finite set U, whereas c j c u is the j-th k-by-N element of the finite set U. Similarly, c for v. c*********************************************************************** c c The variables declared below correspond to the above variables as c follows: c _ c pbvec = p c _ c pbmat = diagonal elements of P c _ c qbvec = q c _ c qbmat = diagonal elements of Q c _ c dbmat = D c u(j,t,nu) = u(j,t) for the nu-th element in U c x(j,t,nu) = x(j,t) " " " " " U c v(i,t,nu) = v(i,t) for the nu-th element in V c y(i,t,nu) = y(i,t) " " " " " V c ue(j,nu) = endpoint vector u for nu-th element c xe(j,nu) = endpoint vector u for nu-th element c ve(i,nu) = endpoint vector v for nu-th element c ye(i,nu) = endpoint vector v for nu-th element c pevec = p c e c pvec = p c pemat = diagonal elements of P c e c pmat = diagonal elements of P c qvec = q c qevec = q c e c qmat = diagonal elements of Q c qemat = diagonal elements of Q c e c dmat = D c cmat = C c cvec = c c cemat = C c e c cevec = c c e c n = N c nvfe = n c nufe = m c udim = k c vdim = l c xydim = dimension of state space for x and y c uedim = k c e c vedim = l c e c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c VARIABLES PASSED THRU SUBROUTINE CALL c integer n,nvfe,nufe integer udim,xydim,vdim integer uedim,vedim c real*8 pbvec(nufe),pbmat(nufe,nufe),qbvec(nvfe),qbmat(nvfe,nvfe) real*8 dbmat(nvfe,nufe) c real*8 u(udim,n,nufe),x(xydim,n,nufe) real*8 ue(uedim,nufe),xe(xydim,nufe) real*8 v(vdim,n,nvfe),y(xydim,n,nvfe) real*8 ve(vedim,nvfe),ye(xydim,nvfe) c real*8 pevec(uedim),pvec(udim),pemat(uedim) real*8 pmat(udim),qvec(vdim),qevec(vedim) real*8 qmat(vdim),qemat(vedim),dmat(vdim,udim) real*8 cmat(vdim,xydim),cvec(xydim),cemat(vedim,xydim) real*8 cevec(xydim) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLES c real*8 sum1,sum2 integer i,ii,j,jj,k,kk,t c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c do 170 j=1,nufe c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c _ c COMPUTE p (pbvec) c sum1 = 0.0d0 do 5 k=1,uedim 5 sum1 = sum1 + pevec(k)*ue(k,j) do 20 k=1,udim sum2 = 0.0d0 do 10 t=1,n 10 sum2 = sum2 + u(k,t,j) sum1 = sum1 + pvec(k)*sum2 20 continue pbvec(j) = sum1 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c _ c COMPUTE P (pbmat) c do 80 jj=1,nufe sum1 = 0.0d0 do 40 k=1,uedim sum1 = sum1 + ue(k,j)*pemat(k)*ue(k,jj) 40 continue do 70 t=1,n do 60 k=1,udim sum1 = sum1 + u(k,t,j)*pmat(k)*u(k,t,jj) 60 continue 70 continue pbmat(j,jj) = sum1 80 continue c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c _ c COMPUTE D (dbmat) c do 160 i=1,nvfe sum1 = 0.0d0 do 82 k=1,vdim sum2 = 0.0d0 do 81 kk=1,udim 81 sum2 = sum2 + dmat(k,kk)*u(kk,1,j) sum1 = sum1 + v(k,1,i)*sum2 82 continue do 84 kk=1,xydim sum2 = 0.0d0 do 83 k=1,vdim 83 sum2 = sum2 + cmat(k,kk)*v(k,1,i) sum1 = sum1 + xe(kk,j)*(sum2+cvec(kk)) 84 continue do 130 t=2,n do 100 k=1,vdim sum2 = 0.0d0 do 90 kk=1,udim 90 sum2 = sum2 + dmat(k,kk)*u(kk,t,j) sum1 = sum1 + v(k,t,i)*sum2 100 continue do 120 kk=1,xydim sum2 = 0.0d0 do 110 k=1,vdim 110 sum2 = sum2 + cmat(k,kk)*v(k,t,i) sum1 = sum1 + x(kk,t-1,j)*(sum2+cvec(kk)) 120 continue 130 continue do 150 kk=1,xydim sum2 = 0.0d0 do 140 k=1,vedim 140 sum2 = sum2 + cemat(k,kk)*ve(k,i) sum1 = sum1 + x(kk,n,j)*(sum2+cevec(kk)) 150 continue dbmat(i,j) = sum1 160 continue c - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 continue c c c do 260 i=1,nvfe c - - - - - - - - - - - - - - - - - - - - - - - - - - - c _ c COMPUTE q (qbvec) c sum1 = 0.0d0 do 175 k=1,vedim 175 sum1 = sum1 + qevec(k)*ve(k,i) do 190 k=1,vdim sum2 = 0.0d0 do 180 t=1,n 180 sum2 = sum2 + v(k,t,i) sum1 = sum1 + qvec(k)*sum2 190 continue qbvec(i) = sum1 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c _ c COMPUTE Q (qbmat) c do 250 ii=1,nvfe sum1 = 0.0d0 do 210 k=1,vedim sum1 = sum1 + ve(k,i)*qemat(k)*ve(k,ii) 210 continue do 240 t=1,n do 230 k=1,vdim sum1 = sum1 + v(k,t,i)*qmat(k)*v(k,t,ii) 230 continue 240 continue qbmat(i,ii) = sum1 250 continue c c - - - - - - - - - - - - - - - - - - - - - - - - - - - 260 continue c c c======================================================================= c c return c END OF subroutine 'COEFF' end :::::::::::::: pdfgm.f :::::::::::::: c The following code is a slightly modified version of Stephen E. Wright's c original code. SUBROUTINE DFGM(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, $ ALPHA,BETA,BESTU,BESTV,EPS,EPSNRM,INFORM,NU, $ WORK,LENWRK,UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME, $ MAXUFE,MAXVFE,NUMAX,DEPS) C C C----------------------------------------------------------------------- C/////////////////////////////////////////////////////////////////////// C/////////////////////////////////////////////////////////////////////// C----------------------------------------------------------------------- C Variables passed thru subroutine call C INTEGER LENWRK INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME INTEGER MAXUFE,MAXVFE,NUMAX INTEGER BESTU,BESTV INTEGER INFORM,NU 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) DOUBLE PRECISION ALPHA(MAXUFE),BETA(MAXVFE) 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 positive c e e c entries on the diagonal. (Certain aspects of the routines written by S.E. c Wright allow for some diagonal entries to be zero.) Also associated with c the problem are the maps 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 The method employed is user-dependent. This routine simply acts as the c primary driver of the actual algorithm. It is expected that the method c used will be similar in structure to the Finite-Generation Method. For c this reason, the data structures have been defined so as to keep track c of "finite sets" of primal and dual vectors. The structure of the algorithm c 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 "numax". c Step 0. (a) Initialization routine (INIT): compute the function c values for the initial guesses and the corresponding c alternative elements giving the max and min defining c these values; compute the trajectories x and y of the c initial guesses and the alternative elements. c (b) Set nu=1. c Step 1. Creation of finite sets (UVSETS): set up the finite sets c of primal and dual vectors which will be used in the direction- c finding procedure, as well as their trajectories; designate c the positions in the U and V arrays in which the direction c vectors, new "best" vectors, and new "alternative" will be c placed; initialize the vectors of barycentric coordinates c which give a convex combination of the elements in the c finite sets. c Step 2. Direction-finding routine (DIRECT): use the finite U and V c sets to generate a "direction" point toward which to linesearch c Step 3. (a) Linesearch routine (LINE): linesearch between the old c "best" elements and the new "direction" elements to obtain c the new "best" elements (and their trajectories), as well c as the corresponding function values and "alternative" c elements (along with their trajectories). c (b) Compute the duality gap EPSILON and normalized duality gap. c Decision Step: If the normalized duality gap is less than 'DEPS', then c stop: the pair (u,ue) and the pair (v,ve) are both EPSILON- c optimal. Otherwise, if nu>>>>>> INITIALIZATION <<<<<<<< C QUIT = .FALSE. C C Compute the initial Finite-Element Sets, including the C Alternative Elements U(ALTU),V(ALTV) which give the C respective maximum and minimum defining the function C values for the Initial Guesses V(BESTV),U(BESTU), and C the Function Values for the Initial Guesses. C CALL INIT(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ 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,NUFE,NVFE, $ GRDU,GRDV,CGU,CGV, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME,MAXUFE,MAXVFE, $ WORK,LENWRK,QUIT) C 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 C C Number Of Iterations: NU = 0 C C C Count the Number of Iterations in each C.G. Cycle: UTIME = 0 VTIME = 0 C C C Pointers indicating positions of Direction Vectors: UDIR = 0 VDIR = 0 C C//////////////////////////// C>>>>>>> Begin Loop <<<<<<<< C>>>>>>> ========== <<<<<<<< C//////////////////////////// C 1000 CONTINUE C IF ((EPSNRM.GT.DEPS).AND.(NU.LT.NUMAX).AND.(.NOT. QUIT)) THEN C C i.e. While ( EPS > DEPS and NU < NUMAX and NOT QUIT ) C ===== --- ---- === -- ----- === --- ---- C Do the following: C == === ========= C NU = NU + 1 C C - - - - - - - - - - - - - - - - - - - - - - - - - C Set up the Finite Sets; Designate the Positions in C which the Direction Vectors should be put, and the C positions for the New 'Best' Guesses and their C "Alternative" Elements. C CALL 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 Linesearch between the "Old Best" elements and the new C "Direction" elements to get the new "Best" elements, their C trajectories, and the corresponding "Alternative" elements. C CALL LINE(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ OLDU,OLDV,UDIR,VDIR,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 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 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 IF ( NU .GE. NUMAX ) THEN C Maximum number of iterations exceeded ! INFORM = 1 ELSE C A User-supplied Exit Criterion has been satisfied ! INFORM = 2 ENDIF C C======================================================================= C RETURN C End of subroutine 'DFGM' END :::::::::::::: pdynfgm.f :::::::::::::: c The following code is a slightly modified version of Stephen E. Wright's c original code. subroutine dynfgm(zzd,nzdlen) c c This routine calls the routine INPUT to get the data for the DYNFGM c program, sets up the working array ZZD, (possibly) obtains the initial c guesses for the vectors U and V , and then calls the routine PROX. c======================================================================= c Variables passed into DYNFGM. c integer nzdlen double precision zzd(nzdlen) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c Local variables c integer n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13, $ n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n24, $ n25,n26,n27,n28,n29,n30,n31,n32,n33,n34,n35, $ n36,n37,n38,n39,n40,n41,n42,n43,n44,n45,n46, $ n47,ncheck integer udim,uedim,xydim,vdim,vedim,mumax,numax,n integer maxufe,maxvfe,mode integer propt(3,5),out c double precision epsmin,dval c logical initu,initv,keepu,keepv c character*60 name,ofile character*64 rsfile c integer lnsize,sadsiz double precision dummy 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. common /prdata/ propt 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 C 3. Locally needed variable of LINE which must be saved double precision gap common /epslin/ gap C c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c call input(zzd,dval,udim,uedim,xydim,vdim,vedim,mumax,numax, $ epsmin,n,name,ofile,propt,initu,initv,mode,nzdlen) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c The following two parameters have been assigned values here so as to c avoid having to change the input routine for the time being. In theory, c they should, in fact, be read in with the rest of the control parameters. keepu = .false. keepv = .false. c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c if ( keepu ) then maxufe = (numax+1)*(mode+2) + 1 else maxufe = 16 endif if ( keepv ) then maxvfe = (numax+1)*(mode+2) + 1 else maxvfe = 16 endif c c zzd( n1 : n1+n*udim*maxufe-1 ) = U (i.e. the finite element set) n1 = 1 c c zzd( n2 : n2+n*vdim*maxvfe-1 ) = V (i.e. the finite element set) n2 = n1 + n*udim*maxufe c c zzd( n3 : n3+n*xydim*maxufe-1 ) = x n3 = n2 + n*vdim*maxvfe c c zzd( n4 : n4+n*xydim*maxvfe-1 ) = y n4 = n3 + n*xydim*maxufe c c zzd( n5 : n5+uedim*maxufe-1 ) = ue (i.e. the endpoints for U ) n5 = n4 + n*xydim*maxvfe c c zzd( n6 : n6+vedim*maxvfe-1 ) = ve (i.e. the endpoints for V ) n6 = n5 + uedim*maxufe c c zzd( n7 : n7+xydim*maxufe-1 ) = xe (i.e. the endpoints for x ) n7 = n6 + vedim*maxvfe c c zzd( n8 : n8+xydim*maxvfe-1 ) = ye (i.e. the endpoints for y ) n8 = n7 + xydim*maxufe c c zzd( n9 : n9+udim-1 ) = uplus (i.e. the upper bounds for u ) n9 = n8 + xydim*maxvfe c c zzd( n10 : n10+udim-1 ) = uminus (i.e. the lower bounds for u ) n10 = n9 + udim c c zzd( n11 : n11+uedim-1 ) = ueplus (i.e. the upper bounds for ue ) n11 = n10 + udim c c zzd( n12 : n12+uedim-1 ) = uemnus (i.e. the lower bounds for ue ) n12 = n11 + uedim c c zzd( n13 : n13+vdim-1 ) = vplus (i.e. the upper bounds for v ) n13 = n12 + uedim c c zzd( n14 : n14+vdim-1 ) = vminus (i.e. the lower bounds for v ) n14 = n13 + vdim c c zzd( n15 : n15+vedim-1 ) = veplus (i.e. the upper bounds for ve ) n15 = n14 + vdim c c zzd( n16 : n16+vedim-1 ) = veminus (i.e. the lower bounds for ve ) n16 = n15 + vedim c c zzd( n17 : n17+udim-1 ) = pmat n17 = n16 + vedim c c zzd( n18 : n18+udim*n-1 ) = pvec n18 = n17 + udim c c zzd( n19 : n19+vdim-1 ) = qmat n19 = n18 + udim*n c c zzd( n20 : n20+vdim*n-1 ) = qvec n20 = n19 + vdim c c zzd( n21 : n21+vdim*udim-1 ) = dmat n21 = n20 + vdim*n c c zzd( n22 : n22+uedim-1 ) = pemat n22 = n21 + vdim*udim c c zzd( n23 : n23+uedim-1 ) = pevec n23 = n22 + uedim c c zzd( n24 : n24+vedim-1 ) = qemat n24 = n23 + uedim c c zzd( n25 : n25+vedim-1 ) = qevec n25 = n24 + vedim c c zzd( n26 : n26+xydim*xydim-1 ) = amat n26 = n25 + vedim c c zzd( n27 : n27+xydim*udim-1 ) = bmat n27 = n26 + xydim*xydim c c zzd( n28 : n28+xydim-1 ) = bvec n28 = n27 + xydim*udim c c zzd( n29 : n29+xydim*vdim-1 ) = cmat n29 = n28 + xydim c c zzd( n30 : n30+xydim-1 ) = cvec n30 = n29 + xydim*vdim c c zzd( n31 : n31+xydim*uedim-1 ) = bemat n31 = n30 + xydim c c zzd( n32 : n32+xydim-1 ) = bevec n32 = n31 + xydim*uedim c c zzd( n33 : n33+xydim*vedim-1 ) = cemat n33 = n32 + xydim c c zzd( n34 : n34+xydim-1 ) = cevec n34 = n33 + xydim*vedim c c zzd( n35 : n35+maxufe)-1 ) = xi (i.e coeff's giving cvx combo of u's) n35 = n34 + xydim c c zzd( n36 : n36+maxvfe-1 )= eta (i.e coeff's giving cvx combo of v's) n36 = n35 + maxufe c c zzd( n37 : n37+ udim -1 ) = pd (augmented version of pmat) n37 = n36 + maxvfe c c zzd( n38 : n38+ n*udim -1 ) = pv (augmented version of pvec) n38 = n37 + udim c c zzd( n39 : n39+ vdim -1 ) = qd (augmented version of qmat) n39 = n38 + n*udim c c zzd( n40 : n40+ n*vdim -1 ) = qv (augmented version of qvec) n40 = n39 + vdim c c zzd( n41 : n41+ uedim -1 ) = ped (augmented version of pemat) n41 = n40 + n*vdim c c zzd( n42 : n42+ uedim -1 ) = pev (augmented version of pevec) n42 = n41 + uedim c c zzd( n43 : n43+ vedim -1 ) = qed (augmented version of qemat) n43 = n42 + uedim c c zzd( n44 : n44+ vedim -1 ) = qev (augmented version of qevec) n44 = n43 + vedim c c zzd( n45 : n45+ maxufe -1 ) = alpha (array of primal function values) n45 = n44 + vedim c c zzd( n46 : n46+ maxvfe -1 ) = beta (array of dual function values) n46 = n45 + maxufe c c LNSIZE is amount of space needed for LNSRCH lnsize = 8*max0(n*vdim+vedim,n*udim+uedim) c SADSIZ is amount of space needed for SADDLE sadsiz = 0 c sadsiz = 43 + 31*maxufe + 30*maxvfe + 8*maxufe*maxvfe c $ + 4*maxufe*maxufe + 5*maxvfe*maxvfe c c zzd( n47 : n47+ max(lnsize,sadsiz)-1 ) = work n47 = n46 + maxvfe c ncheck = n47 + max0(lnsize,sadsiz) c if ( ncheck .gt. nzdlen ) then open(unit=17,file=ofile) write(17,'("ERROR IN PROBLEM-SIZE FOR : ",a)') name write(17,'("''DYNFGM'' REQUIRES MORE MEMORY THAN ",$)') write(17,'("ALLOCATED !!")') write(17,'("AMOUNT OF ADDITIONAL DOUBLE-PRECISION MEMORY",$)') write(17,'(" NEEDED IS: ",$)') write(17,*) ncheck-nzdlen close(unit=17) stop endif c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c do 5 i=n35,ncheck zzd(i) = 0.0d0 5 continue c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c if ( .not. initu ) then c CALCULATE THE INITIAL GUESS FOR 'u' do 15 i=1,udim if ( dabs(zzd(n9+i-1)) .lt. dabs(zzd(n10-1+i)) ) then dummy = zzd(n9+i-1) else dummy = zzd(n10-1+i) endif do 10 j=1,n 10 zzd(n1-1+(j-1)*udim+i) = dummy 15 continue do 17 i=1,uedim if ( dabs(zzd(n11-1+i)) .lt. dabs(zzd(n12-1+i)) ) then zzd(n5+i-1) = zzd(n11-1+i) else zzd(n5+i-1) = zzd(n12-1+i) endif 17 continue endif c c if ( .not. initv ) then c CALCULATE THE INITIAL GUESS FOR 'v' do 25 i=1,vdim if ( dabs(zzd(n13-1+i)) .lt. dabs(zzd(n14-1+i)) ) then dummy = zzd(n13-1+i) else dummy = zzd(n14-1+i) endif do 20 j=1,n 20 zzd(n2-1+(j-1)*vdim+i) = dummy 25 continue do 27 i=1,vedim if ( dabs(zzd(n15-1+i)) .lt. dabs(zzd(n16-1+i)) ) then zzd(n6+i-1) = zzd(n15-1+i) else zzd(n6+i-1) = zzd(n16-1+i) endif 27 continue endif c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c propt(1,5) = 6 propt(2,5) = 12 propt(3,5) = 13 do 191 j=1,3 out = propt(j,5) if ( propt(j,1) .gt. 0) then if ( j .eq. 2 ) then open (unit=out,file=ofile) else if ( j .eq. 3 ) then i=index(ofile,' ') if ( i .eq. 0 ) then rsfile = ofile//'.res' else rsfile = ofile(1:i-1)//'.res' endif open (unit=out,file=rsfile) endif endif if ( propt(j,1) .gt. 1) then if ( mod(propt(j,3),2) .eq. 1) then c I.E. IF propt(j,3) is ODD ! if ( propt(j,2) .ge. 2 ) then write(out,'("Name of Problem : ",a,/)') name write(out,'("Dimensions :")') write(out,'(" primal vector = ",i6)') udim write(out,'(" primal endpoint = ",i6)') uedim write(out,'(" dual vector = ",i6)') vdim write(out,'(" dual endpoint = ",i6)') vedim write(out,'(" trajectories = ",i6)') xydim write(out,'(" number of time periods = ",i6,/)') n write(out,'("Stopping Criteria :")') write(out,'(" normalized duality gap = ",D9.2)') epsmin write(out,'(" number of iterations :")') write(out,'(" outer loop = ",i6)') mumax write(out,'(" inner loop = ",i6,/)') numax write(out,'("Memory Allocation :")') write(out,'(" total memory available = ",i8)') nzdlen write(out,'(" amount of memory used = ",i8)') ncheck else write(out,'(a)') name write(out,*) udim,uedim,vdim,vedim,xydim,n write(out,*) epsmin,mumax,numax endif endif endif 191 continue c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c call prox(zzd(n1),zzd(n2),zzd(n3),zzd(n4),zzd(n5),zzd(n6), $ zzd(n7),zzd(n8),zzd(n9),zzd(n10),zzd(n11),zzd(n12), $ zzd(n13),zzd(n14),zzd(n15),zzd(n16),zzd(n17),zzd(n18), $ zzd(n19),zzd(n20),zzd(n21),zzd(n22),zzd(n23),zzd(n24), $ zzd(n25),zzd(n26),zzd(n27),zzd(n28),zzd(n29),zzd(n30), $ zzd(n31),zzd(n32),zzd(n33),zzd(n34),zzd(n35),zzd(n36), $ zzd(n37),zzd(n38),zzd(n39),zzd(n40),zzd(n41),zzd(n42), $ zzd(n43),zzd(n44),zzd(n45),zzd(n46),dval, $ zzd(n47), nzdlen-n47 ,udim,uedim,vdim,vedim,xydim,n, $ maxufe,maxvfe,mumax,numax,epsmin) c c======================================================================= c return c END OF routine 'DYNFGM' end :::::::::::::: pfgval.f :::::::::::::: subroutine fgval(uopt,vopt,ueopt,veopt,alpha,beta, $ u,v,x,y,ue,ve,xe,ye, $ uupr,ulow,vlow,vupr,ueupr,uelow,veupr,velow, $ pm,pv,qm,qv,dm,pme,pve,qme,qve, $ bm,bv,cm,cv,bme,bve,cme,cve, $ ndu,ndue,nxy,ndv,ndve,ntime) c***************************************************************************** c c Computes 'vopt', 'uopt', 'beta' and 'alpha' from 'u', 'v', 'x' and 'y' by c c N 1 c alpha = sum { p *u(t) + -u(t)*Pu(t) - c*x(t-1) } c t=1 t 2 c c 1 c + p *u(0) + -u(0)*P u(0) - c *x(N) c e 2 e e c c +--- ---+ c maximum | / v(t)*[q - Cx(t-1) - Du(t)]\ | c over | N \ t / | c / v-<=v(t)<=v+ \ | sum | 1 | | c \ for all t=1,N / | t=1 / - - v(t)*Qv(t) \ | c + | and | | \ 2 / | c /v-<=v(N+1)<=v+ \ | 1 | c \ e e / | + v(N+1)*[q - C x(N)] - - v(N+1)*Q v(N+1) | c | e e 2 e | c +--- ---+ c c where 'vopt' is chosen to be an element v which gives the max defining alpha, c c N 1 c beta = sum { q *v(t) - -v(t)*Qv(t) - b*y(t+1) } c t=1 t 2 c c 1 c + q *v(N+1) - -v(N+1)*Q v(N+1) - b *y(1) c e 2 e e c c +--- T T ---+ c | / u(t)*[p - B y(t+1) - D v(t)] \ | c | N \ t / | c minimum | sum | 1 | | c over | t=1 / + - u(t)*Pu(t) \ | c / u-<=u(t)<=u+ \ | \ 2 / | c \ for all t=1,N / | | c + < and > | T 1 T | c / u-= optimal v(t) c veopt <--> optimal v(N+1) c uopt(,t)<--> optimal u(t) c ueopt <--> optimal u(0) c beta <--> beta c alpha <--> alpha c vlow <--> v- c vupr <--> v+ c velow <--> ve- c veupr <--> ve+ c ulow <--> u- c uupr <--> u+ c uelow <--> ue- c ueupr <--> ue+ c qv(,t) <--> q c t c pv(,t) <--> p c t c cm <--> C c bm <--> B c u(,t) <--> u(t) for t=1,N c ue <--> u(0) c v(,t) <--> v(t) for t=1,N c ve <--> v(N+1) c x(t) <--> x(t) for t=1,N c xe <--> x(0) c y(t) <--> y(t) for t=1,N c ye <--> y(N+1) c dm <--> D c qm <--> Q c pm <--> P c qve <--> qe c pve <--> pe c cme <--> Ce c bme <--> Be c qme <--> Qe c pme <--> Pe c ntime <--> N c c----------------------------------------------------------------------- c/////////////////////////////////////////////////////////////////////// c----------------------------------------------------------------------- c VARIABLES PASSED THRU SUBROUTINE CALL c c dimensions: integer ndu,ndue,nxy,ndv,ndve,ntime c c input variables: double precision u(ndu,ntime),v(ndv,ntime),x(nxy,ntime), $ y(nxy,ntime),ue(ndue),ve(ndue),xe(nxy),ye(nxy) c c data: double precision uupr(ndu),ulow(ndu),vlow(ndv),vupr(ndv), $ ueupr(ndue),uelow(ndue),veupr(ndve),velow(ndve), $ pm(ndu),pv(ndu,ntime), $ qm(ndv),qv(ndv,ntime),dm(ndv,ndu), $ pme(ndue),pve(ndue),qme(ndve),qve(ndve), $ bm(nxy,ndu),bv(nxy),cm(ndv,nxy),cv(nxy), $ bme(nxy,ndue),bve(nxy),cme(ndve,nxy),cve(nxy) c c output variables: double precision uopt(ndu,ntime),vopt(ndv,ntime), $ ueopt(ndue),veopt(ndve),alpha,beta c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLES c double precision asum,bsum,sum,sum1 c c----------------------------------------------------------------------- c/////////////////////////////////////////////////////////////////////// c----------------------------------------------------------------------- c c COMPUTE alpha AND vopt,veopt c asum = 0.0d0 c do 10 i=1,ndu asum = asum + pv(i,1)*u(i,1) + pm(i)*u(i,1)*u(i,1)/2.0d0 10 continue c do 15 i=1,nxy 15 asum = asum - cv(i)*xe(i) c c do 100 i=1,ndv c sum = qv(i,1) do 110 j=1,nxy 110 sum = sum - cm(i,j)*xe(j) do 120 j=1,ndu 120 sum = sum - dm(i,j)*u(j,1) c if ( sum .le. vlow(i)*qm(i) ) then sum1 = vlow(i) else if ( sum .gt. vupr(i)*qm(i) ) then sum1 = vupr(i) else c ! qm(i)<>0.0 and v-<= sum/qm(i) <=v+ sum1 = sum/qm(i) endif c vopt(i,1) = sum1 c asum = asum + sum1*sum - qm(i)*sum1*sum1/2.0d0 c 100 continue c c - - - - - - - - - - c do 250 k=2,ntime do 201 i=1,ndu asum = asum + pv(i,k)*u(i,k) + pm(i)*u(i,k)*u(i,k)/2.0d0 201 continue c do 205 i=1,nxy 205 asum = asum - cv(i)*x(i,k-1) c c do 200 i=1,ndv c sum = qv(i,k) do 210 j=1,nxy 210 sum = sum - cm(i,j)*x(j,k-1) do 220 j=1,ndu 220 sum = sum - dm(i,j)*u(j,k) c if ( sum .le. vlow(i)*qm(i) ) then sum1 = vlow(i) else if ( sum .gt. vupr(i)*qm(i) ) then sum1 = vupr(i) else c ! qm(i)<>0.0 and v-<= sum/qm(i) <=v+ sum1 = sum/qm(i) endif c vopt(i,k) = sum1 c asum = asum + sum1*sum - qm(i)*sum1*sum1/2.0d0 c 200 continue c 250 continue c c - - - - - - - - - - c do 301 i=1,ndue asum = asum + pve(i)*ue(i) + pme(i)*ue(i)*ue(i)/2.0d0 301 continue c do 305 i=1,nxy 305 asum = asum - cve(i)*x(i,ntime) c c do 300 i=1,ndve c c sum = qve(i) do 310 j=1,nxy 310 sum = sum - cme(i,j)*x(j,ntime) c if ( sum .le. velow(i)*qme(i) ) then sum1 = velow(i) else if ( sum .gt. veupr(i)*qme(i) ) then sum1 = veupr(i) else c ! qme(i)<>0.0 and ve-<= sum/qme(i) <=ve+ sum1 = sum/qme(i) endif c veopt(i) = sum1 c asum = asum + sum1*sum - qme(i)*sum1*sum1/2.0d0 c 300 continue c alpha = asum c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c COMPUTE beta AND uopt,ueopt c bsum = 0.0d0 c do 401 i=1,ndv bsum = bsum + qv(i,ntime)*v(i,ntime) bsum = bsum - qm(i)*v(i,ntime)*v(i,ntime)/2.0d0 401 continue c do 405 i=1,nxy 405 bsum = bsum - bv(i)*ye(i) c c do 400 i=1,ndu c sum = -pv(i,ntime) do 410 j=1,nxy 410 sum = sum + bm(j,i)*ye(j) do 420 j=1,ndv 420 sum = sum + dm(j,i)*v(j,ntime) c if ( sum .le. ulow(i)*pm(i) ) then sum1 = ulow(i) else if ( sum .gt. uupr(i)*pm(i) ) then sum1 = uupr(i) else c ! pm(i)<>0.0 and u-<= sum/pm(i) <=u+ sum1 = sum/pm(i) endif c uopt(i,ntime) = sum1 c bsum = bsum - sum1*sum + pm(i)*sum1*sum1/2.0d0 c 400 continue c c - - - - - - - - - - c do 550 k=1,ntime-1 do 501 i=1,ndv bsum = bsum + qv(i,k)*v(i,k) - qm(i)*v(i,k)*v(i,k)/2.0d0 501 continue c do 505 i=1,nxy 505 bsum = bsum - bv(i)*y(i,k+1) c c do 500 i=1,ndu c sum = -pv(i,k) do 510 j=1,nxy 510 sum = sum + bm(j,i)*y(j,k+1) do 520 j=1,ndv 520 sum = sum + dm(j,i)*v(j,k) c if ( sum .le. ulow(i)*pm(i) ) then sum1 = ulow(i) else if ( sum .gt. uupr(i)*pm(i) ) then sum1 = uupr(i) else c ! pm(i)<>0.0 and u-<= sum/pm(i) <=u+ sum1 = sum/pm(i) endif c uopt(i,k) = sum1 c bsum = bsum - sum1*sum + pm(i)*sum1*sum1/2.0d0 c 500 continue c 550 continue c c - - - - - - - - - - c do 601 i=1,ndve bsum = bsum + qve(i)*ve(i) - qme(i)*ve(i)*ve(i)/2.0d0 601 continue c do 605 i=1,nxy 605 bsum = bsum - bve(i)*y(i,1) c c do 600 i=1,ndue c sum = -pve(i) do 610 j=1,nxy 610 sum = sum + bme(j,i)*y(j,1) c if ( sum .le. uelow(i)*pme(i) ) then sum1 = uelow(i) else if ( sum .gt. ueupr(i)*pme(i) ) then sum1 = ueupr(i) else c ! pme(i)<>0.0 and ue-<= sum/pme(i) <=ue+ sum1 = sum/pme(i) endif c ueopt(i) = sum1 c bsum = bsum - sum1*sum + pme(i)*sum1*sum1/2.0d0 c 600 continue c beta = bsum c c======================================================================= c return c END OF subroutine 'FGVAL' end :::::::::::::: pgdt7.f :::::::::::::: Program gdt external rands, randu c PARAMETER ( mudim = 20, : muedim = 20, : mxydim = 20, : mvdim = 20, : mvedim = 20 ) 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 INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM C DOUBLE PRECISION UUPR(mUDIM) DOUBLE PRECISION UEUPR(mUEDIM) DOUBLE PRECISION VUPR(mVDIM) DOUBLE PRECISION VEUPR(mVEDIM) C C c======================================================================= c integer numax,n,propt(3,5),mode c double precision deps c character*60 name,ifile,ofile c logical initu,initv logical cont c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLE c double precision wt character*90 data integer lpropt(12), start c c======================================================================= c read(*,'(30x,a)') name read(name,'(9x,a)') data read(data,*) start print*, 'start=', start c read(*,'(30x,a)') data read(data,*) cont c read(*,'(30x,a)') data read(data,*) udim c read(*,'(30x,a)') data read(data,*) uedim c read(*,'(30x,a)') data read(data,*) xydim c read(*,'(30x,a)') data read(data,*) vdim c read(*,'(30x,a)') data read(data,*) vedim c read(*,'(30x,a)') data read(data,*) n c read(*,'(30x,a)') data read(data,*) mumax c read(*,'(30x,a)') data read(data,*) numax c read(*,'(30x,a)') data read(data,*) deps c read(*,'(30x,a)') data read(data,*) mode c read(*,'(30x,a)') data read(data,*) initu c read(*,'(30x,a)') data read(data,*) initv c read(*,'(30x,a)') ifile c read(*,'(30x,a)') ofile c read(*,'(30x,a)') data i=index(data,'(') read(data(i:),1000) lpropt 1000 format( 3(1x,4(i1,1x)) ) kk=0 do 1110 j=1,3 do 1110 k=1,4 kk=kk+1 propt(j,k) = lpropt(kk) 1110 continue c wt=rands(start) c c======================================================================= c open(3, file = ifile) c c generate uupr: do 100 i = 1, udim wt = 10. * randu() - 5. uupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 100 continue c c generate ulow: do 101 i = 1, udim wt = uupr(i) - 4.5 * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 101 continue c c generate ueupr: do 103 i = 1, uedim wt = 10. * randu() - 5. ueupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 103 continue c c generate uelow: do 104 i = 1, uedim wt = ueupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 104 continue c c generate vupr: do 105 i = 1, vdim wt = 10. * randu() - 5. vupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 105 continue c c generate vlow: do 106 i = 1, vdim wt = vupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 106 continue c c generate veupr: do 107 i = 1, vedim wt = 10. * randu() - 5. veupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 107 continue c c generate velow: do 108 i = 1, vedim wt = veupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 108 continue c c generate pmat: do 110 i = 1, udim wt = 2. * randu() + 1.d0 if (wt .lt. 2.) wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 110 continue c c generate pvec: do 111 i = 1, udim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 111 continue c c generate qmat: do 112 i = 1, vdim wt = 2. * randu() + 1.d0 if (wt .lt. 2.) wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 112 continue c c generate qvec: do 113 i = 1, vdim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 113 continue c c generate dmat: do 120 i = 1, vdim do 120 j = 1, udim wt = 20. * randu() - 10. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 120 continue c c generate pemat: do 130 i = 1, uedim wt = 2. * randu() + 1.d0 if (wt .lt. 2.) wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 130 continue c c generate pevec: do 131 i = 1, uedim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 131 continue c c generate qemat: do 132 i = 1, vedim wt = 2. * randu() + 1.d0 if (wt .lt. 2.) wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 132 continue c c generate qevec: do 133 i = 1, vedim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 133 continue c c generate amat: do 140 i = 1, xydim do 140 j = 1, xydim wt = 1.8 * randu() - 0.9 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 140 continue c c generate bmat: do 141 i = 1, xydim do 141 j = 1, udim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 141 continue c c generate bvec: do 142 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 142 continue c c generate cmat: do 143 i = 1, vdim do 143 j = 1, xydim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 143 continue c c generate cvec: do 144 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 144 continue c c c generate bemat: do 145 i = 1, xydim do 145 j = 1, uedim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 145 continue c c generate bevec: do 146 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 146 continue c c generate cemat: do 147 i = 1, vedim do 147 j = 1, xydim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 147 continue c c generate cevec: do 148 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 148 continue stop end ********************************************************************** real function rands(start) integer start real randu integer l,c,m parameter(l=29,c=217,m=1024) integer seed save seed data seed /0/ seed=mod(abs(start),m) entry randu( ) seed=mod(seed*l+c,m) randu=real(seed)/m end ************************************************************************* :::::::::::::: pgdt8.f :::::::::::::: Program gdt external rands, randu c PARAMETER ( mudim = 20, : muedim = 20, : mxydim = 20, : mvdim = 20, : mvedim = 20 ) 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 INTEGER UDIM,UEDIM,VDIM,VEDIM,XYDIM C DOUBLE PRECISION UUPR(mUDIM) DOUBLE PRECISION UEUPR(mUEDIM) DOUBLE PRECISION VUPR(mVDIM) DOUBLE PRECISION VEUPR(mVEDIM) C C c======================================================================= c integer numax,n,propt(3,5),mode c double precision deps c character*60 name,ifile,ofile c logical initu,initv logical cont c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLE c double precision wt character*90 data integer lpropt(12), start c c======================================================================= c read(*,'(30x,a)') name read(name,'(9x,a)') data read(data,*) start print*, 'start=', start c read(*,'(30x,a)') data read(data,*) cont c read(*,'(30x,a)') data read(data,*) udim c read(*,'(30x,a)') data read(data,*) uedim c read(*,'(30x,a)') data read(data,*) xydim c read(*,'(30x,a)') data read(data,*) vdim c read(*,'(30x,a)') data read(data,*) vedim c read(*,'(30x,a)') data read(data,*) n c read(*,'(30x,a)') data read(data,*) mumax c read(*,'(30x,a)') data read(data,*) numax c read(*,'(30x,a)') data read(data,*) deps c read(*,'(30x,a)') data read(data,*) mode c read(*,'(30x,a)') data read(data,*) initu c read(*,'(30x,a)') data read(data,*) initv c read(*,'(30x,a)') ifile c read(*,'(30x,a)') ofile c read(*,'(30x,a)') data i=index(data,'(') read(data(i:),1000) lpropt 1000 format( 3(1x,4(i1,1x)) ) kk=0 do 1110 j=1,3 do 1110 k=1,4 kk=kk+1 propt(j,k) = lpropt(kk) 1110 continue c wt=rands(start) c c======================================================================= c open(3, file = ifile) c c generate uupr: do 100 i = 1, udim wt = 10. * randu() - 5. uupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 100 continue c c generate ulow: do 101 i = 1, udim wt = uupr(i) - 4.5 * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 101 continue c c generate ueupr: do 103 i = 1, uedim wt = 10. * randu() - 5. ueupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 103 continue c c generate uelow: do 104 i = 1, uedim wt = ueupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 104 continue c c generate vupr: do 105 i = 1, vdim wt = 10. * randu() - 5. vupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 105 continue c c generate vlow: do 106 i = 1, vdim wt = vupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 106 continue c c generate veupr: do 107 i = 1, vedim wt = 10. * randu() - 5. veupr(i) = wt if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 107 continue c c generate velow: do 108 i = 1, vedim wt = veupr(i) - 4. * randu() if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 108 continue c c generate pmat: do 110 i = 1, udim wt = 2. * randu() if (wt .lt. 1.) wt = 0. wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 110 continue c c generate pvec: do 111 i = 1, udim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 111 continue c c generate qmat: do 112 i = 1, vdim wt = 2. * randu() if (wt .lt. 1.) wt = 0. wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 112 continue c c generate qvec: do 113 i = 1, vdim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 113 continue c c generate dmat: do 120 i = 1, vdim do 120 j = 1, udim wt = 20. * randu() - 10. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 120 continue c c generate pemat: do 130 i = 1, uedim wt = 2. * randu() if (wt .lt. 1.) wt = 0. wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 130 continue c c generate pevec: do 131 i = 1, uedim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 131 continue c c generate qemat: do 132 i = 1, vedim wt = 2. * randu() if (wt .lt. 1.) wt = 0. wt = 0. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 132 continue c c generate qevec: do 133 i = 1, vedim wt = 4. * randu() - 2. if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 133 continue c c generate amat: do 140 i = 1, xydim do 140 j = 1, xydim wt = 1.8 * randu() - 0.9 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 140 continue c c generate bmat: do 141 i = 1, xydim do 141 j = 1, udim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 141 continue c c generate bvec: do 142 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 142 continue c c generate cmat: do 143 i = 1, vdim do 143 j = 1, xydim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 143 continue c c generate cvec: do 144 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 144 continue c c c generate bemat: do 145 i = 1, xydim do 145 j = 1, uedim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 145 continue c c generate bevec: do 146 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 146 continue c c generate cemat: do 147 i = 1, vedim do 147 j = 1, xydim wt = .4 * randu() - .2 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 147 continue c c generate cevec: do 148 i = 1, xydim wt = randu() - .5 if (wt .lt. 0.) then write (3, '(30x, e17.10)') wt else write (3, '(30x, e16.10)') wt endif 148 continue stop end ********************************************************************** real function rands(start) integer start real randu integer l,c,m parameter(l=29,c=217,m=1024) integer seed save seed data seed /0/ seed=mod(abs(start),m) entry randu( ) seed=mod(seed*l+c,m) randu=real(seed)/m end ************************************************************************* :::::::::::::: pinit.f :::::::::::::: c The following code is a slightly modified version of Stephen E. Wright's c original code. SUBROUTINE INIT(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ 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,NUFE,NVFE, $ GRDU,GRDV,CGU,CGV, $ UDIM,VDIM,XYDIM,UEDIM,VEDIM,TIME,MAXUFE,MAXVFE, $ WORK,LENWRK,QUIT) C*********************************************************************** C This routine is required by the DYNFGM program to do the following: C a. Build the initial finite-element sets and compute the C corresponding trajectories. In particular, it should find C the vectors vdot in argmax{L(u^,v):v in V} and C udot in argmin{L(u,v^):u in U} , where (u^,v^) is the initial C "best" guess inputted to the algorithm. C b. Return (i) the positions in the U and V arrays of the initial C guesses (u^,v^) and also of the elements (udot,vdot); C (ii) the function values for various elements of the C U and V arrays, including AT LEAST the function values C which correspond to u^ and v^, respectively; C (iii) the total number of elements currently in the C U and V arrays. C c. 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 INIT: C U,V,UE,VE -- finite sets of decision vectors; C X,Y,XE,YE -- the corresponding trajectories; 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 BESTU,BESTV -- pointers indicating the positions of the "best" C elements in the finite sets; C ALTU,ALTV -- pointers indicating the positions of the "alternative" C elements in the finite sets which give the max and min C defining the function values for the "best" elements. C QUIT -- Flag whose value is .FALSE. when INIT is called. If it C is decided that the algorithm should be immediately C terminated (due to some criterion specified within INIT C by the user), then INIT should change QUIT to be .TRUE. C and print out any appropriate information. C NOTE: Upon entering INIT, the first element of the U array (i.e. the C sub-array of dimension UDIM by TIME which is pointed to by C the expression U(1,1,1)) contains the initial primal guess and C the first element of the UE array contains the corresponding C endpoint control. If this element is left unchanged, BESTU C should be set equal to 1 before exiting. Otherwise this element C (and its endpoint control) must be copied to another position C in the set and BESTU given the value of the new position. 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 ALPHA(MAXUFE),BETA(MAXVFE) INTEGER BESTU,BESTV,ALTU,ALTV INTEGER GRDU,GRDV,CGU,CGV INTEGER NUFE,NVFE LOGICAL QUIT C C data (values not to be changed in INIT) 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 dimensions for the above arrays C C workspace for use in INIT and other routines called during the C algorithm (NOTE: values placed in this array cannot be saved between C successive calls to INIT.) 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 Stephen E. Wright IN August 1988 C --- --------- ---- --- ------- -- ================= -- =========== C AT IIASA, Laxenburg Austria. 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 3. Variable to be used in LINE, and is initialized in INIT DOUBLE PRECISION EPS COMMON /EPSLIN/ EPS C C Local variables: DOUBLE PRECISION EPSNRM INTEGER OUT C----------------------------------------------------------------------- C BESTU = 1 BESTV = 1 ALTU = 2 ALTV = 2 GRDU = 3 GRDV = 3 CGU = 4 CGV = 4 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Call routine to calculate the trajectories X and Y for U , V C 1 1 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Compute dual element ALTV and function value ALPHA for U = BESTU C 1 C and ALTU and BETA for V = BESTV . C 1 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 trajectories for ALTU and ALTV C CALL TRAJ(U(1,1,ALTU),UE(1,ALTU),V(1,1,ALTV),VE(1,ALTV), $ X(1,1,ALTU),XE(1,ALTU),Y(1,1,ALTV),YE(1,ALTV), $ AMAT,BMAT,CMAT,BVEC,CVEC,BEMAT,CEMAT,BEVEC,CEVEC, $ UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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 4001 FORMAT(4x,'INNER 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,'(4x,"INNER ITERATION 0:")') ENDIF IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(/," u vector:")') WRITE(OUT,'(6(E12.2))') ((U(J,K,BESTU),J=1,UDIM),K=1,TIME) WRITE(OUT,'(" ue vector:")') WRITE(OUT,'(6(E12.2))') (UE(J,BESTU), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" x vector:")') WRITE(OUT,'(6(E12.2))') ((X(J,K,BESTU),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(" xe vector:")') WRITE(OUT,'(6(E12.2))') (XE(J,BESTU), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(34x,"primal objective value:",G20.14)') ALPHA(BESTU) IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(/," v vector:")') WRITE(OUT,'(6(E12.2))') ((V(I,K,BESTV),I=1,VDIM),K=1,TIME) WRITE(OUT,'(" ve vector:")') WRITE(OUT,'(6(E12.2))') (VE(J,BESTV), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(" y vector:")') WRITE(OUT,'(6(E12.2))') ((Y(J,K,BESTV),J=1,XYDIM), : K=1,TIME) WRITE(OUT,'(" ye vector:")') WRITE(OUT,'(6(E12.2))') (YE(J,BESTV), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT, : '(36x,"dual objective value:",G20.14)') BETA(BESTV) 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,BESTU),J=1,UDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (UE(J,BESTU), J=1,UEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((X(J,K,BESTU),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(G26.19)') (XE(J,BESTU), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2) WRITE(OUT,'(G26.19)') ALPHA(BESTU) IF ( MOD(PROPT(JJ,4),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((V(I,K,BESTV),I=1,VDIM), K=1,TIME) WRITE(OUT,'(G26.19)') (VE(J,BESTV), J=1,VEDIM) IF ( MOD(PROPT(JJ,2),2) .EQ. 1 ) THEN WRITE(OUT,'(G26.19)') ((Y(J,K,BESTV),J=1,XYDIM),K=1,TIME) WRITE(OUT,'(6(E12.2))') (YE(J,BESTV), J=1,XYDIM) ENDIF ENDIF IF ( PROPT(JJ,3) .GE. 2 ) WRITE(OUT,'(G26.19)') BETA(BESTV) ENDIF ENDIF 191 CONTINUE C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C Initial size of finite-element sets NUFE = 2 NVFE = 2 C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C End of subroutine 'INIT' RETURN END :::::::::::::: pinput.f :::::::::::::: c The following code is a slightly modified version of Stephen E. Wright's c original code. subroutine input(zzd,dval,udim,uedim,xydim,vdim,vedim,mumax, $ numax,deps,n,name,ofile,propt,initu,initv, $ mode,nzdlen) c c======================================================================= c VARIABLES PASSED THRU SUBROUTINE CALL c integer nzdlen double precision zzd(nzdlen),dval c integer udim,uedim,xydim,vdim,vedim integer mumax,numax,n,propt(3,5),mode double precision deps character*60 name,ifile,ofile logical initu,initv c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c LOCAL VARIABLES c double precision h integer nuinit,nvinit,nuenit,nvenit,nupls,numns,nuepls,nuemns, $ nvpls,nvmns,nvepls,nvemns,npm,npv,nqm,nqv,ndm,npem, $ npev,nqem,nqev,nam,nbm,ncm,ncv,nbem,nbev,ncem,ncev, $ nshm,nmhm,nwork,ncheck integer maxufe,maxvfe logical cont c c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c READ THE PARAMETERS CONTROLLING THE PROGRAM FROM THE STANDARD INPUT c call inctrl(udim,uedim,xydim,vdim,vedim,mumax,numax,n,propt, $ deps,name,ifile,ofile,initu,initv,mode,cont) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c COMPUTE INDICES TO PLACE DATA IN 'zzd' c maxufe = 16 maxvfe = 16 c nuinit = 1 nvinit = nuinit + n*udim*maxufe nuenit = nvinit + vdim*n*maxvfe + xydim*n*(maxufe+maxvfe) nvenit = nuenit + uedim*maxufe nupls = nvenit + vedim*maxufe + xydim*(maxufe+maxvfe) numns = nupls + udim nuepls = numns + udim nuemns = nuepls + uedim nvpls = nuemns + uedim nvmns = nvpls + vdim nvepls = nvmns + vdim nvemns = nvepls + vedim npm = nvemns + vedim npv = npm + udim nqm = npv + udim*n nqv = nqm + vdim ndm = nqv + vdim*n npem = ndm + vdim*udim npev = npem + uedim nqem = npev + uedim nqev = nqem + vedim nam = nqev + vedim nbm = nam + xydim*xydim nbv = nbm + xydim*udim ncm = nbv + xydim ncv = ncm + xydim*vdim nbem = ncv + xydim nbev = nbem + xydim*uedim ncem = nbev + xydim ncev = ncem + xydim*vedim nshm = ncev + xydim nmhm = nshm + xydim*xydim nwork = nmhm + xydim*xydim ncheck = nwork + 2*xydim*xydim + udim + vdim c if ( ncheck .gt. nzdlen ) then open(unit=17,file=ofile) write(17,'("ERROR IN PROBLEM-SIZE FOR : ",a)') name write(17,'("INPUT ROUTINE REQUIRES MORE MEMORY THAN ",$)') write(17,'("ALLOCATED !!")') write(17,'("AMOUNT OF ADDITIONAL DOUBLE-PRECISION MEMORY",$)') write(17,'(" NEEDED IS: ",$)') write(17,*) ncheck-nzdlen write(17,*) 'ncheck,nzdlen',ncheck,nzdlen close(unit=17) stop endif c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c READ IN THE VECTOR AND MATRIX DATA FROM THE FILE 'ifile' c call indata(zzd,nuinit,nvinit,nuenit,nvenit,nupls,numns,nuepls, $ nuemns,nvpls,nvmns,nvepls,nvemns,npem,npev,nqem, $ nqev,nbem,nbev,ncem,ncev,npm,npv,nqm,nqv,ndm,nam, $ nbm,nbv,ncm,ncv,udim,uedim,xydim,vdim,vedim,mumax, $ numax,n,nzdlen,ifile,initu,initv) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c if ( cont ) then c c COMPUTE THE MATRICES 'S , M ' FROM THE MATRIX ' A '. c h h h = 1.0d0/dble(n) call shexp(zzd(nshm),zzd(nmhm),h,zzd(nam),zzd(nwork),xydim) c c - - - - - - - - - - - - - - - - c c USE ' S , M ' TO CONVERT ' p, P, q, Q, D, A, b, B, c, C ' c h h c c INTO ' p , P , q , Q , D , A , b , B , c , C ' and 'd '. c h h h h h h h h h h h c call hcomp(zzd(npv),zzd(npm),zzd(nqv),zzd(nqm),zzd(ndm),dval, $ zzd(nam),zzd(nbv),zzd(nbm),zzd(ncv),zzd(ncm), h, $ zzd(nshm),zzd(nmhm),zzd(nwork),udim,vdim,xydim) CC This additive constant should be added for EACH time period. CC Steve's original manual (and codes) missed this point. dval = dval * dble(n) c c else c The problem is a discrete-time problem, so the additive c constant should be zero. dval = 0.0d0 endif c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c COPY THE COMPONENTS OF P , p , Q , q TO P , p , Q , q FOR t=2,n. c 1 1 1 1 t t t t c do 100 k=1,n-1 do 10 i=0,udim-1 zzd(npv+k*udim+i) = zzd(npv+i) 10 continue do 20 i=0,vdim-1 zzd(nqv+k*vdim+i) = zzd(nqv+i) 20 continue 100 continue c======================================================================= c return c END OF subroutine 'INPUT' end subroutine inctrl(udim,uedim,xydim,vdim,vedim,mumax,numax, $ n,propt,deps,name,ifile,ofile,initu,initv, $ mode,cont) c c======================================================================= c c VARIABLES PASSED THRU SUBROUTINE CALL c integer udim,uedim,xydim,vdim,vedim integer mumax,numax,n,propt(3,5),mode c double precision deps c character*60 name,ifile,ofile c logical initu,initv logical cont c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLE c character*90 data integer lpropt(12) c c======================================================================= c read(*,'(30x,a)') name c read(*,'(30x,a)') data read(data,*) cont c read(*,'(30x,a)') data read(data,*) udim c read(*,'(30x,a)') data read(data,*) uedim c read(*,'(30x,a)') data read(data,*) xydim c read(*,'(30x,a)') data read(data,*) vdim c read(*,'(30x,a)') data read(data,*) vedim c read(*,'(30x,a)') data read(data,*) n c read(*,'(30x,a)') data read(data,*) mumax c read(*,'(30x,a)') data read(data,*) numax c read(*,'(30x,a)') data read(data,*) deps c read(*,'(30x,a)') data read(data,*) mode c read(*,'(30x,a)') data read(data,*) initu c read(*,'(30x,a)') data read(data,*) initv c read(*,'(30x,a)') ifile c read(*,'(30x,a)') ofile c read(*,'(30x,a)') data i=index(data,'(') read(data(i:),100) lpropt 100 format( 3(1x,4(i1,1x)) ) kk=0 do 10 j=1,3 do 10 k=1,4 kk=kk+1 propt(j,k) = lpropt(kk) 10 continue c c c======================================================================= c return c END OF subroutine 'INCTRL' end subroutine indata(zzd,nuinit,nvinit,nuenit,nvenit,nupls,numns, $ nuepls,nuemns,nvpls,nvmns,nvepls,nvemns,npem,npev, $ nqem,nqev,nbem,nbev,ncem,ncev,npm,npv,nqm,nqv,ndm, $ nam,nbm,nbv,ncm,ncv,udim,uedim,xydim,vdim,vedim, $ mumax,numax,n,nzdlen,ifile,initu,initv) c---------------------------------------- c VARIABLES PASSED THROUGH SUBROUTINE CALL: integer nzdlen double precision zzd(nzdlen) c integer udim,uedim,xydim,vdim,vedim,mumax,numax,n character*60 ifile logical initu,initv c integer nuinit,nvinit,nuenit,nvenit,nupls,numns,nuepls,nuemns, $ nvpls,nvmns,nvepls,nvemns,npem,npev,nqem,nqev, $ nbem,nbev,ncem,ncev,npm,npv,nqm,nqv,ndm,nam,nbm,nbv, $ ncm,ncv c---------------------------------------- c LOCAL VARIABLES: character*60 data c c----------------------------------------------------------------------- c/////////////////////////////////////////////////////////////////////// c----------------------------------------------------------------------- c open(unit=17,file=ifile) c c-------------------- if ( initu ) then c READ IN INITIAL GUESS FOR u c AS ( do j=1,n ( do i=1,udim u(i,j) ) ) do 10 i=0,n-1 do 10 j=nuinit,nuinit+udim-1 read(17,'(30x,a)') data read(data,*) zzd( i*udim + j ) 10 continue do 15 i=nuenit,nuenit+uedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 15 continue endif c c-------------------- if ( initv ) then c READ IN INITIAL GUESS FOR v c AS ( do j=1,n ( do i=1,vdim v(i,j) ) ) do 20 i=0,n-1 do 20 j=nvinit,nvinit+vdim-1 read(17,'(30x,a)') data read(data,*) zzd( i*vdim + j ) 20 continue do 25 i=nvenit,nvenit+vedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 25 continue endif c c-------------------- c READ IN uplus,uminus,ueplus,uemnus c do 30 i=nupls,nupls+udim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 30 continue do 33 i=numns,numns+udim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 33 continue do 36 i=nuepls,nuepls+uedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 36 continue do 39 i=nuemns,nuemns+uedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 39 continue c-------------------- c READ IN vplus,vminus,veplus,vemnus c do 40 i=nvpls,nvpls+vdim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 40 continue do 43 i=nvmns,nvmns+vdim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 43 continue do 46 i=nvepls,nvepls+vedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 46 continue do 49 i=nvemns,nvemns+vedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 49 continue c c-------------------- c READ IN P,p,Q,q,D c (NOTE: P and Q are assumed to be diagonal. The input for c them is therefore only the diagonal elements) c do 50 i=npm,npm+udim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 50 continue do 51 i=npv,npv+udim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 51 continue do 52 i=nqm,nqm+vdim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 52 continue do 53 i=nqv,nqv+vdim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 53 continue c c D is read in as: ( do i=1,vdim ( do j=1,udim read d(i,j) ) ) do 59 i=ndm,ndm+vdim-1 do 59 j=0,udim-1 read(17,'(30x,a)') data read(data,*) zzd( j*vdim + i ) 59 continue c c-------------------- c READ IN P ,p ,Q ,q c e e e e c (NOTE: P and Q are assumed to be diagonal, hence only c e e c the diagonal elements are to be given) c do 60 i=npem,npem+uedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 60 continue do 63 i=npev,npev+uedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 63 continue do 66 i=nqem,nqem+vedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 66 continue do 69 i=nqev,nqev+vedim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 69 continue c c-------------------- c READ IN A,B,b,C,c,B ,b ,C ,c c e e e e c c A is read in as: ( do i=1,xydim ( do j=1,xydim read a(i,j) ) do 70 i=nam,nam+xydim-1 do 70 j=0,xydim-1 read(17,'(30x,a)') data read(data,*) zzd( j*xydim + i ) 70 continue c c B is read in as: ( do i=1,xydim ( do j=1,udim read b(i,j) ) do 80 i=nbm,nbm+xydim-1 do 80 j=0,udim-1 read(17,'(30x,a)') data read(data,*) zzd( j*xydim + i ) 80 continue c do 83 i=nbv,nbv+xydim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 83 continue c c C is read in as: ( do i=1,vdim ( do j=1,xydim read c(i,j) ) do 86 i=ncm,ncm+vdim-1 do 86 j=0,xydim-1 read(17,'(30x,a)') data read(data,*) zzd( j*vdim + i ) 86 continue c do 87 i=ncv,ncv+xydim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 87 continue c c B is read in as: ( do i=1,xydim ( do j=1,uedim read b(i,j) ) c e e do 90 i=nbem,nbem+xydim-1 do 90 j=0,uedim-1 read(17,'(30x,a)') data read(data,*) zzd( j*xydim + i ) 90 continue c do 93 i=nbev,nbev+xydim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 93 continue c c C is read in as: ( do i=1,vedim ( do j=1,xydim read c(i,j) ) c e e do 96 i=ncem,ncem+vedim-1 do 96 j=0,xydim-1 read(17,'(30x,a)') data read(data,*) zzd( j*vedim + i ) 96 continue c do 99 i=ncev,ncev+xydim-1 read(17,'(30x,a)') data read(data,*) zzd(i) 99 continue c c-------------------- c close(unit=17) c c======================================================================= c return c END OF subroutine 'INDATA' end c************************************************************************ c* SUBROUTINE SHEXP * c* ---------- ===== * c* * c* This routine takes as input the matrix A, the scalar h and the dim- * c* ension n and returns the matrices S and M. The calculation per- * c* formed is: * c* * c* 2 3 4 2 * c* S = (h /2!)I + (h /3!)A +(h /4!)A + ... * c* * c* M = hI + AS * c* * c* The power series is calculated out to the fifteenth term, which * c* will give "near-" double precision accuracy for n<=4, h<=0.5, * c* and all entries of A in [-1.0,1.0]. (A,S,M are n-by-n) * c* * c* Note: the routine requires a double-precision array of length * c* 2*n*n for workspace, as well as S,M, and A. * c************************************************************************ c subroutine shexp(s,m,h,a,w,n) c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c VARIABLES PASSED THRU SUBROUTINE CALL c integer n double precision s(n,n),m(n,n),h,a(n,n),w(n,n,2) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLES c double precision mmm,sum1,cc integer i,j,ow,nw c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c 2 c INITIALIZE : S = (h /2!)I , W = I c 1 c do 10 i=1,n do 10 j=1,n s(i,j) = 0.0d0 w(i,j,1) = 0.0d0 10 continue mmm = 2.0d0 cc = h*h/mmm do 20 i=1,n w(i,i,1) = 1.0d0 s(i,i) = cc 20 continue ow = 1 nw = 2 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c COMPUTE S FROM A AND h c ii-1 c NOTE: In iteration ii the "old" W (W ) is A c ow c ii c and we compute the "new" W (W ) as W = A = AW . c nw nw ow c do 65 ii=1,15 mmm = mmm + 1.0d0 cc = cc * h/mmm do 60 i=1,n do 50 j=1,n sum1 = 0.0d0 do 40 k=1,n sum1 = sum1 + w(i,k,ow)*a(k,j) 40 continue w(i,j,nw) = sum1 s(i,j) = s(i,j) + cc*sum1 50 continue 60 continue j = ow ow = nw nw = ow 65 continue c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c COMPUTE M FROM A, S, h. c do 90 i=1,n do 80 j=1,n sum1 = 0.0d0 do 70 k=1,n 70 sum1 = sum1 + a(i,k)*s(k,j) m(i,j) = sum1 80 continue m(i,i) = m(i,i) + h 90 continue c c======================================================================= c return c END OF subroutine 'SHEXP' end c************************************************************************ c* SUBROUTINE HCOMP * c* * c* This routine takes as input the matrices S,M,P,Q,D,A,B,C, the * c* vectors p,q,b,c, the scalar h and the dimensions k,l,n and returns * c* the vectors p , q , b , c , and the matrices P , Q , D , A , B , C * c* h h h h h h h h h h * c* in place of the originals, as well as the "residual" constant d . * c* h * c* The transformation is as follows: * c* T T * c* p = hp - B S c , q = hq - CSb * c* h * c* T * c* P = hP , Q = hQ , D = hD + CSB , d = c Sb * c* h h h h * c* T * c* A = I + MA , B = MB , b = Mb , C = CM , c = M c, * c* h h h h h * c* * c* where I is the identity matrix. * c* * c* * c* Note: in the following routine, k=udim, l=vdim, n=xydim. * c************************************************************************ c///////////////////////////////////////////////////////////////////////// c// NOTE : THE MATRICES ' P, Q ' ARE ASSUMED TO BE DIAGONAL AND // c// // c// THEREFORE EACH IS GIVEN BY A 1-DIM'L ARRAY REPRESENTING // c// ITS DIAGONAL ELEMENTS // c///////////////////////////////////////////////////////////////////////// subroutine hcomp(pvec,pmat,qvec,qmat,dmat,dval,amat,bvec,bmat, $ cvec,cmat,h,s,m,work,udim,vdim,xydim) c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c VARIABLES PASSED THRU SUBROUTINE CALL c integer udim,xydim,vdim c double precision pvec(udim),qvec(vdim) double precision pmat(udim),qmat(vdim) double precision dmat(vdim,udim),dval double precision amat(xydim,xydim),bmat(xydim,udim) double precision bvec(xydim),cmat(vdim,xydim),cvec(xydim),h double precision s(xydim,xydim),m(xydim,xydim) double precision work(xydim+udim+vdim) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c LOCAL VARIABLES c double precision sum,sum1 c c======================================================================= c/////////////////////////////////////////////////////////////////////// c======================================================================= c c do 30 i=1,udim sum = 0.0d0 do 20 j=1,xydim sum1 = 0.0d0 do 10 k=1,xydim sum1 = sum1 + s(j,k)*bmat(k,i) 10 continue sum = sum + cvec(j)*sum1 20 continue pvec(i) = h*pvec(i) - sum pmat(i) = h*pmat(i) 30 continue c do 60 i=1,vdim sum = 0.0d0 do 50 j=1,xydim sum1 = 0.0d0 do 40 k=1,xydim sum1 = sum1 + s(k,j)*cmat(i,k) 40 continue sum = sum + bvec(j)*sum1 50 continue qvec(i) = h*qvec(i) - sum qmat(i) = h*qmat(i) 60 continue c do 100 ii=1,udim do 90 i=1,vdim sum = 0.0d0 do 80 j=1,xydim sum1 = 0.0d0 do 70 k=1,xydim sum1 = sum1 + s(k,j)*cmat(i,k) 70 continue sum = sum + bmat(j,ii)*sum1 80 continue dmat(i,ii) = h*dmat(i,ii) + sum 90 continue 100 continue c sum = 0.0d0 do 102 i=1,vdim sum1 = 0.0d0 do 101 j=1,udim sum1 = sum1 + s(i,j)*bvec(j) 101 continue sum = sum + cvec(i)*sum1 102 continue dval = -sum c c do 140 j=1,udim sum =0.0d0 do 110 i=1,xydim work(i) = bmat(i,j) sum = sum + m(1,i)*work(i) 110 continue bmat(1,j) = sum do 130 k=2,xydim sum = 0.0d0 do 120 i=1,xydim 120 sum = sum + m(k,i)*work(i) bmat(k,j) = sum 130 continue 140 continue c c do 145 j=1,xydim sum =0.0d0 do 115 i=1,xydim work(i) = amat(i,j) sum = sum + m(1,i)*work(i) 115 continue amat(1,j) = sum do 135 k=2,xydim sum = 0.0d0 do 125 i=1,xydim 125 sum = sum + m(k,i)*work(i) amat(k,j) = sum 135 continue amat(j,j) = amat(j,j) + 1.0d0 145 continue c c do 240 j=1,vdim sum =0.0d0 do 210 i=1,xydim work(i) = cmat(j,i) sum = sum + m(i,1)*work(i) 210 continue cmat(j,1) = sum do 230 k=2,xydim sum = 0.0d0 do 220 i=1,xydim 220 sum = sum + m(i,k)*work(i) cmat(j,k) = sum 230 continue 240 continue c c sum =0.0d0 do 215 i=1,xydim work(i) = cvec(i) sum = sum + m(i,1)*work(i) 215 continue cvec(1) = sum do 235 k=2,xydim sum = 0.0d0 do 225 i=1,xydim 225 sum = sum + m(i,k)*work(i) cvec(k) = sum 235 continue c c sum =0.0d0 do 315 i=1,xydim work(i) = bvec(i) sum = sum + m(1,i)*work(i) 315 continue bvec(1) = sum do 335 k=2,xydim sum = 0.0d0 do 325 i=1,xydim 325 sum = sum + m(k,i)*work(i) bvec(k) = sum 335 continue c c c======================================================================= c return c END OF subroutine 'HCOMP' end :::::::::::::: pline.f :::::::::::::: SUBROUTINE LINE(U,V,X,Y,UE,VE,XE,YE,ALPHA,BETA, $ OLDU,OLDV,UDIR,VDIR,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 C This routine is required by the DYNFGM program to do the following: C a. Perform a linesearch between the old U-vector OLDU and the C direction U-vector UDIR to obtain the new "best" guess BESTU C and its trajectory BESTX. C b. Perform a linesearch between the old V-vector OLDV and the C direction V-vector VDIR to obtain the new "bset" guess BESTV C and its trajectory BESTY. C c. Return the function values for BESTU and BESTV, as well as the C "alternative" elements ALTV and ALTU which give the respective C max and min defining these function values. 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 [ For a full description of this data, see the documentation in C the routine 'DFGM'. ] C C Variables whose values may be changed within LINE: C U,V,UE,VE -- finite sets of decision vectors; C X,Y,XE,YE -- the corresponding trajectories; C ALPHA,BETA -- the function values associated with the new best C guesses BESTU and BESTV. C QUIT -- Flag whose value is .FALSE. when LINE is called. If it C is decided that the algorithm should be immediately C terminated (due to some criterion specified within C LINE by the user), then LINE 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 ALPHA(MAXUFE),BETA(MAXVFE) LOGICAL QUIT C C Pointers (indicating positions in the finite sets) C 1. UDIR,VDIR. UDIR indicates the position in the U array held by C the U-vector giving the new search "direction". Similarly C for VDIR. C 2. BESTU,BESTV. Indicate the positions in which to place the new C best vectors. C 3. ALTU,ALTV. Indicate the positions in which to place the C "alternative" elements corresponding to the new best vectors. C 4. OLDU,OLDV. OLDU indicates the position in the U array held by the C previous best U-vector. LINE is expected to perform a line- C search between this vector and the direction vector described C above. OLDV serves a similarly purpose. INTEGER OLDU,OLDV INTEGER BESTU,BESTV INTEGER ALTU,ALTV INTEGER UDIR,VDIR C C Data (values not to be changed in LINE) 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 LINE and other routines called during the C algorithm (NOTE: values placed in this array cannot be saved between C successive calls to LINE.) 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 Stephen E. Wright IN August 1988 C --- --------- ---- --- ------- -- ================= -- =========== C AT IIASA, Laxenburg Austria. 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. Locally needed variable which must be saved: initialized in INIT DOUBLE PRECISION EPS COMMON /EPSLIN/ EPS C C Local Variables: DOUBLE PRECISION EPSNRM,RATIO,USTEP,VSTEP,FVAL,SUM INTEGER NOTIFY C C----------------------------------------------------------------------- c DO 181 JJ=1,3 IF (PROPT(JJ,1).EQ.3) THEN WRITE(PROPT(JJ,5),'(/"LINE SEARCH IN PRIMAL VARIABLES:"/)') ENDIF 181 CONTINUE C C Linesearch between the Old BESTU and UDIR to get the New BESTU, C the Best Primal Value ALPHA = f(BESTU), and the C corresponding Dual Element ALTV. C C CALL LNSRCH('min',U(1,1,UDIR),UE(1,UDIR),X(1,1,UDIR),XE(1,UDIR), : U(1,1,OLDU),UE(1,OLDU),X(1,1,OLDU),XE(1,OLDU), : PVEC,PMAT,QVEC,QMAT,PEVEC,PEMAT,QEVEC,QEMAT,DMAT, : CMAT,CVEC,CEMAT,CEVEC,VLOW,VUPR,VELOW,VEUPR, : USTEP,V(1,1,ALTV),VE(1,ALTV),FVAL,NOTIFY, : UDIM,UEDIM,VDIM,VEDIM,XYDIM,TIME,WORK,LENWRK) IF ( NOTIFY .GT. 0 ) THEN WRITE(*,*) 'MEMORY ALLOCATION INSUFFICIENT FOR: lnsrch' STOP ENDIF C ALPHA(BESTU) = FVAL C Set BESTU = USTEP*UDIR + (1-USTEP)*BESTU C (So that BESTU represents the "Best So Far"). C C DO 30 J=1,TIME DO 10 I=1,UDIM SUM = (1.0D0-USTEP)*U(I,J,OLDU) 10 U(I,J,BESTU) = SUM + USTEP*U(I,J,UDIR) DO 20 I=1,XYDIM SUM = (1.0D0-USTEP)*X(I,J,OLDU) 20 X(I,J,BESTU) = SUM + USTEP*X(I,J,UDIR) 30 CONTINUE DO 40 I=1,UEDIM SUM = (1.0D0-USTEP)*UE(I,OLDU) 40 UE(I,BESTU) = SUM + USTEP*UE(I,UDIR) DO 50 I=1,XYDIM SUM = (1.0D0-USTEP)*XE(I,OLDU) 50 XE(I,BESTU) = SUM + USTEP*XE(I,UDIR) C C DO 182 JJ=1,3 IF (PROPT(JJ,1).EQ.3) THEN WRITE(PROPT(JJ,5),'(/"LINE SEARCH IN DUAL VARIABLES:"/)') ENDIF 182 CONTINUE C C Linesearch between BESTV and VDIR to get the New BESTV, C the Best Dual Value BETA = g(BESTV), and the C corresponding Primal Element ALTV. C C CALL LNSRCH('max',V(1,1,VDIR),VE(1,VDIR),Y(1,1,VDIR),YE(1,VDIR), : V(1,1,OLDV),VE(1,OLDV),Y(1,1,OLDV),YE(1,OLDV), : QVEC,QMAT,PVEC,PMAT,QEVEC,QEMAT,PEVEC,PEMAT,DMAT, : BMAT,BVEC,BEMAT,BEVEC,ULOW,UUPR,UELOW,UEUPR, : VSTEP,U(1,1,ALTU),UE(1,ALTU),FVAL,NOTIFY, : VDIM,VEDIM,UDIM,UEDIM,XYDIM,TIME,WORK,LENWRK) IF ( NOTIFY .GT. 0 ) THEN WRITE(*,*) 'MEMORY ALLOCATION INSUFFICIENT FOR: lnsrch' STOP ENDIF C BETA(BESTV) = FVAL C Set BESTV = VSTEP*VDIR + (1-VSTEP)*BESTV C (So that BESTV represents the "Best So Far"). C DO 80 J=1,TIME DO 60 I=1,VDIM SUM = (1.0D0-VSTEP)*V(I,J,OLDV) 60 V(I,J,BESTV) = SUM + VSTEP*V(I,J,VDIR) DO 70 I=1,XYDIM SUM = (1.0D0-VSTEP)*Y(I,J,OLDV) 70 Y(I,J,BESTV) = SUM + VSTEP*Y(I,J,VDIR) 80 CONTINUE DO 90 I=1,VEDIM SUM = (1.0D0-VSTEP)*VE(I,OLDV) 90 VE(I,BESTV) = SUM + VSTEP*VE(I,VDIR) DO 100 I=1,XYDIM SUM = (1.0D0-VSTEP)*YE(I,OLDV) 100 YE(I,BESTV) = SUM + VSTEP*YE(I,VDIR) C C C Compute the Trajectories X,Y for ALTU,ALTV . C CALL TRAJ(U(1,1,ALTU),UE(1,ALTU), : V(1,1,ALTV),VE(1,ALTV), : X(1,1,ALTU),XE(1,ALTU), : Y(1,1,ALTV),YE(1,ALTV), : AMAT,BMAT,CMAT,BVEC,CVEC, : BEMAT,CEMAT,BEVEC,CEVEC, : UDIM,UEDIM,XYDIM,VDIM,VEDIM,TIME) C C - - - - - - - - - - - - - - - - - - - - - - - - - C C Compute the Duality Gap for BESTU,BESTV. C RATIO = (ALPHA(BESTU) - BETA(BESTV))/EPS EPS = ALPHA(BESTU) - BETA(BESTV) IF ( DABS(ALPHA(BESTU)) .GT. 1.0D0 ) THEN EPSNRM = EPS/DABS(ALPHA(BESTU)) ELSE EPSNRM = EPS ENDIF C 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 RETURN C End of subroutine 'LINE' END :::::::::::::: plnprep.f :::::::::::::: subroutine 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, # zlft,zrt,zmid,brk1,brk2,slope,btil,dtil, # zelft,zert,zemid,ebrk1,ebrk2,eslope,betil,detil, # wdim,zdim,wedim,zedim,xydim,ntime,type) 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 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) double precision btil(zdim,ntime),dtil(zdim,ntime),ct,atv,atm double precision betil(zedim),detil(zedim) c c dummy variables and counters: double precision sum1,sum2 integer i,j c c*********************************************************************** c This subroutine takes the given input data and creates the data for a c reformulated saddle-point problem. If TYPE='min' then the output data is c ct = sum [a *w2 + (1/2)w2 *A w2 - c*x2 ] c t=1,N t t t t t-1 c c + a *w2 + (1/2)w2 *A w2 - c *x2 , c e e e e e e N c c atv = sum [ (a + A w2 )*(w1 - w2 ) - c*(x1 - x2 ) ] c t=1,N t t t t t-1 t-1 c c + (a + A w2 )*(w1 - w2 ) - c *(x1 - x2 ) , c e e e e e e N N c c atm = sum [ (w1 - w2 )*A (w1 - w2 ) ] + (w1 - w2 )*A (w1 - w2 ) , c t=1,N t t t t e e e e e c c dtil = D(w2 - w1 ) + C(x2 - x1 ) for t=1,...,N , c t t t t-1 t-1 c c detil = C (x2 - x1 ) , c e N N c c btil = b - Dw2 - Cx2 for t=1,...N , c t t t t-1 c c betil = b - C x2 . c e e N c c If TYPE='max' then the output data is c c ct = sum [ -a *w2 + (1/2)w2 *A w2 + c*x2 ] c t=1,N t t t t t-1 c c - a *w2 + (1/2)w2 *A w2 + c *x2 , c e e e e e e N c c atv = sum [ (-a + A w2 )*(w1 - w2 ) + c*(x1 - x2 ) ] c t=1,N t t t t t-1 t-1 c c + (-a + A w2 )*(w1 - w2 ) + c *(x1 - x2 ) , c e e e e e e N N c c atm = sum [ (w1 - w2 )*A (w1 - w2 ) ] + (w1 - w2 )*A (w1 - w2 ) , c t=1,N t t t t e e e e e c c T T c dtil = -D (w2 - w1 ) - C (x2 - x1 ) for t=1,...,N , c t t t t-1 t-1 c c T c detil = -C (x2 - x1 ) , c e N N c c T T c btil = -b + D w2 + C x2 for t=1,...N , c t t t t-1 c c T c betil = -b + C x2 . c e e N c c c The resulting problem has the form c c minimax ct + atv*lambda + (1/2)lambda*atm*lambda c c + sum [ btil *z - (1/2)z *bm z + (z *dtil )*lambda ] c t=1,N t t t t t t c c + betil*z - (1/2)z *bem z + (z *detil)*lambda c e e e e c c subject to lambda in [0,1] and zlo <= z <= zhi , zelo<= z <=zehi . c e c The primal problem is then c minimize ct + atv*lambda + (1/2)lambda*atm*lambda c + sum [ sum rho (lambda;i) ] + sum rho (lambda;i) c t=1,N i=1,zdim t i=1,zedim e c subject to 0<= lambda <=1 , c where c rho (lambda;i) = c t c max [btil (i)+dtil (i)*lambda -(1/2)z (i)*bm (i)]*z (i) c zlo(i)<=z (i)<=zhi(i) t t t t c t c c = [btil (i) + dtil (i)*lambda - (1/2)zz (lambda;i)*bm (i)]*zz (lambda;i) c t t t t c c / c < zlo(i) if btil (i)+lambda*dtil (i) <= bm (i)*zlo(i) c > t t c and zz (lambda;i) = < zhi(i) if btil (i)+lambda*dtil (i) >= bm (i)*zhi(i) c t > t t c < [btil (i)+lambda*dtil (i)]/bm otherwise c \ t t c c and rho and zz are similarly defined. c e e c c c Note that zz (lambda;i) is a piecewise linear function in lambda which is c t c monotone increasing and can be represented by (at most) six pieces of data, c / c i.e. < zlft (i) if lambda <= brk1 (i) c > t t c zz(lambda;i) = < zrt (i) if lambda >= brk2 (i) c > t t c < zmid (i) + lambda*slope (i) otherwise. c \ t t c c The values of brk1,brk2,zlft,zrt,zmid and slope are also computed here. c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Compute CT, ATV, ATM. ct = 0.0d0 atv = 0.0d0 atm = 0.0d0 if ( type .ne. 'max' ) then c !! First for the term with t=1. do 01 j=1,wdim ct = ct + (av(j,1)+w2(j,1)*am(j)/2.0d0)*w2(j,1) atv = atv + (av(j,1)+am(j)*w2(j,1)) * (w1(j,1)-w2(j,1)) atm = atm + (w1(j,1)-w2(j,1))*am(j)*(w1(j,1)-w2(j,1)) 01 continue do 02 j=1,xydim ct = ct - cv(j)*xye2(j) atv = atv - cv(j) * (xye1(j)-xye2(j)) 02 continue c !! Now do the terms for t=2,...,N. do 06 k=2,ntime do 04 j=1,wdim ct = ct + (av(j,k)+w2(j,k)*am(j)/2.0d0)*w2(j,k) atv = atv + (av(j,k)+am(j)*w2(j,k)) * (w1(j,k)-w2(j,k)) atm = atm + (w1(j,k)-w2(j,k))*am(j)*(w1(j,k)-w2(j,k)) 04 continue do 05 j=1,xydim ct = ct - cv(j)*xy2(j,k-1) atv = atv - cv(j) * (xy1(j,k-1)-xy2(j,k-1)) 05 continue 06 continue c !! Finally, do the endpoint term. do 07 j=1,wedim ct = ct + (aev(j)+we2(j)*aem(j)/2.0d0)*we2(j) atv = atv + (aev(j)+aem(j)*we2(j)) * (we1(j)-we2(j)) atm = atm + (we1(j)-we2(j))*aem(j)*(we1(j)-we2(j)) 07 continue do 10 j=1,xydim ct = ct - cev(j)*xy2(j,ntime) atv = atv - cev(j) * (xy1(j,ntime)-xy2(j,ntime)) 10 continue else c !! First do the endpoint term. do 11 j=1,wedim ct = ct + (-aev(j)+we2(j)*aem(j)/2.0d0)*we2(j) atv = atv + (-aev(j)+aem(j)*we2(j)) * (we1(j)-we2(j)) atm = atm + (we1(j)-we2(j))*aem(j)*(we1(j)-we2(j)) 11 continue do 12 j=1,xydim ct = ct + cev(j)*xy2(j,1) atv = atv + cev(j) * (xy1(j,1)-xy2(j,1)) 12 continue c !! Now do the terms for t=1,...,N-1. do 16 k=1,ntime-1 do 14 j=1,wdim ct = ct + (-av(j,k)+w2(j,k)*am(j)/2.0d0)*w2(j,k) atv =atv + (-av(j,k)+am(j)*w2(j,k)) * (w1(j,k)-w2(j,k)) atm = atm + (w1(j,k)-w2(j,k))*am(j)*(w1(j,k)-w2(j,k)) 14 continue do 15 j=1,xydim ct = ct + cv(j)*xy2(j,k+1) atv = atv + cv(j) * (xy1(j,k+1)-xy2(j,k+1)) 15 continue 16 continue c !! Finally, do the term with t=N. k = ntime do 18 j=1,wdim ct = ct + (-av(j,k)+w2(j,k)*am(j)/2.0d0)*w2(j,k) atv =atv + (-av(j,k)+am(j)*w2(j,k)) * (w1(j,k)-w2(j,k)) atm = atm + (w1(j,k)-w2(j,k))*am(j)*(w1(j,k)-w2(j,k)) 18 continue do 19 j=1,xydim ct = ct + cv(j)*xye2(j) atv = atv + cv(j) * (xye1(j)-xye2(j)) 19 continue endif c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c if ( type .ne. 'max' ) then c First we compute the values of DTIL and BTIL for t=1. c Set sum1 = D (w2 -w1 ) + C (x2 - x1 ) c 1 1 0 0 c c Set sum2 = b - D w2 - C x2 c 1 1 0 do 30 i=1,zdim sum1 = 0.0d0 sum2 = bv(i,1) do 20 j=1,wdim sum1 = sum1 + dm(i+(j-1)*zdim)*(w2(j,1)-w1(j,1)) sum2 = sum2 - dm(i+(j-1)*zdim)*w2(j,1) 20 continue do 22 j=1,xydim sum1 = sum1 + cm(i+(j-1)*zdim)*(xye2(j)-xye1(j)) sum2 = sum2 - cm(i+(j-1)*zdim)*xye2(j) 22 continue dtil(i,1) = sum1 btil(i,1) = sum2 c Now compute BRK1,BRK2,ZLFT,ZRT,ZMID, and SLOPE call lndata(sum1,sum2,bm(i),zhi(i),zlo(i), # brk1(i,1),brk2(i,1),slope(i,1), # zlft(i,1),zrt(i,1),zmid(i,1)) c 30 continue c c Next we compute the values of DTIL and BTIL for t=1,...,N. c Set sum1 = D (w2 -w1 ) + C (x2 - x1 ) c t t t-1 t-1 c c Set sum2 = b - D w2 - C x2 c t t t-1 do 41 k=2,ntime do 40 i=1,zdim sum1 = 0.0d0 sum2 = bv(i,k) do 31 j=1,wdim sum1 = sum1 + dm(i+(j-1)*zdim)*(w2(j,k)-w1(j,k)) sum2 = sum2 - dm(i+(j-1)*zdim)*w2(j,k) 31 continue do 32 j=1,xydim sum1 = sum1 + cm(i+(j-1)*zdim)*(xy2(j,k-1)-xy1(j,k-1)) sum2 = sum2 - cm(i+(j-1)*zdim)*xy2(j,k-1) 32 continue dtil(i,k) = sum1 btil(i,k) = sum2 c Now compute BRK1,BRK2,ZLFT,ZRT,ZMID, and SLOPE call lndata(sum1,sum2,bm(i),zhi(i),zlo(i), # brk1(i,k),brk2(i,k),slope(i,k), # zlft(i,k),zrt(i,k),zmid(i,k)) c 40 continue 41 continue c c Finally we compute the values of DETIL and BETIL . c Set sum1 = C (x2 - x1 ) c e e e c c Set sum2 = b - C x2 c e e e k = ntime do 49 i=1,zedim sum1 = 0.0d0 sum2 = bev(i) do 45 j=1,xydim sum1 = sum1 + cem(i+(j-1)*zedim)*(xy2(j,k)-xy1(j,k)) sum2 = sum2 - cem(i+(j-1)*zedim)*xy2(j,k) 45 continue detil(i) = sum1 betil(i) = sum2 c Now compute EBRK1,EBRK2,ZELFT,ZERT,ZEMID, and ESLOPE call lndata(sum1,sum2,bem(i),zehi(i),zelo(i), # ebrk1(i),ebrk2(i),eslope(i), # zelft(i),zert(i),zemid(i)) c 49 continue else c The problem is a maximization. c T T c Set sum1 = - C (x2 -x1 ) , sum2 = -b + C x2 c e 1 1 e e 1 c do 130 i=1,zedim sum1 = 0.0d0 sum2 = -bev(i) do 122 j=1,xydim sum1 = sum1 - cem(j+(i-1)*xydim)*(xy2(j,1)-xy1(j,1)) sum2 = sum2 + cem(j+(i-1)*xydim)*xy2(j,1) 122 continue detil(i) = sum1 betil(i) = sum2 c Now compute EBRK1,EBRK2,ZELFT,ZERT,ZEMID, and ESLOPE call lndata(sum1,sum2,bem(i),zehi(i),zelo(i), # ebrk1(i),ebrk2(i),eslope(i), # zelft(i),zert(i),zemid(i)) c 130 continue c c Now we compute DTIL, BTIL for t=1,N-1. c c T T c Set sum1 = -D (w2 -w1 ) - C (x2 - x1 ) c t t t+1 t+1 c c T T c Set sum2 = -b + D w2 + C x2 c t t t+1 do 141 k=1,ntime-1 do 140 i=1,zdim sum1 = 0.0d0 sum2 = -bv(i,k) do 131 j=1,wdim sum1 = sum1 - dm(j+(i-1)*wdim)*(w2(j,k)-w1(j,k)) sum2 = sum2 + dm(j+(i-1)*wdim)*w2(j,k) 131 continue do 132 j=1,xydim sum1 = sum1 - cm(j+(i-1)*xydim)*(xy2(j,k+1)-xy1(j,k+1)) sum2 = sum2 + cm(j+(i-1)*xydim)*xy2(j,k+1) 132 continue dtil(i,k) = sum1 btil(i,k) = sum2 c Now compute BRK1,BRK2,ZLFT,ZRT,ZMID, and SLOPE call lndata(sum1,sum2,bm(i),zhi(i),zlo(i), # brk1(i,k),brk2(i,k),slope(i,k), # zlft(i,k),zrt(i,k),zmid(i,k)) c 140 continue 141 continue c c Finally we compute the values of DTIL and BTIL for t=N c c T T c Set sum1 = -D (w2 -w1 ) - C (x2 - x1 ) c N N N+1 N+1 c c c T T c Set sum2 = -b + D w2 + C x2 c N N N+1 k = ntime do 149 i=1,zdim sum1 = 0.0d0 sum2 = -bv(i,k) do 144 j=1,wdim sum1 = sum1 - dm(j+(i-1)*wdim)*(w2(j,k)-w1(j,k)) sum2 = sum2 + dm(j+(i-1)*wdim)*w2(j,k) 144 continue do 147 j=1,xydim sum1 = sum1 - cm(j+(i-1)*xydim)*(xye2(j)-xye1(j)) sum2 = sum2 + cm(j+(i-1)*xydim)*xye2(j) 147 continue dtil(i,ntime) = sum1 btil(i,ntime) = sum2 c Now compute BRK1,BRK2,ZLFT,ZRT,ZMID, and SLOPE call lndata(sum1,sum2,bm(i),zhi(i),zlo(i), # brk1(i,k),brk2(i,k),slope(i,k), # zlft(i,k),zrt(i,k),zmid(i,k)) c 149 continue c c END OF 'max' vs 'min' endif c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c END OF SUBROUTINE 'lnprep' return end ############################################### End of 92007-1.for