117 lines
4.2 KiB
FortranFixed
117 lines
4.2 KiB
FortranFixed
subroutine cpnonsyssolver(funcnleq1,fmin_funcnleq1,
|
|
& f1dim_funcnleq1,x0min,x0ori,xp,x0max,fp,
|
|
& nunknowns,iwhichsolver)
|
|
implicit none
|
|
integer nunknowns,iwhichsolver
|
|
double precision x0min(nunknowns),x0ori(nunknowns),
|
|
& xp(nunknowns),x0max(nunknowns),fp(nunknowns)
|
|
!-------- Specified values ---------------------------------------
|
|
!funcnleq1: the subroutine that calculates the functional values of the
|
|
! the nonlinear system in the following form:
|
|
! funcnleq1(nunknowns,xp,fp,fsqsum)
|
|
!fmin_funcnleq1: the subroutine that calls funcnleq1 and returns fsqsum (half
|
|
! of the sum of the squared functional values of the nonlinear system)
|
|
! fmin_funcnleq1(nunknowns,xp,fsqsum)
|
|
!f1dim_funcnleq1: a function subroutine that returns fsqsum
|
|
! f1dim_funcnleq1(xp)
|
|
! nunknowns: The number of unknowns to be solved
|
|
! x0ori(1:nunknowns): initial guess for the unknowns
|
|
! x0min(1:nunknowns): lower bound of the solution
|
|
! x0max(1:nunknowns): upper bound of the solution
|
|
! --------- Calculated values -------------------------------------
|
|
! fp(1:nunknowns): function values at the last step of iteration
|
|
! xp(1:nunknowns): final solutions
|
|
! iwhichsolver:
|
|
! =1 solved by plain fixed point method 1
|
|
! =2 solved by fixed point method 2
|
|
! =3 solved by fixed point method 3
|
|
! =4 solved by fixed point method 4
|
|
! =6 solved by broydn
|
|
! =7 Solved by multiobjective minimization.
|
|
! =-9999 Best approximation returned. Solution may not be accurate.
|
|
! --------- Local variables ---------------------------------------
|
|
double precision x0(nunknowns),TOLF,stpmax,scldstpmax,
|
|
& sum,tb,tp,xb(nunknowns),fb(nunknowns),fsqsum,
|
|
& f1dim_funcnleq1
|
|
integer i,irepeat,maxrepeats,IERR,notfound
|
|
intrinsic dble
|
|
parameter(maxrepeats=100,notfound=-9999,TOLF=1.0d-10)
|
|
external funcnleq1,fmin_funcnleq1,f1dim_funcnleq1
|
|
!-------------------------------------------------------------------
|
|
stpmax=0.0d0
|
|
sum=0.0d0
|
|
do i=1, nunknowns
|
|
x0(i)=x0ori(i)
|
|
sum=sum+x0ori(i)*x0ori(i)
|
|
stpmax=stpmax+
|
|
& (x0min(i)-x0max(i))*(x0min(i)-x0max(i))
|
|
enddo
|
|
stpmax=dsqrt(stpmax)/4.0d0
|
|
scldstpmax=stpmax/dmax1(dsqrt(sum),dble(nunknowns))
|
|
! In Numerical Recipes, scldstpmax (STPMX) is 100
|
|
scldstpmax=dmax1(100.0d0,scldstpmax)
|
|
iwhichsolver=notfound
|
|
do irepeat=1,maxrepeats
|
|
call cpfixedpoint(funcnleq1,x0min,x0,xp,
|
|
& x0max,fp,nunknowns,TOLF,stpmax,iwhichsolver)
|
|
if(iwhichsolver.ne.notfound)return
|
|
tp=dabs(fp(1))
|
|
xb(1)=xp(1)
|
|
do i=2,nunknowns
|
|
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
|
|
xb(i)=xp(i)
|
|
enddo
|
|
call cpbroydn(x0min,xb,x0max,scldstpmax,nunknowns,
|
|
& fb,funcnleq1,TOLF,IERR)
|
|
call funcnleq1(nunknowns,xb,fb,fsqsum)
|
|
tb=dabs(fb(1))
|
|
do i=2,nunknowns
|
|
if(dabs(fb(i)).gt.tb)tb=dabs(fb(i))
|
|
enddo
|
|
do i=1,nunknowns
|
|
if(xb(i).lt.x0min(i).or.xb(i).gt.x0max(i))then
|
|
tb=1.0d+100
|
|
endif
|
|
enddo
|
|
if(tb.lt.tp)then
|
|
do i=1,nunknowns
|
|
xp(i)=xb(i)
|
|
fp(i)=fb(i)
|
|
enddo
|
|
if(tb.lt.TOLF)then
|
|
iwhichsolver=6
|
|
return
|
|
endif
|
|
endif
|
|
fsqsum=0.0d0
|
|
do i=1,nunknowns
|
|
fsqsum=fsqsum+fp(i)*fp(i)
|
|
enddo
|
|
tp=fsqsum
|
|
call cpnongradopt(nunknowns,fmin_funcnleq1,
|
|
& f1dim_funcnleq1,xp,x0min,x0max,TOLF,fsqsum)
|
|
if(dabs(tp-fsqsum).gt.TOLF)then
|
|
call cpRepeatCompassSearch(nunknowns,xp,fsqsum,
|
|
& x0min,x0max,fmin_funcnleq1,f1dim_funcnleq1,
|
|
& TOLF)
|
|
endif
|
|
call funcnleq1(nunknowns,xp,fp,fsqsum)
|
|
tp=dabs(fp(1))
|
|
do i=2,nunknowns
|
|
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
|
|
enddo
|
|
if(tp.lt.TOLF)then
|
|
iwhichsolver=7
|
|
return
|
|
endif
|
|
IERR=0
|
|
do i=1,nunknowns
|
|
if(dabs(xp(i)-x0(i)).gt.TOLF)IERR=1
|
|
enddo
|
|
if(IERR.eq.0)return
|
|
do i=1,nunknowns
|
|
x0(i)=xp(i)
|
|
enddo
|
|
enddo
|
|
end subroutine cpnonsyssolver
|