From: IN%"kjudd@bucky.Stanford.EDU" 29-JUN-1995 13:38:46.42 To: IN%"judd@HOOVER.STANFORD.EDU" CC: Subj: Dear Colleague, Following are three files which I promised to deliver on request in my 1992 JET paper on Projection Methods. I suggest that you separate the three files, juddread.me, juddgroh.for, and JUDD.IN. The first is a general information file, the second is the FORTRAN code, and the third is a sample input file for the FORTRAN program. If you have any problems, please tell me about them. Ken Judd 415/723-5866 judd@hoover.stanford.edu ############################################################################ ############################################################################ FILE ONE: JUDDREAD.ME ############################################################################ README FILE FOR JUDDGROH.FOR **************************************************************************** **************************************************************************** FORTRAN PROGRAM: JUDDGROH.FOR **************************************************************************** The file juddgroh.for contains all the necessary FORTRAN programs for computing the examples in my JET 1992 paper. Just compile this file and it will read input values from JUDD.IN and produce the output files discussed below. Part of the program file consists of the nonlinear equation solver program which I call HYBJUD. It is a modified version of the nonlinear equation solvers in MINPACK. You may find this nonlinear equation solver to be of independent use to you. Warnings I have tried to stay within FORTRAN77 standards. I do use DO-ENDDO and IF-THEN-ELSEIF constructs, but most compilers support these extensions, and they are part of FORTRAN90. I used the Lahey compilers. The only compiler-specific statement I used was the clock call subroutine, TIME. They appear only in the main program and should be either commented out or replaced by your compiler's clock call. You can find these instances by searching for 'WARNING ON TIME'. **************************************************************************** **************************************************************************** INPUT FILE: JUDD.IN **************************************************************************** Iopt switches whether you use a numerical or analytical Jacobian: Iopt = 1 chooses the analytical Jacobian Iopt = 2 chooses the numerical Jacobian Tol sets the convergence criterion for the nonlinear equation solver. The smaller, the tighter the criterion. Choose values between .01 and .00000001. nchk: The number of capital stocks (uniformly distributed) used to check the accuracy mchk: The number of productivity shocks used to check accuracy ncases: The number of (nk,nt,ngrid,mgrid) combinations inputted (<11) nk: Number of basis elements used in the capital dimension (<10) nt: Number of basis elements in productivity direction(<10) ngrid: Number of capital points used in projection equations (<50) mgrid: Number of productivity values used in projections (<30) (If you violate these bounds, the program may crash) numrho: number of rho's (persistence parameter in productivity process) (<5) list of rho's numsig number of sigma's (variance of productivity shocks) (<5) list of sigma's numgam number of gamma's (risk aversion parameter, negative) (<7) list of gamma's ibasis, igrid ibasis determines whether we use Chebyshev polynomials or ordinary polynomials: 1 chooses Chebyshev 2 chooses ordinary polynomials igrid determines whether we use the Chebyshev zeroes as the grid points, or a uniform grid: 1 chooses Chebyshev zeroes 2 chooses uniform grid capmin, capmax capmin is the minimum capital examined, capmax the maximum. 0 < capmin < 1 < capmax. ithmethod, thmin, thmax ithmethod=1 says that the range of productivity shocks is determined by the maximum possible range given rho and three standard deviation productivity shocks; thmin and thmax are ignored. ithmethod=2 sets the productivity shock range be [thmin,thmax]. Use this to implement Taylor-Uhlig style settings. alpha, beta alpha is capital share; beta is discount factor incase, lcase, step incase is the first (nk,nt,ngrid,mgrid) case solved lcase is the last case step is the step between cases Therefore, we start with incase, then go to incase+step, incase+2*step, etc. until we get past lcase. Below is an input file which works fine. Keep it in case you mess up your JUDD.IN file: 1 .00001 iopt tol 20 10 nchk mchk 8 ncases 2 4 4 7 10 10 15 15 nk 2 3 3 5 6 6 10 10 mt 2 4 7 20 10 25 15 25 ngrd 2 3 5 12 6 12 10 30 mgrd 2 number of sigma's .01 .03 .04 sigma 3 number of rho's .8 .5 .3 rho 8 number of gamma's -15.0 -10.0 -7.0 -5.0 -2.0 -.9 -.5 -.1 gamma 1 1 ibasis igrid .333333 2.0 capmin capmax 1 .999 1.001 ithmethod thmin thmax .33333333333 .9500000 alpha beta 1 4 1 incase lcase step In this file, the program first chooses nk=2, mt=2, ngrid=2, mgrid=2, and then solves for the growth model for only two sigma's, .01 and .03, three rho's, eight gamma's, chooses Chebyshev polynomial basis and Chebyshev zeroes grid, a capital range of [.3,2.0], the wide range of productivity shocks, alpha = .3, and beta=.95. This is a total of 48 growth models which are solved. The last line -- 1 4 1 -- causes these 48 growth models to be solved also for the (nk,nt,ngrid,mgrid) case of (4,3,4,3), (4,3,7,5), and (7,5,20,12). The other cases are read in only, but not executed. **************************************************************************** **************************************************************************** OUTPUT FILES **************************************************************************** JUDD.ERR This file reports the L-p errors of the Euler equation. For example, consider the following two lines: 4 3 7 5 -15.00 0.80 0.01 -3.17 -3.67 -3.79 -3.61 -3.80 This says that for an approximation where (nk,nt,ngrid,mgrid) = (4,3,7,5), and for gamma = -15, rho = .8, sigma = .01, the base-10 log of the maximum error is -3.17, the log of the L-2 norm is -3.67, and the log of the L-1 norm is -3.79. The last two numbers are the maximum error and log L-2 norm over a smaller range around the deterministic steady state. See the error table in the paper for a full discussion. JUDD.CFS This file reports the coefficients of the solution in the form i j a-sub-ij as in the example 1 1 0.16303764690142520000 2 1 0.04700083186822695000 1 2 0.02138972620641763000 2 2 0.00590755192547805100 JUDD.POL This file reports consumption, output, and savings at various choices of capital stock and productivity levels. JUDD.TIM 4 3 7 5 -15.00 0.80 0.01 08:45:48.21 08:45:49.25 12 2 1 0.7E+01 0.4E+01 The first three numbers are gamma, rho, and sigma. The first time is when the nonlinear equation solver is called, and the second time is when the program returns from that call. The next three integers are the number of function calls, the number of Jacobian calls, and the nature of the exit from SNSQE. The last two numbers are the condition numbers of the problem at the initial guess and at the solution. If you have any questions, call me at (415)723-5866. ############################################################################ ############################################################################ FILE TWO: JUDDGROH.FOR ############################################################################ * JUDDGROH.FOR * August 6, 1993 version * This program computes the equilibrium policy functions for an economy * with an autoregressive productive shock. We compute a Jacobian to help * the nonlinear equation solver. * * NK is the number of capital stock functions used, and NT is the number of * productivity shock functions. N is their product, and is the total * number of parameters in the policy function which we must compute. * * This version does a loop over NK,NT,NGRID, and MGRID. * It also computes a variety of measures of the error. * PARAMETER (NKMX = 10, NTMX = 10) PARAMETER (NMX = 100, LWA = 16000) PARAMETER (MXNGRD = 50, MXMGRD = 30, MXN = NKMX, MXM = NTMX) PARAMETER (PI = 3.14159265359D0) IMPLICIT REAL*8(A-H,O-Z) REAL*8 X(NMX), FVEC(NMX), WA(LWA), TOL INTEGER NKMAT(11),NTMAT(11),NGRDMT(11),MGRDMT(11) REAL*8 SAV(MXNGRD,MXMGRD) REAL*8 C(MXNGRD,MXMGRD),PROD(MXNGRD,MXMGRD) REAL*8 PHIM(MXN,MXNGRD),PSIM(MXM,MXMGRD) REAL*8 CAP(MXNGRD),EPS(20), & GHERMX(20),GHERMW(20), & THETA(MXMGRD) REAL*8 SIGMAT(5),RHOMAT(5),GAMMAT(7) REAL*8 PHIC(MXN),PSIT(MXM),DPHIC(MXN) REAL*8 WR(NMX),WI(NMX),RGSTOR(NMX),EIGVEC(2,2) REAL*8 FJAC(NMX,NMX) INTEGER IRGSTOR(NMX) CHARACTER*11 XTIMIN,XTIMOUT * COMMON /KOUNT/ IFCNCT, JACCT COMMON /BASIS/ IBASIS COMMON /LENPAR/ NK,NT COMMON /TECH/ AF,ALPHA COMMON /TASTE/ BETA,GAMMA COMMON /BLOCK1/ C,CAP,CAPSS,CAPMIN,CAPMAX COMMON /PHIBLK/ PHIM,PSIM,EPS,RHO,SIG,THETA COMMON /BLOCK2/ SAV,NGRID,MGRID,PROD COMMON /BLOCK3/ CAPBAR,CAPRNG,THBAR,THRNG,THMMIN,THMMAX EXTERNAL F,JAC * * We will use a Gauss-Hermite quadrature formula for computing expectations * in the stochastic Euler equation. * GHERMX are the evaluated points and GHERMW are the weights. * DATA GHERMX/-2.93063742,-1.98165676,-1.15719371,-.38118699, & .38118699,1.15719371,1.98165676,2.93063742, & 12*0.000/ DATA GHERMW/.00019960,.01707798,.20780233,.66114701, & .66114701,.20780233,.01707798,.00019960, & 12*0.000/ * * Declare input unit number IIN = 10 * We first input critical parameters OPEN(UNIT=IIN,STATUS='OLD',FILE='JUDD.IN') * IOPT and TOL are parameters for the nonlinear equation solver. * IOPT=1 says that a Jacobian routine is available. * TOL is the relative precision which HYBRD attains before stopping. READ(IIN,*)IOPT,TOL * NK is the dimension of the approximation in the capital dimension, * NT in the theta dimension. N is the number of unknown coefficients. * NGRID is the number of capital stocks at which we compute the residuals * used in the quadrature computations of projections. * MGRID is the number of productivity states which are used. * NCHK is the number of capital stocks used in final check. * MCHK is the number of theta values used in final check. READ(IIN,*)NCHK,MCHK READ(IIN,*) NCASES READ(IIN,*) (NKMAT(I),I=1,NCASES) READ(IIN,*) (NTMAT(I),I=1,NCASES) READ(IIN,*) (NGRDMT(I),I=1,NCASES) READ(IIN,*) (MGRDMT(I),I=1,NCASES) READ(IIN,*) NSIG READ(IIN,*) (SIGMAT(I),I=1,NSIG) READ(IIN,*) NRHO READ(IIN,*) (RHOMAT(I),I=1,NRHO) READ(IIN,*) NGAM READ(IIN,*) (GAMMAT(I),I=1,NGAM) * CAPMIN and CAPMAX are the minimum and maximum capital stocks * examined by the program. THMIN and THMAX are the extreme * values of the productivity shocks. ITHMETH tells us how to * compute the range of examined productivity levels. READ(IIN,*) IBASIS,IGRID READ(IIN,*) CAPMIN,CAPMAX READ(IIN,*) ITHMETH,THMMIN,THMMAX * ALPHA is capital share, BETA the discount factor, GAMMA the relative risk * aversion parameter, RHO the autocorrelation in productivity, and SIGMA * the standard deviation in the productivity shock. READ(IIN,*) ALPHA,BETA READ(IIN,*) INCASE,LCASE,ISTEP * * NPRINT is a parameter for HYBRD; zero suppresses printing. * NPRINT = 0 * AF is a scaling parameter factor for the production function, chosen * so that the deterministic steady state is 1.0 AF=(1.0D0/BETA - 1.0D0)/ALPHA CAPSS=1. CONSSS = GNP(1.0D0,1.0D0) *----------------------------------------------------------------------- * If you want to compute a monthly version, make the following * conversions. c BETA = BETA*(1.0D0/12.0D0) c AF=(1.0D0/BETA - 1.0D0)/ALPHA c RHO = RHO*(1.0D0/12.0D0) c SIG = SIG/DSQRT(12.0D0) * If you want to compute a quarterly version, make the following * conversions. c BETA = BETA*(1.0D0/4.0D0) c AF=(1.0D0/BETA - 1.0D0)/ALPHA c RHO = RHO*(1.0D0/4.0D0) c SIG = SIG/DSQRT(4.0D0) *----------------------------------------------------------------------- * CAPRNG and CAPBAR are set so that calls to Chebyshev polynomials * are between -1 and 1. CAPRNG=(CAPMAX-CAPMIN)/2.0D0 CAPBAR=CAPRNG+CAPMIN OPEN(UNIT=12,STATUS='UNKNOWN',FILE='JUDD.ERR') OPEN(UNIT=13,STATUS='UNKNOWN',FILE='JUDD.TIM') OPEN(UNIT=15,STATUS='UNKNOWN',FILE='JUDD.POL') OPEN(UNIT=14,STATUS='UNKNOWN',FILE='JUDD.CFS') WRITE(12,*)'ALPHA = ',ALPHA,' AF = ',AF WRITE(12,*)'BETA = ',BETA WRITE(12,*)'CAPMIN = ',CAPMIN,' CAPMAX = ',CAPMAX WRITE(12,*)'NCHK =', NCHK, 'MCHK = ', MCHK WRITE(12,*)'TOL =',TOL WRITE(13,*)'ALPHA = ',ALPHA,' AF = ',AF WRITE(13,*)'BETA = ',BETA WRITE(13,*)'CAPMIN = ',CAPMIN,' CAPMAX = ',CAPMAX WRITE(13,*)'NCHK =', NCHK, 'MCHK = ', MCHK WRITE(13,*)'TOL =',TOL WRITE(14,*)'ALPHA = ',ALPHA,' AF = ',AF WRITE(14,*)'BETA = ',BETA WRITE(14,*)'CAPMIN = ',CAPMIN,' CAPMAX = ',CAPMAX WRITE(14,*)'NCHK =', NCHK, 'MCHK = ', MCHK WRITE(14,*)'TOL =',TOL WRITE(15,*)'ALPHA = ',ALPHA,' AF = ',AF WRITE(15,*)'BETA = ',BETA WRITE(15,*)'CAPMIN = ',CAPMIN,' CAPMAX = ',CAPMAX WRITE(15,*)'NCHK =', NCHK, 'MCHK = ', MCHK WRITE(15,*)'TOL =',TOL * We compute the linear initial guess discussed in the paper. CALL PHI(1.0D0,3,PHIC,CAPMIN,CAPMAX) X1 = PHIC(1) X2 = PHIC(2) CALL PHI(0.0D0,3,PHIC,CAPMIN,CAPMAX) Y1 = PHIC(1) Y2 = PHIC(2) X2INIT = CONSSS/(X2 - Y2*X1/Y1) X1INIT = -X2INIT*Y2/Y1 DO 988 ICASE = INCASE,LCASE,ISTEP NK = NKMAT(ICASE) NT = NTMAT(ICASE) N = NK*NT NGRID = NGRDMT(ICASE) MGRID = MGRDMT(ICASE) WRITE(12,'(1X,4I6)')NK,NT,NGRID,MGRID WRITE(13,'(1X,4I6)')NK,NT,NGRID,MGRID DO 9999 IGAM=1,NGAM DO 9999 IRHO=1,NRHO DO 9999 ISIG=1,NSIG GAMMA=GAMMAT(IGAM) SIG=SIGMAT(ISIG) RHO=RHOMAT(IRHO) *----------------------------------------------------------------------- * We can also compute the steady-state linearization initial guess. * To activate this guess and use the one above, uncomment these * statements. C BBEIG = GAMMA-GAMMA/BETA+AF*ALPHA*(ALPHA-1.0D0)*AF*BETA C BBEIG = BBEIG/GAMMA C CCEIG = -(1.0D0+AF)*AF*ALPHA*(ALPHA-1.0D0)/GAMMA C SLOPE = (-BBEIG+DSQRT(BBEIG*BBEIG-4.0D0*CCEIG))/2.0D0 C CALL PHID(1.0D0,3,PHIC,DPHIC,CAPMIN,CAPMAX) C X2EIG = SLOPE/DPHIC(2) C X1EIG = AF-X2EIG*PHIC(2) C X1INIT = X1EIG C X2INIT = X2EIG *----------------------------------------------------------------------- * Use a linear starting guess for each of * the gamma-rho-sig combinations. c IHOT = 0 c if(isig.eq.1.and.irho.eq.1.and.igam.eq.1) then IHOT = 0 X(1) = X1INIT X(2) = X2INIT DO 3171 I = 3,N 3171 X(I) = 0.0D0 c endif * We compute the values of the productivity innovations which we will * consider when we compute the conditional expected utility term in * the stochastic Euler equation. DO 915 J=1,8 EPS(J)=DEXP(GHERMX(J)*SIG*DSQRT(2.0D0)) 915 CONTINUE EPSMIN = DEXP(-3.0D0*SIG*DSQRT(2.0D0)) EPSMAX = DEXP(3.0D0*SIG*DSQRT(2.0D0)) * We next compute the capital stocks which we will consider. Here we use * the Chebyshev nodes. IF(IGRID.EQ.1) THEN DO 9161 J=1,NGRID CAP(J)=CAPBAR-CAPRNG*COS(PI*(2*J-1)/(2*NGRID)) 9161 CONTINUE ELSEIF(IGRID.EQ.2) THEN DO 9162 J=1,NGRID CAP(J)=CAPMIN+(J-1)*(CAPMAX-CAPMIN)/(NGRID-1) 9162 CONTINUE ENDIF * THMMIN and THMMAX are the extreme values of the productivity * which will be called. THRNG and THBAR are parameters of the * PSI function. These parameters are set so as to keep the * argument of Chebyshev subroutine between -1 and +1. * If ITHMETH is 1 then we use EPSMIN and EPSMAX to compute * the range of THETA; otherwise, we read in THMMIN and THMMAX. IF(ITHMETH.EQ.1) THEN THMMIN = EPSMIN**(1.0D0/(1.0D0-RHO))*.9999D0 THMMAX = EPSMAX**(1.0D0/(1.0D0-RHO))*1.0001D0 ENDIF THRNG = (THMMAX-THMMIN)/2.0D0 THBAR = THRNG+THMMIN * We compute the productivity state values which we will examine * when we compute projections. DO 917 J=1,MGRID THETA(J)=THBAR-THRNG*COS(PI*(2*J-1)/(2*MGRID)) 917 CONTINUE * We compute output, integration weights, and basis function values * at the grid points, storing the results in matrices for easy * future access. This is done to minimize calls to PHI, PHID, PHIDD DO 918 J=1,NGRID DO 918 L=1,MGRID PROD(J,L)=GNP(CAP(J),THETA(L)) 918 CONTINUE DO 919 J=1,NGRID CALL PHI(CAP(J),NK,PHIC,CAPMIN,CAPMAX) DO 919 IK=1,NK PHIM(IK,J)=PHIC(IK) 919 CONTINUE DO 920 L=1,MGRID CALL PHI(THETA(L),NT,PSIT,THMMIN,THMMAX) DO 920 IT=1,NT PSIM(IT,L)=PSIT(IT) 920 CONTINUE * Compute the initial condition number. LDFJAC = N CALL JAC (N, X, FVEC, FJAC, LDFJAC, IFLAG) MATZ = 0 CALL RG(N,N,FJAC,WR,WI,MATZ,EIGVEC,IRGSTOR,RGSTOR,IERR) EMAX = 0 EMIN = 1.0D100 DO I=1,N EMAG = DSQRT(WR(I)*WR(I)+WI(I)*WI(I)) IF(EMAG.GT.EMAX) EMAX = EMAG IF(EMAG.LT.EMIN) EMIN = EMAG ENDDO CONDINIT = EMAX/EMIN IFCNCT = 0 JACCT = 0 * WARNING ON TIME * The following call to TIME is compiler specific. You should * either replace this call with your compiler's subroutine, * or comment out this call and uncomment the XTIMIN assignment. CALL TIME(XTIMIN) C XTIMIN = 'NO TIME ' WRITE(*,*) GAMMA,RHO,SIG,XTIMIN * * SOLVE NONLINEAR EQUATIONS * CALL HYBJUD (F, N, X, FVEC, FJAC, IHOT, TOL, INFO, WA, LWA) * WARNING ON TIME * The following call to TIME is compiler specific. You should * either replace this call with your compiler's subroutine, * or comment out this call and uncomment the XTIMIN assignment. CALL TIME(XTIMOUT) c XTIMIN = 'NO TIME ' WRITE(*,*) XTIMOUT,INFO,IFCNCT,JACCT LDFJAC = N CALL JAC (N, X, FVEC, FJAC, LDFJAC, IFLAG) MATZ = 0 CALL RG(N,N,FJAC,WR,WI,MATZ,EIGVEC,IRGSTOR,RGSTOR,IERR) EMAX = 0 EMIN = 1.0D100 DO I=1,N EMAG = DSQRT(WR(I)*WR(I)+WI(I)*WI(I)) IF(EMAG.GT.EMAX) EMAX = EMAG IF(EMAG.LT.EMIN) EMIN = EMAG ENDDO CONDIT = EMAX/EMIN WRITE(*,'(2E9.1)') CONDIT,CONDINIT WRITE(13,'(3F6.2,2(3X,A11),3I5,2E8.1)') GAMMA,RHO,SIG,XTIMIN, &XTIMOUT,IFCNCT,JACCT,INFO,CONDIT,CONDINIT * * PRINT RESULTS * IF (INFO .NE. 1) WRITE (*,810) INFO,GAMMA,RHO,SIG, & NK,NT,NGRID,MGRID 810 FORMAT (' INFO =', I3, 3F7.3, 4I5) WRITE(14,*)GAMMA,RHO,SIG WRITE(14,*)THMMIN,THMMAX,THRNG,THBAR WRITE(14,*)CAPMIN,CAPMAX,CAPRNG,CAPBAR DO 561,I=1,N IT=(I-1)/NK IK=I-IT*NK IT=IT+1 WRITE(14,59)IK,IT,X(I) 561 CONTINUE 59 FORMAT(1X,2I6,F25.20) CALL F(N, X, FVEC, IFLAG) WRITE(15,*) WRITE(15,*)'GAMMA = ',GAMMA, ' RHO = ',RHO WRITE(15,*)'SIGMA = ',SIG WRITE(15,*)'THMMIN = ',THMMIN,' THMMAX = ',THMMAX WRITE(15,*)'CAPMIN = ',CAPMIN,' CAPMAX = ',CAPMAX WRITE(15,*)' THETA CAP CONS SAV PROD ' DO 56,IT=1,MGRID DO 56,IK=1,NGRID WRITE(15,58)THETA(IT),CAP(IK),C(IK,IT),SAV(IK,IT),PROD(IK,IT) 58 FORMAT(1X,2F5.2,3F10.7) 56 CONTINUE CALL FCHK(N , X, NCHK, MCHK) 9999 CONTINUE 988 CONTINUE STOP 800 FORMAT (' INITIAL GUESS', /, 2F12.6) END ************************************************************************ SUBROUTINE F (N, X, FVEC, IFLAG) PARAMETER (NKMX = 10, NTMX = 10) PARAMETER (MXNGRD = 50, MXMGRD = 30, MXN = NKMX, MXM = NTMX) IMPLICIT REAL*8(A-H,O-Z) REAL*8 X(N), FVEC(N) REAL*8 SAV(MXNGRD,MXMGRD),A(MXNGRD,MXMGRD) REAL*8 C(MXNGRD,MXMGRD),PROD(MXNGRD,MXMGRD) REAL*8 PHIM(MXN,MXNGRD),PSIM(MXM,MXMGRD) REAL*8 CAP(MXNGRD),EPS(20), & THETA(MXMGRD) REAL*8 AMPK2(20),PSIM2(20,MXMGRD), & PHIM2(MXNGRD),C2(20), & GHERMX(20),GHERMW(20) REAL*8 PHIC(MXN),PSIT(MXM),DPHIC(MXN) * COMMON /KOUNT/ IFCNCT, JACCT COMMON /LENPAR/ NK,NT COMMON /TECH/ AF,ALPHA COMMON /TASTE/ BETA,GAMMA COMMON /BLOCK1/ C,CAP,CAPSS,CAPMIN,CAPMAX COMMON /PHIBLK/ PHIM,PSIM,EPS,RHO,SIG,THETA COMMON /BLOCK2/ SAV,NGRID,MGRID,PROD COMMON /BLOCK3/ CAPBAR,CAPRNG,THBAR,THRNG,THMMIN,THMMAX DATA GHERMX/-2.93063742,-1.98165676,-1.15719371,-.38118699, & .38118699,1.15719371,1.98165676,2.93063742, & 12*0.000/ DATA GHERMW/.00019960,.01707798,.20780233,.66114701, & .66114701,.20780233,.01707798,.00019960, & 12*0.000/ IFCNCT = IFCNCT + 1 * * -------------------------------------------------------------------- * F computes projections of the Euler equation error for a guess of the * consumption policy function. X is the vector of coefficients for * the basis functions on which we base our guesses. A nonlinear equation * package finds the X which sets all projections equal to zero. * --------------------------------------------------------------------- * * First we must decode X, a vector, to get the two-dimensional set of * coefficients, A(I,J), of our two-dimensional policy function. FVEC * is the vector of projections of the Euler equation errors. DO 15,I=1,N IT=(I-1)/NK IK=I-IT*NK IT=IT+1 A(IK,IT)=X(I) FVEC(I)=0.0D0 15 CONTINUE * We have picked the points, (CAP(I), THETA(J)), at which we compute the * Euler equation errors. First, using the coefficients given by X, compute the * consumption, CC, at that point, then savings, SAVV, and next * period's capital stock, CAP2. With this we can compute the * deviation from the Euler equation. DO 30,IT=1,MGRID DO 30,IK=1,NGRID CC=0. * Using the policy function expressed in A, we compute consumption * at the point represented by (IT,IK). DO 31,LK=1,NK DO 31,LT=1,NT CC=CC+A(LK,LT)*PHIM(LK,IK)*PSIM(LT,IT) 31 CONTINUE * Now compute savings and next period's capital stock. SAVV=PROD(IK,IT)-CC CAP2=CAP(IK)+SAVV C(IK,IT)=CC SAV(IK,IT)=SAVV * We compute the possible values of the tomorrow's productivity state, * and marginal product of capital, AMPK2(I) one value for each of the four * values of the innovation which we examine. DO 303 IT1=1,8 THETA2=THETA(IT)**RHO*EPS(IT1) AMPK2(IT1)=1.0D0+GNPK(CAP2,THETA2) * We compute the values of the basis functions at tomorrow's * possible capital and productivity values, to be used below to * compute tomorrow's consumption levels. CALL PHI(THETA2,NT,PSIT,THMMIN,THMMAX) DO 303 JT=1,NT PSIM2(JT,IT1)=PSIT(JT) 303 CONTINUE CALL PHI(CAP2,NK,PHIC,CAPMIN,CAPMAX) DO 304 JK=1,NK PHIM2(JK)=PHIC(JK) 304 CONTINUE * We compute the Euler equation error, VAL, at the point (IK,IT). VAL=0.0D0 DO 302 IT1=1,8 * We have to compute the next period's possible consumption values, * C2, one for each possible productivity innovation. C2(IT1)=0.D0 DO 301 JK=1,NK DO 301 JT=1,NT C2(IT1)=C2(IT1)+A(JK,JT)*PHIM2(JK)*PSIM2(JT,IT1) 301 CONTINUE VAL=VAL+UTILP(C2(IT1))*AMPK2(IT1)*GHERMW(IT1) 302 CONTINUE C The following is the usual form of the Euler equation C VAL=UTILP(CC)-BETA*VAL/DSQRT(3.14159D0) C Instead, we used the consumption equivalent form, as discussed C in the paper. VAL = CC - (BETA*VAL/DSQRT(3.14159D0))**(1.0D0/GAMMA) * We now update our running computation of each projection. DO 60 JK=1,NK DO 60 JT=1,NT J = JK + (JT-1)*NK FVEC(J)=FVEC(J)+VAL*PHIM(JK,IK)*PSIM(JT,IT) 60 CONTINUE 30 CONTINUE C WRITE(*,'(1X,4F14.6)') FVEC(1), FVEC(2), FVEC(3),FVEC(4) RETURN END ************************************************************************ * JAC * This program computes the gradients of F with respect to the components * of X. * SUBROUTINE JAC (N, X, FVEC, FJAC, LDFJAC, IFLAG) PARAMETER (NKMX = 10, NTMX = 10) PARAMETER (MXNGRD = 50, MXMGRD = 30, MXN = NKMX, MXM = NTMX) IMPLICIT REAL*8(A-H,O-Z) REAL*8 X(N), FVEC(N), FJAC(LDFJAC,N),GRAD(NKMX,NTMX) REAL*8 SAV(MXNGRD,MXMGRD),A(MXNGRD,MXMGRD) REAL*8 C(MXNGRD,MXMGRD),PROD(MXNGRD,MXMGRD) REAL*8 PHIM(MXN,MXNGRD),PSIM(MXM,MXMGRD) REAL*8 CAP(MXNGRD),EPS(20), & THETA(MXMGRD) REAL*8 AMPK2(20),AMPPK2(20),PSIM2(20,MXMGRD), & PHIM2(MXNGRD), & GHERMX(20),GHERMW(20),PHIXM2(MXNGRD) REAL*8 PHIC(MXN),PSIT(MXM),DPHIC(MXN) * COMMON /KOUNT/ IFCNCT, JACCT COMMON /LENPAR/ NK,NT COMMON /TECH/ AF,ALPHA COMMON /TASTE/ BETA,GAMMA COMMON /BLOCK1/ C,CAP,CAPSS,CAPMIN,CAPMAX COMMON /PHIBLK/ PHIM,PSIM,EPS,RHO,SIG,THETA COMMON /BLOCK2/ SAV,NGRID,MGRID,PROD COMMON /BLOCK3/ CAPBAR,CAPRNG,THBAR,THRNG,THMMIN,THMMAX DATA GHERMX/-2.93063742,-1.98165676,-1.15719371,-.38118699, & .38118699,1.15719371,1.98165676,2.93063742, & 12*0.000/ DATA GHERMW/.00019960,.01707798,.20780233,.66114701, & .66114701,.20780233,.01707798,.00019960, & 12*0.000/ JACCT = JACCT + 1 * We first duplicate the steps in F, as well as initialize the Jacobian * matrix, FJAC, to zero. DO 15,I=1,N IT=(I-1)/NK IK=I-IT*NK IT=IT+1 A(IK,IT)=X(I) DO 15 J=1,N FJAC(I,J)=0.D0 15 CONTINUE DO 30,IT=1,MGRID DO 30,IK=1,NGRID CC=0. DO 31,LK=1,NK DO 31,LT=1,NT CC=CC+A(LK,LT)*PHIM(LK,IK)*PSIM(LT,IT) GRAD(LK,LT)=0.0D0 31 CONTINUE SAVV=PROD(IK,IT)-CC CAP2=CAP(IK)+SAVV C(IK,IT)=CC SAV(IK,IT)=SAVV DO 303 IT1=1,8 THETA2=THETA(IT)**RHO*EPS(IT1) AMPK2(IT1)=1.0+GNPK(CAP2,THETA2) AMPPK2(IT1)=GNPKK(CAP2,THETA2) CALL PHI(THETA2,NT,PSIT,THMMIN,THMMAX) DO 303 JT=1,NT PSIM2(JT,IT1)=PSIT(JT) 303 CONTINUE CALL PHID(CAP2,NK,PHIC,DPHIC,CAPMIN,CAPMAX) DO 304 JK=1,NK PHIM2(JK)=PHIC(JK) PHIXM2(JK)=DPHIC(JK) 304 CONTINUE VAL=0.0 DO 302 IT1=1,8 C2=0.D0 CP2=0.D0 DO 301 JK=1,NK DO 301 JT=1,NT C2=C2+A(JK,JT)*PHIM2(JK)*PSIM2(JT,IT1) CP2=CP2+A(JK,JT)*PHIXM2(JK)*PSIM2(JT,IT1) 301 CONTINUE VAL=VAL+UTILP(C2)*AMPK2(IT1)*GHERMW(IT1) * We now now compute the derivative of VAL at (IK,IT) with respect to X, * beginning with the sum of terms due to the expected marginal value * of future utility, the RHS of the stochastic Euler equation. UC2 = UTILP(C2) UCC2 = UTILPP(C2) DO 305 JK=1,NK DO 305 JT=1,NT GRAD(JK,JT)=GRAD(JK,JT)+ & (-UC2*AMPPK2(IT1)*PHIM(JK,IK)*PSIM(JT,IT)+ & UCC2*(-PHIM(JK,IK)*PSIM(JT,IT)*CP2 & + PHIM2(JK)*PSIM2(JT,IT1))*AMPK2(IT1))*GHERMW(IT1) 305 CONTINUE 302 CONTINUE C VAL = CC - (BETA*VAL/DSQRT(3.14159D0))**(1.0D0/GAMMA) UCC = UTILPP(CC) DO 51 JK=1,NK DO 51 JT=1,NT GRAD(JK,JT)= PHIM(JK,IK)*PSIM(JT,IT) & -(BETA/DSQRT(3.14159D0))**(1.0D0/GAMMA) & *GRAD(JK,JT)*VAL**(1.0D0/GAMMA-1.0D0)/GAMMA 51 CONTINUE * We now update our running computation of the Jacobian of the projections. DO 60 JK=1,NK DO 60 JT=1,NT J = JK + (JT-1)*NK DO 61 NJK=1,NK DO 61 NJT=1,NT JJAC=NJK+(NJT-1)*NK FJAC(J,JJAC)=FJAC(J,JJAC)+GRAD(NJK,NJT) & *PHIM(JK,IK)*PSIM(JT,IT) 61 CONTINUE 60 CONTINUE 30 CONTINUE C WRITE(*,*) FJAC(1,1), FJAC(5,5), FJAC(7,7) RETURN END ************************************************************************ * FCHK checks the solution at capital stocks and productivity * levels which were not used in computing the solution. This * is an accuracy check. It is the same as F except for the * state variables checked. SUBROUTINE FCHK (N, X, NCHK, MCHK) PARAMETER (NKMX = 10, NTMX = 10) PARAMETER (MXNGRD = 50, MXMGRD = 30, MXN = NKMX, MXM = NTMX) IMPLICIT REAL*8(A-H,O-Z) REAL*8 X(N) REAL*8 SAV(MXNGRD,MXMGRD),A(MXNGRD,MXMGRD) REAL*8 C(MXNGRD,MXMGRD),PROD(MXNGRD,MXMGRD) REAL*8 PHIM(MXN,MXNGRD),PSIM(MXM,MXMGRD) REAL*8 CAP(MXNGRD),EPS(20),THETA(MXMGRD) REAL*8 AMPK2(20),PSIM2(20,MXMGRD), & PHIM2(MXNGRD),C2(20), & GHERMX(20),GHERMW(20) REAL*8 PHIC(MXN),PSIT(MXM),DPHIC(MXN) * COMMON /LENPAR/ NK,NT COMMON /TECH/ AF,ALPHA COMMON /TASTE/ BETA,GAMMA COMMON /BLOCK1/ C,CAP,CAPSS,CAPMIN,CAPMAX COMMON /PHIBLK/ PHIM,PSIM,EPS,RHO,SIG,THETA COMMON /BLOCK2/ SAV,NGRID,MGRID,PROD COMMON /BLOCK3/ CAPBAR,CAPRNG,THBAR,THRNG,THMMIN,THMMAX DATA GHERMX/-2.93063742,-1.98165676,-1.15719371,-.38118699, & .38118699,1.15719371,1.98165676,2.93063742, & 12*0.000/ DATA GHERMW/.00019960,.01707798,.20780233,.66114701, & .66114701,.20780233,.01707798,.00019960, & 12*0.000/ IF(NCHK.EQ.0) RETURN DO 15,I=1,N IT=(I-1)/NK IK=I-IT*NK IT=IT+1 A(IK,IT)=X(I) 15 CONTINUE ERRL2 = 0.0D0 ERRL1 = 0.0D0 ERRLINF = 0.0D0 ERRINR = 0.0D0 ERININF = 0.0D0 NINNER = 0 DO 30,IT=1,MCHK DO 30,IK=1,NCHK CAPCHK = CAPMIN+(IK-1)*(CAPMAX-CAPMIN)/(NCHK-1) THCHK = THMMIN+(IT-1)*(THMMAX-THMMIN)/(MCHK-1) OUTPUT = GNP(CAPCHK,THCHK) CALL PHI(CAPCHK,NK,PHIC,CAPMIN,CAPMAX) CALL PHI(THCHK,NT,PSIT,THMMIN,THMMAX) CC=0. DO 31,LK=1,NK DO 31,LT=1,NT CC=CC+A(LK,LT)*PHIC(LK)*PSIT(LT) 31 CONTINUE SAVV=OUTPUT-CC CAP2=CAPCHK+SAVV DO 303 IT1=1,8 THETA2=THCHK**RHO*EPS(IT1) AMPK2(IT1)=1.0D0+GNPK(CAP2,THETA2) CALL PHI(THETA2,NT,PSIT,THMMIN,THMMAX) DO 303 JT=1,NT PSIM2(JT,IT1)=PSIT(JT) 303 CONTINUE CALL PHI(CAP2,NK,PHIC,CAPMIN,CAPMAX) DO 304 JK=1,NK PHIM2(JK)=PHIC(JK) 304 CONTINUE VAL=0.0D0 DO 302 IT1=1,8 C2(IT1)=0.D0 DO 301 JK=1,NK DO 301 JT=1,NT C2(IT1)=C2(IT1)+A(JK,JT)*PHIM2(JK)*PSIM2(JT,IT1) 301 CONTINUE VAL=VAL+UTILP(C2(IT1))*AMPK2(IT1)*GHERMW(IT1) 302 CONTINUE VAL = CC - (BETA*VAL/DSQRT(3.14159D0))**(1.0D0/GAMMA) C WRITE(12,'(1X,4F12.6,F15.11)') CAPCHK,THCHK,CC,VAL,VAL/CC ABSERR = DABS(VAL/CC) ERRL1 = ERRL1 + ABSERR ERRL2 = ERRL2 + ABSERR*ABSERR THMNIN = .5D0+.5D0*THMMIN THMXIN = .5D0+.5D0*THMMAX IF(CAPCHK.GT.0.75D0.AND.CAPCHK.LT.1.25D0 &.AND.THCHK.GT.THMNIN.AND.THCHK.LT.THMXIN) THEN NINNER = NINNER + 1 ERRINR = ERRINR + ABSERR*ABSERR IF(ABSERR.GT.ERININF) ERININF = ABSERR ENDIF IF(ABSERR.GT.ERRLINF) THEN ERRLINF = ABSERR CAPBIG = CAPCHK THBIG = THCHK ENDIF 30 CONTINUE ERRL1 = ERRL1/(MCHK*NCHK) ERRL2 = (ERRL2/(MCHK*NCHK))**(.5D0) ERRINR = (ERRINR/NINNER)**(.5D0) ERRLINF = DLOG(ERRLINF)/DLOG(10.D0) ERRL2 = DLOG(ERRL2)/DLOG(10.D0) ERRL1 = DLOG(ERRL1)/DLOG(10.D0) ERININF = DLOG(ERININF)/DLOG(10.D0) ERRINR = DLOG(ERRINR)/DLOG(10.D0) C WRITE(12,'(1X,4I6)')NK,NT,NGRID,MGRID WRITE(12,'(1X,8F7.2)')GAMMA,RHO,SIG,ERRLINF,ERRL2,ERRL1, & ERININF,ERRINR WRITE(*,'(1X,4I6)')NK,NT,NGRID,MGRID WRITE(*,'(1X,8F7.2)')GAMMA,RHO,SIG,ERRLINF,ERRL2,ERRL1, & ERININF,ERRINR RETURN END ************************************************************************ SUBROUTINE PHI(ST,N,PHIM,A,B) REAL*8 A,B,ST,PHIM(N) COMMON /BASIS/ IBASIS IF(IBASIS.EQ.1) THEN * PHI is the Chebyshev function XX = -1.0D0 + 2.0D0*(ST-A)/(B-A) PHIM(1) = 1.0D0 PHIM(2) = XX DO 1,I=3,N PHIM(I) = 2.0D0*XX*PHIM(I-1)-PHIM(I-2) 1 CONTINUE RETURN ELSEIF(IBASIS.EQ.2) THEN * PHI is the ordinary polynomial XX = ST PHIM(1) = 1.0D0 PHIM(2) = XX DO 2,I=3,N PHIM(I) = XX*PHIM(I-1) 2 CONTINUE RETURN ENDIF END ************************************************************************ SUBROUTINE PHID(ST,N,PHIM,DPHIM,A,B) REAL*8 A,B,ST,PHIM(N),DPHIM(N) COMMON /BASIS/IBASIS IF(IBASIS.EQ.1) THEN * PHID computes the Chebyshev function and its derivative. XX = -1.0D0 + 2.0D0*(ST-A)/(B-A) PHIM(1) = 1.0D0 DPHIM(1) = 0.0D0 PHIM(2) = XX DPHIM(2) = 1.0D0 DO 1,I=3,N PHIM(I) = 2.0D0*XX*PHIM(I-1)-PHIM(I-2) 1 DPHIM(I) = 2.0D0*PHIM(I-1)+2.0D0*XX*DPHIM(I-1)-DPHIM(I-2) DO 11,I=2,N 11 DPHIM(I) = 2.0D0*DPHIM(I)/(B-A) RETURN ELSEIF(IBASIS.EQ.2) THEN * PHID computes the ordinary polynomial and its derivative. XX = ST PHIM(1) = 1.0D0 DPHIM(1) = 0.0D0 PHIM(2) = XX DPHIM(2) = 1.0D0 DO 2,I=3,N PHIM(I) = XX*PHIM(I-1) 2 DPHIM(I) = PHIM(I-1)+XX*DPHIM(I-1) ENDIF RETURN END ************************************************************************ * Following are function routines for utility, production, and their * derivatives. These functions must be defined for all values of * consumption and capital stock. So we take standard definitions * but take quadratic Taylor expansions around small values for arguments * which are negative and only slightly positive. ************************************************************************ REAL*8 FUNCTION GNP(CAP,TH) REAL*8 CAP,TH,AF,ALPHA COMMON /TECH/ AF,ALPHA IF(CAP.GT.0.001D0) THEN GNP = AF*CAP**ALPHA*TH ELSE GNP = .001D0**ALPHA*AF*TH GNP = GNP + (CAP-.001D0)*AF*ALPHA*0.001D0**(ALPHA-1.0D0)*TH GNP = GNP + ALPHA*(ALPHA-1.0D0)*0.001D0**(ALPHA-2.0D0)*TH & *(CAP-0.001D0)**2/2.0D0 ENDIF RETURN END ************************************************************************ REAL*8 FUNCTION GNPK(CAP,TH) REAL*8 CAP,TH,AF,ALPHA COMMON /TECH/ AF,ALPHA IF(CAP.GT.0.001D0) THEN GNPK = AF*CAP**(ALPHA-1.0D0)*ALPHA*TH ELSE GNPK = AF*ALPHA*0.001D0**(ALPHA-1.0D0)*TH GNPK = GNPK + ALPHA*(ALPHA-1.0D0)*0.001D0**(ALPHA-2.0D0)*TH & *(CAP-0.001D0) ENDIF RETURN END ************************************************************************ REAL*8 FUNCTION UTIL(CONS) REAL*8 CONS,GAMMA,BETA COMMON /TASTE/ BETA,GAMMA IF(CONS.GT.0.01D0) THEN UTIL = CONS**(1.0D0+GAMMA)/(1.0D0+GAMMA) ELSE UTIL = .01D0**(1.0D0+GAMMA)/(1.0D0+GAMMA) UTIL = UTIL + (CONS-.01D0)*0.01D0**GAMMA UTIL = UTIL + GAMMA*0.01D0**(GAMMA-1.0D0) & *(CONS-0.01D0)**2/2.0D0 ENDIF RETURN END ************************************************************************ REAL*8 FUNCTION UTILP(CONS) REAL*8 CONS,GAMMA,BETA COMMON /TASTE/ BETA,GAMMA IF(CONS.GT.0.01D0) THEN UTILP = CONS**GAMMA ELSE UTILP = 0.01D0**GAMMA UTILP = UTILP + GAMMA*0.01D0**(GAMMA-1.0D0) & *(CONS-0.01D0) ENDIF RETURN END ************************************************************************ REAL*8 FUNCTION UTILPP(CONS) REAL*8 CONS,GAMMA,BETA COMMON /TASTE/ BETA,GAMMA IF(CONS.GT.0.01D0) THEN UTILPP = GAMMA*CONS**(-1.0D0+GAMMA) ELSE UTILPP = GAMMA*0.01D0**(GAMMA-1.0D0) ENDIF RETURN END ************************************************************************ REAL*8 FUNCTION GNPKK(CAP,TH) REAL*8 CAP,TH,AF,ALPHA COMMON /TECH/ AF,ALPHA IF(CAP.GT.0.001D0) THEN GNPKK = ALPHA*(ALPHA-1.0D0)*AF*CAP**(ALPHA-2.0D0)*TH ELSE GNPKK = ALPHA*(ALPHA-1.0D0)*0.001D0**(ALPHA-2.0D0)*TH ENDIF RETURN END ************************************************************************ ************************************************************************ C PROGRAM HYBJUD C implicit double precision (a-h,o-z) ************************************************************************ * N is the maximal number of unknowns C PARAMETER (N = 10) ************************************************************************ * HYBRD parameters and arrays C parameter (lwa = (3*n*n+13*n)/2) C double precision x(n),fvec(n),fjac(n,n) C integer wa(lwa) ************************************************************************ * Name of function C external fcn ************************************************************************ * Name input file C open(unit = 8, file = 'hyb.in', status = 'unknown') * Read input C read(8,*) tol C read(8,*) ihot2 * Set initial guess C x(1) = 0.9d0 C x(2) = 0.9d0 C x(3) = 1.9d0 * Ihot = 0 initially C ihot = 0 * Call nonlinear equation solver C call HYBJUD(fcn,3,x,fvec,FJAC,IHOT,tol,info,wa,lwa) C write(*,*) x C write(*,*) info C write(*,*) C x(1) = .95d0*x(1) C x(2) = .99d0*x(2) C x(3) = 1.05d0*x(3) C ihot = ihot2 C call HYBJUD(fcn,3,x,fvec,FJAC,IHOT,tol,info,wa,lwa) C write(*,*) x C write(*,*) info C stop C end C subroutine fcn(n,x,fvec,iflag) C implicit double precision (a-h,o-z) * Dimension arguments C integer n,iflag C double precision x(n),fvec(n) * Declare auxiliary arrays * * Declare common variables and blocks * * Compute fvec C fvec(1) = x(1)**2 + x(2)**2 + x(3)**2 - 14.0d0 C fvec(2) = (x(2)-1.0d0)**3 - 1.0d0 C fvec(3) = x(1) - x(2) + (x(3)-2.0d0)**4 C write(*,'(1x,6f10.6)') x(1), x(2), x(3), C & fvec(1),fvec(2),fvec(3) C return C end subroutine HYBJUD(fcn,n,x,fvec,FJAC,IHOT,tol,info,wa,lwa) integer n,info,lwa implicit double precision (a-h,o-z) double precision tol double precision x(n),fvec(n),wa(lwa),FJAC(N,N) external fcn c ********** c c subroutine HYBJUD c c the purpose of HYBJUD is to find a zero of a system of c n nonlinear functions in n variables by a modification c of the powell hybrid method. this is done by using the c more general nonlinear equation solver hybrd. the user c must provide a subroutine which calculates the functions. c the jacobian is then calculated by a forward-difference c approximation. c c the subroutine statement is c c subroutine HYBJUD(fcn,n,x,fvec,tol,info,wa,lwa) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(n,x,fvec,iflag) c integer n,iflag c double precision x(n),fvec(n) c ---------- c calculate the functions at x and c return this vector in fvec. c --------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of HYBJUD. c in this case set iflag to a negative integer. c c n is a positive integer input variable set to the number c of functions and variables. c c x is an array of length n. on input x must contain c an initial estimate of the solution vector. on output x c contains the final estimate of the solution vector. c c fvec is an output array of length n which contains c the functions evaluated at the output x. c c tol is a nonnegative input variable. termination occurs c when the algorithm estimates that the relative error c between x and the solution is at most tol. c c info is an integer output variable. if the user has c terminated execution, info is set to the (negative) c value of iflag. see description of fcn. otherwise, c info is set as follows. c c info = 0 improper input parameters. c c info = 1 algorithm estimates that the relative error c between x and the solution is at most tol. c c info = 2 number of calls to fcn has reached or exceeded c 200*(n+1). c c info = 3 tol is too small. no further improvement in c the approximate solution x is possible. c c info = 4 iteration is not making good progress. c c wa is a work array of length lwa. c c lwa is a positive integer input variable not less than c (n*(3*n+13))/2. c c subprograms called c c user-supplied ...... fcn c c minpack-supplied ... hybrd c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer index,j,lr,maxfev,ml,mode,mu,nfev,nprint double precision epsfcn,factor,one,xtol,zero data factor,one,zero /1.0d2,1.0d0,0.0d0/ info = 0 c c check the input parameters for errors. c if (n .le. 0 .or. tol .lt. zero .or. lwa .lt. (n*(3*n + 13))/2) * go to 20 c c call hybrd. c maxfev = 200*(n + 1) xtol = tol ml = n - 1 mu = n - 1 epsfcn = zero mode = 2 do 10 j = 1, n wa(j) = one 10 continue nprint = 0 lr = (n*(n + 1))/2 index = 6*n + lr call hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,wa(1),mode, * factor,nprint,info,nfev,FJAC,n,wa(6*n+1),lr, * wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1),IHOT) if (info .eq. 5) info = 4 IJ = 0 DO I=1,N DO J=1,N IJ = IJ + 1 SUM = 0.0D0 JJ = J DO L=1,J SUM = SUM + FJAC(I,L)*WA(6*N+1+JJ-1) JJ = JJ + N - L ENDDO WA(INDEX+IJ) = SUM ENDDO ENDDO IJ = 0 DO I = 1,N DO J = 1,N IJ = IJ+1 FJAC(I,J) = WA(INDEX+IJ) ENDDO ENDDO 20 continue return c c last card of subroutine HYBJUD. c end subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag, * mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr, * qtf,wa1,wa2,wa3,wa4,IHOT) integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr implicit double precision (a-h,o-z) double precision xtol,epsfcn,factor double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr), * qtf(n),wa1(n),wa2(n),wa3(n),wa4(n) external fcn c ********** c c subroutine hybrd c c the purpose of hybrd is to find a zero of a system of c n nonlinear functions in n variables by a modification c of the powell hybrid method. the user must provide a c subroutine which calculates the functions. the jacobian is c then calculated by a forward-difference approximation. c c the subroutine statement is c c subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn, c diag,mode,factor,nprint,info,nfev,fjac, c ldfjac,r,lr,qtf,wa1,wa2,wa3,wa4) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(n,x,fvec,iflag) c integer n,iflag c double precision x(n),fvec(n) c ---------- c calculate the functions at x and c return this vector in fvec. c --------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of hybrd. c in this case set iflag to a negative integer. c c n is a positive integer input variable set to the number c of functions and variables. c c x is an array of length n. on input x must contain c an initial estimate of the solution vector. on output x c contains the final estimate of the solution vector. c c fvec is an output array of length n which contains c the functions evaluated at the output x. c c xtol is a nonnegative input variable. termination c occurs when the relative error between two consecutive c iterates is at most xtol. c c maxfev is a positive integer input variable. termination c occurs when the number of calls to fcn is at least maxfev c by the end of an iteration. c c ml is a nonnegative integer input variable which specifies c the number of subdiagonals within the band of the c jacobian matrix. if the jacobian is not banded, set c ml to at least n - 1. c c mu is a nonnegative integer input variable which specifies c the number of superdiagonals within the band of the c jacobian matrix. if the jacobian is not banded, set c mu to at least n - 1. c c epsfcn is an input variable used in determining a suitable c step length for the forward-difference approximation. this c approximation assumes that the relative errors in the c functions are of the order of epsfcn. if epsfcn is less c than the machine precision, it is assumed that the relative c errors in the functions are of the order of the machine c precision. c c diag is an array of length n. if mode = 1 (see c below), diag is internally set. if mode = 2, diag c must contain positive entries that serve as c multiplicative scale factors for the variables. c c mode is an integer input variable. if mode = 1, the c variables will be scaled internally. if mode = 2, c the scaling is specified by the input diag. other c values of mode are equivalent to mode = 1. c c factor is a positive input variable used in determining the c initial step bound. this bound is set to the product of c factor and the euclidean norm of diag*x if nonzero, or else c to factor itself. in most cases factor should lie in the c interval (.1,100.). 100. is a generally recommended value. c c nprint is an integer input variable that enables controlled c printing of iterates if it is positive. in this case, c fcn is called with iflag = 0 at the beginning of the first c iteration and every nprint iterations thereafter and c immediately prior to return, with x and fvec available c for printing. if nprint is not positive, no special calls c of fcn with iflag = 0 are made. c c info is an integer output variable. if the user has c terminated execution, info is set to the (negative) c value of iflag. see description of fcn. otherwise, c info is set as follows. c c info = 0 improper input parameters. c c info = 1 relative error between two consecutive iterates c is at most xtol. c c info = 2 number of calls to fcn has reached or exceeded c maxfev. c c info = 3 xtol is too small. no further improvement in c the approximate solution x is possible. c c info = 4 iteration is not making good progress, as c measured by the improvement from the last c five jacobian evaluations. c c info = 5 iteration is not making good progress, as c measured by the improvement from the last c ten iterations. c c nfev is an integer output variable set to the number of c calls to fcn. c c fjac is an output n by n array which contains the c orthogonal matrix q produced by the qr factorization c of the final approximate jacobian. c c ldfjac is a positive integer input variable not less than n c which specifies the leading dimension of the array fjac. c c r is an output array of length lr which contains the c upper triangular matrix produced by the qr factorization c of the final approximate jacobian, stored rowwise. c c lr is a positive integer input variable not less than c (n*(n+1))/2. c c qtf is an output array of length n which contains c the vector (q transpose)*fvec. c c wa1, wa2, wa3, and wa4 are work arrays of length n. c c subprograms called c c user-supplied ...... fcn c c minpack-supplied ... dogleg,dpmpar,enorm,fdjac1, c qform,qrfac,r1mpyq,r1updt c c fortran-supplied ... dabs,dmax1,dmin1,min0,mod c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,iflag,iter,j,jm1,l,msum,ncfail,ncsuc,nslow1,nslow2 integer iwa(1) logical jeval,sing double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm, * prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm, * zero double precision dpmpar,enorm data one,p1,p5,p001,p0001,zero * /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c info = 0 iflag = 0 nfev = 0 c c check the input parameters for errors. c if (n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0 * .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero * .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) go to 300 if (mode .ne. 2) go to 20 do 10 j = 1, n if (diag(j) .le. zero) go to 300 10 continue 20 continue c c evaluate the function at the starting point c and calculate its norm. c iflag = 1 call fcn(n,x,fvec,iflag) nfev = 1 if (iflag .lt. 0) go to 300 fnorm = enorm(n,fvec) c c determine the number of calls to fcn needed to compute c the jacobian matrix. c msum = min0(ml+mu+1,n) c c initialize iteration counter and monitors. c iter = 1 ncsuc = 0 ncfail = 0 nslow1 = 0 nslow2 = 0 c c beginning of the outer loop. c 30 continue jeval = .true. c c calculate the jacobian matrix. c iflag = 2 C C IF IHOT=1, THEN WE USE EXISTING FJAC. OTHERWISE, COMPUTE NEW ONE. C IF(IHOT.EQ.0) THEN call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1, * wa2) nfev = nfev + msum if (iflag .lt. 0) go to 300 ELSE IHOT = 0 ENDIF c c compute the qr factorization of the jacobian. c call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3) c c on the first iteration and if mode is 1, scale according c to the norms of the columns of the initial jacobian. c if (iter .ne. 1) go to 70 if (mode .eq. 2) go to 50 do 40 j = 1, n diag(j) = wa2(j) if (wa2(j) .eq. zero) diag(j) = one 40 continue 50 continue c c on the first iteration, calculate the norm of the scaled x c and initialize the step bound delta. c do 60 j = 1, n wa3(j) = diag(j)*x(j) 60 continue xnorm = enorm(n,wa3) delta = factor*xnorm if (delta .eq. zero) delta = factor 70 continue c c form (q transpose)*fvec and store in qtf. c do 80 i = 1, n qtf(i) = fvec(i) 80 continue do 120 j = 1, n if (fjac(j,j) .eq. zero) go to 110 sum = zero do 90 i = j, n sum = sum + fjac(i,j)*qtf(i) 90 continue temp = -sum/fjac(j,j) do 100 i = j, n qtf(i) = qtf(i) + fjac(i,j)*temp 100 continue 110 continue 120 continue c c copy the triangular factor of the qr factorization into r. c sing = .false. do 150 j = 1, n l = j jm1 = j - 1 if (jm1 .lt. 1) go to 140 do 130 i = 1, jm1 r(l) = fjac(i,j) l = l + n - i 130 continue 140 continue r(l) = wa1(j) if (wa1(j) .eq. zero) sing = .true. 150 continue c c accumulate the orthogonal factor in fjac. c call qform(n,n,fjac,ldfjac,wa1) c c rescale if necessary. c if (mode .eq. 2) go to 170 do 160 j = 1, n diag(j) = dmax1(diag(j),wa2(j)) 160 continue 170 continue c c beginning of the inner loop. c 180 continue c c if requested, call fcn to enable printing of iterates. c if (nprint .le. 0) go to 190 iflag = 0 if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag) if (iflag .lt. 0) go to 300 190 continue c c determine the direction p. c call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3) c c store the direction p and x + p. calculate the norm of p. c do 200 j = 1, n wa1(j) = -wa1(j) wa2(j) = x(j) + wa1(j) wa3(j) = diag(j)*wa1(j) 200 continue pnorm = enorm(n,wa3) c c on the first iteration, adjust the initial step bound. c if (iter .eq. 1) delta = dmin1(delta,pnorm) c c evaluate the function at x + p and calculate its norm. c iflag = 1 call fcn(n,wa2,wa4,iflag) nfev = nfev + 1 if (iflag .lt. 0) go to 300 fnorm1 = enorm(n,wa4) c c compute the scaled actual reduction. c actred = -one if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 c c compute the scaled predicted reduction. c l = 1 do 220 i = 1, n sum = zero do 210 j = i, n sum = sum + r(l)*wa1(j) l = l + 1 210 continue wa3(i) = qtf(i) + sum 220 continue temp = enorm(n,wa3) prered = zero if (temp .lt. fnorm) prered = one - (temp/fnorm)**2 c c compute the ratio of the actual to the predicted c reduction. c ratio = zero if (prered .gt. zero) ratio = actred/prered c c update the step bound. c if (ratio .ge. p1) go to 230 ncsuc = 0 ncfail = ncfail + 1 delta = p5*delta go to 240 230 continue ncfail = 0 ncsuc = ncsuc + 1 if (ratio .ge. p5 .or. ncsuc .gt. 1) * delta = dmax1(delta,pnorm/p5) if (dabs(ratio-one) .le. p1) delta = pnorm/p5 240 continue c c test for successful iteration. c if (ratio .lt. p0001) go to 260 c c successful iteration. update x, fvec, and their norms. c do 250 j = 1, n x(j) = wa2(j) wa2(j) = diag(j)*x(j) fvec(j) = wa4(j) 250 continue xnorm = enorm(n,wa2) fnorm = fnorm1 iter = iter + 1 260 continue c c determine the progress of the iteration. c nslow1 = nslow1 + 1 if (actred .ge. p001) nslow1 = 0 if (jeval) nslow2 = nslow2 + 1 if (actred .ge. p1) nslow2 = 0 c c test for convergence. c if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1 if (info .ne. 0) go to 300 c c tests for termination and stringent tolerances. c if (nfev .ge. maxfev) info = 2 if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3 if (nslow2 .eq. 5) info = 4 if (nslow1 .eq. 10) info = 5 if (info .ne. 0) go to 300 c c criterion for recalculating jacobian approximation c by forward differences. c if (ncfail .eq. 2) go to 290 c c calculate the rank one modification to the jacobian c and update qtf if necessary. c do 280 j = 1, n sum = zero do 270 i = 1, n sum = sum + fjac(i,j)*wa4(i) 270 continue wa2(j) = (sum - wa3(j))/pnorm wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm) if (ratio .ge. p0001) qtf(j) = sum 280 continue c c compute the qr factorization of the updated jacobian. c call r1updt(n,n,r,lr,wa1,wa2,wa3,sing) call r1mpyq(n,n,fjac,ldfjac,wa2,wa3) call r1mpyq(1,n,qtf,1,wa2,wa3) c c end of the inner loop. c jeval = .false. go to 180 290 continue c c end of the outer loop. c go to 30 300 continue c c termination, either normal or user imposed. c if (iflag .lt. 0) info = iflag iflag = 0 if (nprint .gt. 0) call fcn(n,x,fvec,iflag) return c c last card of subroutine hybrd. c end subroutine qform(m,n,q,ldq,wa) integer m,n,ldq implicit double precision (a-h,o-z) double precision q(ldq,m),wa(m) c ********** c c subroutine qform c c this subroutine proceeds from the computed qr factorization of c an m by n matrix a to accumulate the m by m orthogonal matrix c q from its factored form. c c the subroutine statement is c c subroutine qform(m,n,q,ldq,wa) c c where c c m is a positive integer input variable set to the number c of rows of a and the order of q. c c n is a positive integer input variable set to the number c of columns of a. c c q is an m by m array. on input the full lower trapezoid in c the first min(m,n) columns of q contains the factored form. c on output q has been accumulated into a square matrix. c c ldq is a positive integer input variable not less than m c which specifies the leading dimension of the array q. c c wa is a work array of length m. c c subprograms called c c fortran-supplied ... min0 c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,jm1,k,l,minmn,np1 double precision one,sum,temp,zero data one,zero /1.0d0,0.0d0/ c c zero out upper triangle of q in the first min(m,n) columns. c minmn = min0(m,n) if (minmn .lt. 2) go to 30 do 20 j = 2, minmn jm1 = j - 1 do 10 i = 1, jm1 q(i,j) = zero 10 continue 20 continue 30 continue c c initialize remaining columns to those of the identity matrix. c np1 = n + 1 if (m .lt. np1) go to 60 do 50 j = np1, m do 40 i = 1, m q(i,j) = zero 40 continue q(j,j) = one 50 continue 60 continue c c accumulate q from its factored form. c do 120 l = 1, minmn k = minmn - l + 1 do 70 i = k, m wa(i) = q(i,k) q(i,k) = zero 70 continue q(k,k) = one if (wa(k) .eq. zero) go to 110 do 100 j = k, m sum = zero do 80 i = k, m sum = sum + q(i,j)*wa(i) 80 continue temp = sum/wa(k) do 90 i = k, m q(i,j) = q(i,j) - temp*wa(i) 90 continue 100 continue 110 continue 120 continue return c c last card of subroutine qform. c end subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) integer m,n,lda,lipvt integer ipvt(lipvt) logical pivot implicit double precision (a-h,o-z) double precision a(lda,n),rdiag(n),acnorm(n),wa(n) c ********** c c subroutine qrfac c c this subroutine uses householder transformations with column c pivoting (optional) to compute a qr factorization of the c m by n matrix a. that is, qrfac determines an orthogonal c matrix q, a permutation matrix p, and an upper trapezoidal c matrix r with diagonal elements of nonincreasing magnitude, c such that a*p = q*r. the householder transformation for c column k, k = 1,2,...,min(m,n), is of the form c c t c i - (1/u(k))*u*u c c where u has zeros in the first k-1 positions. the form of c this transformation and the method of pivoting first c appeared in the corresponding linpack subroutine. c c the subroutine statement is c c subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) c c where c c m is a positive integer input variable set to the number c of rows of a. c c n is a positive integer input variable set to the number c of columns of a. c c a is an m by n array. on input a contains the matrix for c which the qr factorization is to be computed. on output c the strict upper trapezoidal part of a contains the strict c upper trapezoidal part of r, and the lower trapezoidal c part of a contains a factored form of q (the non-trivial c elements of the u vectors described above). c c lda is a positive integer input variable not less than m c which specifies the leading dimension of the array a. c c pivot is a logical input variable. if pivot is set true, c then column pivoting is enforced. if pivot is set false, c then no column pivoting is done. c c ipvt is an integer output array of length lipvt. ipvt c defines the permutation matrix p such that a*p = q*r. c column j of p is column ipvt(j) of the identity matrix. c if pivot is false, ipvt is not referenced. c c lipvt is a positive integer input variable. if pivot is false, c then lipvt may be as small as 1. if pivot is true, then c lipvt must be at least n. c c rdiag is an output array of length n which contains the c diagonal elements of r. c c acnorm is an output array of length n which contains the c norms of the corresponding columns of the input matrix a. c if this information is not needed, then acnorm can coincide c with rdiag. c c wa is a work array of length n. if pivot is false, then wa c can coincide with rdiag. c c subprograms called c c minpack-supplied ... dpmpar,enorm c c fortran-supplied ... dmax1,dsqrt,min0 c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,jp1,k,kmax,minmn double precision ajnorm,epsmch,one,p05,sum,temp,zero double precision dpmpar,enorm data one,p05,zero /1.0d0,5.0d-2,0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c c compute the initial column norms and initialize several arrays. c do 10 j = 1, n acnorm(j) = enorm(m,a(1,j)) rdiag(j) = acnorm(j) wa(j) = rdiag(j) if (pivot) ipvt(j) = j 10 continue c c reduce a to r with householder transformations. c minmn = min0(m,n) do 110 j = 1, minmn if (.not.pivot) go to 40 c c bring the column of largest norm into the pivot position. c kmax = j do 20 k = j, n if (rdiag(k) .gt. rdiag(kmax)) kmax = k 20 continue if (kmax .eq. j) go to 40 do 30 i = 1, m temp = a(i,j) a(i,j) = a(i,kmax) a(i,kmax) = temp 30 continue rdiag(kmax) = rdiag(j) wa(kmax) = wa(j) k = ipvt(j) ipvt(j) = ipvt(kmax) ipvt(kmax) = k 40 continue c c compute the householder transformation to reduce the c j-th column of a to a multiple of the j-th unit vector. c ajnorm = enorm(m-j+1,a(j,j)) if (ajnorm .eq. zero) go to 100 if (a(j,j) .lt. zero) ajnorm = -ajnorm do 50 i = j, m a(i,j) = a(i,j)/ajnorm 50 continue a(j,j) = a(j,j) + one c c apply the transformation to the remaining columns c and update the norms. c jp1 = j + 1 if (n .lt. jp1) go to 100 do 90 k = jp1, n sum = zero do 60 i = j, m sum = sum + a(i,j)*a(i,k) 60 continue temp = sum/a(j,j) do 70 i = j, m a(i,k) = a(i,k) - temp*a(i,j) 70 continue if (.not.pivot .or. rdiag(k) .eq. zero) go to 80 temp = a(j,k)/rdiag(k) rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2)) if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80 rdiag(k) = enorm(m-j,a(jp1,k)) wa(k) = rdiag(k) 80 continue 90 continue 100 continue rdiag(j) = -ajnorm 110 continue return c c last card of subroutine qrfac. c end subroutine r1mpyq(m,n,a,lda,v,w) integer m,n,lda implicit double precision (a-h,o-z) double precision a(lda,n),v(n),w(n) c ********** c c subroutine r1mpyq c c given an m by n matrix a, this subroutine computes a*q where c q is the product of 2*(n - 1) transformations c c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1) c c and gv(i), gw(i) are givens rotations in the (i,n) plane which c eliminate elements in the i-th and n-th planes, respectively. c q itself is not given, rather the information to recover the c gv, gw rotations is supplied. c c the subroutine statement is c c subroutine r1mpyq(m,n,a,lda,v,w) c c where c c m is a positive integer input variable set to the number c of rows of a. c c n is a positive integer input variable set to the number c of columns of a. c c a is an m by n array. on input a must contain the matrix c to be postmultiplied by the orthogonal matrix q c described above. on output a*q has replaced a. c c lda is a positive integer input variable not less than m c which specifies the leading dimension of the array a. c c v is an input array of length n. v(i) must contain the c information necessary to recover the givens rotation gv(i) c described above. c c w is an input array of length n. w(i) must contain the c information necessary to recover the givens rotation gw(i) c described above. c c subroutines called c c fortran-supplied ... dabs,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,nmj,nm1 double precision cos,one,sin,temp data one /1.0d0/ c c apply the first set of givens rotations to a. c nm1 = n - 1 if (nm1 .lt. 1) go to 50 do 20 nmj = 1, nm1 j = n - nmj if (dabs(v(j)) .gt. one) cos = one/v(j) if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2) if (dabs(v(j)) .le. one) sin = v(j) if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2) do 10 i = 1, m temp = cos*a(i,j) - sin*a(i,n) a(i,n) = sin*a(i,j) + cos*a(i,n) a(i,j) = temp 10 continue 20 continue c c apply the second set of givens rotations to a. c do 40 j = 1, nm1 if (dabs(w(j)) .gt. one) cos = one/w(j) if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2) if (dabs(w(j)) .le. one) sin = w(j) if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2) do 30 i = 1, m temp = cos*a(i,j) + sin*a(i,n) a(i,n) = -sin*a(i,j) + cos*a(i,n) a(i,j) = temp 30 continue 40 continue 50 continue return c c last card of subroutine r1mpyq. c end subroutine r1updt(m,n,s,ls,u,v,w,sing) integer m,n,ls logical sing implicit double precision (a-h,o-z) double precision s(ls),u(m),v(n),w(m) c ********** c c subroutine r1updt c c given an m by n lower trapezoidal matrix s, an m-vector u, c and an n-vector v, the problem is to determine an c orthogonal matrix q such that c c t c (s + u*v )*q c c is again lower trapezoidal. c c this subroutine determines q as the product of 2*(n - 1) c transformations c c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1) c c where gv(i), gw(i) are givens rotations in the (i,n) plane c which eliminate elements in the i-th and n-th planes, c respectively. q itself is not accumulated, rather the c information to recover the gv, gw rotations is returned. c c the subroutine statement is c c subroutine r1updt(m,n,s,ls,u,v,w,sing) c c where c c m is a positive integer input variable set to the number c of rows of s. c c n is a positive integer input variable set to the number c of columns of s. n must not exceed m. c c s is an array of length ls. on input s must contain the lower c trapezoidal matrix s stored by columns. on output s contains c the lower trapezoidal matrix produced as described above. c c ls is a positive integer input variable not less than c (n*(2*m-n+1))/2. c c u is an input array of length m which must contain the c vector u. c c v is an array of length n. on input v must contain the vector c v. on output v(i) contains the information necessary to c recover the givens rotation gv(i) described above. c c w is an output array of length m. w(i) contains information c necessary to recover the givens rotation gw(i) described c above. c c sing is a logical output variable. sing is set true if any c of the diagonal elements of the output s are zero. otherwise c sing is set false. c c subprograms called c c minpack-supplied ... dpmpar c c fortran-supplied ... dabs,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more, c john l. nazareth c c ********** integer i,j,jj,l,nmj,nm1 double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp, * zero double precision dpmpar data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/ c c giant is the largest magnitude. c giant = dpmpar(3) c c initialize the diagonal element pointer. c jj = (n*(2*m - n + 1))/2 - (m - n) c c move the nontrivial part of the last column of s into w. c l = jj do 10 i = n, m w(i) = s(l) l = l + 1 10 continue c c rotate the vector v into a multiple of the n-th unit vector c in such a way that a spike is introduced into w. c nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 nmj = 1, nm1 j = n - nmj jj = jj - (m - j + 1) w(j) = zero if (v(j) .eq. zero) go to 50 c c determine a givens rotation which eliminates the c j-th element of v. c if (dabs(v(n)) .ge. dabs(v(j))) go to 20 cotan = v(n)/v(j) sin = p5/dsqrt(p25+p25*cotan**2) cos = sin*cotan tau = one if (dabs(cos)*giant .gt. one) tau = one/cos go to 30 20 continue tan = v(j)/v(n) cos = p5/dsqrt(p25+p25*tan**2) sin = cos*tan tau = sin 30 continue c c apply the transformation to v and store the information c necessary to recover the givens rotation. c v(n) = sin*v(j) + cos*v(n) v(j) = tau c c apply the transformation to s and extend the spike in w. c l = jj do 40 i = j, m temp = cos*s(l) - sin*w(i) w(i) = sin*s(l) + cos*w(i) s(l) = temp l = l + 1 40 continue 50 continue 60 continue 70 continue c c add the spike from the rank 1 update to w. c do 80 i = 1, m w(i) = w(i) + v(n)*u(i) 80 continue c c eliminate the spike. c sing = .false. if (nm1 .lt. 1) go to 140 do 130 j = 1, nm1 if (w(j) .eq. zero) go to 120 c c determine a givens rotation which eliminates the c j-th element of the spike. c if (dabs(s(jj)) .ge. dabs(w(j))) go to 90 cotan = s(jj)/w(j) sin = p5/dsqrt(p25+p25*cotan**2) cos = sin*cotan tau = one if (dabs(cos)*giant .gt. one) tau = one/cos go to 100 90 continue tan = w(j)/s(jj) cos = p5/dsqrt(p25+p25*tan**2) sin = cos*tan tau = sin 100 continue c c apply the transformation to s and reduce the spike in w. c l = jj do 110 i = j, m temp = cos*s(l) + sin*w(i) w(i) = -sin*s(l) + cos*w(i) s(l) = temp l = l + 1 110 continue c c store the information necessary to recover the c givens rotation. c w(j) = tau 120 continue c c test for zero diagonal elements in the output s. c if (s(jj) .eq. zero) sing = .true. jj = jj + (m - j + 1) 130 continue 140 continue c c move w back into the last column of the output s. c l = jj do 150 i = n, m s(l) = w(i) l = l + 1 150 continue if (s(jj) .eq. zero) sing = .true. return c c last card of subroutine r1updt. c end subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2) integer n,lr implicit double precision (a-h,o-z) double precision delta double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n) c ********** c c subroutine dogleg c c given an m by n matrix a, an n by n nonsingular diagonal c matrix d, an m-vector b, and a positive number delta, the c problem is to determine the convex combination x of the c gauss-newton and scaled gradient directions that minimizes c (a*x - b) in the least squares sense, subject to the c restriction that the euclidean norm of d*x be at most delta. c c this subroutine completes the solution of the problem c if it is provided with the necessary information from the c qr factorization of a. that is, if a = q*r, where q has c orthogonal columns and r is an upper triangular matrix, c then dogleg expects the full upper triangle of r and c the first n components of (q transpose)*b. c c the subroutine statement is c c subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2) c c where c c n is a positive integer input variable set to the order of r. c c r is an input array of length lr which must contain the upper c triangular matrix r stored by rows. c c lr is a positive integer input variable not less than c (n*(n+1))/2. c c diag is an input array of length n which must contain the c diagonal elements of the matrix d. c c qtb is an input array of length n which must contain the first c n elements of the vector (q transpose)*b. c c delta is a positive input variable which specifies an upper c bound on the euclidean norm of d*x. c c x is an output array of length n which contains the desired c convex combination of the gauss-newton direction and the c scaled gradient direction. c c wa1 and wa2 are work arrays of length n. c c subprograms called c c minpack-supplied ... dpmpar,enorm c c fortran-supplied ... dabs,dmax1,dmin1,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,jj,jp1,k,l double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum, * temp,zero double precision dpmpar,enorm data one,zero /1.0d0,0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c c first, calculate the gauss-newton direction. c jj = (n*(n + 1))/2 + 1 do 50 k = 1, n j = n - k + 1 jp1 = j + 1 jj = jj - k l = jj + 1 sum = zero if (n .lt. jp1) go to 20 do 10 i = jp1, n sum = sum + r(l)*x(i) l = l + 1 10 continue 20 continue temp = r(jj) if (temp .ne. zero) go to 40 l = j do 30 i = 1, j temp = dmax1(temp,dabs(r(l))) l = l + n - i 30 continue temp = epsmch*temp if (temp .eq. zero) temp = epsmch 40 continue x(j) = (qtb(j) - sum)/temp 50 continue c c test whether the gauss-newton direction is acceptable. c do 60 j = 1, n wa1(j) = zero wa2(j) = diag(j)*x(j) 60 continue qnorm = enorm(n,wa2) if (qnorm .le. delta) go to 140 c c the gauss-newton direction is not acceptable. c next, calculate the scaled gradient direction. c l = 1 do 80 j = 1, n temp = qtb(j) do 70 i = j, n wa1(i) = wa1(i) + r(l)*temp l = l + 1 70 continue wa1(j) = wa1(j)/diag(j) 80 continue c c calculate the norm of the scaled gradient and test for c the special case in which the scaled gradient is zero. c gnorm = enorm(n,wa1) sgnorm = zero alpha = delta/qnorm if (gnorm .eq. zero) go to 120 c c calculate the point along the scaled gradient c at which the quadratic is minimized. c do 90 j = 1, n wa1(j) = (wa1(j)/gnorm)/diag(j) 90 continue l = 1 do 110 j = 1, n sum = zero do 100 i = j, n sum = sum + r(l)*wa1(i) l = l + 1 100 continue wa2(j) = sum 110 continue temp = enorm(n,wa2) sgnorm = (gnorm/temp)/temp c c test whether the scaled gradient direction is acceptable. c alpha = zero if (sgnorm .ge. delta) go to 120 c c the scaled gradient direction is not acceptable. c finally, calculate the point along the dogleg c at which the quadratic is minimized. c bnorm = enorm(n,qtb) temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta) temp = temp - (delta/qnorm)*(sgnorm/delta)**2 * + dsqrt((temp-(delta/qnorm))**2 * +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2)) alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp 120 continue c c form appropriate convex combination of the gauss-newton c direction and the scaled gradient direction. c temp = (one - alpha)*dmin1(sgnorm,delta) do 130 j = 1, n x(j) = temp*wa1(j) + alpha*x(j) 130 continue 140 continue return c c last card of subroutine dogleg. c end double precision function dpmpar(i) integer i c ********** c c function dpmpar c c This function provides double precision machine parameters c when the appropriate set of data statements is activated (by c removing the c from column 1) and all other data statements are c rendered inactive. Most of the parameter values were obtained c from the corresponding Bell Laboratories Port Library function. c c The function statement is c c double precision function dpmpar(i) c c where c c i is an integer input variable set to 1, 2, or 3 which c selects the desired machine parameter. If the machine has c t base b digits and its smallest and largest exponents are c emin and emax, respectively, then these parameters are c c dpmpar(1) = b**(1 - t), the machine precision, c c dpmpar(2) = b**(emin - 1), the smallest magnitude, c c dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude. c c Argonne National Laboratory. MINPACK Project. June 1983. c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More c c ********** integer mcheps(4) integer minmag(4) integer maxmag(4) double precision dmach(3) equivalence (dmach(1),mcheps(1)) equivalence (dmach(2),minmag(1)) equivalence (dmach(3),maxmag(1)) c c Machine constants for PC's. c data dmach(1)/0.00000000001d0/ data dmach(2)/1.0d-300/ data dmach(3)/1.0d300/ c data mcheps(1),mcheps(2) / 9472, 0 / c data minmag(1),minmag(2) / 128, 0 / c data maxmag(1),maxmag(2) / -32769, -1 / c dpmpar = dmach(i) return c c Last card of function dpmpar. c end double precision function enorm(n,x) integer n double precision x(n) c ********** c c function enorm c c given an n-vector x, this function calculates the c euclidean norm of x. c c the euclidean norm is computed by accumulating the sum of c squares in three different sums. the sums of squares for the c small and large components are scaled so that no overflows c occur. non-destructive underflows are permitted. underflows c and overflows do not occur in the computation of the unscaled c sum of squares for the intermediate components. c the definitions of small, intermediate and large components c depend on two constants, rdwarf and rgiant. the main c restrictions on these constants are that rdwarf**2 not c underflow and rgiant**2 not overflow. the constants c given here are suitable for every known computer. c c the function statement is c c double precision function enorm(n,x) c c where c c n is a positive integer input variable. c c x is an input array of length n. c c subprograms called c c fortran-supplied ... dabs,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs, * x1max,x3max,zero data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/ s1 = zero s2 = zero s3 = zero x1max = zero x3max = zero floatn = n agiant = rgiant/floatn do 90 i = 1, n xabs = dabs(x(i)) if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70 if (xabs .le. rdwarf) go to 30 c c sum for large components. c if (xabs .le. x1max) go to 10 s1 = one + s1*(x1max/xabs)**2 x1max = xabs go to 20 10 continue s1 = s1 + (xabs/x1max)**2 20 continue go to 60 30 continue c c sum for small components. c if (xabs .le. x3max) go to 40 s3 = one + s3*(x3max/xabs)**2 x3max = xabs go to 50 40 continue if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2 50 continue 60 continue go to 80 70 continue c c sum for intermediate components. c s2 = s2 + xabs**2 80 continue 90 continue c c calculation of norm. c if (s1 .eq. zero) go to 100 enorm = x1max*dsqrt(s1+(s2/x1max)/x1max) go to 130 100 continue if (s2 .eq. zero) go to 110 if (s2 .ge. x3max) * enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3))) if (s2 .lt. x3max) * enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3))) go to 120 110 continue enorm = x3max*dsqrt(s3) 120 continue 130 continue return c c last card of function enorm. c end subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn, * wa1,wa2) integer n,ldfjac,iflag,ml,mu implicit double precision (a-h,o-z) double precision epsfcn double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n) c ********** c c subroutine fdjac1 c c this subroutine computes a forward-difference approximation c to the n by n jacobian matrix associated with a specified c problem of n functions in n variables. if the jacobian has c a banded form, then function evaluations are saved by only c approximating the nonzero terms. c c the subroutine statement is c c subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn, c wa1,wa2) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(n,x,fvec,iflag) c integer n,iflag c double precision x(n),fvec(n) c ---------- c calculate the functions at x and c return this vector in fvec. c ---------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of fdjac1. c in this case set iflag to a negative integer. c c n is a positive integer input variable set to the number c of functions and variables. c c x is an input array of length n. c c fvec is an input array of length n which must contain the c functions evaluated at x. c c fjac is an output n by n array which contains the c approximation to the jacobian matrix evaluated at x. c c ldfjac is a positive integer input variable not less than n c which specifies the leading dimension of the array fjac. c c iflag is an integer variable which can be used to terminate c the execution of fdjac1. see description of fcn. c c ml is a nonnegative integer input variable which specifies c the number of subdiagonals within the band of the c jacobian matrix. if the jacobian is not banded, set c ml to at least n - 1. c c epsfcn is an input variable used in determining a suitable c step length for the forward-difference approximation. this c approximation assumes that the relative errors in the c functions are of the order of epsfcn. if epsfcn is less c than the machine precision, it is assumed that the relative c errors in the functions are of the order of the machine c precision. c c mu is a nonnegative integer input variable which specifies c the number of superdiagonals within the band of the c jacobian matrix. if the jacobian is not banded, set c mu to at least n - 1. c c wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at c least n, then the jacobian is considered dense, and wa2 is c not referenced. c c subprograms called c c minpack-supplied ... dpmpar c c fortran-supplied ... dabs,dmax1,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,k,msum double precision eps,epsmch,h,temp,zero double precision dpmpar data zero /0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c eps = dsqrt(dmax1(epsfcn,epsmch)) msum = ml + mu + 1 if (msum .lt. n) go to 40 c c computation of dense approximate jacobian. c do 20 j = 1, n temp = x(j) h = eps*dabs(temp) C C I have changed the difference formula to avoid h becoming too small. C C if (h .eq. zero) h = eps IF (DABS(H).LT.EPS) H = EPS x(j) = temp + h call fcn(n,x,wa1,iflag) if (iflag .lt. 0) go to 30 x(j) = temp do 10 i = 1, n fjac(i,j) = (wa1(i) - fvec(i))/h 10 continue 20 continue 30 continue go to 110 40 continue c c computation of banded approximate jacobian. c do 90 k = 1, msum do 60 j = k, n, msum wa2(j) = x(j) h = eps*dabs(wa2(j)) C C I have changed the difference formula to avoid h becoming too small. C C if (h .eq. zero) h = eps IF (DABS(H).LT.EPS) H = EPS x(j) = wa2(j) + h 60 continue call fcn(n,x,wa1,iflag) if (iflag .lt. 0) go to 100 do 80 j = k, n, msum x(j) = wa2(j) h = eps*dabs(wa2(j)) if (h .eq. zero) h = eps do 70 i = 1, n fjac(i,j) = zero if (i .ge. j - mu .and. i .le. j + ml) * fjac(i,j) = (wa1(i) - fvec(i))/h 70 continue 80 continue 90 continue 100 continue 110 continue return c c last card of subroutine fdjac1. c end ************************************************************************ ************************************************************************ *** from netlib, Thu Oct 31 01:57:13 EST 1991 *** subroutine rg(nm,n,a,wr,wi,matz,z,iv1,fv1,ierr) c integer n,nm,is1,is2,ierr,matz double precision a(nm,n),wr(n),wi(n),z(nm,n),fv1(n) integer iv1(n) c c this subroutine calls the recommended sequence of c subroutines from the eigensystem subroutine package (eispack) c to find the eigenvalues and eigenvectors (if desired) c of a real general matrix. c c on input c c nm must be set to the row dimension of the two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix a. c c a contains the real general matrix. c c matz is an integer variable set equal to zero if c only eigenvalues are desired. otherwise it is set to c any non-zero integer for both eigenvalues and eigenvectors. c c on output c c wr and wi contain the real and imaginary parts, c respectively, of the eigenvalues. complex conjugate c pairs of eigenvalues appear consecutively with the c eigenvalue having the positive imaginary part first. c c z contains the real and imaginary parts of the eigenvectors c if matz is not zero. if the j-th eigenvalue is real, the c j-th column of z contains its eigenvector. if the j-th c eigenvalue is complex with positive imaginary part, the c j-th and (j+1)-th columns of z contain the real and c imaginary parts of its eigenvector. the conjugate of this c vector is the eigenvector for the conjugate eigenvalue. c c ierr is an integer output variable set equal to an error c completion code described in the documentation for hqr c and hqr2. the normal completion code is zero. c c iv1 and fv1 are temporary storage arrays. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (n .le. nm) go to 10 ierr = 10 * n go to 50 c 10 call balanc(nm,n,a,is1,is2,fv1) call elmhes(nm,n,is1,is2,a,iv1) if (matz .ne. 0) go to 20 c .......... find eigenvalues only .......... call hqr(nm,n,is1,is2,a,wr,wi,ierr) go to 50 c .......... find both eigenvalues and eigenvectors .......... 20 call eltran(nm,n,is1,is2,a,iv1,z) call hqr2(nm,n,is1,is2,a,wr,wi,z,ierr) if (ierr .ne. 0) go to 50 call balbak(nm,n,is1,is2,fv1,n,z) 50 return end subroutine balanc(nm,n,a,low,igh,scale) c integer i,j,k,l,m,n,jj,nm,igh,low,iexc double precision a(nm,n),scale(n) double precision c,f,g,r,s,b2,radix logical noconv c c this subroutine is a translation of the algol procedure balance, c num. math. 13, 293-304(1969) by parlett and reinsch. c handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). c c this subroutine balances a real matrix and isolates c eigenvalues whenever possible. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c a contains the input matrix to be balanced. c c on output c c a contains the balanced matrix. c c low and igh are two integers such that a(i,j) c is equal to zero if c (1) i is greater than j and c (2) j=1,...,low-1 or i=igh+1,...,n. c c scale contains information determining the c permutations and scaling factors used. c c suppose that the principal submatrix in rows low through igh c has been balanced, that p(j) denotes the index interchanged c with j during the permutation step, and that the elements c of the diagonal matrix used are denoted by d(i,j). then c scale(j) = p(j), for j = 1,...,low-1 c = d(j,j), j = low,...,igh c = p(j) j = igh+1,...,n. c the order in which the interchanges are made is n to igh+1, c then 1 to low-1. c c note that 1 is returned for igh if igh is zero formally. c c the algol procedure exc contained in balance appears in c balanc in line. (note that the algol roles of identifiers c k,l have been reversed.) c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c radix = 16.0d0 c b2 = radix * radix k = 1 l = n go to 100 c .......... in-line procedure for row and c column exchange .......... 20 scale(m) = j if (j .eq. m) go to 50 c do 30 i = 1, l f = a(i,j) a(i,j) = a(i,m) a(i,m) = f 30 continue c do 40 i = k, n f = a(j,i) a(j,i) = a(m,i) a(m,i) = f 40 continue c 50 go to (80,130), iexc c .......... search for rows isolating an eigenvalue c and push them down .......... 80 if (l .eq. 1) go to 280 l = l - 1 c .......... for j=l step -1 until 1 do -- .......... 100 do 120 jj = 1, l j = l + 1 - jj c do 110 i = 1, l if (i .eq. j) go to 110 if (a(j,i) .ne. 0.0d0) go to 120 110 continue c m = l iexc = 1 go to 20 120 continue c go to 140 c .......... search for columns isolating an eigenvalue c and push them left .......... 130 k = k + 1 c 140 do 170 j = k, l c do 150 i = k, l if (i .eq. j) go to 150 if (a(i,j) .ne. 0.0d0) go to 170 150 continue c m = k iexc = 2 go to 20 170 continue c .......... now balance the submatrix in rows k to l .......... do 180 i = k, l 180 scale(i) = 1.0d0 c .......... iterative loop for norm reduction .......... 190 noconv = .false. c do 270 i = k, l c = 0.0d0 r = 0.0d0 c do 200 j = k, l if (j .eq. i) go to 200 c = c + dabs(a(j,i)) r = r + dabs(a(i,j)) 200 continue c .......... guard against zero c or r due to underflow .......... if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 g = r / radix f = 1.0d0 s = c + r 210 if (c .ge. g) go to 220 f = f * radix c = c * b2 go to 210 220 g = r * radix 230 if (c .lt. g) go to 240 f = f / radix c = c / b2 go to 230 c .......... now balance .......... 240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 g = 1.0d0 / f scale(i) = scale(i) * f noconv = .true. c do 250 j = k, n 250 a(i,j) = a(i,j) * g c do 260 j = 1, l 260 a(j,i) = a(j,i) * f c 270 continue c if (noconv) go to 190 c 280 low = k igh = l return end subroutine balbak(nm,n,low,igh,scale,m,z) c integer i,j,k,m,n,ii,nm,igh,low double precision scale(n),z(nm,m) double precision s c c this subroutine is a translation of the algol procedure balbak, c num. math. 13, 293-304(1969) by parlett and reinsch. c handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). c c this subroutine forms the eigenvectors of a real general c matrix by back transforming those of the corresponding c balanced matrix determined by balanc. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by balanc. c c scale contains information determining the permutations c and scaling factors used by balanc. c c m is the number of columns of z to be back transformed. c c z contains the real and imaginary parts of the eigen- c vectors to be back transformed in its first m columns. c c on output c c z contains the real and imaginary parts of the c transformed eigenvectors in its first m columns. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (m .eq. 0) go to 200 if (igh .eq. low) go to 120 c do 110 i = low, igh s = scale(i) c .......... left hand eigenvectors are back transformed c if the foregoing statement is replaced by c s=1.0d0/scale(i). .......... do 100 j = 1, m 100 z(i,j) = z(i,j) * s c 110 continue c ......... for i=low-1 step -1 until 1, c igh+1 step 1 until n do -- .......... 120 do 140 ii = 1, n i = ii if (i .ge. low .and. i .le. igh) go to 140 if (i .lt. low) i = low - ii k = scale(i) if (k .eq. i) go to 140 c do 130 j = 1, m s = z(i,j) z(i,j) = z(k,j) z(k,j) = s 130 continue c 140 continue c 200 return end subroutine elmhes(nm,n,low,igh,a,int) c integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1 double precision a(nm,n) double precision x,y integer int(igh) c c this subroutine is a translation of the algol procedure elmhes, c num. math. 12, 349-368(1968) by martin and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). c c given a real general matrix, this subroutine c reduces a submatrix situated in rows and columns c low through igh to upper hessenberg form by c stabilized elementary similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c a contains the input matrix. c c on output c c a contains the hessenberg matrix. the multipliers c which were used in the reduction are stored in the c remaining triangle under the hessenberg matrix. c c int contains information on the rows and columns c interchanged in the reduction. c only elements low through igh are used. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c la = igh - 1 kp1 = low + 1 if (la .lt. kp1) go to 200 c do 180 m = kp1, la mm1 = m - 1 x = 0.0d0 i = m c do 100 j = m, igh if (dabs(a(j,mm1)) .le. dabs(x)) go to 100 x = a(j,mm1) i = j 100 continue c int(m) = i if (i .eq. m) go to 130 c .......... interchange rows and columns of a .......... do 110 j = mm1, n y = a(i,j) a(i,j) = a(m,j) a(m,j) = y 110 continue c do 120 j = 1, igh y = a(j,i) a(j,i) = a(j,m) a(j,m) = y 120 continue c .......... end interchange .......... 130 if (x .eq. 0.0d0) go to 180 mp1 = m + 1 c do 160 i = mp1, igh y = a(i,mm1) if (y .eq. 0.0d0) go to 160 y = y / x a(i,mm1) = y c do 140 j = m, n 140 a(i,j) = a(i,j) - y * a(m,j) c do 150 j = 1, igh 150 a(j,m) = a(j,m) + y * a(j,i) c 160 continue c 180 continue c 200 return end subroutine eltran(nm,n,low,igh,a,int,z) c integer i,j,n,kl,mm,mp,nm,igh,low,mp1 double precision a(nm,igh),z(nm,n) integer int(igh) c c this subroutine is a translation of the algol procedure elmtrans, c num. math. 16, 181-204(1970) by peters and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). c c this subroutine accumulates the stabilized elementary c similarity transformations used in the reduction of a c real general matrix to upper hessenberg form by elmhes. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c a contains the multipliers which were used in the c reduction by elmhes in its lower triangle c below the subdiagonal. c c int contains information on the rows and columns c interchanged in the reduction by elmhes. c only elements low through igh are used. c c on output c c z contains the transformation matrix produced in the c reduction by elmhes. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c c .......... initialize z to identity matrix .......... do 80 j = 1, n c do 60 i = 1, n 60 z(i,j) = 0.0d0 c z(j,j) = 1.0d0 80 continue c kl = igh - low - 1 if (kl .lt. 1) go to 200 c .......... for mp=igh-1 step -1 until low+1 do -- .......... do 140 mm = 1, kl mp = igh - mm mp1 = mp + 1 c do 100 i = mp1, igh 100 z(i,mp) = a(i,mp-1) c i = int(mp) if (i .eq. mp) go to 140 c do 130 j = mp, igh z(mp,j) = z(i,j) z(i,j) = 0.0d0 130 continue c z(i,mp) = 1.0d0 140 continue c 200 return end subroutine hqr(nm,n,low,igh,h,wr,wi,ierr) C RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG) c integer i,j,k,l,m,n,en,ll,mm,na,nm,igh,itn,its,low,mp2,enm2,ierr double precision h(nm,n),wr(n),wi(n) double precision p,q,r,s,t,w,x,y,zz,norm,tst1,tst2 logical notlas c c this subroutine is a translation of the algol procedure hqr, c num. math. 14, 219-231(1970) by martin, peters, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 359-371(1971). c c this subroutine finds the eigenvalues of a real c upper hessenberg matrix by the qr method. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c h contains the upper hessenberg matrix. information about c the transformations used in the reduction to hessenberg c form by elmhes or orthes, if performed, is stored c in the remaining triangle under the hessenberg matrix. c c on output c c h has been destroyed. therefore, it must be saved c before calling hqr if subsequent calculation and c back transformation of eigenvectors is to be performed. c c wr and wi contain the real and imaginary parts, c respectively, of the eigenvalues. the eigenvalues c are unordered except that complex conjugate pairs c of values appear consecutively with the eigenvalue c having the positive imaginary part first. if an c error exit is made, the eigenvalues should be correct c for indices ierr+1,...,n. c c ierr is set to c zero for normal return, c j if the limit of 30*n iterations is exhausted c while the j-th eigenvalue is being sought. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated september 1989. c c ------------------------------------------------------------------ c ierr = 0 norm = 0.0d0 k = 1 c .......... store roots isolated by balanc c and compute matrix norm .......... do 50 i = 1, n c do 40 j = k, n 40 norm = norm + dabs(h(i,j)) c k = i if (i .ge. low .and. i .le. igh) go to 50 wr(i) = h(i,i) wi(i) = 0.0d0 50 continue c en = igh t = 0.0d0 itn = 30*n c .......... search for next eigenvalues .......... 60 if (en .lt. low) go to 1001 its = 0 na = en - 1 enm2 = na - 1 c .......... look for single small sub-diagonal element c for l=en step -1 until low do -- .......... 70 do 80 ll = low, en l = en + low - ll if (l .eq. low) go to 100 s = dabs(h(l-1,l-1)) + dabs(h(l,l)) if (s .eq. 0.0d0) s = norm tst1 = s tst2 = tst1 + dabs(h(l,l-1)) if (tst2 .eq. tst1) go to 100 80 continue c .......... form shift .......... 100 x = h(en,en) if (l .eq. en) go to 270 y = h(na,na) w = h(en,na) * h(na,en) if (l .eq. na) go to 280 if (itn .eq. 0) go to 1000 if (its .ne. 10 .and. its .ne. 20) go to 130 c .......... form exceptional shift .......... t = t + x c do 120 i = low, en 120 h(i,i) = h(i,i) - x c s = dabs(h(en,na)) + dabs(h(na,enm2)) x = 0.75d0 * s y = x w = -0.4375d0 * s * s 130 its = its + 1 itn = itn - 1 c .......... look for two consecutive small c sub-diagonal elements. c for m=en-2 step -1 until l do -- .......... do 140 mm = l, enm2 m = enm2 + l - mm zz = h(m,m) r = x - zz s = y - zz p = (r * s - w) / h(m+1,m) + h(m,m+1) q = h(m+1,m+1) - zz - r - s r = h(m+2,m+1) s = dabs(p) + dabs(q) + dabs(r) p = p / s q = q / s r = r / s if (m .eq. l) go to 150 tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) if (tst2 .eq. tst1) go to 150 140 continue c 150 mp2 = m + 2 c do 160 i = mp2, en h(i,i-2) = 0.0d0 if (i .eq. mp2) go to 160 h(i,i-3) = 0.0d0 160 continue c .......... double qr step involving rows l to en and c columns m to en .......... do 260 k = m, na notlas = k .ne. na if (k .eq. m) go to 170 p = h(k,k-1) q = h(k+1,k-1) r = 0.0d0 if (notlas) r = h(k+2,k-1) x = dabs(p) + dabs(q) + dabs(r) if (x .eq. 0.0d0) go to 260 p = p / x q = q / x r = r / x 170 s = dsign(dsqrt(p*p+q*q+r*r),p) if (k .eq. m) go to 180 h(k,k-1) = -s * x go to 190 180 if (l .ne. m) h(k,k-1) = -h(k,k-1) 190 p = p + s x = p / s y = q / s zz = r / s q = q / p r = r / p if (notlas) go to 225 c .......... row modification .......... do 200 j = k, EN p = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y 200 continue c j = min0(en,k+3) c .......... column modification .......... do 210 i = L, j p = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q 210 continue go to 255 225 continue c .......... row modification .......... do 230 j = k, EN p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y h(k+2,j) = h(k+2,j) - p * zz 230 continue c j = min0(en,k+3) c .......... column modification .......... do 240 i = L, j p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q h(i,k+2) = h(i,k+2) - p * r 240 continue 255 continue c 260 continue c go to 70 c .......... one root found .......... 270 wr(en) = x + t wi(en) = 0.0d0 en = na go to 60 c .......... two roots found .......... 280 p = (y - x) / 2.0d0 q = p * p + w zz = dsqrt(dabs(q)) x = x + t if (q .lt. 0.0d0) go to 320 c .......... real pair .......... zz = p + dsign(zz,p) wr(na) = x + zz wr(en) = wr(na) if (zz .ne. 0.0d0) wr(en) = x - w / zz wi(na) = 0.0d0 wi(en) = 0.0d0 go to 330 c .......... complex pair .......... 320 wr(na) = x + p wr(en) = x + p wi(na) = zz wi(en) = -zz 330 en = enm2 go to 60 c .......... set error -- all eigenvalues have not c converged after 30*n iterations .......... 1000 ierr = en 1001 return end subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr) c integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn, x igh,itn,its,low,mp2,enm2,ierr double precision h(nm,n),wr(n),wi(n),z(nm,n) double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2 logical notlas c c this subroutine is a translation of the algol procedure hqr2, c num. math. 16, 181-204(1970) by peters and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a real upper hessenberg matrix by the qr method. the c eigenvectors of a real general matrix can also be found c if elmhes and eltran or orthes and ortran have c been used to reduce this general matrix to hessenberg form c and to accumulate the similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c low and igh are integers determined by the balancing c subroutine balanc. if balanc has not been used, c set low=1, igh=n. c c h contains the upper hessenberg matrix. c c z contains the transformation matrix produced by eltran c after the reduction by elmhes, or by ortran after the c reduction by orthes, if performed. if the eigenvectors c of the hessenberg matrix are desired, z must contain the c identity matrix. c c on output c c h has been destroyed. c c wr and wi contain the real and imaginary parts, c respectively, of the eigenvalues. the eigenvalues c are unordered except that complex conjugate pairs c of values appear consecutively with the eigenvalue c having the positive imaginary part first. if an c error exit is made, the eigenvalues should be correct c for indices ierr+1,...,n. c c z contains the real and imaginary parts of the eigenvectors. c if the i-th eigenvalue is real, the i-th column of z c contains its eigenvector. if the i-th eigenvalue is complex c with positive imaginary part, the i-th and (i+1)-th c columns of z contain the real and imaginary parts of its c eigenvector. the eigenvectors are unnormalized. if an c error exit is made, none of the eigenvectors has been found. c c ierr is set to c zero for normal return, c j if the limit of 30*n iterations is exhausted c while the j-th eigenvalue is being sought. c c calls cdiv for complex division. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 norm = 0.0d0 k = 1 c .......... store roots isolated by balanc c and compute matrix norm .......... do 50 i = 1, n c do 40 j = k, n 40 norm = norm + dabs(h(i,j)) c k = i if (i .ge. low .and. i .le. igh) go to 50 wr(i) = h(i,i) wi(i) = 0.0d0 50 continue c en = igh t = 0.0d0 itn = 30*n c .......... search for next eigenvalues .......... 60 if (en .lt. low) go to 340 its = 0 na = en - 1 enm2 = na - 1 c .......... look for single small sub-diagonal element c for l=en step -1 until low do -- .......... 70 do 80 ll = low, en l = en + low - ll if (l .eq. low) go to 100 s = dabs(h(l-1,l-1)) + dabs(h(l,l)) if (s .eq. 0.0d0) s = norm tst1 = s tst2 = tst1 + dabs(h(l,l-1)) if (tst2 .eq. tst1) go to 100 80 continue c .......... form shift .......... 100 x = h(en,en) if (l .eq. en) go to 270 y = h(na,na) w = h(en,na) * h(na,en) if (l .eq. na) go to 280 if (itn .eq. 0) go to 1000 if (its .ne. 10 .and. its .ne. 20) go to 130 c .......... form exceptional shift .......... t = t + x c do 120 i = low, en 120 h(i,i) = h(i,i) - x c s = dabs(h(en,na)) + dabs(h(na,enm2)) x = 0.75d0 * s y = x w = -0.4375d0 * s * s 130 its = its + 1 itn = itn - 1 c .......... look for two consecutive small c sub-diagonal elements. c for m=en-2 step -1 until l do -- .......... do 140 mm = l, enm2 m = enm2 + l - mm zz = h(m,m) r = x - zz s = y - zz p = (r * s - w) / h(m+1,m) + h(m,m+1) q = h(m+1,m+1) - zz - r - s r = h(m+2,m+1) s = dabs(p) + dabs(q) + dabs(r) p = p / s q = q / s r = r / s if (m .eq. l) go to 150 tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) if (tst2 .eq. tst1) go to 150 140 continue c 150 mp2 = m + 2 c do 160 i = mp2, en h(i,i-2) = 0.0d0 if (i .eq. mp2) go to 160 h(i,i-3) = 0.0d0 160 continue c .......... double qr step involving rows l to en and c columns m to en .......... do 260 k = m, na notlas = k .ne. na if (k .eq. m) go to 170 p = h(k,k-1) q = h(k+1,k-1) r = 0.0d0 if (notlas) r = h(k+2,k-1) x = dabs(p) + dabs(q) + dabs(r) if (x .eq. 0.0d0) go to 260 p = p / x q = q / x r = r / x 170 s = dsign(dsqrt(p*p+q*q+r*r),p) if (k .eq. m) go to 180 h(k,k-1) = -s * x go to 190 180 if (l .ne. m) h(k,k-1) = -h(k,k-1) 190 p = p + s x = p / s y = q / s zz = r / s q = q / p r = r / p if (notlas) go to 225 c .......... row modification .......... do 200 j = k, n p = h(k,j) + q * h(k+1,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y 200 continue c j = min0(en,k+3) c .......... column modification .......... do 210 i = 1, j p = x * h(i,k) + y * h(i,k+1) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q 210 continue c .......... accumulate transformations .......... do 220 i = low, igh p = x * z(i,k) + y * z(i,k+1) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q 220 continue go to 255 225 continue c .......... row modification .......... do 230 j = k, n p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) h(k,j) = h(k,j) - p * x h(k+1,j) = h(k+1,j) - p * y h(k+2,j) = h(k+2,j) - p * zz 230 continue c j = min0(en,k+3) c .......... column modification .......... do 240 i = 1, j p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) h(i,k) = h(i,k) - p h(i,k+1) = h(i,k+1) - p * q h(i,k+2) = h(i,k+2) - p * r 240 continue c .......... accumulate transformations .......... do 250 i = low, igh p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) z(i,k) = z(i,k) - p z(i,k+1) = z(i,k+1) - p * q z(i,k+2) = z(i,k+2) - p * r 250 continue 255 continue c 260 continue c go to 70 c .......... one root found .......... 270 h(en,en) = x + t wr(en) = h(en,en) wi(en) = 0.0d0 en = na go to 60 c .......... two roots found .......... 280 p = (y - x) / 2.0d0 q = p * p + w zz = dsqrt(dabs(q)) h(en,en) = x + t x = h(en,en) h(na,na) = y + t if (q .lt. 0.0d0) go to 320 c .......... real pair .......... zz = p + dsign(zz,p) wr(na) = x + zz wr(en) = wr(na) if (zz .ne. 0.0d0) wr(en) = x - w / zz wi(na) = 0.0d0 wi(en) = 0.0d0 x = h(en,na) s = dabs(x) + dabs(zz) p = x / s q = zz / s r = dsqrt(p*p+q*q) p = p / r q = q / r c .......... row modification .......... do 290 j = na, n zz = h(na,j) h(na,j) = q * zz + p * h(en,j) h(en,j) = q * h(en,j) - p * zz 290 continue c .......... column modification .......... do 300 i = 1, en zz = h(i,na) h(i,na) = q * zz + p * h(i,en) h(i,en) = q * h(i,en) - p * zz 300 continue c .......... accumulate transformations .......... do 310 i = low, igh zz = z(i,na) z(i,na) = q * zz + p * z(i,en) z(i,en) = q * z(i,en) - p * zz 310 continue c go to 330 c .......... complex pair .......... 320 wr(na) = x + p wr(en) = x + p wi(na) = zz wi(en) = -zz 330 en = enm2 go to 60 c .......... all roots found. backsubstitute to find c vectors of upper triangular form .......... 340 if (norm .eq. 0.0d0) go to 1001 c .......... for en=n step -1 until 1 do -- .......... do 800 nn = 1, n en = n + 1 - nn p = wr(en) q = wi(en) na = en - 1 if (q) 710, 600, 800 c .......... real vector .......... 600 m = en h(en,en) = 1.0d0 if (na .eq. 0) go to 800 c .......... for i=en-1 step -1 until 1 do -- .......... do 700 ii = 1, na i = en - ii w = h(i,i) - p r = 0.0d0 c do 610 j = m, en 610 r = r + h(i,j) * h(j,en) c if (wi(i) .ge. 0.0d0) go to 630 zz = w s = r go to 700 630 m = i if (wi(i) .ne. 0.0d0) go to 640 t = w if (t .ne. 0.0d0) go to 635 tst1 = norm t = tst1 632 t = 0.01d0 * t tst2 = norm + t if (tst2 .gt. tst1) go to 632 635 h(i,en) = -r / t go to 680 c .......... solve real equations .......... 640 x = h(i,i+1) y = h(i+1,i) q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) t = (x * s - zz * r) / q h(i,en) = t if (dabs(x) .le. dabs(zz)) go to 650 h(i+1,en) = (-r - w * t) / x go to 680 650 h(i+1,en) = (-s - y * t) / zz c c .......... overflow control .......... 680 t = dabs(h(i,en)) if (t .eq. 0.0d0) go to 700 tst1 = t tst2 = tst1 + 1.0d0/tst1 if (tst2 .gt. tst1) go to 700 do 690 j = i, en h(j,en) = h(j,en)/t 690 continue c 700 continue c .......... end real vector .......... go to 800 c .......... complex vector .......... 710 m = na c .......... last vector component chosen imaginary so that c eigenvector matrix is triangular .......... if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720 h(na,na) = q / h(en,na) h(na,en) = -(h(en,en) - p) / h(en,na) go to 730 720 call cdiv(0.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) 730 h(en,na) = 0.0d0 h(en,en) = 1.0d0 enm2 = na - 1 if (enm2 .eq. 0) go to 800 c .......... for i=en-2 step -1 until 1 do -- .......... do 795 ii = 1, enm2 i = na - ii w = h(i,i) - p ra = 0.0d0 sa = 0.0d0 c do 760 j = m, en ra = ra + h(i,j) * h(j,na) sa = sa + h(i,j) * h(j,en) 760 continue c if (wi(i) .ge. 0.0d0) go to 770 zz = w r = ra s = sa go to 795 770 m = i if (wi(i) .ne. 0.0d0) go to 780 call cdiv(-ra,-sa,w,q,h(i,na),h(i,en)) go to 790 c .......... solve complex equations .......... 780 x = h(i,i+1) y = h(i+1,i) vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q vi = (wr(i) - p) * 2.0d0 * q if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784 tst1 = norm * (dabs(w) + dabs(q) + dabs(x) x + dabs(y) + dabs(zz)) vr = tst1 783 vr = 0.01d0 * vr tst2 = tst1 + vr if (tst2 .gt. tst1) go to 783 784 call cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi, x h(i,na),h(i,en)) if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785 h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x go to 790 785 call cdiv(-r-y*h(i,na),-s-y*h(i,en),zz,q, x h(i+1,na),h(i+1,en)) c c .......... overflow control .......... 790 t = dmax1(dabs(h(i,na)), dabs(h(i,en))) if (t .eq. 0.0d0) go to 795 tst1 = t tst2 = tst1 + 1.0d0/tst1 if (tst2 .gt. tst1) go to 795 do 792 j = i, en h(j,na) = h(j,na)/t h(j,en) = h(j,en)/t 792 continue c 795 continue c .......... end complex vector .......... 800 continue c .......... end back substitution. c