**** options **** jprobl=1 square 2 notch cut out of square 5 squeezed square 6 quarter circle 7 quarter of circular fin waveguide 1,1 11 quarter of circular fin waveguide general jguess=1 the trivial guess 0. (all cases) 6 an excellent guess for the notch 8 guess based on arg(wt(i)), beta=0 9 assumes x(1),...,x(2*n-4) set in main.f nrml = 1 0 -> 0, wn -> wtn, side 1 -> circle of side_t 1 3 optimize vertex errors, starting from nrml=1 jeqn = 2 re, im wn = wtn 5 sigma(j) = sigmat(j), 1<=j<=n-1 **** program outline **** main sincp sincp2 scrss n2f sfneq sinteg sstpar selim strans sray ode sfode smtov scirc sle2 snewt2 sinter scrss smapc srinv srarg scplot ----------------------------------------------------------- scirc find radius and orientation from derivatives scplot plot polygon scrss search for parameter in normalizing transformation selim solve for 3 betas sexprl relative exponential sfneq define nonlinear system sfode define ode system sincp initialize problem and parameters sincp2 find auxiliary information from initial data sinteg given parameters, find computed polygon sinter nonlinear system for intersection of circles sle2 2 by 2 linear systems smapc find center of Mobius image of circle smobs apply implicitly defined Mobius transformation smtov find vertices from midpoint data snewt2 newton's method srarg compute relative angle about center sray integrate ode to midpoints of sides srinv zero-tolerant divide sstpar store new guess at parameters strans transform to constrained variables ----------------------------------------------------------- current target w wt vertices wc wtc centers wdir wtdir orientations side j connects w(j-1) and w(j) boundary is traced counterclockwise around polygon wdir = 1 for convex polygons miscellaneous common variables: epsfcn estimate of error in computing sfneq common variables no longer used: s punp ptrace hybeps c n number of sides c wtdir orientations (convex=1, concave=-1) c w0 image of 0 inside polygon c wt vertices c wtm point on side wt(n) -- wt(1) c wtc centers of sides c ****** conformal mapping of circular arc polygons ****** c Petter Bjorstad and Eric Grosse 7 Apr 1986 c Siam J Sci Stat Comp Apr 1987? c c CAUTION. This is an experimental research code. We certainly c cannot promise you it will always converge quickly! Please c let us know what troubles you have; although there are no c immediate plans to turn it into good mathematical software, c we would certainly welcome any improvements. c c The routine N2F called here is availble by asking netlib: c send n2f ivset from port c real x(16) c common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) common /plotcm/ punp, ptrace, mtrace, ntrace, odet integer punp, ptrace, mtrace, ntrace complex odet(100,10) common /method/ nrml, jeqn integer nrml, jeqn common /wavegd/ cfws, cfwt real cfws, cfwt c integer m, neqn, nl2iv(2000) c 82*m real nl2v(3000) c 105 + m*(neqn+2*m+17) + 2*neqn complex z2,wz external sfode, sfneq real y(6), r1, r2, relerr, abserr, awork(226) integer iwork(5),iflag integer i, j, k, nfun, niter c integer mxfcal, mxiter, afctol, rfctol, dltfdj, sctol integer xctol, delta0, dltfdc, lmax0, covprt, rdreq parameter ( mxfcal=17, mxiter=18, afctol=31, rfctol=32 ) parameter ( dltfdj=43, sctol=37, xctol=33, delta0=44 ) parameter ( dltfdc=42, lmax0=35, covprt=14, rdreq=57 ) c cfws=1.1 cfwt=1 c from converged 1,1 x( 1) = -4.412843E-01 x( 2) = -9.410260E-01 x( 3) = 1.542326E+00 x( 4) = -2.231437E+00 x( 5) = 2.489092E-02 x( 6) = -2.438678E-01 call sincp(x,2,1,1,2) m = 2*n-4 c if(jeqn.eq.1)then neqn = 2*n-1 else if(jeqn.eq.2)then neqn = 2*n-2 else if(jeqn.eq.3)then neqn = 3*n-2 else if(jeqn.eq.4)then neqn = n-1 else if(jeqn.eq.5)then neqn = n-1 else write(6,1119) jeqn 1119 format(' equation schema ',i5,' is not implemented') stop endif c it = 0 call ivset(1,nl2iv,2000,3000,nl2v) nl2iv(mxfcal) = maxfun nl2iv(mxiter) = 50 c nl2iv(covprt) = 0 c nl2iv(rdreq ) = 0 nl2v(lmax0 ) = dmax nl2v(afctol) = epsfcn**2 nl2v(rfctol) = epsfcn nl2v(sctol ) = epsfcn nl2v(dltfdj) = sqrt(epsfcn) nl2v(xctol ) = max(0.1*eps**(0.5),10*epsfcn) nl2v(delta0) = sqrt(epsfcn) nl2v(dltfdc) = (epsfcn)**(1./3.) call n2f(neqn,m,x,sfneq,nl2iv,2000,3000,nl2v,nl2iv,nl2v,sfode) c final function evaluation to be sure common block values are c up to date call sfneq(neqn,m,x,i,nl2v,nl2iv,nl2v,sfode) c write(6,*) 'hessian=' c k=nl2iv(56)-1 c do 569 i = 1, m c write(6,*) (nl2v(j+k),j=1,i) c k = k + i c 569 continue do 10 i = 1, m write (6,1003) i, x(i) 1003 format (' x(', i2, ') =', 1pe14.6) 10 continue abserr = 0 do 15 i = 1, n abserr = max(abserr,abs(w(i)-wt(i))) 15 continue write(6,1005) abserr 1005 format('max abs vertex err =',1pe10.2) do 20 i = 1, n write (6,1004) i, beta(i) 1004 format ('beta(', i2, ') =', 1pe14.6) 20 continue call resist(n,x) c stop end C subroutine scplot(wb, wa, wd) C integer wd(10) C complex wb(10), wa(10) C end Compliments of netlib Fri Apr 19 13:21:51 EST 1985 c modified to flag finite difference calls to calcr by ifidif=1 SUBROUTINE N2F(N, P, X, CALCR, IV, LIV, LV, V, 1 UIPARM, URPARM, UFPARM) common /scplt/ ifidif C C *** MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY.. C *** THIS AMOUNTS TO N2G WITHOUT THE SUBROUTINE PARAMETER CALCJ. C C *** PARAMETERS *** C INTEGER N, P, LIV, LV C/6 INTEGER IV(LIV), UIPARM(1) REAL X(P), V(LV), URPARM(1) C/7 C INTEGER IV(LIV), UIPARM(*) C REAL X(P), V(LV), URPARM(*) C/ EXTERNAL CALCR, UFPARM C C----------------------------- DISCUSSION ---------------------------- C C THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO CALL C RN2G. C THE PARAMETERS FOR N2F ARE THE SAME AS THOSE FOR N2G C (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING C CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, N2F COMPUTES C AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE C V(DLTFDJ) BELOW. N2F USES FUNCTION VALUES ONLY WHEN COMPUT- C THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS C THAT N2G MAY USE). TO DO SO, N2F SETS IV(COVREQ) TO MINUS C ITS ABSOLUTE VALUE. THUS V(DELTA0) IS NEVER REFERENCED AND ONLY C V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC). C THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO- C BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION C COUNT IV(NFCALL), BUT ARE RECORDED IN IV(NGCALL) INSTEAD. C C V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE C FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN- C VOLVING X(I), THE STEP SIZE FIRST TRIED IS C V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)), C WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF C THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN C SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE- C LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF. C DEFAULT = MACHEP**0.5. C C *** REFERENCE *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. C SOFTWARE, VOL. 7, NO. 3. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL IVSET, RN2G, N2RDP C C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C RN2G... CARRIES OUT OPTIMIZATION ITERATIONS. C N2RDP... PRINTS REGRESSION DIAGNOSTICS. C C *** LOCAL VARIABLES *** C INTEGER D1, DK, DR1, I, IV1, J1K, K, N1, N2, NF, NG, RD1, R1, RN REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO C C *** IV AND V COMPONENTS *** C INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL, 1 NGCALL, NGCOV, R, REGD, REGD0, TOOBIG, VNEED C/6 DATA COVREQ/15/, D/27/, DINIT/38/, DLTFDJ/43/, J/70/, MODE/35/, 1 NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, NGCOV/53/, 2 R/61/, REGD/67/, REGD0/82/, TOOBIG/2/, VNEED/4/ C/7 C PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35, C 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, C 2 R=61, REGD=67, REGD0=82, TOOBIG=2, VNEED=4) C/ DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V) IV(COVREQ) = -IABS(IV(COVREQ)) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+2) CALL RN2G(X, V, IV, LIV, LV, N, N, N1, N2, P, V, V, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + N IV(J) = IV(REGD0) + N IV(NEXTV) = IV(J) + N*P IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RN = R1 + N - 1 RD1 = IV(REGD0) C 20 CALL RN2G(V(D1), V(DR1), IV, LIV, LV, N, N, N1, N2, P, V(R1), 1 V(RD1), V, X) IF (IV(1)-2) 30, 50, 100 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) ifidif = 0 CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** C C *** INITIALIZE D IF NECESSARY *** C 50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO) 1 CALL V7SCP(P, V(D1), ONE) C J1K = DR1 DK = D1 NG = IV(NGCALL) - 1 IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1 DO 90 K = 1, P XK = X(K) H = V(DLTFDJ) * AMAX1( ABS(XK), ONE/V(DK)) H0 = H DK = DK + 1 60 X(K) = XK + H NF = IV(NFGCAL) ifidif = 0 CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM) NG = NG + 1 IF (NF .GT. 0) GO TO 70 H = NEGPT5 * H IF ( ABS(H/H0) .GE. HLIM) GO TO 60 IV(TOOBIG) = 1 IV(NGCALL) = NG GO TO 20 70 X(K) = XK IV(NGCALL) = NG DO 80 I = R1, RN V(J1K) = (V(J1K) - V(I)) / H J1K = J1K + 1 80 CONTINUE 90 CONTINUE GO TO 20 C 100 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 CALL N2RDP(IV, LIV, LV, N, V(RD1), V) C 999 RETURN C C *** LAST CARD OF N2F FOLLOWS *** END c name changed for use inside sinteg Compliments of netlib Fri Apr 19 13:21:51 EST 1985 c modified to flag finite difference calls to calcr by ifidif=1 SUBROUTINE N2F2(N, P, X, CALCR, IV, LIV, LV, V, 1 UIPARM, URPARM, UFPARM) common /scplt/ ifidif C C *** MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY.. C *** THIS AMOUNTS TO N2G WITHOUT THE SUBROUTINE PARAMETER CALCJ. C C *** PARAMETERS *** C INTEGER N, P, LIV, LV C/6 INTEGER IV(LIV), UIPARM(1) REAL X(P), V(LV), URPARM(1) C/7 C INTEGER IV(LIV), UIPARM(*) C REAL X(P), V(LV), URPARM(*) C/ EXTERNAL CALCR, UFPARM C C----------------------------- DISCUSSION ---------------------------- C C THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO CALL C RN2G. C THE PARAMETERS FOR N2F ARE THE SAME AS THOSE FOR N2G C (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING C CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, N2F COMPUTES C AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE C V(DLTFDJ) BELOW. N2F USES FUNCTION VALUES ONLY WHEN COMPUT- C THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS C THAT N2G MAY USE). TO DO SO, N2F SETS IV(COVREQ) TO MINUS C ITS ABSOLUTE VALUE. THUS V(DELTA0) IS NEVER REFERENCED AND ONLY C V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC). C THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO- C BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION C COUNT IV(NFCALL), BUT ARE RECORDED IN IV(NGCALL) INSTEAD. C C V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE C FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN- C VOLVING X(I), THE STEP SIZE FIRST TRIED IS C V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)), C WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF C THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN C SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE- C LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF. C DEFAULT = MACHEP**0.5. C C *** REFERENCE *** C C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. C SOFTWARE, VOL. 7, NO. 3. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL IVSET, RN2G, N2RDP C C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C RN2G... CARRIES OUT OPTIMIZATION ITERATIONS. C N2RDP... PRINTS REGRESSION DIAGNOSTICS. C C *** LOCAL VARIABLES *** C INTEGER D1, DK, DR1, I, IV1, J1K, K, N1, N2, NF, NG, RD1, R1, RN REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO C C *** IV AND V COMPONENTS *** C INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL, 1 NGCALL, NGCOV, R, REGD, REGD0, TOOBIG, VNEED C/6 DATA COVREQ/15/, D/27/, DINIT/38/, DLTFDJ/43/, J/70/, MODE/35/, 1 NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, NGCOV/53/, 2 R/61/, REGD/67/, REGD0/82/, TOOBIG/2/, VNEED/4/ C/7 C PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35, C 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, C 2 R=61, REGD=67, REGD0=82, TOOBIG=2, VNEED=4) C/ DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/ C C--------------------------------- BODY ------------------------------ C write(6,*) 'n2f2{' IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V) IV(COVREQ) = -IABS(IV(COVREQ)) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+2) CALL RN2G(X, V, IV, LIV, LV, N, N, N1, N2, P, V, V, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + N IV(J) = IV(REGD0) + N IV(NEXTV) = IV(J) + N*P IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RN = R1 + N - 1 RD1 = IV(REGD0) C 20 CALL RN2G(V(D1), V(DR1), IV, LIV, LV, N, N, N1, N2, P, V(R1), 1 V(RD1), V, X) IF (IV(1)-2) 30, 50, 100 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) ifidif = 0 CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** C C *** INITIALIZE D IF NECESSARY *** C 50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO) 1 CALL V7SCP(P, V(D1), ONE) C J1K = DR1 DK = D1 NG = IV(NGCALL) - 1 IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1 DO 90 K = 1, P XK = X(K) H = V(DLTFDJ) * AMAX1( ABS(XK), ONE/V(DK)) H0 = H DK = DK + 1 60 X(K) = XK + H NF = IV(NFGCAL) ifidif = 0 CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM) NG = NG + 1 IF (NF .GT. 0) GO TO 70 H = NEGPT5 * H IF ( ABS(H/H0) .GE. HLIM) GO TO 60 IV(TOOBIG) = 1 IV(NGCALL) = NG GO TO 20 70 X(K) = XK IV(NGCALL) = NG DO 80 I = R1, RN V(J1K) = (V(J1K) - V(I)) / H J1K = J1K + 1 80 CONTINUE 90 CONTINUE GO TO 20 C 100 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 CALL N2RDP(IV, LIV, LV, N, V(RD1), V) C 999 continue write(6,*) 'n2f2}' return END subroutine resist ( n, x ) integer n real x(25), ang(25), ang2(25) complex z(25) integer i,nm real k, r complex a,b,c,d,chi c call stransa(x,z,ang,n) write(6,1001)(i,ang(i),i=1,n) 1001 format('tau(',i2,')=',f10.6) write(6,1002)(i,z(i),i=1,n) 1002 format('z(',i2,')=',1p2e16.8) c do 190 i=1,n-3 do 180 j=i+2,n-1 a = z(i) b = z(i+1) c = z(j) d = z(j+1) chi = (b-a)*(d-c)/((c-b)*(a-d)) k = real(chi) k = 1 - 2*(k+sqrt(k*k-k)) r = 2*ellip(k)/ellip(sqrt(1-k*k)) write(6,1004) i, j, r 180 continue 190 continue 1004 format('resistance from side',i2,' to side ',i2,' = ',1pe16.8) return end c------------------------------------------------- real function ellip ( k ) c complete elliptic integral c from Hart, Computer Approximations, 154 real eps, k, a, b, t, pi pi = 4*atan(1.) eps = (r1mach(4)) a = 1 b = sqrt(1-k*k) 10 continue t = (a+b)/2 b = sqrt(a*b) a = t if(abs(a-b).gt.eps*(a+b)) goto 10 ellip = pi/(2*a) return end function sarc ( z1, z2, c, wdir ) c c arc length between points z1 and z2 on circle centered at c c wdir = 1 counterclockwise from z1 to z2 c -1 clockwise c real sarc complex z1, z2, c, c1, t integer wdir real r, s, pi pi = 4.*atan(1.) r = .5*(abs(z1-c)+abs(z2-c)) s = abs(z1-z2) c phi = angle counterclockwise from z1-c to z2-c t = z1 - c t = (z2-c)/t phi = atan2(aimag(t),real(t)) if(phi.lt.0) phi = phi + 2*pi if(phi.eq.0.0)then t = c-z1 c1 = .5*(z1+z2) + 30*max(abs(z1),abs(z2))*t/abs(t) t = z1 - c1 t = (z2-c1)/t phi = atan2(aimag(t),real(t)) if(phi.lt.0) phi = phi + 2*pi endif s = min(1.,max(-1.,.5*s/r)) if( (wdir.eq. 1 .and. phi.le.pi) .or. + (wdir.eq.-1 .and. phi.gt.pi) )then sarc = 2 * r * asin(s) else sarc = 2 * r * (pi-asin(s)) end if return end function oarc ( z1, z2, c) c c arc length between points z1 and z2 on circle centered at c c defined in some magic way ?!?! c real oarc complex z1, z2, c real pi, arg1, arg2, t, t0 pi=4*atan(1.e0) arg1= atan2(aimag(z1-c), real(z1-c)) arg2= atan2(aimag(z2-c), real(z2-c)) t=arg1-arg2 if(abs(z1-z2).le.0.1e0*abs(z2-c)) t=asin(aimag((z1-z2)/(z2-c))) t=mod(t+3*pi,2*pi)-pi if(abs(t).lt.0.01)goto 100 oarc=abs(z2-c)*t return 100 t0=t t=t*t oarc=abs(z1-z2)/sqrt(1-(t/(3*4))*(1-(t/(5*6))*(1-(t/(7*8))*(1- $ (t/(9*10))*(1-(t/(11*12))*(1-(t/(13*14)))))))) if(t0.lt.0) oarc = - oarc return end subroutine scirc(z, w, w1, w2, c, r, orient) c integer orient real r complex z, w, w1, w2, c c c this routine computes the radius, and the c orientation of the arc from derivative information. c real phidot, w11r, w11i, w11s, w22r, w22i complex z1, w11, w22 c z1 = z* cmplx(0.e0, 1.e0) w11 = w1*z1 w11r = real(w11) w11i = aimag(w11) w11s = w11r**2+w11i**2 w22 = -z*(z*w2+w1) w22r = real(w22) w22i = aimag(w22) phidot = w11r*w22i-w11i*w22r orient = sign(1.e0,phidot) r = abs(phidot)/w11s r = 1.e0/r c = w + cmplx(0.e0,1.e0)*orient*r*w11 r = r * sqrt(w11s) return end subroutine scplot(wb, wa, wd) c c trace boundary of circular arc polygon c integer wd(10) complex wb(10), wa(10) common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) common /scplt/ ifidif integer i c if(ifidif.eq.1) return do 10 i = 1, n write (6,1002) wb(i), wa(i), wd(i) 1002 format ('* ',1p4e15.7,i3) 10 continue return end subroutine scrss(wm,w1,wn,aa) c complex wm,w1,wn,aa c c this routine finds aa , the parameter in an intermediate c mobius transformation that enforces the circle to c circle condition. c common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c real a,b,c,d,dmb,cma,cl,dnom real x,y c a = real(wm/wn) b = aimag(wm/wn) c = real(w1/wn) d = aimag(w1/wn) c dnom = (a-1.e0)*(a+1.e0) + b*b if( abs(dnom).lt.radeps) go to 10 a = (a-1.e0)/dnom b = b/dnom dnom = (c-1.e0)*(c+1.e0) + d*d if( abs(dnom).lt.radeps) go to 20 c = (c-1.e0)/dnom d = d/dnom c dmb = d-b cma = c-a dnom = dmb*dmb+cma*cma x = (cma*(cma-2.e0*b*dmb)+dmb*dmb*(2.e0*a-1.e0))/dnom y = 2.e0*cma*(b*cma+(1.e0-a)*dmb)/dnom aa= cmplx(x,y) return 10 aa= cmplx(1.e0,0.e0) cl = (a-1.e0)/b dnom = (c-1.e0)*(c+1.e0) + d*d if( abs(dnom).lt.radeps) return c = (c-1.e0)/dnom d = d/dnom x = (cl*(cl-2.e0*d)+2.e0*c-1.e0)/(cl*cl+1.e0) y = 2.e0*cl*(d*cl-c+1.e0)/(cl*cl+1.e0) aa= cmplx(x,y) return 20 cl = (c-1.e0)/d x = (cl*(cl-2.e0*b)+2.e0*a-1.e0)/(cl*cl+1.e0) y = 2.e0*cl*(b*cl-a+1.e0)/(cl*cl+1.e0) aa= cmplx(x,y) return end subroutine selim(n,a,bet,gam) c integer n,i real a(1),bet(1),gam(1) c c eliminates algebraic sideconditions. c c input. c a(i) i=1,n -vertex coordinates in halfplane. c (on real axis.) c bet(i) i=1,n-3 -first n-3 beta-values. c gam(i) i=1,n -given constants. c output. c bet(n-2),bet(n-1) and bet(n). c real a1,a2,a3,b1,b2,b3 c b1=0.e0 b2=0.e0 b3=0.e0 bet(n-2)=0.e0 bet(n-1)=0.e0 bet(n) =0.e0 c c compute right hand side. c do 10 i=1,n b1=b1+bet(i) b2=b2+2.e0*bet(i)*a(i)+gam(i) b3=b3+a(i)*(bet(i)*a(i)+gam(i)) 10 continue a1=a(n-2) a2=a(n-1) a3=a(n) b1=b1-(bet(n)+bet(n-1)+bet(n-2)) b2=b2-2.e0*(bet(n)*a3+bet(n-1)*a2+bet(n-2)*a1) b3=b3-(bet(n)*a3*a3+bet(n-1)*a2*a2+bet(n-2)*a1*a1) c c start elimination. c b1=-b1 b2=-(b2+2.e0*b1*a1) b3=-(b3+.5e0*b2*(a1+a2)+b1*a1*a1) b3=b3/((a3-a1)*(a3-a2)) b2=.5e0*(b2-2.e0*b3*(a3-a1))/(a2-a1) b1=b1-b2-b3 bet(n-2)=b1 bet(n-1)=b2 bet(n) =b3 return end function sexprl (z) c complex version of Wayne's Aug 80 fnlib routine eric 3 Apr 81 c simplified error handling c undeclare abs, cexp -> exp 18 apr 85 c c evaluate (cexp(z)-1)/z . for small abs(z), we use the taylor c series. we could instead use the expression c sexprl(z) = (exp(x)*exp(i*y)-1)/z c = (x*exprel(x) * (1 - 2* sin(y/2)**2) - 2* sin(y/2)**2 c + i* sin(y)*(1+x*exprel(x))) / z c complex z, sexprl real rbnd, sqeps, drlz, r1mach real r, xn, xln, alneps data nterms, rbnd, sqeps / 0, 2*0.0e0 / c if (nterms.ne.0) go to 10 alneps = log(r1mach(3)) xn = 3.72 - 0.3*alneps xln = log((xn+1.0)/1.36) nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5 rbnd = r1mach(3) sqeps = 0.0 c 10 r = abs(z) if (r.le.0.5) go to 20 c sexprl = (exp(z) - 1.0) / z drlz= real(z) if( abs(drlz).lt.sqeps.and.abs(sexprl).lt.sqeps) goto 101 return c 20 sexprl = (1.0, 0.0) if (r.lt.rbnd) return c sexprl = (0.0, 0.0) do 30 i=1,nterms sexprl = 1.0 + sexprl*z/ float(nterms+2-i) 30 continue return c 101 write(6,1001) 1001 format('zexprl answer lt half prec, z near ni2pi') stop end subroutine sfneq(neqn, m, x, nf, f, nl2iv, nl2v, foo) c integer m, neqn, nl2iv(1) real x(m), f(neqn), nl2v(1) external foo c c this routine defines the system of nonlinear equations. c common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) common /method/ nrml, jeqn integer nrml, jeqn c common /sfail/ ifail integer ifail c integer i, k, is it=it+1 ifail=0 call sinteg(x) if(ifail.eq.0)goto 100 nf=0 return 100 continue if(jeqn.eq.1)then do 111 i = 1, n-1 k = 2*i f(k-1) = real(wt(i)-w(i)) f(k) = aimag(wt(i)-w(i)) 111 continue do 112 i=1,1 f(2*n-2+i) = arct(i) - arcc(i) 112 continue else if(jeqn.eq.2)then do 121 i = 1, n-1 k = 2*i f(k-1) = real(wt(i)-w(i)) f(k) = aimag(wt(i)-w(i)) 121 continue else if(jeqn.eq.3)then do 131 i = 1, n-1 k = 2*i f(k-1) = real(wt(i)-w(i)) f(k) = aimag(wt(i)-w(i)) 131 continue do 132 i=1,n f(2*n-2+i) = arct(i) - arcc(i) 132 continue else if(jeqn.eq.4)then do 141 i = 1, n-1 f(i) = abs(wt(i)-w(i)) 141 continue else if(jeqn.eq.5)then do 151 i = 1, n-1 f(i) = arct(i) - arcc(i) 151 continue end if call scplot(w,wc,wdir) return end subroutine sfode(r,y,yp) c real r,y(6),yp(6) c c routine called by ode. direct calculation. c for the schwartzian , third order equation. c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd,zk,ci,zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c common/node/ nfun integer nfun c integer i real ar, ai, br, bi real denom real z2r, z2i, csumr, csumi, c2r, c2i, zkr, zki real cr, cii, zdr, zdi, c1r, c1i c c zk is cut in circle-halfplane map. c zd is exp(%i*teta), teta direction of integration. c z(i) , i=1,2, n. contains the vertices in the c unit disk. (singular points.) c c-------------- c for debugging if(r.lt.0)then write(6,*)" zk=cmplx",zk write(6,*)" zd=cmplx",zd write(6,*)" z2=",z2 write(6,*)" ci=cmplx",ci write(6,*)" c2=",c2 do 8881 i=1,n write(6,*)" cc(",i,")=cmplx",cc(i) write(6,*)" gamma(",i,")=",gamma(i) write(6,*)" beta(",i,")=",beta(i) 8881 continue return end if c-------------- nfun = nfun + 1 zkr = real(zk) zki = aimag(zk) zdr = real(zd) zdi = aimag(zd) c cy2= cmplx(y(3),y(4)) c cy3= cmplx(y(5),y(6)) z2r=r*zdr z2i=r*zdi c csum = zero csumr = 0.e0 csumi = 0.e0 c c2 = 1.e0/(zk-z2) denom = (zkr-z2r)**2 + (zki-z2i)**2 c2r = (zkr-z2r) / denom c2i = - (zki-z2i) / denom c c = ci*c2*(zk+z2) ar = zkr+z2r ai = zki+z2i cii = c2r*ar-c2i*ai cr = -(c2i*ar+c2r*ai) c c2 = 8.e0*(c2*c2*zd*zk)**2 ar = c2r*c2r-c2i*c2i ai = 2*c2r*c2i br = ar*zdr-ai*zdi bi = ai*zdr+ar*zdi ar = br*zkr-bi*zki ai = bi*zkr+br*zki c2r = 8*(ar*ar-ai*ai) c2i = 8*(2*ar*ai) do 10 i = 1,n c c1 = .5e0/(c-cc(i)) ar = cr- real(cc(i)) denom = ar**2+cii**2 c1r = .5e0*ar/denom c1i = -.5e0*cii/denom c csum = csum + c1*(c1*gamma(i)+beta(i)) ar = c1r*gamma(i)+beta(i) ai = c1i*gamma(i) csumr = csumr + c1r*ar-c1i*ai csumi = csumi + c1i*ar+c1r*ai 10 continue c c1 = 1.5e0*cy3**2/cy2-c2*cy2*csum ar = 1.5e0*(y(5)*y(5)-y(6)*y(6)) ai = 1.5e0*2*y(5)*y(6) denom = y(3)**2+y(4)**2 c1r = (ar*y(3)+ai*y(4))/denom c1i = (ai*y(3)-ar*y(4))/denom ar = c2r*y(3)-c2i*y(4) ai = c2i*y(3)+c2r*y(4) c1r = c1r - (ar*csumr-ai*csumi) c1i = c1i - (ai*csumr+ar*csumi) yp(1)=y(3) yp(2)=y(4) yp(3)=y(5) yp(4)=y(6) yp(5)=c1r yp(6)=c1i c write(6,*)'r,yp6=',r,yp(6) return end subroutine sincp ( x, jprobl, jguess, inrml, ieqn ) c c initialize parameters describing target circular arc polygon c real x(1) integer jprobl, jguess common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) common /plotcm/ punp, ptrace, mtrace, ntrace, odet integer punp, ptrace, mtrace, ntrace complex odet(100,10) common /method/ nrml, jeqn integer nrml, jeqn integer i, j, k, m real cfwa, cfws, cfwt,t,t1, h real th(11) complex wtm, u common /wavegd/ cfws, cfwt c pi = 4.e0*atan(1.e0) eps = r1mach(4) iv(1) = 4 ntrace = 2 nrml = inrml jeqn = ieqn c ---- tolerances ---- radeps = 100*eps odeeps = 2000*eps epsfcn = 10*odeeps maxfun = 800 dmax = 1.e-2 hybeps = 10*odeeps c ci = cmplx(0.e0, 1.e0) zero = cmplx(0.e0, 0.e0) zk = cmplx( cos(0.33333e0), sin(0.33333e0)) c c ----------- particular regions ---------------- c c n number of sides c wtdir orientations (convex=1, concave=-1) c w0 image of 0 inside polygon c wt vertices c wtm point on side wt(n) -- wt(1) c wtc centers of sides c go to (10,20,30,40,50,60,70,30,30,30,710), jprobl c 10 n = 4 write(6,2010) 2010 format(' square example') do 15 i = 1, n wtdir(i) = 1 15 continue w0 = zero wt(1) = cmplx( .5e0, .5e0) wt(2) = cmplx(-.5e0, .5e0) wt(3) = cmplx(-.5e0, -.5e0) wt(4) = cmplx( .5e0, -.5e0) wtc(1) = cmplx(-1.e30, 0.) wtc(2) = cmplx( 0., -1.e30) wtc(3) = cmplx( 1.e30, 0.) wtc(4) = cmplx( 0., 1.e30) wtm = cmplx( .5e0, .0e0) c the following is needed to avoid "stiffness" in example 1,1,1,5 zk = cmplx( cos(0.11111e0), sin(0.11111e0)) go to 99 c 20 n = 5 write(6,2020) 2020 format(' notch example') do 25 i = 1, n wtdir(i) = 1 25 continue wtdir(1)=-1 w0 = zero wt(1) = cmplx( .0e0, .5e0) wt(2) = cmplx(-.5e0, .5e0) wt(3) = cmplx(-.5e0, -.5e0) wt(4) = cmplx( .5e0, -.5e0) wt(5) = cmplx( .5e0, .0e0) wtc(1) = cmplx( .5e0, .5e0) wtc(2) = cmplx( 0., -1.e30) wtc(3) = cmplx( 1.e30, 0.) wtc(4) = cmplx( 0., 1.e30) wtc(5) = cmplx(-1.e30, 0.) t = .25e0 * (2.e0-sqrt(2.e0)) wtm = cmplx(t,t) go to 99 c c 30 stop 40 stop 50 n = 4 write(6,2060) 2060 format(' sqeezed square example') wtdir(1) = -1 wtdir(2) = 1 wtdir(3) = -1 wtdir(4) = 1 w0 = zero wt(1) = cmplx( .5e0, .5e0) wt(2) = cmplx(-.5e0, .5e0) wt(3) = cmplx(-.5e0, -.5e0) wt(4) = cmplx( .5e0, -.5e0) wtc(1) = cmplx(4., 0.) wtc(2) = cmplx(0., -4.) wtc(3) = cmplx(-4., 0.) wtc(4) = cmplx(0., 4.) wtm = cmplx( real(wtc(1))-abs(wt(1)-wtc(1)), 0.e0) go to 99 60 n = 3 write(6,2050) 2050 format(' quarter circle example') wtdir(1) = 1 wtdir(2) = 1 wtdir(3) = 1 w0 = cmplx(.4e0,.4e0) wt(1) = cmplx( 0.e0, 1.e0) - w0 wt(2) = cmplx( 0.e0, 0.e0) - w0 wt(3) = cmplx( 1.e0, 0.e0) - w0 wtc(1) = cmplx(0., 0.) - w0 wtc(2) = cmplx(1.e30, 0.) - w0 wtc(3) = cmplx(0., 1.e30) - w0 wtm = cmplx( 1/sqrt(2.e0), 1/sqrt(2.e0) ) - w0 w0 = 0 go to 99 c 70 n = 5 cfws = 1 cfwt = 1 c jump here from 710 71 continue do 72 i = 1, n wtdir(i) = 1 72 continue write(6,2070) cfws, cfwt 2070 format(' quarter circular fin waveguide example' + /' 1/s, 1/t = ',2f15.2) cfwa = 1.e0 cfws = 1.e0/cfws cfwt = 1.e0/cfwt w0 = cmplx( cfwt/200+cfwa/2, cfws/4 ) wt(1)= cmplx( cfwt/2, cfwa ) - w0 wt(2)= cmplx( cfwt/2, cfws/2 ) - w0 wt(3)= cmplx( 0.e0, cfws/2 ) - w0 wt(4)= 0 - w0 wt(5)= cmplx( cfwt/2+cfwa, 0. ) - w0 wtc(1)=cmplx(cfwt/2,0.) - w0 wtm=cmplx( cfwt/2, 0.e0 ) + cfwa*cmplx(1/sqrt(2.),1/sqrt(2.)) - w0 w0 = 0 do 75 i = 1, 4 u = wt(i+1) - wt(i) u = u / abs(u) u = wt(i) - u*(real(u)*real(wt(i))+aimag(u)*aimag(wt(i))) wtc(i+1) = -1.e10 * u 75 continue go to 99 710 n = 5 go to 71 c 1 2 3 4 5 6 7 8 9 99 go to (100,110,120,130,140,150,160,170,200),jguess 100 continue m = 2*n-4 do 105 i=1,m x(i) = 0.e0 105 continue go to 200 110 stop x(1) = 2.3 x(2) = -1.3 x(3) = 0. x(4) = 1.3 go to 200 120 stop 130 stop 140 stop 150 continue x( 1) = 2.343616E+00 x( 2) = -1.331171E+00 x( 3) = -6.384445E-04 x( 4) = 1.332062E+00 x( 5) = 1.343304E+00 x( 6) = -6.599826E-01 go to 200 160 stop 170 continue m = 2*n-4 do 171 i=1,m x(i) = 0.e0 171 continue m = n-1 th(1) = 0 th(n+1) = 2*pi do 173 i=1,m u = (wt(i)-w0)/(wt(n)-w0) th(i+1) = atan2(aimag(u),real(u)) if(th(i+1).lt.0) th(i+1) = th(i+1) + 2*pi 173 continue do 175 i=1,m x(i) = log( (th(i+1)-th(i))/(th(i+2)-th(i+1)) ) 175 continue goto 200 c 200 call sincp2(x,wtm) write(6,1002) n 1002 format('* ',i3) call scplot(wt,wtc,wtdir) do 210 i = 1, n write(6,*) 'wt,wtc=', wt(i), wtc(i) 210 continue write(6,1001) radeps,odeeps 1001 format( x /' radeps = ',1pe11.3 x /' odeeps = ',1pe11.3 x ) return end subroutine sincp2 ( x, wtm ) c real x(1) complex wtm common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) common /method/ nrml, jeqn integer nrml, jeqn integer i real sarc, phij, phii complex bear c a = cmplx(1.e0,0.e0) b = zero c = zero d = a winit1 = a winit2 = 2.e0*(1.e0/zk-ci) c c compute normalizing constant for true polygon c (this is the last time wtm is used) if(nrml.eq.1.or.nrml.eq.3)then call scrss(wtm,wt(1),wt(n),cnr) else if(nrml.eq.2)then else stop end if do 350 i = 1, n j = mod(i,n)+1 bear = wtc(i)-wt(i) bear = bear / abs(bear) phii = atan2(aimag(bear),real(bear)) if(wtdir(i).ne.1) phii = phii + pi bear = wtc(j)-wt(i) bear = bear / abs(bear) phij = atan2(aimag(bear),real(bear)) if(wtdir(j).ne.1) phij = phij + pi alpha(i) = mod(5*pi+phii-phij,2*pi) gamma(i) = 1.e0-(alpha(i)/pi)**2 350 continue do 360 i = 1, n j = mod(i-2+n,n)+1 radt(i) = abs(wt(i)-wtc(i)) arct(i) = sarc(wt(j),wt(i),wtc(i),wtdir(i)) 360 continue return end subroutine sinteg(x) c real x(1) c c this routine is the main driver defining the nonlinear c system of equations. c common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 c common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c common /capcm3/ wm complex wm(10) common/node/ nfun integer nfun common /plotcm/ punp, ptrace, mtrace, ntrace, odet integer punp, ptrace, mtrace, ntrace complex odet(100,10) common /method/ nrml, jeqn integer nrml, jeqn c common /sfail/ ifail integer ifail c c integer i, is, j real marg(10), zarg(10) real sarc complex mc, wd(2, 10) complex aa, mapwn c c ...for nrml=3 optimization of normalization by n2f2... integer nl2iv(300) real nl2v(300) external werr integer rdreq, lmax0 complex copt, dopt parameter ( rdreq=57, lmax0=35 ) c c c write(6,1010) it 1010 format(/' -------------------------- ',i5) nfun = 0 call sstpar(x) do 10 i = 1, n zarg(i)=atan2(aimag(z(i)), real(z(i))) if(zarg(i).lt.0) zarg(i) = zarg(i) + 2*pi 10 continue marg(1) = zarg(1)/2.e0 is = n-1 do 20 i = 2, is marg(i) = zarg(i-1)+(zarg(i)-zarg(i-1))/2.e0 20 continue marg(n) = zarg(n-1)+(2.e0*pi-zarg(n-1))/2.e0 winit2 = (a/d**2)*(d*winit2-2.e0*c*winit1**2) winit1 = (a/d)*winit1 c write(6,1019) winit1, winit2 1019 format(' (sinteg) winit1 = ',2e15.7 + /' winit2 = ',2e15.7 ) mtrace = n do 30 i = 1, n call sray(marg(i),wm(i),wd(1,i),wd(2,i), + ntrace,odet(1,i)) if(ifail.ne.0) return 30 continue c write(6,2) nfun 2 format(1x,'ode used ',i7,2x,'function evaluations') call smtov(wm, wd) if(ifail.ne.0) return c c **** normalization **** c.... if(nrml.eq.1.or.nrml.eq.3)then c a0 = 1./w(n) c a1 = aa c a2 = cnr c a3 = 1./wt(n) call scrss(wm(1),w(1),w(n),aa) if( (abs(aa-1).le.radeps) .and. (abs(cnr-1).le.radeps) ) then a = 1./w(n) b = zero c = zero d = 1./wt(n) else if( (abs(aa-1).le.radeps) )then a = (1.-aa)/w(n) b = zero c = 1./(w(n)*wt(n)) d = aa/wt(n) else if( (abs(cnr-1).le.radeps) )then a = aa/w(n) b = zero c = 1./(w(n)*wt(n)) d = (aa-1.)/wt(n) else a = aa*(1.-cnr)/w(n) b = zero c = (aa-cnr)/(w(n)*wt(n)) d = (cnr/wt(n))*(1.-aa) end if c.... else if(nrml.eq.2)then a = 1. b = zero aa = wt(1)*wt(n)*(w(n)-w(1)) c = ( wt(1)*w(n) - w(1)*wt(n) ) / aa d = w(1)*w(n)*(wt(n)-wt(1)) / aa c.... else write(6,1038) nrml 1038 format(' normalization ',i5,' is not implemented') stop end if c.... if(nrml.eq.3)then copt = 0 call ivset(1,nl2iv,300,300,nl2v) nl2iv(rdreq) = 0 nl2v(lmax0) = .2 call n2f2(2*(n-1),2,copt,werr, $ nl2iv,300,300,nl2v,nl2iv,nl2v,werr) dopt = 1 - copt*w(n) c = copt*a+dopt*c d = dopt*d end if c.... mapwn = a*w(n)/(c*w(n)+d) c write(6,1040) wm(1),w(1),w(n),aa,cnr,mapwn 1040 format(' wm =',2e13.4 + /' w1 =',2e13.4 + /' wn =',2e13.4 + /' aa =',2e13.4 + /' cnr =',2e13.4 + /' mapwn =',2e13.4) c write(6,1039) a, c, d 1039 format(' a =',2e13.4 + /' c =',2e13.4 + /' d =',2e13.4) c c apply the mobius transformation c c write(6,1041) (w(i),wc(i),i=1,n) c do 44 i = 1, mtrace c do 43 j = 1, ntrace c odet(j,i) = a*odet(j,i)/(c*odet(j,i)+d) c 43 continue c 44 continue do 45 i = 1, n w(i) = a*w(i)/(c*w(i)+d) mc = a*wc(i)/(c*wc(i)+d) call smapc(wc(i), abs(wm(i)-wc(i)), a, b, c, d, wc(i)) wm(i) = a*wm(i)/(c*wm(i)+d) if (abs(wm(i)-wc(i)) .lt. abs(mc-wc(i))) wdir(i) = -wdir(i) 45 continue c write(6,1041) (w(i),wc(i),i=1,n) 1041 format(' w, wc =',10(/1x,4e13.4)) c c restore initialization values to naive guess, c for debugging c c a = cmplx(1.e0,0.e0) c b = zero c c = zero c d = a c winit1 = a c winit2 = 2.e0*(1.e0/zk-ci) c c compute radius of final circular arcs. c do 50 i = 1, n radc(i) = abs(wc(i)-w(i)) 50 continue c c compute arc-length of final circular arcs. c do 60 i = 1, n is = mod(i-2+n,n)+1 arcc(i) = sarc(w(is),w(i),wc(i),wdir(i)) 60 continue return end c*************************************************************** subroutine werr ( nn, mm, copt, nf, f, nl2iv, nl2v, foo ) common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c complex copt, dopt, ws, lc, ld integer i, k, nn, mm, nf, nl2iv(1) real f(nn), nl2v(1) external foo c c starting normalization is u = a x / ( c x + d ). c follow this by something of the form v = u / ( copt u + dopt ) c (We optimize copt rather than c directly in a hope to deal c with a better conditioned problem.) c dopt = 1 - copt*w(n) lc = copt*a+dopt*c ld = dopt*d do 190 i=1,n-1 k = 2*i ws = a*w(i)/(lc*w(i)+ld) f(k-1) = real(ws-wt(i)) f(k) = aimag(ws-wt(i)) 190 continue end subroutine sinter ( m, n, s, ijac, idj, jac, f, p ) c c this routine defines the equation for the intersection c of two circles, based on an expansion around the c midpoint data. c real s(1), jac(idj,1), f(1), p(1) integer n, ijac, idj complex wm, wn, ws1, ws2, ci, e2 complex sexprl, e real r1, r2, sr ci= cmplx(0.e0,1.e0) if(ijac.eq.1) goto 200 c wm= cmplx(p(1),p(2)) wn= cmplx(p(3),p(4)) r1=p(5) e=ci*s(1)/r1 e=sexprl(e) e=s(1)*ci*e ws1=wm+wn*e wm= cmplx(p(6),p(7)) wn= cmplx(p(8),p(9)) r2=p(10) e=ci*s(2)/r2 e=sexprl(e) e=s(2)*ci*e ws2=wm+wn*e f(1)= real(ws1)- real(ws2) f(2)=aimag(ws1)-aimag(ws2) goto 900 c 200 continue sr=s(1)/p(5) e2= cmplx(p(3),p(4))* cmplx( cos(sr), sin(sr)) jac(1,1)=-aimag(e2) jac(2,1)= real(e2) sr=s(2)/p(10) e2= cmplx(p(8),p(9))* cmplx( cos(sr), sin(sr)) jac(1,2)=aimag(e2) jac(2,2)=- real(e2) 900 continue return end subroutine smapc(c0, r, a, b, c, d, c3) c real r complex c0, a, b, c, d, c3 c c given the center c0 and radius r of a circle c and the coefficients of a mobius transformation (a*z+b)/(c*z+d), c compute the center c3 of the image of the circle. c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c real e, srinv complex c1, c2, b1, b2 c if (abs(c) .ge. eps) goto 10 c3 = (a*c0+b)/d goto 40 10 c1 = c0+d/c e = abs(c1) if (e .ge. eps) goto 20 c3 = a/c goto 40 20 if ( abs(e-r) .ge. 10*eps) goto 30 c3 = cmplx(srinv(0.0e0), 0.0e0) goto 40 30 b1 = (1.0e0-r/e)*c1 b2 = (r/e+1.0e0)*c1 if(b1*b2*c.eq.0)write(6,*)"smapc",r,e,c1,b1,b2,c c2 = .5e0*(1.e0/b1+1.e0/b2) c3 = a/c+((b*c-a*d)/(c*c))*c2 40 continue return end subroutine smtov(wm, wd) c complex wm(10), wd(2, 10) c c solve for vertices from midpoint information c common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c common /sfail/ ifail integer ifail c integer i, j, k, is, maxstp, iprint real ws(2) real tol, tr1(2), tr2(2,2), tr3(12), parm(10) real u, v, marg(10), zarg(10) real t, b1, b2, e real theta, sarc, oarc real det, a11, a12, a21, a22 complex w1, w2, ni, nj, zm, em, sexprl external sinter c do 10 i = 1, n zarg(i)=atan2(aimag(z(i)), real(z(i))) if(zarg(i).lt.0) zarg(i) = zarg(i) + 2*pi 10 continue marg(1) = zarg(1)/2.e0 is = n-1 do 20 i = 2, is marg(i) = zarg(i-1)+(zarg(i)-zarg(i-1))/2.e0 20 continue marg(n) = zarg(n-1)+(2.e0*pi-zarg(n-1))/2.e0 do 30 i = 1, n zm = cmplx( cos(marg(i)), sin(marg(i))) call scirc(zm, wm(i), wd(1, i), wd(2, i), wc(i), radc(i), + wdir(i)) c write(6,1025) i,wm(i),i,wd(1,i),i,wd(2,i),i,wc(i), c + i,radc(i),i,wdir(i) 1025 format( + /' wm( ',i2,') = ',1p2e13.4 + /' wd(1,',i2,') = ',1p2e13.4 + /' wd(2,',i2,') = ',1p2e13.4 + /' wc( ',i2,') = ',1p2e13.4 + /' radc(',i2,') = ',1p1e13.4 + /' wdir(',i2,') = ',i3 ) 30 continue c c find intersection points c do 90 k = 1, n i = k j = mod(k, n)+1 if(abs(wc(i)-wc(j)).le.(radc(i)+radc(j))) goto 39 write(6,1039) 1039 format(' sides do not intersect!') ifail=1 return 39 if (radc(j) .le. radc(i)) goto 40 i = j j = k 40 if (radeps*radc(j) .le. 1.) goto 50 c c two lines c b1 = ( real(wm(i)))* real(wd(1,i)) 1 +(aimag(wm(i)))*aimag(wd(1,i)) b2 = ( real(wm(j)))* real(wd(1,j)) 1 +(aimag(wm(j)))*aimag(wd(1,j)) a11 = real(wd(1, i)) a12 = aimag(wd(1, i)) a21 = real(wd(1, j)) a22 = aimag(wd(1, j)) det = (a11)*a22-(a12)*a21 w1=cmplx((b1*a22-a12*b2)/det,(-(a11*b1-b2*a21))/det) w2=cmplx(1.e+15, 0.e0) goto 70 50 if (radeps*radc(i) .le. 1.) goto 60 c c circle and line c b1 = real(wd(1,i)) b2 = aimag(wd(1,i)) e = sqrt(b1**2+b2**2) u = ( real(wd(1, i))*( real(wm(i)-wc(j))) 1 +aimag(wd(1, i))*(aimag(wm(i)-wc(j))))/e v = sqrt((radc(j)-u)*(radc(j)+u)) w1 = wc(j)+ cmplx(u, v)*wd(1, i)/e w2 = wc(j)+ cmplx(u, -v)*wd(1, i)/e goto 70 c c two circles c 60 continue if(radeps*radc(i).gt.radc(j) .or. 1 radeps*radc(j).gt.radc(i) ) then write(6,1802) 1802 format(" circles with wildly differing radii") ifail=3 return endif b1 = ( real(wc(j)))- real(wc(i)) b2 = (aimag(wc(j)))-aimag(wc(i)) e = sqrt(b1**2+b2**2) b1 = radc(j) b2 = radc(i) u = (e**2-b2**2+b1**2)/(2.e0*e) v = sqrt((b1-u)*(b1+u)) w1 = wc(j)+ cmplx(u, v)*(wc(i)-wc(j))/e w2 = wc(j)+ cmplx(u, -v)*(wc(i)-wc(j))/e c c choose intersection point with correct interior angle c 70 i = k j = mod(k, n)+1 ni = w1-wc(i) nj = w1-wc(j) if(abs(ni).eq.0.0.or.abs(nj).eq.0.0) goto 801 ni = cmplx( float(wdir(i)), 0.e0)*ni/abs(ni) nj = cmplx( float(wdir(j)), 0.e0)*nj/abs(nj) theta=atan2(aimag(nj),real(nj))-atan2(aimag(ni),real(ni)) if ( abs(mod(alpha(k)-theta+9.e0*pi, 2.e0*pi)-pi) 1 .ge. 1.e0) go to 80 w(k) = w1 goto 85 80 w(k) = w2 85 continue zm= cmplx( cos(marg(i)), sin(marg(i))) ni=wdir(i)*wd(1,i)*zm ni=ni/abs(ni) zm= cmplx( cos(marg(j)), sin(marg(j))) nj=wdir(j)*wd(1,j)*zm nj=nj/abs(nj) c write(6,1031)i,ni,j,nj,w(i) 1031 format( + /' snewt2' + /' i=',i2 + ,' ni=',0p2f7.3 + /' j=',i2 + ,' nj=',0p2f7.3 + /' w before =',0p2e12.4 ) parm(1)= real(wm(i)) parm(2)=aimag(wm(i)) parm(3)= real(ni) parm(4)=aimag(ni) parm(5)=radc(i) parm(6)= real(wm(j)) parm(7)=aimag(wm(j)) parm(8)= real(nj) parm(9)=aimag(nj) parm(10)=radc(j) tol=5*eps*abs(w(k)) maxstp=10 iprint=0 c**** 20 Apr 85 added wdir c write(6,1803) w(k),wm(i),wc(i),wm(j),wc(j) c1803 format('for',1p10e10.2) c ws(1)=sarc(w(k),wm(i),wc(i),1) c ws(2)=sarc(w(k),wm(j),wc(j),1) c write(6,1804) ws(1), ws(2) c1804 format('new',1p10e10.2) ws(1)=oarc(w(k),wm(i),wc(i)) ws(2)=oarc(w(k),wm(j),wc(j)) c write(6,1805) ws(1), ws(2) c1805 format('old',1p10e10.2) call snewt2(2,2,ws,tr1,tol,maxstp,1.0e0,iprint, + 2,tr2,sinter,parm,tr3) em=ci*ws(1)*sexprl(ci*ws(1)/radc(i)) w(k)=wm(i)+ni*em c write(6,1035)w(k) 1035 format(' w after =',2e12.4 ) 90 continue return 801 continue write(6,1801) 1801 format(' smtov: entire side mapped onto a point') ifail=2 return end subroutine snewt2(m,n,x,f,tol,maxstp,dmax, + iprint,idj,jac,fun,parm,w) c integer m,n,idj,maxstp,iprint real x(1),f(1),jac(idj,1) real dmax real w(1) real parm(1) external fun c c This is a simple Newton code for solving a system of nonlinear c equations. If the number of equations m exceeds the number of c unknowns n then a nonlinear least squares solution is attempted. c c dmax is the maximum 2-norm of any step that will be attempted. c c The user must supply a subroutine fun, declared external in the c calling program , with the following callingsequence: c c fun(m,n,x,ijac,idj,jac,f,parm) c c if ijac=0 : f should return the function value at x. c if ijac=1 : jac should return the jacobian Df(i)/Dx(j). c In this case f is input and already contains c the function value at x. (from a previous call) c c This code will return when the 2-norm of f is less than tol or c when maxstp steps have been tried. (Each involving one jacobian c evaluation and at most 3 function evaluations.) c c The code may also give up and in this case some information c about the function is written out. c c w is a workspace containing at least 5*n+m elements. c parm is a real array in which the user can pass information c to fun. It is never referenced by Newton. c c Local variables. c integer i1,i2,i3,i4,i5 real fnrm c if(maxstp.lt.1) go to 100 if(m.lt.n) go to 200 c c assign workspace. c i1 = n+1 i2 = i1+m i3 = i2+n i4 = i3+n i5 = i4+n c nfun = 0 njac = 0 c call snwt2(m,n,x,f,tol,maxstp,dmax,iprint, + idj,jac,fun,parm,nfun,njac,fnrm, + w(1),w(i1),w(i2),w(i3),w(i4),w(i5)) c if(iprint.le.0) return write(6,1) fnrm write(6,2) nfun write(6,3) njac return 100 write(6,4) return 200 write(6,5) return 1 format(1x,'on exit from newton the norm of f is,'e12.3) 2 format(1x,'number of function evaluations were ,'i6) 3 format(1x,'number of jacobian evaluations were ,'i6) 4 format(1x,'newton returned because maxstp was ,'i6) 5 format(1x,'newton returned because m is less than n, m,n=',2i4) end subroutine snwt2(m,n,x,f,tol,maxstp,dmax,iprint, + idj,jac,fun,parm,nfun,njac,fnrm, + xnew,fnew,stepn,qraux,jpvt,stepg) c integer m,n,maxstp,iprint,idj,nfun,njac real tol,dmax,fnrm integer jpvt(1) real x(1),f(1) real jac(idj,1) real xnew(1),fnew(1),stepn(1),qraux(1),stepg(1) real parm(1) c c this routine must always be called from the c driver routine newtn. c c local variables. c integer i,j,it,info real fnrm1,fnrm2,fnrm3,fmin real dec1,dec2 real a,b,c,t real alfa real sdot,snrm2 real tmp(10,10),xtmp(10),ytmp(10),ztmp(10) c c dec1 is used to test for decrease. c dec2 is used to test for decrease if dec1 cannot be satisfied. c alfa is steplenght parameter. c dec1 = .2e0 dec2 = .9e0 call scopy(n,x,1,xnew,1) nfun = nfun + 1 call fun(m,n,x,0,idj,jac,f,parm) fnrm = snrm2(m,f,1) if(fnrm.lt.tol) go to 600 fmin = fnrm c c start main iteration loop. c do 400 it=1,maxstp if(iprint.le.0) go to 210 write(6,10) it,fnrm write(6,11)(x(i),i=1,n) write(6,12)(f(i),i=1,m) 210 njac = njac + 1 call fun(m,n,x,1,idj,jac,f,parm) c c if iprint is 3 or larger then diagnostic info c about the jacobian can be (computed and) printed c here. c if(iprint.le.2) go to 220 write(6,18) write(6,19)((jac(i,j),j=1,n),i=1,m) do 215 i=1,n do 215 j=1,n tmp(i,j)=jac(i,j) 215 continue do 216 i=1,n ytmp(i)=f(i) 216 continue call sgeco(tmp,10,n,xtmp,rcond,ztmp) call sgesl(tmp,10,n,xtmp,ytmp,0) write(6,719) rcond write(6,718)(ztmp(i),i=1,n) write(6,718)(ytmp(i),i=1,n) 719 format(1x,'rcond=',e12.3) 718 format(1x,4e12.3) 220 do 230 i = 1,n jpvt(i) = 0 230 continue call sqrdc(jac,idj,m,n,qraux,jpvt,stepg,1) call sqrsl(jac,idj,m,n,qraux,f,f,f,stepg,f,f,100,info) if(info.ne.0) go to 520 do 240 i = 1,n stepn(jpvt(i)) = -stepg(i) 240 continue if(iprint.le.2) go to 250 write(6,20)(stepn(i),i=1,n) c c test if newton step is too large. c 250 fnrm1 = snrm2(n,stepn,1) if(fnrm1.lt.dmax) go to 258 if(iprint.ge.3) write(6,22) fnrm1,dmax c c the gradient of the squared norm of f is computed and c normalized to lenght dmax. The step is taken to be c between the newton step and this c gradient direction. (See below for details.) c do 256 j = 1,n t = 0.e0 do 255 i = 1,m t = t - f(i)*jac(i,j) 255 continue stepg(j) = t 256 continue c t = dmax/fnrm1 call sscal(n,1.e0/fnrm1,stepn,1) fnrm1 = snrm2(n,stepg,1) call sscal(n,-1.e0/fnrm1,stepg,1) fnrm1 = sdot(n,stepn,1,stepg,1) if(iprint.ge.2) write(6,17) fnrm1 c do 257 i = 1,n c stepn(i) = t*stepn(i) + (1.e0-t)*stepg(i) c 257 continue call sscal(n,dmax,stepn,1) c c try a unit step. c 258 call saxpy(n,1.e0,stepn,1,xnew,1) if(iprint.ge.3) write(6,21)(xnew(i),i=1,n) nfun = nfun + 1 call fun(m,n,xnew,0,idj,jac,fnew,parm) fnrm1 = snrm2(m,fnew,1) if(iprint.le.1) go to 260 write(6,13) fnrm1 write(6,11) (xnew(i),i=1,n) write(6,12) (fnew(i),i=1,m) 260 fmin = min(fmin,fnrm1) if(fnrm1.lt.dec1*fnrm) go to 350 if(fnrm1.gt.fmin) go to 270 alfa = 1.e0 call scopy(m,fnew,1,f,1) c c try half a unit step. c 270 call saxpy(n,-.5e0,stepn,1,xnew,1) nfun = nfun + 1 call fun(m,n,xnew,0,idj,jac,fnew,parm) fnrm2 = snrm2(m,fnew,1) if(iprint.le.1) go to 280 write(6,14) fnrm2 write(6,11)(xnew(i),i=1,n) write(6,12)(fnew(i),i=1,m) 280 fmin = min(fmin,fnrm2) if(fnrm2.lt.dec1*fnrm) go to 350 if(fnrm2.gt.fmin) go to 290 alfa = .5e0 call scopy(m,fnew,1,f,1) c c compute the quadratic passing through fnrm, fnrm2 and fnrm1. c 290 a = 2.e0*(fnrm+fnrm1-2.e0*fnrm2) b = 4.e0*fnrm2-3.e0*fnrm-fnrm1 c = fnrm if(a.le.0.e0) go to 310 t = -b/(2.e0*a)-.5e0 c c try a step to the minimum of a quadratic model. c call saxpy(n,t,stepn,1,xnew,1) t = t + .5e0 nfun = nfun + 1 call fun(m,n,xnew,0,idj,jac,fnew,parm) fnrm3 = snrm2(m,fnew,1) if(iprint.le.1) go to 300 write(6,15) a,b,c write(6,16) fnrm3,t write(6,11)(xnew(i),i=1,n) write(6,12)(fnew(i),i=1,m) 300 fmin = min(fmin,fnrm3) if(fnrm3.lt..5e0*(dec1+dec2)*fnrm) go to 350 if(fnrm3.gt.fmin) go to 310 alfa = t call scopy(m,fnew,1,f,1) 310 if(fmin.gt.dec2*fnrm) go to 530 call scopy(n,x,1,xnew,1) call saxpy(n,alfa,stepn,1,xnew,1) go to 360 c c update and prepare for a new iteration. c 350 call scopy(m,fnew,1,f,1) 360 call scopy(n,xnew,1,x,1) fnrm = fmin if(fnrm.lt.tol) go to 600 400 continue write(6,1) return 520 write(6,3) return 530 write(6,4) fnrm write(6,11)(x(i),i=1,n) write(6,7)(stepn(i),i=1,n) write(6,5) fnrm,fnrm2,fnrm1 write(6,6) alfa,fmin write(6,11)(xnew(i),i=1,n) write(6,12)(f(i),i=1,m) return 600 continue return 1 format(1x,'maxstp iterations in newtn') 3 format(1x,'jacobian has rank less than n') 4 format(1x,'Newton gave up, norm of f is:'e12.4) 5 format(1x,'Norm of f in x+alfa*stepn , alfa=0,1/2,1, is :'3e12.3) 6 format(1x,'With step= ',e12.4,1x,',the norm is :',e12.3) 7 format(1x,'s= ',1x,5e14.6,4(/5x,5e14.6)) 10 format(/1x,'At the start of iteration',i3,4x,'Norm of f= ',e14.6) 11 format(1x,'x= ',1x,5e14.6,4(/5x,5e14.6)) 12 format(1x,'f= ',1x,5e14.6,4(/5x,5e14.6)) 13 format(1x,'Unit step,',1x,'Norm of f=',e12.4) 14 format(1x,'Half step,',1x,'Norm of f=',e12.4) 15 format(1x,'Quadratic model coeff a,b,c=',1x,3e12.4) 16 format(1x,'Min step,',1x,'Norm of f=',e12.4,2x,'alfa=',e12.4) 17 format(1x,'angle between newton and gradient direction,',e12.3) 18 format(1x,'jacobian matrix.') 19 format(10(5x,5e12.3/)) 20 format(1x,'newton step,'/,5(5x,5e12.3/)) 21 format(1x,'newton step gives x=,'/,5(5x,5e12.3/)) 22 format(1x,'step too large, norm of step,dmax= ',2e12.3) end subroutine sray(theta,y0,y1,y2,ntrace,trace) c c trace radial ray in unit disk to unnormalized polygon c integer ntrace real theta complex y0, y1, y2, trace(ntrace) c common /initc/ a, b, c, d, winit1, winit2 complex a, b, c, d, winit1, winit2 c common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun integer maxfun real radeps,odeeps,hybeps,epsfcn,dmax c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c common /savex/ xxx real xxx(10) common/node/ nfun integer nfun c external sfode, sfode4 integer iflag, iwork(50), odeit, odemit, step real y(6), relerr, abserr, r1, r2, rwork(600) real yeps(6), ymag(6) common /sfail/ ifail integer ifail integer istate, jt, index, order, limits(7) real dt, dtmax data odemit /10/ c zd = cmplx( cos(theta), sin(theta)) y(1) = 0.e0 y(2) = 0.e0 y(3) = real(winit1*zd) y(4) = aimag(winit1*zd) y(5) = real(winit2*zd**2) y(6) = aimag(winit2*zd**2) r1 = 0.e0 relerr = odeeps abserr = 1.e-5*relerr iflag = 1 odeit = 1 trace(1) = cmplx(y(1),y(2)) do 50 step=2,ntrace r2 = (step-1.0e0)/(ntrace-1.0e0) c c istate=1 c jt=2 c call lsoda(sfode4,6,y,r1,r2,1,relerr,abserr,1,istate, c $ 0,rwork,200,iwork,26,sfode4,jt) c if(istate.ne.2) write(6,*)"lsoda returned",istate c... c index=0 c dt=.001 c dtmax=.1 c do 10 i=1,6 c yeps(i) = relerr c 10 continue c 11 call stride(6,sfode,sfode,r1,r2,y,ymag,index,dt,dtmax, c $ order,yeps,101,limits,iwork,rwork) c if(index.ge.1) goto 11 cc... c goto 40 21 call ode(sfode,6,y,r1,r2,relerr,abserr,iflag,rwork,iwork) if(iflag.eq.2) goto 40 if(iflag.ne.3) goto 22 write(6,1021) relerr,abserr 1021 format(" ode increased relerr,abserr to ",1p2e11.2) odeit=odeit+1 if(odeit.le.odemit)goto 21 ifail=1 return c 22 if(iflag.ne.4) goto 23 write(6,1022) 1022 format(" ode complained that it had used more than 500 steps") odeit=odeit+1 if(odeit.le.odemit)goto 21 ifail=1 return c 23 if(iflag.ne.5) goto 24 write(6,1023) r1, r2 1023 format(" ode complained equations are stiff r=",2f8.4) c------ m = 2*n-4 write(6,8024) (i,xxx(i),i=1,m) 8024 format(' x(',i2,')=',f22.16) y(1) = 0.e0 y(2) = 0.e0 y(3) = real(winit1*zd) y(4) = aimag(winit1*zd) y(5) = real(winit2*zd**2) y(6) = aimag(winit2*zd**2) do 8023 i=1,6 write(6,*)" y(",i,")=",y(i) 8023 continue call sfode(-1,y,yp) stop c------ odeit=odeit+1 if(odeit.le.odemit)goto 21 ifail=1 return c 24 write(6,1) iflag 1 format(1x,'ode returned iflag = ',1x,i8/) ifail=1 return c 40 continue trace(step) = cmplx(y(1),y(2)) 50 continue c y0 = cmplx(y(1), y(2)) y1 = cmplx(y(3), y(4))/zd y2 = cmplx(y(5), y(6))/zd**2 return end function srinv(x) c real x c if ( abs(x) .ge. 1.e-15) goto 1 srinv = sign(1.e+15, x) goto 2 1 srinv = 1.0e0/x 2 continue return end subroutine sstpar(x) c real x(1) c c given guesses of vertex preimages and n-3 of the c auxiliary parameters, solve for remaining 3 parameters c common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero complex wt(10), cnr, wtc(10), z(10), cc(10), w(10) complex zd, zk, ci, zero, wc(10), w0 c common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir, 1 wdir, iv, radt, radc, arct, arcc integer n, it, wtdir(10), wdir(10), iv(1) real s(6) real alpha(10), gamma(10), beta(10), pi, eps real radt(10), radc(10), arct(10), arcc(10) c common /savex/ xxx real xxx(10) integer i, m real a(10) call strans(x, z, n) m = 2*n-4 do 20 i = 1, n cc(i)= ci*(zk+z(i))/(zk-z(i)) a(i) = real(cc(i)) xxx(i) = x(i) 20 continue do 30 i = n, m beta(i-n+1) = x(i) 30 continue call selim(n, a, beta, gamma) return end subroutine strans(x,z,n) integer n real x(1), a(25) complex z(n) if(n.gt.25)then write(6,*)'strans a dimension too samll' stop end if call stransa(x,z,a,n) return end subroutine stransa(x,z,a,n) c c given optimization parameters x, back transform to c get angles a and points z on unit circle c c x(i)=log((teta(i)-teta(i-1))/(teta(i+1)-teta(i))) c for i=1,2.. n-1. c teta(0)=0,teta(n)=twopi. c z(i)=exp(%i*teta(i)). c integer n real x(1), a(n) complex z(n) integer i,nm real y,sum,twopi twopi=8.e0*atan(1.e0) y=1.e0 sum=y nm=n-1 do 10 i=1,nm y=y/exp(x(i)) sum=sum + y 10 continue y=twopi/sum sum=y a(1)=sum z(1)= cmplx( cos(y), sin(y)) do 20 i=2,nm y=y/exp(x(i-1)) sum=sum+y a(i)=sum z(i)= cmplx( cos(sum), sin(sum)) 20 continue a(n)=0.e0 z(n)= cmplx(1.e0,0.e0) return end c From jcb Thu Feb 6 15:16 EST 1986 SUBROUTINE STRIDE(N,F,JAC,X,XPRINT,Y,YMAG,INDEX,STP,STPMAX,M,EPS, + METHOD,LIMITS,INTGRS,REALS) COMMENT OOO: C C S T R I D E C C STABLE RUNGE-KUTTA INTEGRATOR FOR DIFFERENTIAL EQUATIONS C C J C BUTCHER, K BURRAGE, F H CHIPMAN C C MARCH 1979 C C C THIS VERSION OF STRIDE CONTAINS THE CHANGES THAT WERE SPECIFIED C IN F.CHIPMANS LETTER TO P.W.GAFFNEY ON SEPTEMBER 21,1979. C COMMENT 020: DECLARATIONS C DIMENSION Y(1),YMAG(1),EPS(1),LIMITS(1),INTGRS(1),REALS(1) INTEGER YTEMP,YSAVE,FTOTAL,FVECT,YVECT,FSPACE,YSPACE,EFACTS, +T,TINV,XI,DERIV,YRES,W,YLAST,FACTS LOGICAL CONVGD,DIVGD,EXCESS,REVAL,OK C COMMENT 040: IF NOT FIRST CALL UNPACK REALS AND INTGRS C IF(INDEX.EQ.0)GO TO 120 GO TO 080 C COMMENT 060: STATEMENTS REPRESENTING THE 'PAUSE' ROUTINE C 060 IF(INDEX .GT. 1) LIMITS(INDEX)=LIMITS(INDEX)-1 IF (INDEX .GT. 2) LIMITS(1)=LIMITS(1)-1 IF (INDEX.LT.3) GO TO 070 IF (INTGRS(LIMMEM+INDEX).GT.0 .AND. + LIMITS(INDEX)+INTGRS(LIMMEM+INDEX).LT.0) GOTO 070 GOTO 100 C THEN 070 IF(INDEX.GE.3)LIMITS(INDEX)=INTGRS(LIMMEM+INDEX) IF(INDEX.EQ.1)RETURN INTGRS(1) = MLAST INTGRS(3) = MSTATE INTGRS(4) = MTRAN INTGRS(5) = NJAC INTGRS(6) = NSTEP INTGRS(7) = 0 IF(REVAL)INTGRS(7) = 1 REALS(2) = EGROW REALS(3) = ERREST REALS(4) = H REALS(5) = HFACTS REALS(6) = HLAST REALS(7) = XRES REALS(8) = XLAST REALS(9) = XJAC RETURN 080 IF(INDEX.GT.1)GO TO 085 MAXM = M IF(MAXM.GT.15)MAXM=15 INTGRS(2) = MAXM 085 URND = REALS(1) MAXM = INTGRS(2) MAXI=7 LIMMEM = MAXI + MAXM + 1 IPERM = LIMMEM + 7 EFACTS = 9 FVECT = EFACTS + INTGRS(MAXI+MAXM+1) YVECT = FVECT + MAXM YSPACE = YVECT + MAXM - N FSPACE = YSPACE + (MAXM + 2)*N T = FSPACE + (MAXM + 2)*N + N - MAXM TINV = T + MAXM**2 XI = TINV + MAXM**2 + MAXM DERIV = XI + MAXM*(MAXM + 1)/2 YRES = DERIV + N YTEMP = YRES + N YSAVE = YTEMP + N W = YSAVE + N FTOTAL = W + N YLAST = FTOTAL + N FACTS = YLAST + N - N MERPST=METHOD / 100 MTYPE=METHOD / 10 MAXITS=METHOD-10*MTYPE MTYPE=MTYPE-10*MERPST IF(INDEX.NE.1)GO TO 090 GO TO 095 C THEN 090 MLAST = INTGRS(1) MSTATE = INTGRS(3) MTRAN = INTGRS(4) NJAC = INTGRS(5) NSTEP = INTGRS(6) REVAL = INTGRS(7).EQ.1 EGROW = REALS(2) ERREST = REALS(3) H = REALS(4) HFACTS = REALS(5) HLAST = REALS(6) XRES = REALS(7) XLAST = REALS(8) XJAC = REALS(9) C FI 095 CONTINUE C FI 100 DO 105 I=1,7 IF(LIMITS(I).GE.0)GO TO 101 GO TO 105 C THEN 101 INTGRS(LIMMEM+I)=LIMITS(I) LIMITS(I) = -1 C FI 105 CONTINUE C OD DO 110 I=1,2 110 IF(INTGRS(LIMMEM+I).GT.0.AND.LIMITS(I)+ + INTGRS(LIMMEM+I).LT.0)INDEX=0 C OD IF ( INDEX .GT. 3 ) STP=0 IF(INDEX.EQ.0)RETURN GO TO (140,860,915,760,760,820,885),INDEX C COMMENT 120: INITIALIZE CONSTANTS AND STARTING VALUES C 120 MAXI=7 TEMP = 1.E0 URND=1.E0 125 URND = URND/2.E0 IF (1.E0 + URND/2.E0 .GT. TEMP) GO TO 125 REALS(1) = URND INTGRS(MAXI + 1) = 0 INTGRS(MAXI + 2) = INTGRS(MAXI + 1) + 1 INTGRS(MAXI + 3) = INTGRS(MAXI + 2) + 2 INTGRS(MAXI + 4) = INTGRS(MAXI + 3) + 2 INTGRS(MAXI + 5) = INTGRS(MAXI + 4) + 3 INTGRS(MAXI + 6) = INTGRS(MAXI + 5) + 3 INTGRS(MAXI + 7) = INTGRS(MAXI + 6) + 4 INTGRS(MAXI + 8) = INTGRS(MAXI + 7) + 4 INTGRS(MAXI + 9) = INTGRS(MAXI + 8) + 4 INTGRS(MAXI + 10) = INTGRS(MAXI + 9) + 5 INTGRS(MAXI + 11) = INTGRS(MAXI + 10) + 5 INTGRS(MAXI + 12) = INTGRS(MAXI + 11) + 5 INTGRS(MAXI + 13) = INTGRS(MAXI + 12) + 6 INTGRS(MAXI + 14) = INTGRS(MAXI + 13) + 6 INTGRS(MAXI + 15) = INTGRS(MAXI + 14) + 6 INTGRS(MAXI + 16) = INTGRS(MAXI + 15) + 6 DO 130 I=1 , N 130 YMAG(I)=0.E0 C OD M = 10 DO 135 I=1,7 135 LIMITS(I) = 0 LIMITS(2) = 1 STPMAX=0.E0 STP=0.E0 INDEX=1 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 140 H=STP XLAST=X DO 145 I=1 , N 145 IF( ABS(Y(I)).GT.YMAG(I) ) YMAG(I)=ABS(Y(I)) C OD XRES=0.E0 DO 150 J=1 , N 150 REALS(YRES+J)=0.E0 C OD IF(H.EQ.0.E0) GOTO 155 GOTO 175 C THEN 155 ENORM=0.E0 DO 160 I=1 , N 160 IF( YMAG(I).GT.EPS(I)*ENORM ) ENORM=YMAG(I)/EPS(I) C OD IF(ENORM.NE.0.E0) GOTO 165 GOTO 175 C THEN 165 CALL F(X,Y,REALS(YTEMP+1)) DO 170 I = 1 , N 170 IF(REALS(YTEMP+I).NE.0.E0.AND. (H.EQ.0.E0 + .OR.H*ABS(REALS(YTEMP+I)).GT.1.25E0*EPS(I)*ENORM-YMAG(I))) + H=(1.25E0*EPS(I)*ENORM-YMAG(I))/ABS(REALS(YTEMP+I)) C OD C FI C FI 175 IF( H.EQ.0.E0 ) H=XPRINT-X IF(ABS(H).GT.STPMAX .AND. STPMAX.GT.0.E0) H=STPMAX IF(H*(XPRINT-X).LT.0.E0) H=-H MLAST = 0 M=1 REVAL = .TRUE. NJAC=0 STP=0.E0 MSTATE = 0 MTRAN=0 HFACTS=0.E0 HLAST=0.E0 ERREST=1.E0 XJAC=X C COMMENT 180: COMPUTE LAGUERRE ZEROS AND RELATED CONSTANTS C REALS(XI+1) =1 REALS(EFACTS+1)=0.5E0 DO 205 M=2 , MAXM DO 205 I=1 , M K= M * (M-1) / 2 IF (I.EQ.1) GO TO 185 IF (I.EQ.M) GO TO 190 TEMP= (REALS(XI+K-M+I) + REALS(XI+K-M+I+1))/2.E0 GO TO 195 185 TEMP= REALS(XI+K-M+2) /2.E0 GO TO 195 190 TEMP= REALS(XI+K)*(1.E0+1.E0/(M-1.9E0)) 195 P= 1 PPRIME= 0.E0 DO 200 J= 1 , M PPRIME= PPRIME - P 200 P= P + TEMP * PPRIME/J C OD DELTA=P/PPRIME TEMP=TEMP-DELTA IF( (DELTA/TEMP) ** 2 .GE. .01E0 * URND )GO TO 195 REALS(XI+K+I) = TEMP 205 IF( I.LE.INTGRS(MAXI+M+1)-INTGRS(MAXI+M) ) + REALS(EFACTS+INTGRS(MAXI+M)+I)=ABS(PPRIME*TEMP)/(M+1) C OD C OD C COMMENT 220: PERFORM STEPS C M = 1 NSTEP=0 C DO 220 NSTEP=NSTEP+1 IF( H*(XPRINT-X).LT.0.E0 ) GO TO 225 GO TO 230 C THEN 225 MSTATE=0 H=-H MLAST=0 C FI 230 MROOTS=M*(M-1) / 2 REVAL = REVAL.AND.MTYPE .GT. 1 TEMP=HFACTS/H IF(.NOT.REVAL.AND.MTYPE.GT.1) + REVAL=TEMP.GT.1.5625E0.OR.TEMP.LT.0.64E0.OR. NSTEP.GE.NJAC + 20 IF(MTYPE.LE.1)HFACTS=0.E0 IF (REVAL) GOTO 240 GOTO 360 C THEN C COMMENT 240: RE-EVALUATE JACOBIAN IF APPROPRIATE C 240 IF( MTYPE.EQ.2 ) GO TO 241 GO TO 245 C THEN 241 CALL JAC(X,Y,REALS(FACTS+N+1)) GO TO 265 C ELSE 245 SQURND=SQRT(URND) CALL F(X,Y,REALS(W+1)) D = 0.E0 DO 250 I = 1 , N 250 D = D + ABS(REALS(W+I))/EPS(I) C OD D = D * ABS(H) * URND * 1.E-3 DO 260 J = 1 , N DELTA = D * EPS(J) + SQURND * YMAG(J) TEMP = Y(J) Y(J) = TEMP + DELTA DELTA = Y(J) - TEMP CALL F(X,Y,REALS(DERIV+1)) DO 255 I = 1 , N REALS(FACTS+N*J+I) = 0.E0 255 IF(DELTA .NE. 0.E0) REALS(FACTS+N*J+I)= + (REALS(DERIV + I)-REALS(W + I))/DELTA 260 Y(J) = TEMP C OD C OD C FI 265 XJAC=X NJAC=NSTEP REVAL=.FALSE. C COMMENT 280: FACTORIZE ( IDENTITY - H * JACOBIAN ) C HFACTS= H DO 290 I= 1 , N DO 285 J= 1 , N 285 REALS(FACTS+N*J+I)=-H*REALS(FACTS+N*J+I) C OD REALS(FACTS+N*I+I)=1+REALS(FACTS+N*I+I) 290 INTGRS(IPERM+I)= I C OD DO 350 K= 1 , N IF(.NOT.REVAL) GOTO 291 GOTO 350 C THEN 291 PIVOT= 0.E0 IPIVOT= 0 IF(K .GT. 1) GOTO 295 GOTO 330 C THEN 295 DO 325 I =1,N TEMP= REALS(FACTS+N*K+I) IF(I .LT. K) GOTO 296 GOTO 310 C THEN 296 IF(I.EQ.1)GO TO 305 IM1 = I-1 DO 300 J=1 , IM1 300 TEMP= TEMP-REALS(FACTS+N*J+I)*REALS(FACTS+N*K+J) C OD 305 REALS(FACTS+N*K+I)= TEMP/REALS(FACTS+N*I+I) GO TO 325 C ELSE 310 KM1 = K-1 IF(K.EQ.1)GO TO 320 DO 315 J=1 , KM1 315 TEMP= TEMP-REALS(FACTS+N*J+I)*REALS(FACTS+N*K+J) C OD 320 REALS(FACTS+N*K+I)= TEMP IF( ABS(TEMP/EPS(INTGRS(IPERM+I))) .GT. PIVOT ) + GOTO 321 GOTO 325 C THEN 321 PIVOT= ABS(TEMP/EPS(INTGRS(IPERM+I))) IPIVOT= I C FI C FI 325 CONTINUE C OD GO TO 350 C ELSE 330 DO 335 I= 1 , N IF( ABS(REALS(FACTS+N+I)/EPS(I)) .GT. PIVOT ) + GOTO 331 GOTO 335 C THEN 331 PIVOT= ABS(REALS(FACTS+N+I)/EPS(I)) IPIVOT= I 335 CONTINUE C FI C OD IF( IPIVOT .EQ. 0 ) GOTO 336 GOTO 340 C THEN 336 H = 0.9E0*H MSTATE = 1 GO TO 240 C ELSE 340 IF( IPIVOT .NE. K )GO TO 341 GOTO 350 C THEN 341 J= INTGRS(IPERM+IPIVOT) INTGRS(IPERM+IPIVOT)= INTGRS(IPERM+K) INTGRS(IPERM+K)= J DO 345 J= 1 , N TEMP = REALS(FACTS+N*J+IPIVOT) REALS(FACTS+N*J+IPIVOT) = REALS(FACTS+N*J+K) 345 REALS(FACTS+N*J+K) = TEMP C OD C FI C FI C FI 350 CONTINUE C FI C OD C FI C COMMENT 360: COMPUTE STARTING VALUES FOR ITERATIONS C 360 IF( STP .NE. 0.E0 .OR. H .NE. HLAST )GO TO 361 GOTO 440 C THEN 361 IF( MLAST.GT.M ) MLAST=M IF( MLAST.EQ.0 )GO TO 365 GOTO 375 C THEN 365 MP1 = M+1 DO 370 I=1 , MP1 DO 370 J=1 , N 370 REALS(YSPACE+N*I+J)=0.E0 C OD C OD GO TO 440 C ELSE 375 MP1 = M+1 DO 425 KP1 = 1,MP1 K = KP1 - 1 D = 0.E0 IF(K.NE.0)D=H*REALS(XI+MROOTS+K) D = (STP + D)/HLAST P=1.E0 PPRIME=0.E0 TEMP=0.E0 IF( K.EQ.0 )GO TO 380 GOTO 390 C THEN 380 DO 385 J=1,N REALS(YSPACE+N*M+N+J)=0.E0 385 REALS(YTEMP+J)=0.E0 C OD GO TO 400 C ELSE 390 DO 395 J=1 , N 395 REALS(YSPACE+N*K+J)=-REALS(YTEMP+J) C OD C FI 400 DO 425 I=1 , MLAST PPRIME=PPRIME-P TEMP=(D/I)*PPRIME P=P+TEMP IF( K.EQ.0 )GO TO 405 GOTO 415 C THEN 405 DO 410 J = 1,N 410 REALS(YTEMP+J)=REALS(YTEMP+J)-TEMP * + REALS(FSPACE+N*I+J) C OD GO TO 425 C ELSE 415 DO 420 J=1 , N 420 REALS(YSPACE+N*K+J)= + REALS(YSPACE+N*K+J)-TEMP*REALS(FSPACE+N*I+J) C OD C FI C OD C OD C FI C FI 425 CONTINUE C COMMENT 440: GENERATE TRANSFORMATION MATRICES T AND TINV C WHEN NECESSARY C 440 IF( M .NE. MTRAN )GO TO 441 GOTO 460 C THEN 441 MTRAN = M DO 455 I= 1 , M P= 1.E0 PPRIME= 0.E0 TEMP= REALS(XI+MROOTS + I) REALS(T+MAXM+I)= P IF(M.EQ.1)GO TO 450 DO 445 J= 2 , M PPRIME= PPRIME - P P= P + TEMP*PPRIME/(J-1) 445 REALS(T+MAXM*J+I)= P C OD 450 DO 455 K= 1 , M 455 REALS(TINV+MAXM*I+K)= TEMP*REALS(T+MAXM*K+I)/(M*P)**2 C OD C OD C FI C COMMENT 460: CARRY OUT ITERATIONS C 460 NITS = 0 CONVGD = .FALSE. DIVGD = .FALSE. EXCESS = .FALSE. 465 NITS = NITS + 1 IF(.NOT.(CONVGD.OR.DIVGD.OR.EXCESS))GO TO 480 GOTO 740 C THEN C COMMENT 480: EVALUATE DERIVATIVES AND TRANSFORM USING TINV C 480 ISAVE=INTGRS(MAXI+M+1)-INTGRS(MAXI+M) DO 485 J =1 , N 485 REALS(YSAVE+J)=REALS(YSPACE+N*ISAVE+J) C OD MP1 = M+1 DO 495 I=1 , MP1 DO 490 J=1 , N 490 REALS(YTEMP+J)= Y(J) + REALS(YSPACE+N*I+J) C OD TEMP = X IF (I .LE. M) TEMP = X+H*REALS(XI+MROOTS+I) CALL F(TEMP,REALS(YTEMP+1),REALS(DERIV+1)) DO 495 K=1 , N 495 REALS(FSPACE+N*I+K)= H*REALS(DERIV+K) C OD C OD DO 510 J=1 , N DO 505 I=1 , M TEMP=0.E0 TEMPF=0.E0 DO 500 K=1 , M TEMP=TEMP + REALS(TINV+MAXM*K+I)*REALS(YSPACE+N*K+J) 500 TEMPF=TEMPF+ REALS(TINV+MAXM*K+I)*REALS(FSPACE+N*K+J) C OD REALS(YVECT+I)= TEMP 505 REALS(FVECT+I)= TEMPF C OD DO 510 I=1 , M REALS(YSPACE+N*I+J)= REALS(YVECT+I) 510 REALS(FSPACE+N*I+J)= REALS(FVECT+I) C OD C OD C COMMENT 520: CHOOSE ITERATION OPTION C IF( MTYPE.LT.2 ) GO TO 540 GOTO 560 C THEN C COMMENT 540: FIXED POINT ITERATION C 540 DO 550 J=1 , N TEMP=REALS(FSPACE+N+J) REALS(YSPACE+N+J)=TEMP IF(M.EQ.1)GO TO 550 DO 545 I=2 , M REALS(YSPACE+N*I+J)= + REALS(FSPACE+N*I+J)-REALS(FSPACE+N*I-N+J) 545 TEMP=TEMP + REALS(FSPACE+N*I+J) C OD 550 REALS(YSPACE+N*M+N+J)=REALS(FSPACE+N*M+N+J) - TEMP C OD GO TO 660 C ELSE C COMMENT 560: NEWTON ITERATION C 560 REL= 2.E0/(1.E0+H/HFACTS) MP1 = M+1 DO 640 I=1 , MP1 DO 570 J=1 , N K=INTGRS(IPERM+J) TEMP= REALS(YSPACE+N*I+K)-REALS(FSPACE+N*I+K) IF(I .NE. 1) TEMP = TEMP+REALS(YLAST+K) IF(J.EQ.1)GO TO 570 JM1 = J-1 DO 565 K=1 , JM1 565 TEMP=TEMP - REALS(FACTS+N*K+J) * REALS(YTEMP+K) C OD 570 REALS(YTEMP+J)= TEMP/REALS(FACTS+N*J+J) C OD DO 640 NP1MJ = 1 , N J = N + 1 -NP1MJ TEMP=REALS(YTEMP+J) IF(J.EQ.N)GO TO 580 JP1 = J + 1 DO 575 K = JP1 , N 575 TEMP = TEMP - REALS(FACTS+N*K+J) * REALS(YTEMP+K) C OD 580 REALS(YTEMP+J)=TEMP IF(I.LE.M) GO TO 585 GO TO 600 C THEN 585 IF(I.EQ.1) GO TO 590 GO TO 595 C THEN 590 REALS(FTOTAL+J)=REALS(FSPACE+N+J) GO TO 600 C ELSE 595 REALS(FTOTAL+J)=REALS(FSPACE+N*I+J)+ + REALS(FTOTAL+J) C FI C FI 600 IF(I.EQ.1 .OR. I.EQ.M+1) GOTO 601 GOTO 605 C THEN 601 REALS(W+J) = TEMP GOTO 610 C ELSE 605 REALS(W+J)=TEMP+REALS(W+J) C FI 610 IF(I.LE.M) GOTO 611 GOTO 625 C THEN 611 IF(I.EQ.M) GOTO 615 GOTO 620 C THEN 615 REALS(YLAST+J)=REALS(FTOTAL+J) GOTO 625 C ELSE 620 REALS(YLAST+J)=REALS(FSPACE+N*I+J) + -REALS(W+J) C FI C FI 625 REALS(YSPACE+N*I+J)=REALS(YSPACE+N*I+J) + -REL*REALS(W+J) IF(I.LE.M) GOTO 630 GOTO 640 C THEN 630 IF(I.EQ.1) GOTO631 GOTO 635 C THEN 631 REALS(FSPACE+N*I+J)=REALS(YSPACE+N*I+J) GOTO 640 C ELSE 635 REALS(FSPACE+N*I+J)=REALS(YSPACE+N*I+J) + +REALS(FSPACE+N*I-N+J) C FI 640 CONTINUE C FI C OD C OD C FI 540 C COMMENT 660: RETRANSFORM USING T C 660 DO 675 J=1 , N DO 670 I=1 , M TEMP = REALS(YSPACE+N+J) IF(M.EQ.1)GO TO 670 DO 665 K=2 , M 665 TEMP= TEMP + REALS(T+MAXM*K+I)*REALS(YSPACE+N*K+J) 670 REALS(YVECT+I)= TEMP C OD C OD DO 675 I=1 , M 675 REALS(YSPACE+N*I+J)=REALS(YVECT+I) C OD C OD C COMMENT 680: TEST CONVERGENCE C ITLIM = MAXITS + M - MLAST IF( NITS.GT.1 ) XNORM=ENORM ENORM=0.E0 DO 685 J=1 , N TEMP= + ABS(REALS(YSPACE+N*ISAVE+J)-REALS(YSAVE+J))/EPS(J) 685 IF( TEMP .GT. ENORM ) ENORM=TEMP C OD IF( MERPST.EQ.0 ) + ENORM=ENORM/(ABS(H)*REALS(XI+MROOTS+ISAVE)) IF(MTYPE.EQ.0 .OR. NITS.EQ.1) GOTO 690 GOTO 695 C THEN 690 DIVGD=.FALSE. GOTO 700 C ELSE 695 DIVGD=ENORM .GT. XNORM + 1.0E0 C FI 700 IF(DIVGD) GOTO 701 GOTO 705 C THEN 701 CONVGD=.FALSE. EXCESS=.FALSE. GOTO 725 C ELSE 705 IF(ITLIM.EQ.NITS .AND. MTYPE.EQ.0) GOTO 706 GOTO 710 C THEN 706 CONVGD=.TRUE. GOTO 720 C ELSE 710 IF(NITS.EQ.1) GOTO 711 GOTO 715 C THEN 711 CONVGD=ENORM.LT.1.E0 GOTO 720 C ELSE 715 CONVGD=ENORM**2.LT.XNORM.OR.ENORM.LT.1.E0 C FI C FI 720 IF(.NOT.CONVGD) EXCESS=ITLIM.EQ.NITS C FI 725 GO TO 465 C FI 480 C COMMENT 740: CONVERGENCE FAILURE: TAKE REMEDIAL ACTION C 740 IF( .NOT. CONVGD ) GO TO 741 GOTO 780 C THEN 741 STP=H*REALS(XI+MROOTS+ISAVE) IF( MTYPE.LT.2.OR.X.EQ.XJAC )GO TO 745 GOTO 755 C THEN 745 IF (DIVGD) GOTO 746 GOTO 750 C THEN 746 H=0.9E0*H*XNORM**(1.E0-1.E0/ITLIM)/ENORM GOTO 755 C ELSE 750 H=0.9E0*H*ENORM**(-1.E0/ITLIM) C FI C FI 755 IF( DIVGD ) MLAST=0 REVAL = MTYPE.GT.1 MSTATE = 0 INDEX = 5 IF(DIVGD)INDEX = 4 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 760 GO TO 975 C ELSE C COMMENT 780: CONVERGENCE SUCCESS: SELECT APPROPRIATE ORDER C 780 DO 785 J = 1 , N 785 REALS(FSPACE+N*M+N+J) = REALS(YSPACE+N*M+N+J) C OD STP=0.E0 DO 810 MP1MI = 1 , M I = M + 1 - MP1MI IF (STP .NE. 0.E0 )GO TO 815 C THEN C ESCAPE FROM DO LOOP C ELSE ENORM = 0.E0 DO 790 J=1 , N 790 IF(ABS(REALS(FSPACE+N*I+N+J)).GT.EPS(J)*ENORM) + ENORM = ABS(REALS(FSPACE+N*I+N+J))/EPS(J) C OD TEMP = REALS(EFACTS+INTGRS(MAXI+I+1))+1 LBD=INTGRS(MAXI+I+1)-INTGRS(MAXI+I) DO 805 MAXMK =1 , LBD K = LBD + 1 - MAXMK IF( REALS(EFACTS+INTGRS(MAXI+I)+K) .GT. TEMP + .OR.STP.NE.0.E0)GO TO 810 C THEN C ESCAPE FROM LOOP C ELSE TEMP = REALS(EFACTS+INTGRS(MAXI+I)+K) C = TEMP IF( MERPST .EQ. 0 ) C = C/ABS(H*REALS(XI+MROOTS+K)) OK=I.EQ.M.AND.K .EQ. INTGRS(MAXI+I + 1) + -INTGRS(MAXI+I) IF( OK )GO TO 795 GOTO 800 C THEN 795 EGROW = 1.5E0 IF(C*ENORM .LT. 1.3E0*ERREST) EGROW=C*ENORM/ERREST IF(EGROW .GT. 1.4E0) MSTATE = 0 IF( EGROW .LT. 1.E0 ) EGROW = 1.E0 C FI 800 IF( ENORM*C.LT.1.E0 ) GO TO 801 GOTO 805 C THEN 801 STP = H*REALS(XI+MROOTS+K) MNEW = I IF( .NOT.OK ) MSTATE = 0 C FI C FI C OD 805 CONTINUE C FI 810 MROOTS=MROOTS-I+1 C OD 815 IF(STP.EQ.0.E0) GO TO 816 GOTO 825 C THEN 816 STP=H*REALS(XI+(M*(M-1)/2)+ISAVE) MSTATE=0 REVAL=MTYPE.GT.0 INDEX = 6 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 820 GO TO 840 C ELSE 825 M=MNEW C FI C COMMENT 840: INTERPOLATE OUTPUT VALUE IF REQUIRED C 840 IF(STP.NE.0.E0) GO TO 841 GOTO 880 C THEN 841 IF((X+1.01*STP-XPRINT)*(XPRINT-X).GT.0.E0)GO TO 842 GOTO 880 C THEN 842 DO 845 I=1 , N 845 REALS(YTEMP+I)=Y(I) D=(XPRINT-X)/H P=D PPRIME=0.E0 DO 850 K=1 , M PPRIME=PPRIME-P P=P+D*PPRIME/K DO 850 I=1 , N 850 REALS(YTEMP+I)=REALS(YTEMP+I)- + PPRIME*REALS(FSPACE+N*K+I)/K C OD C OD DO 855 I=1,N TEMP = Y(I) Y(I) = REALS(YTEMP+I) 855 REALS(YTEMP+I) = TEMP C OD TEMP = X X = XPRINT XPRINT=XPRINT+(XPRINT-XLAST) XLAST=TEMP INDEX=2 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 860 TEMP = XLAST XLAST = X X = TEMP DO 865 I=1,N 865 Y(I) = REALS(YTEMP+I) C OD GO TO 841 C FI C FI C COMMENT 880: COMPLETE COMPUTATION C 880 IF((XPRINT-X)*H.LT.0.E0)GO TO 881 GOTO 890 C THEN 881 STP = H*REALS(XI+(M*(M-1)/2)+INTGRS(MAXI+M+1)- + INTGRS(MAXI+M)) MSTATE = 0 INDEX = 7 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 885 CONTINUE C FI 890 IF(STP.NE.0.E0) GO TO 891 GOTO 920 C THEN 891 D=STP/H P=D PPRIME=0.E0 DO 895 K=1 , M PPRIME=PPRIME-P P=P+D*PPRIME/K DO 895 I=1 , N 895 REALS(YRES+I)= + REALS(YRES+I) - PPRIME*REALS(FSPACE+N*K+I)/K C OD C OD XRES=XRES+STP TEMP=X X=X+XRES XRES=XRES-(X-TEMP) DO 900 I=1 , N TEMP=Y(I) Y(I)=TEMP+REALS(YRES+I) REALS(YRES+I)=REALS(YRES+I)-(Y(I)-TEMP) YMAG(I)= .9E0* YMAG(I) 900 IF( ABS(Y(I)).GT.YMAG(I) ) YMAG(I) = ABS(Y(I)) C OD INDEX=3 GO TO 060 C GO TO AND RETURN FROM PAUSE ROUTINE 915 CONTINUE C FI C COMMENT 920: SPECIFY ORDER AND H FOR NEXT STEP C 920 IF(H.LE.0.E0) D=-STPMAX IF(H.GT.0.E0) D=STPMAX IF(STPMAX .EQ. 0.E0) D=2.E2*H IF(STP.NE.0.E0 .AND. ABS(D).GT.5.E0*ABS(STP)) + D=5.E0*STP MNEW=M IF( MSTATE .EQ. 4) GO TO 925 GOTO 935 C THEN 925 MNEW = M+1 DO 930 J=1,N 930 REALS(FSPACE+N*(M+2)+J)= + (REALS(FSPACE+N*(M+2)+J)-REALS(FSPACE+N*M+N+J))*H/STP C OD C FI 935 MLAST = M HLAST = H MROOTS = MNEW*(MNEW+1) / 2 SCORE = 0.E0 DO 950 MNP1MI = 1, MNEW I = MNEW + 1 - MNP1MI MROOTS = MROOTS-I XIM = REALS(XI+MROOTS+INTGRS(MAXI+I+1)-INTGRS(MAXI+I)) ENORM = 0.E0 HERR = D/XIM DO 940 J = 1 , N 940 IF(ABS(REALS(FSPACE+N*I+N+J)).GT.EPS(J)*ENORM) + ENORM = ABS(REALS(FSPACE+N*I+N+J))/EPS(J) C OD ENORM = REALS(EFACTS+INTGRS(MAXI+I+1))*ENORM IF( MERPST.EQ.0 ) ENORM=ENORM/(ABS(H)*XIM) TEMP =1.4E0 IF( I .EQ. MLAST ) TEMP = 2.E0 IF(I.EQ.MLAST.AND.MSTATE.GE.2)TEMP=1.2E0*EGROW TEMP = TEMP * ENORM IF( TEMP.GT.(H/HERR)**(I+MERPST) ) + HERR=H/TEMP**(1.E0/(I+MERPST)) TEMP=ABS(HERR)*XIM/(I+1) IF( TEMP .GT. SCORE ) GO TO 945 GOTO 950 C THEN 945 IF( SCORE.GT.0.E0 ) MSTATE = 0 SCORE=TEMP HNEW=HERR M=I ERREST = ENORM*(HNEW/H)**(M+MERPST) C FI 950 CONTINUE C OD IF(MSTATE.NE.3.OR.(M.LT.MAXM.AND.EGROW.LT.1.2E0)) + GO TO 955 GOTO 960 C THEN 955 MSTATE = MSTATE + 1 IF(MSTATE .EQ. 5)MSTATE=2 C FI 960 IF(MSTATE .EQ. 4) GO TO 961 GOTO 970 C THEN 961 DO 965 J=1,N 965 REALS(FSPACE+N*(M+2)+J)=REALS(FSPACE+N*M+N+J) C OD GO TO 975 C ELSE 970 H=HNEW C FI C FI 740 975 GO TO 220 C OD END