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