************************************************************************
* SUBROUTINE PBEQU              ALL SYSTEMS                   05/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR SOLUTION OF DENSE SYSTEMS OF NONLINEAR
* EQUATIONS USING THE BROYDEN QUASI-NEWTON METHOD.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RO  AF(N)  VECTOR CONTAINING VALUES OF THE APPROXIMATED
*         FUNCTIONS.
*  II  IPAR(2)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*  RI  RPAR(7)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PBEQ.
*      RPAR(7)  TRUST REGION STEPSIZE.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN
*         SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NREM  NUMBER OF CONSTRAINT DELETIONS.
*  IO  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PBEQ  SOLUTION OF DENSE NONLINEAR SYSTEMS OF EQUATIONS BY THE
*         BROYDEN QUASI-NEWTON METHOD.
*
* EXTERNAL SUBROUTINES :
*  SE  FUN  COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL FUN(N,KA,X,FA) WHERE N IS THE NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(N) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE
*         APPROXIMATED FUNCTION.
*
      SUBROUTINE PBEQU(N,X,AF,IPAR,RPAR,F,IPRNT,ITERM)
      DOUBLE PRECISION F
      INTEGER          IPRNT,ITERM,N
      DOUBLE PRECISION AF(*),RPAR(7),X(*)
      INTEGER          IPAR(2)
      INTEGER          LAFO,LAG,LG,LGA,LGO,LH,LXO,LS
      INTEGER          NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON           /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      DOUBLE PRECISION RA(:)
      ALLOCATABLE RA
      ALLOCATE (RA((N+6)*N+N*(N+1)/2))
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGA = 1
      LAG = LGA + N
      LG = LAG + N * N
      LS = LG + N
      LXO = LS + N
      LGO = LXO + N
      LAFO = LGO + N
      LH = LAFO + N
      CALL PBEQ(N,X,RA(LGA),RA(LAG),RA(LG),RA(LS),RA(LXO),RA(LGO),
     +          RA(LH),AF,RA(LAFO),RPAR(1),RPAR(2),RPAR(3),RPAR(4),
     +          RPAR(5),RPAR(7),F,IPAR(1),IPAR(2),IPRNT,ITERM)
      DEALLOCATE(RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PBEQ                   ALL SYSTEMS                95/12/01
* PORTABILITY : ALL SYSTEMS
* 95/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF SPARSE NONLINEAR SYSTEMS OF EQUATIONS BY THE NEWTON
* METHOD USING THE PRECONDITIONED SMOOTHED CGS SUBALGORITHM FOR
* ITERATIVE SOLUTION OF LINEARIZED SYSTEMS.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RA  GA(N)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RA  AG(N*N)  RECTANGULAR MATRIX WHICH IS USED FOR THE DIRECTION
*        VECTOR DETERMINATION.
*  RA  G(N)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(N)  DIRECTION VECTOR.
*  RA  XO(N)  AUXILIARY VECTOR.
*  RA  GO(N)  AUXILIARY VECTOR.
*  RA  H(N*(N+1)/2)  AUXILIARY MATRIX.
*  RO  AF(N)  VECTOR WHOSE ELEMENTS ARE VALUES OF THE APPROXIMATED
*         FUNCTIONS.
*  RA  AFO(N)  AUXILIARY VECTOR.
*  RI  XMAX  MAXIMUM STEPSIZE.
*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
*  RI  TOLB  TOLERANCE FOR THE FUNCTION FALUE.
*  RI  TOLG  TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RI  XDEL  TRUST REGION STEPSIZE.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  MIT  MAXIMUN NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUN NUMBER OF FUNCTION EVALUATIONS.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN
*         SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE.
*  IO  ITERM  CAUSE OF TERMINATION.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NREM  NUMBER OF CONSTRAINT DELETIONS.
*  IO  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PA1SQ1  COMPUTATION OF THE VALUE AND THE GRADIENT OF THE
*         OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES
*         OF THE APPROXIMATED FUNCTIONS.
*  S   PNSTEP  VECTOR ON THE BOUNDARY OF THE TRUST REGION.
*  S   PS0G01  TRUST REGION STEPSIZE USING ONLY FUNCTION VALUES.
*  S   PYFUT1  TEST ON TERMINATION.
*  S   PYTRUD  COMPUTATION OF THE VARIABLE VECTOR AND THE GRADIENT
*         DIFFERENCES.
*  S   MXDCMM  MATRIX-VECTOR PRODUCT. RECTANGULAR MATRIX IS STORED
*         COLUMNWISE.
*  S   MXDPRB  BACK SUBSTITUTION USING THE CHOLESKI DECOMPOSITION
*         OBTAINED BY MXDPRF.
*  S   MXDPRC  CORRECTION OF THE CHOLESKI DECOMPOSITION OBTAINED BY
*         MXDPRF.
*  S   MXDPRM  MATRIX-VECTOR PRODUCT USING THE CHOLESKI DECOMPOSITION
*         OBTAINED BY MXDPRF.
*  S   MXDRMM  MATRIX-VECTOR PRODUCT. RECTANGULAR MATRIX IS STORED
*         ROWWISE.
*  S   MXDRQF  QR DECOMPOSITION OF ROWWISE STORED RECTANGULAR MATRIX
*         USING HOUSEHOLDER TRANSFORMATIONS WITHOUT PIVOTING.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
* EXTERNAL SUBROUTINES :
*  SE  FUN  COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
*         CALLING SEQUENCE: CALL FUN(N,KA,X,FA) WHERE N IS A NUMBER
*         OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION,
*         X(N) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE
*         APPROXIMATED FUNCTION.
*
* METHOD :
* BROYDEN GOOD QUASI-NEWTON METHOD.
*
      SUBROUTINE PBEQ(N,X,GA,AG,G,S,XO,GO,H,AF,AFO,XMAX,TOLX,TOLF,
     +                TOLB,TOLG,XDEL,F,MIT,MFV,IPRNT,ITERM)
      DOUBLE PRECISION F,GMAX,TOLB,TOLF,TOLG,TOLX,XMAX,XDEL
      INTEGER          IPRNT,ITERM,MFV,MIT,N
      DOUBLE PRECISION AF(*),AFO(*),AG(*),G(*),GA(*),GO(*),H(*),S(*),
     +                 X(*),XO(*)
      DOUBLE PRECISION BET1,BET2,GAM1,GAM2,DMAX,EPS4,EPS5,ETA0,ETA2,
     +                 ETA9,FMIN,FO,FP,GNORM,P,PO,PP,R,RMAX,RO,SNORM,
     +                 B1,B2,B3,D3,S1,XDELO,UMAX,T,MXVDOT
      INTEGER          I,INF,IREST,ITERD,ITERS,NRED,IPOM1,IPOM2,ITERH,
     +                 KD,KIT,LD,MRED,MTESF,MTESX,NTESF,NTESX,LDS,
     +                 IDEC,INITS,KTERS,IEST,ITES,MAXST,IRES1,IRES2,
     +                 ISYS,MFG,INITD,MES1,MES2,MES3,MOT3,KDA,IDIR
      INTEGER          NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON           /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (ABS(IPRNT).GT.1) WRITE (6,'(1X,''ENTRY TO PBEQ :'')')
*
*     INITIATION
*
      NRES=0
      NDEC=0
      NREM=0
      NADD=0
      NIT=0
      NFV=0
      NFG=0
      NFH=0
      ISYS = 0
      IEST = 1
      ITES = 1
      MTESX = 2
      MTESF = 2
      INITS = 1
      INITD = 1
      ITERM = 0
      ITERD = 0
      ITERS = 2
      ITERH = 0
      KTERS = 0
      IREST = 1
      IRES1=999
      IRES2 = 0
      MRED = 5
      IDEC = 0
      IDIR = 0
      EPS4=0.10D 0
      EPS5=0.90D 0
      ETA0 = 1.0D-15
      ETA2 = 1.0D-18
      ETA9 = 1.0D 60
      BET1=0.05D 0
      BET2=0.75D 0
      GAM1=2.0D 0
      GAM2=1.0D 6
      FMIN = 0.0D 0
      IF (TOLX.LE.0.0D0) TOLX = 1.0D-16
      IF (TOLF.LE.0.0D0) TOLF = 1.0D-16
      IF (TOLB.LE.0.0D0) TOLB = 1.0D-16
      IF (TOLG.LE.0.0D0) TOLG = 1.0D-16
      IF (XMAX.LE.0.0D0) XMAX = 1.0D 16
      IF (MIT.LE.0) MIT = 1000
      IF (MFV.LE.0) MFV = 3000
      MFG = 2*MIT
      MES1 = 3
      MES2 = 2
      MES3 = 1
      MOT3 = 0
      KDA = 0
      KD =  0
      LD = -1
      KIT = 0
      FO = FMIN
      GMAX = ETA9
      DMAX = ETA9
      T=0.0D 0
      FO=FMIN
      IPOM1=0
      IPOM2=1
*
*     -----------------
*     MODEL DESCRIPTION
*     -----------------
*
      CALL PA1SQ1(N,X,F,AF,GA,AG,G,ETA0,KDA,KD,LD,NFV,NFG)
*
*     ------------------------
*     END OF MODEL DESCRIPTION
*     ------------------------
*
11020 CONTINUE
      IF (ABS(IPRNT).GT.1)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F ='', G16.9)') NIT,NFV,NFG,F
      CALL PYFUT1(N,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
     & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1,
     & IRES2,IREST,ITERS,ITERM)
      IF(ITERM.NE.0) GO TO 11090
