subroutine SWATAN (N5,U5,Y5) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 2001-06-18 SWCOMP Krogh Replaced ". and." with " .and." C>> 1996-04-15 SWCOMP Krogh Changed single common block name to SWCOM. C>> 1995-11-21 SWCOMP Krogh Removed multiple entries for C conversion. C>> 1994-11-11 SWCOMP Krogh Declared all vars. C>> 1994-10-31 SWCOMP Krogh Changes to use M77CON C>> 1987-12-07 SWCOMP Lawson Initial code. c--S replaces "?": ?WCOMP,?WATAN,?WASIN,?WACOS,?WATN2,?WSUM,?WDIF, c--& ?WSQRT,?WEXP,?WSIN,?WCOS,?WTAN,?WSINH,?WCOSH,?WTANH,?WSET,?WSUM1, c--& ?WDIF1,?WPRO1,?WQUO1,?WLOG,?WPWRI,?WCHN,?WPRO,?WQUO,?PASCL, c--& ?WRCHN,?WACSI,?WCOM c C File SWCOMP [- Wengert derivative COMPutation] c The file SWCOMP contains a set of program units to perform c computation on (N+1)-tuples representing the value of a quantity c and its first N derivatives with respect to a single variable. c In the comments we use T as the name of the ultimate independent c variable. c The parameter NMAX must be set in program units c SWATAN, SWSUM, and SWPRO c to establish the largest order of derivative the package c will be able to handle. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c C. L. Lawson, JPL, 1971. C References: c 1. Wengert, R. E., A simple automatic derivative evaluation c program, Comm. ACM, 1, Aug 1964, 463-464. c 2. C. L. Lawson, Computing Derivatives using W-arithmetic and c U-arithmetic., JPL Appl Math TM 289, Sept 1971. c Revised by CLL for Fortran 77 in Jan 1987. Deleted subr WSTART. c Put a first time flag in WPRO. Now user does not need to make c any initialization call. c 9/18/87, CLL. Added new entry names SWASIN & SWACOS. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NMAX Setting NMAX in the parameter statement determines c the highest order derivative this program unit will be able c to handle. The following arrays have dimensions depending c on NMAX: S1(), S2(), S3(), S4(), C() c c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c ERROR HANDLING c c All detected error conditions are essentially fatal for the c requested operation. We do not stop; rather, we issue an error c message and return. By use of the the MATH77 library subroutine, c ERMSET, the user can change the error processing action to cause a c STOP following the error message. c c Error Called c No. name Explanation c c 1 SWASIN Infinite deriv when arg = -1 or +1 c 1 SWACOS Infinite deriv when arg = -1 or +1 c 2 SWSQRT Infinite deriv when arg = 0 c 3 SWQUO1 Zero divisor c 4 SWPWRI U**M is infinite when U = 0 and M < 0 c 5 SWPRO Require dimension NMAX .ge. NDERIV c 6 SWQUO Require dimension NMAX .ge. NDERIV c 7 SWQUO Zero divisor c 8 SPASCL Require N .ge. 0 c 9 SWRCHN Require U(2) .ne. 0. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c PROGRAM UNITS, entry NAMES, and calls c c All of these subroutine and entry names may be called directly by c a user, however, it is expected that SWCHN, SWRCHN, and SPASCL c would only be used to augment the package with new functions. c C Subroutine names: Calls to: c c subroutine SWATAN (N5,U5,Y5) SWPWRI,SWCHN c subroutine SWASIN (N16,U16,Y16) SWACSI c subroutine SWACOS (N16,U16,Y16) SWACSI c subroutine SWACSI (N16,U16,Y16,ACOSIN) SWSQRT,SWPWRI,SWPRO1,SWCHN c subroutine SWATN2 (N9,Y9,X9,A9) SWSUM,SWDIF,SWPRO,SWQUO,SWCHN c c subroutine SWSUM (NDERIV,U,V,W) None c subroutine SWDIF (N2,U2,V2,W2) None c subroutine SWSQRT (N15,U15,Y15) SWCHN c subroutine SWEXP (N4,U4,Y4) SWCHN c subroutine SWSIN (N7,U7,Y7) SWCHN c subroutine SWCOS (N7,U7,Y7) SWCHN c subroutine SWTAN (N7,U7,Y7) SWCHN c subroutine SWSINH (N7,U7,Y7) SWCHN c subroutine SWCOSH (N7,U7,Y7) SWCHN c subroutine SWTANH (N7,U7,Y7) SWCHN, SWPRO, SWQUO c subroutine SWSET (N10,UVAL,UDER,U10) None c subroutine SWSUM1 (N11,C11,U11,Y11) None c subroutine SWDIF1 (N12,C12,U12,Y12) None c subroutine SWPRO1 (N13,C13,U13,Y13) None c subroutine SWQUO1 (N14,C14,U14,Y14) SWQUO c c subroutine SWLOG (NDERIV,U,Y) SWCHN c subroutine SWPWRI (NDERIV,MPWR,U,Y) SWCHN c c subroutine SWCHN (NDERIV,U,F) SWPRO c c subroutine SWPRO (NDERIV,U,V,W) SPASCL c subroutine SWQUO (NDERIV,U,V,W) SPASCL c c subroutine SPASCL (N ,C) None c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Dependencies between program units c c The progran units are SWATAN,SWSUM,SWLOG,SWCHN,SWPRO,SPASCL c SWATAN has calls into SWSUM,SWLOG,SWCHN,SWPRO c SWSUM has calls into SWLOG,SWCHN,SWPRO c SWLOG has calls into SWCHN c SWCHN has calls into SWPRO c SWPRO has calls into SPASCL c SPASCL has no calls c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c integer N5 real U5(*),Y5(*) integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) integer I, NS real S1(NMAXP1) c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C COMPUTE.. Y5=ATAN(U5) and derivs w.r.t. T C Y5(1)= ATAN (U5(1)) IF(N5 .eq. 0) return C C deriv of Y w.r.t. U IS 1./(1. + U**2) C START BY CONSTRUCTING 1.+U**2 AND ITS derivs IN S1(). C S1(1)=1.0E0 + U5(1)**2 NS=N5-1 IF(NS .eq. 0) go to 52 S1(2)= U5(1)+U5(1) IF(NS .eq. 1) go to 52 S1(3)= 2.0E0 IF(NS .eq. 2) go to 52 DO 50 I=3,NS 50 S1(I+1) = 0.0E0 C C NOW S1() = 1.+U**2 and derivs w.r.t. U . C COMPUTE 1./S1 and derivs w.r.t. U. STORE RESULT STARTING IN Y(2) C 52 call SWPWRI (NS,-1,S1,Y5(2)) C C Convert Y from derivs w.r.t. U TO derivs w.r.t. T call SWCHN(N5,U5,Y5) return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWASIN(N16,U16,Y16) C C COMPUTES.. Y = ASIN(U) and derivs w.r.t. T C integer N16 real U16(*), Y16(*) c Y16(1)= ASIN (U16(1)) if (N16 .ne. 0) call SWACSI(N16,U16,Y16,.FALSE.) return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWACOS(N16,U16,Y16) C C COMPUTES.. Y = ACOS(U) and derivs w.r.t. T C integer N16 real U16(*), Y16(*) c Y16(1)= ACOS (U16(1)) if (N16 .ne. 0) call SWACSI(N16,U16,Y16,.TRUE.) return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWACSI(N16,U16,Y16,ACOSIN) c Common code for SWACOS and SWASIN -- ACOSIN - .true. for SWACOS integer N16 real U16(*), Y16(*) logical ACOSIN c integer I, NS integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1),S2(NMAXP1) c C C Deriv of Y w.r.t. U is (+,-)1.0/sqrt(1.0 - U**2) C Start by constructing 1.0 - U**2 and its derivs in S1(). C S1(1)=1.0E0 - U16(1)**2 if(S1(1) .eq. 0.0E0) then c Error condition. if(ACOSIN) then call ERMSG('SWACOS',1,0, * 'Deriv of ACOS(x) is infinite at x = -1 or +1','.') else call ERMSG('SWASIN',1,0, * 'Deriv of ASIN(x) is infinite at x = -1 or +1','.') endif return endif c NS=N16-1 IF(NS .eq. 0) go to 66 S1(2)= -2.0E0 * U16(1) IF(NS .eq. 1) go to 66 S1(3)= -2.0E0 IF(NS .eq. 2) go to 66 DO 64 I=3,NS 64 S1(I+1) = 0.0E0 C C Now S1() = 1.0 - U**2 and derivs w.r.t. U. c Compute S2() = sqrt(1.0 - U**2) and derivs w.r.t. U. C Then change sign if doing Arc Cosine. c 66 continue call SWSQRT(NS,S1,S2) if(ACOSIN) call SWPRO1(NS,-1.0E0,S2,S2) c C Compute 1.0/S2 and derivs w.r.t. U. Store result starting in Y(2) c call SWPWRI (NS,-1,S2,Y16(2)) C C Convert Y from derivs w.r.t. U to derivs w.r.t. T call SWCHN(N16,U16,Y16) return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWATN2(N9,Y9,X9,A9) C C COMPUTES.. A = ATAN2(Y,X) = ATAN(Y/X) and derivs w.r.t. T C integer N9 real Y9(*),X9(*),A9(*) integer I, NM1, NP1 real BIG, XX, YY integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1),S2(NMAXP1),S3(NMAXP1),S4(NMAXP1) c A9(1)= ATAN2(Y9(1),X9(1)) IF(N9 - 1) 90,92,94 90 return C SPECIAL FOR N9 = 1 92 continue BIG= max( ABS(X9(1)), ABS(Y9(1))) YY =Y9(1)/BIG XX =X9(1)/BIG A9(2)=(XX*(Y9(2)/BIG) - YY*(X9(2)/BIG) )/(XX**2 + YY**2) return C GENERAL CASE N9 .ge. 2 94 continue NP1=N9+1 BIG= max( ABS(X9(1)), ABS(Y9(1))) DO 96 I=1,NP1 S1(I)=X9(I)/BIG 96 S2(I)= Y9(I)/BIG NM1=N9-1 call SWPRO (NM1,S1,S2(2),S3) call SWPRO (NM1,S2,S1(2),S4) call SWDIF (NM1,S3,S4,S3) call SWPRO (NM1,S1 ,S1 ,S1 ) call SWPRO (NM1,S2,S2,S2) call SWSUM (NM1,S1 ,S2,S1 ) call SWQUO (NM1,S3,S1 ,A9(2)) return end subroutine SWSUM (NDERIV,U,V,W) C C COMPUTE.. W=U+V and derivs w.r.t. T C integer NDERIV real U(*),V(*),W(*) integer I, NP1 c NP1=NDERIV+1 DO 10 I=1,NP1 10 W(I)=U(I)+V(I) return end C C C subroutine SWDIF (N2,U2,V2,W2) C C COMPUTE.. W2=U2-V2 and derivs w.r.t. T integer N2 real U2(*),V2(*),W2(*) integer I, NP1 c NP1= N2+1 DO 20 I=1,NP1 20 W2(I)=U2(I)-V2(I) return end C C C subroutine SWSQRT (N15,U15,Y15) C C COMPUTE.. Y= SQRT(U) and derivs w.r.t. T C integer N15 real U15(*),Y15(*) integer I real FAC, SCALE, SQRTU integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1) c SQRTU = SQRT(U15(1)) if(N15 .eq. 0) go to 152 if(SQRTU .eq. 0.0E0) then call ERMSG('SWSQRT',2,0, * 'Deriv of sqrt(x) is infinite at x = 0','.') do 155 I = 2,N15+1 Y15(I) = 0.0E0 155 continue return endif if(N15 .gt. 1) go to 154 Y15(2)=U15(2)/(SQRTU+SQRTU) 152 Y15(1)=SQRTU return C GENERAL CASE FOR N .ge. 2 154 continue SCALE=1.0E0/U15(1) Y15(1)=1.0E0 S1(1) =1.0E0 FAC=0.5E0 DO 156 I=1,N15 Y15(I+1)=Y15(I)*FAC FAC=FAC-1.0E0 156 S1(I+1)=U15(I+1)*SCALE call SWCHN(N15,S1,Y15) Y15(1)=SQRTU DO 158 I=1,N15 158 Y15(I+1)=Y15(I+1)*SQRTU return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C subroutine SWEXP (N4,U4,Y4) C C COMPUTE.. Y4=EXP(U4) and derivs w.r.t. T C integer N4 real U4(*),Y4(*) integer I C Y4(1)=EXP(U4(1)) IF(N4 .eq. 0) return DO 40 I=1,N4 40 Y4(I+1)=Y4(I) call SWCHN(N4,U4,Y4) return end C C C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWSIN(N7,U7,Y7) C C COMPUTE.. Y7=SIN(U7) and derivs w.r.t. T C integer N7 real U7(*),Y7(*) integer I C Y7(1) = SIN(U7(1)) IF(N7 .eq. 0) return Y7(2) = COS(U7(1)) if(N7 .eq.1) then Y7(2) = U7(2) * Y7(2) else DO 70 I=2,N7 Y7(I+1)=-Y7(I-1) 70 continue call SWCHN(N7,U7,Y7) end if return end C C C subroutine SWCOS (N7,U7,Y7) C C COMPUTE.. Y7=COS(U7) and derivs w.r.t. T C integer N7 real U7(*),Y7(*) integer I C Y7(1)= COS (U7(1)) IF(N7 .eq. 0) return Y7(2)= -SIN(U7(1)) if(N7 .eq.1) then Y7(2) = U7(2) * Y7(2) else DO 70 I=2,N7 Y7(I+1)=-Y7(I-1) 70 continue call SWCHN(N7,U7,Y7) end if return end c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWTAN (N7,U7,Y7) C C COMPUTE.. Y7=TAN(U7) and derivs w.r.t. T C The first deriv is 1/cos(u)**2 c integer N7 real U7(*),Y7(*) integer I, N7M1 real CS, SN integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1) c if(N7 .eq. 0) then Y7(1) = tan(U7(1)) return end if SN = SIN(U7(1)) CS = COS(U7(1)) Y7(1)= SN/CS c c Compute Cos(U) and its derivs w.r.t. U in S1(2..N+1) c S1(2) = CS S1(3) = -SN DO 76 I=4,N7+1 76 S1(I) = -S1(I-2) N7M1 = N7-1 c Convert to derivs w.r.t. t call SWCHN(N7M1,U7,S1(2)) c Compute Cos(U)**2 & derivs w.r.t. t call SWPRO(N7M1, S1(2), S1(2), S1(2)) c Divide first deriv of U by Cos(U)**2, c both with derivs w.r.t. t call SWQUO(N7M1, U7(2), S1(2), Y7(2)) return end c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWSINH(N7,U7,Y7) C C COMPUTE.. Y7=SINH(U7) and derivs w.r.t. T C integer N7 real U7(*),Y7(*) integer I C Y7(1) = SINH(U7(1)) IF(N7 .eq. 0) return Y7(2) = COSH(U7(1)) if(N7 .eq.1) then Y7(2) = Y7(2) * U7(2) else DO 80 I=2,N7 Y7(I+1) = Y7(I-1) 80 CONTINUE call SWCHN(N7,U7,Y7) end if return end C C C subroutine SWCOSH (N7,U7,Y7) C C COMPUTE.. Y7=COSH(U7) and derivs w.r.t. T C integer N7 real U7(*),Y7(*) integer I C Y7(1)= COSH (U7(1)) IF(N7 .eq. 0) return Y7(2)= SINH(U7(1)) if(N7 .eq.1) then Y7(2) = Y7(2) * U7(2) else DO 80 I=2,N7 Y7(I+1) = Y7(I-1) 80 CONTINUE call SWCHN(N7,U7,Y7) end if return end C C C subroutine SWTANH (N7,U7,Y7) C C COMPUTE.. Y7=TANH(U7) and derivs w.r.t. T C The first deriv is 1/cosh(u)**2 c integer N7 real U7(*),Y7(*) integer I, N7M1 real CH, SH integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1) C if(N7 .eq. 0) then Y7(1) = tanh(U7(1)) return endif SH = SINH(U7(1)) CH = COSH(U7(1)) Y7(1)= SH/CH c c Compute Cosh(U) and its derivs w.r.t. U in S1(2..N+1) c S1(2) = CH S1(3) = SH DO 94 I=4,N7+1 94 S1(I) = S1(I-2) N7M1 = N7-1 c Convert to derivs w.r.t. t call SWCHN(N7M1,U7,S1(2)) c Compute Cosh(U)**2 & derivs w.r.t. t call SWPRO(N7M1, S1(2), S1(2), S1(2)) c Divide first deriv of U by Cosh(U)**2, c both with derivs w.r.t. t call SWQUO(N7M1, U7(2), S1(2), Y7(2)) return end C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWSET (N10,UVAL,UDER,U10) C C INITIALIZE THE ARRAY U10().. SET VALUE=UVAL, 1-ST DERIV = UDER, C AND HIGHER derivs = ZERO integer N10 real U10(*), UVAL, UDER integer I C IF(N10-1) 104,102,100 100 DO 101 I=2, N10 101 U10(I+1)= 0.0E0 102 U10(2)=UDER 104 U10(1)=UVAL return end C C C subroutine SWSUM1 (N11,C11,U11,Y11) C C COMPUTE.. Y=C+U WHERE C IS A SCALAR C integer N11 real C11, U11(*), Y11(*) integer I C Y11(1)=C11+U11(1) DO 110 I=1,N11 110 Y11(I+1)=U11(I+1) return end C C C subroutine SWDIF1 (N12,C12,U12,Y12) C C COMPUTE.. Y=C-U WHERE C IS A SCALAR C integer N12 real C12, U12(*), Y12(*) integer I C Y12(1)=C12-U12(1) DO 120 I=1,N12 120 Y12(I+1)=-U12(I+1) return end C C C subroutine SWPRO1(N13,C13,U13,Y13) C C COMPUTE.. Y=C*U WHERE C IS A SCALAR C integer N13 real C13, U13(*), Y13(*) integer I C DO 130 I=1,N13+1 130 Y13(I) = C13*U13(I) return end C C C subroutine SWQUO1(N14,C14,U14,Y14) C C COMPUTE.. Y=C/U WHERE C IS A SCALAR C integer N14 real C14, U14(*), Y14(*) integer I integer NMAX, NMAXP1 parameter (NMAX=10, NMAXP1=NMAX+1) real S1(NMAXP1) C if(U14(1) .eq. 0.0E0) then call ERMSG('SWQUO1',3,0, 'Zero divisor.','.') return endif S1(1) = C14 do 160 I = 2, N14+1 S1(I) = 0.0E0 160 continue call SWQUO (N14,S1,U14,Y14) return end subroutine SWLOG (NDERIV,U,Y) C C. L. LAWSON, JPL, 1971 NOV 19 C C COMPUTE.. Y= LOG BASE E OF U and derivs w.r.t. T C ERROR STOP IN this subr IF U(1) .le. 0. C integer II, NDERIV real U(*),Y(*) real FAC, UINV C Y(1)= LOG(U(1)) if (NDERIV .le. 1) then if (NDERIV .eq. 1) Y(2)= U(2)/U(1) return end if UINV=1.0E0/U(1) Y(2)=UINV FAC=-UINV DO 120 II = 2, NDERIV Y(II+1)=Y(II)*FAC FAC=FAC-UINV 120 continue C C Y() NOW CONTAINS VALUE and derivs w.r.t. U C USE CHAIN RULE TO convert TO derivs w.r.t. T C call SWCHN(NDERIV,U,Y) return end C subroutine SWPWRI (NDERIV,MPWR,U,Y) C C COMPUTE.. Y = U**MPWR and derivs w.r.t. T C MPWR is an integer, POS.,NEG.,OR ZERO. C MPWR does not depend ON T. C U may depend on T C Issues an error message if both (U(1) .eq. 0.) and (MPWR .lt. 0) C If MPWR .eq. 0 then set Y(1) = 1. and all derivs = 0. C integer I, II, II1, II2, J, M, MPWR, NDERIV real U(*),Y(*) logical SHORT real D, FAC, UINV C II2=NDERIV SHORT = .FALSE. M=MPWR IF(M ) 110,35,45 C HERE MPWR .eq. 0 35 Y(1)= 1.0E0 IF(NDERIV .le. 0) return DO 40 J=1,NDERIV 40 Y(J+1) = 0.0E0 return C HERE MPWR IS POSITIVE 45 IF(U(1) .ne. 0.0E0) go to 50 Y(1)=0.0E0 IF(NDERIV .eq. 0) return IF(M .le. NDERIV) go to 46 II=1 go to 130 C 46 DO 47 I=1,NDERIV 47 Y(I+1) = 0.0E0 FAC=1.0E0 D =0.0E0 DO 48 I=1,M D=D+1.0E0 48 FAC=FAC*D Y(M+1) = FAC go to 150 C 50 IF( M .ge. NDERIV) go to 112 II2=M SHORT = .TRUE. go to 112 C C HERE MPWR .lt 0 110 continue if(U(1) .eq. 0.0E0) then call IERM1('SWPWRI',4,0, * 'U**M is infinite when U = 0. and M < 0', * 'M',M,'.') do 111 I = 1,NDERIV+1 Y(I) = 0.0E0 111 continue return endif 112 continue Y(1)=U(1)**M if(NDERIV .eq. 0) return UINV=1.0E0/U(1) FAC=UINV*M II1=1 DO 120 II=II1,II2 Y(II+1)=Y(II)*FAC FAC=FAC-UINV 120 continue if(.NOT. SHORT) go to 150 II=II2+1 C C SET HIGHER derivs TO ZERO. C 130 continue DO 140 J=II,NDERIV 140 Y(J+1)=0.0E0 C C Y() NOW CONTAINS VALUE and derivs w.r.t. U C USE CHAIN RULE TO convert TO derivs w.r.t. T C 150 call SWCHN(NDERIV,U,Y) return end subroutine SWCHN(NDERIV,U,Y) C c Implements the chain rule. c Given values of y(u) and its derivs with respect to u in Y(), and c u(x) and its derivs with respect to x in U(), this subr replaces c the contents of Y() with y(u(x)) and its derivs with respect to x. c C NDERIV HIGHEST ORDER DERIVATIVE c C (U(I),I=1,NDERIV+1) INPUT.. value of u and derivs w.r.t. x C C (Y(I),I=1,NDERIV+1) INPUT.. value of y and derivs w.r.t. u C OUTPUT.. (Y(I+1),I=1,NDERIV) will be C replaced by derivs w.r.t. x integer I, L, NDERIV, NP1 real U(*), Y(*) C if(NDERIV .eq. 0) return if(U(2) .ne. 1.) go to 20 C U(2) .eq. 1., See if higher derivs are 0. if(NDERIV .eq. 1) return DO 10 I=2,NDERIV if(U(I+1) .ne. 0.) go to 50 10 continue return C Test if all derivs of U are 0. 20 DO 30 I=1,NDERIV if(U(I+1) .ne. 0.) go to 50 30 continue C DO 40 I=1,NDERIV 40 Y(I+1)=0. return C C Here for the general case when U is not T and not constant. C 50 NP1=NDERIV+1 DO 60 L=0,NDERIV-1 60 call SWPRO (L, Y(NP1-L), U(2), Y(NP1-L)) return end subroutine SWRCHN(NDERIV, U, Y) C c Implements the "reverse" chain rule. c Given values of y(u(x)) and its derivs with respect to x in Y(), c and u(x) and its derivs with respect to x in U(), this subr replac c the contents of Y() with y(u) and its derivs with respect to u. c Requires u'(x) .ne. 0. c c Note that this subr can be used to compute a representation of c the inverse function of a given function. This transformation c is called "reversion" of a power series. To do this let y() be c the inverse function of u() in a neighborhood of a point, x0. c Thus y(u(x)) = x for all x in a neighborhood of x0. c Then the value of y() and its derivs with respect to x, evaluated c at x0 are (x0, 1.0, 0.0, ..., 0.0). c If Y() contains these values on entry, and U() contains u(x) and c its derivs w.r.t. x, evaluated at x0, then on return Y() will c contain y(u) and its derivs w.r.t. u, evaluated at u0 = u(x0). c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NDERIV [in] Highest order derivative to be considered. C C (U(I),I=1,NDERIV+1) [in] Contains the value of u(x) and its c first NDERIV derivatives with respect to x, evaluated at x0. c Require U(2) .ne. 0. c C (Y(I),I=1,NDERIV+1) [inout] On entry, Y() must contain y(u(x)) c and its first NDERIV derivatives w.r.t. x, evaluated at x0. c On return, Y() will contain y(u) and its first NDERIV c derivatives w.r.t. u, evaluated at u0 = u(x0). c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - integer NDERIV, NP1, L real U(*), Y(*) c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - NP1=NDERIV+1 c if( U(2) .eq. 0.0E0 .and. NP1 .ge. 2) then call ERMSG('SWRCHN',9,0,'U(2) is zero.','.') return endif c do 60 L=2, NP1 CALL SWQUO (NP1 - L, Y(L), U(2), Y(L)) 60 continue return end subroutine SWPRO(NDERIV,U,V,W) C COMPUTE W(1)=U(1)*V(1) AND DERIVATIVES . C NDERIV HIGHEST ORDER DERIVATIVE TO BE COMPUTED C NDERIV .ge. 0 C (U(I),I=1,NDERIV+1) INPUT ARRAY.. C U(1)= VALUE C U(I+1)= I-TH DERIVATIVE OF U(1) C (V(I),I=1,NDERIV+1) INPUT ARRAY.. SAME FORMAT AS U(). C C (W(I),I=1,NDERIV+1) OUTPUT ARRAY.. C W(1)= U(1)*V(1) C W(I+1)=Ith deriv of W(1) c It is permissible for W() to occupy the same storage locations as c U() and/or V(). c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NMAX Setting NMAX in the parameter statement determines c the highest order derivative this program unit will be able c to handle. The dimension of the internal saved array C() c will be set as a function of NMAX. C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - integer I, ID2, IMID, IP1, J, JC, NDERIV, NP1 real U(*),V(*),W(*) integer NMAX, NCOEF parameter(NMAX = 10, NCOEF = 1 + (NMAX*(NMAX+1))/2 ) real C(NCOEF) real TEMP logical FIRST common / SWCOM / C save FIRST, /SWCOM/ data FIRST /.TRUE. / C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(FIRST) then call SPASCL(NMAX, C) FIRST = .FALSE. endif if( NDERIV .gt. NMAX) then call IERM1('SWPRO',5,0, * 'Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV', * 'NMAX',NMAX,',') call IERV1('NDERIV',NDERIV,'.') return endif NP1=NDERIV+1 c JC = (NDERIV*(NDERIV-1))/2 c DO 20 I=NP1,1,-1 c TEMP=0.0E0 c IP1=I+1 c DO 10 J=1,I c TEMP = TEMP + U(J)*V(IP1-J)*C(JC+J) c 10 continue c W(I) = TEMP c JC=JC- I+2 c 20 continue c c The following code does the same as the above commented-out c code but reduces the number of multiplies. c Treats the cases of NDERIV = 0, 1, 2, 3, and 4 special for c efficiency, since it is expected that most usage will involve c small values of NDERIV. In particular the chain rule subroutine, c _WCHN, will call this subr with a sequence of values of NDERIV c going down to 0. c Treats NDERIV > 4 in a general way, but takes account of c the fact that the C() terms in the "do 10" loop above are c symmetric about the middle term and the first and last term c are each = 1. c go to (41,42,43,44,45), NP1 c c Fortran 77 specifies that a Computed Go To drops through is the c index is out of range. Thus we arrive here if NP1 < 1 or > 5. c NP1 < 1 would be an error. We choose not to use time to test for c it however. c NP1 > 5 is valid. In that case we do the following loop for I c running down from NP1 to 6, and then do the special statements c for 5, 4, 3, 2, and 1. c JC = (NDERIV*(NDERIV-1))/2 DO 20 I=NP1,6,-1 TEMP = U(1) * V(I) + U(I) * V(1) IP1=I+1 ID2 = I/2 DO 10 J=2,ID2 TEMP = TEMP + C(JC+J) * (U(J)*V(IP1-J) + U(IP1-J)*V(J)) 10 continue if(2*ID2 .eq. I) then W(I) = TEMP else IMID = ID2 + 1 W(I) = TEMP + C(JC+IMID) * U(IMID) * V(IMID) endif JC=JC- I+2 20 continue c 45 W(5) = U(1)*V(5) + 4.0E0*(U(2)*V(4) + U(4)*V(2)) + * 6.0E0*U(3)*V(3) + U(5)*V(1) 44 W(4) = U(1) * V(4) + 3.0E0*(U(2)*V(3) + U(3) * V(2)) + U(4)*V(1) 43 W(3) = U(1) * V(3) + 2.0E0 * U(2) * V(2) + U(3) * V(1) 42 W(2) = U(1) * V(2) + U(2) * V(1) 41 W(1) = U(1) * V(1) return end c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine SWQUO (NDERIV, U ,V ,W ) C C COMPUTE.. W = U / V and derivs w.r.t. T c c V() must be distinct in storage from W(), but U() can be the same c array as W(). c C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - integer NDERIV real U(*),V(*),W(*) integer I, ID2, IMID, IP1, J, JC, NP1 integer NMAX, NCOEF parameter(NMAX = 10, NCOEF = 1 + (NMAX*(NMAX+1))/2 ) real C(NCOEF) real FAC, TEMP logical FIRST common / SWCOM / C save FIRST, /SWCOM/ data FIRST /.TRUE. / C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(FIRST) THEN CALL SPASCL(NMAX, C) FIRST = .FALSE. END IF if( NDERIV .GT. NMAX) then call IERM1('SWQUO',6,0, * 'Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV', * 'NMAX',NMAX,',') call IERV1('NDERIV',NDERIV,'.') return endif c if( V(1) .eq. 0.0E0) then call ERMSG('SWQUO',7,0,'Zero divisor.','.') return endif c FAC = 1.0E0/V(1) W(1) = U(1) * FAC NP1=NDERIV+1 JC = 0 c DO 40 I=2,NP1 c TEMP = U(I) c IP1=I+1 c DO 30 J=2,I c TEMP = TEMP - V(J)*W(IP1-J)*C(JC+J) c 30 CONTINUE c W(I) = TEMP * FAC c JC=JC + I-1 c 40 CONTINUE DO 40 I= 2, NP1 TEMP = U(I) - W(1) * V(I) IP1=I+1 ID2 = I/2 DO 30 J=2,ID2 TEMP = TEMP - C(JC+J) * (W(J)*V(IP1-J) + W(IP1-J)*V(J)) 30 continue if(2*ID2 .ne. I) then IMID = ID2 + 1 TEMP = TEMP - C(JC+IMID) * W(IMID) * V(IMID) endif W(I) = TEMP * FAC JC=JC + I-1 40 continue return end subroutine SPASCL(N ,C) C C.L.LAWSON,JPL, 1969 DEC 10 C C N MAXIMUM ORDER DERIVATIVE TO BE C COMPUTED USING W-ARITHMETIC C C (C(I),I=1,NC) FIRST N+1 ROWS OF PASCAL TRIANGLE PACKED C WITH 1-ST DIAGONAL OMITTED AFTER 1-ST ROW. C NC= (N*(N+1)/2) + 1 C EXAMPLE.. IF N=3 then TRIANGLE IS 1 C 1 1 C 1 2 1 C 1 3 3 1 C C WHICH WILL BE STORED AS.. 1,1,2,1,3,3,1 C integer I, J, K, L, N logical DONE real C(*) save DONE data DONE / .false. / C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if (DONE) return DONE = .true. if(N .lt. 0) then call IERM1('SPASCL',8,0, 'Require N .ge. 0','N',N,'.') return endif C(1)=1.0E0 if(N .eq. 0) return C(2)=1.0E0 if(N .eq. 1) return K=3 DO 520 I=1,N-1 L=I+1 DO 510 J=1,I C(K)=C(K-L)+C(K-L+1) K=K+1 510 continue C(K)=1.0E0 K=K+1 520 continue return end