share/hompack/fortran/tangns.f

   1       SUBROUTINE TANGNS(RHOLEN,Y,YP,TZ,YPOLD,A,QR,LENQR,PIVOT,
   2      &     PP,RHOVEC,WORK,NFE,N,IFLAG,PAR,IPAR)
   3 C
   4 C THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP,
   5 C AND THEN CALCULATES THE (UNIT) TANGENT VECTOR AND THE NEWTON STEP
   6 C USING A PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
   7 C
   8 C ON INPUT:
   9 C
  10 C RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT
  11 C    (A, X, LAMBDA) IS TO BE COMPUTED.  IF  RHOLEN >= 0  THE NORM IS NOT
  12 C    COMPUTED AND  RHOLEN  IS NOT CHANGED.
  13 C
  14 C Y(1:N+1) = CURRENT POINT (X(S), LAMBDA(S)).
  15 C
  16 C YPOLD(1:N+1) = UNIT TANGENT VECTOR AT PREVIOUS POINT ON THE ZERO
  17 C    CURVE OF THE HOMOTOPY MAP.
  18 C
  19 C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
  20 C
  21 C QR(1:LENQR), PIVOT(1:N+2), PP(1:N), RHOVEC(1:N), WORK(1:8*(N+1)+LENQR)
  22 C    ARE WORK ARRAYS USED FOR THE JACOBIAN MATRIX AND CONJUGATE GRADIENT
  23 C    ITERATION.
  24 C
  25 C LENQR = LENGTH OF THE ONE-DIMENSIONAL ARRAY  QR  USED TO CONTAIN THE
  26 C    N X N SYMMETRIC JACOBIAN MATRIX WITH RESPECT TO X IN PACKED SKYLINE
  27 C    STORAGE FORMAT.
  28 C
  29 C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY
  30 C    FUNCTION EVALUATIONS.
  31 C
  32 C N = DIMENSION OF X.
  33 C
  34 C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
  35 C
  36 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  37 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  38 C    RHO, RHOJS.
  39 C
  40 C ON OUTPUT:
  41 C
  42 C RHOLEN = ||RHO(A, X(S), LAMBDA(S)|| IF  RHOLEN < 0  ON INPUT.
  43 C    OTHERWISE  RHOLEN  IS UNCHANGED.
  44 C
  45 C Y, YPOLD, A, N  ARE UNCHANGED.
  46 C
  47 C YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF
  48 C    D(HOMOTOPY MAP)/DS = 0  AT  Y(S) = (X(S), LAMBDA(S)) .
  49 C
  50 C TZ(1:N+1) = THE NEWTON STEP = -(PSEUDO INVERSE OF  (D RHO(A,Y(S))/DX ,
  51 C    D RHO(A,Y(S))/D LAMBDA)) * RHO(A,Y(S)) .
  52 C
  53 C NFE  HAS BEEN INCRMENTED BY 1.
  54 C
  55 C IFLAG  IS UNCHANGED, UNLESS THE PRECONDITIONED CONJUGATE GRADIENT
  56 C    ITERATION FAILS TO CONVERGE, IN WHICH CASE THE TANGENT AND NEWTON
  57 C    STEP VECTORS ARE NOT COMPUTED AND  TANGNS  RETURNS WITH  IFLAG = 4 .
  58 C
  59 C
  60 C CALLS  F (OR  RHO ), FJACS (OR  RHOJS ), PCGDS , PCGNS , AND THE BLAS
  61 C    ROUTINES  DAXPY , DCOPY , DDOT , DNRM2 , DSCAL .
  62 C
  63       DOUBLE PRECISION DDOT,DNRM2,LAMBDA,RHOLEN,SIGMA,YPNORM
  64       INTEGER IFLAG,J,LENQR,N,NFE,NP1,NP2,N2P3,N3P4,N4P5
  65 C
  66 C *****  ARRAY DECLARATIONS.  *****
  67 C
  68       DOUBLE PRECISION Y(N+1),YP(N+1),TZ(N+1),YPOLD(N+1),A(N),PAR(1)
  69       INTEGER IPAR(1)
  70 C
  71 C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
  72       DOUBLE PRECISION QR(LENQR),PP(N),RHOVEC(N),
  73      &     WORK(8*(N+1)+LENQR)
  74       INTEGER PIVOT(N+2)
  75 C
  76 C *****  END OF DIMENSIONAL INFORMATION.  *****
  77 C
  78 C
  79       NP1=N+1
  80       NP2=N+2
  81       N2P3=2*N+3
  82       N3P4=3*N+4
  83       N4P5=4*N+5
  84       NFE=NFE+1
  85 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
  86 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
  87 C
  88 C COMPUTE THE JACOBIAN MATRIX, STORE IT IN  [QR | -PP] .
  89 C
  90       LAMBDA=Y(NP1)
  91       IF (IFLAG .EQ. -2) THEN
  92 C
  93 C  [QR | -PP] = ( D RHO(A,X,LAMBDA)/DX , D RHO(A,X,LAMBDA)/D LAMBDA ) ,
  94 C  RHOVEC = RHO(A,X,LAMBDA) .
  95 C
  96         CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR)
  97         CALL RHO(A,LAMBDA,Y,RHOVEC,PAR,IPAR)
  98       ELSE
  99         CALL F(Y,PP)
 100         CALL DCOPY(N,Y,1,RHOVEC,1)
 101         CALL DAXPY(N,-1.0D0,A,1,RHOVEC,1)
 102         IF (IFLAG .EQ. 0) THEN
 103 C
 104 C      [QR | -PP] = ( I - LAMBDA*DF(X) , A - F(X) ) ,
 105 C      RHOVEC = X - A + LAMBDA*(A - F(X)) .
 106 C
 107           CALL DAXPY(N,-1.0D0,A,1,PP,1)
 108           CALL FJACS(Y,QR,LENQR,PIVOT)
 109           CALL DSCAL(LENQR,-LAMBDA,QR,1)
 110           DO 120 J=1,N
 111             QR(PIVOT(J))=QR(PIVOT(J)) + 1.0
 112 120       CONTINUE
 113           CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
 114         ELSE
 115 C
 116 C   [QR | -PP] = ( LAMBDA*DF(X) + (1 - LAMBDA)*I , F(X) - X + A ) ,
 117 C   RHOVEC = X - A + LAMBDA*(F(X) - X + A) .
 118 C
 119           CALL DSCAL(N,-1.0D0,PP,1)
 120           CALL DAXPY(N,1.0D0,RHOVEC,1,PP,1)
 121           CALL FJACS(Y,QR,LENQR,PIVOT)
 122           CALL DSCAL(LENQR,LAMBDA,QR,1)
 123           SIGMA=1.0 - LAMBDA
 124           DO 170 J=1,N
 125             QR(PIVOT(J))=QR(PIVOT(J)) + SIGMA
 126 170       CONTINUE
 127           CALL DAXPY(N,-LAMBDA,PP,1,RHOVEC,1)
 128         ENDIF
 129       ENDIF
 130 C
 131 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
 132 C COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED.
 133       IF (RHOLEN .LT. 0.0) RHOLEN=DNRM2(N,RHOVEC,1)
 134 C
 135 C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, BY A
 136 C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.  THIS IS DONE BY SOLVING
 137 C SEVERAL AUXILLARY SYSTEMS, WHOSE PREVIOUS SOLUTIONS HAVE BEEN LEFT IN
 138 C WORK(1:N+1) AND WORK(N+2:2*N+2).
 139       CALL DCOPY(2*NP1,WORK,1,WORK(N3P4),1)
 140       CALL DCOPY(NP1,YPOLD,1,YP,1)
 141       CALL PCGDS(N,QR,LENQR,PIVOT,PP,YP,WORK(N2P3),IFLAG)
 142       IF (IFLAG .GT. 0) RETURN
 143       CALL DCOPY(2*NP1,WORK(N3P4),1,WORK,1)
 144       YPNORM=DNRM2(NP1,YP,1)
 145       CALL DSCAL(NP1,1.0/YPNORM,YP,1)
 146       IF (DDOT(NP1,YP,1,YPOLD,1) .LT. 0.0)
 147      &     CALL DSCAL(NP1,-1.0D0,YP,1)
 148 C YP  IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION.
 149 C
 150 C COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)).
 151 C V IS GIVEN BY  P - (P,Q)Q  , WHERE P IS ANY SOLUTION OF
 152 C [D RHO] V = -RHO  AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO].
 153 C
 154       CALL DSCAL(2*NP1,0.0D0,WORK(N3P4),1)
 155       CALL DCOPY(NP1,YPOLD,1,TZ,1)
 156       CALL PCGNS(N,QR,LENQR,PIVOT,PP,RHOVEC,TZ,WORK(N2P3),IFLAG)
 157       IF (IFLAG .GT. 0) RETURN
 158 C TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q
 159 C IN THE KERNEL(THE TANGENT).
 160       SIGMA=DDOT(NP1,TZ,1,YP,1)
 161       CALL DAXPY(NP1,-SIGMA,YP,1,TZ,1)
 162 C
 163 C TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (X(S), LAMBDA(S)).
 164 C
 165       RETURN
 166       END