c OPTPACK was developed to test an algorithm presented in the paper c "Analysis and implementation of a dual algorithm for constrained c optimization," Journal of Optimization Theory and Applications, 79 (1993), c pp. 427-462. I have tested the software using the test problems c described in that paper as well other other research type problems, however, c I make no guarantees concerning the code. c c There are comment statements at the start of the 3 main routines c (bmin for bound constraints, lmin for bound + linear inequality constraints, c and min for bound + linear inequality + nonlinear constraints). c I can provide pointers, but the user is pretty much on his own. c c Bill Hager c c ------------------------------------------------------------------- c / / / ____/ William W. Hager c / / / / 358 Little Hall c / / / ____/ Department of Mathematics c /__/ / / UNIVERSITY OF FLORIDA c _______/ __/ Gainesville FL 32611-8105 c c Phone : (352) 392-0281 x 244 E-mail : hager@math.ufl.edu c FAX : (352) 392-6254 c http://www.math.ufl.edu/~hager c ------------------------------------------------------------------- C ________________________________________________________ C | | C |SOLVE AN OPTIMIZATION PROBLEM WITH NONLINEAR CONSTRAINTS| C | (ALL REAL VARIABLE ARE DOUBLE PRECISION) | C | | C | INPUT: | C | | C | X --STARTING GUESS (LENGTH AT LEAST N) | C | | C | ST --INITIAL STEP SIZE (WHEN ST = 0., PROGRAM | C | SELECTS INITIAL STEP | C | | C | IN --DESCRIBES THE STARTING POINT | C | =0, INFEASIBLE STARTING POINT, ALL | C | CONSTRAINTS NONBINDING INITIALLY | C | =1, INFEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS ARE GIVEN BY VALUES OF IA | C | =2, INFEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS DETERMINED BY X | C | =3, INFEASIBLE STARTING POINT, USE K, IA, | C | AND U AS GIVEN | C | =4, FEASIBLE STARTING POINT, ALL | C | CONSTRAINTS NONBINDING INITIALLY | C | =5, FEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS ARE GIVEN BY VALUES OF IA | C | =6, FEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS DETERMINED BY X | C | =7, FEASIBLE STARTING POINT, USE K, IA, | C | AND U AS GIVEN | C | | C | TL --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | LM --MAXIMUM NUMBER OF ITERATIONS | C | | C | PNI --INPUT PENALTY PARAMETER (TRY PNI= 0 FIRST,| C | IF A POSITIVE PENALTY IS REQUIRED, PROGRAM| C | STOPS, CALL ERR(IO) TO GET ERROR MESSAGE) | C | | C | A --COEFFICIENT MATRIX FOR LINEAR CONSTRAINTS | C | (AT LEAST L+NL ROWS) | C | | C | LA --LEADING DIMENSION OF A ARRAY | C | | C | LG --LEADING DIMENSION OF DG ARRAY | C | | C | L --NUMBER OF LINEAR EQUALITY CONSTRAINTS | C | | C | NL --NUMBER OF NONLINEAR EQUALITY CONSTRAINTS | C | | C | N --NUMBER OF UNKNOWNS | C | | C | B --RIGHT SIDE VECTOR FOR LINEAR CONSTRAINTS | C | | C | BL --LOWER BOUNDS FOR COMPONENTS OF X | C | | C | BU --UPPER BOUNDS FOR COMPONENTS OF X | C | | C | CT --CUTOFF FOR DECIDING IF A CONSTRAINT IS | C | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE| C | OF NORM OF TRUE SOLUTION) | C | | C | VLF --NAME OF ROUTINE TO EVALUATE COST FUNTION | C | | C | VLG --NAME OF ROUTINE TO EVALUATE CONSTRAINTS | C | | C | GRF --NAME OF ROUTINE TO EVALUATE COST GRADIENT | C | | C | GRG --NAME OF ROUTINE TO EVALUATE CONSTRAINT | C | GRADIENT | C | | C | W --WORK ARRAY (LENGTH AT LEAST 5M+3N IF | C | PENALTY NOT INCREASED, LENGTH AT LEAST | C | .5NL(NL+5)+L+3N OTHERWISE | C | | C | | C | OUTPUT: | C | | C | X --SOLUTION | C | | C | LE --MULTIPLIERS FOR EQUALITY CONSTRAINTS | C | (LENGTH AT LEAST L+NL) | C | | C | LI --MULTIPLIERS FOR INEQUALITY CONSTRAINTS | C | (LENGTH AT LEAST N) | C | | C | U --CHOLESKY FACTORIZATION ASSOCIATED WITH | C | CONSTRAINT PROJECTION (LENGTH AT LEAST | C | .5M(M+1) WHERE M = L + NL) | C | | C | IA --= 0 IF CONSTRAINT IS NONBINDING | C | = 1 IF CONSTRAINT IS AT UPPER BOUND | C | =-1 IF CONSTRAINT IS AT LOWER BOUND | C | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) | C | | C | K --NUMBER OF FREE VARIABLES | C | | C | F --VALUE OF COST AT OPTIMAL POINT | C | | C | G --VALUE OF CONSTRAINTS AT OPTIMAL POINT | C | (LENGTH AT LEAST NL) | C | | C | DF --DERIVATIVE OF COST AT OPTIMAL POINT | C | (LENGTH AT LEAST N) | C | | C | DG --DERIVATIVE OF CONSTRAINTS AT OPTIMAL POINT | C | | C | E --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | IT --NUMBER OF ITERATIONS | C | | C | ST --FINAL STEP SIZE | C | | C | PN --FINAL VALUE OF PENALTY | C | | C | IO --ERROR FLAG (CALL ERR(IO) TO SEE IF AN | C | ERROR WAS DETECTED) | C |________________________________________________________| C subroutine min(x,le,li,u,ia,k,f,g,df,dg,e,it,st,pn,io, 1 in,tl,lm,pni,a,la,lg,l,nl,n,b,bl,bu,ct,vlf,vlg,grf,grg,w) integer ia(1),i,in,io,ip,is,it,j,k,l,la,lg,lm integer m,m1,m2,m3,m4,m5,m6,m7,n,nl,n1,n2,n3,n4 real*8 a(la,1),b(1),bl(1),bu(1),df(1),dg(lg,1),g(1),le(1),li(1) real*8 u(1),w(1),x(1),big,ct,e,f,kt,pn,pni,r,s,st,t,tl,vlf external grf,grg,vlf,vlg m = l + nl if ( m .le. n ) goto 10 io = 8 return 10 if ( m .gt. 0 ) goto 20 call bmin(x,ia,k,f,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vlf,grf,w) return 20 if ( nl .gt. 0 ) goto 25 call lmin(x,le,li,u,ia,k,f,df,e,it,st,io, & in,tl,lm,a,la,l,n,b,bl,bu,ct,vlf,grf,w) return 25 m1 = m + 1 m2 = m + m1 m3 = m + m2 m4 = m + m3 m5 = n + m4 m6 = m + m5 + n m7 = m6 - 1 n1 = .5*(nl*nl+nl) + 1 n2 = n1 + nl n3 = n2 + l n4 = n3 + nl pn = pni c c largest floating point number c big = 1.d50 if ( in .lt. 4 ) goto 110 if ( in .eq. 7 ) goto 130 k = n - l if ( in .eq. 6 ) goto 60 if ( in .eq. 5 ) goto 40 do 30 i = 1,n 30 ia(i) = 0 goto 90 40 do 50 i = 1,n if ( ia(i) .eq. 0 ) goto 50 k = k - 1 if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 50 continue goto 90 60 do 80 i = 1,n ia(i) = 0 if ( x(i) .le. bl(i) ) goto 70 if ( x(i) .lt. bu(i) ) goto 80 k = k - 1 x(i) = bu(i) ia(i) = 1 goto 80 70 k = k - 1 x(i) = bl(i) ia(i) = -1 80 continue 90 if ( k .ge. 0 ) goto 100 io = 1 return 100 call mat(u,ia,a,la,l,n,0) call fac(u,0,l,io) if ( io .eq. 0 ) goto 130 io = 1 return 110 call feas(li,u,ia,k,io,in,x,a,la,l,n,b,bl,bu,w) do 120 i = 1,n 120 x(i) = li(i) if ( io .gt. 0 ) return 130 f = vlf(x) call vlg(g,x) call grf(df,x) call grg(dg,x) it = 0 io = 0 call erg(t,g,nl) 140 is = 0 ip = 0 if ( k .lt. nl ) ip = 1 if ( k .ge. nl ) k = k - nl call xpn(u,l,m) call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m) call err(io) if ( io .gt. 0 ) goto 420 call mul(le,li,s,u,ia,a,la,m,n,df,l,nl,w(m2)) e = s + t if ( e .le. tl ) return kt = .95*e goto 190 150 call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m) call err(io) if ( io .gt. 0 ) goto 420 call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2)) e = s + t if ( e .gt. tl ) goto 180 160 do 170 i = 1,m 170 le(i) = w(i+m) return 180 if ( t .le. s ) goto 320 190 is = is + 1 if ( is .lt. 40 ) goto 200 io = 6 return 200 if ( l .eq. 0 ) goto 230 do 210 i = 1,l 210 w(i) = b(i) do 220 j = 1,n s = x(j) do 220 i = 1,l 220 w(i) = w(i) - s*a(i,j) 230 do 240 i = 1,n 240 w(i+m7) = ia(i) if ( ip .eq. 1 ) goto 250 call nf1(li,u,ia,k,io,x,a,la,m,n,w,bl,bu,w(m1),w(m2),w(m2+n)) call err(io) if ( io .gt. 0 ) goto 420 if ( io .eq. 0 ) goto 290 io = 0 goto 300 250 do 260 i = 1,m w(i+m) = 0. 260 if ( i .gt. l ) w(i+m) = 1. call nf2(li,u,ia,k,io,x,a,la,m,n,w,w(m1),bl,bu,w(m2),w(m3),w(m5)) if ( io .eq. 0 ) goto 280 do 270 i = 1,n j = w(i+m7) call mod(u,ia,k,io,a,la,m,i,j,w) if ( io .gt. 0 ) return 270 continue goto 420 280 ip = 0 290 call lsq(li,u,ia,k,io,x,a,la,m,n,bl,bu,w,w(m1),w(m2)) call err(io) if ( io .gt. 0 ) goto 420 300 do 310 i = 1,n 310 w(i) = x(i) call arm(x,li,w,u,ia,k,io,t,w(m6),a,la,m,n,vlg,g,nl,tl) f = vlf(x) call grf(df,x) call grg(dg,x) if ( io .eq. 0 ) goto 150 return 320 r = .95*s is = -1 do 330 i = 1,m 330 le(i) = w(i+m) 340 is = is + 1 if ( is .le. 3 ) goto 350 if ( e .le. kt ) goto 390 goto 420 350 j = k call cgn(x,e,i,st,io,.1*tl,j,lm-it,n,is,vlf,vlg,grf,grg, 1 nl,pn,a,la,m,bl,bu,ct,u,ia,k,le(l+1),f,g,df,dg,lg, 2 w,w(m1),w(m2),w(m3),w(m4)) call err(io) it = it + i if ( io .gt. 0 ) return call erg(t,g,nl) call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m) call err(io) if ( io .gt. 0 ) return call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2)) e = s + t if ( e .le. tl ) goto 160 if ( it .lt. lm ) goto 360 io = 9 goto 160 360 if ( s .gt. r ) goto 380 r = .95*s is = 0 do 370 i = 1,m 370 le(i) = w(i+m) if ( 4.*t .lt. s ) goto 340 if ( e .gt. kt ) goto 410 kt = .95*e goto 190 380 if ( t .gt. .5*kt ) goto 420 if ( 4.*t .lt. s ) goto 340 if ( e .gt. kt ) goto 410 390 kt = .95*e is = 0 do 400 i = 1,m 400 le(i) = w(i+m) goto 190 410 if ( t .le. .25*kt ) goto 340 420 call shk(u,l,m) is = 0 if ( ip .eq. 0 ) k = k + nl if ( it .gt. 0 ) pn = 5.*pn if ( pn .gt. 0. ) goto 430 io = 7 return 430 j = k call cgp(x,e,i,st,io,.1*tl,j,lm-it,n,vlf,vlg,grf,grg, 1 l,nl,pn,a,la,bl,bu,ct,u,ia,k,le(l+1),f,g,df,dg,lg, 2 w,w(n1),w(n2),w(n3),w(n4)) call err(io) it = it + i if ( io .gt. 0 ) return if ( it .lt. lm ) goto 440 io = 9 return 440 call erg(t,g,nl) call xpn(u,l,m) ip = 0 if ( k .lt. nl ) ip = 1 if ( k .ge. nl ) k = k - nl call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m) call err(io) if ( io .gt. 0 ) return call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2)) r = e e = s + t if ( e .le. tl ) goto 160 if ( e .le. .6*kt ) goto 450 if ( r .le. t ) goto 460 call shk(u,l,m) if ( ip .eq. 0 ) k = k + nl goto 430 450 pn = pn*.4 460 do 470 i = 1,m 470 le(i) = w(i+m) kt = e goto 190 end C ________________________________________________________ C | | C | SOLVE AN OPTIMIZATION PROBLEM WITH UPPER AND LOWER | C | BOUND CONSTRAINTS | C | (ALL REAL VARIABLE ARE DOUBLE PRECISION) | C | | C | INPUT: | C | | C | X --STARTING GUESS (LENGTH AT LEAST N) | C | | C | ST --INITIAL STEP SIZE (WHEN ST = 0., PROGRAM | C | SELECTS INITIAL STEP | C | | C | IN --DESCRIBES THE STARTING POINT | C | =0, ALL CONSTRAINTS NONBINDING INITIALLY | C | =1, BINDING CONSTRAINTS GIVEN BY IA | C | =2, BINDING CONSTRAINTS DETERMINED BY X | C | =3, USE K, IA, AND W AS GIVEN | C | | C | TL --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | LM --MAXIMUM NUMBER OF ITERATIONS | C | | C | N --NUMBER OF UNKNOWNS | C | | C | BL --LOWER BOUNDS FOR COMPONENTS OF X | C | | C | BU --UPPER BOUNDS FOR COMPONENTS OF X | C | | C | CT --CUTOFF FOR DECIDING IF A CONSTRAINT IS | C | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE| C | OF NORM OF TRUE SOLUTION) | C | | C | VL --NAME OF ROUTINE TO EVALUATE COST FUNTION | C | | C | GR --NAME OF ROUTINE TO EVALUATE COST GRADIENT | C | | C | W --WORK ARRAY (LENGTH AT LEAST 3N) | C | | C | OUTPUT: | C | | C | X --SOLUTION | C | | C | IA --= 0 IF CONSTRAINT IS NONBINDING | C | = 1 IF CONSTRAINT IS AT UPPER BOUND | C | =-1 IF CONSTRAINT IS AT LOWER BOUND | C | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) | C | | C | K --NUMBER OF FREE VARIABLES | C | | C | VF --VALUE OF COST AT OPTIMAL POINT | C | | C | DF --DERIVATIVE OF COST AT OPTIMAL POINT | C | (LENGTH AT LEAST N) | C | | C | E --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | IT --NUMBER OF ITERATIONS | C | | C | ST --FINAL STEP SIZE | C | | C | IO --ERROR FLAG (CALL ERR(IO) TO SEE IF AN | C | ERROR WAS DETECTED) | C |________________________________________________________| C subroutine bmin 1 (x,ia,k,vf,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vl,gr,w) real*8 a(1),b(1),bl(1),bu(1),c(1),df(1),x(1),u(1),w(1) real*8 ct,e,st,t,tl,vf,vl integer ia(1),i,in,io,it,j,k,lm,n external gr,vl it = 0 io = 0 if ( in .eq. 3 ) goto 80 k = n do 10 i = 1,n if ( bl(i) .lt. bu(i) ) goto 10 t = bl(i) bl(i) = bu(i) bu(i) = t 10 continue if ( in .eq. 2 ) goto 50 if ( in .eq. 1 ) goto 30 do 20 i = 1,n 20 ia(i) = 0 goto 80 30 do 40 i = 1,n if ( ia(i) .eq. 0 ) goto 40 k = k - 1 if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 40 continue goto 80 50 do 70 i = 1,n ia(i) = 0 if ( x(i) .lt. bu(i) ) goto 60 x(i) = bu(i) ia(i) = 1 k = k - 1 goto 70 60 if ( x(i) .gt. bl(i) ) goto 70 x(i) = bl(i) ia(i) = -1 k = k - 1 70 continue 80 call gr(df,x) vf = vl(x) 90 j = k call cgl(x,e,i,st,io,tl,j,lm-it,n,vl,gr, 1 a,1,0,bl,bu,ct,u,ia,k,vf,df,b,c,w) it = it + i if ( io .gt. 0 ) return if ( e .le. tl ) return if ( it .lt. lm ) goto 90 io = 9 return end C ________________________________________________________ C | | C | SOLVE AN OPTIMIZATION PROBLEM WITH LINEAR CONSTRAINTS | C | (ALL REAL VARIABLE ARE DOUBLE PRECISION) | C | | C | INPUT: | C | | C | X --STARTING GUESS (LENGTH AT LEAST N) | C | | C | ST --INITIAL STEP SIZE (WHEN ST = 0., PROGRAM | C | SELECTS INITIAL STEP | C | | C | IN --DESCRIBES THE STARTING POINT | C | =0, INFEASIBLE STARTING POINT, ALL | C | CONSTRAINTS NONBINDING INITIALLY | C | =1, INFEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS ARE GIVEN BY VALUES OF IA | C | =2, INFEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS DETERMINED BY X | C | =3, INFEASIBLE STARTING POINT, USE K, IA, | C | AND W AS GIVEN | C | =4, FEASIBLE STARTING POINT, ALL | C | CONSTRAINTS NONBINDING INITIALLY | C | =5, FEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS ARE GIVEN BY VALUES OF IA | C | =6, FEASIBLE STARTING POINT, BINDING | C | CONSTRAINTS DETERMINED BY X | C | =7, FEASIBLE STARTING POINT, USE K, IA, | C | AND U AS GIVEN | C | | C | TL --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | LM --MAXIMUM NUMBER OF ITERATIONS | C | | C | A --COEFFICIENT MATRIX FOR LINEAR CONSTRAINTS | C | (AT LEAST M ROWS) | C | | C | LA --LEADING DIMENSION OF A ARRAY | C | | C | M --NUMBER OF EQUALITY CONSTRAINTS | C | | C | N --NUMBER OF UNKNOWNS | C | | C | B --RIGHT SIDE VECTOR FOR LINEAR CONSTRAINTS | C | (LENGTH AT LEAST M) | C | | C | BL --LOWER BOUNDS FOR COMPONENTS OF X | C | | C | BU --UPPER BOUNDS FOR COMPONENTS OF X | C | | C | CT --CUTOFF FOR DECIDING IF A CONSTRAINT IS | C | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE| C | OF NORM OF TRUE SOLUTION) | C | | C | VL --NAME OF ROUTINE TO EVALUATE COST FUNTION | C | | C | GR --NAME OF ROUTINE TO EVALUATE COST GRADIENT | C | | C | W --WORK ARRAY(LENGTH AT LEAST MAX(4M+2N,M+3N)| C | | C | OUTPUT: | C | | C | X --SOLUTION | C | | C | LE --MULTIPLIERS FOR EQUALITY CONSTRAINTS | C | | C | LI --MULTIPLIERS FOR INEQUALITY CONSTRAINTS | C | | C | U --CHOLESKY FACTORIZATION ASSOCIATED WITH | C | CONSTRAINT PROJECTION | C | (LENGTH AT LEAST .5M(M+1)) | C | | C | IA --= 0 IF CONSTRAINT IS NONBINDING | C | = 1 IF CONSTRAINT IS AT UPPER BOUND | C | =-1 IF CONSTRAINT IS AT LOWER BOUND | C | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) | C | | C | K --NUMBER OF FREE VARIABLES | C | | C | VF --VALUE OF COST AT OPTIMAL POINT | C | | C | DF --DERIVATIVE OF COST AT OPTIMAL POINT | C | (LENGTH AT LEAST N) | C | | C | E --ERROR TOLERANCE (INFINITY NORM OF KUHN- | C | TUCKER ERROR + CONSTRAINT VIOLATION) | C | | C | IT --NUMBER OF ITERATIONS | C | | C | ST --FINAL STEP SIZE | C | | C | IO --ERROR FLAG (CALL ERR(IO) TO SEE IF AN | C | ERROR WAS DETECTED) | C |________________________________________________________| C subroutine lmin(x,le,li,u,ia,k,vf,df,e,it,st,io, 1 in,tl,lm,a,la,m,n,b,bl,bu,ct,vl,gr,w) real*8 a(1),b(1),bl(1),bu(1),df(1),le(1),li(1),u(1),w(1),x(1),z(1) real*8 ct,e,tl,st,vf,vl integer ia(1),i,in,io,it,k,la,lm,m,m1,n external gr,vl m1 = m + 1 if ( m .le. n ) goto 10 io = 8 return 10 if ( m .gt. 0 ) goto 20 call bmin(x,ia,k,vf,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vl,gr,w) io = 0 return 20 if ( in .lt. 4 ) goto 100 if ( in .eq. 7 ) goto 120 k = n - m if ( in .eq. 6 ) goto 60 if ( in .eq. 5 ) goto 40 do 30 i = 1,n 30 ia(i) = 0 goto 90 40 do 50 i = 1,n if ( ia(i) .eq. 0 ) goto 50 k = k - 1 if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 50 continue if ( k .ge. 0 ) goto 90 io = 8 return 60 do 80 i = 1,n ia(i) = 0 if ( x(i) .le. bl(i) ) goto 70 if ( x(i) .lt. bu(i) ) goto 80 k = k - 1 x(i) = bu(i) ia(i) = 1 goto 80 70 k = k - 1 x(i) = bl(i) ia(i) = -1 80 continue if ( k .ge. 0 ) goto 90 io = 8 return 90 call mat(u,ia,a,la,m,n,0) call fac(u,0,m,io) if ( io .gt. 0 ) return goto 120 100 do 110 i = 1,n 110 li(i) = x(i) call feas(x,u,ia,k,io,in,li,a,la,m,n,b,bl,bu,w) if ( io .gt. 0 ) return 120 call gr(df,x) vf = vl(x) it = 0 io = 0 130 j = k c c Projected cg step c call cgl(x,e,i,st,io,tl,j,lm-it,n,vl,gr, 1 a,la,m,bl,bu,ct,u,ia,k,vf,df,w,le,w(m1)) it = it + i if ( io .eq. 1 ) return if ( io .gt. 1 ) goto 140 if ( it .ge. lm ) goto 140 if ( e .gt. tl ) goto 130 140 call mul(le,li,e,u,ia,a,la,m,n,df,m,0,z) return end subroutine arm(x,y,z,u,iy,k,io,t,iz,a,la,m,n,vlg,g,nl,tl) real*8 a(1),g(1),iz(1),u(1),x(1),y(1),z(1),r,s,t,tl integer iy(1),i,io,j,k,m,n,nl io = 0 call vlg(g,y) call erg(s,g,nl) if ( s .le. tl ) goto 10 if ( s .gt. .5*t ) goto 30 10 t = s do 20 i = 1,n 20 x(i) = y(i) return 30 do 40 i = 1,n if ( iy(i) .eq. 0 ) goto 40 if ( iy(i) .eq. iz(i) ) goto 40 call mod(u,iy,k,io,a,la,m,i,0,x) if ( io .eq. 1 ) return 40 continue r = .5 j = 0 50 j = j + 1 if ( j .gt. 20 ) goto 70 do 60 i = 1,n 60 x(i) = r*y(i) + (1-r)*z(i) call vlg(g,x) call erg(s,g,nl) r = .5*r if ( s .gt. t*(1-r) ) goto 50 t = s return 70 io = 6 t = s return end subroutine cgl(x,e,it,st,io,tl,lm1,lm2,n,vl,gr, 1 aa,la,m,bl,bu,ct,u,ia,k,vf,df,bb,cc,h) integer ia(1),i,ib,io,it,j,k,lm1,lm2,m,n,na,nb,nc,nd real*8 h(n,1),x(1),y(50),z(50) real*8 aa(1),bb(1),bl(1),bu(1),cc(1),df(1),u(1),vf real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0 real*8 da,db,e,f,f0,f1,fa,fb,fc,g,l3,p,q,r,s,st,tl,v,w real*8 fv,fd external gr,vl data a1/.1d0/,a2/.9d0/,a3/5.d0/,a4/.2d0/,a5/10.d0/,a6/.9d0/ data a7/.3d0/ c c set big = largest floating point number c big = 1.d50 a8 = a3 + .01d0 it = 0 io = 0 f = vf do 10 i = 1,n 10 h(i,3) = df(i) l3 = 1./dlog(a3) call pre(h(1,2),e,u,ia,k,i,lm1,io,h(1,3),aa,la,m,n,bb,cc) if ( e .le. tl ) goto 690 if ( io .gt. 0 ) goto 690 if ( lm1 .eq. 0 ) return a = st if ( a .gt. 0. ) goto 30 do 20 i = 1,n 20 if ( dabs(x(i)) .gt. a ) a = dabs(x(i)) a = .01*a/e if ( a .eq. 0. ) a = 1. 30 g = 0. do 40 i = 1,n h(i,1) = -h(i,2) 40 g = g + h(i,2)*h(i,3) if ( g .lt. 0. ) goto 650 d = -g 50 call cut(s,ib,x,h,ia,bl,bu,ct,n,big) na = 0 nb = 0 nc = 0 nd = 0 if ( s .le. 0. ) goto 760 if ( a .gt. s ) a = s fa = fv(a,x,h,n,vl) c0 = a f0 = fa ny = 2 y(1) = 0. z(1) = f y(2) = a z(2) = fa v = a1*d w = a2*d iq = 0 if ( fa .le. f ) goto 60 c = a b = 0. a = 0. fc = fa fb = f fa = f goto 70 60 c = 0. b = 0. fc = f fb = f iq = 1 70 q = (d+(f-f0)/c0)/c0 if ( q .lt. 0. ) goto 90 q = a p = fa 80 nd = nd + 1 if ( nd .gt. 25 ) goto 640 if ( s .eq. q ) goto 710 q = a3*q if ( q .gt. s ) q = s p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .lt. w*q ) goto 80 goto 270 90 q = .5*d/q if ( q .lt. s ) goto 110 if ( c0 .lt. s ) goto 100 f1 = f0 c1 = c0 q = s goto 140 100 q = s 110 if ( q .lt. .01*c0 ) q = .01*c0 p = fv(q,x,h,n,vl) if ( p .le. f0 ) goto 120 f1 = f0 c1 = c0 f0 = p c0 = q goto 130 120 f1 = p c1 = q 130 call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) 140 if ( a .eq. 0. ) goto 150 if ( fa-f .ge. v*a ) goto 170 if ( fa-f .lt. w*a ) goto 220 goto 290 150 q = c0 if ( c1 .lt. q ) q = c1 160 na = na + 1 if ( na .gt. 25 ) goto 660 q = a4*q p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .ge. v*q ) goto 160 goto 260 170 if ( c0 .gt. c1 ) goto 210 if ( f0-f .gt. v*c0 ) goto 190 if ( f0-f .ge. w*c0 ) goto 350 if ( c1 .le. a5*c0 ) goto 350 r = dlog(c1/c0) r = .999*dexp(-r/idint(r*l3+.999)) q = c1 180 q = q*r if ( q .lt. c0 ) goto 350 p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) na = na + 1 if ( p-f .gt. v*q ) goto 180 goto 350 190 q = c0 200 na = na + 1 if ( na .gt. 25 ) goto 660 q = a4*q p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .ge. v*q ) goto 200 goto 260 210 q = a goto 200 220 if ( c0 .lt. c1 ) goto 320 if ( f0-f .ge. v*c0 ) goto 240 if ( f0-f .ge. w*c0 ) goto 260 q = c0 p = f0 230 nd = nd + 1 if ( nd .gt. 25 ) goto 640 if ( s .eq. q ) goto 710 q = a3*q if ( q .gt. s ) q = s p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .lt. w*q ) goto 230 goto 260 240 if ( c0 .le. a5*c1 ) goto 260 r = dlog(c0/c1) r = 1.001*dexp(r/idint(r*l3+.999)) q = a 250 q = q*r if ( q .gt. c0 ) goto 260 nd = nd + 1 p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .lt. w*q ) goto 250 260 if ( iq .eq. 1 ) goto 350 270 if ( b .eq. 0. ) goto 290 if ( c .eq. 0. ) goto 280 v = c - a w = a - b r = 1./v f0 = 1./w p = fc - fa q = fb - fa e = p*r + q*f0 if ( dsign(e,c-b) .ne. e ) goto 350 if ( e .eq. 0. ) goto 350 q = (p*r)*w - (q*f0)*v q = a - .5*q/e goto 300 280 r = 1./a f0 = 1./b p = r*(fa-f) - d q = f0*(fb-f) - d e = a - b v = (r*p-f0*q)/e w = (a*q*f0-b*p*r)/e v = w*w-3.*v*d if ( v .lt. 0. ) v = 0. v = dsqrt(v) if ( w+v .eq. 0. ) goto 350 q = -d/(w+v) if ( q .le. 0. ) goto 350 goto 300 290 if ( iq .eq. 1 ) goto 350 q = (d+(f-fa)/a)/a if ( q .ge. 0. ) goto 350 q = .5*d/q 300 if ( q .gt. s ) q = s do 310 i = 1,ny 310 if ( y(i) .eq. q ) goto 350 p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) goto 350 320 continue if ( f0-f .gt. v*c0 ) goto 330 if ( f0-f .gt. w*c0 ) goto 350 330 q = a p = fa 340 nd = nd + 1 if ( nd .gt. 25 ) goto 640 if ( s .eq. q ) goto 710 q = a3*q if ( q .gt. s ) q = s p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .lt. w*q ) goto 340 goto 260 350 da = fd(a,x,h,n,gr) if ( da .gt. a6*g ) goto 440 if ( s .eq. a ) goto 780 if ( da .ge. 0. ) goto 590 r = a q = 0. do 360 i = 1,ny if ( y(i) .gt. a ) goto 400 if ( y(i) .le. q ) goto 360 if ( y(i) .eq. a ) goto 360 q = y(i) 360 continue if ( a .le. a8*q ) goto 590 q = a p = fa 370 nd = nd + 1 if ( nd .gt. 25 ) goto 640 if ( s .eq. q ) goto 710 q = a3*q if ( q .gt. s ) q = s p = fv(q,x,h,n,vl) f1 = fa call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p .lt. f1 ) goto 370 if ( a .gt. r ) goto 390 do 380 i = 1,n 380 h(i,2) = x(i) + a*h(i,1) goto 590 390 da = fd(a,x,h,n,gr) if ( da .gt. a6*g ) goto 440 goto 590 400 q = y(i) do 410 j = i,ny if ( y(j) .le. a ) goto 410 if ( y(j) .lt. q ) q = y(j) 410 continue if ( q .le. a5*a ) goto 590 f0 = dlog(q/a) f0 = 1.001*dexp(f0/idint(f0*l3+.999)) v = a 420 v = v*f0 if ( v .ge. q ) goto 350 p = fv(v,x,h,n,vl) f1 = fa call ins(v,p,a,b,c,fa,fb,fc,ny,y,z) if ( p .lt. f1 ) goto 420 if ( a .gt. r ) goto 350 do 430 i = 1,n 430 h(i,2) = x(i) + a*h(i,1) goto 590 440 b = 0. j = 1 i = j 450 i = i + 1 if ( i .gt. ny ) goto 460 if ( y(i) .ge. a ) goto 450 if ( y(i) .lt. b ) goto 450 b = y(i) j = i goto 450 460 fb = z(j) db = d if ( b .ne. 0. ) db = fd(b,x,h,n,gr) 470 w = 2.*dabs(b-a) call cub(c,a,b,fa,fb,da,db) nc = 1 goto 510 480 w = .5*w if ( w .lt. dabs(c0-c) ) goto 580 if ( c0 .lt. c ) goto 490 if ( d0 .ge. d ) goto 500 goto 580 490 if ( d0 .gt. d ) goto 580 500 call cub(c,c,c0,f,f0,d,d0) nc = nc + 1 if ( nc .gt. 30 ) goto 630 510 r = dmax1(a,b) s = dmin1(a,b) if ( c .gt. r ) goto 520 if ( c .gt. s ) goto 530 c = s + (s-c) s = .5*(a+b) if ( c .gt. s ) c = s goto 530 520 c = r - (c-r) s = .5*(a+b) if ( c .lt. s ) c = s 530 c0 = a f0 = fa d0 = da d = fd(c,x,h,n,gr) f = fv(c,x,h,n,vl) if ( f .lt. fa ) goto 540 b = c fb = f db = d goto 480 540 if ( c .lt. a ) goto 570 if ( d .lt. 0. ) goto 560 550 b = a fb = fa db = da 560 a = c fa = f da = d if ( d .gt. a6*g ) goto 480 goto 590 570 if ( d .lt. 0. ) goto 550 goto 560 580 c = .5*(a+b) nb = nb + 1 w = dabs(b-a) goto 530 590 do 600 i = 1,n 600 x(i) = h(i,2) i = lm1 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc) it = it + 1 f = fa d = da a = a7*a if ( io .gt. 0 ) goto 690 if ( e .le. tl ) goto 690 if ( it .ge. lm2 ) goto 690 if ( it .ge. lm1 ) goto 690 if ( i .lt. lm1 ) goto 30 r = 0. do 610 i = 1,n 610 r = r + h(i,2)*h(i,3) if ( r .lt. 0. ) goto 650 s = r/g g = r d = 0. do 620 i = 1,n h(i,1) = -h(i,2) + s*h(i,1) 620 d = d + h(i,1)*h(i,3) goto 50 630 if ( d .lt. g ) goto 590 io = 4 return 640 io = 3 return 650 io = 4 return 660 q = q*a3**25 nd = 0 670 nd = nd + 1 if ( nd .gt. 25 ) goto 680 q = a3*q if ( q .ge. s ) goto 680 p = fv(q,x,h,n,vl) call ins(q,p,a,b,c,fa,fb,fc,ny,y,z) if ( p-f .gt. v*q ) goto 670 goto 140 680 io = 5 return 690 st = a vf = f do 700 i = 1,n 700 df(i) = h(i,3) return 710 call mod(u,ia,k,io,aa,la,m,ib,1,bb) do 720 i = 1,n 720 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) a = q*a7 730 call gr(h(1,3),x) 740 f = p 750 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc) it = it + 1 if ( it .ge. lm2 ) goto 690 if ( it .ge. lm1 ) goto 690 if ( io .gt. 0 ) goto 690 if ( e .le. tl ) goto 690 goto 30 760 call mod(u,ia,k,io,aa,la,m,ib,1,bb) q = -s p = fv(q,x,h,n,vl) do 770 i = 1,n 770 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) goto 730 780 call mod(u,ia,k,io,aa,la,m,ib,1,bb) do 790 i = 1,n 790 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) p = fa goto 740 end subroutine cgn(x,e,it,st,io,tl,lm1,lm2,n,is,vlf,vlg,grf,grg,nl,pn, 1 aa,la,m,bl,bu,ct,u,ia,k,le,vf,vg,df,dg,lg,gg,bb,cc,dd,h) integer ia(1),i,ib,io,is,it,j,k,la,lg,lm1,lm2,m,n,nl,na,nb,nc,nd real*8 df(1),dg(lg,1),gg(1),h(n,1),le(1),vg(1),x(1),y(50),z(50) real*8 aa(1),bb(1),bl(1),bu(1),cc(1),dd(1),u(1) real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0,da,db real*8 e,f,f0,f1,fa,fb,fc,g,l3,p,p2,pn,q,r,s,st,tl,v,vf,vv,w real*8 pv,pd,pd2 external grf,grg,vlf,vlg data a1/.1d0/,a2/.9d0/,a3/5.d0/,a4/.2d0/,a5/10.d0/,a6/.9d0/ data a7/.3d0/ c c set big = largest floating point number c big = 1.d50 a8 = a3 + .01d0 p2 = .5*pn it = 0 io = 0 f = vf if ( is .eq. 0 ) goto 40 do 10 i = 1,nl dd(i) = vg(i) 10 f = f + vg(i)*le(i) do 30 j = 1,n s = df(j) do 20 i = 1,nl 20 s = s + le(i)*dg(i,j) h(j,3) = s 30 continue goto 80 40 do 50 i = 1,nl dd(i) = 0. bb(i) = le(i) + pn*vg(i) 50 f = f + vg(i)*le(i) + p2*vg(i)**2 do 70 j = 1,n s = df(j) do 60 i = 1,nl 60 s = s + bb(i)*dg(i,j) h(j,3) = s 70 continue 80 l3 = 1./dlog(a3) call pre(h(1,2),e,u,ia,k,i,lm1,io,h(1,3),aa,la,m,n,bb,cc) if ( e .le. tl ) goto 770 if ( io .gt. 0 ) goto 770 if ( lm1 .eq. 0 ) return a = st if ( a .gt. 0. ) goto 100 do 90 i = 1,n 90 if ( dabs(x(i)) .gt. a ) a = dabs(x(i)) a = .01*a/e if ( a .eq. 0. ) a = 1. 100 g = 0. do 110 i = 1,n h(i,1) = -h(i,2) 110 g = g + h(i,2)*h(i,3) if ( g .lt. 0. ) goto 730 d = -g 120 call cut(s,ib,x,h,ia,bl,bu,ct,n,big) na = 0 nb = 0 nc = 0 nd = 0 if ( s .le. 0. ) goto 860 if ( a .gt. s ) a = s p = pv(a,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) fa = p c0 = a f0 = fa ny = 2 y(1) = 0. z(1) = f y(2) = a z(2) = fa v = a1*d w = a2*d iq = 0 if ( fa .le. f ) goto 130 c = a b = 0. a = 0. fc = fa fb = f fa = f goto 150 130 c = 0. b = 0. fc = f fb = f iq = 1 vf = vv do 140 i = 1,nl 140 vg(i) = gg(i) 150 q = (d+(f-f0)/c0)/c0 if ( q .lt. 0. ) goto 170 q = a 160 nd = nd + 1 if ( nd .gt. 25 ) goto 720 if ( s .eq. q ) goto 780 q = a3*q if ( q .gt. s ) q = s p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 160 goto 350 170 q = .5*d/q if ( q .lt. s ) goto 190 if ( c0 .lt. s ) goto 180 f1 = f0 c1 = c0 q = s goto 220 180 q = s 190 if ( q .lt. .01*c0 ) q = .01*c0 p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) if ( p .le. f0 ) goto 200 f1 = f0 c1 = c0 f0 = p c0 = q goto 210 200 f1 = p c1 = q 210 call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) 220 if ( a .eq. 0. ) goto 230 if ( fa-f .ge. v*a ) goto 250 if ( fa-f .lt. w*a ) goto 300 goto 370 230 q = c0 if ( c1 .lt. q ) q = c1 240 na = na + 1 if ( na .gt. 25 ) goto 740 q = a4*q p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .ge. v*q ) goto 240 goto 340 250 if ( c0 .gt. c1 ) goto 290 if ( f0-f .gt. v*c0 ) goto 270 if ( f0-f .ge. w*c0 ) goto 430 if ( c1 .le. a5*c0 ) goto 430 r = dlog(c1/c0) r = .999*dexp(-r/idint(r*l3+.999)) q = c1 260 q = q*r if ( q .lt. c0 ) goto 430 p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) na = na + 1 if ( p-f .gt. v*q ) goto 260 goto 430 270 q = c0 280 na = na + 1 if ( na .gt. 25 ) goto 740 q = a4*q p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .ge. v*q ) goto 280 goto 340 290 q = a goto 280 300 if ( c0 .lt. c1 ) goto 400 if ( f0-f .ge. v*c0 ) goto 320 if ( f0-f .ge. w*c0 ) goto 340 q = c0 p = f0 310 nd = nd + 1 if ( nd .gt. 25 ) goto 720 if ( s .eq. q ) goto 780 q = a3*q if ( q .gt. s ) q = s p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 310 goto 340 320 if ( c0 .le. a5*c1 ) goto 340 r = dlog(c0/c1) r = 1.001*dexp(r/idint(r*l3+.999)) q = a 330 q = q*r if ( q .gt. c0 ) goto 340 nd = nd + 1 p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 330 340 if ( iq .eq. 1 ) goto 430 350 if ( b .eq. 0. ) goto 370 if ( c .eq. 0. ) goto 360 v = c - a w = a - b r = 1./v f0 = 1./w p = fc - fa q = fb - fa e = p*r + q*f0 if ( dsign(e,c-b) .ne. e ) goto 430 if ( e .eq. 0. ) goto 430 q = (p*r)*w - (q*f0)*v q = a - .5*q/e goto 380 360 r = 1./a f0 = 1./b p = r*(fa-f) - d q = f0*(fb-f) - d e = a - b v = (r*p-f0*q)/e w = (a*q*f0-b*p*r)/e v = w*w-3.*v*d if ( v .lt. 0. ) v = 0. v = dsqrt(v) if ( w+v .eq. 0. ) goto 430 q = -d/(w+v) if ( q .le. 0. ) goto 430 goto 380 370 if ( iq .eq. 1 ) goto 430 q = (d+(f-fa)/a)/a if ( q .ge. 0. ) goto 430 q = .5*d/q 380 if ( q .gt. s ) q = s do 390 i = 1,ny 390 if ( y(i) .eq. q ) goto 430 p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) goto 430 400 continue if ( f0-f .gt. v*c0 ) goto 410 if ( f0-f .gt. w*c0 ) goto 430 410 q = a p = fa 420 nd = nd + 1 if ( nd .gt. 25 ) goto 720 if ( s .eq. q ) goto 780 q = a3*q if ( q .gt. s ) q = s p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 420 goto 340 430 da = pd(a,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb) if ( da .gt. a6*g ) goto 520 if ( s .eq. a ) goto 880 if ( da .ge. 0. ) goto 670 r = a q = 0. do 440 i = 1,ny if ( y(i) .gt. a ) goto 480 if ( y(i) .le. q ) goto 440 if ( y(i) .eq. a ) goto 440 q = y(i) 440 continue if ( a .le. a8*q ) goto 670 q = a p = fa 450 nd = nd + 1 if ( nd .gt. 25 ) goto 720 if ( s .eq. q ) goto 780 q = a3*q if ( q .gt. s ) q = s p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) f1 = fa call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p .lt. f1 ) goto 450 if ( a .gt. r ) goto 470 do 460 i = 1,n 460 h(i,2) = x(i) + a*h(i,1) goto 670 470 da = pd(a,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb) if ( da .gt. a6*g ) goto 520 goto 670 480 q = y(i) do 490 j = i,ny if ( y(j) .le. a ) goto 490 if ( y(j) .lt. q ) q = y(j) 490 continue if ( q .le. a5*a ) goto 670 f0 = dlog(q/a) f0 = 1.001*dexp(f0/idint(f0*l3+.999)) v = a 500 v = v*f0 if ( v .ge. q ) goto 430 p = pv(v,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) f1 = fa call igs(v,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p .lt. f1 ) goto 500 if ( a .gt. r ) goto 430 do 510 i = 1,n 510 h(i,2) = x(i) + a*h(i,1) goto 670 520 b = 0. j = 1 i = j 530 i = i + 1 if ( i .gt. ny ) goto 540 if ( y(i) .ge. a ) goto 530 if ( y(i) .lt. b ) goto 530 b = y(i) j = i goto 530 540 fb = z(j) db = d if ( b .ne. 0. ) 1 db = pd2(b,x,h,df,dg,gg,dd,nl,le,vlg,grf,grg,pn,n,lg,bb) 550 w = 2.*dabs(b-a) call cub(c,a,b,fa,fb,da,db) nc = 1 goto 590 560 w = .5*w if ( w .lt. dabs(c0-c) ) goto 660 if ( c0 .lt. c ) goto 570 if ( d0 .ge. d ) goto 580 goto 660 570 if ( d0 .gt. d ) goto 660 580 call cub(c,c,c0,f,f0,d,d0) nc = nc + 1 if ( nc .gt. 30 ) goto 710 590 r = dmax1(a,b) s = dmin1(a,b) if ( c .gt. r ) goto 600 if ( c .gt. s ) goto 610 c = s + (s-c) s = .5*(a+b) if ( c .gt. s ) c = s goto 610 600 c = r - (c-r) s = .5*(a+b) if ( c .lt. s ) c = s 610 c0 = a f0 = fa d0 = da f = pv(c,x,h,vf,vg,dd,nl,le,vlf,vlg,p2,n) d = pd(c,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb) if ( f .lt. fa ) goto 620 b = c fb = f db = d goto 560 620 if ( c .lt. a ) goto 650 if ( d .lt. 0. ) goto 640 630 b = a fb = fa db = da 640 a = c fa = f da = d if ( d .gt. a6*g ) goto 560 goto 670 650 if ( d .lt. 0. ) goto 630 goto 640 660 c = .5*(a+b) nb = nb + 1 w = dabs(b-a) goto 610 670 do 680 i = 1,n 680 x(i) = h(i,2) i = lm1 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc) it = it + 1 f = fa d = da a = a7*a if ( io .gt. 0 ) goto 770 if ( e .le. tl ) goto 770 if ( it .ge. lm2 ) goto 770 if ( it .ge. lm1 ) goto 770 if ( i .lt. lm1 ) goto 100 r = 0. do 690 i = 1,n 690 r = r + h(i,2)*h(i,3) if ( r .lt. 0. ) goto 730 s = r/g g = r d = 0. do 700 i = 1,n h(i,1) = -h(i,2) + s*h(i,1) 700 d = d + h(i,1)*h(i,3) goto 120 710 if ( d .lt. g ) goto 670 io = 4 return 720 io = 3 return 730 io = 4 return 740 q = q*a3**25 nd = 0 750 nd = nd + 1 if ( nd .gt. 25 ) goto 760 q = a3*q if ( q .ge. s ) goto 760 p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .gt. v*q ) goto 750 goto 220 760 io = 5 return 770 st = a return 780 call mod(u,ia,k,io,aa,la,m,ib,1,bb) do 790 i = 1,n 790 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) a = q*a7 800 call grf(df,x) call grg(dg,x) 810 do 820 i = 1,nl 820 bb(i) = le(i) + pn*(vg(i)-dd(i)) do 840 j = 1,n s = df(j) do 830 i = 1,nl 830 s = s + bb(i)*dg(i,j) h(j,3) = s 840 continue f = p 850 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc) it = it + 1 if ( it .ge. lm2 ) goto 770 if ( it .ge. lm1 ) goto 770 if ( io .gt. 0 ) goto 770 if ( e .le. tl ) goto 770 goto 100 860 call mod(u,ia,k,io,aa,la,m,ib,1,bb) q = -s p = pv(q,x,h,vf,vg,dd,nl,le,vlf,vlg,p2,n) do 870 i = 1,n 870 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) goto 800 880 call mod(u,ia,k,io,aa,la,m,ib,1,bb) do 890 i = 1,n 890 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) p = fa goto 810 end subroutine cgp(x,e,it,st,io,tl,lm1,lm2,n,vlf,vlg,grf,grg,l,nl,pn, 1 aa,la,bl,bu,ct,u,ia,k,le,vf,vg,df,dg,lg,uu,gg,bb,cc,h) integer ia(1),i,ib,io,it,j,k,l,la,lg,lm1,lm2,n,nl,na,nb,nc,nd real*8 df(1),dg(lg,1),gg(1),h(n,1),le(1),vg(1),x(1),y(50),z(50) real*8 aa(1),bb(1),bl(1),bu(1),cc(1),u(1),uu(1) real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0,da,db real*8 e,f,f0,f1,fa,fb,fc,g,l3,p,p2,pn,q,r,s,st,tl,v,vf,vv,w real*8 qv,qd,qd2 external grf,grg,vlf,vlg data a1/.1d0/,a2/.9d0/,a3/5.d0/,a4/.2d0/,a5/10.d0/,a6/.9d0/ data a7/.3d0/ c c set big = largest floating point number c i = l if ( l .eq. 0 ) i = 1 call fab(uu,u,ia,io,aa,la,dg,lg,pn,n,i,l,nl,bb,cc) if ( io .gt. 0 ) return big = 1.d50 a8 = a3 + .01d0 p2 = .5*pn it = 0 io = 0 f = vf do 10 i = 1,nl cc(i) = le(i) + pn*vg(i) 10 f = f + (le(i)+p2*vg(i))*vg(i) do 30 j = 1,n s = df(j) do 20 i = 1,nl 20 s = s + cc(i)*dg(i,j) h(j,3) = s 30 continue l3 = 1./dlog(a3) call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc) if ( e .le. tl ) goto 720 if ( io .gt. 0 ) goto 720 if ( lm1 .eq. 0 ) return a = st if ( a .gt. 0. ) goto 50 do 40 i = 1,n 40 if ( dabs(x(i)) .gt. a ) a = dabs(x(i)) a = .01*a/e if ( a .eq. 0. ) a = 1. 50 g = 0. do 60 i = 1,n h(i,1) = -h(i,2) 60 g = g + h(i,2)*h(i,3) if ( g .lt. 0. ) goto 680 d = -g 70 call cut(s,ib,x,h,ia,bl,bu,ct,n,big) na = 0 nb = 0 nc = 0 nd = 0 if ( s .le. 0. ) goto 790 if ( a .gt. s ) a = s p = qv(a,x,h,vv,gg,nl,le,vlf,vlg,p2,n) fa = p c0 = a f0 = fa ny = 2 y(1) = 0. z(1) = f y(2) = a z(2) = fa v = a1*d w = a2*d iq = 0 if ( fa .le. f ) goto 80 c = a b = 0. a = 0. fc = fa fb = f fa = f goto 100 80 c = 0. b = 0. fc = f fb = f iq = 1 vf = vv do 90 i = 1,nl 90 vg(i) = gg(i) 100 q = (d+(f-f0)/c0)/c0 if ( q .lt. 0. ) goto 120 q = a 110 nd = nd + 1 if ( nd .gt. 25 ) goto 670 if ( s .eq. q ) goto 730 q = a3*q if ( q .gt. s ) q = s p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 110 goto 300 120 q = .5*d/q if ( q .lt. s ) goto 140 if ( c0 .lt. s ) goto 130 f1 = f0 c1 = c0 q = s goto 170 130 q = s 140 if ( q .lt. .01*c0 ) q = .01*c0 p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) if ( p .le. f0 ) goto 150 f1 = f0 c1 = c0 f0 = p c0 = q goto 160 150 f1 = p c1 = q 160 call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) 170 if ( a .eq. 0. ) goto 180 if ( fa-f .ge. v*a ) goto 200 if ( fa-f .lt. w*a ) goto 250 goto 320 180 q = c0 if ( c1 .lt. q ) q = c1 190 na = na + 1 if ( na .gt. 25 ) goto 690 q = a4*q p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .ge. v*q ) goto 190 goto 290 200 if ( c0 .gt. c1 ) goto 240 if ( f0-f .gt. v*c0 ) goto 220 if ( f0-f .ge. w*c0 ) goto 380 if ( c1 .le. a5*c0 ) goto 380 r = dlog(c1/c0) r = .999*dexp(-r/idint(r*l3+.999)) q = c1 210 q = q*r if ( q .lt. c0 ) goto 380 p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) na = na + 1 if ( p-f .gt. v*q ) goto 210 goto 380 220 q = c0 230 na = na + 1 if ( na .gt. 25 ) goto 690 q = a4*q p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .ge. v*q ) goto 230 goto 290 240 q = a goto 230 250 if ( c0 .lt. c1 ) goto 350 if ( f0-f .ge. v*c0 ) goto 270 if ( f0-f .ge. w*c0 ) goto 290 q = c0 p = f0 260 nd = nd + 1 if ( nd .gt. 25 ) goto 670 if ( s .eq. q ) goto 730 q = a3*q if ( q .gt. s ) q = s p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 260 goto 290 270 if ( c0 .le. a5*c1 ) goto 290 r = dlog(c0/c1) r = 1.001*dexp(r/idint(r*l3+.999)) q = a 280 q = q*r if ( q .gt. c0 ) goto 290 nd = nd + 1 p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 280 290 if ( iq .eq. 1 ) goto 380 300 if ( b .eq. 0. ) goto 320 if ( c .eq. 0. ) goto 310 v = c - a w = a - b r = 1./v f0 = 1./w p = fc - fa q = fb - fa e = p*r + q*f0 if ( dsign(e,c-b) .ne. e ) goto 380 if ( e .eq. 0. ) goto 380 q = (p*r)*w - (q*f0)*v q = a - .5*q/e goto 330 310 r = 1./a f0 = 1./b p = r*(fa-f) - d q = f0*(fb-f) - d e = a - b v = (r*p-f0*q)/e w = (a*q*f0-b*p*r)/e v = w*w-3.*v*d if ( v .lt. 0. ) v = 0. v = dsqrt(v) if ( w+v .eq. 0. ) goto 380 q = -d/(w+v) if ( q .le. 0. ) goto 380 goto 330 320 if ( iq .eq. 1 ) goto 380 q = (d+(f-fa)/a)/a if ( q .ge. 0. ) goto 380 q = .5*d/q 330 if ( q .gt. s ) q = s do 340 i = 1,ny 340 if ( y(i) .eq. q ) goto 380 p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) goto 380 350 continue if ( f0-f .gt. v*c0 ) goto 360 if ( f0-f .gt. w*c0 ) goto 380 360 q = a p = fa 370 nd = nd + 1 if ( nd .gt. 25 ) goto 670 if ( s .eq. q ) goto 730 q = a3*q if ( q .gt. s ) q = s p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .lt. w*q ) goto 370 goto 290 380 da = qd(a,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc) if ( da .gt. a6*g ) goto 470 if ( s .eq. a ) goto 800 if ( da .ge. 0. ) goto 620 r = a q = 0. do 390 i = 1,ny if ( y(i) .gt. a ) goto 430 if ( y(i) .le. q ) goto 390 if ( y(i) .eq. a ) goto 390 q = y(i) 390 continue if ( a .le. a8*q ) goto 620 q = a p = fa 400 nd = nd + 1 if ( nd .gt. 25 ) goto 670 if ( s .eq. q ) goto 730 q = a3*q if ( q .gt. s ) q = s p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) f1 = fa call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p .lt. f1 ) goto 400 if ( a .gt. r ) goto 420 do 410 i = 1,n 410 h(i,2) = x(i) + a*h(i,1) goto 620 420 da = qd(a,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc) if ( da .gt. a6*g ) goto 470 goto 620 430 q = y(i) do 440 j = i,ny if ( y(j) .le. a ) goto 440 if ( y(j) .lt. q ) q = y(j) 440 continue if ( q .le. a5*a ) goto 620 f0 = dlog(q/a) f0 = 1.001*dexp(f0/idint(f0*l3+.999)) v = a 450 v = v*f0 if ( v .ge. q ) goto 380 p = qv(v,x,h,vv,gg,nl,le,vlf,vlg,p2,n) f1 = fa call igs(v,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p .lt. f1 ) goto 450 if ( a .gt. r ) goto 380 do 460 i = 1,n 460 h(i,2) = x(i) + a*h(i,1) goto 620 470 b = 0. j = 1 i = j 480 i = i + 1 if ( i .gt. ny ) goto 490 if ( y(i) .ge. a ) goto 480 if ( y(i) .lt. b ) goto 480 b = y(i) j = i goto 480 490 fb = z(j) db = d if ( b .ne. 0. ) 1 db = qd2(b,x,h,df,dg,gg,nl,le,vlg,grf,grg,pn,n,lg,cc) 500 w = 2.*dabs(b-a) call cub(c,a,b,fa,fb,da,db) nc = 1 goto 540 510 w = .5*w if ( w .lt. dabs(c0-c) ) goto 610 if ( c0 .lt. c ) goto 520 if ( d0 .ge. d ) goto 530 goto 610 520 if ( d0 .gt. d ) goto 610 530 call cub(c,c,c0,f,f0,d,d0) nc = nc + 1 if ( nc .gt. 30 ) goto 660 540 r = dmax1(a,b) s = dmin1(a,b) if ( c .gt. r ) goto 550 if ( c .gt. s ) goto 560 c = s + (s-c) s = .5*(a+b) if ( c .gt. s ) c = s goto 560 550 c = r - (c-r) s = .5*(a+b) if ( c .lt. s ) c = s 560 c0 = a f0 = fa d0 = da f = qv(c,x,h,vf,vg,nl,le,vlf,vlg,p2,n) d = qd(c,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc) if ( f .lt. fa ) goto 570 b = c fb = f db = d goto 510 570 if ( c .lt. a ) goto 600 if ( d .lt. 0. ) goto 590 580 b = a fb = fa db = da 590 a = c fa = f da = d if ( d .gt. a6*g ) goto 510 goto 620 600 if ( d .lt. 0. ) goto 580 goto 590 610 c = .5*(a+b) nb = nb + 1 w = dabs(b-a) goto 560 620 do 630 i = 1,n 630 x(i) = h(i,2) i = lm1 call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc) it = it + 1 f = fa d = da a = a7*a if ( io .gt. 0 ) goto 720 if ( e .le. tl ) goto 720 if ( it .ge. lm2 ) goto 720 if ( it .ge. lm1 ) goto 720 if ( i .lt. lm1 ) goto 50 r = 0. do 640 i = 1,n 640 r = r + h(i,2)*h(i,3) if ( r .lt. 0. ) goto 680 s = r/g g = r d = 0. do 650 i = 1,n h(i,1) = -h(i,2) + s*h(i,1) 650 d = d + h(i,1)*h(i,3) goto 70 660 if ( d .lt. g ) goto 620 io = 4 return 670 io = 3 return 680 io = 4 return 690 q = q*a3**25 nd = 0 700 nd = nd + 1 if ( nd .gt. 25 ) goto 710 q = a3*q if ( q .ge. s ) goto 710 p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n) call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z) if ( p-f .gt. v*q ) goto 700 goto 170 710 io = 5 return 720 st = a return 730 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc) if ( io .gt. 0 ) goto 720 call mod(u,ia,k,io,aa,la,l,ib,1,bb) do 740 i = 1,n 740 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) a = q*a7 745 call grf(df,x) call grg(dg,x) 746 do 750 i = 1,nl 750 cc(i) = le(i) + pn*vg(i) do 770 j = 1,n s = df(j) do 760 i = 1,nl 760 s = s + cc(i)*dg(i,j) h(j,3) = s 770 continue f = p 780 call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc) it = it + 1 if ( it .ge. lm2 ) goto 720 if ( it .ge. lm1 ) goto 720 if ( io .gt. 0 ) goto 720 if ( e .le. tl ) goto 720 goto 50 790 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc) if ( io .gt. 0 ) goto 720 call mod(u,ia,k,io,aa,la,l,ib,1,bb) q = -s p = qv(q,x,h,vf,vg,nl,le,vlf,vlg,p2,n) do 795 i = 1,n 795 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) goto 745 800 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc) if ( io .gt. 0 ) goto 720 call mod(u,ia,k,io,aa,la,l,ib,1,bb) do 810 i = 1,n 810 x(i) = h(i,2) x(ib) = bl(ib) if ( ia(ib) .gt. 0 ) x(ib) = bu(ib) goto 746 end subroutine cub(x,a,b,c,d,e,f) real*8 a,b,c,d,e,f,g,v,w,x,y,z g = b - a v = e + f - 3*(d-c)/g w = v*v-e*f if ( w .lt. 0. ) w = 0. w = dsign(dsqrt(w),g) y = e + v z = f + v if ( dsign(y,g) .ne. y ) goto 30 if ( dsign(z,g) .ne. z ) goto 20 if ( z .eq. 0. ) goto 20 10 x = b - g*f/(z+w) return 20 if ( c .lt. d ) x = a if ( c .ge. d ) x = b return 30 if ( dsign(z,g) .ne. z ) goto 40 if ( dabs(e) .gt. dabs(f) ) goto 10 40 x = a + g*e/(y-w) return end subroutine cut(s,ib,x,d,ia,bl,bu,ct,n,big) real*8 bl(1),bu(1),d(1),x(1),big,ct,s,t integer ia(1),i,ib,k s = big k = 0 do 20 i = 1,n if ( ia(i) .ne. 0 ) goto 20 if ( d(i) .eq. 0. ) goto 20 if ( d(i) .gt. 0. ) goto 10 t = x(i) - bl(i) if ( t .le. ct ) k = 1 t = -t/d(i) if ( t .gt. s ) goto 20 s = t ib = -i goto 20 10 t = bu(i) - x(i) if ( t .le. ct ) k = 1 t = t/d(i) if ( t .gt. s ) goto 20 s = t ib = i 20 continue if ( s .lt. 0. ) return if ( k .eq. 1 ) goto 30 s = dabs(s) return 30 s = -s return end subroutine erg(t,g,nl) real*8 g(1),t integer i,nl t = 0. do 10 i = 1,nl 10 t = dmax1(dabs(g(i)),t) return end subroutine err(io) integer io goto (10,20,30,40,50,60,70,80,90), io write(6,*) 'no detected errors' return 10 write(6,*) 'dependent constraint gradients' return 20 write(6,*) 'problem is infeasible' return 30 write(6,*) 'function decreases with no minimum' return 40 write(6,*) 'inconsistent gradient values' return 50 write(6,*) 'unable to satisfy Armijo condition' return 60 write(6,*) 'cannot reduce constraint error below gradient error' return 70 write(6,*) 'input penalty is zero but a positive value is needed' write(6,*) 'for convergence in this problem' return 80 write(6,*) 'more constraints imposed than there are unknowns' return 90 write(6,*) 'iterations at limit' return end subroutine fab(v,u,ia,io,a,la,dg,lg,pn,n,ll,l,nl,b,w) integer ia(1),i,io,j,k,l,ll,la,lg,nl,n,o real*8 a(la,1),dg(lg,1),b(1),u(1),v(1),w(ll,1),pn,t do 10 j = 1,n do 10 i = 1,nl 10 a(l+i,j) = dg(i,j) j = (nl*nl+nl)/2 do 20 i = 1,j 20 v(i) = 0. o = 1 do 30 i = 1,nl v(o) = 1./pn o = o + nl - i + 1 30 continue do 60 j = 1,n if ( ia(j) .ne. 0 ) goto 60 o = 0 do 50 i = 1,nl t = dg(i,j) do 40 k = i,nl o = o + 1 v(o) = v(o) + t*dg(k,j) 40 continue 50 continue 60 continue if ( l .eq. 0 ) goto 140 do 70 j = 1,nl do 70 i = 1,l 70 w(i,j) = 0. do 100 j = 1,n if ( ia(j) .ne. 0 ) goto 100 do 90 i = 1,nl t = dg(i,j) do 80 k = 1,l w(k,i) = w(k,i) + t*a(k,j) 80 continue 90 continue 100 continue o = 0 do 130 j = 1,nl do 110 i = 1,l 110 b(i) = w(i,j) call sol(b,u,l,b) do 130 i = j,nl t = 0. do 120 k = 1,l 120 t = t + b(k)*w(k,i) o = o + 1 v(o) = v(o) - t 130 continue 140 call fac(v,0,nl,io) return end subroutine fac(a,p,n,io) real*8 a(1),s,t integer e,f,g,h,i,io,j,k,l,m,n,o,p,q io = 0 if ( n .eq. 0 ) return q = ((n+n+1-p)*p)/2 l = p if ( l .lt. 1 ) l = 1 h = n k = 1 10 s = a(k) if ( s .le. 0. ) goto 50 if ( k .gt. q ) s = dsqrt(s) a(k) = s if ( h .eq. 1 ) return c -------------------------- c |*** save pivot entry ***| c -------------------------- e = k + l f = k + h - 1 do 20 i = e,f 20 a(i) = a(i)/s o = l if ( l .gt. 1 ) l = l - 1 k = k + h g = k h = h - 1 m = h j = 0 30 if ( o .gt. 0 ) o = o - 1 e = g + o j = j - m m = m - 1 f = g + m t = a(g+j) c --------------------------- c |*** eliminate by rows ***| c --------------------------- do 40 i = e,f 40 a(i) = a(i) - t*a(i+j) g = f + 1 if ( m .gt. 0 ) goto 30 goto 10 50 io = 1 return end subroutine fat(u,io,in,ia,a,la,l,n,g,h,z,dg,lg,nl,m) integer ia(1),i,in,io,j,k,l,lg,m,n,nl,o,p real*8 a(la,1),dg(lg,1),g(1),h(1),u(1),z(1),t c c normalize rows of g c do 10 i = 1,nl 10 h(i) = 0. do 20 j = 1,n do 20 i = 1,nl 20 if ( dabs(dg(i,j)) .gt. h(i) ) h(i) = dabs(dg(i,j)) do 30 i = 1,nl if ( h(i) .eq. 0. ) goto 120 h(i) = 1./h(i) 30 z(i) = h(i) do 40 j = 1,n do 40 i = 1,nl 40 a(i+l,j) = h(i)*dg(i,j) do 50 i = 1,nl 50 h(i) = -h(i)*g(i) o = l if ( l .eq. 0 ) goto 110 p = o do 70 k = 1,l do 60 i = 1,nl o = o + 1 60 u(o) = 0. p = p - 1 o = o + p 70 continue do 100 j = 1,n if ( ia(j) .ne. 0 ) goto 100 o = l p = o do 90 k = 1,l t = a(k,j) do 80 i = 1,nl o = o + 1 80 u(o) = u(o) + t*a(i+l,j) p = p - 1 o = o + p 90 continue 100 continue 110 call mat(u(o+1),ia,a(l+1,1),la,nl,n,in) call fac(u,l,m,io) return 120 io = 1 return end subroutine fea(x,u,ia,k,io,y,a,la,m,n,b,bl,bu,c,d,g) integer ia(1),h,i,io,it,j,k,l,la,m,n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),u(1) double precision x(1),y(1),b1,b2,big,e,f,q,r,s,t,v,w c c set big = largest floating point number c big = 1.d50 io = 0 e = 1.d-3 f = .999d0 if ( m .gt. 0 ) goto 30 do 20 i = 1,n ia(i) = 0 x(i) = y(i) if ( x(i) .lt. bu(i) ) goto 10 x(i) = bu(i) ia(i) = 1 goto 20 10 if ( x(i) .gt. bl(i) ) goto 20 x(i) = bl(i) ia(i) = -1 20 continue return 30 call sol(b,u,m,b) do 50 j = 1,n if ( ia(j) .ne. 0 ) goto 50 t = y(j) do 40 i = 1,m 40 t = t + b(i)*a(i,j) x(j) = t 50 continue l = 0 it = -1 do 70 i = 1,n d(i) = 0. if ( iabs(ia(i)) .eq. 1 ) goto 70 if ( x(i) .lt. bl(i) ) goto 60 if ( x(i) .le. bu(i) ) goto 70 ia(i) = 2 l = l + 1 goto 70 60 ia(i) = -2 l = l + 1 70 continue b2 = 0. if ( l .gt. 0 ) goto 190 do 80 i = 1,n if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 80 continue return 90 r = 0. s = 0. t = big h = 0 if ( it .eq. 0 ) q = 0. if ( it .gt. 0 ) q = b2/b1 do 120 i = 1,n if ( iabs(ia(i)) .eq. 1 ) goto 120 d(i) = q*d(i) - g(i) s = s + d(i)*g(i) w = bu(i) if ( d(i) .gt. 0. ) goto 100 if ( d(i) .eq. 0. ) goto 120 if ( ia(i) .lt. 0 ) goto 120 if ( ia(i) .eq. 0 ) w = bl(i) if ( ia(i) .gt. 0 ) r = r + d(i)**2 goto 110 100 if ( ia(i) .gt. 0 ) goto 120 if ( ia(i) .eq. 0 ) goto 110 r = r + d(i)**2 w = bl(i) 110 v = (w-x(i))/d(i) if ( v .ge. t ) goto 120 t = v h = i 120 continue s = -s/r if ( t .le. s ) goto 130 if ( h .eq. 0 ) goto 170 if ( ia(h) .eq. 0 ) goto 170 130 s = t if ( d(h) .gt. 0. ) goto 140 if ( ia(h) .gt. 0 ) goto 150 h = -h goto 160 140 if ( ia(h) .eq. 0 ) goto 160 h = -h 150 l = l - 1 160 call mod(u,ia,k,io,a,la,m,h,1,b) if ( io .gt. 0 ) return it = -1 170 do 180 i = 1,n if ( iabs(ia(i)) .eq. 1 ) goto 180 x(i) = x(i) + s*d(i) 180 continue if ( l .eq. 0 ) goto 390 c c *** compute the projected gradient *** c 190 do 200 i = 1,m 200 c(i) = 0. do 220 j = 1,n if ( iabs(ia(j)) .lt. 2 ) goto 220 if ( ia(j) .gt. 0 ) t = e*bl(j) + f*bu(j) - x(j) if ( ia(j) .lt. 0 ) t = e*bu(j) + f*bl(j) - x(j) do 210 i = 1,m 210 c(i) = c(i) + t*a(i,j) 220 continue call sol(b,u,m,c) it = it + 1 b1 = b2 b2 = 0. r = 0. s = 0. do 250 j = 1,n t = 0. if ( ia(j) .eq. 2 ) t = x(j) - e*bl(j) - f*bu(j) if ( ia(j) .eq.-2 ) t = x(j) - e*bu(j) - f*bl(j) do 230 i = 1,m 230 t = t + b(i)*a(i,j) if ( iabs(ia(j)) .eq. 1 ) goto 240 b2 = b2 + t*t g(j) = t r = dmax1(dabs(t),r) goto 250 240 t = t*ia(j) if ( t .le. s ) goto 250 s = t h = j 250 continue if ( s .gt. 10.*r ) goto 260 if ( b2 .eq. 0. ) goto 290 if ( it .ge. l ) goto 290 if ( it .ge. k ) goto 290 goto 90 260 x(h) = bl(h) if ( ia(h) .gt. 0 ) x(h) = bu(h) call mod(u,ia,k,io,a,la,m,h,0,b) if ( io .gt. 0 ) return call sol(b,u,m,c) it = 0 b2 = 0. do 280 j = 1,n if ( iabs(ia(j)) .eq. 1 ) goto 280 t = 0. if ( ia(j) .eq. 2 ) t = x(j) - e*bl(j) - f*bu(j) if ( ia(j) .eq.-2 ) t = x(j) - e*bu(j) - f*bl(j) do 270 i = 1,m 270 t = t + b(i)*a(i,j) b2 = b2 + t*t g(j) = t 280 continue if ( b2 .eq. 0. ) goto 190 goto 90 290 if ( k .ge. l ) goto 300 io = 2 return 300 do 310 i = 1,m 310 b(i) = 0. do 320 i = 1,n if ( ia(i) .eq. 0 ) g(i) = x(i) if ( ia(i) .eq. 1 ) g(i) = bu(i) if ( ia(i) .eq.-1 ) g(i) = bl(i) if ( iabs(ia(i)) .eq. 2 ) g(i) = x(i) 320 continue do 340 j = 1,n if ( iabs(ia(j)) .lt. 2 ) goto 340 call mod(u,ia,k,io,a,la,m,j,ia(j),b) if ( io .gt. 0 ) return if ( ia(j) .lt. 0 ) t = x(j) - bl(j) if ( ia(j) .gt. 0 ) t = x(j) - bu(j) do 330 i = 1,m 330 b(i) = b(i) + t*a(i,j) 340 continue call sol(b,u,m,b) do 380 j = 1,n if ( ia(j) .ne. 0 ) goto 380 t = x(j) do 350 i = 1,m 350 t = t + b(i)*a(i,j) x(j) = t if ( ia(j) .ne. 0 ) goto 380 if ( t .lt. bl(j) ) goto 360 if ( t .le. bu(j) ) goto 380 360 io = 2 do 370 i = 1,n 370 x(i) = g(i) return 380 continue 390 do 400 i = 1,m 400 c(i) = 0. do 420 j = 1,n if ( ia(j) .ne. 0 ) goto 420 x(j) = x(j) - y(j) t = x(j) do 410 i = 1,m 410 c(i) = c(i) + t*a(i,j) 420 continue 430 call sol(b,u,m,c) 440 s = 1. h = 0 do 480 j = 1,n if ( ia(j) .ne. 0 ) goto 480 t = 0. do 450 i = 1,m 450 t = t + b(i)*a(i,j) g(j) = t r = t + y(j) if ( r .lt. bu(j) ) goto 460 r = (bu(j)-x(j)-y(j))/(t-x(j)) if ( r .ge. s ) goto 470 s = r h = j goto 470 460 if ( r .gt. bl(j) ) goto 470 r = (bl(j)-x(j)-y(j))/(t-x(j)) if ( r .ge. s ) goto 470 s = r h = -j 470 g(j) = t 480 continue if ( h .eq. 0 ) goto 510 call mod(u,ia,k,io,a,la,m,h,1,b) if ( io .gt. 0 ) return if ( ia(h) .gt. 0 ) t = bu(h) - y(h) if ( ia(h) .lt. 0 ) t = bl(h) - y(h) do 490 i = 1,m 490 c(i) = c(i) - t*a(i,h) call sol(b,u,m,c) do 500 i = 1,n 500 if ( ia(i) .eq. 0 ) x(i) = x(i) + s*(g(i)-x(i)) goto 530 510 do 520 i = 1,n 520 if ( ia(i) .eq. 0 ) x(i) = g(i) 530 s = 0. do 550 j = 1,n if ( ia(j) .eq. 0 ) goto 550 if ( ia(j) .eq. 1 ) t = bu(j) - y(j) if ( ia(j) .eq.-1 ) t = bl(j) - y(j) do 540 i = 1,m 540 t = t - b(i)*a(i,j) t = t*ia(j) if ( t .le. s ) goto 550 s = t l = j 550 continue if ( s .eq. 0 ) goto 570 if ( ia(l) .eq. 1 ) t = bu(l) - y(l) if ( ia(l) .eq.-1 ) t = bl(l) - y(l) x(l) = t call mod(u,ia,k,io,a,la,m,l,0,b) if ( io .gt. 0 ) return do 560 i = 1,m 560 c(i) = c(i) + t*a(i,l) goto 430 570 if ( h .gt. 0 ) goto 440 do 580 i = 1,n if ( ia(i) .eq. 0 ) x(i) = x(i) + y(i) if ( ia(i) .eq. 1 ) x(i) = bu(i) 580 if ( ia(i) .eq.-1 ) x(i) = bl(i) return end c c in = 0 -- all constraints are nonbinding initially c in = 1 -- binding constraints are given by values of ia c in = 2 -- binding constraints determined by x, bl, and bu c in = 3 -- use k,ia, and u as given c c LENGTH OF W AT LEAST 2M+2N IF NF1 INVOKED C 4M+2N IF NF2 INVOKED c subroutine feas(x,u,ia,k,io,in,y,a,la,m,n,b,bl,bu,w) integer ia(1),i,in,io,ip,j,k,l,la,m,n real*8 a(la,1),b(1),bl(1),bu(1),u(1),x(1),y(1),w(1),t ip = 0 if ( in .eq. 3 ) goto 100 k = n - m do 10 i = 1,n if ( bl(i) .lt. bu(i) ) goto 10 t = bl(i) bl(i) = bu(i) bu(i) = t 10 continue if ( in .eq. 2 ) goto 50 if ( in .eq. 1 ) goto 30 do 20 i = 1,n 20 ia(i) = 0 goto 90 30 do 40 i = 1,n 40 if ( ia(i) .ne. 0 ) k = k - 1 if ( k .ge. 0 ) goto 90 goto 80 50 do 70 i = 1,n ia(i) = 0 if ( y(i) .ne. bu(i) ) goto 60 ia(i) = 1 k = k - 1 goto 70 60 if ( y(i) .ne. bl(i) ) goto 70 ia(i) = -1 k = k - 1 70 continue if ( k .ge. 0 ) goto 90 80 k = k + m ip = 1 90 call mat(u,ia,a,la,m,n,ip) call fac(u,0,m,io) if ( io .gt. 0 ) return 100 if ( m .eq. 0 ) goto 130 do 110 i = 1,m 110 w(i) = b(i) do 120 j = 1,n t = y(j) do 120 i = 1,m 120 w(i) = w(i) - t*a(i,j) 130 i = m + 1 j = i + m if ( ip .eq. 1 ) goto 140 l = j + n call nf1(x,u,ia,k,io,y,a,la,m,n,w,bl,bu,w(i),w(j),w(l)) if ( m .gt. 0 ) goto 160 return 140 do 150 l = 1,m 150 w(l+m) = 1. l = j + m call nf2(x,u,ia,k,io,y,a,la,m,n,w,w(i),bl,bu,w(j),w(l),w(l+m+n)) 160 if ( io .gt. 0 ) return call lsq(x,u,ia,k,io,y,a,la,m,n,bl,bu,w,w(i),w(j)) return end double precision function fv(a,x,h,n,vl) integer i,n real*8 h(n,1),x(1),a,vl do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) fv = vl(h(1,2)) return end c double precision function pv(a,x,h,f,g,d,nl,le,vlf,vlg,p2,n) integer i,nl,n real*8 d(1),g(1),h(n,1),le(1),x(1),a,f,p2,t,vlf do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) f = vlf(h(1,2)) call vlg(g,h(1,2)) t = f do 20 i = 1,nl 20 t = t + le(i)*g(i) + p2*(g(i)-d(i))**2 pv = t return end c double precision function qv(a,x,h,f,g,nl,le,vlf,vlg,p2,n) integer i,nl,n real*8 g(1),h(n,1),le(1),x(1),a,f,p2,t,vlf do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) f = vlf(h(1,2)) call vlg(g,h(1,2)) t = f do 20 i = 1,nl 20 t = t + (le(i)+p2*g(i))*g(i) qv = t return end c double precision function fd(a,x,h,n,gr) real*8 h(n,1),x(1),a,d do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) call gr(h(1,3),h(1,2)) d = 0. do 20 i = 1,n 20 d = d + h(i,1)*h(i,3) fd = d return end c double precision function pd(a,x,h,df,dg,g,d,nl,le, 1 grf,grg,pn,n,lg,b) integer i,nl,lg,n real*8 b(1),d(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) call grf(df,h(1,2)) call grg(dg,h(1,2)) do 20 i = 1,nl 20 b(i) = le(i) + pn*(g(i)-d(i)) s = 0. do 40 j = 1,n t = df(j) do 30 i = 1,nl 30 t = t + b(i)*dg(i,j) h(j,3) = t s = s + t*h(j,1) 40 continue pd = s return end c double precision function qd(a,x,h,df,dg,g,nl,le, 1 grf,grg,pn,n,lg,b) integer i,nl,lg,n real*8 b(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) call grf(df,h(1,2)) call grg(dg,h(1,2)) do 20 i = 1,nl 20 b(i) = le(i) + pn*g(i) s = 0. do 40 j = 1,n t = df(j) do 30 i = 1,nl 30 t = t + b(i)*dg(i,j) h(j,3) = t s = s + t*h(j,1) 40 continue qd = s return end c double precision function pd2(a,x,h,df,dg,g,d,nl,le, 1 vlg,grf,grg,pn,n,lg,b) integer i,nl,lg,n real*8 b(1),d(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) call grf(df,h(1,2)) call grg(dg,h(1,2)) call vlg(g,h(1,2)) do 20 i = 1,nl 20 b(i) = le(i) + pn*(g(i)-d(i)) s = 0. do 40 j = 1,n t = df(j) do 30 i = 1,nl 30 t = t + b(i)*dg(i,j) h(j,3) = t s = s + t*h(j,1) 40 continue pd2 = s return end c double precision function qd2(a,x,h,df,dg,g,nl,le, 1 vlg,grf,grg,pn,n,lg,b) integer i,nl,lg,n real*8 b(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t do 10 i = 1,n 10 h(i,2) = x(i) + a*h(i,1) call grf(df,h(1,2)) call grg(dg,h(1,2)) call vlg(g,h(1,2)) do 20 i = 1,nl 20 b(i) = le(i) + pn*g(i) s = 0. do 40 j = 1,n t = df(j) do 30 i = 1,nl 30 t = t + b(i)*dg(i,j) h(j,3) = t s = s + t*h(j,1) 40 continue qd2 = s return end subroutine igs(s,f,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,j,y,z) real*8 gg(1),vg(1),y(1),z(1),a,b,c,f,fa,fb,fc,s,vf,vv integer j j = j + 1 y(j) = s z(j) = f if ( f .le. fa ) goto 20 if ( f .le. fb ) goto 10 if ( f .gt. fc ) return c = s fc = f return 10 c = b b = s fc = fb fb = f return 20 c = b b = a a = s fc = fb fb = fa fa = f vf = vv do 30 i = 1,nl 30 vg(i) = gg(i) return end subroutine ins(s,f,a,b,c,fa,fb,fc,j,y,z) real*8 a,b,c,f,fa,fb,fc,s,y(1),z(1) integer j j = j + 1 y(j) = s z(j) = f if ( f .le. fa ) goto 20 if ( f .le. fb ) goto 10 if ( f .gt. fc ) return c = s fc = f return 10 c = b b = s fc = fb fb = f return 20 c = b b = a a = s fc = fb fb = fa fa = f return end subroutine lsq(x,u,ia,k,io,y,a,la,m,n,bl,bu,b,c,g) integer ia(1),h,i,io,j,k,la,m,n double precision a(la,1),b(1),bl(1),bu(1),c(1),g(1),u(1) double precision x(1),y(1),r,s,t io = 0 if ( m .gt. 0 ) goto 30 do 20 i = 1,n x(i) = y(i) if ( x(i) .lt. bl(i) ) goto 10 if ( x(i) .le. bu(i) ) goto 20 x(i) = bu(i) if ( ia(i) .eq. 0 ) k = k - 1 ia(i) = 1 goto 20 10 x(i) = bl(i) if ( ia(i) .eq. 0 ) k = k - 1 ia(i) = -1 20 continue return 30 do 40 i = 1,m 40 c(i) = 0. do 60 j = 1,n if ( ia(j) .ne. 0 ) goto 60 x(j) = x(j) - y(j) t = x(j) do 50 i = 1,m 50 c(i) = c(i) + t*a(i,j) 60 continue 70 call sol(b,u,m,c) 80 s = 1. h = 0 do 120 j = 1,n if ( ia(j) .ne. 0 ) goto 120 t = 0. do 90 i = 1,m 90 t = t + b(i)*a(i,j) g(j) = t r = t + y(j) if ( r .lt. bu(j) ) goto 100 r = (bu(j)-x(j)-y(j))/(t-x(j)) if ( r .ge. s ) goto 110 s = r h = j goto 110 100 if ( r .gt. bl(j) ) goto 110 r = (bl(j)-x(j)-y(j))/(t-x(j)) if ( r .ge. s ) goto 110 s = r h = -j 110 g(j) = t 120 continue if ( h .eq. 0 ) goto 150 call mod(u,ia,k,io,a,la,m,h,1,b) if ( io .gt. 0 ) return if ( ia(h) .gt. 0 ) t = bu(h) - y(h) if ( ia(h) .lt. 0 ) t = bl(h) - y(h) do 130 i = 1,m 130 c(i) = c(i) - t*a(i,h) call sol(b,u,m,c) do 140 i = 1,n 140 if ( ia(i) .eq. 0 ) x(i) = x(i) + s*(g(i)-x(i)) goto 80 150 do 160 i = 1,n 160 if ( ia(i) .eq. 0 ) x(i) = g(i) 170 s = 0. do 190 j = 1,n if ( ia(j) .eq. 0 ) goto 190 if ( ia(j) .eq. 1 ) t = bu(j) - y(j) if ( ia(j) .eq.-1 ) t = bl(j) - y(j) do 180 i = 1,m 180 t = t - b(i)*a(i,j) t = t*ia(j) if ( t .le. s ) goto 190 s = t l = j 190 continue if ( s .eq. 0 ) goto 210 if ( ia(l) .eq. 1 ) t = bu(l) - y(l) if ( ia(l) .eq.-1 ) t = bl(l) - y(l) x(l) = t call mod(u,ia,k,io,a,la,m,l,0,b) if ( io .gt. 0 ) return do 200 i = 1,m 200 c(i) = c(i) + t*a(i,l) goto 70 210 if ( h .gt. 0 ) goto 80 do 220 i = 1,n if ( ia(i) .eq. 0 ) x(i) = x(i) + y(i) if ( ia(i) .eq. 1 ) x(i) = bu(i) 220 if ( ia(i) .eq.-1 ) x(i) = bl(i) return end subroutine mat(u,ia,a,la,m,n,in) integer ia(1),i,in,j,k,l,la,m,n real*8 a(la,1),u(1),t if ( m .eq. 0 ) return j = (m*m+m)/2 do 10 i = 1,j 10 u(i) = 0. if ( in .eq. 0 ) goto 30 l = 1 do 20 i = 1,m u(l) = 1. l = l + m - i + 1 20 continue 30 do 60 j = 1,n if ( ia(j) .ne. 0 ) goto 60 l = 0 do 50 i = 1,m t = a(i,j) do 40 k = i,m l = l + 1 u(l) = u(l) + t*a(k,j) 40 continue 50 continue 60 continue return end subroutine mob(v,io,u,ia,k,a,la,n,l,nl,ib,g,b,c) integer ia(1),g,h,i,ib,io,j,k,l,la,n,nl real*8 a(la,1),b(1),c(1),u(1),v(1),s,t h = iabs(ib) if ( l .gt. 0 ) goto 20 do 10 i = 1,nl 10 c(i) = a(i+l,h) goto 110 20 do 30 i = 1,l 30 b(i) = a(i,h) call sol(b,u,l,b) t = 0. do 40 i = 1,l 40 t = t - a(i,h)*b(i) s = t + 1. if ( s .gt. 0. ) goto 50 io = 1 return 50 do 60 i = 1,nl 60 c(i) = a(i+l,h) do 90 j = 1,n if ( ia(j) .ne. 0 ) goto 90 t = 0. do 70 i = 1,l 70 t = t - a(i,j)*b(i) do 80 i = 1,nl 80 c(i) = c(i) + t*a(i+l,j) 90 continue s = 1./dsqrt(s) do 100 i = 1,nl 100 c(i) = s*c(i) 110 h = ib call mod(v,ia,k,io,a,la,nl,h,g,c) return end subroutine mod(u,ia,o,io,a,la,n,q,h,v) integer ia(1),h,i,j,k,l,la,m,n,nm1,o,p,q real*8 a(la,1),v(1),u(1),c,r,s,t io = 0 p = h if ( p .ne. 0 ) goto 10 if ( ia(q) .eq. 0 ) return ia(q) = 0 o = o + 1 goto 60 10 if ( p .lt. 3 ) goto 40 if ( p .gt. 3 ) goto 30 do 20 i = 1,n 20 v(i) = 0. v(q) = 1. o = o - 1 goto 80 30 if ( p .eq. 5 ) p = 0 goto 80 40 if ( q .gt. 0 ) goto 50 p = -p q = -q 50 i = ia(q) if ( p .lt. 0 ) ia(q) = -1 if ( p .gt. 0 ) ia(q) = 1 if ( iabs(i) .eq. 1 ) return o = o - 1 if ( o .ge. 0 ) goto 60 io = 1 return 60 if ( n .eq. 0 ) return do 70 i = 1,n 70 v(i) = a(i,q) 80 k = 0 m = n j = 1 90 t = v(j)/u(j+k) v(j) = t if ( j .eq. n ) goto 110 j = j + 1 do 100 i = j,n 100 v(i) = v(i) - t*u(i+k) m = m - 1 k = k + m goto 90 110 r = v(n) l = (n*n+n)/2 + 1 nm1 = n - 1 do 170 k = 1,nm1 m = l - 2 l = m - k s = r j = n - k c = v(j) if ( dabs(c) .gt. dabs(s) ) goto 130 if ( s .eq. 0. ) goto 120 r = dabs(s)*dsqrt(1.+(c/s)**2) goto 140 120 c = 1. s = 0. goto 150 130 r = dabs(c)*dsqrt(1.+(s/c)**2) 140 c = c/r s = s/r c c process rows j and j+1 of U c 150 v(j+1) = -s*u(l) u(l) = c*u(l) l = l + 1 do 160 i = l,m t = u(i+k) u(i+k) = c*t - s*u(i) u(i) = s*t + c*u(i) 160 continue 170 continue if ( p .ne. 0 ) r = 1 - r*r if ( p .eq. 0 ) r = 1 + r*r if ( r .gt. 0. ) goto 180 io = 1 return 180 r = dsqrt(r) do 190 i = 1,n 190 u(i) = r*u(i) c c restore upper triangular form c m = 0 k = n do 250 j = 2,n l = m + 1 m = m + k k = k - 1 c = u(l) s = v(j) if ( dabs(c) .gt. dabs(s) ) goto 210 if ( s .eq. 0. ) goto 200 r = dabs(s)*dsqrt(1.+(c/s)**2) goto 220 200 c = 1. s = 0. r = 0. goto 230 210 r = dabs(c)*dsqrt(1.+(s/c)**2) 220 c = c/r s = s/r 230 u(l) = r l = l + 1 do 240 i = l,m t = u(i+k) u(i+k) = c*t - s*u(i) u(i) = s*t + c*u(i) 240 continue 250 continue return end subroutine mul(le,li,e,u,ia,a,la,m,n,df,l,nl,z) integer ia(1),i,j,l,la,m,n,nl real*8 a(la,1),df(1),le(1),li(1),u(1),z(1),e,t do 10 i = 1,m 10 le(i) = 0. do 30 j = 1,n if ( ia(j) .ne. 0 ) goto 30 t = df(j) do 20 i = 1,m 20 le(i) = le(i) - t*a(i,j) 30 continue call sol(le,u,m,le) e = 0. do 60 j = 1,n li(j) = 0. t = df(j) do 40 i = 1,m 40 t = t + le(i)*a(i,j) if ( ia(j) .eq. 0 ) goto 50 if ( ia(j) .gt. 0 ) t = -t li(j) = t if ( t .ge. 0. ) goto 60 li(j) = 0. 50 e = dmax1(e,dabs(t)) 60 continue if ( nl .eq. 0 ) return do 70 i = 1,nl 70 le(i+l) = le(i+l)*z(i) return end subroutine mup(p,e,u,ia,le,a,la,n,dg,lg,pn,g,df,l,nl,y,z) integer ia(1),i,j,l,la,lg,n,nl real*8 a(la,1),dg(lg,1),df(1),g(1),le(1),p(1),u(1),y(1),z(1) real*8 e,pn,t do 10 i = 1,nl 10 z(i) = le(i) + pn*g(i) if ( l .eq. 0 ) goto 90 do 20 i = 1,l 20 p(i) = 0. do 50 j = 1,n t = df(j) do 30 i = 1,nl 30 t = t + z(i)*dg(i,j) y(j) = t if ( ia(j) .ne. 0 ) goto 50 do 40 i = 1,l 40 p(i) = p(i) - t*a(i,j) 50 continue call sol(p,u,l,p) e = 0. do 80 j = 1,n t = y(j) do 60 i = 1,l 60 t = t + p(i)*a(i,j) if ( ia(j) .eq. 0 ) goto 70 if ( ia(j) .gt. 0 ) t = -t if ( t .ge. 0. ) goto 80 70 e = dmax1(e,dabs(t)) 80 continue return 90 e = 0. do 120 j = 1,n t = df(j) do 100 i = 1,nl 100 t = t + z(i)*dg(i,j) if ( ia(j) .eq. 0 ) goto 110 if ( ia(j) .gt. 0 ) t = -t if ( t .ge. 0. ) goto 120 110 e = dmax1(e,dabs(t)) 120 continue return end subroutine nf1(x,u,ia,k,io,y,a,la,m,n,b,bl,bu,c,d,g) integer ia(1),h,i,io,it,j,k,l,la,m,n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),u(1) double precision x(1),y(1),big,e,f,p,q,r,s,t c c set big = largest floating point number c big = 1.d50 io = 0 if ( m .gt. 0 ) goto 30 do 20 i = 1,n x(i) = y(i) if ( x(i) .lt. bl(i) ) goto 10 if ( x(i) .le. bu(i) ) goto 20 x(i) = bu(i) if ( ia(i) .eq. 0 ) k = k - 1 ia(i) = 1 goto 20 10 x(i) = bl(i) if ( ia(i) .eq. 0 ) k = k - 1 ia(i) = -1 20 continue return 30 call sol(b,u,m,b) do 50 j = 1,n if ( ia(j) .ne. 0 ) goto 50 t = y(j) do 40 i = 1,m 40 t = t + b(i)*a(i,j) x(j) = t 50 continue l = 0 do 70 i = 1,n d(i) = 0. if ( iabs(ia(i)) .eq. 1 ) goto 70 if ( x(i) .lt. bl(i) ) goto 60 if ( x(i) .le. bu(i) ) goto 70 ia(i) = 2 l = l + 1 goto 70 60 ia(i) = -2 l = l + 1 70 continue f = 0. it = 0 if ( l .gt. 0 ) goto 230 io = -1 80 do 90 i = 1,n if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 90 continue return 100 r = 0. s = 0. t = big h = 0 it = it + 1 if ( it .eq. 1 ) q = 0. if ( it .gt. 1 ) q = f/e do 190 i = 1,n if ( iabs(ia(i)) .eq. 1 ) goto 190 d(i) = q*d(i) - g(i) s = s + d(i)*g(i) if ( ia(i) .ne. 0 ) goto 110 if ( d(i) .eq. 0. ) goto 190 if ( d(i) .lt. 0. ) goto 120 goto 170 110 r = r + d(i)**2 if ( ia(i) .gt. 0 ) goto 150 if ( x(i) .ge. bl(i) ) goto 140 if ( d(i) .le. 0. ) goto 190 120 p = (bl(i)-x(i))/d(i) 130 if ( p .ge. t ) goto 190 t = p h = -i goto 190 140 p = 0. goto 130 150 if ( x(i) .gt. bu(i) ) goto 160 p = 0. goto 180 160 if ( d(i) .ge. 0. ) goto 190 170 p = (bu(i)-x(i))/d(i) 180 if ( p .ge. t ) goto 190 t = p h = i 190 continue if ( r .ne. 0. ) s = -s/r j = iabs(h) if ( t .le. s ) goto 200 if ( ia(j) .eq. 0 ) goto 210 200 s = t if ( ia(j) .ne. 0 ) l = l - 1 call mod(u,ia,k,io,a,la,m,h,1,b) if ( io .gt. 0 ) return it = 0 210 do 220 i = 1,n if ( iabs(ia(i)) .eq. 1 ) goto 220 x(i) = x(i) + s*d(i) 220 continue if ( l .eq. 0 ) goto 80 c c *** compute the projected gradient *** c 230 do 240 i = 1,m 240 c(i) = 0. do 260 j = 1,n if ( iabs(ia(j)) .lt. 2 ) goto 260 if ( ia(j) .gt. 0 ) t = bu(j) - x(j) if ( ia(j) .lt. 0 ) t = bl(j) - x(j) do 250 i = 1,m 250 c(i) = c(i) + t*a(i,j) 260 continue call sol(b,u,m,c) e = f f = 0. r = 0. s = 0. do 290 j = 1,n t = 0. if ( ia(j) .eq. 2 ) t = x(j) - bu(j) if ( ia(j) .eq.-2 ) t = x(j) - bl(j) do 270 i = 1,m 270 t = t + b(i)*a(i,j) if ( iabs(ia(j)) .eq. 1 ) goto 280 f = f + t*t g(j) = t r = dmax1(dabs(t),r) goto 290 280 t = t*ia(j) if ( t .le. s ) goto 290 s = t h = j 290 continue 300 if ( s .gt. 10.*r ) goto 310 if ( r .eq. 0. ) goto 360 if ( it .ge. l ) goto 360 if ( it .ge. k ) goto 360 goto 100 310 x(h) = bl(h) if ( ia(h) .gt. 0 ) x(h) = bu(h) call mod(u,ia,k,io,a,la,m,h,0,b) if ( io .gt. 0 ) return call sol(b,u,m,c) it = 0 f = 0. r = 0. do 330 j = 1,n if ( iabs(ia(j)) .eq. 1 ) goto 330 t = 0. if ( ia(j) .eq. 2 ) t = x(j) - bu(j) if ( ia(j) .eq.-2 ) t = x(j) - bl(j) do 320 i = 1,m 320 t = t + b(i)*a(i,j) f = f + t*t g(j) = t r = dmax1(dabs(t),r) 330 continue if ( r .gt. 0. ) goto 100 do 350 j = 1,n if ( iabs(ia(j)) .ne. 1 ) goto 350 t = 0. do 340 i = 1,m 340 t = t + b(i)*a(i,j) t = t*ia(j) if ( t .le. s ) goto 350 s = t h = j 350 continue goto 300 360 io = 2 return end subroutine nf2(x,u,ia,k,io,y,a,la,m,n,b,ib,bl,bu,c,d,g) integer ia(1),h,i,io,it,j,k,l,la,m,n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),ib(1) double precision u(1),x(1),y(1),big,e,f,p,q,r,s,t c c set big = largest floating point number c big = 1.d50 io = 0 do 10 i = 1,n 10 x(i) = y(i) if ( m .eq. 0 ) goto 360 l = 0 do 20 i = 1,m 20 if ( ib(i) .ne. 0 ) l = l + 1 if ( l .eq. 0 ) goto 360 do 30 i = 1,m if ( ib(i) .ne. 0. ) ib(i) = dsign(1.d0,b(i)) 30 if ( ib(i) .eq. 0. ) b(i) = 0. f = 0. it = 0 goto 220 40 t = big it = it + 1 if ( it .eq. 1 ) q = 0. if ( it .gt. 1 ) q = f/e do 60 i = 1,n d(i) = 0. if ( ia(i) .ne. 0 ) goto 60 d(i) = q*d(i) - g(i) if ( d(i) .gt. 0. ) goto 50 if ( d(i) .eq. 0. ) goto 60 p = (bl(i)-x(i))/d(i) if ( p .ge. t ) goto 60 t = p h = -i goto 60 50 p = (bu(i)-x(i))/d(i) if ( p .ge. t ) goto 60 t = p h = i 60 continue r = 0. s = 0. do 70 i = 1,m c(i) = 0. if ( ib(i) .eq. 0 ) goto 70 j = i + n d(j) = q*d(j) - g(j) r = r + d(j)**2 s = s + d(j)*b(i) 70 continue if ( r .ne. 0. ) s = -s/r do 90 j = 1,n r = d(j) if ( r .eq. 0. ) goto 90 do 80 i = 1,m 80 c(i) = c(i) + r*a(i,j) 90 continue r = big do 130 i = 1,m if ( ib(i) .eq. 0. ) goto 130 if ( ib(i) .gt. 0. ) goto 100 if ( b(i) .lt. 0. ) goto 110 p = 0. goto 120 100 if ( b(i) .gt. 0. ) goto 110 p = 0. goto 120 110 if ( c(i) .eq. 0. ) goto 130 p = b(i)/c(i) if ( p .le. 0. ) goto 130 120 if ( p .ge. r ) goto 130 r = p j = i 130 continue if ( t .lt. r ) goto 160 do 140 i = 1,n 140 x(i) = x(i) + r*d(i) do 150 i = 1,m 150 if ( ib(i) .ne. 0. ) b(i) = b(i) - r*c(i) ib(j) = 0. b(j) = 0. call mod(u,ia,k,io,a,la,m,j,3,c) if ( io .gt. 0 ) return l = l - 1 if ( l .eq. 0 ) goto 360 it = 0 goto 220 160 if ( t .le. s ) goto 190 do 170 i = 1,n 170 x(i) = x(i) + s*d(i) do 180 i = 1,m 180 if ( ib(i) .ne. 0. ) b(i) = b(i) - s*c(i) goto 220 190 do 200 i = 1,n 200 x(i) = x(i) + t*d(i) do 210 i = 1,m 210 if ( ib(i) .ne. 0. ) b(i) = b(i) - t*c(i) call mod(u,ia,k,io,a,la,m,h,1,c) it = 0 c c *** compute the projected gradient *** c 220 call sol(c,u,m,b) 230 e = f f = 0. r = 0. s = 0. do 260 j = 1,n t = 0. do 240 i = 1,m 240 t = t - c(i)*a(i,j) if ( iabs(ia(j)) .eq. 1 ) goto 250 f = f + t*t g(j) = t r = dmax1(dabs(t),r) goto 260 250 t = t*ia(j) if ( t .le. s ) goto 260 s = t h = j 260 continue 270 if ( s .gt. 10.*r ) goto 300 if ( r .eq. 0. ) goto 350 if ( it .ge. l ) goto 350 if ( it .ge. k ) goto 350 280 do 290 i = 1,m if ( ib(i) .eq. 0. ) goto 290 t = b(i) - c(i) f = f + t*t g(i+n) = t 290 continue goto 40 300 x(h) = bl(h) if ( ia(h) .gt. 0 ) x(h) = bu(h) call mod(u,ia,k,io,a,la,m,h,0,c) if ( io .gt. 0 ) return call sol(c,u,m,b) it = 0 f = 0. r = 0. do 320 j = 1,n if ( ia(j) .ne. 0 ) goto 320 t = 0. do 310 i = 1,m 310 t = t - c(i)*a(i,j) f = f + t*t g(j) = t r = dmax1(dabs(t),r) 320 continue if ( r .gt. 0. ) goto 280 s = 0. do 340 j = 1,n if ( ia(j) .eq. 0 ) goto 340 t = 0. do 330 i = 1,m 330 t = t - c(i)*a(i,j) t = t*ia(j) if ( t .le. s ) goto 340 s = t h = j 340 continue goto 270 350 io = 2 return 360 do 370 i = 1,n if ( ia(i) .gt. 0 ) x(i) = bu(i) if ( ia(i) .lt. 0 ) x(i) = bl(i) 370 continue return end subroutine pre(p,e,u,ia,k,l,lm,io,g,a,la,m,n,b,c) integer ia(1),i,io,j,k,l,la,lm,m,n real*8 a(la,1),b(1),c(1),g(1),p(1),u(1),e,s,t if ( m .gt. 0 ) goto 60 s = 0. e = s do 20 j = 1,n t = g(j) if ( ia(j) .eq. 0 ) goto 10 if ( ia(j) .lt. 0 ) t = -t if ( t .le. s ) goto 20 s = t l = j goto 20 10 e = dmax1(e,dabs(t)) 20 continue if ( s .gt. 10.*e ) goto 30 l = 0 goto 40 30 i = ia(l) call mod(u,ia,k,io,a,la,m,l,0,b) if ( i .lt. 0 ) l = -l lm = lm + 1 40 do 50 j = 1,n p(j) = 0. if ( ia(j) .eq. 0 ) p(j) = g(j) 50 continue e = dmax1(e,s) return 60 do 70 i = 1,m 70 c(i) = 0. do 90 j = 1,n if ( ia(j) .ne. 0 ) goto 90 t = g(j) do 80 i = 1,m 80 c(i) = c(i) - t*a(i,j) 90 continue call sol(b,u,m,c) s = 0. e = s do 120 j = 1,n t = g(j) do 100 i = 1,m 100 t = t + b(i)*a(i,j) if ( ia(j) .eq. 0 ) goto 110 p(j) = 0. if ( ia(j) .lt. 0 ) t = -t if ( t .le. s ) goto 120 l = j s = t goto 120 110 p(j) = t e = dmax1(e,dabs(t)) 120 continue if ( s .gt. 10.*e ) goto 130 l = 0 goto 190 130 t = g(l) do 140 i = 1,m 140 c(i) = c(i) - t*a(i,l) i = ia(l) call mod(u,ia,k,io,a,la,m,l,0,b) if ( i .lt. 0 ) l = -l if ( io .gt. 0 ) return lm = lm + 1 call sol(b,u,m,c) 150 do 180 j = 1,n if ( ia(j) .eq. 0 ) goto 160 p(j) = 0. goto 180 160 t = g(j) do 170 i = 1,m 170 t = t + b(i)*a(i,j) p(j) = t 180 continue 190 e = dmax1(e,s) return end subroutine pro(p,u,ia,a,la,m,n,b) integer ia(1),i,j,la,m,n real*8 a(la,1),b(1),p(1),u(1),t if ( m .gt. 0 ) goto 20 do 10 j = 1,n if ( ia(j) .ne. 0 ) p(j) = 0. 10 continue return 20 do 30 i = 1,m 30 b(i) = 0. do 50 j = 1,n if ( ia(j) .ne. 0 ) goto 50 t = p(j) do 40 i = 1,m 40 b(i) = b(i) - t*a(i,j) 50 continue call sol(b,u,m,b) do 80 j = 1,n t = 0. if ( ia(j) .ne. 0 ) goto 70 t = p(j) do 60 i = 1,m 60 t = t + b(i)*a(i,j) 70 p(j) = t 80 continue return end subroutine prp(p,e,u,v,ia,k,lm,io,g,a,la,l,nl,n,b,c) integer ia(1),h,i,ib,io,j,k,l,la,lm,nl,n real*8 a(la,1),b(1),c(1),g(1),p(1),u(1),v(1),e,t call pre(p,e,u,ia,k,ib,lm,io,g,a,la,l,n,b,c) h = iabs(ib) if ( ib .eq. 0 ) goto 10 call mob(v,io,u,ia,k,a,la,n,l,nl,ib,5,b,c) if ( io .gt. 0 ) return 10 do 20 i = 1,nl 20 c(i) = 0. do 40 j = 1,n t = p(j) if ( t .eq. 0. ) goto 40 do 30 i = 1,nl 30 c(i) = c(i) + t*a(i+l,j) 40 continue call sol(c,v,nl,c) if ( l .eq. 0 ) goto 60 do 50 i = 1,l 50 b(i) = 0. 60 do 90 j = 1,n if ( ia(j) .ne. 0 ) goto 90 t = 0. do 70 i = 1,nl 70 t = t + a(i+l,j)*c(i) p(j) = p(j) - t if ( l .eq. 0 ) goto 90 do 80 i = 1,l 80 b(i) = b(i) + t*a(i,j) 90 continue if ( l .eq. 0 ) goto 120 call sol(b,u,l,b) do 110 j = 1,n if ( ia(j) .ne. 0 ) goto 110 t = 0. do 100 i = 1,l 100 t = t + b(i)*a(i,j) p(j) = p(j) + t 110 continue 120 if ( h .eq. 0 ) return if ( ib .lt. 0 ) goto 130 if ( p(h) .ge. 0. ) return goto 140 130 if ( p(h) .le. 0. ) return 140 lm = lm - 1 call mob(v,io,u,ia,k,a,la,n,l,nl,ib,4,b,c) if ( io .gt. 0 ) return call mod(u,ia,k,io,a,la,l,ib,1,b) h = 0 if ( l .gt. 0 ) goto 160 do 150 j = 1,n p(j) = 0. 150 if ( ia(j) .eq. 0 ) p(j) = g(j) goto 10 160 do 170 i = 1,l 170 b(i) = 0. do 190 j = 1,n if ( ia(j) .ne. 0 ) goto 190 t = g(j) do 180 i = 1,l 180 b(i) = b(i) - t*a(i,j) 190 continue call sol(b,u,l,b) do 220 j = 1,n t = 0. if ( ia(j) .ne. 0 ) goto 210 t = g(j) do 200 i = 1,l 200 t = t + b(i)*a(i,j) 210 p(j) = t 220 continue goto 10 end subroutine shk(a,m,n) real*8 a(1) integer g,h,i,j,k,l,m,n if ( m .eq. 0 ) return g = n - m j = 0 l = 0 h = m 10 k = l + 1 l = l + h h = h - 1 do 20 i = k,l 20 a(i) = a(i+j) j = j + g if ( k .lt. l ) goto 10 return end subroutine sol(x,a,n,b) real*8 a(1),b(1),x(1),t integer i,j,k,l,n do 10 i = 1,n 10 x(i) = b(i) l = 0 k = 1 c ----------------------------- c |*** forward elimination ***| c ----------------------------- 20 if ( k .eq. n ) goto 40 t = x(k)/a(k+l) x(k) = t j = l l = l + n - k k = k + 1 if ( t .eq. 0. ) goto 20 do 30 i = k,n 30 x(i) = x(i) - t*a(i+j) goto 20 c ----------------------------------- c |*** back substitution by rows ***| c ----------------------------------- 40 x(n) = x(n)/a(k+l)**2 50 if ( k .eq. 1 ) return j = k k = k - 1 l = l + k - n t = x(k) do 60 i = j,n 60 t = t - x(i)*a(i+l) x(k) = t/a(k+l) goto 50 end subroutine xpn(a,m,n) real*8 a(1) integer g,h,i,j,k,l,m,n if ( m .eq. 0 ) return g = (m*m+m)/2 h = g l = n - m j = l*m k = 0 20 if ( g .le. 1 ) return j = j - l do 30 i = g,h 30 a(i+j) = a(i) h = g - 1 k = k + 1 g = h - k goto 20 end