SUBROUTINE TANGNS(RHOLEN,Y,YP,TZ,YPOLD,A,QR,LENQR,PIVOT,
     &     PP,RHOVEC,WORK,NFE,N,IFLAG,PAR,IPAR)
C
C THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP,
C AND THEN CALCULATES THE (UNIT) TANGENT VECTOR AND THE NEWTON STEP
C USING A PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
C
C ON INPUT:
C
C RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT
C    (A, X, LAMBDA) 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 (X(S), LAMBDA(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:LENQR), PIVOT(1:N+2), PP(1:N), RHOVEC(1:N), WORK(1:8*(N+1)+LENQR)
C    ARE WORK ARRAYS USED FOR THE JACOBIAN MATRIX AND CONJUGATE GRADIENT
C    ITERATION.
C
C LENQR = LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR  USED TO CONTAIN THE
C    N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X IN PACKED SKYLINE
C    STORAGE FORMAT.
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, RHOJS.
C
C ON OUTPUT:
C
C RHOLEN = ||RHO(A, X(S), LAMBDA(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) = (X(S), LAMBDA(S)) .
C
C TZ(1:N+1) = THE NEWTON STEP = -(PSEUDO INVERSE OF  (D RHO(A,Y(S))/DX ,
C    D RHO(A,Y(S))/D LAMBDA)) * RHO(A,Y(S)) .
C
C NFE  HAS BEEN INCRMENTED BY 1.
C
C IFLAG  IS UNCHANGED, UNLESS THE PRECONDITIONED CONJUGATE GRADIENT
C    ITERATION FAILS TO CONVERGE, IN WHICH CASE THE TANGENT AND NEWTON 
C    STEP VECTORS ARE NOT COMPUTED AND  TANGNS  RETURNS WITH  IFLAG = 4 .
C
C
C CALLS  F (OR  RHO ), FJACS (OR  RHOJS ), PCGDS , PCGNS , AND THE BLAS
C    ROUTINES  DAXPY , DCOPY , DDOT , DNRM2 , DSCAL .
C
      DOUBLE PRECISION DDOT,DNRM2,LAMBDA,RHOLEN,SIGMA,YPNORM
      INTEGER IFLAG,J,LENQR,N,NFE,NP1,NP2,N2P3,N3P4,N4P5
C
C *****  ARRAY DECLARATIONS.  *****
C
      DOUBLE PRECISION Y(N+1),YP(N+1),TZ(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),RHOVEC(N),
     &     WORK(8*(N+1)+LENQR)
      INTEGER PIVOT(N+2)
C
C *****  END OF DIMENSIONAL INFORMATION.  *****
C
C
      NP1=N+1
      NP2=N+2
      N2P3=2*N+3
      N3P4=3*N+4
      N4P5=4*N+5
      NFE=NFE+1
C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C
C COMPUTE THE JACOBIAN MATRIX, STORE IT IN  [QR | -PP] .
C
      LAMBDA=Y(NP1)
      IF (IFLAG .EQ. -2) THEN
C
C  [QR | -PP] = ( D RHO(A,X,LAMBDA)/DX , D RHO(A,X,LAMBDA)/D LAMBDA ) ,
C  RHOVEC = RHO(A,X,LAMBDA) .
C
        CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR)
        CALL RHO(A,LAMBDA,Y,RHOVEC,PAR,IPAR)
      ELSE
        CALL F(Y,PP)
        CALL DCOPY(N,Y,1,RHOVEC,1)
        CALL DAXPY(N,-1.0D0,A,1,RHOVEC,1)
        IF (IFLAG .EQ. 0) THEN
C
C      [QR | -PP] = ( I - LAMBDA*DF(X) , A - F(X) ) ,
C      RHOVEC = X - A + LAMBDA*(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
          CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
        ELSE
C
C   [QR | -PP] = ( LAMBDA*DF(X) + (1 - LAMBDA)*I , F(X) - X + A ) ,
C   RHOVEC = X - A + LAMBDA*(F(X) - X + A) .
C
          CALL DSCAL(N,-1.0D0,PP,1)
          CALL DAXPY(N,1.0D0,RHOVEC,1,PP,1)
          CALL FJACS(Y,QR,LENQR,PIVOT)
          CALL DSCAL(LENQR,LAMBDA,QR,1)
          SIGMA=1.0 - LAMBDA
          DO 170 J=1,N
            QR(PIVOT(J))=QR(PIVOT(J)) + SIGMA
170       CONTINUE
          CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
        ENDIF
      ENDIF
C
C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
C COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED.
      IF (RHOLEN .LT. 0.0) RHOLEN=DNRM2(N,RHOVEC,1)
C
C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, BY A
C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.  THIS IS DONE BY SOLVING
C SEVERAL AUXILLARY SYSTEMS, WHOSE PREVIOUS SOLUTIONS HAVE BEEN LEFT IN
C WORK(1:N+1) AND WORK(N+2:2*N+2).
      CALL DCOPY(2*NP1,WORK,1,WORK(N3P4),1)
      CALL DCOPY(NP1,YPOLD,1,YP,1)
      CALL PCGDS(N,QR,LENQR,PIVOT,PP,YP,WORK(N2P3),IFLAG)
      IF (IFLAG .GT. 0) RETURN
      CALL DCOPY(2*NP1,WORK(N3P4),1,WORK,1)
      YPNORM=DNRM2(NP1,YP,1)
      CALL DSCAL(NP1,1.0/YPNORM,YP,1)
      IF (DDOT(NP1,YP,1,YPOLD,1) .LT. 0.0) 
     &     CALL DSCAL(NP1,-1.0D0,YP,1)
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
      CALL DSCAL(2*NP1,0.0D0,WORK(N3P4),1)
      CALL DCOPY(NP1,YPOLD,1,TZ,1)
      CALL PCGNS(N,QR,LENQR,PIVOT,PP,RHOVEC,TZ,WORK(N2P3),IFLAG)
      IF (IFLAG .GT. 0) RETURN
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)
      CALL DAXPY(NP1,-SIGMA,YP,1,TZ,1)
C
C TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (X(S), LAMBDA(S)).
C
      RETURN
      END