11030 CONTINUE
      IF (IREST.LE.0) GO TO 12420
      KD=1
      CALL PA1SQ1(N,X,F,AF,GA,AG,G,ETA0,KDA,KD,LD,NFV,NFG)
      IDEC=0
      IF (IDIR.EQ.0) THEN
      IF (KIT.LT.NIT) THEN
      NRES=NRES+1
      KIT = NIT
      ELSE
      ITERM=-10
      IF (ITERS.LT.0) ITERM=ITERS-5
      IF (GMAX.LE.1.0D 2*TOLG) ITERM=-ITERM
      ENDIF
      ENDIF
12420 CONTINUE
      IF(ITERM.NE.0) GO TO 11090
*
*     -----------------------
*     DIRECTION DETERMINATION
*     TEMPLATE : UDGGA1
*     -----------------------
*
      IF (IDEC.LT.0) THEN
      IDEC=1
      INF=0
      ENDIF
      IF (IDEC.EQ.0) THEN
      ELSEIF(IDEC.EQ.1) THEN
      ELSE
      ITERD=-1
      GO TO 12530
      ENDIF
      IF (LD.LE.0) THEN
      CALL MXDCMM(N,N,AG,AF,G)
      IF (IDEC.EQ.1) THEN
      CALL MXVCOP(N,G,GO)
      CALL MXDPRM(N,H,G,-1)
      ENDIF
      ELSE IF (IDEC.EQ.1) THEN
      IF (IPOM1.EQ.1) THEN
      CALL MXVCOP(N,G,GO)
      CALL MXDPRB(N,H,GO,1)
      ELSE
      CALL MXDCMM(N,N,AG,AF,GO)
      ENDIF
      ENDIF
      B2=MXVDOT(N,G,G)
      GNORM=SQRT(B2)
      IF (KD.LE.0) THEN
      IF (GNORM.LE.TOLG) THEN
      ITERM=4
      GO TO 11090
      ENDIF
      ELSE
      IF (GNORM.LE.0.0D 0) THEN
      ITERM=4
      GO TO 11090
      ENDIF
      ENDIF
      INITD=MAX(ABS(INITD),1)
      IF (IDEC.EQ.0) THEN
      CALL MXDRMM(N,N,AG,G,AFO)
      B1=MXVDOT(N,AFO,AFO)
      ELSE
      CALL MXVCOP(N,G,S)
      CALL MXDPRM(N,H,S,1)
      B1=MXVDOT(N,S,S)
      ENDIF
      IF (XDEL.LE.0.0D 0) THEN
