SUBROUTINE HFUN1P(QDG,LAMBDA,X,
     $ PDG,CL,COEF,RHO,
     $ DRHOX,DRHOL,XDGM1,XDG,
     $ G,DG,PXDGM1,PXDG,
     $ F,DF,XX,TRM,
     $ DTRM,CLX,DXNP1,
     $ N,MMAXT,IDEG,
     $ NUMT,KDEG)
C
C  HFUN1P  EVALUATES THE CONTINUATION EQUATION "RHO".
C
C  NOTE THAT:
C    DRHOX IS THE "REALIFICATION" OF DCRHOX, WHERE
C    DCRHOX DENOTES THE (COMPLEX) PARTIAL
C    DERIVATIVE MATRIX OF THE CONTINUATION SYSTEM
C    WITH RESPECT TO X,  AND
C    DRHOL IS THE "REALIFICATION" OF DCRHOL, WHERE
C    DCRHOL DENOTES THE (COMPLEX) PARTIAL
C    DERIVATIVE MATRIX OF THE CONTINUATION SYSTEM
C    WITH RESPECT TO LAMBDA. THUS
C      DRHOX(2J-1,2K-1) = DCRHOX(1,J,K)
C      DRHOX(2J  ,2K  ) = DCRHOX(1,J,K)
C      DRHOX(2J-1,2K  ) =-DCRHOX(2,J,K)
C      DRHOX(2J  ,2K-1) = DCRHOX(2,J,K)
C      DRHOL(2J-1,N2P1) = DCRHOL(1,J)
C      DRHOL(2J  ,N2P1) = DCRHOL(2,J)
C       RHO(2J-1)      = CRHO(1,J)
C       RHO(2J  )      = CRHO(2,J)
C    WHERE CRHO DENOTES THE (COMPLEX) CONTINUATION SYSTEM,
C    THE INITIAL "1" OR "2" DENOTES REAL OR IMAGINARY PARTS,
C    RESPECTIVELY, "J" INDEXES THE EQUATION, "K" INDEXES THE PARTIAL
C    DERIVATIVE, AND NEITHER DCRHOX NOR DCRHOL ARE PROGRAM VARIABLES.
C
C  ON INPUT:
C
C    QDG  IS THE "RANDOM" PARAMETER "A".
C
C    LAMBDA  IS THE CONTINUATION PARAMETER.
C
C    X    IS THE INDEPENDENT VARIABLE.
C
C    PDG  IS ONE OF THE PARAMETERS THAT DEFINES G (SEE SUBROUTINE
C         GFUNP).
C
C    CL   IS ONE OF THE PARAMETERS THAT DEFINES F (SEE SUBROUTINE
C         FFUNP).
C
C    COEF  IS ONE OF THE PARAMETERS THAT DEFINES F (SEE SUBROUTINE
C         FFUNP).
C
C  ON OUTPUT:
C
C    RHO    IS THE HOMOTOPY.
C
C    DRHOX  CONTAINS THE PARTIAL DERIVATIVES OF RHO WITH RESPECT
C         TO X.
C
C    DRHOL  CONTAINS THE PARTIAL DERIVATIVES OF RHO WITH RESPECT
C         TO LAMBDA.
C
C  THE FOLLOWING ARE VARIABLES WHOSE WORKSPACE IS PASSED FROM HFUNP:
C    XDGM1
C    XDG
C    G
C    DG
C    PXDGM1
C    PXDG
C    F
C    DF
C    XX
C    TRM
C    DTRM
C    CLX
C    DXNP1
C    N
C    MMAXT
C    IDEG
C    NUMT
C    KDEG
C
C  OTHER VARIABLES:
C    ONEML
C
C  SUBROUTINES:  GFUNP, FFUNP.
C
C DECLARATION OF INPUT, WORKSPACE, AND OUTPUT:
      INTEGER N,MMAXT,IDEG,
     $ NUMT,KDEG
      DOUBLE PRECISION QDG,LAMBDA,X,
     $ PDG,CL,COEF,RHO,
     $ DRHOX,DRHOL,XDGM1,XDG,
     $ G,DG,PXDGM1,PXDG,
     $ F,DF,XX,TRM,
     $ DTRM,CLX,DXNP1
      DIMENSION IDEG(N),NUMT(N),KDEG(N,N+1,MMAXT)
      DIMENSION QDG(2,N),X(2,N),
     $  PDG(2,N),CL(2,N+1),COEF(N,MMAXT),RHO(2*N),
     $  DRHOX(2*N,2*N),DRHOL(2*N),XDGM1(2,N),XDG(2,N),
     $  G(2,N),DG(2,N),PXDGM1(2,N),PXDG(2,N),
     $  F(2,N), DF(2,N,N+1),XX(2,N,N+1,MMAXT),TRM(2,N,MMAXT),
     $  DTRM(2,N,N+1,MMAXT),CLX(2,N),DXNP1(2,N)
C
C DECLARATION OF VARIABLES:
      INTEGER J,J2,J2M1,K,K2,K2M1
      DOUBLE PRECISION ONEML
C
      CALL GFUNP(N,IDEG,PDG,QDG,X,XDGM1,XDG,PXDGM1,PXDG,G,DG)
      CALL FFUNP(N,NUMT,MMAXT,KDEG,COEF,CL,X,
     $  XX,TRM,DTRM,CLX,DXNP1,F,DF)
      ONEML=1.0 - LAMBDA
      DO 30 J=1,N
          J2=2*J
          J2M1=J2-1
          DO 20 K=1,N
              K2=2*K
              K2M1=K2-1
              DRHOX(J2M1,K2M1)= LAMBDA*DF(1,J,K)
              DRHOX(J2  ,K2  )= DRHOX(J2M1,K2M1)
              DRHOX(J2  ,K2M1)= LAMBDA*DF(2,J,K)
              DRHOX(J2M1,K2  )=-DRHOX(J2  ,K2M1)
  20      CONTINUE
          DRHOX(J2M1,J2M1)= DRHOX(J2M1,J2M1) + ONEML*DG(1,J)
          DRHOX(J2  ,J2  )= DRHOX(J2M1,J2M1)
          DRHOX(J2  ,J2M1)= DRHOX(J2  ,J2M1) + ONEML*DG(2,J)
          DRHOX(J2M1,J2  )=-DRHOX(J2  ,J2M1)
          DRHOL(J2M1)     =   F(1,J)      -        G(1,J)
          DRHOL(J2)       =   F(2,J)      -        G(2,J)
          RHO(J2M1)      = LAMBDA*F(1,J) + ONEML* G(1,J)
          RHO(J2  )      = LAMBDA*F(2,J) + ONEML* G(2,J)
  30  CONTINUE
      RETURN
      END