SUBROUTINE TANGQF(Y,YP,YPOLD,A,QT,R,W,S,T,N,IFLAG,NFE,PAR,IPAR)
C
C SUBROUTINE  TANGQF  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 IN ADDITION, THE MATRIX AUG IS UPDATED SO THAT THE LAST ROW IS
C YP  INSTEAD OF  YPOLD  ON RETURN.
C
C
C ON INPUT:
C
C Y(1:N+1) = CURRENT POINT (LAMBDA(S), X(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 W(1:N+1), S(1:N+1), T(1:N+1)  ARE WORK ARRAYS.
C
C N  IS THE DIMENSION OF X, WHERE  Y=(LAMBDA(S),X(S)).
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, RHOJAC.
C
C
C ON OUTPUT:
C
C Y, YPOLD, A, N  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) = (LAMBDA(S), X(S)).
C
C QT(1:N+1,1:N+1)  CONTAINS  Q TRANSPOSE  OF THE QR FACTORIZATION OF
C    THE JACOBIAN MATRIX OF RHO EVALUATED AT  Y  AUGMENTED BY
C    YP TRANSPOSE.
C
C R(1:(N+1)*(N+2)/2)  CONTAINS THE UPPER TRIANGLE (STORED BY ROWS)
C    OF THE  R  PART OF THE QR FACTORIZATION OF THE AUGMENTED JACOBIAN
C    MATRIX.
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), FJAC
C    (OR RHOJAC, IF IFLAG = -2), R1UPQF (WHICH IS AN ENTRY POINT OF
C    UPQRQF), QRFAQF, QRSLQF.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION DNRM2
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION LAMBDA, ONE, YPNRM
        INTEGER I, J, JP1, NP1
C
C     SCALAR ARGUMENTS
C
        INTEGER N, IFLAG, NFE
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION Y(N+1), YP(N+1), YPOLD(N+1), A(N),
     $    QT(N+1,N+1), R((N+1)*(N+2)/2), W(N+1), S(N+1), T(N+1),PAR(1)
        INTEGER IPAR(1)
C
C ***** END OF DECLARATIONS *****
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
        ONE = 1.0
        NFE = NFE + 1
        NP1 = N + 1
        LAMBDA = Y(1)
C
C ***** DEFINE THE AUGMENTED JACOBIAN MATRIX *****
C
C QT = AUG.
C
        IF (IFLAG .EQ. -2) THEN
C
C         CURVE TRACKING PROBLEM:
C         D(RHO) = (D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX).
C
          DO 10 J = 1,NP1
            CALL RHOJAC(A,LAMBDA,Y(2),QT(1,J),J,PAR,IPAR)
  10      CONTINUE
        ELSE IF (IFLAG .EQ. -1) THEN
C
C         ZERO FINDING PROBLEM:
C         D(RHO) = (F(X) - X + A, LAMBDA*DF(X) + (1-LAMBDA)*I)
C
          CALL F(Y(2),QT(1,1))
          DO 20 I=1,N
            QT(I,1) = A(I) - Y(I+1) + QT(I,1)
  20      CONTINUE
          DO 30 J= 1,N
            JP1 = J+1
            CALL FJAC(Y(2),QT(1,JP1),J)
            CALL DSCAL(N,LAMBDA,QT(1,JP1),1)
            QT(J,JP1) = 1.0 - LAMBDA + QT(J,JP1)
  30      CONTINUE
        ELSE
C
C         FIXED POINT PROBLEM:
C         D(RHO) = (A - F(X), I - LAMBDA*DF(X)).
C
          CALL F(Y(2),QT(1,1))
          CALL DSCAL(N,-ONE,QT(1,1),1)
          CALL DAXPY(N,ONE,A,1,QT(1,1),1)
          DO 50 J=1,N
            JP1 = J+1
            CALL FJAC(Y(2),QT(1,JP1),J)
            CALL DSCAL(N,-LAMBDA,QT(1,JP1),1)
            QT(J,JP1) = 1.0 + QT(J,JP1)
  50      CONTINUE
        END IF
C
C     DEFINE LAST ROW OF QT  = YPOLD.
C
        CALL DCOPY(NP1,YPOLD,1,QT(NP1,1),NP1)
C
C ***** END OF DEFINITION OF AUGMENTED JACOBIAN MATRIX *****
C
C                                          T
C ***** SOLVE SYSTEM  AUG*YPT = (0,...,0,1)  *****
C
C FACTOR MATRIX.
C
        CALL QRFAQF(QT,R,NP1,IFLAG)
C
C IF MATRIX IS SINGULAR, THEN QUIT.
C
C       IF (IFLAG .EQ. 4) RETURN
C
C ELSE SOLVE SYSTEM  R*YP = QT*(0,...,0,1)  FOR YP.
C
        DO 70 J=1,N
          YP(J) = 0.0
  70    CONTINUE
        YP(NP1) = 1.0
        CALL QRSLQF(QT,R,YP,W,NP1)
C
C     COMPUTE UNIT VECTOR.
C
        YPNRM = 1.0/DNRM2(NP1,YP,1)
        CALL DSCAL(NP1,YPNRM,YP,1)
C
C ***** SYSTEM SOLVED *****
C
C ***** UPDATE AUGMENTED SYSTEM SO THAT LAST ROW IS YP *****
C
C S=YP-YPOLD,  T = QT*E(NP1).
C
        CALL DCOPY(NP1,YP,1,S,1)
        CALL DAXPY(NP1,-ONE,YPOLD,1,S,1)
        CALL DCOPY(NP1,QT(1,NP1),1,T,1)
        CALL R1UPQF(NP1,S,T,QT,R,W)
C
        RETURN
C
C ***** END OF SUBROUTINE TANGQF *****
        END