share/hompack/fortran/fodeds.f

   1       SUBROUTINE FODEDS(S,Y,YP,YPOLD,A,QR,LENQR,PIVOT,PP,WORK,
   2      $     NFE,N,IFLAG,PAR,IPAR)
   3 C
   4 C SUBROUTINE  FODEDS  IS USED BY SUBROUTINE  STEPDS  TO SPECIFY THE
   5 C ORDINARY DIFFERENTIAL EQUATION  DY/DS = G(S,Y) , WHOSE SOLUTION
   6 C IS THE ZERO CURVE OF THE HOMOTOPY MAP.  S = ARC LENGTH,
   7 C YP = DY/DS, AND  Y(S) = (X(S), LAMBDA(S)) .
   8 C
   9 C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
  10 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
  11 C    RHOA, RHOJS.
  12 C
  13       DOUBLE PRECISION DDOT,DNRM2,LAMBDA,S,TEMP,YPNORM
  14       INTEGER IFLAG,J,LENQR,N,NFE,NP1
  15 C
  16 C *****  ARRAY DECLARATIONS.  *****
  17 C
  18       DOUBLE PRECISION Y(N+1),YP(N+1),YPOLD(N+1),A(N),PAR(1)
  19       INTEGER IPAR(1)
  20 C
  21 C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL.
  22       DOUBLE PRECISION QR(LENQR),PP(N),WORK(6*(N+1)+LENQR)
  23       INTEGER PIVOT(N+2)
  24 C
  25 C *****  END OF DIMENSIONAL INFORMATION.  *****
  26 C
  27 C
  28       NP1=N+1
  29       NFE=NFE+1
  30 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS.
  31       LAMBDA=Y(NP1)
  32 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
  33 C
  34 C COMPUTE THE JACOBIAN MATRIX, STORE IT IN  [QR | -PP] .
  35 C
  36       IF (IFLAG .EQ. -2) THEN
  37 C
  38 C  [QR | -PP] = [ D RHO(A,X,LAMBDA)/DX | D RHO(A,X,LAMBDA)/D LAMBDA ]  .
  39 C
  40         CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR)
  41 C  PP = - (D RHO(A,X,LAMBDA)/D LAMBDA) .
  42 C
  43       ELSE
  44         CALL F(Y,PP)
  45         IF (IFLAG .EQ. 0) THEN
  46 C
  47 C      [QR | -PP] = [ I - LAMBDA*DF(X) | A - F(X) ]  .
  48 C
  49           CALL DAXPY(N,-1.0D0,A,1,PP,1)
  50           CALL FJACS(Y,QR,LENQR,PIVOT)
  51           CALL DSCAL(LENQR,-LAMBDA,QR,1)
  52           DO 120 J=1,N
  53             QR(PIVOT(J))=QR(PIVOT(J)) + 1.0
  54 120       CONTINUE
  55         ELSE
  56 C
  57 C   [QR | -PP] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] .
  58 C
  59           CALL DSCAL(N,-1.0D0,PP,1)
  60           CALL DAXPY(N,1.0D0,Y,1,PP,1)
  61           CALL DAXPY(N,-1.0D0,A,1,PP,1)
  62           CALL FJACS(Y,QR,LENQR,PIVOT)
  63           CALL DSCAL(LENQR,LAMBDA,QR,1)
  64           TEMP=1.0 - LAMBDA
  65           DO 170 J=1,N
  66             QR(PIVOT(J))=QR(PIVOT(J)) + TEMP
  67 170       CONTINUE
  68         ENDIF
  69       ENDIF
  70 C
  71 C   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
  72       CALL DCOPY(NP1,YPOLD,1,YP,1)
  73 C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, USING A
  74 C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
  75       CALL PCGDS(N,QR,LENQR,PIVOT,PP,YP,WORK,IFLAG)
  76       IF (IFLAG .GT. 0) RETURN
  77 C
  78 C NORMALIZE TANGENT VECTOR YP.
  79       YPNORM=DNRM2(NP1,YP,1)
  80       CALL DSCAL(NP1,1.0/YPNORM,YP,1)
  81 C
  82 C CHOOSE UNIT TANGENT VECTOR DIRECTION TO MAINTAIN CONTINUITY.
  83       IF (DDOT(NP1,YP,1,YPOLD,1) .LT. 0.0)
  84      $     CALL DSCAL(NP1,-1.0D0,YP,1)
  85 C
  86 C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN  YPOLD .
  87       CALL DCOPY(NP1,YP,1,YPOLD,1)
  88 C
  89       RETURN
  90       END