SUBROUTINE FODE(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG, $ PAR,IPAR) C C SUBROUTINE FODE IS USED BY SUBROUTINE STEPS TO SPECIFY THE C ORDINARY DIFFERENTIAL EQUATION DY/DS = G(S,Y) , WHOSE SOLUTION C IS THE ZERO CURVE OF THE HOMOTOPY MAP. S = ARC LENGTH, C YP = DY/DS, AND Y(S) = (LAMBDA(S), X(S)) . 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 RHOA, RHOJAC. C DOUBLE PRECISION DDOT,DNRM2,S,SUM,YPNORM INTEGER I,IERR,IFLAG,IK,J,K,KP1,KPIV,LW,N,NFE,NP1 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+1),ALPHA(N),TZ(N+1) INTEGER PIVOT(N+1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C C NP1=N+1 NFE=NFE+1 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. C * * * * * * * * * * * * * * * * * C C COMPUTE THE JACOBIAN MATRIX, STORE IT 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 DO 30 K=1,NP1 CALL RHOJAC(A,Y(1),Y(2),QR(1,K),K,PAR,IPAR) 30 CONTINUE ELSE CALL F(Y(2),TZ) IF (IFLAG .EQ. 0) THEN C C QR = ( A - F(X), I - LAMBDA*DF(X) ) . C DO 100 J=1,N 100 QR(J,1)=A(J)-TZ(J) DO 120 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 110 J=1,N 110 QR(J,KP1)=-Y(1)*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 140 DO 150 J=1,N 150 QR(J,1)=TZ(J)-Y(J+1)+A(J) DO 170 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 160 J=1,N 160 QR(J,KP1)=Y(1)*TZ(J) 170 QR(K,KP1)=1.0-Y(1)+QR(K,KP1) ENDIF ENDIF C C * * * * * * * * * * * * * * * * * C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM. 210 CALL DCPOSE(N,N,QR,ALPHA,PIVOT,IERR,TZ,YP) IF (IERR .EQ. 0) GO TO 220 IFLAG=4 RETURN C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS. 220 TZ(NP1)=1.0 DO 240 LW=1,N I=NP1-LW IK=I+1 SUM=0.0 DO 230 J=IK,NP1 230 SUM=SUM+QR(I,J)*TZ(J) 240 TZ(I)=-SUM/ALPHA(I) YPNORM=DNRM2(NP1,TZ,1) DO 260 K=1,NP1 KPIV=PIVOT(K) 260 YP(KPIV)=TZ(K)/YPNORM IF (DDOT(NP1,YP,1,YPOLD,1) .GE. 0.0) GO TO 280 DO 270 I=1,NP1 270 YP(I)=-YP(I) C C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN YPOLD . 280 DO 290 I=1,NP1 290 YPOLD(I)=YP(I) RETURN END