SUBROUTINE TANGNF(RHOLEN,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT, $ NFE,N,IFLAG,PAR,IPAR) C C THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP, C COMPUTES A QR DECOMPOSITION OF THAT MATRIX, AND THEN CALCULATES THE C (UNIT) TANGENT VECTOR AND THE NEWTON STEP. C C ON INPUT: C C RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT C (A, LAMBDA, X) IS TO BE COMPUTED. IF RHOLEN >= 0 THE NORM IS NOT C COMPUTED AND RHOLEN IS NOT CHANGED. C C Y(1:N+1) = CURRENT POINT (LAMBDA(S), X(S)). C C YPOLD(1:N+1) = UNIT TANGENT VECTOR AT PREVIOUS POINT ON THE ZERO C CURVE OF THE HOMOTOPY MAP. C C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. C C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) ARE WORK ARRAYS C USED FOR THE QR FACTORIZATION. C C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY C FUNCTION EVALUATIONS. C C N = DIMENSION OF X. C C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. 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 RHO, RHOJAC. C C ON OUTPUT: C C RHOLEN = ||RHO(A, LAMBDA(S), X(S)|| IF RHOLEN < 0 ON INPUT. C OTHERWISE RHOLEN IS UNCHANGED. C C Y, YPOLD, A, N ARE UNCHANGED. C C YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF C D(HOMOTOPY MAP)/DS = 0 AT Y(S) = (LAMBDA(S), X(S)) . C C TZ = THE NEWTON STEP = -(PSEUDO INVERSE OF (D RHO(A,Y(S))/D LAMBDA , C D RHO(A,Y(S))/DX)) * RHO(A,Y(S)) . C C NFE HAS BEEN INCRMENTED BY 1. C C IFLAG IS UNCHANGED, UNLESS THE QR FACTORIZATION DETECTS A RANK < N, C IN WHICH CASE THE TANGENT AND NEWTON STEP VECTORS ARE NOT COMPUTED C AND TANGNF RETURNS WITH IFLAG = 4 . C C C CALLS DDOT , DNRM2 , F (OR RHO ), FJAC (OR RHOJAC ). C DOUBLE PRECISION ALPHAK,BETA,DDOT,DNRM2,LAMBDA,QRKK,RHOLEN, $ SIGMA,SUM,YPNORM INTEGER I,IFLAG,J,JBAR,K,KP1,N,NFE,NP1,NP2 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+2),ALPHA(N),TZ(N+1) INTEGER PIVOT(N+1) C C ***** END OF DIMENSIONAL INFORMATION. ***** C C LAMBDA=Y(1) NP1=N+1 NP2=N+2 NFE=NFE+1 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. C * * * * * * * * * * * * * * * * * C C COMPUTE THE JACOBIAN MATRIX, STORE IT AND HOMOTOPY MAP 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 RHO(A,LAMBDA,X) ) . C DO 30 K=1,NP1 CALL RHOJAC(A,LAMBDA,Y(2),QR(1,K),K,PAR,IPAR) 30 CONTINUE CALL RHO(A,LAMBDA,Y(2),QR(1,NP2),PAR,IPAR) ELSE CALL F(Y(2),TZ) IF (IFLAG .EQ. 0) THEN C C QR = ( A - F(X), I - LAMBDA*DF(X) , C X - A + LAMBDA*(A - F(X)) ) . C DO 100 J=1,N SIGMA=A(J) BETA=SIGMA-TZ(J) QR(J,1)=BETA 100 QR(J,NP2)=Y(J+1)-SIGMA+LAMBDA*BETA DO 120 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 110 J=1,N 110 QR(J,KP1)=-LAMBDA*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 X - A + LAMBDA*(F(X) - X + A) ) . C 140 DO 150 J=1,N SIGMA=Y(J+1)-A(J) BETA=TZ(J)-SIGMA QR(J,1)=BETA 150 QR(J,NP2)=SIGMA+LAMBDA*BETA DO 170 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 160 J=1,N 160 QR(J,KP1)=LAMBDA*TZ(J) 170 QR(K,KP1)=1.0-LAMBDA+QR(K,KP1) ENDIF ENDIF C C * * * * * * * * * * * * * * * * * C COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED. IF (RHOLEN .LT. 0.0) RHOLEN=DNRM2(N,QR(1,NP2),1) C C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM. C C THE FOLLOWING CODE IS A MODIFICATION OF THE ALGOL PROCEDURE C DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, C NUMER. MATH. 7 (1965) 269-276. C DO 220 J=1,NP1 YP(J)=DDOT(N,QR(1,J),1,QR(1,J),1) 220 PIVOT(J)=J DO 300 K=1,N SIGMA=YP(K) JBAR=K KP1=K+1 DO 240 J=KP1,NP1 IF (SIGMA .GE. YP(J)) GO TO 240 SIGMA=YP(J) JBAR=J 240 CONTINUE IF (JBAR .EQ. K) GO TO 260 I=PIVOT(K) PIVOT(K)=PIVOT(JBAR) PIVOT(JBAR)=I YP(JBAR)=YP(K) YP(K)=SIGMA DO 250 I=1,N SIGMA=QR(I,K) QR(I,K)=QR(I,JBAR) QR(I,JBAR)=SIGMA 250 CONTINUE C END OF COLUMN INTERCHANGE. 260 SIGMA=DDOT(N-K+1,QR(K,K),1,QR(K,K),1) IF (SIGMA .EQ. 0.0) THEN IFLAG=4 RETURN ENDIF 270 IF (K .EQ. N) GO TO 300 QRKK=QR(K,K) ALPHAK=-SQRT(SIGMA) IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK ALPHA(K)=ALPHAK BETA=1.0/(SIGMA-QRKK*ALPHAK) QR(K,K)=QRKK-ALPHAK DO 290 J=KP1,NP2 SIGMA=BETA*DDOT(N-K+1,QR(K,K),1,QR(K,J),1) DO 280 I=K,N QR(I,J)=QR(I,J)-QR(I,K)*SIGMA 280 CONTINUE IF (J .LT. NP2) YP(J)=YP(J)-QR(K,J)**2 290 CONTINUE 300 CONTINUE ALPHA(N)=QR(N,N) C C C COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS. TZ(NP1)=1.0 DO 340 I=N,1,-1 SUM=0.0 DO 330 J=I+1,NP1 330 SUM=SUM+QR(I,J)*TZ(J) 340 TZ(I)=-SUM/ALPHA(I) YPNORM=DNRM2(NP1,TZ,1) DO 360 K=1,NP1 360 YP(PIVOT(K))=TZ(K)/YPNORM IF (DDOT(NP1,YP,1,YPOLD,1) .GE. 0.0) GO TO 380 DO 370 I=1,NP1 370 YP(I)=-YP(I) C YP IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION. C C COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)). C V IS GIVEN BY P - (P,Q)Q , WHERE P IS ANY SOLUTION OF C [D RHO] V = -RHO AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO]. C 380 DO 440 I=N,1,-1 SUM=QR(I,NP1)+QR(I,NP2) DO 430 J=I+1,N 430 SUM=SUM+QR(I,J)*ALPHA(J) 440 ALPHA(I)=-SUM/ALPHA(I) DO 450 K=1,N 450 TZ(PIVOT(K))=ALPHA(K) TZ(PIVOT(NP1))=1.0 C TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q C IN THE KERNEL(THE TANGENT). SIGMA=DDOT(NP1,TZ,1,YP,1) DO 470 J=1,NP1 TZ(J)=TZ(J)-SIGMA*YP(J) 470 CONTINUE C TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (LAMBDA(S), X(S)). RETURN END