{"id":272,"date":"2021-03-15T12:46:05","date_gmt":"2021-03-15T16:46:05","guid":{"rendered":"https:\/\/people.clas.ufl.edu\/hager\/?page_id=272"},"modified":"2026-03-19T08:18:25","modified_gmt":"2026-03-19T12:18:25","slug":"optpack","status":"publish","type":"page","link":"https:\/\/people.clas.ufl.edu\/hager\/optpack\/","title":{"rendered":"OPTPACK"},"content":{"rendered":"\r\n
\r\n\t
\r\n\t\t
\r\n\t\t\t
\r\n\t\t\t\t\n

OPTPACK<\/h1>\n\n\n\n

c OPTPACK was developed to test an algorithm presented in the paper
\nc “Analysis and implementation of a dual algorithm for constrained
\nc optimization,” Journal of Optimization Theory and Applications, 79 (1993),
\nc pp. 427-462. I have tested the software using the test problems
\nc described in that paper as well other other research type problems, however,
\nc I make no guarantees concerning the code.
\nc
\nc There are comment statements at the start of the 3 main routines
\nc (bmin for bound constraints, lmin for bound + linear inequality constraints,
\nc and min for bound + linear inequality + nonlinear constraints).
\nc I can provide pointers, but the user is pretty much on his own.
\nc
\nc Bill Hager
\nc
\nc ——————————————————————-
\nc \/ \/ \/ ____\/ William W. Hager
\nc \/ \/ \/ \/ 358 Little Hall
\nc \/ \/ \/ ____\/ Department of Mathematics
\nc \/__\/ \/ \/ UNIVERSITY OF FLORIDA
\nc _______\/ __\/ Gainesville FL 32611-8105
\nc
\nc Phone : (352) 392-0281 x 244 E-mail : hager@math.ufl.edu
\nc FAX : (352) 392-6254
\nc http:\/\/www.math.ufl.edu\/~hager
\nc ——————————————————————-<\/p>\n\n\n\n\n\n

C ________________________________________________________
\nC | |
\nC |SOLVE AN OPTIMIZATION PROBLEM WITH NONLINEAR CONSTRAINTS|
\nC | (ALL REAL VARIABLE ARE DOUBLE PRECISION) |
\nC | |
\nC | INPUT: |
\nC | |
\nC | X –STARTING GUESS (LENGTH AT LEAST N) |
\nC | |
\nC | ST –INITIAL STEP SIZE (WHEN ST = 0., PROGRAM |
\nC | SELECTS INITIAL STEP |
\nC | |
\nC | IN –DESCRIBES THE STARTING POINT |
\nC | =0, INFEASIBLE STARTING POINT, ALL |
\nC | CONSTRAINTS NONBINDING INITIALLY |
\nC | =1, INFEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS ARE GIVEN BY VALUES OF IA |
\nC | =2, INFEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS DETERMINED BY X |
\nC | =3, INFEASIBLE STARTING POINT, USE K, IA, |
\nC | AND U AS GIVEN |
\nC | =4, FEASIBLE STARTING POINT, ALL |
\nC | CONSTRAINTS NONBINDING INITIALLY |
\nC | =5, FEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS ARE GIVEN BY VALUES OF IA |
\nC | =6, FEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS DETERMINED BY X |
\nC | =7, FEASIBLE STARTING POINT, USE K, IA, |
\nC | AND U AS GIVEN |
\nC | |
\nC | TL –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | LM –MAXIMUM NUMBER OF ITERATIONS |
\nC | |
\nC | PNI –INPUT PENALTY PARAMETER (TRY PNI= 0 FIRST,|
\nC | IF A POSITIVE PENALTY IS REQUIRED, PROGRAM|
\nC | STOPS, CALL ERR(IO) TO GET ERROR MESSAGE) |
\nC | |
\nC | A –COEFFICIENT MATRIX FOR LINEAR CONSTRAINTS |
\nC | (AT LEAST L+NL ROWS) |
\nC | |
\nC | LA –LEADING DIMENSION OF A ARRAY |
\nC | |
\nC | LG –LEADING DIMENSION OF DG ARRAY |
\nC | |
\nC | L –NUMBER OF LINEAR EQUALITY CONSTRAINTS |
\nC | |
\nC | NL –NUMBER OF NONLINEAR EQUALITY CONSTRAINTS |
\nC | |
\nC | N –NUMBER OF UNKNOWNS |
\nC | |
\nC | B –RIGHT SIDE VECTOR FOR LINEAR CONSTRAINTS |
\nC | |
\nC | BL –LOWER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | BU –UPPER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | CT –CUTOFF FOR DECIDING IF A CONSTRAINT IS |
\nC | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE|
\nC | OF NORM OF TRUE SOLUTION) |
\nC | |
\nC | VLF –NAME OF ROUTINE TO EVALUATE COST FUNTION |
\nC | |
\nC | VLG –NAME OF ROUTINE TO EVALUATE CONSTRAINTS |
\nC | |
\nC | GRF –NAME OF ROUTINE TO EVALUATE COST GRADIENT |
\nC | |
\nC | GRG –NAME OF ROUTINE TO EVALUATE CONSTRAINT |
\nC | GRADIENT |
\nC | |
\nC | W –WORK ARRAY (LENGTH AT LEAST 5M+3N IF |
\nC | PENALTY NOT INCREASED, LENGTH AT LEAST |
\nC | .5NL(NL+5)+L+3N OTHERWISE |
\nC | |
\nC | |
\nC | OUTPUT: |
\nC | |
\nC | X –SOLUTION |
\nC | |
\nC | LE –MULTIPLIERS FOR EQUALITY CONSTRAINTS |
\nC | (LENGTH AT LEAST L+NL) |
\nC | |
\nC | LI –MULTIPLIERS FOR INEQUALITY CONSTRAINTS |
\nC | (LENGTH AT LEAST N) |
\nC | |
\nC | U –CHOLESKY FACTORIZATION ASSOCIATED WITH |
\nC | CONSTRAINT PROJECTION (LENGTH AT LEAST |
\nC | .5M(M+1) WHERE M = L + NL) |
\nC | |
\nC | IA –= 0 IF CONSTRAINT IS NONBINDING |
\nC | = 1 IF CONSTRAINT IS AT UPPER BOUND |
\nC | =-1 IF CONSTRAINT IS AT LOWER BOUND |
\nC | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) |
\nC | |
\nC | K –NUMBER OF FREE VARIABLES |
\nC | |
\nC | F –VALUE OF COST AT OPTIMAL POINT |
\nC | |
\nC | G –VALUE OF CONSTRAINTS AT OPTIMAL POINT |
\nC | (LENGTH AT LEAST NL) |
\nC | |
\nC | DF –DERIVATIVE OF COST AT OPTIMAL POINT |
\nC | (LENGTH AT LEAST N) |
\nC | |
\nC | DG –DERIVATIVE OF CONSTRAINTS AT OPTIMAL POINT |
\nC | |
\nC | E –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | IT –NUMBER OF ITERATIONS |
\nC | |
\nC | ST –FINAL STEP SIZE |
\nC | |
\nC | PN –FINAL VALUE OF PENALTY |
\nC | |
\nC | IO –ERROR FLAG (CALL ERR(IO) TO SEE IF AN |
\nC | ERROR WAS DETECTED) |
\nC |________________________________________________________|
\nC
\n subroutine min(x,le,li,u,ia,k,f,g,df,dg,e,it,st,pn,io,
\n 1 in,tl,lm,pni,a,la,lg,l,nl,n,b,bl,bu,ct,vlf,vlg,grf,grg,w)
\n integer ia(1),i,in,io,ip,is,it,j,k,l,la,lg,lm
\n integer m,m1,m2,m3,m4,m5,m6,m7,n,nl,n1,n2,n3,n4
\n real*8 a(la,1),b(1),bl(1),bu(1),df(1),dg(lg,1),g(1),le(1),li(1)
\n real*8 u(1),w(1),x(1),big,ct,e,f,kt,pn,pni,r,s,st,t,tl,vlf
\n external grf,grg,vlf,vlg
\n m = l + nl
\n if ( m .le. n ) goto 10
\n io = 8
\n return
\n10 if ( m .gt. 0 ) goto 20
\n call bmin(x,ia,k,f,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vlf,grf,w)
\n return
\n20 if ( nl .gt. 0 ) goto 25
\n call lmin(x,le,li,u,ia,k,f,df,e,it,st,io,
\n & in,tl,lm,a,la,l,n,b,bl,bu,ct,vlf,grf,w)
\n return
\n25 m1 = m + 1
\n m2 = m + m1
\n m3 = m + m2
\n m4 = m + m3
\n m5 = n + m4
\n m6 = m + m5 + n
\n m7 = m6 – 1
\n n1 = .5*(nl*nl+nl) + 1
\n n2 = n1 + nl
\n n3 = n2 + l
\n n4 = n3 + nl
\n pn = pni
\nc
\nc largest floating point number
\nc
\n big = 1.d50
\n if ( in .lt. 4 ) goto 110
\n if ( in .eq. 7 ) goto 130
\n k = n – l
\n if ( in .eq. 6 ) goto 60
\n if ( in .eq. 5 ) goto 40
\n do 30 i = 1,n
\n30 ia(i) = 0
\n goto 90
\n40 do 50 i = 1,n
\n if ( ia(i) .eq. 0 ) goto 50
\n k = k – 1
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n50 continue
\n goto 90
\n60 do 80 i = 1,n
\n ia(i) = 0
\n if ( x(i) .le. bl(i) ) goto 70
\n if ( x(i) .lt. bu(i) ) goto 80
\n k = k – 1
\n x(i) = bu(i)
\n ia(i) = 1
\n goto 80
\n70 k = k – 1
\n x(i) = bl(i)
\n ia(i) = -1
\n80 continue
\n90 if ( k .ge. 0 ) goto 100
\n io = 1
\n return
\n100 call mat(u,ia,a,la,l,n,0)
\n call fac(u,0,l,io)
\n if ( io .eq. 0 ) goto 130
\n io = 1
\n return
\n110 call feas(li,u,ia,k,io,in,x,a,la,l,n,b,bl,bu,w)
\n do 120 i = 1,n
\n120 x(i) = li(i)
\n if ( io .gt. 0 ) return
\n130 f = vlf(x)
\n call vlg(g,x)
\n call grf(df,x)
\n call grg(dg,x)
\n it = 0
\n io = 0
\n call erg(t,g,nl)
\n140 is = 0
\n ip = 0
\n if ( k .lt. nl ) ip = 1
\n if ( k .ge. nl ) k = k – nl
\n call xpn(u,l,m)
\n call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m)
\n call err(io)
\n if ( io .gt. 0 ) goto 420
\n call mul(le,li,s,u,ia,a,la,m,n,df,l,nl,w(m2))
\n e = s + t
\n if ( e .le. tl ) return
\n kt = .95*e
\n goto 190
\n150 call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m)
\n call err(io)
\n if ( io .gt. 0 ) goto 420
\n call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2))
\n e = s + t
\n if ( e .gt. tl ) goto 180
\n160 do 170 i = 1,m
\n170 le(i) = w(i+m)
\n return
\n180 if ( t .le. s ) goto 320
\n190 is = is + 1
\n if ( is .lt. 40 ) goto 200
\n io = 6
\n return
\n200 if ( l .eq. 0 ) goto 230
\n do 210 i = 1,l
\n210 w(i) = b(i)
\n do 220 j = 1,n
\n s = x(j)
\n do 220 i = 1,l
\n220 w(i) = w(i) – s*a(i,j)
\n230 do 240 i = 1,n
\n240 w(i+m7) = ia(i)
\n if ( ip .eq. 1 ) goto 250
\n call nf1(li,u,ia,k,io,x,a,la,m,n,w,bl,bu,w(m1),w(m2),w(m2+n))
\n call err(io)
\n if ( io .gt. 0 ) goto 420
\n if ( io .eq. 0 ) goto 290
\n io = 0
\n goto 300
\n250 do 260 i = 1,m
\n w(i+m) = 0.
\n260 if ( i .gt. l ) w(i+m) = 1.
\n call nf2(li,u,ia,k,io,x,a,la,m,n,w,w(m1),bl,bu,w(m2),w(m3),w(m5))
\n if ( io .eq. 0 ) goto 280
\n do 270 i = 1,n
\n j = w(i+m7)
\n call mod(u,ia,k,io,a,la,m,i,j,w)
\n if ( io .gt. 0 ) return
\n270 continue
\n goto 420
\n280 ip = 0
\n290 call lsq(li,u,ia,k,io,x,a,la,m,n,bl,bu,w,w(m1),w(m2))
\n call err(io)
\n if ( io .gt. 0 ) goto 420
\n300 do 310 i = 1,n
\n310 w(i) = x(i)
\n call arm(x,li,w,u,ia,k,io,t,w(m6),a,la,m,n,vlg,g,nl,tl)
\n f = vlf(x)
\n call grf(df,x)
\n call grg(dg,x)
\n if ( io .eq. 0 ) goto 150
\n return
\n320 r = .95*s
\n is = -1
\n do 330 i = 1,m
\n330 le(i) = w(i+m)
\n340 is = is + 1
\n if ( is .le. 3 ) goto 350
\n if ( e .le. kt ) goto 390
\n goto 420
\n350 j = k
\n call cgn(x,e,i,st,io,.1*tl,j,lm-it,n,is,vlf,vlg,grf,grg,
\n 1 nl,pn,a,la,m,bl,bu,ct,u,ia,k,le(l+1),f,g,df,dg,lg,
\n 2 w,w(m1),w(m2),w(m3),w(m4))
\n call err(io)
\n it = it + i
\n if ( io .gt. 0 ) return
\n call erg(t,g,nl)
\n call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m)
\n call err(io)
\n if ( io .gt. 0 ) return
\n call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2))
\n e = s + t
\n if ( e .le. tl ) goto 160
\n if ( it .lt. lm ) goto 360
\n io = 9
\n goto 160
\n360 if ( s .gt. r ) goto 380
\n r = .95*s
\n is = 0
\n do 370 i = 1,m
\n370 le(i) = w(i+m)
\n if ( 4.*t .lt. s ) goto 340
\n if ( e .gt. kt ) goto 410
\n kt = .95*e
\n goto 190
\n380 if ( t .gt. .5*kt ) goto 420
\n if ( 4.*t .lt. s ) goto 340
\n if ( e .gt. kt ) goto 410
\n390 kt = .95*e
\n is = 0
\n do 400 i = 1,m
\n400 le(i) = w(i+m)
\n goto 190
\n410 if ( t .le. .25*kt ) goto 340
\n420 call shk(u,l,m)
\n is = 0
\n if ( ip .eq. 0 ) k = k + nl
\n if ( it .gt. 0 ) pn = 5.*pn
\n if ( pn .gt. 0. ) goto 430
\n io = 7
\n return
\n430 j = k
\n call cgp(x,e,i,st,io,.1*tl,j,lm-it,n,vlf,vlg,grf,grg,
\n 1 l,nl,pn,a,la,bl,bu,ct,u,ia,k,le(l+1),f,g,df,dg,lg,
\n 2 w,w(n1),w(n2),w(n3),w(n4))
\n call err(io)
\n it = it + i
\n if ( io .gt. 0 ) return
\n if ( it .lt. lm ) goto 440
\n io = 9
\n return
\n440 call erg(t,g,nl)
\n call xpn(u,l,m)
\n ip = 0
\n if ( k .lt. nl ) ip = 1
\n if ( k .ge. nl ) k = k – nl
\n call fat(u,io,ip,ia,a,la,l,n,g,w(l+1),w(m2),dg,lg,nl,m)
\n call err(io)
\n if ( io .gt. 0 ) return
\n call mul(w(m1),li,s,u,ia,a,la,m,n,df,l,nl,w(m2))
\n r = e
\n e = s + t
\n if ( e .le. tl ) goto 160
\n if ( e .le. .6*kt ) goto 450
\n if ( r .le. t ) goto 460
\n call shk(u,l,m)
\n if ( ip .eq. 0 ) k = k + nl
\n goto 430
\n450 pn = pn*.4
\n460 do 470 i = 1,m
\n470 le(i) = w(i+m)
\n kt = e
\n goto 190
\n end
\nC ________________________________________________________
\nC | |
\nC | SOLVE AN OPTIMIZATION PROBLEM WITH UPPER AND LOWER |
\nC | BOUND CONSTRAINTS |
\nC | (ALL REAL VARIABLE ARE DOUBLE PRECISION) |
\nC | |
\nC | INPUT: |
\nC | |
\nC | X –STARTING GUESS (LENGTH AT LEAST N) |
\nC | |
\nC | ST –INITIAL STEP SIZE (WHEN ST = 0., PROGRAM |
\nC | SELECTS INITIAL STEP |
\nC | |
\nC | IN –DESCRIBES THE STARTING POINT |
\nC | =0, ALL CONSTRAINTS NONBINDING INITIALLY |
\nC | =1, BINDING CONSTRAINTS GIVEN BY IA |
\nC | =2, BINDING CONSTRAINTS DETERMINED BY X |
\nC | =3, USE K, IA, AND W AS GIVEN |
\nC | |
\nC | TL –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | LM –MAXIMUM NUMBER OF ITERATIONS |
\nC | |
\nC | N –NUMBER OF UNKNOWNS |
\nC | |
\nC | BL –LOWER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | BU –UPPER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | CT –CUTOFF FOR DECIDING IF A CONSTRAINT IS |
\nC | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE|
\nC | OF NORM OF TRUE SOLUTION) |
\nC | |
\nC | VL –NAME OF ROUTINE TO EVALUATE COST FUNTION |
\nC | |
\nC | GR –NAME OF ROUTINE TO EVALUATE COST GRADIENT |
\nC | |
\nC | W –WORK ARRAY (LENGTH AT LEAST 3N) |
\nC | |
\nC | OUTPUT: |
\nC | |
\nC | X –SOLUTION |
\nC | |
\nC | IA –= 0 IF CONSTRAINT IS NONBINDING |
\nC | = 1 IF CONSTRAINT IS AT UPPER BOUND |
\nC | =-1 IF CONSTRAINT IS AT LOWER BOUND |
\nC | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) |
\nC | |
\nC | K –NUMBER OF FREE VARIABLES |
\nC | |
\nC | VF –VALUE OF COST AT OPTIMAL POINT |
\nC | |
\nC | DF –DERIVATIVE OF COST AT OPTIMAL POINT |
\nC | (LENGTH AT LEAST N) |
\nC | |
\nC | E –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | IT –NUMBER OF ITERATIONS |
\nC | |
\nC | ST –FINAL STEP SIZE |
\nC | |
\nC | IO –ERROR FLAG (CALL ERR(IO) TO SEE IF AN |
\nC | ERROR WAS DETECTED) |
\nC |________________________________________________________|
\nC
\n subroutine bmin
\n 1 (x,ia,k,vf,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vl,gr,w)
\n real*8 a(1),b(1),bl(1),bu(1),c(1),df(1),x(1),u(1),w(1)
\n real*8 ct,e,st,t,tl,vf,vl
\n integer ia(1),i,in,io,it,j,k,lm,n
\n external gr,vl
\n it = 0
\n io = 0
\n if ( in .eq. 3 ) goto 80
\n k = n
\n do 10 i = 1,n
\n if ( bl(i) .lt. bu(i) ) goto 10
\n t = bl(i)
\n bl(i) = bu(i)
\n bu(i) = t
\n10 continue
\n if ( in .eq. 2 ) goto 50
\n if ( in .eq. 1 ) goto 30
\n do 20 i = 1,n
\n20 ia(i) = 0
\n goto 80
\n30 do 40 i = 1,n
\n if ( ia(i) .eq. 0 ) goto 40
\n k = k – 1
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n40 continue
\n goto 80
\n50 do 70 i = 1,n
\n ia(i) = 0
\n if ( x(i) .lt. bu(i) ) goto 60
\n x(i) = bu(i)
\n ia(i) = 1
\n k = k – 1
\n goto 70
\n60 if ( x(i) .gt. bl(i) ) goto 70
\n x(i) = bl(i)
\n ia(i) = -1
\n k = k – 1
\n70 continue
\n80 call gr(df,x)
\n vf = vl(x)
\n90 j = k
\n call cgl(x,e,i,st,io,tl,j,lm-it,n,vl,gr,
\n 1 a,1,0,bl,bu,ct,u,ia,k,vf,df,b,c,w)
\n it = it + i
\n if ( io .gt. 0 ) return
\n if ( e .le. tl ) return
\n if ( it .lt. lm ) goto 90
\n io = 9
\n return
\n end
\nC ________________________________________________________
\nC | |
\nC | SOLVE AN OPTIMIZATION PROBLEM WITH LINEAR CONSTRAINTS |
\nC | (ALL REAL VARIABLE ARE DOUBLE PRECISION) |
\nC | |
\nC | INPUT: |
\nC | |
\nC | X –STARTING GUESS (LENGTH AT LEAST N) |
\nC | |
\nC | ST –INITIAL STEP SIZE (WHEN ST = 0., PROGRAM |
\nC | SELECTS INITIAL STEP |
\nC | |
\nC | IN –DESCRIBES THE STARTING POINT |
\nC | =0, INFEASIBLE STARTING POINT, ALL |
\nC | CONSTRAINTS NONBINDING INITIALLY |
\nC | =1, INFEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS ARE GIVEN BY VALUES OF IA |
\nC | =2, INFEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS DETERMINED BY X |
\nC | =3, INFEASIBLE STARTING POINT, USE K, IA, |
\nC | AND W AS GIVEN |
\nC | =4, FEASIBLE STARTING POINT, ALL |
\nC | CONSTRAINTS NONBINDING INITIALLY |
\nC | =5, FEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS ARE GIVEN BY VALUES OF IA |
\nC | =6, FEASIBLE STARTING POINT, BINDING |
\nC | CONSTRAINTS DETERMINED BY X |
\nC | =7, FEASIBLE STARTING POINT, USE K, IA, |
\nC | AND U AS GIVEN |
\nC | |
\nC | TL –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | LM –MAXIMUM NUMBER OF ITERATIONS |
\nC | |
\nC | A –COEFFICIENT MATRIX FOR LINEAR CONSTRAINTS |
\nC | (AT LEAST M ROWS) |
\nC | |
\nC | LA –LEADING DIMENSION OF A ARRAY |
\nC | |
\nC | M –NUMBER OF EQUALITY CONSTRAINTS |
\nC | |
\nC | N –NUMBER OF UNKNOWNS |
\nC | |
\nC | B –RIGHT SIDE VECTOR FOR LINEAR CONSTRAINTS |
\nC | (LENGTH AT LEAST M) |
\nC | |
\nC | BL –LOWER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | BU –UPPER BOUNDS FOR COMPONENTS OF X |
\nC | |
\nC | CT –CUTOFF FOR DECIDING IF A CONSTRAINT IS |
\nC | BINDING (FOR EXAMPLE, .0001 TIMES ESTIMATE|
\nC | OF NORM OF TRUE SOLUTION) |
\nC | |
\nC | VL –NAME OF ROUTINE TO EVALUATE COST FUNTION |
\nC | |
\nC | GR –NAME OF ROUTINE TO EVALUATE COST GRADIENT |
\nC | |
\nC | W –WORK ARRAY(LENGTH AT LEAST MAX(4M+2N,M+3N)|
\nC | |
\nC | OUTPUT: |
\nC | |
\nC | X –SOLUTION |
\nC | |
\nC | LE –MULTIPLIERS FOR EQUALITY CONSTRAINTS |
\nC | |
\nC | LI –MULTIPLIERS FOR INEQUALITY CONSTRAINTS |
\nC | |
\nC | U –CHOLESKY FACTORIZATION ASSOCIATED WITH |
\nC | CONSTRAINT PROJECTION |
\nC | (LENGTH AT LEAST .5M(M+1)) |
\nC | |
\nC | IA –= 0 IF CONSTRAINT IS NONBINDING |
\nC | = 1 IF CONSTRAINT IS AT UPPER BOUND |
\nC | =-1 IF CONSTRAINT IS AT LOWER BOUND |
\nC | (INTEGER ARRAY WITH AT LEAST N ELEMENTS) |
\nC | |
\nC | K –NUMBER OF FREE VARIABLES |
\nC | |
\nC | VF –VALUE OF COST AT OPTIMAL POINT |
\nC | |
\nC | DF –DERIVATIVE OF COST AT OPTIMAL POINT |
\nC | (LENGTH AT LEAST N) |
\nC | |
\nC | E –ERROR TOLERANCE (INFINITY NORM OF KUHN- |
\nC | TUCKER ERROR + CONSTRAINT VIOLATION) |
\nC | |
\nC | IT –NUMBER OF ITERATIONS |
\nC | |
\nC | ST –FINAL STEP SIZE |
\nC | |
\nC | IO –ERROR FLAG (CALL ERR(IO) TO SEE IF AN |
\nC | ERROR WAS DETECTED) |
\nC |________________________________________________________|
\nC
\n subroutine lmin(x,le,li,u,ia,k,vf,df,e,it,st,io,
\n 1 in,tl,lm,a,la,m,n,b,bl,bu,ct,vl,gr,w)
\n real*8 a(1),b(1),bl(1),bu(1),df(1),le(1),li(1),u(1),w(1),x(1),z(1)
\n real*8 ct,e,tl,st,vf,vl
\n integer ia(1),i,in,io,it,k,la,lm,m,m1,n
\n external gr,vl
\n m1 = m + 1
\n if ( m .le. n ) goto 10
\n io = 8
\n return
\n10 if ( m .gt. 0 ) goto 20
\n call bmin(x,ia,k,vf,df,e,it,st,io,in,tl,lm,n,bl,bu,ct,vl,gr,w)
\n io = 0
\n return
\n20 if ( in .lt. 4 ) goto 100
\n if ( in .eq. 7 ) goto 120
\n k = n – m
\n if ( in .eq. 6 ) goto 60
\n if ( in .eq. 5 ) goto 40
\n do 30 i = 1,n
\n30 ia(i) = 0
\n goto 90
\n40 do 50 i = 1,n
\n if ( ia(i) .eq. 0 ) goto 50
\n k = k – 1
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n50 continue
\n if ( k .ge. 0 ) goto 90
\n io = 8
\n return
\n60 do 80 i = 1,n
\n ia(i) = 0
\n if ( x(i) .le. bl(i) ) goto 70
\n if ( x(i) .lt. bu(i) ) goto 80
\n k = k – 1
\n x(i) = bu(i)
\n ia(i) = 1
\n goto 80
\n70 k = k – 1
\n x(i) = bl(i)
\n ia(i) = -1
\n80 continue
\n if ( k .ge. 0 ) goto 90
\n io = 8
\n return
\n90 call mat(u,ia,a,la,m,n,0)
\n call fac(u,0,m,io)
\n if ( io .gt. 0 ) return
\n goto 120
\n100 do 110 i = 1,n
\n110 li(i) = x(i)
\n call feas(x,u,ia,k,io,in,li,a,la,m,n,b,bl,bu,w)
\n if ( io .gt. 0 ) return
\n120 call gr(df,x)
\n vf = vl(x)
\n it = 0
\n io = 0
\n130 j = k
\nc
\nc Projected cg step
\nc
\n call cgl(x,e,i,st,io,tl,j,lm-it,n,vl,gr,
\n 1 a,la,m,bl,bu,ct,u,ia,k,vf,df,w,le,w(m1))
\n it = it + i
\n if ( io .eq. 1 ) return
\n if ( io .gt. 1 ) goto 140
\n if ( it .ge. lm ) goto 140
\n if ( e .gt. tl ) goto 130
\n140 call mul(le,li,e,u,ia,a,la,m,n,df,m,0,z)
\n return
\n end
\n subroutine arm(x,y,z,u,iy,k,io,t,iz,a,la,m,n,vlg,g,nl,tl)
\n real*8 a(1),g(1),iz(1),u(1),x(1),y(1),z(1),r,s,t,tl
\n integer iy(1),i,io,j,k,m,n,nl
\n io = 0
\n call vlg(g,y)
\n call erg(s,g,nl)
\n if ( s .le. tl ) goto 10
\n if ( s .gt. .5*t ) goto 30
\n10 t = s
\n do 20 i = 1,n
\n20 x(i) = y(i)
\n return
\n30 do 40 i = 1,n
\n if ( iy(i) .eq. 0 ) goto 40
\n if ( iy(i) .eq. iz(i) ) goto 40
\n call mod(u,iy,k,io,a,la,m,i,0,x)
\n if ( io .eq. 1 ) return
\n40 continue
\n r = .5
\n j = 0
\n50 j = j + 1
\n if ( j .gt. 20 ) goto 70
\n do 60 i = 1,n
\n60 x(i) = r*y(i) + (1-r)*z(i)
\n call vlg(g,x)
\n call erg(s,g,nl)
\n r = .5*r
\n if ( s .gt. t*(1-r) ) goto 50
\n t = s
\n return
\n70 io = 6
\n t = s
\n return
\n end
\n subroutine cgl(x,e,it,st,io,tl,lm1,lm2,n,vl,gr,
\n 1 aa,la,m,bl,bu,ct,u,ia,k,vf,df,bb,cc,h)
\n integer ia(1),i,ib,io,it,j,k,lm1,lm2,m,n,na,nb,nc,nd
\n real*8 h(n,1),x(1),y(50),z(50)
\n real*8 aa(1),bb(1),bl(1),bu(1),cc(1),df(1),u(1),vf
\n real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0
\n real*8 da,db,e,f,f0,f1,fa,fb,fc,g,l3,p,q,r,s,st,tl,v,w
\n real*8 fv,fd
\n external gr,vl
\n data a1\/.1d0\/,a2\/.9d0\/,a3\/5.d0\/,a4\/.2d0\/,a5\/10.d0\/,a6\/.9d0\/
\n data a7\/.3d0\/
\nc
\nc set big = largest floating point number
\nc
\n big = 1.d50
\n a8 = a3 + .01d0
\n it = 0
\n io = 0
\n f = vf
\n do 10 i = 1,n
\n10 h(i,3) = df(i)
\n l3 = 1.\/dlog(a3)
\n call pre(h(1,2),e,u,ia,k,i,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n if ( e .le. tl ) goto 690
\n if ( io .gt. 0 ) goto 690
\n if ( lm1 .eq. 0 ) return
\n a = st
\n if ( a .gt. 0. ) goto 30
\n do 20 i = 1,n
\n20 if ( dabs(x(i)) .gt. a ) a = dabs(x(i))
\n a = .01*a\/e
\n if ( a .eq. 0. ) a = 1.
\n30 g = 0.
\n do 40 i = 1,n
\n h(i,1) = -h(i,2)
\n40 g = g + h(i,2)*h(i,3)
\n if ( g .lt. 0. ) goto 650
\n d = -g
\n50 call cut(s,ib,x,h,ia,bl,bu,ct,n,big)
\n na = 0
\n nb = 0
\n nc = 0
\n nd = 0
\n if ( s .le. 0. ) goto 760
\n if ( a .gt. s ) a = s
\n fa = fv(a,x,h,n,vl)
\n c0 = a
\n f0 = fa
\n ny = 2
\n y(1) = 0.
\n z(1) = f
\n y(2) = a
\n z(2) = fa
\n v = a1*d
\n w = a2*d
\n iq = 0
\n if ( fa .le. f ) goto 60
\n c = a
\n b = 0.
\n a = 0.
\n fc = fa
\n fb = f
\n fa = f
\n goto 70
\n60 c = 0.
\n b = 0.
\n fc = f
\n fb = f
\n iq = 1
\n70 q = (d+(f-f0)\/c0)\/c0
\n if ( q .lt. 0. ) goto 90
\n q = a
\n p = fa
\n80 nd = nd + 1
\n if ( nd .gt. 25 ) goto 640
\n if ( s .eq. q ) goto 710
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .lt. w*q ) goto 80
\n goto 270
\n90 q = .5*d\/q
\n if ( q .lt. s ) goto 110
\n if ( c0 .lt. s ) goto 100
\n f1 = f0
\n c1 = c0
\n q = s
\n goto 140
\n100 q = s
\n110 if ( q .lt. .01*c0 ) q = .01*c0
\n p = fv(q,x,h,n,vl)
\n if ( p .le. f0 ) goto 120
\n f1 = f0
\n c1 = c0
\n f0 = p
\n c0 = q
\n goto 130
\n120 f1 = p
\n c1 = q
\n130 call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n140 if ( a .eq. 0. ) goto 150
\n if ( fa-f .ge. v*a ) goto 170
\n if ( fa-f .lt. w*a ) goto 220
\n goto 290
\n150 q = c0
\n if ( c1 .lt. q ) q = c1
\n160 na = na + 1
\n if ( na .gt. 25 ) goto 660
\n q = a4*q
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .ge. v*q ) goto 160
\n goto 260
\n170 if ( c0 .gt. c1 ) goto 210
\n if ( f0-f .gt. v*c0 ) goto 190
\n if ( f0-f .ge. w*c0 ) goto 350
\n if ( c1 .le. a5*c0 ) goto 350
\n r = dlog(c1\/c0)
\n r = .999*dexp(-r\/idint(r*l3+.999))
\n q = c1
\n180 q = q*r
\n if ( q .lt. c0 ) goto 350
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n na = na + 1
\n if ( p-f .gt. v*q ) goto 180
\n goto 350
\n190 q = c0
\n200 na = na + 1
\n if ( na .gt. 25 ) goto 660
\n q = a4*q
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .ge. v*q ) goto 200
\n goto 260
\n210 q = a
\n goto 200
\n220 if ( c0 .lt. c1 ) goto 320
\n if ( f0-f .ge. v*c0 ) goto 240
\n if ( f0-f .ge. w*c0 ) goto 260
\n q = c0
\n p = f0
\n230 nd = nd + 1
\n if ( nd .gt. 25 ) goto 640
\n if ( s .eq. q ) goto 710
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .lt. w*q ) goto 230
\n goto 260
\n240 if ( c0 .le. a5*c1 ) goto 260
\n r = dlog(c0\/c1)
\n r = 1.001*dexp(r\/idint(r*l3+.999))
\n q = a
\n250 q = q*r
\n if ( q .gt. c0 ) goto 260
\n nd = nd + 1
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .lt. w*q ) goto 250
\n260 if ( iq .eq. 1 ) goto 350
\n270 if ( b .eq. 0. ) goto 290
\n if ( c .eq. 0. ) goto 280
\n v = c – a
\n w = a – b
\n r = 1.\/v
\n f0 = 1.\/w
\n p = fc – fa
\n q = fb – fa
\n e = p*r + q*f0
\n if ( dsign(e,c-b) .ne. e ) goto 350
\n if ( e .eq. 0. ) goto 350
\n q = (p*r)*w – (q*f0)*v
\n q = a – .5*q\/e
\n goto 300
\n280 r = 1.\/a
\n f0 = 1.\/b
\n p = r*(fa-f) – d
\n q = f0*(fb-f) – d
\n e = a – b
\n v = (r*p-f0*q)\/e
\n w = (a*q*f0-b*p*r)\/e
\n v = w*w-3.*v*d
\n if ( v .lt. 0. ) v = 0.
\n v = dsqrt(v)
\n if ( w+v .eq. 0. ) goto 350
\n q = -d\/(w+v)
\n if ( q .le. 0. ) goto 350
\n goto 300
\n290 if ( iq .eq. 1 ) goto 350
\n q = (d+(f-fa)\/a)\/a
\n if ( q .ge. 0. ) goto 350
\n q = .5*d\/q
\n300 if ( q .gt. s ) q = s
\n do 310 i = 1,ny
\n310 if ( y(i) .eq. q ) goto 350
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n goto 350
\n320 continue
\n if ( f0-f .gt. v*c0 ) goto 330
\n if ( f0-f .gt. w*c0 ) goto 350
\n330 q = a
\n p = fa
\n340 nd = nd + 1
\n if ( nd .gt. 25 ) goto 640
\n if ( s .eq. q ) goto 710
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .lt. w*q ) goto 340
\n goto 260
\n350 da = fd(a,x,h,n,gr)
\n if ( da .gt. a6*g ) goto 440
\n if ( s .eq. a ) goto 780
\n if ( da .ge. 0. ) goto 590
\n r = a
\n q = 0.
\n do 360 i = 1,ny
\n if ( y(i) .gt. a ) goto 400
\n if ( y(i) .le. q ) goto 360
\n if ( y(i) .eq. a ) goto 360
\n q = y(i)
\n360 continue
\n if ( a .le. a8*q ) goto 590
\n q = a
\n p = fa
\n370 nd = nd + 1
\n if ( nd .gt. 25 ) goto 640
\n if ( s .eq. q ) goto 710
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = fv(q,x,h,n,vl)
\n f1 = fa
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p .lt. f1 ) goto 370
\n if ( a .gt. r ) goto 390
\n do 380 i = 1,n
\n380 h(i,2) = x(i) + a*h(i,1)
\n goto 590
\n390 da = fd(a,x,h,n,gr)
\n if ( da .gt. a6*g ) goto 440
\n goto 590
\n400 q = y(i)
\n do 410 j = i,ny
\n if ( y(j) .le. a ) goto 410
\n if ( y(j) .lt. q ) q = y(j)
\n410 continue
\n if ( q .le. a5*a ) goto 590
\n f0 = dlog(q\/a)
\n f0 = 1.001*dexp(f0\/idint(f0*l3+.999))
\n v = a
\n420 v = v*f0
\n if ( v .ge. q ) goto 350
\n p = fv(v,x,h,n,vl)
\n f1 = fa
\n call ins(v,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p .lt. f1 ) goto 420
\n if ( a .gt. r ) goto 350
\n do 430 i = 1,n
\n430 h(i,2) = x(i) + a*h(i,1)
\n goto 590
\n440 b = 0.
\n j = 1
\n i = j
\n450 i = i + 1
\n if ( i .gt. ny ) goto 460
\n if ( y(i) .ge. a ) goto 450
\n if ( y(i) .lt. b ) goto 450
\n b = y(i)
\n j = i
\n goto 450
\n460 fb = z(j)
\n db = d
\n if ( b .ne. 0. ) db = fd(b,x,h,n,gr)
\n470 w = 2.*dabs(b-a)
\n call cub(c,a,b,fa,fb,da,db)
\n nc = 1
\n goto 510
\n480 w = .5*w
\n if ( w .lt. dabs(c0-c) ) goto 580
\n if ( c0 .lt. c ) goto 490
\n if ( d0 .ge. d ) goto 500
\n goto 580
\n490 if ( d0 .gt. d ) goto 580
\n500 call cub(c,c,c0,f,f0,d,d0)
\n nc = nc + 1
\n if ( nc .gt. 30 ) goto 630
\n510 r = dmax1(a,b)
\n s = dmin1(a,b)
\n if ( c .gt. r ) goto 520
\n if ( c .gt. s ) goto 530
\n c = s + (s-c)
\n s = .5*(a+b)
\n if ( c .gt. s ) c = s
\n goto 530
\n520 c = r – (c-r)
\n s = .5*(a+b)
\n if ( c .lt. s ) c = s
\n530 c0 = a
\n f0 = fa
\n d0 = da
\n d = fd(c,x,h,n,gr)
\n f = fv(c,x,h,n,vl)
\n if ( f .lt. fa ) goto 540
\n b = c
\n fb = f
\n db = d
\n goto 480
\n540 if ( c .lt. a ) goto 570
\n if ( d .lt. 0. ) goto 560
\n550 b = a
\n fb = fa
\n db = da
\n560 a = c
\n fa = f
\n da = d
\n if ( d .gt. a6*g ) goto 480
\n goto 590
\n570 if ( d .lt. 0. ) goto 550
\n goto 560
\n580 c = .5*(a+b)
\n nb = nb + 1
\n w = dabs(b-a)
\n goto 530
\n590 do 600 i = 1,n
\n600 x(i) = h(i,2)
\n i = lm1
\n call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n it = it + 1
\n f = fa
\n d = da
\n a = a7*a
\n if ( io .gt. 0 ) goto 690
\n if ( e .le. tl ) goto 690
\n if ( it .ge. lm2 ) goto 690
\n if ( it .ge. lm1 ) goto 690
\n if ( i .lt. lm1 ) goto 30
\n r = 0.
\n do 610 i = 1,n
\n610 r = r + h(i,2)*h(i,3)
\n if ( r .lt. 0. ) goto 650
\n s = r\/g
\n g = r
\n d = 0.
\n do 620 i = 1,n
\n h(i,1) = -h(i,2) + s*h(i,1)
\n620 d = d + h(i,1)*h(i,3)
\n goto 50
\n630 if ( d .lt. g ) goto 590
\n io = 4
\n return
\n640 io = 3
\n return
\n650 io = 4
\n return
\n660 q = q*a3**25
\n nd = 0
\n670 nd = nd + 1
\n if ( nd .gt. 25 ) goto 680
\n q = a3*q
\n if ( q .ge. s ) goto 680
\n p = fv(q,x,h,n,vl)
\n call ins(q,p,a,b,c,fa,fb,fc,ny,y,z)
\n if ( p-f .gt. v*q ) goto 670
\n goto 140
\n680 io = 5
\n return
\n690 st = a
\n vf = f
\n do 700 i = 1,n
\n700 df(i) = h(i,3)
\n return
\n710 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n do 720 i = 1,n
\n720 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n a = q*a7
\n730 call gr(h(1,3),x)
\n740 f = p
\n750 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n it = it + 1
\n if ( it .ge. lm2 ) goto 690
\n if ( it .ge. lm1 ) goto 690
\n if ( io .gt. 0 ) goto 690
\n if ( e .le. tl ) goto 690
\n goto 30
\n760 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n q = -s
\n p = fv(q,x,h,n,vl)
\n do 770 i = 1,n
\n770 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n goto 730
\n780 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n do 790 i = 1,n
\n790 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n p = fa
\n goto 740
\n end
\n subroutine cgn(x,e,it,st,io,tl,lm1,lm2,n,is,vlf,vlg,grf,grg,nl,pn,
\n 1 aa,la,m,bl,bu,ct,u,ia,k,le,vf,vg,df,dg,lg,gg,bb,cc,dd,h)
\n integer ia(1),i,ib,io,is,it,j,k,la,lg,lm1,lm2,m,n,nl,na,nb,nc,nd
\n real*8 df(1),dg(lg,1),gg(1),h(n,1),le(1),vg(1),x(1),y(50),z(50)
\n real*8 aa(1),bb(1),bl(1),bu(1),cc(1),dd(1),u(1)
\n real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0,da,db
\n real*8 e,f,f0,f1,fa,fb,fc,g,l3,p,p2,pn,q,r,s,st,tl,v,vf,vv,w
\n real*8 pv,pd,pd2
\n external grf,grg,vlf,vlg
\n data a1\/.1d0\/,a2\/.9d0\/,a3\/5.d0\/,a4\/.2d0\/,a5\/10.d0\/,a6\/.9d0\/
\n data a7\/.3d0\/
\nc
\nc set big = largest floating point number
\nc
\n big = 1.d50
\n a8 = a3 + .01d0
\n p2 = .5*pn
\n it = 0
\n io = 0
\n f = vf
\n if ( is .eq. 0 ) goto 40
\n do 10 i = 1,nl
\n dd(i) = vg(i)
\n10 f = f + vg(i)*le(i)
\n do 30 j = 1,n
\n s = df(j)
\n do 20 i = 1,nl
\n20 s = s + le(i)*dg(i,j)
\n h(j,3) = s
\n30 continue
\n goto 80
\n40 do 50 i = 1,nl
\n dd(i) = 0.
\n bb(i) = le(i) + pn*vg(i)
\n50 f = f + vg(i)*le(i) + p2*vg(i)**2
\n do 70 j = 1,n
\n s = df(j)
\n do 60 i = 1,nl
\n60 s = s + bb(i)*dg(i,j)
\n h(j,3) = s
\n70 continue
\n80 l3 = 1.\/dlog(a3)
\n call pre(h(1,2),e,u,ia,k,i,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n if ( e .le. tl ) goto 770
\n if ( io .gt. 0 ) goto 770
\n if ( lm1 .eq. 0 ) return
\n a = st
\n if ( a .gt. 0. ) goto 100
\n do 90 i = 1,n
\n90 if ( dabs(x(i)) .gt. a ) a = dabs(x(i))
\n a = .01*a\/e
\n if ( a .eq. 0. ) a = 1.
\n100 g = 0.
\n do 110 i = 1,n
\n h(i,1) = -h(i,2)
\n110 g = g + h(i,2)*h(i,3)
\n if ( g .lt. 0. ) goto 730
\n d = -g
\n120 call cut(s,ib,x,h,ia,bl,bu,ct,n,big)
\n na = 0
\n nb = 0
\n nc = 0
\n nd = 0
\n if ( s .le. 0. ) goto 860
\n if ( a .gt. s ) a = s
\n p = pv(a,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n fa = p
\n c0 = a
\n f0 = fa
\n ny = 2
\n y(1) = 0.
\n z(1) = f
\n y(2) = a
\n z(2) = fa
\n v = a1*d
\n w = a2*d
\n iq = 0
\n if ( fa .le. f ) goto 130
\n c = a
\n b = 0.
\n a = 0.
\n fc = fa
\n fb = f
\n fa = f
\n goto 150
\n130 c = 0.
\n b = 0.
\n fc = f
\n fb = f
\n iq = 1
\n vf = vv
\n do 140 i = 1,nl
\n140 vg(i) = gg(i)
\n150 q = (d+(f-f0)\/c0)\/c0
\n if ( q .lt. 0. ) goto 170
\n q = a
\n160 nd = nd + 1
\n if ( nd .gt. 25 ) goto 720
\n if ( s .eq. q ) goto 780
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 160
\n goto 350
\n170 q = .5*d\/q
\n if ( q .lt. s ) goto 190
\n if ( c0 .lt. s ) goto 180
\n f1 = f0
\n c1 = c0
\n q = s
\n goto 220
\n180 q = s
\n190 if ( q .lt. .01*c0 ) q = .01*c0
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n if ( p .le. f0 ) goto 200
\n f1 = f0
\n c1 = c0
\n f0 = p
\n c0 = q
\n goto 210
\n200 f1 = p
\n c1 = q
\n210 call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n220 if ( a .eq. 0. ) goto 230
\n if ( fa-f .ge. v*a ) goto 250
\n if ( fa-f .lt. w*a ) goto 300
\n goto 370
\n230 q = c0
\n if ( c1 .lt. q ) q = c1
\n240 na = na + 1
\n if ( na .gt. 25 ) goto 740
\n q = a4*q
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .ge. v*q ) goto 240
\n goto 340
\n250 if ( c0 .gt. c1 ) goto 290
\n if ( f0-f .gt. v*c0 ) goto 270
\n if ( f0-f .ge. w*c0 ) goto 430
\n if ( c1 .le. a5*c0 ) goto 430
\n r = dlog(c1\/c0)
\n r = .999*dexp(-r\/idint(r*l3+.999))
\n q = c1
\n260 q = q*r
\n if ( q .lt. c0 ) goto 430
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n na = na + 1
\n if ( p-f .gt. v*q ) goto 260
\n goto 430
\n270 q = c0
\n280 na = na + 1
\n if ( na .gt. 25 ) goto 740
\n q = a4*q
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .ge. v*q ) goto 280
\n goto 340
\n290 q = a
\n goto 280
\n300 if ( c0 .lt. c1 ) goto 400
\n if ( f0-f .ge. v*c0 ) goto 320
\n if ( f0-f .ge. w*c0 ) goto 340
\n q = c0
\n p = f0
\n310 nd = nd + 1
\n if ( nd .gt. 25 ) goto 720
\n if ( s .eq. q ) goto 780
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 310
\n goto 340
\n320 if ( c0 .le. a5*c1 ) goto 340
\n r = dlog(c0\/c1)
\n r = 1.001*dexp(r\/idint(r*l3+.999))
\n q = a
\n330 q = q*r
\n if ( q .gt. c0 ) goto 340
\n nd = nd + 1
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 330
\n340 if ( iq .eq. 1 ) goto 430
\n350 if ( b .eq. 0. ) goto 370
\n if ( c .eq. 0. ) goto 360
\n v = c – a
\n w = a – b
\n r = 1.\/v
\n f0 = 1.\/w
\n p = fc – fa
\n q = fb – fa
\n e = p*r + q*f0
\n if ( dsign(e,c-b) .ne. e ) goto 430
\n if ( e .eq. 0. ) goto 430
\n q = (p*r)*w – (q*f0)*v
\n q = a – .5*q\/e
\n goto 380
\n360 r = 1.\/a
\n f0 = 1.\/b
\n p = r*(fa-f) – d
\n q = f0*(fb-f) – d
\n e = a – b
\n v = (r*p-f0*q)\/e
\n w = (a*q*f0-b*p*r)\/e
\n v = w*w-3.*v*d
\n if ( v .lt. 0. ) v = 0.
\n v = dsqrt(v)
\n if ( w+v .eq. 0. ) goto 430
\n q = -d\/(w+v)
\n if ( q .le. 0. ) goto 430
\n goto 380
\n370 if ( iq .eq. 1 ) goto 430
\n q = (d+(f-fa)\/a)\/a
\n if ( q .ge. 0. ) goto 430
\n q = .5*d\/q
\n380 if ( q .gt. s ) q = s
\n do 390 i = 1,ny
\n390 if ( y(i) .eq. q ) goto 430
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n goto 430
\n400 continue
\n if ( f0-f .gt. v*c0 ) goto 410
\n if ( f0-f .gt. w*c0 ) goto 430
\n410 q = a
\n p = fa
\n420 nd = nd + 1
\n if ( nd .gt. 25 ) goto 720
\n if ( s .eq. q ) goto 780
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 420
\n goto 340
\n430 da = pd(a,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb)
\n if ( da .gt. a6*g ) goto 520
\n if ( s .eq. a ) goto 880
\n if ( da .ge. 0. ) goto 670
\n r = a
\n q = 0.
\n do 440 i = 1,ny
\n if ( y(i) .gt. a ) goto 480
\n if ( y(i) .le. q ) goto 440
\n if ( y(i) .eq. a ) goto 440
\n q = y(i)
\n440 continue
\n if ( a .le. a8*q ) goto 670
\n q = a
\n p = fa
\n450 nd = nd + 1
\n if ( nd .gt. 25 ) goto 720
\n if ( s .eq. q ) goto 780
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n f1 = fa
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p .lt. f1 ) goto 450
\n if ( a .gt. r ) goto 470
\n do 460 i = 1,n
\n460 h(i,2) = x(i) + a*h(i,1)
\n goto 670
\n470 da = pd(a,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb)
\n if ( da .gt. a6*g ) goto 520
\n goto 670
\n480 q = y(i)
\n do 490 j = i,ny
\n if ( y(j) .le. a ) goto 490
\n if ( y(j) .lt. q ) q = y(j)
\n490 continue
\n if ( q .le. a5*a ) goto 670
\n f0 = dlog(q\/a)
\n f0 = 1.001*dexp(f0\/idint(f0*l3+.999))
\n v = a
\n500 v = v*f0
\n if ( v .ge. q ) goto 430
\n p = pv(v,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n f1 = fa
\n call igs(v,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p .lt. f1 ) goto 500
\n if ( a .gt. r ) goto 430
\n do 510 i = 1,n
\n510 h(i,2) = x(i) + a*h(i,1)
\n goto 670
\n520 b = 0.
\n j = 1
\n i = j
\n530 i = i + 1
\n if ( i .gt. ny ) goto 540
\n if ( y(i) .ge. a ) goto 530
\n if ( y(i) .lt. b ) goto 530
\n b = y(i)
\n j = i
\n goto 530
\n540 fb = z(j)
\n db = d
\n if ( b .ne. 0. )
\n 1 db = pd2(b,x,h,df,dg,gg,dd,nl,le,vlg,grf,grg,pn,n,lg,bb)
\n550 w = 2.*dabs(b-a)
\n call cub(c,a,b,fa,fb,da,db)
\n nc = 1
\n goto 590
\n560 w = .5*w
\n if ( w .lt. dabs(c0-c) ) goto 660
\n if ( c0 .lt. c ) goto 570
\n if ( d0 .ge. d ) goto 580
\n goto 660
\n570 if ( d0 .gt. d ) goto 660
\n580 call cub(c,c,c0,f,f0,d,d0)
\n nc = nc + 1
\n if ( nc .gt. 30 ) goto 710
\n590 r = dmax1(a,b)
\n s = dmin1(a,b)
\n if ( c .gt. r ) goto 600
\n if ( c .gt. s ) goto 610
\n c = s + (s-c)
\n s = .5*(a+b)
\n if ( c .gt. s ) c = s
\n goto 610
\n600 c = r – (c-r)
\n s = .5*(a+b)
\n if ( c .lt. s ) c = s
\n610 c0 = a
\n f0 = fa
\n d0 = da
\n f = pv(c,x,h,vf,vg,dd,nl,le,vlf,vlg,p2,n)
\n d = pd(c,x,h,df,dg,vg,dd,nl,le,grf,grg,pn,n,lg,bb)
\n if ( f .lt. fa ) goto 620
\n b = c
\n fb = f
\n db = d
\n goto 560
\n620 if ( c .lt. a ) goto 650
\n if ( d .lt. 0. ) goto 640
\n630 b = a
\n fb = fa
\n db = da
\n640 a = c
\n fa = f
\n da = d
\n if ( d .gt. a6*g ) goto 560
\n goto 670
\n650 if ( d .lt. 0. ) goto 630
\n goto 640
\n660 c = .5*(a+b)
\n nb = nb + 1
\n w = dabs(b-a)
\n goto 610
\n670 do 680 i = 1,n
\n680 x(i) = h(i,2)
\n i = lm1
\n call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n it = it + 1
\n f = fa
\n d = da
\n a = a7*a
\n if ( io .gt. 0 ) goto 770
\n if ( e .le. tl ) goto 770
\n if ( it .ge. lm2 ) goto 770
\n if ( it .ge. lm1 ) goto 770
\n if ( i .lt. lm1 ) goto 100
\n r = 0.
\n do 690 i = 1,n
\n690 r = r + h(i,2)*h(i,3)
\n if ( r .lt. 0. ) goto 730
\n s = r\/g
\n g = r
\n d = 0.
\n do 700 i = 1,n
\n h(i,1) = -h(i,2) + s*h(i,1)
\n700 d = d + h(i,1)*h(i,3)
\n goto 120
\n710 if ( d .lt. g ) goto 670
\n io = 4
\n return
\n720 io = 3
\n return
\n730 io = 4
\n return
\n740 q = q*a3**25
\n nd = 0
\n750 nd = nd + 1
\n if ( nd .gt. 25 ) goto 760
\n q = a3*q
\n if ( q .ge. s ) goto 760
\n p = pv(q,x,h,vv,gg,dd,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .gt. v*q ) goto 750
\n goto 220
\n760 io = 5
\n return
\n770 st = a
\n return
\n780 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n do 790 i = 1,n
\n790 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n a = q*a7
\n800 call grf(df,x)
\n call grg(dg,x)
\n810 do 820 i = 1,nl
\n820 bb(i) = le(i) + pn*(vg(i)-dd(i))
\n do 840 j = 1,n
\n s = df(j)
\n do 830 i = 1,nl
\n830 s = s + bb(i)*dg(i,j)
\n h(j,3) = s
\n840 continue
\n f = p
\n850 call pre(h(1,2),e,u,ia,k,j,lm1,io,h(1,3),aa,la,m,n,bb,cc)
\n it = it + 1
\n if ( it .ge. lm2 ) goto 770
\n if ( it .ge. lm1 ) goto 770
\n if ( io .gt. 0 ) goto 770
\n if ( e .le. tl ) goto 770
\n goto 100
\n860 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n q = -s
\n p = pv(q,x,h,vf,vg,dd,nl,le,vlf,vlg,p2,n)
\n do 870 i = 1,n
\n870 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n goto 800
\n880 call mod(u,ia,k,io,aa,la,m,ib,1,bb)
\n do 890 i = 1,n
\n890 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n p = fa
\n goto 810
\n end
\n subroutine cgp(x,e,it,st,io,tl,lm1,lm2,n,vlf,vlg,grf,grg,l,nl,pn,
\n 1 aa,la,bl,bu,ct,u,ia,k,le,vf,vg,df,dg,lg,uu,gg,bb,cc,h)
\n integer ia(1),i,ib,io,it,j,k,l,la,lg,lm1,lm2,n,nl,na,nb,nc,nd
\n real*8 df(1),dg(lg,1),gg(1),h(n,1),le(1),vg(1),x(1),y(50),z(50)
\n real*8 aa(1),bb(1),bl(1),bu(1),cc(1),u(1),uu(1)
\n real*8 a1,a2,a3,a4,a5,a6,a7,a8,a,b,big,c,c0,c1,ct,d,d0,da,db
\n real*8 e,f,f0,f1,fa,fb,fc,g,l3,p,p2,pn,q,r,s,st,tl,v,vf,vv,w
\n real*8 qv,qd,qd2
\n external grf,grg,vlf,vlg
\n data a1\/.1d0\/,a2\/.9d0\/,a3\/5.d0\/,a4\/.2d0\/,a5\/10.d0\/,a6\/.9d0\/
\n data a7\/.3d0\/
\nc
\nc set big = largest floating point number
\nc
\n i = l
\n if ( l .eq. 0 ) i = 1
\n call fab(uu,u,ia,io,aa,la,dg,lg,pn,n,i,l,nl,bb,cc)
\n if ( io .gt. 0 ) return
\n big = 1.d50
\n a8 = a3 + .01d0
\n p2 = .5*pn
\n it = 0
\n io = 0
\n f = vf
\n do 10 i = 1,nl
\n cc(i) = le(i) + pn*vg(i)
\n10 f = f + (le(i)+p2*vg(i))*vg(i)
\n do 30 j = 1,n
\n s = df(j)
\n do 20 i = 1,nl
\n20 s = s + cc(i)*dg(i,j)
\n h(j,3) = s
\n30 continue
\n l3 = 1.\/dlog(a3)
\n call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc)
\n if ( e .le. tl ) goto 720
\n if ( io .gt. 0 ) goto 720
\n if ( lm1 .eq. 0 ) return
\n a = st
\n if ( a .gt. 0. ) goto 50
\n do 40 i = 1,n
\n40 if ( dabs(x(i)) .gt. a ) a = dabs(x(i))
\n a = .01*a\/e
\n if ( a .eq. 0. ) a = 1.
\n50 g = 0.
\n do 60 i = 1,n
\n h(i,1) = -h(i,2)
\n60 g = g + h(i,2)*h(i,3)
\n if ( g .lt. 0. ) goto 680
\n d = -g
\n70 call cut(s,ib,x,h,ia,bl,bu,ct,n,big)
\n na = 0
\n nb = 0
\n nc = 0
\n nd = 0
\n if ( s .le. 0. ) goto 790
\n if ( a .gt. s ) a = s
\n p = qv(a,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n fa = p
\n c0 = a
\n f0 = fa
\n ny = 2
\n y(1) = 0.
\n z(1) = f
\n y(2) = a
\n z(2) = fa
\n v = a1*d
\n w = a2*d
\n iq = 0
\n if ( fa .le. f ) goto 80
\n c = a
\n b = 0.
\n a = 0.
\n fc = fa
\n fb = f
\n fa = f
\n goto 100
\n80 c = 0.
\n b = 0.
\n fc = f
\n fb = f
\n iq = 1
\n vf = vv
\n do 90 i = 1,nl
\n90 vg(i) = gg(i)
\n100 q = (d+(f-f0)\/c0)\/c0
\n if ( q .lt. 0. ) goto 120
\n q = a
\n110 nd = nd + 1
\n if ( nd .gt. 25 ) goto 670
\n if ( s .eq. q ) goto 730
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 110
\n goto 300
\n120 q = .5*d\/q
\n if ( q .lt. s ) goto 140
\n if ( c0 .lt. s ) goto 130
\n f1 = f0
\n c1 = c0
\n q = s
\n goto 170
\n130 q = s
\n140 if ( q .lt. .01*c0 ) q = .01*c0
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n if ( p .le. f0 ) goto 150
\n f1 = f0
\n c1 = c0
\n f0 = p
\n c0 = q
\n goto 160
\n150 f1 = p
\n c1 = q
\n160 call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n170 if ( a .eq. 0. ) goto 180
\n if ( fa-f .ge. v*a ) goto 200
\n if ( fa-f .lt. w*a ) goto 250
\n goto 320
\n180 q = c0
\n if ( c1 .lt. q ) q = c1
\n190 na = na + 1
\n if ( na .gt. 25 ) goto 690
\n q = a4*q
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .ge. v*q ) goto 190
\n goto 290
\n200 if ( c0 .gt. c1 ) goto 240
\n if ( f0-f .gt. v*c0 ) goto 220
\n if ( f0-f .ge. w*c0 ) goto 380
\n if ( c1 .le. a5*c0 ) goto 380
\n r = dlog(c1\/c0)
\n r = .999*dexp(-r\/idint(r*l3+.999))
\n q = c1
\n210 q = q*r
\n if ( q .lt. c0 ) goto 380
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n na = na + 1
\n if ( p-f .gt. v*q ) goto 210
\n goto 380
\n220 q = c0
\n230 na = na + 1
\n if ( na .gt. 25 ) goto 690
\n q = a4*q
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .ge. v*q ) goto 230
\n goto 290
\n240 q = a
\n goto 230
\n250 if ( c0 .lt. c1 ) goto 350
\n if ( f0-f .ge. v*c0 ) goto 270
\n if ( f0-f .ge. w*c0 ) goto 290
\n q = c0
\n p = f0
\n260 nd = nd + 1
\n if ( nd .gt. 25 ) goto 670
\n if ( s .eq. q ) goto 730
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 260
\n goto 290
\n270 if ( c0 .le. a5*c1 ) goto 290
\n r = dlog(c0\/c1)
\n r = 1.001*dexp(r\/idint(r*l3+.999))
\n q = a
\n280 q = q*r
\n if ( q .gt. c0 ) goto 290
\n nd = nd + 1
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 280
\n290 if ( iq .eq. 1 ) goto 380
\n300 if ( b .eq. 0. ) goto 320
\n if ( c .eq. 0. ) goto 310
\n v = c – a
\n w = a – b
\n r = 1.\/v
\n f0 = 1.\/w
\n p = fc – fa
\n q = fb – fa
\n e = p*r + q*f0
\n if ( dsign(e,c-b) .ne. e ) goto 380
\n if ( e .eq. 0. ) goto 380
\n q = (p*r)*w – (q*f0)*v
\n q = a – .5*q\/e
\n goto 330
\n310 r = 1.\/a
\n f0 = 1.\/b
\n p = r*(fa-f) – d
\n q = f0*(fb-f) – d
\n e = a – b
\n v = (r*p-f0*q)\/e
\n w = (a*q*f0-b*p*r)\/e
\n v = w*w-3.*v*d
\n if ( v .lt. 0. ) v = 0.
\n v = dsqrt(v)
\n if ( w+v .eq. 0. ) goto 380
\n q = -d\/(w+v)
\n if ( q .le. 0. ) goto 380
\n goto 330
\n320 if ( iq .eq. 1 ) goto 380
\n q = (d+(f-fa)\/a)\/a
\n if ( q .ge. 0. ) goto 380
\n q = .5*d\/q
\n330 if ( q .gt. s ) q = s
\n do 340 i = 1,ny
\n340 if ( y(i) .eq. q ) goto 380
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n goto 380
\n350 continue
\n if ( f0-f .gt. v*c0 ) goto 360
\n if ( f0-f .gt. w*c0 ) goto 380
\n360 q = a
\n p = fa
\n370 nd = nd + 1
\n if ( nd .gt. 25 ) goto 670
\n if ( s .eq. q ) goto 730
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .lt. w*q ) goto 370
\n goto 290
\n380 da = qd(a,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc)
\n if ( da .gt. a6*g ) goto 470
\n if ( s .eq. a ) goto 800
\n if ( da .ge. 0. ) goto 620
\n r = a
\n q = 0.
\n do 390 i = 1,ny
\n if ( y(i) .gt. a ) goto 430
\n if ( y(i) .le. q ) goto 390
\n if ( y(i) .eq. a ) goto 390
\n q = y(i)
\n390 continue
\n if ( a .le. a8*q ) goto 620
\n q = a
\n p = fa
\n400 nd = nd + 1
\n if ( nd .gt. 25 ) goto 670
\n if ( s .eq. q ) goto 730
\n q = a3*q
\n if ( q .gt. s ) q = s
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n f1 = fa
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p .lt. f1 ) goto 400
\n if ( a .gt. r ) goto 420
\n do 410 i = 1,n
\n410 h(i,2) = x(i) + a*h(i,1)
\n goto 620
\n420 da = qd(a,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc)
\n if ( da .gt. a6*g ) goto 470
\n goto 620
\n430 q = y(i)
\n do 440 j = i,ny
\n if ( y(j) .le. a ) goto 440
\n if ( y(j) .lt. q ) q = y(j)
\n440 continue
\n if ( q .le. a5*a ) goto 620
\n f0 = dlog(q\/a)
\n f0 = 1.001*dexp(f0\/idint(f0*l3+.999))
\n v = a
\n450 v = v*f0
\n if ( v .ge. q ) goto 380
\n p = qv(v,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n f1 = fa
\n call igs(v,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p .lt. f1 ) goto 450
\n if ( a .gt. r ) goto 380
\n do 460 i = 1,n
\n460 h(i,2) = x(i) + a*h(i,1)
\n goto 620
\n470 b = 0.
\n j = 1
\n i = j
\n480 i = i + 1
\n if ( i .gt. ny ) goto 490
\n if ( y(i) .ge. a ) goto 480
\n if ( y(i) .lt. b ) goto 480
\n b = y(i)
\n j = i
\n goto 480
\n490 fb = z(j)
\n db = d
\n if ( b .ne. 0. )
\n 1 db = qd2(b,x,h,df,dg,gg,nl,le,vlg,grf,grg,pn,n,lg,cc)
\n500 w = 2.*dabs(b-a)
\n call cub(c,a,b,fa,fb,da,db)
\n nc = 1
\n goto 540
\n510 w = .5*w
\n if ( w .lt. dabs(c0-c) ) goto 610
\n if ( c0 .lt. c ) goto 520
\n if ( d0 .ge. d ) goto 530
\n goto 610
\n520 if ( d0 .gt. d ) goto 610
\n530 call cub(c,c,c0,f,f0,d,d0)
\n nc = nc + 1
\n if ( nc .gt. 30 ) goto 660
\n540 r = dmax1(a,b)
\n s = dmin1(a,b)
\n if ( c .gt. r ) goto 550
\n if ( c .gt. s ) goto 560
\n c = s + (s-c)
\n s = .5*(a+b)
\n if ( c .gt. s ) c = s
\n goto 560
\n550 c = r – (c-r)
\n s = .5*(a+b)
\n if ( c .lt. s ) c = s
\n560 c0 = a
\n f0 = fa
\n d0 = da
\n f = qv(c,x,h,vf,vg,nl,le,vlf,vlg,p2,n)
\n d = qd(c,x,h,df,dg,vg,nl,le,grf,grg,pn,n,lg,cc)
\n if ( f .lt. fa ) goto 570
\n b = c
\n fb = f
\n db = d
\n goto 510
\n570 if ( c .lt. a ) goto 600
\n if ( d .lt. 0. ) goto 590
\n580 b = a
\n fb = fa
\n db = da
\n590 a = c
\n fa = f
\n da = d
\n if ( d .gt. a6*g ) goto 510
\n goto 620
\n600 if ( d .lt. 0. ) goto 580
\n goto 590
\n610 c = .5*(a+b)
\n nb = nb + 1
\n w = dabs(b-a)
\n goto 560
\n620 do 630 i = 1,n
\n630 x(i) = h(i,2)
\n i = lm1
\n call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc)
\n it = it + 1
\n f = fa
\n d = da
\n a = a7*a
\n if ( io .gt. 0 ) goto 720
\n if ( e .le. tl ) goto 720
\n if ( it .ge. lm2 ) goto 720
\n if ( it .ge. lm1 ) goto 720
\n if ( i .lt. lm1 ) goto 50
\n r = 0.
\n do 640 i = 1,n
\n640 r = r + h(i,2)*h(i,3)
\n if ( r .lt. 0. ) goto 680
\n s = r\/g
\n g = r
\n d = 0.
\n do 650 i = 1,n
\n h(i,1) = -h(i,2) + s*h(i,1)
\n650 d = d + h(i,1)*h(i,3)
\n goto 70
\n660 if ( d .lt. g ) goto 620
\n io = 4
\n return
\n670 io = 3
\n return
\n680 io = 4
\n return
\n690 q = q*a3**25
\n nd = 0
\n700 nd = nd + 1
\n if ( nd .gt. 25 ) goto 710
\n q = a3*q
\n if ( q .ge. s ) goto 710
\n p = qv(q,x,h,vv,gg,nl,le,vlf,vlg,p2,n)
\n call igs(q,p,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,ny,y,z)
\n if ( p-f .gt. v*q ) goto 700
\n goto 170
\n710 io = 5
\n return
\n720 st = a
\n return
\n730 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc)
\n if ( io .gt. 0 ) goto 720
\n call mod(u,ia,k,io,aa,la,l,ib,1,bb)
\n do 740 i = 1,n
\n740 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n a = q*a7
\n745 call grf(df,x)
\n call grg(dg,x)
\n746 do 750 i = 1,nl
\n750 cc(i) = le(i) + pn*vg(i)
\n do 770 j = 1,n
\n s = df(j)
\n do 760 i = 1,nl
\n760 s = s + cc(i)*dg(i,j)
\n h(j,3) = s
\n770 continue
\n f = p
\n780 call prp(h(1,2),e,u,uu,ia,k,lm1,io,h(1,3),aa,la,l,nl,n,bb,cc)
\n it = it + 1
\n if ( it .ge. lm2 ) goto 720
\n if ( it .ge. lm1 ) goto 720
\n if ( io .gt. 0 ) goto 720
\n if ( e .le. tl ) goto 720
\n goto 50
\n790 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc)
\n if ( io .gt. 0 ) goto 720
\n call mod(u,ia,k,io,aa,la,l,ib,1,bb)
\n q = -s
\n p = qv(q,x,h,vf,vg,nl,le,vlf,vlg,p2,n)
\n do 795 i = 1,n
\n795 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n goto 745
\n800 call mob(uu,io,u,ia,k,aa,la,n,l,nl,ib,4,bb,cc)
\n if ( io .gt. 0 ) goto 720
\n call mod(u,ia,k,io,aa,la,l,ib,1,bb)
\n do 810 i = 1,n
\n810 x(i) = h(i,2)
\n x(ib) = bl(ib)
\n if ( ia(ib) .gt. 0 ) x(ib) = bu(ib)
\n goto 746
\n end
\n subroutine cub(x,a,b,c,d,e,f)
\n real*8 a,b,c,d,e,f,g,v,w,x,y,z
\n g = b – a
\n v = e + f – 3*(d-c)\/g
\n w = v*v-e*f
\n if ( w .lt. 0. ) w = 0.
\n w = dsign(dsqrt(w),g)
\n y = e + v
\n z = f + v
\n if ( dsign(y,g) .ne. y ) goto 30
\n if ( dsign(z,g) .ne. z ) goto 20
\n if ( z .eq. 0. ) goto 20
\n10 x = b – g*f\/(z+w)
\n return
\n20 if ( c .lt. d ) x = a
\n if ( c .ge. d ) x = b
\n return
\n30 if ( dsign(z,g) .ne. z ) goto 40
\n if ( dabs(e) .gt. dabs(f) ) goto 10
\n40 x = a + g*e\/(y-w)
\n return
\n end
\n subroutine cut(s,ib,x,d,ia,bl,bu,ct,n,big)
\n real*8 bl(1),bu(1),d(1),x(1),big,ct,s,t
\n integer ia(1),i,ib,k
\n s = big
\n k = 0
\n do 20 i = 1,n
\n if ( ia(i) .ne. 0 ) goto 20
\n if ( d(i) .eq. 0. ) goto 20
\n if ( d(i) .gt. 0. ) goto 10
\n t = x(i) – bl(i)
\n if ( t .le. ct ) k = 1
\n t = -t\/d(i)
\n if ( t .gt. s ) goto 20
\n s = t
\n ib = -i
\n goto 20
\n10 t = bu(i) – x(i)
\n if ( t .le. ct ) k = 1
\n t = t\/d(i)
\n if ( t .gt. s ) goto 20
\n s = t
\n ib = i
\n20 continue
\n if ( s .lt. 0. ) return
\n if ( k .eq. 1 ) goto 30
\n s = dabs(s)
\n return
\n30 s = -s
\n return
\n end
\n subroutine erg(t,g,nl)
\n real*8 g(1),t
\n integer i,nl
\n t = 0.
\n do 10 i = 1,nl
\n10 t = dmax1(dabs(g(i)),t)
\n return
\n end
\n subroutine err(io)
\n integer io
\n goto (10,20,30,40,50,60,70,80,90), io
\n write(6,*) ‘no detected errors’
\n return
\n10 write(6,*) ‘dependent constraint gradients’
\n return
\n20 write(6,*) ‘problem is infeasible’
\n return
\n30 write(6,*) ‘function decreases with no minimum’
\n return
\n40 write(6,*) ‘inconsistent gradient values’
\n return
\n50 write(6,*) ‘unable to satisfy Armijo condition’
\n return
\n60 write(6,*) ‘cannot reduce constraint error below gradient error’
\n return
\n70 write(6,*) ‘input penalty is zero but a positive value is needed’
\n write(6,*) ‘for convergence in this problem’
\n return
\n80 write(6,*) ‘more constraints imposed than there are unknowns’
\n return
\n90 write(6,*) ‘iterations at limit’
\n return
\n end
\n subroutine fab(v,u,ia,io,a,la,dg,lg,pn,n,ll,l,nl,b,w)
\n integer ia(1),i,io,j,k,l,ll,la,lg,nl,n,o
\n real*8 a(la,1),dg(lg,1),b(1),u(1),v(1),w(ll,1),pn,t
\n do 10 j = 1,n
\n do 10 i = 1,nl
\n10 a(l+i,j) = dg(i,j)
\n j = (nl*nl+nl)\/2
\n do 20 i = 1,j
\n20 v(i) = 0.
\n o = 1
\n do 30 i = 1,nl
\n v(o) = 1.\/pn
\n o = o + nl – i + 1
\n30 continue
\n do 60 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 60
\n o = 0
\n do 50 i = 1,nl
\n t = dg(i,j)
\n do 40 k = i,nl
\n o = o + 1
\n v(o) = v(o) + t*dg(k,j)
\n40 continue
\n50 continue
\n60 continue
\n if ( l .eq. 0 ) goto 140
\n do 70 j = 1,nl
\n do 70 i = 1,l
\n70 w(i,j) = 0.
\n do 100 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 100
\n do 90 i = 1,nl
\n t = dg(i,j)
\n do 80 k = 1,l
\n w(k,i) = w(k,i) + t*a(k,j)
\n80 continue
\n90 continue
\n100 continue
\n o = 0
\n do 130 j = 1,nl
\n do 110 i = 1,l
\n110 b(i) = w(i,j)
\n call sol(b,u,l,b)
\n do 130 i = j,nl
\n t = 0.
\n do 120 k = 1,l
\n120 t = t + b(k)*w(k,i)
\n o = o + 1
\n v(o) = v(o) – t
\n130 continue
\n140 call fac(v,0,nl,io)
\n return
\n end
\n subroutine fac(a,p,n,io)
\n real*8 a(1),s,t
\n integer e,f,g,h,i,io,j,k,l,m,n,o,p,q
\n io = 0
\n if ( n .eq. 0 ) return
\n q = ((n+n+1-p)*p)\/2
\n l = p
\n if ( l .lt. 1 ) l = 1
\n h = n
\n k = 1
\n10 s = a(k)
\n if ( s .le. 0. ) goto 50
\n if ( k .gt. q ) s = dsqrt(s)
\n a(k) = s
\n if ( h .eq. 1 ) return
\nc ————————–
\nc |*** save pivot entry ***|
\nc ————————–
\n e = k + l
\n f = k + h – 1
\n do 20 i = e,f
\n20 a(i) = a(i)\/s
\n o = l
\n if ( l .gt. 1 ) l = l – 1
\n k = k + h
\n g = k
\n h = h – 1
\n m = h
\n j = 0
\n30 if ( o .gt. 0 ) o = o – 1
\n e = g + o
\n j = j – m
\n m = m – 1
\n f = g + m
\n t = a(g+j)
\nc —————————
\nc |*** eliminate by rows ***|
\nc —————————
\n do 40 i = e,f
\n40 a(i) = a(i) – t*a(i+j)
\n g = f + 1
\n if ( m .gt. 0 ) goto 30
\n goto 10
\n50 io = 1
\n return
\n end
\n subroutine fat(u,io,in,ia,a,la,l,n,g,h,z,dg,lg,nl,m)
\n integer ia(1),i,in,io,j,k,l,lg,m,n,nl,o,p
\n real*8 a(la,1),dg(lg,1),g(1),h(1),u(1),z(1),t
\nc
\nc normalize rows of g
\nc
\n do 10 i = 1,nl
\n10 h(i) = 0.
\n do 20 j = 1,n
\n do 20 i = 1,nl
\n20 if ( dabs(dg(i,j)) .gt. h(i) ) h(i) = dabs(dg(i,j))
\n do 30 i = 1,nl
\n if ( h(i) .eq. 0. ) goto 120
\n h(i) = 1.\/h(i)
\n30 z(i) = h(i)
\n do 40 j = 1,n
\n do 40 i = 1,nl
\n40 a(i+l,j) = h(i)*dg(i,j)
\n do 50 i = 1,nl
\n50 h(i) = -h(i)*g(i)
\n o = l
\n if ( l .eq. 0 ) goto 110
\n p = o
\n do 70 k = 1,l
\n do 60 i = 1,nl
\n o = o + 1
\n60 u(o) = 0.
\n p = p – 1
\n o = o + p
\n70 continue
\n do 100 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 100
\n o = l
\n p = o
\n do 90 k = 1,l
\n t = a(k,j)
\n do 80 i = 1,nl
\n o = o + 1
\n80 u(o) = u(o) + t*a(i+l,j)
\n p = p – 1
\n o = o + p
\n90 continue
\n100 continue
\n110 call mat(u(o+1),ia,a(l+1,1),la,nl,n,in)
\n call fac(u,l,m,io)
\n return
\n120 io = 1
\n return
\n end
\n subroutine fea(x,u,ia,k,io,y,a,la,m,n,b,bl,bu,c,d,g)
\n integer ia(1),h,i,io,it,j,k,l,la,m,n
\n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),u(1)
\n double precision x(1),y(1),b1,b2,big,e,f,q,r,s,t,v,w
\nc
\nc set big = largest floating point number
\nc
\n big = 1.d50
\n io = 0
\n e = 1.d-3
\n f = .999d0
\n if ( m .gt. 0 ) goto 30
\n do 20 i = 1,n
\n ia(i) = 0
\n x(i) = y(i)
\n if ( x(i) .lt. bu(i) ) goto 10
\n x(i) = bu(i)
\n ia(i) = 1
\n goto 20
\n10 if ( x(i) .gt. bl(i) ) goto 20
\n x(i) = bl(i)
\n ia(i) = -1
\n20 continue
\n return
\n30 call sol(b,u,m,b)
\n do 50 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 50
\n t = y(j)
\n do 40 i = 1,m
\n40 t = t + b(i)*a(i,j)
\n x(j) = t
\n50 continue
\n l = 0
\n it = -1
\n do 70 i = 1,n
\n d(i) = 0.
\n if ( iabs(ia(i)) .eq. 1 ) goto 70
\n if ( x(i) .lt. bl(i) ) goto 60
\n if ( x(i) .le. bu(i) ) goto 70
\n ia(i) = 2
\n l = l + 1
\n goto 70
\n60 ia(i) = -2
\n l = l + 1
\n70 continue
\n b2 = 0.
\n if ( l .gt. 0 ) goto 190
\n do 80 i = 1,n
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n80 continue
\n return
\n90 r = 0.
\n s = 0.
\n t = big
\n h = 0
\n if ( it .eq. 0 ) q = 0.
\n if ( it .gt. 0 ) q = b2\/b1
\n do 120 i = 1,n
\n if ( iabs(ia(i)) .eq. 1 ) goto 120
\n d(i) = q*d(i) – g(i)
\n s = s + d(i)*g(i)
\n w = bu(i)
\n if ( d(i) .gt. 0. ) goto 100
\n if ( d(i) .eq. 0. ) goto 120
\n if ( ia(i) .lt. 0 ) goto 120
\n if ( ia(i) .eq. 0 ) w = bl(i)
\n if ( ia(i) .gt. 0 ) r = r + d(i)**2
\n goto 110
\n100 if ( ia(i) .gt. 0 ) goto 120
\n if ( ia(i) .eq. 0 ) goto 110
\n r = r + d(i)**2
\n w = bl(i)
\n110 v = (w-x(i))\/d(i)
\n if ( v .ge. t ) goto 120
\n t = v
\n h = i
\n120 continue
\n s = -s\/r
\n if ( t .le. s ) goto 130
\n if ( h .eq. 0 ) goto 170
\n if ( ia(h) .eq. 0 ) goto 170
\n130 s = t
\n if ( d(h) .gt. 0. ) goto 140
\n if ( ia(h) .gt. 0 ) goto 150
\n h = -h
\n goto 160
\n140 if ( ia(h) .eq. 0 ) goto 160
\n h = -h
\n150 l = l – 1
\n160 call mod(u,ia,k,io,a,la,m,h,1,b)
\n if ( io .gt. 0 ) return
\n it = -1
\n170 do 180 i = 1,n
\n if ( iabs(ia(i)) .eq. 1 ) goto 180
\n x(i) = x(i) + s*d(i)
\n180 continue
\n if ( l .eq. 0 ) goto 390
\nc
\nc *** compute the projected gradient ***
\nc
\n190 do 200 i = 1,m
\n200 c(i) = 0.
\n do 220 j = 1,n
\n if ( iabs(ia(j)) .lt. 2 ) goto 220
\n if ( ia(j) .gt. 0 ) t = e*bl(j) + f*bu(j) – x(j)
\n if ( ia(j) .lt. 0 ) t = e*bu(j) + f*bl(j) – x(j)
\n do 210 i = 1,m
\n210 c(i) = c(i) + t*a(i,j)
\n220 continue
\n call sol(b,u,m,c)
\n it = it + 1
\n b1 = b2
\n b2 = 0.
\n r = 0.
\n s = 0.
\n do 250 j = 1,n
\n t = 0.
\n if ( ia(j) .eq. 2 ) t = x(j) – e*bl(j) – f*bu(j)
\n if ( ia(j) .eq.-2 ) t = x(j) – e*bu(j) – f*bl(j)
\n do 230 i = 1,m
\n230 t = t + b(i)*a(i,j)
\n if ( iabs(ia(j)) .eq. 1 ) goto 240
\n b2 = b2 + t*t
\n g(j) = t
\n r = dmax1(dabs(t),r)
\n goto 250
\n240 t = t*ia(j)
\n if ( t .le. s ) goto 250
\n s = t
\n h = j
\n250 continue
\n if ( s .gt. 10.*r ) goto 260
\n if ( b2 .eq. 0. ) goto 290
\n if ( it .ge. l ) goto 290
\n if ( it .ge. k ) goto 290
\n goto 90
\n260 x(h) = bl(h)
\n if ( ia(h) .gt. 0 ) x(h) = bu(h)
\n call mod(u,ia,k,io,a,la,m,h,0,b)
\n if ( io .gt. 0 ) return
\n call sol(b,u,m,c)
\n it = 0
\n b2 = 0.
\n do 280 j = 1,n
\n if ( iabs(ia(j)) .eq. 1 ) goto 280
\n t = 0.
\n if ( ia(j) .eq. 2 ) t = x(j) – e*bl(j) – f*bu(j)
\n if ( ia(j) .eq.-2 ) t = x(j) – e*bu(j) – f*bl(j)
\n do 270 i = 1,m
\n270 t = t + b(i)*a(i,j)
\n b2 = b2 + t*t
\n g(j) = t
\n280 continue
\n if ( b2 .eq. 0. ) goto 190
\n goto 90
\n290 if ( k .ge. l ) goto 300
\n io = 2
\n return
\n300 do 310 i = 1,m
\n310 b(i) = 0.
\n do 320 i = 1,n
\n if ( ia(i) .eq. 0 ) g(i) = x(i)
\n if ( ia(i) .eq. 1 ) g(i) = bu(i)
\n if ( ia(i) .eq.-1 ) g(i) = bl(i)
\n if ( iabs(ia(i)) .eq. 2 ) g(i) = x(i)
\n320 continue
\n do 340 j = 1,n
\n if ( iabs(ia(j)) .lt. 2 ) goto 340
\n call mod(u,ia,k,io,a,la,m,j,ia(j),b)
\n if ( io .gt. 0 ) return
\n if ( ia(j) .lt. 0 ) t = x(j) – bl(j)
\n if ( ia(j) .gt. 0 ) t = x(j) – bu(j)
\n do 330 i = 1,m
\n330 b(i) = b(i) + t*a(i,j)
\n340 continue
\n call sol(b,u,m,b)
\n do 380 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 380
\n t = x(j)
\n do 350 i = 1,m
\n350 t = t + b(i)*a(i,j)
\n x(j) = t
\n if ( ia(j) .ne. 0 ) goto 380
\n if ( t .lt. bl(j) ) goto 360
\n if ( t .le. bu(j) ) goto 380
\n360 io = 2
\n do 370 i = 1,n
\n370 x(i) = g(i)
\n return
\n380 continue
\n390 do 400 i = 1,m
\n400 c(i) = 0.
\n do 420 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 420
\n x(j) = x(j) – y(j)
\n t = x(j)
\n do 410 i = 1,m
\n410 c(i) = c(i) + t*a(i,j)
\n420 continue
\n430 call sol(b,u,m,c)
\n440 s = 1.
\n h = 0
\n do 480 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 480
\n t = 0.
\n do 450 i = 1,m
\n450 t = t + b(i)*a(i,j)
\n g(j) = t
\n r = t + y(j)
\n if ( r .lt. bu(j) ) goto 460
\n r = (bu(j)-x(j)-y(j))\/(t-x(j))
\n if ( r .ge. s ) goto 470
\n s = r
\n h = j
\n goto 470
\n460 if ( r .gt. bl(j) ) goto 470
\n r = (bl(j)-x(j)-y(j))\/(t-x(j))
\n if ( r .ge. s ) goto 470
\n s = r
\n h = -j
\n470 g(j) = t
\n480 continue
\n if ( h .eq. 0 ) goto 510
\n call mod(u,ia,k,io,a,la,m,h,1,b)
\n if ( io .gt. 0 ) return
\n if ( ia(h) .gt. 0 ) t = bu(h) – y(h)
\n if ( ia(h) .lt. 0 ) t = bl(h) – y(h)
\n do 490 i = 1,m
\n490 c(i) = c(i) – t*a(i,h)
\n call sol(b,u,m,c)
\n do 500 i = 1,n
\n500 if ( ia(i) .eq. 0 ) x(i) = x(i) + s*(g(i)-x(i))
\n goto 530
\n510 do 520 i = 1,n
\n520 if ( ia(i) .eq. 0 ) x(i) = g(i)
\n530 s = 0.
\n do 550 j = 1,n
\n if ( ia(j) .eq. 0 ) goto 550
\n if ( ia(j) .eq. 1 ) t = bu(j) – y(j)
\n if ( ia(j) .eq.-1 ) t = bl(j) – y(j)
\n do 540 i = 1,m
\n540 t = t – b(i)*a(i,j)
\n t = t*ia(j)
\n if ( t .le. s ) goto 550
\n s = t
\n l = j
\n550 continue
\n if ( s .eq. 0 ) goto 570
\n if ( ia(l) .eq. 1 ) t = bu(l) – y(l)
\n if ( ia(l) .eq.-1 ) t = bl(l) – y(l)
\n x(l) = t
\n call mod(u,ia,k,io,a,la,m,l,0,b)
\n if ( io .gt. 0 ) return
\n do 560 i = 1,m
\n560 c(i) = c(i) + t*a(i,l)
\n goto 430
\n570 if ( h .gt. 0 ) goto 440
\n do 580 i = 1,n
\n if ( ia(i) .eq. 0 ) x(i) = x(i) + y(i)
\n if ( ia(i) .eq. 1 ) x(i) = bu(i)
\n580 if ( ia(i) .eq.-1 ) x(i) = bl(i)
\n return
\n end
\nc
\nc in = 0 — all constraints are nonbinding initially
\nc in = 1 — binding constraints are given by values of ia
\nc in = 2 — binding constraints determined by x, bl, and bu
\nc in = 3 — use k,ia, and u as given
\nc
\nc LENGTH OF W AT LEAST 2M+2N IF NF1 INVOKED
\nC 4M+2N IF NF2 INVOKED
\nc
\n subroutine feas(x,u,ia,k,io,in,y,a,la,m,n,b,bl,bu,w)
\n integer ia(1),i,in,io,ip,j,k,l,la,m,n
\n real*8 a(la,1),b(1),bl(1),bu(1),u(1),x(1),y(1),w(1),t
\n ip = 0
\n if ( in .eq. 3 ) goto 100
\n k = n – m
\n do 10 i = 1,n
\n if ( bl(i) .lt. bu(i) ) goto 10
\n t = bl(i)
\n bl(i) = bu(i)
\n bu(i) = t
\n10 continue
\n if ( in .eq. 2 ) goto 50
\n if ( in .eq. 1 ) goto 30
\n do 20 i = 1,n
\n20 ia(i) = 0
\n goto 90
\n30 do 40 i = 1,n
\n40 if ( ia(i) .ne. 0 ) k = k – 1
\n if ( k .ge. 0 ) goto 90
\n goto 80
\n50 do 70 i = 1,n
\n ia(i) = 0
\n if ( y(i) .ne. bu(i) ) goto 60
\n ia(i) = 1
\n k = k – 1
\n goto 70
\n60 if ( y(i) .ne. bl(i) ) goto 70
\n ia(i) = -1
\n k = k – 1
\n70 continue
\n if ( k .ge. 0 ) goto 90
\n80 k = k + m
\n ip = 1
\n90 call mat(u,ia,a,la,m,n,ip)
\n call fac(u,0,m,io)
\n if ( io .gt. 0 ) return
\n100 if ( m .eq. 0 ) goto 130
\n do 110 i = 1,m
\n110 w(i) = b(i)
\n do 120 j = 1,n
\n t = y(j)
\n do 120 i = 1,m
\n120 w(i) = w(i) – t*a(i,j)
\n130 i = m + 1
\n j = i + m
\n if ( ip .eq. 1 ) goto 140
\n l = j + n
\n call nf1(x,u,ia,k,io,y,a,la,m,n,w,bl,bu,w(i),w(j),w(l))
\n if ( m .gt. 0 ) goto 160
\n return
\n140 do 150 l = 1,m
\n150 w(l+m) = 1.
\n l = j + m
\n 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))
\n160 if ( io .gt. 0 ) return
\n call lsq(x,u,ia,k,io,y,a,la,m,n,bl,bu,w,w(i),w(j))
\n return
\n end
\n double precision function fv(a,x,h,n,vl)
\n integer i,n
\n real*8 h(n,1),x(1),a,vl
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n fv = vl(h(1,2))
\n return
\n end
\nc
\n double precision function pv(a,x,h,f,g,d,nl,le,vlf,vlg,p2,n)
\n integer i,nl,n
\n real*8 d(1),g(1),h(n,1),le(1),x(1),a,f,p2,t,vlf
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n f = vlf(h(1,2))
\n call vlg(g,h(1,2))
\n t = f
\n do 20 i = 1,nl
\n20 t = t + le(i)*g(i) + p2*(g(i)-d(i))**2
\n pv = t
\n return
\n end
\nc
\n double precision function qv(a,x,h,f,g,nl,le,vlf,vlg,p2,n)
\n integer i,nl,n
\n real*8 g(1),h(n,1),le(1),x(1),a,f,p2,t,vlf
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n f = vlf(h(1,2))
\n call vlg(g,h(1,2))
\n t = f
\n do 20 i = 1,nl
\n20 t = t + (le(i)+p2*g(i))*g(i)
\n qv = t
\n return
\n end
\nc
\n double precision function fd(a,x,h,n,gr)
\n real*8 h(n,1),x(1),a,d
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n call gr(h(1,3),h(1,2))
\n d = 0.
\n do 20 i = 1,n
\n20 d = d + h(i,1)*h(i,3)
\n fd = d
\n return
\n end
\nc
\n double precision function pd(a,x,h,df,dg,g,d,nl,le,
\n 1 grf,grg,pn,n,lg,b)
\n integer i,nl,lg,n
\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
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n call grf(df,h(1,2))
\n call grg(dg,h(1,2))
\n do 20 i = 1,nl
\n20 b(i) = le(i) + pn*(g(i)-d(i))
\n s = 0.
\n do 40 j = 1,n
\n t = df(j)
\n do 30 i = 1,nl
\n30 t = t + b(i)*dg(i,j)
\n h(j,3) = t
\n s = s + t*h(j,1)
\n40 continue
\n pd = s
\n return
\n end
\nc
\n double precision function qd(a,x,h,df,dg,g,nl,le,
\n 1 grf,grg,pn,n,lg,b)
\n integer i,nl,lg,n
\n real*8 b(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n call grf(df,h(1,2))
\n call grg(dg,h(1,2))
\n do 20 i = 1,nl
\n20 b(i) = le(i) + pn*g(i)
\n s = 0.
\n do 40 j = 1,n
\n t = df(j)
\n do 30 i = 1,nl
\n30 t = t + b(i)*dg(i,j)
\n h(j,3) = t
\n s = s + t*h(j,1)
\n40 continue
\n qd = s
\n return
\n end
\nc
\n double precision function pd2(a,x,h,df,dg,g,d,nl,le,
\n 1 vlg,grf,grg,pn,n,lg,b)
\n integer i,nl,lg,n
\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
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n call grf(df,h(1,2))
\n call grg(dg,h(1,2))
\n call vlg(g,h(1,2))
\n do 20 i = 1,nl
\n20 b(i) = le(i) + pn*(g(i)-d(i))
\n s = 0.
\n do 40 j = 1,n
\n t = df(j)
\n do 30 i = 1,nl
\n30 t = t + b(i)*dg(i,j)
\n h(j,3) = t
\n s = s + t*h(j,1)
\n40 continue
\n pd2 = s
\n return
\n end
\nc
\n double precision function qd2(a,x,h,df,dg,g,nl,le,
\n 1 vlg,grf,grg,pn,n,lg,b)
\n integer i,nl,lg,n
\n real*8 b(1),df(1),dg(lg,1),g(1),h(n,1),le(1),x(1),a,pn,s,t
\n do 10 i = 1,n
\n10 h(i,2) = x(i) + a*h(i,1)
\n call grf(df,h(1,2))
\n call grg(dg,h(1,2))
\n call vlg(g,h(1,2))
\n do 20 i = 1,nl
\n20 b(i) = le(i) + pn*g(i)
\n s = 0.
\n do 40 j = 1,n
\n t = df(j)
\n do 30 i = 1,nl
\n30 t = t + b(i)*dg(i,j)
\n h(j,3) = t
\n s = s + t*h(j,1)
\n40 continue
\n qd2 = s
\n return
\n end
\n subroutine igs(s,f,a,b,c,fa,fb,fc,vf,vv,vg,gg,nl,j,y,z)
\n real*8 gg(1),vg(1),y(1),z(1),a,b,c,f,fa,fb,fc,s,vf,vv
\n integer j
\n j = j + 1
\n y(j) = s
\n z(j) = f
\n if ( f .le. fa ) goto 20
\n if ( f .le. fb ) goto 10
\n if ( f .gt. fc ) return
\n c = s
\n fc = f
\n return
\n10 c = b
\n b = s
\n fc = fb
\n fb = f
\n return
\n20 c = b
\n b = a
\n a = s
\n fc = fb
\n fb = fa
\n fa = f
\n vf = vv
\n do 30 i = 1,nl
\n30 vg(i) = gg(i)
\n return
\n end
\n subroutine ins(s,f,a,b,c,fa,fb,fc,j,y,z)
\n real*8 a,b,c,f,fa,fb,fc,s,y(1),z(1)
\n integer j
\n j = j + 1
\n y(j) = s
\n z(j) = f
\n if ( f .le. fa ) goto 20
\n if ( f .le. fb ) goto 10
\n if ( f .gt. fc ) return
\n c = s
\n fc = f
\n return
\n10 c = b
\n b = s
\n fc = fb
\n fb = f
\n return
\n20 c = b
\n b = a
\n a = s
\n fc = fb
\n fb = fa
\n fa = f
\n return
\n end
\n subroutine lsq(x,u,ia,k,io,y,a,la,m,n,bl,bu,b,c,g)
\n integer ia(1),h,i,io,j,k,la,m,n
\n double precision a(la,1),b(1),bl(1),bu(1),c(1),g(1),u(1)
\n double precision x(1),y(1),r,s,t
\n io = 0
\n if ( m .gt. 0 ) goto 30
\n do 20 i = 1,n
\n x(i) = y(i)
\n if ( x(i) .lt. bl(i) ) goto 10
\n if ( x(i) .le. bu(i) ) goto 20
\n x(i) = bu(i)
\n if ( ia(i) .eq. 0 ) k = k – 1
\n ia(i) = 1
\n goto 20
\n10 x(i) = bl(i)
\n if ( ia(i) .eq. 0 ) k = k – 1
\n ia(i) = -1
\n20 continue
\n return
\n30 do 40 i = 1,m
\n40 c(i) = 0.
\n do 60 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 60
\n x(j) = x(j) – y(j)
\n t = x(j)
\n do 50 i = 1,m
\n50 c(i) = c(i) + t*a(i,j)
\n60 continue
\n70 call sol(b,u,m,c)
\n80 s = 1.
\n h = 0
\n do 120 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 120
\n t = 0.
\n do 90 i = 1,m
\n90 t = t + b(i)*a(i,j)
\n g(j) = t
\n r = t + y(j)
\n if ( r .lt. bu(j) ) goto 100
\n r = (bu(j)-x(j)-y(j))\/(t-x(j))
\n if ( r .ge. s ) goto 110
\n s = r
\n h = j
\n goto 110
\n100 if ( r .gt. bl(j) ) goto 110
\n r = (bl(j)-x(j)-y(j))\/(t-x(j))
\n if ( r .ge. s ) goto 110
\n s = r
\n h = -j
\n110 g(j) = t
\n120 continue
\n if ( h .eq. 0 ) goto 150
\n call mod(u,ia,k,io,a,la,m,h,1,b)
\n if ( io .gt. 0 ) return
\n if ( ia(h) .gt. 0 ) t = bu(h) – y(h)
\n if ( ia(h) .lt. 0 ) t = bl(h) – y(h)
\n do 130 i = 1,m
\n130 c(i) = c(i) – t*a(i,h)
\n call sol(b,u,m,c)
\n do 140 i = 1,n
\n140 if ( ia(i) .eq. 0 ) x(i) = x(i) + s*(g(i)-x(i))
\n goto 80
\n150 do 160 i = 1,n
\n160 if ( ia(i) .eq. 0 ) x(i) = g(i)
\n170 s = 0.
\n do 190 j = 1,n
\n if ( ia(j) .eq. 0 ) goto 190
\n if ( ia(j) .eq. 1 ) t = bu(j) – y(j)
\n if ( ia(j) .eq.-1 ) t = bl(j) – y(j)
\n do 180 i = 1,m
\n180 t = t – b(i)*a(i,j)
\n t = t*ia(j)
\n if ( t .le. s ) goto 190
\n s = t
\n l = j
\n190 continue
\n if ( s .eq. 0 ) goto 210
\n if ( ia(l) .eq. 1 ) t = bu(l) – y(l)
\n if ( ia(l) .eq.-1 ) t = bl(l) – y(l)
\n x(l) = t
\n call mod(u,ia,k,io,a,la,m,l,0,b)
\n if ( io .gt. 0 ) return
\n do 200 i = 1,m
\n200 c(i) = c(i) + t*a(i,l)
\n goto 70
\n210 if ( h .gt. 0 ) goto 80
\n do 220 i = 1,n
\n if ( ia(i) .eq. 0 ) x(i) = x(i) + y(i)
\n if ( ia(i) .eq. 1 ) x(i) = bu(i)
\n220 if ( ia(i) .eq.-1 ) x(i) = bl(i)
\n return
\n end
\n subroutine mat(u,ia,a,la,m,n,in)
\n integer ia(1),i,in,j,k,l,la,m,n
\n real*8 a(la,1),u(1),t
\n if ( m .eq. 0 ) return
\n j = (m*m+m)\/2
\n do 10 i = 1,j
\n10 u(i) = 0.
\n if ( in .eq. 0 ) goto 30
\n l = 1
\n do 20 i = 1,m
\n u(l) = 1.
\n l = l + m – i + 1
\n20 continue
\n30 do 60 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 60
\n l = 0
\n do 50 i = 1,m
\n t = a(i,j)
\n do 40 k = i,m
\n l = l + 1
\n u(l) = u(l) + t*a(k,j)
\n40 continue
\n50 continue
\n60 continue
\n return
\n end
\n subroutine mob(v,io,u,ia,k,a,la,n,l,nl,ib,g,b,c)
\n integer ia(1),g,h,i,ib,io,j,k,l,la,n,nl
\n real*8 a(la,1),b(1),c(1),u(1),v(1),s,t
\n h = iabs(ib)
\n if ( l .gt. 0 ) goto 20
\n do 10 i = 1,nl
\n10 c(i) = a(i+l,h)
\n goto 110
\n20 do 30 i = 1,l
\n30 b(i) = a(i,h)
\n call sol(b,u,l,b)
\n t = 0.
\n do 40 i = 1,l
\n40 t = t – a(i,h)*b(i)
\n s = t + 1.
\n if ( s .gt. 0. ) goto 50
\n io = 1
\n return
\n50 do 60 i = 1,nl
\n60 c(i) = a(i+l,h)
\n do 90 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 90
\n t = 0.
\n do 70 i = 1,l
\n70 t = t – a(i,j)*b(i)
\n do 80 i = 1,nl
\n80 c(i) = c(i) + t*a(i+l,j)
\n90 continue
\n s = 1.\/dsqrt(s)
\n do 100 i = 1,nl
\n100 c(i) = s*c(i)
\n110 h = ib
\n call mod(v,ia,k,io,a,la,nl,h,g,c)
\n return
\n end
\n subroutine mod(u,ia,o,io,a,la,n,q,h,v)
\n integer ia(1),h,i,j,k,l,la,m,n,nm1,o,p,q
\n real*8 a(la,1),v(1),u(1),c,r,s,t
\n io = 0
\n p = h
\n if ( p .ne. 0 ) goto 10
\n if ( ia(q) .eq. 0 ) return
\n ia(q) = 0
\n o = o + 1
\n goto 60
\n10 if ( p .lt. 3 ) goto 40
\n if ( p .gt. 3 ) goto 30
\n do 20 i = 1,n
\n20 v(i) = 0.
\n v(q) = 1.
\n o = o – 1
\n goto 80
\n30 if ( p .eq. 5 ) p = 0
\n goto 80
\n40 if ( q .gt. 0 ) goto 50
\n p = -p
\n q = -q
\n50 i = ia(q)
\n if ( p .lt. 0 ) ia(q) = -1
\n if ( p .gt. 0 ) ia(q) = 1
\n if ( iabs(i) .eq. 1 ) return
\n o = o – 1
\n if ( o .ge. 0 ) goto 60
\n io = 1
\n return
\n60 if ( n .eq. 0 ) return
\n do 70 i = 1,n
\n70 v(i) = a(i,q)
\n80 k = 0
\n m = n
\n j = 1
\n90 t = v(j)\/u(j+k)
\n v(j) = t
\n if ( j .eq. n ) goto 110
\n j = j + 1
\n do 100 i = j,n
\n100 v(i) = v(i) – t*u(i+k)
\n m = m – 1
\n k = k + m
\n goto 90
\n110 r = v(n)
\n l = (n*n+n)\/2 + 1
\n nm1 = n – 1
\n do 170 k = 1,nm1
\n m = l – 2
\n l = m – k
\n s = r
\n j = n – k
\n c = v(j)
\n if ( dabs(c) .gt. dabs(s) ) goto 130
\n if ( s .eq. 0. ) goto 120
\n r = dabs(s)*dsqrt(1.+(c\/s)**2)
\n goto 140
\n120 c = 1.
\n s = 0.
\n goto 150
\n130 r = dabs(c)*dsqrt(1.+(s\/c)**2)
\n140 c = c\/r
\n s = s\/r
\nc
\nc process rows j and j+1 of U
\nc
\n150 v(j+1) = -s*u(l)
\n u(l) = c*u(l)
\n l = l + 1
\n do 160 i = l,m
\n t = u(i+k)
\n u(i+k) = c*t – s*u(i)
\n u(i) = s*t + c*u(i)
\n160 continue
\n170 continue
\n if ( p .ne. 0 ) r = 1 – r*r
\n if ( p .eq. 0 ) r = 1 + r*r
\n if ( r .gt. 0. ) goto 180
\n io = 1
\n return
\n180 r = dsqrt(r)
\n do 190 i = 1,n
\n190 u(i) = r*u(i)
\nc
\nc restore upper triangular form
\nc
\n m = 0
\n k = n
\n do 250 j = 2,n
\n l = m + 1
\n m = m + k
\n k = k – 1
\n c = u(l)
\n s = v(j)
\n if ( dabs(c) .gt. dabs(s) ) goto 210
\n if ( s .eq. 0. ) goto 200
\n r = dabs(s)*dsqrt(1.+(c\/s)**2)
\n goto 220
\n200 c = 1.
\n s = 0.
\n r = 0.
\n goto 230
\n210 r = dabs(c)*dsqrt(1.+(s\/c)**2)
\n220 c = c\/r
\n s = s\/r
\n230 u(l) = r
\n l = l + 1
\n do 240 i = l,m
\n t = u(i+k)
\n u(i+k) = c*t – s*u(i)
\n u(i) = s*t + c*u(i)
\n240 continue
\n250 continue
\n return
\n end
\n subroutine mul(le,li,e,u,ia,a,la,m,n,df,l,nl,z)
\n integer ia(1),i,j,l,la,m,n,nl
\n real*8 a(la,1),df(1),le(1),li(1),u(1),z(1),e,t
\n do 10 i = 1,m
\n10 le(i) = 0.
\n do 30 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 30
\n t = df(j)
\n do 20 i = 1,m
\n20 le(i) = le(i) – t*a(i,j)
\n30 continue
\n call sol(le,u,m,le)
\n e = 0.
\n do 60 j = 1,n
\n li(j) = 0.
\n t = df(j)
\n do 40 i = 1,m
\n40 t = t + le(i)*a(i,j)
\n if ( ia(j) .eq. 0 ) goto 50
\n if ( ia(j) .gt. 0 ) t = -t
\n li(j) = t
\n if ( t .ge. 0. ) goto 60
\n li(j) = 0.
\n50 e = dmax1(e,dabs(t))
\n60 continue
\n if ( nl .eq. 0 ) return
\n do 70 i = 1,nl
\n70 le(i+l) = le(i+l)*z(i)
\n return
\n end
\n subroutine mup(p,e,u,ia,le,a,la,n,dg,lg,pn,g,df,l,nl,y,z)
\n integer ia(1),i,j,l,la,lg,n,nl
\n real*8 a(la,1),dg(lg,1),df(1),g(1),le(1),p(1),u(1),y(1),z(1)
\n real*8 e,pn,t
\n do 10 i = 1,nl
\n10 z(i) = le(i) + pn*g(i)
\n if ( l .eq. 0 ) goto 90
\n do 20 i = 1,l
\n20 p(i) = 0.
\n do 50 j = 1,n
\n t = df(j)
\n do 30 i = 1,nl
\n30 t = t + z(i)*dg(i,j)
\n y(j) = t
\n if ( ia(j) .ne. 0 ) goto 50
\n do 40 i = 1,l
\n40 p(i) = p(i) – t*a(i,j)
\n50 continue
\n call sol(p,u,l,p)
\n e = 0.
\n do 80 j = 1,n
\n t = y(j)
\n do 60 i = 1,l
\n60 t = t + p(i)*a(i,j)
\n if ( ia(j) .eq. 0 ) goto 70
\n if ( ia(j) .gt. 0 ) t = -t
\n if ( t .ge. 0. ) goto 80
\n70 e = dmax1(e,dabs(t))
\n80 continue
\n return
\n90 e = 0.
\n do 120 j = 1,n
\n t = df(j)
\n do 100 i = 1,nl
\n100 t = t + z(i)*dg(i,j)
\n if ( ia(j) .eq. 0 ) goto 110
\n if ( ia(j) .gt. 0 ) t = -t
\n if ( t .ge. 0. ) goto 120
\n110 e = dmax1(e,dabs(t))
\n120 continue
\n return
\n end
\n subroutine nf1(x,u,ia,k,io,y,a,la,m,n,b,bl,bu,c,d,g)
\n integer ia(1),h,i,io,it,j,k,l,la,m,n
\n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),u(1)
\n double precision x(1),y(1),big,e,f,p,q,r,s,t
\nc
\nc set big = largest floating point number
\nc
\n big = 1.d50
\n io = 0
\n if ( m .gt. 0 ) goto 30
\n do 20 i = 1,n
\n x(i) = y(i)
\n if ( x(i) .lt. bl(i) ) goto 10
\n if ( x(i) .le. bu(i) ) goto 20
\n x(i) = bu(i)
\n if ( ia(i) .eq. 0 ) k = k – 1
\n ia(i) = 1
\n goto 20
\n10 x(i) = bl(i)
\n if ( ia(i) .eq. 0 ) k = k – 1
\n ia(i) = -1
\n20 continue
\n return
\n30 call sol(b,u,m,b)
\n do 50 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 50
\n t = y(j)
\n do 40 i = 1,m
\n40 t = t + b(i)*a(i,j)
\n x(j) = t
\n50 continue
\n l = 0
\n do 70 i = 1,n
\n d(i) = 0.
\n if ( iabs(ia(i)) .eq. 1 ) goto 70
\n if ( x(i) .lt. bl(i) ) goto 60
\n if ( x(i) .le. bu(i) ) goto 70
\n ia(i) = 2
\n l = l + 1
\n goto 70
\n60 ia(i) = -2
\n l = l + 1
\n70 continue
\n f = 0.
\n it = 0
\n if ( l .gt. 0 ) goto 230
\n io = -1
\n80 do 90 i = 1,n
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n90 continue
\n return
\n100 r = 0.
\n s = 0.
\n t = big
\n h = 0
\n it = it + 1
\n if ( it .eq. 1 ) q = 0.
\n if ( it .gt. 1 ) q = f\/e
\n do 190 i = 1,n
\n if ( iabs(ia(i)) .eq. 1 ) goto 190
\n d(i) = q*d(i) – g(i)
\n s = s + d(i)*g(i)
\n if ( ia(i) .ne. 0 ) goto 110
\n if ( d(i) .eq. 0. ) goto 190
\n if ( d(i) .lt. 0. ) goto 120
\n goto 170
\n110 r = r + d(i)**2
\n if ( ia(i) .gt. 0 ) goto 150
\n if ( x(i) .ge. bl(i) ) goto 140
\n if ( d(i) .le. 0. ) goto 190
\n120 p = (bl(i)-x(i))\/d(i)
\n130 if ( p .ge. t ) goto 190
\n t = p
\n h = -i
\n goto 190
\n140 p = 0.
\n goto 130
\n150 if ( x(i) .gt. bu(i) ) goto 160
\n p = 0.
\n goto 180
\n160 if ( d(i) .ge. 0. ) goto 190
\n170 p = (bu(i)-x(i))\/d(i)
\n180 if ( p .ge. t ) goto 190
\n t = p
\n h = i
\n190 continue
\n if ( r .ne. 0. ) s = -s\/r
\n j = iabs(h)
\n if ( t .le. s ) goto 200
\n if ( ia(j) .eq. 0 ) goto 210
\n200 s = t
\n if ( ia(j) .ne. 0 ) l = l – 1
\n call mod(u,ia,k,io,a,la,m,h,1,b)
\n if ( io .gt. 0 ) return
\n it = 0
\n210 do 220 i = 1,n
\n if ( iabs(ia(i)) .eq. 1 ) goto 220
\n x(i) = x(i) + s*d(i)
\n220 continue
\n if ( l .eq. 0 ) goto 80
\nc
\nc *** compute the projected gradient ***
\nc
\n230 do 240 i = 1,m
\n240 c(i) = 0.
\n do 260 j = 1,n
\n if ( iabs(ia(j)) .lt. 2 ) goto 260
\n if ( ia(j) .gt. 0 ) t = bu(j) – x(j)
\n if ( ia(j) .lt. 0 ) t = bl(j) – x(j)
\n do 250 i = 1,m
\n250 c(i) = c(i) + t*a(i,j)
\n260 continue
\n call sol(b,u,m,c)
\n e = f
\n f = 0.
\n r = 0.
\n s = 0.
\n do 290 j = 1,n
\n t = 0.
\n if ( ia(j) .eq. 2 ) t = x(j) – bu(j)
\n if ( ia(j) .eq.-2 ) t = x(j) – bl(j)
\n do 270 i = 1,m
\n270 t = t + b(i)*a(i,j)
\n if ( iabs(ia(j)) .eq. 1 ) goto 280
\n f = f + t*t
\n g(j) = t
\n r = dmax1(dabs(t),r)
\n goto 290
\n280 t = t*ia(j)
\n if ( t .le. s ) goto 290
\n s = t
\n h = j
\n290 continue
\n300 if ( s .gt. 10.*r ) goto 310
\n if ( r .eq. 0. ) goto 360
\n if ( it .ge. l ) goto 360
\n if ( it .ge. k ) goto 360
\n goto 100
\n310 x(h) = bl(h)
\n if ( ia(h) .gt. 0 ) x(h) = bu(h)
\n call mod(u,ia,k,io,a,la,m,h,0,b)
\n if ( io .gt. 0 ) return
\n call sol(b,u,m,c)
\n it = 0
\n f = 0.
\n r = 0.
\n do 330 j = 1,n
\n if ( iabs(ia(j)) .eq. 1 ) goto 330
\n t = 0.
\n if ( ia(j) .eq. 2 ) t = x(j) – bu(j)
\n if ( ia(j) .eq.-2 ) t = x(j) – bl(j)
\n do 320 i = 1,m
\n320 t = t + b(i)*a(i,j)
\n f = f + t*t
\n g(j) = t
\n r = dmax1(dabs(t),r)
\n330 continue
\n if ( r .gt. 0. ) goto 100
\n do 350 j = 1,n
\n if ( iabs(ia(j)) .ne. 1 ) goto 350
\n t = 0.
\n do 340 i = 1,m
\n340 t = t + b(i)*a(i,j)
\n t = t*ia(j)
\n if ( t .le. s ) goto 350
\n s = t
\n h = j
\n350 continue
\n goto 300
\n360 io = 2
\n return
\n end
\n subroutine nf2(x,u,ia,k,io,y,a,la,m,n,b,ib,bl,bu,c,d,g)
\n integer ia(1),h,i,io,it,j,k,l,la,m,n
\n double precision a(la,1),b(1),bl(1),bu(1),c(1),d(1),g(1),ib(1)
\n double precision u(1),x(1),y(1),big,e,f,p,q,r,s,t
\nc
\nc set big = largest floating point number
\nc
\n big = 1.d50
\n io = 0
\n do 10 i = 1,n
\n10 x(i) = y(i)
\n if ( m .eq. 0 ) goto 360
\n l = 0
\n do 20 i = 1,m
\n20 if ( ib(i) .ne. 0 ) l = l + 1
\n if ( l .eq. 0 ) goto 360
\n do 30 i = 1,m
\n if ( ib(i) .ne. 0. ) ib(i) = dsign(1.d0,b(i))
\n30 if ( ib(i) .eq. 0. ) b(i) = 0.
\n f = 0.
\n it = 0
\n goto 220
\n40 t = big
\n it = it + 1
\n if ( it .eq. 1 ) q = 0.
\n if ( it .gt. 1 ) q = f\/e
\n do 60 i = 1,n
\n d(i) = 0.
\n if ( ia(i) .ne. 0 ) goto 60
\n d(i) = q*d(i) – g(i)
\n if ( d(i) .gt. 0. ) goto 50
\n if ( d(i) .eq. 0. ) goto 60
\n p = (bl(i)-x(i))\/d(i)
\n if ( p .ge. t ) goto 60
\n t = p
\n h = -i
\n goto 60
\n50 p = (bu(i)-x(i))\/d(i)
\n if ( p .ge. t ) goto 60
\n t = p
\n h = i
\n60 continue
\n r = 0.
\n s = 0.
\n do 70 i = 1,m
\n c(i) = 0.
\n if ( ib(i) .eq. 0 ) goto 70
\n j = i + n
\n d(j) = q*d(j) – g(j)
\n r = r + d(j)**2
\n s = s + d(j)*b(i)
\n70 continue
\n if ( r .ne. 0. ) s = -s\/r
\n do 90 j = 1,n
\n r = d(j)
\n if ( r .eq. 0. ) goto 90
\n do 80 i = 1,m
\n80 c(i) = c(i) + r*a(i,j)
\n90 continue
\n r = big
\n do 130 i = 1,m
\n if ( ib(i) .eq. 0. ) goto 130
\n if ( ib(i) .gt. 0. ) goto 100
\n if ( b(i) .lt. 0. ) goto 110
\n p = 0.
\n goto 120
\n100 if ( b(i) .gt. 0. ) goto 110
\n p = 0.
\n goto 120
\n110 if ( c(i) .eq. 0. ) goto 130
\n p = b(i)\/c(i)
\n if ( p .le. 0. ) goto 130
\n120 if ( p .ge. r ) goto 130
\n r = p
\n j = i
\n130 continue
\n if ( t .lt. r ) goto 160
\n do 140 i = 1,n
\n140 x(i) = x(i) + r*d(i)
\n do 150 i = 1,m
\n150 if ( ib(i) .ne. 0. ) b(i) = b(i) – r*c(i)
\n ib(j) = 0.
\n b(j) = 0.
\n call mod(u,ia,k,io,a,la,m,j,3,c)
\n if ( io .gt. 0 ) return
\n l = l – 1
\n if ( l .eq. 0 ) goto 360
\n it = 0
\n goto 220
\n160 if ( t .le. s ) goto 190
\n do 170 i = 1,n
\n170 x(i) = x(i) + s*d(i)
\n do 180 i = 1,m
\n180 if ( ib(i) .ne. 0. ) b(i) = b(i) – s*c(i)
\n goto 220
\n190 do 200 i = 1,n
\n200 x(i) = x(i) + t*d(i)
\n do 210 i = 1,m
\n210 if ( ib(i) .ne. 0. ) b(i) = b(i) – t*c(i)
\n call mod(u,ia,k,io,a,la,m,h,1,c)
\n it = 0
\nc
\nc *** compute the projected gradient ***
\nc
\n220 call sol(c,u,m,b)
\n230 e = f
\n f = 0.
\n r = 0.
\n s = 0.
\n do 260 j = 1,n
\n t = 0.
\n do 240 i = 1,m
\n240 t = t – c(i)*a(i,j)
\n if ( iabs(ia(j)) .eq. 1 ) goto 250
\n f = f + t*t
\n g(j) = t
\n r = dmax1(dabs(t),r)
\n goto 260
\n250 t = t*ia(j)
\n if ( t .le. s ) goto 260
\n s = t
\n h = j
\n260 continue
\n270 if ( s .gt. 10.*r ) goto 300
\n if ( r .eq. 0. ) goto 350
\n if ( it .ge. l ) goto 350
\n if ( it .ge. k ) goto 350
\n280 do 290 i = 1,m
\n if ( ib(i) .eq. 0. ) goto 290
\n t = b(i) – c(i)
\n f = f + t*t
\n g(i+n) = t
\n290 continue
\n goto 40
\n300 x(h) = bl(h)
\n if ( ia(h) .gt. 0 ) x(h) = bu(h)
\n call mod(u,ia,k,io,a,la,m,h,0,c)
\n if ( io .gt. 0 ) return
\n call sol(c,u,m,b)
\n it = 0
\n f = 0.
\n r = 0.
\n do 320 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 320
\n t = 0.
\n do 310 i = 1,m
\n310 t = t – c(i)*a(i,j)
\n f = f + t*t
\n g(j) = t
\n r = dmax1(dabs(t),r)
\n320 continue
\n if ( r .gt. 0. ) goto 280
\n s = 0.
\n do 340 j = 1,n
\n if ( ia(j) .eq. 0 ) goto 340
\n t = 0.
\n do 330 i = 1,m
\n330 t = t – c(i)*a(i,j)
\n t = t*ia(j)
\n if ( t .le. s ) goto 340
\n s = t
\n h = j
\n340 continue
\n goto 270
\n350 io = 2
\n return
\n360 do 370 i = 1,n
\n if ( ia(i) .gt. 0 ) x(i) = bu(i)
\n if ( ia(i) .lt. 0 ) x(i) = bl(i)
\n370 continue
\n return
\n end
\n subroutine pre(p,e,u,ia,k,l,lm,io,g,a,la,m,n,b,c)
\n integer ia(1),i,io,j,k,l,la,lm,m,n
\n real*8 a(la,1),b(1),c(1),g(1),p(1),u(1),e,s,t
\n if ( m .gt. 0 ) goto 60
\n s = 0.
\n e = s
\n do 20 j = 1,n
\n t = g(j)
\n if ( ia(j) .eq. 0 ) goto 10
\n if ( ia(j) .lt. 0 ) t = -t
\n if ( t .le. s ) goto 20
\n s = t
\n l = j
\n goto 20
\n10 e = dmax1(e,dabs(t))
\n20 continue
\n if ( s .gt. 10.*e ) goto 30
\n l = 0
\n goto 40
\n30 i = ia(l)
\n call mod(u,ia,k,io,a,la,m,l,0,b)
\n if ( i .lt. 0 ) l = -l
\n lm = lm + 1
\n40 do 50 j = 1,n
\n p(j) = 0.
\n if ( ia(j) .eq. 0 ) p(j) = g(j)
\n50 continue
\n e = dmax1(e,s)
\n return
\n60 do 70 i = 1,m
\n70 c(i) = 0.
\n do 90 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 90
\n t = g(j)
\n do 80 i = 1,m
\n80 c(i) = c(i) – t*a(i,j)
\n90 continue
\n call sol(b,u,m,c)
\n s = 0.
\n e = s
\n do 120 j = 1,n
\n t = g(j)
\n do 100 i = 1,m
\n100 t = t + b(i)*a(i,j)
\n if ( ia(j) .eq. 0 ) goto 110
\n p(j) = 0.
\n if ( ia(j) .lt. 0 ) t = -t
\n if ( t .le. s ) goto 120
\n l = j
\n s = t
\n goto 120
\n110 p(j) = t
\n e = dmax1(e,dabs(t))
\n120 continue
\n if ( s .gt. 10.*e ) goto 130
\n l = 0
\n goto 190
\n130 t = g(l)
\n do 140 i = 1,m
\n140 c(i) = c(i) – t*a(i,l)
\n i = ia(l)
\n call mod(u,ia,k,io,a,la,m,l,0,b)
\n if ( i .lt. 0 ) l = -l
\n if ( io .gt. 0 ) return
\n lm = lm + 1
\n call sol(b,u,m,c)
\n150 do 180 j = 1,n
\n if ( ia(j) .eq. 0 ) goto 160
\n p(j) = 0.
\n goto 180
\n160 t = g(j)
\n do 170 i = 1,m
\n170 t = t + b(i)*a(i,j)
\n p(j) = t
\n180 continue
\n190 e = dmax1(e,s)
\n return
\n end
\n subroutine pro(p,u,ia,a,la,m,n,b)
\n integer ia(1),i,j,la,m,n
\n real*8 a(la,1),b(1),p(1),u(1),t
\n if ( m .gt. 0 ) goto 20
\n do 10 j = 1,n
\n if ( ia(j) .ne. 0 ) p(j) = 0.
\n10 continue
\n return
\n20 do 30 i = 1,m
\n30 b(i) = 0.
\n do 50 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 50
\n t = p(j)
\n do 40 i = 1,m
\n40 b(i) = b(i) – t*a(i,j)
\n50 continue
\n call sol(b,u,m,b)
\n do 80 j = 1,n
\n t = 0.
\n if ( ia(j) .ne. 0 ) goto 70
\n t = p(j)
\n do 60 i = 1,m
\n60 t = t + b(i)*a(i,j)
\n70 p(j) = t
\n80 continue
\n return
\n end
\n subroutine prp(p,e,u,v,ia,k,lm,io,g,a,la,l,nl,n,b,c)
\n integer ia(1),h,i,ib,io,j,k,l,la,lm,nl,n
\n real*8 a(la,1),b(1),c(1),g(1),p(1),u(1),v(1),e,t
\n call pre(p,e,u,ia,k,ib,lm,io,g,a,la,l,n,b,c)
\n h = iabs(ib)
\n if ( ib .eq. 0 ) goto 10
\n call mob(v,io,u,ia,k,a,la,n,l,nl,ib,5,b,c)
\n if ( io .gt. 0 ) return
\n10 do 20 i = 1,nl
\n20 c(i) = 0.
\n do 40 j = 1,n
\n t = p(j)
\n if ( t .eq. 0. ) goto 40
\n do 30 i = 1,nl
\n30 c(i) = c(i) + t*a(i+l,j)
\n40 continue
\n call sol(c,v,nl,c)
\n if ( l .eq. 0 ) goto 60
\n do 50 i = 1,l
\n50 b(i) = 0.
\n60 do 90 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 90
\n t = 0.
\n do 70 i = 1,nl
\n70 t = t + a(i+l,j)*c(i)
\n p(j) = p(j) – t
\n if ( l .eq. 0 ) goto 90
\n do 80 i = 1,l
\n80 b(i) = b(i) + t*a(i,j)
\n90 continue
\n if ( l .eq. 0 ) goto 120
\n call sol(b,u,l,b)
\n do 110 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 110
\n t = 0.
\n do 100 i = 1,l
\n100 t = t + b(i)*a(i,j)
\n p(j) = p(j) + t
\n110 continue
\n120 if ( h .eq. 0 ) return
\n if ( ib .lt. 0 ) goto 130
\n if ( p(h) .ge. 0. ) return
\n goto 140
\n130 if ( p(h) .le. 0. ) return
\n140 lm = lm – 1
\n call mob(v,io,u,ia,k,a,la,n,l,nl,ib,4,b,c)
\n if ( io .gt. 0 ) return
\n call mod(u,ia,k,io,a,la,l,ib,1,b)
\n h = 0
\n if ( l .gt. 0 ) goto 160
\n do 150 j = 1,n
\n p(j) = 0.
\n150 if ( ia(j) .eq. 0 ) p(j) = g(j)
\n goto 10
\n160 do 170 i = 1,l
\n170 b(i) = 0.
\n do 190 j = 1,n
\n if ( ia(j) .ne. 0 ) goto 190
\n t = g(j)
\n do 180 i = 1,l
\n180 b(i) = b(i) – t*a(i,j)
\n190 continue
\n call sol(b,u,l,b)
\n do 220 j = 1,n
\n t = 0.
\n if ( ia(j) .ne. 0 ) goto 210
\n t = g(j)
\n do 200 i = 1,l
\n200 t = t + b(i)*a(i,j)
\n210 p(j) = t
\n220 continue
\n goto 10
\n end
\n subroutine shk(a,m,n)
\n real*8 a(1)
\n integer g,h,i,j,k,l,m,n
\n if ( m .eq. 0 ) return
\n g = n – m
\n j = 0
\n l = 0
\n h = m
\n10 k = l + 1
\n l = l + h
\n h = h – 1
\n do 20 i = k,l
\n20 a(i) = a(i+j)
\n j = j + g
\n if ( k .lt. l ) goto 10
\n return
\n end
\n subroutine sol(x,a,n,b)
\n real*8 a(1),b(1),x(1),t
\n integer i,j,k,l,n
\n do 10 i = 1,n
\n10 x(i) = b(i)
\n l = 0
\n k = 1
\nc —————————–
\nc |*** forward elimination ***|
\nc —————————–
\n20 if ( k .eq. n ) goto 40
\n t = x(k)\/a(k+l)
\n x(k) = t
\n j = l
\n l = l + n – k
\n k = k + 1
\n if ( t .eq. 0. ) goto 20
\n do 30 i = k,n
\n30 x(i) = x(i) – t*a(i+j)
\n goto 20
\nc ———————————–
\nc |*** back substitution by rows ***|
\nc ———————————–
\n40 x(n) = x(n)\/a(k+l)**2
\n50 if ( k .eq. 1 ) return
\n j = k
\n k = k – 1
\n l = l + k – n
\n t = x(k)
\n do 60 i = j,n
\n60 t = t – x(i)*a(i+l)
\n x(k) = t\/a(k+l)
\n goto 50
\n end
\n subroutine xpn(a,m,n)
\n real*8 a(1)
\n integer g,h,i,j,k,l,m,n
\n if ( m .eq. 0 ) return
\n g = (m*m+m)\/2
\n h = g
\n l = n – m
\n j = l*m
\n k = 0
\n20 if ( g .le. 1 ) return
\n j = j – l
\n do 30 i = g,h
\n30 a(i+j) = a(i)
\n h = g – 1
\n k = k + 1
\n g = h – k
\n goto 20
\n end<\/p>\n\n\n\r\n\t\t\t<\/div>\r\n\t\t<\/div>\r\n\t<\/div>\r\n<\/section>\r\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1075,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"_acf_changed":false,"featured_post":"","footnotes":"","_links_to":"","_links_to_target":""},"class_list":["post-272","page","type-page","status-publish","hentry"],"acf":[],"_links":{"self":[{"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/pages\/272","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/users\/1075"}],"replies":[{"embeddable":true,"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/comments?post=272"}],"version-history":[{"count":2,"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/pages\/272\/revisions"}],"predecessor-version":[{"id":802,"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/pages\/272\/revisions\/802"}],"wp:attachment":[{"href":"https:\/\/people.clas.ufl.edu\/hager\/wp-json\/wp\/v2\/media?parent=272"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}