*
*     INITIAL TRUST REGION BOUND
*
      IF (B1.LE.0.0D 0) THEN
      XDEL=GNORM
      ELSE
      XDEL=(B2/B1)*GNORM
      ENDIF
      IF(IEST.EQ.1) XDEL=MIN(XDEL,4.0D0*(F-FMIN)/GNORM)
      XDEL=MIN(XDEL,XMAX)
      ENDIF
      IF (B1.LE.0.0D0.OR.B2*GNORM.GE.B1*XDEL) THEN
*
*     SCALED STEEPEST DESCENT DIRECTION IS ACCEPTED
*
      CALL MXVSCL(N,-XDEL/GNORM,G,S)
      SNORM=XDEL
      ITERD=3
      GO TO 12520
      ENDIF
      IF (IDEC.EQ.0) THEN
*
*     QR DECOMPOSITION
*
      S1=ETA2
      CALL MXDRQF(N,N,AG,H)
      CALL MXDPRC(N,H,INF,S1)
      IF (IPOM1.EQ.1) THEN
      CALL MXVCOP(N,G,GO)
      CALL MXDPRB(N,H,GO,1)
      ELSE
      CALL MXDCMM(N,N,AG,AF,GO)
      ENDIF
      NDEC=NDEC+1
      IDEC=1
      ENDIF
