      SUBROUTINE TANGNF(RHOLEN,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,
     $                  NFE,N,IFLAG,PAR,IPAR)
C
C THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP,
C COMPUTES A QR DECOMPOSITION OF THAT MATRIX, AND THEN CALCULATES THE
C (UNIT) TANGENT VECTOR AND THE NEWTON STEP.
C
C ON INPUT:
C
C RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT
C    (A, LAMBDA, X) IS TO BE COMPUTED.  IF  RHOLEN >= 0  THE NORM IS NOT
C    COMPUTED AND  RHOLEN  IS NOT CHANGED.
C
C Y(1:N+1) = CURRENT POINT (LAMBDA(S), X(S)).
C
C YPOLD(1:N+1) = UNIT TANGENT VECTOR AT PREVIOUS POINT ON THE ZERO
C    CURVE OF THE HOMOTOPY MAP.
C
C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
C
C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1)  ARE WORK ARRAYS
C    USED FOR THE QR FACTORIZATION.
C
C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY
C    FUNCTION EVALUATIONS.
C
C N = DIMENSION OF X.
C
C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
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 ON OUTPUT:
C
C RHOLEN = ||RHO(A, LAMBDA(S), X(S)|| IF  RHOLEN < 0  ON INPUT.
C    OTHERWISE  RHOLEN  IS UNCHANGED.
C
C Y, YPOLD, A, N  ARE UNCHANGED.
C
C YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF
C    D(HOMOTOPY MAP)/DS = 0  AT  Y(S) = (LAMBDA(S), X(S)) .
C
C TZ = THE NEWTON STEP = -(PSEUDO INVERSE OF  (D RHO(A,Y(S))/D LAMBDA ,
C    D RHO(A,Y(S))/DX)) * RHO(A,Y(S)) .
C
C NFE  HAS BEEN INCRMENTED BY 1.
C
C IFLAG  IS UNCHANGED, UNLESS THE QR FACTORIZATION DETECTS A RANK < N,
C    IN WHICH CASE THE TANGENT AND NEWTON STEP VECTORS ARE NOT COMPUTED
C    AND  TANGNF  RETURNS WITH  IFLAG = 4 .
C
C
C CALLS  DDOT , DNRM2 , F (OR  RHO ), FJAC (OR  RHOJAC ).
C
      DOUBLE PRECISION ALPHAK,BETA,DDOT,DNRM2,LAMBDA,QRKK,RHOLEN,
     $       SIGMA,SUM,YPNORM
      INTEGER I,IFLAG,J,JBAR,K,KP1,N,NFE,NP1,NP2
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(N,N+2),ALPHA(N),TZ(N+1)
      INTEGER PIVOT(N+1)
C
C *****  END OF DIMENSIONAL INFORMATION.  *****
C
C
      LAMBDA=Y(1)
      NP1=N+1
      NP2=N+2
      NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C
C COMPUTE THE JACOBIAN MATRIX, STORE IT AND HOMOTOPY MAP IN QR.
C
      IF (IFLAG .EQ. -2) THEN
C
C  QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX ,
C                                              RHO(A,LAMBDA,X) )  .
C
        DO 30 K=1,NP1
          CALL RHOJAC(A,LAMBDA,Y(2),QR(1,K),K,PAR,IPAR)
30      CONTINUE
        CALL RHO(A,LAMBDA,Y(2),QR(1,NP2),PAR,IPAR)
      ELSE
        CALL F(Y(2),TZ)
        IF (IFLAG .EQ. 0) THEN
C
C      QR = ( A - F(X), I - LAMBDA*DF(X) ,
C                                 X - A + LAMBDA*(A - F(X)) )  .
C
          DO 100 J=1,N
            SIGMA=A(J)
            BETA=SIGMA-TZ(J)
            QR(J,1)=BETA
100       QR(J,NP2)=Y(J+1)-SIGMA+LAMBDA*BETA
          DO 120 K=1,N
            CALL FJAC(Y(2),TZ,K)
            KP1=K+1
            DO 110 J=1,N
110         QR(J,KP1)=-LAMBDA*TZ(J)
120       QR(K,KP1)=1.0+QR(K,KP1)
        ELSE
