SUBROUTINE CXFFT3N(C,MX,MY,MZ,ISIGN) C Parameters: INTEGER MXTRIG PARAMETER(MXTRIG=1000) C Arguments: INTEGER ISIGN,MX,MY,MZ REAL C(*) C C Local scalars: INTEGER INC,IND,JUMP,K,LOT,MXOLD,MYOLD,MZOLD C Local arrays: REAL TRIG1(MXTRIG),TRIG2(MXTRIG),TRIG3(MXTRIG) SAVE TRIG1,TRIG2,TRIG3,MXOLD,MYOLD,MZOLD C*********************************************************************** C Interface routine for "Generalized Prime Factor Algorithm" FFT code C of Temperton. C Interface written by P.J.Flatau C C History: C C*********************************************************************** data mxold, myold, mzold/0,0,0/ if(mx.ne.mxold .or. my.ne.myold .or. mz.ne.mzold) then mxold=mx myold=my mzold=mz CALL SETGPFA(trig1,mX) CALL SETGPFA(trig2,mY) CALL SETGPFA(trig3,mZ) endif C --- first dimension INC=2 JUMP=2*MX LOT=MY*MZ CALL GPFA(C(1),c(2),TRIG1,INC,JUMP,MX,LOT,ISIGN) C --- second dimension INC=MX*2 JUMP=2 LOT=MX C --- one plane at a time DO 10 K=1,MZ ind=1+2*MX*MY*(K-1) CALL GPFA(C(IND),C(IND+1),TRIG2,INC,JUMP, & MY,LOT,ISIGN) 10 CONTINUE C --- third dimension INC=MX*MY*2 JUMP=2 LOT=MX*MY CALL GPFA(C(1),C(2),TRIG3,INC,JUMP,MZ,LOT, & ISIGN) RETURN END SUBROUTINE DIELEC(CDIEL,WAVE,IDVOUT,CFLEPS,CXEPS, & MXCOMP,MXWAVT,NCOMP,E1A,E2A,WVA) IMPLICIT NONE C Arguments: INTEGER MXCOMP,MXWAVT CHARACTER*6 CDIEL CHARACTER*60 CFLEPS(MXCOMP) COMPLEX CXEPS(MXCOMP) INTEGER IDVOUT,NCOMP REAL WAVE REAL E1A(MXWAVT),E2A(MXWAVT), & WVA(MXWAVT) C Local Variables: LOGICAL INIT CHARACTER*80 CDESCR COMPLEX CXI INTEGER I,J,JJ,NWAVT,IWV,IREN,IIMN,IEPS1,IEPS2,ICOL REAL DUM1,DUM2,E1,E2,TEMP REAL XIN(10) EXTERNAL ERRMSG SAVE CXI,IEPS1,INIT,NWAVT *********************************************************************** C Given: C CDIEL = 'TABLES' to read from input tables C 'H2OICE' for H2O ice C 'H2OLIQ' for H2O liquid C WAVE = wavelength (micron) C IDVOUT=output unit number C CFLEPS(1-NCOMP) = names of files containing dielectric data C MXCOMP,MXWAVT = dimensioning information C NCOMP = number of components C NWAV = number of wavelengths C E1A, E2A, WVA = scratch arrays C Returns: C CXEPS(1-NCOMP) = dielectric constant for materials C 1-NCOMP C NOTE: C If CDIEL='TABLES', then it is assumed that file(s) containing C the table(s) will have following format: C line1 = description (CHARACTER*80) C line2 = IWV, IREN, IIMN, IEPS1, IEPS2 C where IWV = column in which wavelengths are tabulated C IREN = col. in which Re(N) is tabulated (0 if not) C IIMN = col. in which Im(N) is tabulated (0 if not) C IEPS1= col. in which Re(EPS) is tabulated (0 if not) C IEPS2= col. in which Im(EPS) is tabulated (0 if not) C line3 = header line (will not be read) C line4... = data, according to columns specified on line2 C wavelengths must be monotonic, either increasing or decreasing C C B.T.Draine, Princeton Univ. Obs. C History: C 90.12.04 (BTD): Rewritten to allow H2OICE and H2OLIQ options. C 90.12.21 (BTD): Correct handling of 'TABLES' option when more than C one wavelength is considered. C 91.05.02 (BTD): Added IDVOUT to argument list for subr. INTERP C 91.09.12 (BTD): Moved declaration of MXCOMP and MWAVT ahead of C declaration of CFLEPS and CXEPS C 93.12.06 (BTD): Added IEPS1,INIT, and NWAVT to SAVE statement C (without this did not run properly on SGI Indigo) C 96.12.16 (BTD): Corrected temperature to T=250K in output statement C for H2OICE option C 98.12.21 (BTD): changed dimension of CFLPAR from CHARACTER*40 to C CHARACTER*60 to allow longer file names C (also changed in reapar.f and DDSCAT.f) C Copyright (C) 1993,1996,1998 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DATA CXI/(0.,1.)/,INIT/.TRUE./ IF(CDIEL.EQ.'TABLES')THEN IF(INIT.OR.NCOMP.GT.1)THEN DO 1000 J=1,NCOMP OPEN(UNIT=3,FILE=CFLEPS(J),STATUS='OLD') C Read header line: READ(3,9000)CDESCR WRITE(IDVOUT,9000)CDESCR C Read line specifying columns for wavelength,Re(n),Im(n),Re(eps),Im(eps): READ(3,*)IWV,IREN,IIMN,IEPS1,IEPS2 ICOL=IWV IF(IREN.GT.ICOL)ICOL=IREN IF(IIMN.GT.ICOL)ICOL=IIMN IF(IEPS1.GT.ICOL)ICOL=IEPS1 IF(IEPS2.GT.ICOL)ICOL=IEPS2 C Skip header line: READ(3,*) DO 0500 I=1,MXWAVT READ(3,*,END=0600)(XIN(JJ),JJ=1,ICOL) WVA(I)=XIN(IWV) IF(IEPS1.GT.0)THEN E1A(I)=XIN(IEPS1) E2A(I)=XIN(IEPS2) ELSE E1A(I)=XIN(IREN) E2A(I)=XIN(IIMN) ENDIF NWAVT=I 0500 CONTINUE C Check whether there is unread data remaining in file READ(3,*,END=0600)(XIN(JJ),JJ=1,ICOL) C If this point is reached, apparently unread data remains, so C issue warning: CALL ERRMSG('WARNING','DIELEC', & 'parameter MXWAVT not large enough to read full dielec file') 0600 CLOSE(3) C C*** Have completed reading in table for this composition. C Now interpolate CALL INTERP(WAVE,E1,E2,WVA,E1A,E2A, & IDVOUT,MXWAVT,NWAVT) IF(IEPS1.GT.0)THEN CXEPS(J)=E1+CXI*E2 ELSE CXEPS(J)=(E1**2-E2**2)+2.*CXI*E1*E2 ENDIF 1000 CONTINUE INIT=.FALSE. ELSE C*** Perform this only if NCOMP=1 and previously initialized: CALL INTERP(WAVE,E1,E2,WVA,E1A,E2A,IDVOUT,MXWAVT,NWAVT) IF(IEPS1.GT.0)THEN CXEPS(1)=E1+CXI*E2 ELSE CXEPS(1)=(E1**2-E2**2)+2.*CXI*E1*E2 ENDIF ENDIF ELSEIF(CDIEL.EQ.'H2OICE')THEN IF(NCOMP.GT.1)CALL ERRMSG('FATAL','DIELE', & 'H2OICE with NCOMP>1??') C For the moment set T=250K TEMP=250. WRITE(CDESCR,FMT='(A)')' H2O ice at T=250K' WRITE(IDVOUT,9000)CDESCR CALL REFICE(0,WAVE,TEMP,E1,E2,DUM1,DUM2) CXEPS(1)=E1**2-E2**2+2.*CXI*E1*E2 RETURN ELSEIF(CDIEL.EQ.'H2OLIQ')THEN WRITE(CDESCR,FMT='(A)')' H2O liquid at T=280K' WRITE(IDVOUT,9000)CDESCR IF(NCOMP.GT.1)CALL ERRMSG('FATAL','DIELE', & 'H2OLIQ with NCOMP>1??') C For the moment set T=280K TEMP=280. CALL REFWAT(0,WAVE,TEMP,E1,E2,DUM1,DUM2) CXEPS(1)=E1**2-E2**2+2.*CXI*E1*E2 RETURN ELSE CALL ERRMSG('FATAL','DIELE','unrecognized value of CDIEL') ENDIF RETURN 9000 FORMAT(A80) END SUBROUTINE DIVIDE(CDIVID,X1,X2,NX,MXARY,ARY) IMPLICIT NONE C Arguments: REAL X1,X2 INTEGER MXARY,NX CHARACTER CDIVID*3 REAL ARY(MXARY) C C Local variables: INTEGER I REAL DELTA C C External Subroutines: EXTERNAL ERRMSG,WRIMSG C C Intrinsic Functions: INTRINSIC REAL C*********************************************************************** C Given: C CDIVID='LIN','INV', or 'LOG' C X1 = lower limit to interval C X2 = upper limit to interval C NX = number of elements C MXARY=dimensioning information for array ARY C C Returns: C ARY(1-NX)=vector of points with C ARY(1)=X1 C ARY(NX)=X2 C ARY(J) values spaced either C linearly (if CDIVID.EQ.'LIN') C linearly in 1/X (if CDIVID.EQ.'INV') C logarithmically (if CDIVID.EQ.'LOG') C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** IF(NX.GT.MXARY)THEN CALL ERRMSG('FATAL','DIVIDE',' NX .GT. MXARY ') ENDIF C IF(X1.GT.X2)THEN C CALL ERRMSG('FATAL','DIVIDE',' X1 .GT. X2 ') C ENDIF IF(NX.LT.1)THEN CALL ERRMSG('FATAL','DIVIDE',' NX .LT. 1 ') ENDIF IF(NX.EQ.1)THEN CALL WRIMSG('DIVIDE',' Only one element initialized ') ARY(1)=X1 ELSE IF(CDIVID.EQ.'LIN')THEN DELTA=(X2-X1)/REAL(NX-1) DO 10 I=1,NX ARY(I)=X1+DELTA*REAL(I-1) 10 CONTINUE ELSEIF(CDIVID.EQ.'INV')THEN DELTA=(1./X2-1./X1)/REAL(NX-1) DO 20 I=1,NX ARY(I)=1./(1./X1+DELTA*REAL(I-1)) 20 CONTINUE ELSEIF(CDIVID.EQ.'LOG')THEN DELTA=ALOG(X2/X1)/REAL(NX-1) DO 30 I=1,NX ARY(I)=EXP(ALOG(X1)+DELTA*REAL(I-1)) 30 CONTINUE ENDIF ENDIF RETURN END SUBROUTINE EIGV (IBAL,A,KA,N,WR,WI,ZR,ZI,IERR) C Arguments: INTEGER IBAL,IERR,KA,N REAL A(KA,N), WR(N), WI(N), ZR(KA,N), ZI(KA,N) C Local variables: INTEGER IGH,J,K,KP1,LOW C----------------------------------------------------------------------- C EIGENVALUES AND EIGENVECTORS OF REAL MATRICES C C This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C----------------------------------------------------------------------- C LOW = 1 IGH = N IF (IBAL .NE. 0) CALL BALANC (KA,N,A,LOW,IGH,ZI) CALL ELMHS0 (KA,N,LOW,IGH,A,WR) CALL ELTRN0 (KA,N,LOW,IGH,A,WR,ZR) CALL HQR2 (KA,N,LOW,IGH,A,WR,WI,ZR,IERR) IF (IERR .NE. 0) RETURN IF (IBAL .NE. 0) CALL BALBAK (KA,N,LOW,IGH,ZI,N,ZR) C DO 30 K = 1,N IF (WI(K)) 30,10,20 10 DO 11 J = 1,N 11 ZI(J,K) = 0.0 GO TO 30 20 KP1 = K + 1 DO 21 J = 1,N ZI(J,K) = ZR(J,KP1) ZR(J,KP1) = ZR(J,K) 21 ZI(J,KP1) = -ZI(J,K) 30 CONTINUE RETURN END SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL A(NM,N),SCALE(N) REAL C,F,G,R,S,B2,RADIX INTEGER IPMPAR C REAL ABS LOGICAL NOCONV C----------------------------------------------------------------------- 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 This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C----------------------------------------------------------------------- C C ********** RADIX IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE BASE OF THE MACHINE FLOATING POINT REPRESENTATION. C RADIX = IPMPAR(4) C C ********** 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.0) 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.0) 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.0 C ********** ITERATIVE LOOP FOR NORM REDUCTION ********** 190 NOCONV = .FALSE. C DO 270 I = K, L C = 0.0 R = 0.0 C DO 200 J = K, L IF (J .EQ. I) GO TO 200 C = C + ABS(A(J,I)) R = R + ABS(A(I,J)) 200 CONTINUE C ********** GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ********** IF (C .EQ. 0.0 .OR. R .EQ. 0.0) GO TO 270 G = R / RADIX F = 1.0 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.95 * S) GO TO 270 G = 1.0 / 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 REAL SCALE(N),Z(NM,M) REAL S C----------------------------------------------------------------------- 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 This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 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.0/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 ELMHS0(NM,N,LOW,IGH,A,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 REAL A(NM,N) REAL X,Y C REAL ABS REAL INT(IGH) C----------------------------------------------------------------------- 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 This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 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.0 I = M C DO 100 J = M, IGH IF (ABS(A(J,MM1)) .LE. ABS(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.0) GO TO 180 MP1 = M + 1 C DO 160 I = MP1, IGH Y = A(I,MM1) IF (Y .EQ. 0.0) 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 ELTRN0(NM,N,LOW,IGH,A,INT,Z) C INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 REAL A(NM,IGH),Z(NM,N) REAL INT(IGH) C----------------------------------------------------------------------- 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 ELMHS0. 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 ELMHS0 IN ITS LOWER TRIANGLE C BELOW THE SUBDIAGONAL, C C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS C INTERCHANGED IN THE REDUCTION BY ELMHS0. 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 ELMHS0. C C This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C----------------------------------------------------------------------- C C ********** INITIALIZE Z TO IDENTITY MATRIX ********** DO 80 I = 1, N C DO 60 J = 1, N 60 Z(I,J) = 0.0 C Z(I,I) = 1.0 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.0 130 CONTINUE C Z(I,MP) = 1.0 140 CONTINUE C 200 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,ITS,LOW,MP2,ENM2,IERR REAL H(NM,N),WR(N),WI(N),Z(NM,N) REAL P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,MACHEP,SPMPAR C REAL SQRT,ABS C INTEGER MIN0 LOGICAL NOTLAS COMPLEX Z3 C COMPLEX CMPLX C REAL REAL,AIMAG C----------------------------------------------------------------------- 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 J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C ARITHMETIC IS REAL EXCEPT FOR THE REPLACEMENT OF THE ALGOL C PROCEDURE CDIV BY COMPLEX DIVISION. C C This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C----------------------------------------------------------------------- C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER. ASSIGN C MACHEP THE VALUE U WHERE U IS THE SMALLEST POSITIVE C FLOATING POINT NUMBER SUCH THAT 1.0 + U .GT. 1.0 . C MACHEP = SPMPAR(1) C C ********** C IERR = 0 NORM = 0.0 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 + ABS(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.0 50 CONTINUE C EN = IGH T = 0.0 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 = ABS(H(L-1,L-1)) + ABS(H(L,L)) IF (S .EQ. 0.0) S = NORM IF (ABS(H(L,L-1)) .LE. MACHEP * S) 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 (ITS .EQ. 30) 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 = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75 * S Y = X W = -0.4375 * S * S 130 ITS = ITS + 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 = ABS(P) + ABS(Q) + ABS(R) P = P / S Q = Q / S R = R / S IF (M .EQ. L) GO TO 150 IF (ABS(H(M,M-1)) * (ABS(Q) + ABS(R)) .LE. MACHEP * ABS(P) X * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0 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.0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X .EQ. 0.0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = SQRT(P*P + Q*Q + R*R) IF (P .LT. 0.0) S = -S 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 C ********** ROW MODIFICATION ********** DO 210 J = K, N P = H(K,J) + Q * H(K+1,J) IF (.NOT. NOTLAS) GO TO 200 P = P + R * H(K+2,J) H(K+2,J) = H(K+2,J) - P * ZZ 200 H(K+1,J) = H(K+1,J) - P * Y H(K,J) = H(K,J) - P * X 210 CONTINUE C J = MIN0(EN,K+3) C ********** COLUMN MODIFICATION ********** DO 230 I = 1, J P = X * H(I,K) + Y * H(I,K+1) IF (.NOT. NOTLAS) GO TO 220 P = P + ZZ * H(I,K+2) H(I,K+2) = H(I,K+2) - P * R 220 H(I,K+1) = H(I,K+1) - P * Q H(I,K) = H(I,K) - P 230 CONTINUE C ********** ACCUMULATE TRANSFORMATIONS ********** DO 250 I = LOW, IGH P = X * Z(I,K) + Y * Z(I,K+1) IF (.NOT. NOTLAS) GO TO 240 P = P + ZZ * Z(I,K+2) Z(I,K+2) = Z(I,K+2) - P * R 240 Z(I,K+1) = Z(I,K+1) - P * Q Z(I,K) = Z(I,K) - P 250 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.0 EN = NA GO TO 60 C ********** TWO ROOTS FOUND ********** 280 P = (Y - X) / 2.0 Q = P * P + W ZZ = SQRT(ABS(Q)) H(EN,EN) = X + T X = H(EN,EN) H(NA,NA) = Y + T IF (Q .LT. 0.0) GO TO 320 C ********** REAL PAIR ********** IF (P .LT. 0.0) ZZ = -ZZ ZZ = P + ZZ WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ .NE. 0.0) WR(EN) = X - W / ZZ WI(NA) = 0.0 WI(EN) = 0.0 X = H(EN,NA) S = ABS(X) + ABS(ZZ) P = X / S Q = ZZ / S R = SQRT(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.0) 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.0 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 = H(I,EN) IF (M .GT. NA) GO TO 620 C DO 610 J = M, NA 610 R = R + H(I,J) * H(J,EN) C 620 IF (WI(I) .GE. 0.0) GO TO 630 ZZ = W S = R GO TO 700 630 M = I IF (WI(I) .NE. 0.0) GO TO 640 T = W IF (W .EQ. 0.0) T = MACHEP * NORM H(I,EN) = -R / T GO TO 700 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 (ABS(X) .LE. ABS(ZZ)) GO TO 650 H(I+1,EN) = (-R - W * T) / X GO TO 700 650 H(I+1,EN) = (-S - Y * T) / ZZ 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 (ABS(H(EN,NA)) .LE. ABS(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 Z3 = CMPLX(0.0,-H(NA,EN)) / CMPLX(H(NA,NA)-P,Q) H(NA,NA) = REAL(Z3) H(NA,EN) = AIMAG(Z3) 730 H(EN,NA) = 0.0 H(EN,EN) = 1.0 ENM2 = NA - 1 IF (ENM2 .EQ. 0) GO TO 800 C ********** FOR I=EN-2 STEP -1 UNTIL 1 DO -- ********** DO 790 II = 1, ENM2 I = NA - II W = H(I,I) - P RA = 0.0 SA = H(I,EN) C DO 760 J = M, NA RA = RA + H(I,J) * H(J,NA) SA = SA + H(I,J) * H(J,EN) 760 CONTINUE C IF (WI(I) .GE. 0.0) GO TO 770 ZZ = W R = RA S = SA GO TO 790 770 M = I IF (WI(I) .NE. 0.0) GO TO 780 Z3 = CMPLX(-RA,-SA) / CMPLX(W,Q) H(I,NA) = REAL(Z3) H(I,EN) = AIMAG(Z3) 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.0 * Q IF (VR .EQ. 0.0 .AND. VI .EQ. 0.0) VR = MACHEP * NORM X * (ABS(W) + ABS(Q) + ABS(X) + ABS(Y) + ABS(ZZ)) Z3 = CMPLX(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA) / CMPLX(VR,VI) H(I,NA) = REAL(Z3) H(I,EN) = AIMAG(Z3) IF (ABS(X) .LE. ABS(ZZ) + ABS(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 Z3 = CMPLX(-R-Y*H(I,NA),-S-Y*H(I,EN)) / CMPLX(ZZ,Q) H(I+1,NA) = REAL(Z3) H(I+1,EN) = AIMAG(Z3) 790 CONTINUE C ********** END COMPLEX VECTOR ********** 800 CONTINUE C ********** END BACK SUBSTITUTION. C VECTORS OF ISOLATED ROOTS ********** DO 840 I = 1, N IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840 C DO 820 J = I, N 820 Z(I,J) = H(I,J) C 840 CONTINUE C ********** MULTIPLY BY TRANSFORMATION MATRIX TO GIVE C VECTORS OF ORIGINAL FULL MATRIX. C FOR J=N STEP -1 UNTIL LOW DO -- ********** DO 880 JJ = LOW, N J = N + LOW - JJ M = MIN0(J,IGH) C DO 880 I = LOW, IGH ZZ = 0.0 C DO 860 K = LOW, M 860 ZZ = ZZ + Z(I,K) * H(K,J) C Z(I,J) = ZZ 880 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 1000 IERR = EN 1001 RETURN END INTEGER FUNCTION IPMPAR(I) C Arguments: INTEGER I C Local variables: INTEGER IMACH(10) C C----------------------------------------------------------------------- C C IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER C THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER C HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ... C C INTEGERS. C C ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM C C SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1. C C IPMPAR(1) = A, THE BASE. C C IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS. C C IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C C IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING C POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE C NONZERO NUMBERS ARE REPRESENTED IN THE FORM C C SIGN (B**E) * (X(1)/B + ... + X(M)/B**M) C C WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M, C X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX. C C IPMPAR(4) = B, THE BASE. C C SINGLE-PRECISION C C IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(7) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C C IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(10) = EMAX, THE LARGEST EXPONENT E. C C----------------------------------------------------------------------- C C TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED, ACTIVATE C THE DATA STATMENTS FOR THE COMPUTER BY REMOVING THE C FROM C COLUMN 1. (ALL THE OTHER DATA STATEMENTS SHOULD HAVE C IN C COLUMN 1.) C C IF DATA STATEMENTS ARE NOT GIVEN FOR THE COMPUTER BEING USED, C THEN THE FORTRAN MANUAL FOR THE COMPUTER NORMALLY GIVES THE C CONSTANTS IPMPAR(1), IPMPAR(2), AND IPMPAR(3) FOR THE INTEGER C ARITHMETIC. HOWEVER, HELP MAY BE NEEDED TO OBTAIN THE CONSTANTS C IPMPAR(4),...,IPMPAR(10) FOR THE SINGLE AND DOUBLE PRECISION C ARITHMETICS. THE SUBROUTINES MACH AND RADIX ARE PROVIDED FOR C THIS PURPOSE. C C----------------------------------------------------------------------- C C IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY C P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES). C IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE C FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES. C C This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C History: C 95.11.16 (BTD) Explicit declaration of all variables. C C----------------------------------------------------------------------- C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 33 / C DATA IMACH( 3) / 8589934591 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -256 / C DATA IMACH( 7) / 255 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) / -256 / C DATA IMACH(10) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -50 / C DATA IMACH(10) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -32754 / C DATA IMACH(10) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C 60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT C ARITHMETIC (NOS OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 48 / C DATA IMACH( 3) / 281474976710655 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -974 / C DATA IMACH( 7) / 1070 / C DATA IMACH( 8) / 95 / C DATA IMACH( 9) / -926 / C DATA IMACH(10) / 1070 / C C MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT C ARITHMETIC (NOS/VE OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -4096 / C DATA IMACH( 7) / 4095 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -4096 / C DATA IMACH(10) / 4095 / C C MACHINE CONSTANTS FOR THE CRAY 1. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -8192 / C DATA IMACH( 7) / 8191 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -8192 / C DATA IMACH(10) / 8191 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 23 / C DATA IMACH( 3) / 8388607 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 38 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 55 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 9000. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -126 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE AMDAHL 470/V6, THE ICL 2900, THE ITEL AS/6, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN, C RM FORTRAN, PROFESSIONAL FORTRAN, AND LAHEY FORTRAN. C DATA IMACH( 1) / 2 / DATA IMACH( 2) / 31 / DATA IMACH( 3) / 2147483647 / DATA IMACH( 4) / 2 / DATA IMACH( 5) / 24 / DATA IMACH( 6) / -125 / DATA IMACH( 7) / 128 / DATA IMACH( 8) / 53 / DATA IMACH( 9) / -1021 / DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE MICROVAX - VMS FORTRAN. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 54 / C DATA IMACH( 9) / -101 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 62 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE SUN 3. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) /-1024 / C DATA IMACH(10) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C IPMPAR = IMACH(I) RETURN END REAL FUNCTION SPMPAR (I) C Arguments: INTEGER I C----------------------------------------------------------------------- C C SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR C THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT C I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE C SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND C ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN C C SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, C C SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, C C SPMPAR(3) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C----------------------------------------------------------------------- C WRITTEN BY C ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN VIRGINIA C C This routine obtained from: C Naval Surface Warfare Center Library of Mathematics Subroutines C (January 1993) C by Alfred H. Morris, Jr. C Naval Surface Warfare Center C Dahlgren Division C Dahlgren VA 22448-5000 C C----------------------------------------------------------------------- INTEGER EMAX,EMIN,IBETA,IPMPAR,M REAL B,BINV,BM1,ONE,W,Z C IF (I .GT. 1) GO TO 10 B = IPMPAR(4) M = IPMPAR(5) SPMPAR = B**(1 - M) RETURN C 10 IF (I .GT. 2) GO TO 20 B = IPMPAR(4) EMIN = IPMPAR(6) ONE = FLOAT(1) BINV = ONE/B W = B**(EMIN + 2) SPMPAR = ((W * BINV) * BINV) * BINV RETURN C 20 IBETA = IPMPAR(4) M = IPMPAR(5) EMAX = IPMPAR(7) C B = IBETA BM1 = IBETA - 1 ONE = FLOAT(1) Z = B**(M - 1) W = ((Z - ONE)*B + BM1)/(B*Z) C Z = B**(EMAX - 2) SPMPAR = ((W * Z) * B) * B RETURN END SUBROUTINE ERRMSG(CSTATS,CSUBRT,CMSGNM) IMPLICIT NONE C Arguments: CHARACTER CMSGNM*(*),CSTATS*(*),CSUBRT*(*) C C Local variables: INTEGER IOERR C C*********************************************************************** C Given: C CSTATS = 'WARNING' or 'FATAL' C CSUBRT = name of subroutine C CMESGN = message C Prints a warning message in a standardized way, C and STOPs if CSTATS='FATAL' C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C History: C 96.11.14 (PJF) Remove "getset" and hardwire "ioerr" C*********************************************************************** DATA IOERR/6/ IF(CSTATS.EQ.'FATAL')THEN WRITE (IOERR,FMT=9000) CSUBRT 9000 FORMAT(/,' >>>>> FATAL ERROR IN PROCEDURE: ',A) WRITE(IOERR,FMT=9010)CMSGNM 9010 FORMAT(1X,A) WRITE(IOERR,FMT=9020) 9020 FORMAT(' >>>>> EXECUTION ABORTED ') STOP ELSEIF(CSTATS.EQ.'WARNING')THEN WRITE(IOERR,FMT=9030)CSUBRT 9030 FORMAT(/,' >>>>> WARNING IN PROCEDURE: ',A) WRITE(IOERR,FMT=9010)CMSGNM 9040 FORMAT(1X,A) ENDIF RETURN END