OPTPACK

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