SUBROUTINE FODEDS(S,Y,YP,YPOLD,A,QR,LENQR,PIVOT,PP,WORK,
     $     NFE,N,IFLAG,PAR,IPAR)
C
C SUBROUTINE  FODEDS  IS USED BY SUBROUTINE  STEPDS  TO SPECIFY THE
C ORDINARY DIFFERENTIAL EQUATION  DY/DS = G(S,Y) , WHOSE SOLUTION
C IS THE ZERO CURVE OF THE HOMOTOPY MAP.  S = ARC LENGTH,
C YP = DY/DS, AND  Y(S) = (X(S), LAMBDA(S)) .
C
C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHOA, RHOJS.
C
      DOUBLE PRECISION DDOT,DNRM2,LAMBDA,S,TEMP,YPNORM
      INTEGER IFLAG,J,LENQR,N,NFE,NP1
C
C *****  ARRAY DECLARATIONS.  *****
C
      DOUBLE PRECISION Y(N+1),YP(N+1),YPOLD(N+1),A(N),PAR(1)
      INTEGER IPAR(1)
C
C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
      DOUBLE PRECISION QR(LENQR),PP(N),WORK(6*(N+1)+LENQR)
      INTEGER PIVOT(N+2)
C
C *****  END OF DIMENSIONAL INFORMATION.  *****
C
C
      NP1=N+1
      NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
      LAMBDA=Y(NP1)
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C
C COMPUTE THE JACOBIAN MATRIX, STORE IT IN  [QR | -PP] .
C
      IF (IFLAG .EQ. -2) THEN
C
C  [QR | -PP] = [ D RHO(A,X,LAMBDA)/DX | D RHO(A,X,LAMBDA)/D LAMBDA ]  .
C
        CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR)
C  PP = - (D RHO(A,X,LAMBDA)/D LAMBDA) .
C
      ELSE
        CALL F(Y,PP)
        IF (IFLAG .EQ. 0) THEN
C
C      [QR | -PP] = [ I - LAMBDA*DF(X) | A - F(X) ]  .
C
          CALL DAXPY(N,-1.0D0,A,1,PP,1)
          CALL FJACS(Y,QR,LENQR,PIVOT)
          CALL DSCAL(LENQR,-LAMBDA,QR,1)
          DO 120 J=1,N
            QR(PIVOT(J))=QR(PIVOT(J)) + 1.0
120       CONTINUE
        ELSE
C
C   [QR | -PP] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] .
C
          CALL DSCAL(N,-1.0D0,PP,1)
          CALL DAXPY(N,1.0D0,Y,1,PP,1)
          CALL DAXPY(N,-1.0D0,A,1,PP,1)
          CALL FJACS(Y,QR,LENQR,PIVOT)
          CALL DSCAL(LENQR,LAMBDA,QR,1)
          TEMP=1.0 - LAMBDA
          DO 170 J=1,N
            QR(PIVOT(J))=QR(PIVOT(J)) + TEMP
170       CONTINUE
        ENDIF
      ENDIF
C
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
      CALL DCOPY(NP1,YPOLD,1,YP,1)
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, USING A
C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
      CALL PCGDS(N,QR,LENQR,PIVOT,PP,YP,WORK,IFLAG)
      IF (IFLAG .GT. 0) RETURN
C
C NORMALIZE TANGENT VECTOR YP.
      YPNORM=DNRM2(NP1,YP,1)
      CALL DSCAL(NP1,1.0/YPNORM,YP,1)
C
C CHOOSE UNIT TANGENT VECTOR DIRECTION TO MAINTAIN CONTINUITY.
      IF (DDOT(NP1,YP,1,YPOLD,1) .LT. 0.0)
     $     CALL DSCAL(NP1,-1.0D0,YP,1)
C
C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN  YPOLD .
      CALL DCOPY(NP1,YP,1,YPOLD,1)
C
      RETURN
      END