share/hompack/fortran/tangnf.f

   1       SUBROUTINE TANGNF(RHOLEN,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,
   2      $                  NFE,N,IFLAG,PAR,IPAR)
   3 C
   4 C THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP,
   5 C COMPUTES A QR DECOMPOSITION OF THAT MATRIX, AND THEN CALCULATES THE
   6 C (UNIT) TANGENT VECTOR AND THE NEWTON STEP.
   7 C
   8 C ON INPUT:
   9 C
  10 C RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT
  11 C    (A, LAMBDA, X) 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 (LAMBDA(S), X(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:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1)  ARE WORK ARRAYS
  22 C    USED FOR THE QR FACTORIZATION.
  23 C
  24 C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY
  25 C    FUNCTION EVALUATIONS.
  26 C
  27 C N = DIMENSION OF X.
  28 C
  29 C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
  30 C
  31 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  32 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  33 C    RHO, RHOJAC.
  34 C
  35 C ON OUTPUT:
  36 C
  37 C RHOLEN = ||RHO(A, LAMBDA(S), X(S)|| IF  RHOLEN < 0  ON INPUT.
  38 C    OTHERWISE  RHOLEN  IS UNCHANGED.
  39 C
  40 C Y, YPOLD, A, N  ARE UNCHANGED.
  41 C
  42 C YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF
  43 C    D(HOMOTOPY MAP)/DS = 0  AT  Y(S) = (LAMBDA(S), X(S)) .
  44 C
  45 C TZ = THE NEWTON STEP = -(PSEUDO INVERSE OF  (D RHO(A,Y(S))/D LAMBDA ,
  46 C    D RHO(A,Y(S))/DX)) * RHO(A,Y(S)) .
  47 C
  48 C NFE  HAS BEEN INCRMENTED BY 1.
  49 C
  50 C IFLAG  IS UNCHANGED, UNLESS THE QR FACTORIZATION DETECTS A RANK < N,
  51 C    IN WHICH CASE THE TANGENT AND NEWTON STEP VECTORS ARE NOT COMPUTED
  52 C    AND  TANGNF  RETURNS WITH  IFLAG = 4 .
  53 C
  54 C
  55 C CALLS  DDOT , DNRM2 , F (OR  RHO ), FJAC (OR  RHOJAC ).
  56 C
  57       DOUBLE PRECISION ALPHAK,BETA,DDOT,DNRM2,LAMBDA,QRKK,RHOLEN,
  58      $       SIGMA,SUM,YPNORM
  59       INTEGER I,IFLAG,J,JBAR,K,KP1,N,NFE,NP1,NP2
  60 C
  61 C *****  ARRAY DECLARATIONS.  *****
  62 C
  63       DOUBLE PRECISION Y(N+1),YP(N+1),YPOLD(N+1),A(N),PAR(1)
  64       INTEGER IPAR(1)
  65 C
  66 C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
  67       DOUBLE PRECISION QR(N,N+2),ALPHA(N),TZ(N+1)
  68       INTEGER PIVOT(N+1)
  69 C
  70 C *****  END OF DIMENSIONAL INFORMATION.  *****
  71 C
  72 C
  73       LAMBDA=Y(1)
  74       NP1=N+1
  75       NP2=N+2
  76       NFE=NFE+1
  77 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
  78 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
  79 C
  80 C COMPUTE THE JACOBIAN MATRIX, STORE IT AND HOMOTOPY MAP IN QR.
  81 C
  82       IF (IFLAG .EQ. -2) THEN
  83 C
  84 C  QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX ,
  85 C                                              RHO(A,LAMBDA,X) )  .
  86 C
  87         DO 30 K=1,NP1
  88           CALL RHOJAC(A,LAMBDA,Y(2),QR(1,K),K,PAR,IPAR)
  89 30      CONTINUE
  90         CALL RHO(A,LAMBDA,Y(2),QR(1,NP2),PAR,IPAR)
  91       ELSE
  92         CALL F(Y(2),TZ)
  93         IF (IFLAG .EQ. 0) THEN
  94 C
  95 C      QR = ( A - F(X), I - LAMBDA*DF(X) ,
  96 C                                 X - A + LAMBDA*(A - F(X)) )  .
  97 C
  98           DO 100 J=1,N
  99             SIGMA=A(J)
 100             BETA=SIGMA-TZ(J)
 101             QR(J,1)=BETA
 102 100       QR(J,NP2)=Y(J+1)-SIGMA+LAMBDA*BETA
 103           DO 120 K=1,N
 104             CALL FJAC(Y(2),TZ,K)
 105             KP1=K+1
 106             DO 110 J=1,N
 107 110         QR(J,KP1)=-LAMBDA*TZ(J)
 108 120       QR(K,KP1)=1.0+QR(K,KP1)
 109         ELSE
 110 C
 111 C   QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I ,
 112 C                                  X - A + LAMBDA*(F(X) - X + A) )  .
 113 C
 114 140       DO 150 J=1,N
 115             SIGMA=Y(J+1)-A(J)
 116             BETA=TZ(J)-SIGMA
 117             QR(J,1)=BETA
 118 150       QR(J,NP2)=SIGMA+LAMBDA*BETA
 119           DO 170 K=1,N
 120             CALL FJAC(Y(2),TZ,K)
 121             KP1=K+1
 122             DO 160 J=1,N
 123 160         QR(J,KP1)=LAMBDA*TZ(J)
 124 170       QR(K,KP1)=1.0-LAMBDA+QR(K,KP1)
 125         ENDIF
 126       ENDIF
 127 C
 128 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
 129 C COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED.
 130       IF (RHOLEN .LT. 0.0) RHOLEN=DNRM2(N,QR(1,NP2),1)
 131 C
 132 C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM.
 133 C
 134 C THE FOLLOWING CODE IS A MODIFICATION OF THE ALGOL PROCEDURE
 135 C DECOMPOSE  IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST
 136 C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS,
 137 C NUMER. MATH. 7 (1965) 269-276.
 138 C
 139       DO 220 J=1,NP1
 140         YP(J)=DDOT(N,QR(1,J),1,QR(1,J),1)
 141 220   PIVOT(J)=J
 142       DO 300 K=1,N
 143         SIGMA=YP(K)
 144         JBAR=K
 145         KP1=K+1
 146         DO 240 J=KP1,NP1
 147           IF (SIGMA .GE. YP(J)) GO TO 240
 148           SIGMA=YP(J)
 149           JBAR=J
 150 240     CONTINUE
 151         IF (JBAR .EQ. K) GO TO 260
 152         I=PIVOT(K)
 153         PIVOT(K)=PIVOT(JBAR)
 154         PIVOT(JBAR)=I
 155         YP(JBAR)=YP(K)
 156         YP(K)=SIGMA
 157         DO 250 I=1,N
 158           SIGMA=QR(I,K)
 159           QR(I,K)=QR(I,JBAR)
 160           QR(I,JBAR)=SIGMA
 161 250     CONTINUE
 162 C   END OF COLUMN INTERCHANGE.
 163 260     SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1)
 164         IF (SIGMA .EQ. 0.0) THEN
 165           IFLAG=4
 166           RETURN
 167         ENDIF
 168 270     IF (K .EQ. N) GO TO 300
 169         QRKK=QR(K,K)
 170         ALPHAK=-SQRT(SIGMA)
 171         IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK
 172         ALPHA(K)=ALPHAK
 173         BETA=1.0/(SIGMA-QRKK*ALPHAK)
 174         QR(K,K)=QRKK-ALPHAK
 175         DO 290 J=KP1,NP2
 176           SIGMA=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1)
 177           DO 280 I=K,N
 178             QR(I,J)=QR(I,J)-QR(I,K)*SIGMA
 179 280       CONTINUE
 180           IF (J .LT. NP2) YP(J)=YP(J)-QR(K,J)**2
 181 290     CONTINUE
 182 300   CONTINUE
 183       ALPHA(N)=QR(N,N)
 184 C
 185 C
 186 C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS.
 187       TZ(NP1)=1.0
 188       DO 340 I=N,1,-1
 189         SUM=0.0
 190         DO 330 J=I+1,NP1
 191 330     SUM=SUM+QR(I,J)*TZ(J)
 192 340   TZ(I)=-SUM/ALPHA(I)
 193       YPNORM=DNRM2(NP1,TZ,1)
 194       DO 360 K=1,NP1
 195 360   YP(PIVOT(K))=TZ(K)/YPNORM
 196       IF (DDOT(NP1,YP,1,YPOLD,1) .GE. 0.0) GO TO 380
 197       DO 370 I=1,NP1
 198 370   YP(I)=-YP(I)
 199 C YP  IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION.
 200 C
 201 C COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)).
 202 C V IS GIVEN BY  P - (P,Q)Q  , WHERE P IS ANY SOLUTION OF
 203 C [D RHO] V = -RHO  AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO].
 204 C
 205 380   DO 440 I=N,1,-1
 206         SUM=QR(I,NP1)+QR(I,NP2)
 207         DO 430 J=I+1,N
 208 430     SUM=SUM+QR(I,J)*ALPHA(J)
 209 440   ALPHA(I)=-SUM/ALPHA(I)
 210       DO 450 K=1,N
 211 450   TZ(PIVOT(K))=ALPHA(K)
 212       TZ(PIVOT(NP1))=1.0
 213 C TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q
 214 C IN THE KERNEL(THE TANGENT).
 215       SIGMA=DDOT(NP1,TZ,1,YP,1)
 216       DO 470 J=1,NP1
 217         TZ(J)=TZ(J)-SIGMA*YP(J)
 218 470   CONTINUE
 219 C TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (LAMBDA(S), X(S)).
 220       RETURN
 221       END