*
*     COMPUTATION OF THE NEWTON DIRECTION
*
      CALL MXDPRB(N,H,GO,-1)
      D3=SQRT(MXVDOT(N,GO,GO))
*
*     COMPUTATION OF THE STEEPEST DESCENT DIRECTION
*
      B2=B2/B1
      SNORM=B2*GNORM
      CALL MXVSCL(N,-B2,G,S)
      CALL MXVNEG(N,GO,GO)
      CALL MXVDIF(N,GO,S,XO)
      B1=MXVDOT(N,S,XO)
      B2=MXVDOT(N,XO,XO)
      IF (B2.LE.1.0D-8*XDEL*XDEL) THEN
*
*     NEWTON AND THE STEEPEST DESCENT DIRECTION ARE
*     APPROXIMATELY EQUAL
*
      CALL MXVCOP(N,GO,S)
      SNORM=D3
      ITERD=1
      ELSE IF (B1.LE.0.0D 0) THEN
*
*     BOUNDARY STEP WITH NEGATIVE INCREMENT
*
      CALL PNSTEP(XDEL,SNORM,-B1,B2,B3)
      CALL MXVDIR(N,-B3,XO,S,S)
      SNORM=XDEL
      ITERD=3
      ELSE IF (D3.LE.XDEL) THEN
*
*     NEWTON DIRECTION IS ACCEPTED
*
      CALL MXVCOP(N,GO,S)
      SNORM=D3
      ITERD=1
      ELSE
*
*     DOUBLE DOGLEG STRATEGY
*
      D3=XDEL/D3
      B3=MXVDOT(N,S,GO)
      D3=MAX(D3,SNORM*SNORM/B3)
      CALL MXVDIR(N,-D3,GO,S,XO)
      B1=SNORM*SNORM-D3*B3
      B2=MXVDOT(N,XO,XO)
      CALL PNSTEP(XDEL,SNORM,-B1,B2,B3)
      CALL MXVDIR(N,-B3,XO,S,S)
      SNORM=XDEL
      ITERD=3
      ENDIF
12520 CONTINUE
      IF (IDEC.EQ.0) THEN
      CALL MXDRMM(N,N,AG,S,AFO)
      PP=MXVDOT(N,AFO,AFO)*0.5D0
      ELSE
      CALL MXVCOP(N,S,GO)
      CALL MXDPRM(N,H,GO,1)
      PP=MXVDOT(N,GO,GO)*0.5D0
      ENDIF
