SUBROUTINE FODEDS(S,Y,YP,YPOLD,A,QR,LENQR,PIVOT,PP,WORK, $ NFE,N,IFLAG,PAR,IPAR) C C SUBROUTINE FODEDS IS USED BY SUBROUTINE STEPDS 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) = (X(S), LAMBDA(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, RHOJS. C DOUBLE PRECISION DDOT,DNRM2,LAMBDA,S,TEMP,YPNORM INTEGER IFLAG,J,LENQR,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(LENQR),PP(N),WORK(6*(N+1)+LENQR) INTEGER PIVOT(N+2) C C ***** END OF DIMENSIONAL INFORMATION. ***** C C NP1=N+1 NFE=NFE+1 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. LAMBDA=Y(NP1) C * * * * * * * * * * * * * * * * * C C COMPUTE THE JACOBIAN MATRIX, STORE IT IN [QR | -PP] . C IF (IFLAG .EQ. -2) THEN C C [QR | -PP] = [ D RHO(A,X,LAMBDA)/DX | D RHO(A,X,LAMBDA)/D LAMBDA ] . C CALL RHOJS(A,LAMBDA,Y,QR,LENQR,PIVOT,PP,PAR,IPAR) C PP = - (D RHO(A,X,LAMBDA)/D LAMBDA) . C ELSE CALL F(Y,PP) IF (IFLAG .EQ. 0) THEN C C [QR | -PP] = [ I - LAMBDA*DF(X) | 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 ELSE C C [QR | -PP] = [ LAMBDA*DF(X) + (1 - LAMBDA)*I | F(X) - X + A ] . C CALL DSCAL(N,-1.0D0,PP,1) CALL DAXPY(N,1.0D0,Y,1,PP,1) CALL DAXPY(N,-1.0D0,A,1,PP,1) CALL FJACS(Y,QR,LENQR,PIVOT) CALL DSCAL(LENQR,LAMBDA,QR,1) TEMP=1.0 - LAMBDA DO 170 J=1,N QR(PIVOT(J))=QR(PIVOT(J)) + TEMP 170 CONTINUE ENDIF ENDIF C C * * * * * * * * * * * * * * * * * CALL DCOPY(NP1,YPOLD,1,YP,1) C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS, USING A C PRECONDITIONED CONJUGATE GRADIENT ALGORITHM. CALL PCGDS(N,QR,LENQR,PIVOT,PP,YP,WORK,IFLAG) IF (IFLAG .GT. 0) RETURN C C NORMALIZE TANGENT VECTOR YP. YPNORM=DNRM2(NP1,YP,1) CALL DSCAL(NP1,1.0/YPNORM,YP,1) C C CHOOSE UNIT TANGENT VECTOR DIRECTION TO MAINTAIN CONTINUITY. IF (DDOT(NP1,YP,1,YPOLD,1) .LT. 0.0) $ CALL DSCAL(NP1,-1.0D0,YP,1) C C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN YPOLD . CALL DCOPY(NP1,YP,1,YPOLD,1) C RETURN END