SUBROUTINE TANGQS(Y,YP,YPOLD,A,QR,PIVOT,PP,RHOVEC,WORK,N,LENQR,
     $     IFLAG,NFE,PAR,IPAR)
C
C SUBROUTINE  TANGQS  COMPUTES THE UNIT TANGENT VECTOR  YP  TO THE
C ZERO CURVE OF THE HOMOTOPY MAP AT  Y  BY GENERATING THE AUGMENTED
C JACOBIAN MATRIX
C
C           --           --
C           |  D(RHO(Y))  |
C     AUG = |        T    |,   WHERE RHO IS THE HOMOTOPY MAP
C           |   YPOLD     |
C           --           --
C
C SOLVING THE SYSTEM
C                                T
C         AUG*YPT = (0,0,...,0,1)    FOR YPT.
C
C AND FINALLY COMPUTING  YP = YPT/||YPT||.
C
C
C ON INPUT:
C
C Y(1:N+1) = CURRENT POINT (X(S), LAMBDA(S)).
C
C YP(1:N+1)  IS UNDEFINED ON INPUT.
C
C YPOLD(1:N+1) = UNIT TANGENT VECTOR AT THE PREVIOUS POINT ON THE
C    ZERO CURVE OF THE HOMOTOPY MAP.
C
C A(1:N)  IS THE PARAMETER VECTOR IN THE HOMOTOPY MAP.
C
C QR(1:LENQR)  IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC
C    JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE
C    STORAGE FORMAT.  LENQR  AND  PIVOT  DESCRIBE THE DATA
C    STRUCTURE IN  QR.  (SEE SUBROUTINE  PCGQS  FOR A DESCRIPTION
C    OF THIS DATA STRUCTURE).
C
C PIVOT(1:N+2)  IS A WORK ARRAY WHOSE FIRST N+1 COMPONENTS CONTAINI
C    THE INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC
C    JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN  QR.
C
C PP(1:N)  IS A WORK ARRAY CONTAINING THE NEGATIVE OF THE LAST COLUMN
C    OF THE JACOBIAN MATRIX  -[D RHO/D LAMBDA].
C
C RHOVEC(1:N+1), IS A WORK ARRAY USED TO CALCULATE THE TANGENT
C    VECTOR.
C
C WORK(1:8*(N+1)+LENQR)  IS A WORK ARRAY USED BY THE CONJUGATE GRADIENT
C    ALGORITHM TO SOLVE LINEAR SYSTEMS.
C
C N  IS THE DIMENSION OF X, WHERE  Y=(X(S),LAMBDA(S)).
C
C LENQR  IS THE LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR.
C
C IFLAG  IS -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
C
C NFE  IS THE NUMBER OF JACOBIAN EVALUATIONS.
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    RHO, RHOJS.
C
C
C ON OUTPUT:
C
C Y, YPOLD, A, N, LENQR  ARE UNCHANGED.
C
C YP(1:N+1) CONTAINS THE NEW UNIT TANGENT VECTOR TO THE ZERO
C    CURVE OF THE HOMOTOPY MAP AT  Y(S) = (X(S),LAMBDA(S)).
C
C IFLAG  = -2, -1, OR 0, (UNCHANGED) ON A NORMAL RETURN.
C        = 4 IF THE AUGMENTED JACOBIAN MATRIX HAS RANK LESS THAN N+1.
C
C NFE  HAS BEEN INCREMENTED BY 1.
C
C
C CALLS  DCOPY, DNRM2, DSCAL, F (OR RHO IF IFLAG = -2), FJACS
C    (OR RHOJS, IF IFLAG = -2), PCGQS.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION DNRM2
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION LAMBDA, ONE, SIGMA, YPNRM
        INTEGER J, NP1, PCGWK, ZU
C
C     SCALAR ARGUMENTS
C
        INTEGER N, LENQR, IFLAG, NFE
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION Y(N+1), YP(N+1), YPOLD(N+1), A(N),
     $    QR(LENQR), PP(N), RHOVEC(N+1), WORK(8*(N+1)+LENQR),PAR(1)
        INTEGER PIVOT(N+2), IPAR(1)
C
C ***** END OF DECLARATIONS *****
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
        ONE = 1.0
        NFE = NFE + 1
        NP1 = N + 1
        LAMBDA = Y(NP1)
        PCGWK = 2*N+3
        ZU = 3*N+4
C
C ***** DEFINE THE AUGMENTED JACOBIAN MATRIX *****
C
C COMPUTE JACOBIAN MATRIX, STORE IT IN [QR|-PP].
C
        IF (IFLAG .EQ. -2) THEN
C
C         CURVE TRACKING PROBLEM.
C
          CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR)
        ELSE IF (IFLAG .EQ. -1) THEN
C
C         ZERO FINDING PROBLEM.
C
          CALL F(Y,PP)
          CALL DSCAL(N,-ONE,PP,1)
          CALL DAXPY(N,ONE,Y,1,PP,1)
          CALL DAXPY(N,-ONE,A,1,PP,1)
          CALL FJACS(Y,QR,LENQR,PIVOT)
          CALL DSCAL(LENQR,LAMBDA,QR,1)
          SIGMA = 1.0-LAMBDA
          DO 10 J=1,N
            QR(PIVOT(J))=QR(PIVOT(J))+SIGMA
  10      CONTINUE
        ELSE
C
C         FIXED POINT PROBLEM
C
          CALL F(Y,PP)
          CALL DAXPY(N,-ONE,A,1,PP,1)
          CALL FJACS(Y,QR,LENQR,PIVOT)
          CALL DSCAL(LENQR,-LAMBDA,QR,1)
          DO 20 J=1,N
            QR(PIVOT(J))=QR(PIVOT(J)) + 1.0
  20      CONTINUE
        ENDIF
C
C ***** END OF DEFINITION OF AUGMENTED JACOBIAN MATRIX *****
C
C                                          T
C ***** SOLVE SYSTEM  AUG*YPT = (0,...,0,1)  *****
C
C INITIALIZE STARTING POINT FOR THE CONJUGATE GRADIENT ALGORITHM
C TO BE THE SOLUTIONS FROM THE PREVIOUS CALL TO  TANGQS.
C
        CALL DCOPY(2*NP1,WORK,1,WORK(ZU),1)
C
C RHOVEC = -(0,...,0,1)**T
C
        DO 30 J=1,N
          RHOVEC(J)=0.0
30      CONTINUE
        RHOVEC(NP1) = -1.0
C
C SOLVE SYSTEM.
C
        CALL PCGQS(N,QR,LENQR,PIVOT,PP,YPOLD,RHOVEC,YP,WORK(PCGWK),
     $      IFLAG)
        IF (IFLAG .GT. 0) RETURN
C
C NORMALIZE THE TANGENT.
C
        YPNRM = 1.0/DNRM2(NP1,YP,1)
        CALL DSCAL(NP1,YPNRM,YP,1)
C
C SAVE SOLUTIONS FROM CONJUGATE GRADIENT ALGORITHM FOR NEXT CALL
C TO TANGQS.
C
        CALL DCOPY(2*NP1,WORK(ZU),1,WORK,1)
C
        RETURN
C
C ***** END OF SUBROUTINE TANGQS *****
        END