12530 CONTINUE
*
*     ------------------------------
*     END OF DIRECTION DETERMINATION
*     ------------------------------
*
      IF(ITERD.LT.0) THEN
        ITERM=ITERD
      ENDIF
      P=MXVDOT(N,S,G)
      IF (INITD.LT.0.AND.P.GT.0.0D 0) THEN
*
*     CHANGE OF THE SIGN
*
        CALL MXVNEG(N,S,S)
        P=-P
      ENDIF
*
*     TEST ON LOCALLY CONSTRAINED STEP
*
      IF(SNORM.LE.0.0D 0) THEN
        IREST=MAX(IREST,1)
      ELSE
        IREST=0
      ENDIF
      IF (IREST.EQ.0) THEN
*
*     PREPARATION OF SIMPLE SEARCH
*
        RMAX=XMAX/SNORM
      ENDIF
      IF(ITERM.NE.0) GO TO 11090
      IF(IREST.NE.0) GO TO 11030
      LDS=LD
11050 CONTINUE
      FP = FO
      RO = 0.0D 0
      FO = F
      PO = P
      CALL MXVCOP(N,X,XO)
      CALL MXVCOP(N,G,GO)
      CALL MXVCOP(N,AF,AFO)
11060 CONTINUE
      CALL PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1,BET2,
     & GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED,
     & KTERS,MES1,MES2,MES3,ISYS)
      IF (ISYS.EQ.0) GO TO 11064
C      GO TO (11064,11062) ISYS+1
11062 CONTINUE
      CALL MXVDIR(N,R,S,XO,X)
      CALL PA1SQ1(N,X,F,AF,GA,AG,G,ETA0,KDA,KD,LD,NFV,NFG)
      GO TO 11060
11064 CONTINUE
      IF (ITERS.LE.0) THEN
      IF (IDIR.EQ.0.OR.MOT3.LE.0) THEN
      R=0.0D 0
      F=FO
      P=PO
      CALL MXVCOP(N,XO,X)
      ENDIF
      IF (MOT3.LE.0) CALL MXVCOP(N,AFO,AF)
      IF(IDIR.EQ.0) THEN
      IREST=MAX(IREST,1)
      LD=LDS
      GO TO 11030
      ELSE IF (MOT3.EQ.0) THEN
      LD=LDS
      GO TO 11030
      ELSE
      ITERD=0
      ENDIF
      ENDIF
      IF (IPOM2.GT.0) THEN
      CALL PUDRV1(R,FO,F,PO,IPOM1,IPOM2,NRED,IREST)
      ENDIF
      IF (IPOM1.EQ.1) KD=1
      IF(KD.GT.LD) THEN
      CALL PA1SQ1(N,X,F,AF,GA,AG,G,ETA0,KDA,KD,LD,NFV,NFG)
      ENDIF
      KD=0
      IF(IPOM1.EQ.2) THEN
      CALL PYTRUF(N,N,X,XO,AF,AFO,R,F,FO,P,PO,DMAX,KD,LD,ITERS)
      CALL PUDBQ1(N,N,H,ETA2,AG,S,XO,AFO,1,ITERH,IDEC,NDEC)
      ELSE
      CALL PYTRUD(N,X,XO,G,GO,R,F,FO,P,PO,DMAX,KD,LD,ITERS)
      ENDIF
      IF(IDIR.EQ.0) THEN
      IF(ITERH.NE.0) THEN
      IREST=MAX(IREST,1)
      ENDIF
      GO TO 11020
      ELSE
      GO TO 11030
      ENDIF
11090 CONTINUE
11100 CONTINUE
      IF (IPRNT.GT.1.OR.IPRNT.LT.0)
     & WRITE(6,'(1X,''EXIT FROM PBEQ :'')')
      IF (IPRNT.NE.0)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F ='', G16.9,2X,''ITERM='',I3)') NIT,NFV,NFG,
     & F,ITERM
      IF (IPRNT.LT.0)
     & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
     & (X(I),I=1,N)
13599 CONTINUE
*
*     -----------------
*     END OF METHOD (1)
*     -----------------
*
10399 CONTINUE
      END