C
C   QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I ,
C                                  X - A + LAMBDA*(F(X) - X + A) )  .
C
140       DO 150 J=1,N
            SIGMA=Y(J+1)-A(J)
            BETA=TZ(J)-SIGMA
            QR(J,1)=BETA
150       QR(J,NP2)=SIGMA+LAMBDA*BETA
          DO 170 K=1,N
            CALL FJAC(Y(2),TZ,K)
            KP1=K+1
            DO 160 J=1,N
160         QR(J,KP1)=LAMBDA*TZ(J)
170       QR(K,KP1)=1.0-LAMBDA+QR(K,KP1)
        ENDIF
      ENDIF
C
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED.
      IF (RHOLEN .LT. 0.0) RHOLEN=DNRM2(N,QR(1,NP2),1)
C
C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM.
C
C THE FOLLOWING CODE IS A MODIFICATION OF THE ALGOL PROCEDURE
C DECOMPOSE  IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST
C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS,
C NUMER. MATH. 7 (1965) 269-276.
C
      DO 220 J=1,NP1
        YP(J)=DDOT(N,QR(1,J),1,QR(1,J),1)
220   PIVOT(J)=J
      DO 300 K=1,N
        SIGMA=YP(K)
        JBAR=K
        KP1=K+1
        DO 240 J=KP1,NP1
          IF (SIGMA .GE. YP(J)) GO TO 240
          SIGMA=YP(J)
          JBAR=J
240     CONTINUE
        IF (JBAR .EQ. K) GO TO 260
        I=PIVOT(K)
        PIVOT(K)=PIVOT(JBAR)
        PIVOT(JBAR)=I
        YP(JBAR)=YP(K)
        YP(K)=SIGMA
        DO 250 I=1,N
          SIGMA=QR(I,K)
          QR(I,K)=QR(I,JBAR)
          QR(I,JBAR)=SIGMA
250     CONTINUE
C   END OF COLUMN INTERCHANGE.
260     SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1)
        IF (SIGMA .EQ. 0.0) THEN
          IFLAG=4
          RETURN
        ENDIF
270     IF (K .EQ. N) GO TO 300
        QRKK=QR(K,K)
        ALPHAK=-SQRT(SIGMA)
        IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK
        ALPHA(K)=ALPHAK
        BETA=1.0/(SIGMA-QRKK*ALPHAK)
        QR(K,K)=QRKK-ALPHAK
        DO 290 J=KP1,NP2
          SIGMA=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1)
          DO 280 I=K,N
            QR(I,J)=QR(I,J)-QR(I,K)*SIGMA
280       CONTINUE
          IF (J .LT. NP2) YP(J)=YP(J)-QR(K,J)**2
290     CONTINUE
300   CONTINUE
      ALPHA(N)=QR(N,N)
C
C
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS.
      TZ(NP1)=1.0
      DO 340 I=N,1,-1
        SUM=0.0
        DO 330 J=I+1,NP1
330     SUM=SUM+QR(I,J)*TZ(J)
340   TZ(I)=-SUM/ALPHA(I)
      YPNORM=DNRM2(NP1,TZ,1)
      DO 360 K=1,NP1
360   YP(PIVOT(K))=TZ(K)/YPNORM
      IF (DDOT(NP1,YP,1,YPOLD,1) .GE. 0.0) GO TO 380
      DO 370 I=1,NP1
370   YP(I)=-YP(I)
C YP  IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION.
C
C COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)).
C V IS GIVEN BY  P - (P,Q)Q  , WHERE P IS ANY SOLUTION OF
C [D RHO] V = -RHO  AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO].
C
380   DO 440 I=N,1,-1
        SUM=QR(I,NP1)+QR(I,NP2)
        DO 430 J=I+1,N
430     SUM=SUM+QR(I,J)*ALPHA(J)
440   ALPHA(I)=-SUM/ALPHA(I)
      DO 450 K=1,N
450   TZ(PIVOT(K))=ALPHA(K)
      TZ(PIVOT(NP1))=1.0
C TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q
C IN THE KERNEL(THE TANGENT).
      SIGMA=DDOT(NP1,TZ,1,YP,1)
      DO 470 J=1,NP1
        TZ(J)=TZ(J)-SIGMA*YP(J)
470   CONTINUE
C TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (LAMBDA(S), X(S)).
      RETURN
      END
