|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 92928 (0x16b00)
Types: TextFile
Names: »taldiscr«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦this⟧ »taldiscr«
(mode list.yes
discrlist=set 400
discrete=set 200
scope day discrete discrlist
head cpu
if listing.yes
o discrlist
discrete=fortran cardmode.yes cond.yes names.yes,
index.no trunc.yes
o c
head cpu
mode list.no
finisb)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
program discrete
c version 1980-11-13
c common adjusted for logical variables
c double presicion removed
c character handling changed
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C-----------------------------------------------------------------------
C DISCRETE - VERSION 1A (MARCH 1976). FOR THE AUTOMATIC ANALYSIS OF
C MULTICOMPONENT EXPONENTIAL DECAY DATA FOR UP TO 5 COMPONENTS.
C MUST BE RUN ON A COMPUTER WITH AT LEAST 7 SIGNIFICANT FIGURES IN
C ITS SINGLE PRECISION FORTRAN ARITHMETIC. OTHERWISE THE DOUBLE
C PRECISION VERSION 1B SHOULD BE USED.
C REFERENCES - S. W. PROVENCHER (1976) WRITE UP OF DISCRETE,
C (1976) IN PREPARATION,
C (1976) J. CHEM. PHYS., VOL. 64, 2772,
C (1976) BIOPHYS. J., VOL. 16, 27-41.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - BLOCKDATA, DATAIO, WEIGHT, FANLYZ, YANLYZ
C WHICH IN TURN CALL - LSTSQR, EVAR, VARF, ETHEOR, PIVOT, PIVOT1,
C ANLERR, FISHER, RESIDU, PLPRES
C-----------------------------------------------------------------------
c----------------------------------------------------------------------
c change of input to read beam foil data in free format
c----------------------------------------------------------------------
external in,out; zone in,out
external freadb,freadr,ifreadi,freadra; logical freadb
LOGICAL LAST, REGINT, LBLOKA
COMMON /BLOK/ RBLOKA(3535), IBLOKA(169), NMAX, NINTMX, IBLOKB(9),
a ibloka 149->169
1 REGINT, LBLOKA(38)
A ADJUSTED FROM 29 TO 38
C
C***********************************************************************
C***********************************************************************
C
C THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES
C THAT YOU MAY HAVE TO MAKE IN THE MAIN PROGRAM. (SEE ALSO THE
C BLOCKDATA SUBPROGRAM.)
C
C***********************************************************************
C IF YOU CHANGE NMAX IN THE BLOCKDATA SUBPROGRAM, YOU MUST READJUST
C THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NMAX -
C
DIMENSION T(333), Y(333), SQRTW(333), YLYFIT(333)
C
C***********************************************************************
C ORDINARILY E SHOULD BE DIMENSIONED E(NMAX,6) IN THE FOLLOWING
C STATEMENT. HOWEVER, IF YOU WILL ALWAYS INPUT REGINT=.TRUE., AND
C IF NMAX IS GREATER THAN 333, YOU CAN SAVE STORAGE BY DIMENSIONING
C IT AS E(1,6) -
C
DIMENSION E(333,6)
C
C***********************************************************************
C ORDINARILY GSE SHOULD BE DIMENSIONED GSE(500,4). HOWEVER, IF YOU
C WILL ALWAYS INPUT REGINT=.FALSE., AND IF NMAX IS LESS THAN 333,
C YOU CAN SAVE A LITTLE STORAGE BY DIMENSIONING IT AS GSE(1,4) IN
C THE FOLLOWING STATEMENT -
C
DIMENSION GSE(500,4)
C
C***********************************************************************
C IF YOU CHANGE NINTMX IN THE BLOCKDATA SUBPROGRAM, YOU MUST
C READJUST THE DIMENSIONS IN THE FOLLOWING STATEMENT TO NINTMX -
C
DIMENSION TSTART(10), NT(10), DELTAT(10)
C
C THIS IS THE END OF ALL POSSIBLE CHANGES YOU MAY HAVE TO MAKE IN
C THE MAIN PROGRAM.
C
C***********************************************************************
C***********************************************************************
C
EQUIVALENCE (E(1,1),GSE(1,1))
C-----------------------------------------------------------------------
C DATAIO IS FOR INPUT AND OUTPUT OF INPUT DATA
100 CALL DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST)
IF (REGINT) GO TO 110
NDIME=NMAX
NDIMG=1
GO TO 200
110 NDIME=1
NDIMG=500
C-----------------------------------------------------------------------
C WEIGHT GENERATES CRUDE STARTING VALUES FOR CONSTRAINED LEAST SQUARES
C FITS TO RAW DATA TO GENERATE WEIGHTS AND DOES INITIAL
C COMPUTATIONS OF QUANTITIES FOR LATER USE.
200 CALL WEIGHT (.FALSE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
C-----------------------------------------------------------------------
C FANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF
C THE RAW DATA AND THE TRANSFORMS TO GET STARTING ESTIMATES
C FOR THE FINAL LEAST SQUARES ANALYSES OF THE RAW DATA.
CALL FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
C-----------------------------------------------------------------------
C WEIGHT USES STARTING VALUES FROM FANLYZ FOR CONSTRAINED LEAST
C SQUARES FITS TO RAW DATA TO GENERATE WEIGHTS FOR THE FINAL
C LEAST SQUARES FITS TO THE RAW DATA.
CALL WEIGHT (.TRUE.,Y,T,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
C-----------------------------------------------------------------------
C YANLYZ IS FOR THE CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS OF
C THE RAW DATA (USING THE STARTING VALUES OBTAINED FROM FANLYZ
C AND THE WEIGHTS FROM WEIGHT).
CALL YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
C-----------------------------------------------------------------------
IF (.NOT.LAST) GO TO 100
STOP
END
SUBROUTINE BLOCKDATA
LOGICAL REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT, PLOTRS,
1 NONNEG, FFAIL, PIVALP, PRFINL
REAL LAMST, LAMF, LAMNMX
LONG LABEL
COMMON /BLOK/ TRBASE(11), GF(11,126), DGFL(11,126), DGRIDF,
1 GFSTL, FHAT(11), DELTA, VARMIN, FHATMX, SUMWTY,
2 SUMWT, ZERO(6), CONVRG, LAMST(5), LAMF(5,5,2), DGRIDE, GESTL,
3 VAR, LAMNMX(2,2), PALPHA(6), PLAM(6), PRECIS, SIGYY, ENPHI,
4 SE(6,6), SIGMAP(11), PCERR(11), VGRID2(252), PLMTRY(6),
5 DELTAP(10), DFCAPL(6,5), SFDE(5), ADUB(11,11), DDSE(5,5)
COMMON /BLOK/ DTOLER(11), F(11,6), DF(11,6), JGEADD, NBINOM(20,5),
1 NVGRID(5), NPERG2(5), MDIMA, IBASE, IWTSGN, MINTER, MLAMMX,
2 NGRID2(5), NLAMMX, IWT, LABEL(20), NINT, NF, N, IWRITE, IREAD,
3 LINEPG, NMAX, NINTMX, MTRY, ITER, JLAMP1, NFTRY(5), NSIGFG,
4 REGINT, FAILED, PRPREL, NOBASE, WTED, REPEAT, PLOTRS, NONNEG,
5 FFAIL(5,2), PIVALP(11), PRFINL
C
C***********************************************************************
C***********************************************************************
C
C THE FOLLOWING INSTRUCTIONS DESCRIBE ALL POSSIBLE NECESSARY CHANGES
C THAT YOU MAY HAVE TO MAKE IN THE BLOCKDATA SUBPROGRAM. (SEE THE
C MAIN PROGRAM ALSO.)
C
C***********************************************************************
C ALL STATEMENTS FOR READING CARDS AND FOR PRINTING OUTPUT ARE OF
C THE FORM READ (IREAD,... AND WRITE(IWRITE,.... THUS, IREAD AND
C IWRITE MAY HAVE TO BE CHANGED TO THE VALUES APPROPRIATE FOR YOUR
C COMPUTER CENTER IN THE FOLLOWING STATEMENT -
C
DATA IREAD/5/,IWRITE/6/
C
C***********************************************************************
C IN THE FOLLOWING STATEMENT, NSIGFG MUST BE SET TO THE NUMBER OF
C SIGNIFICANT FIGURES CARRIED IN THE SINGLE PRECISION FLOATING POINT
C ARITHMETIC OF YOUR COMPUTER (E.G., 7 FOR IBM 360, 8 FOR
C UNIVAC 1108, ETC.) -
C
DATA NSIGFG/11/
C
C***********************************************************************
C IN THE FOLLOWING STATEMENT, LINEPG MUST BE SET TO THE NUMBER OF
C LINES PER PAGE AVAILABLE FOR PRINTED OUTPUT ON YOUR PRINTER -
C
DATA LINEPG/66/
C
C***********************************************************************
C NMAX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM NUMBER OF DATA
C POINTS YOU WILL EVER USE. CORE STORAGE CAN BE SAVED BY KEEPING
C NMAX AS LITTLE ABOVE THIS MAXIMUM AS POSSIBLE. IF YOU CHANGE NMAX
C IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE DIMENSIONS AS
C DIRECTED IN THE MAIN PROGRAM -
C
DATA NMAX/333/
C
C***********************************************************************
C NINTMX MUST BE GREATER THAN OR EQUAL TO THE MAXIMUM VALUE OF NINT
C YOU WILL EVER USE WITH REGINT=.TRUE.. IF YOU CHANGE
C NINTMX IN THE FOLLOWING STATEMENT, YOU MUST ALSO RESET THE
C DIMENSIONS AS DIRECTED IN THE MAIN PROGRAM -
C
DATA NINTMX /10/
C
C THIS IS THE END OF ALL POSSIBLE NECESSARY CHANGES YOU MAY HAVE TO
C MAKE IN THE BLOCK DATA SUBPROGRAM
C
C***********************************************************************
C***********************************************************************
C
DATA NGRID2/0, 20, 12, 10, 10/,MINTER/7/,
1 CONVRG/5.E-5/,MLAMMX/5/, SIGMAP/11*0./
END
C-----------------------------------------------------------------------
C SUBROUTINE DATAIO. FOR INPUT AND OUTPUT OF INITIAL INPUT DATA
C-----------------------------------------------------------------------
C CALLS NO EXTERNAL SUBPROGRAMS.
C-----------------------------------------------------------------------
SUBROUTINE DATAIO (T,Y,SQRTW,NMAX,TSTART,DELTAT,NT,NINTMX,LAST)
external in,out; zone in,out
external freadb,freadr,ifreadi,freadra; logical freadb
LOGICAL REGINT, LAST, NONNEG, NOBASE, PRY, PRPREL, PLOTRS, REPEAT,
1 LBLOKA, LBLOKB, LBLOKC, ABORT, PRFINL
LONG LA(9), LATF(2),LABEL
COMMON /BLOK/ RBLOKA(3535), IBLOKA(113), IWTSGN, IBLOKB, MLAMMX,
1 IBLOKC(5), NLAMMX, IWT, LABEL(20), NINT, IBLOKD, N, IWRITE,
2 IREAD, IBLOKE(3), MTRY, IBLOKF(8), REGINT, LBLOKA, PRPREL,
3 NOBASE, LBLOKB, REPEAT, PLOTRS, NONNEG, LBLOKC(21), PRFINL
DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), TSTART(NINTMX),
1 NT(NINTMX), DELTAT(NINTMX)
LATF(1)='T '
LATF(2)='F '
5010 FORMAT (12A6)
READ (IREAD,5010) LABEL
5020 FORMAT (35H1DISCRETE - VERSION 1A (MARCH 1976))
WRITE (IWRITE,5020)
5222 FORMAT(12A6)
WRITE(IWRITE,5222) LABEL
ABORT=.FALSE.
c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c5100 FORMAT (9A1,3I3)
c READ (IREAD,5100) LA,NLAMMX,IWT,MTRY
read(iread,5100) la
5100 format(9a1)
nlammx=ifreadi(in)
iwt =ifreadi(in)
mtry =ifreadi(in)
c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
DO 101 J=1,9
IF (LA(J).EQ.LATF(1) .OR. LA(J).EQ.LATF(2)) GO TO 101
5101 FORMAT (/43H THE FOLLOWING CARD HAS NO T OR F IN COLUMN,I2)
WRITE (IWRITE,5101) J
ABORT=.TRUE.
101 CONTINUE
IF (NLAMMX.GE.1 .AND. NLAMMX.LE.MLAMMX) GO TO 102
5102 FORMAT (/9H NLAMMX =,I4,21H IS NOT BETWEEN 1 AND,I2)
WRITE (IWRITE,5102) NLAMMX, MLAMMX
ABORT=.TRUE.
102 IF (IWT.EQ.1 .OR. IABS(IWT).EQ.2 .OR. IABS(IWT).EQ.3
1 .OR. IWT.EQ.4) GO TO 103
5103 FORMAT (/1X,I4,29H IS AN IMPROPER VALUE FOR IWT)
WRITE (IWRITE,5103) IWT
ABORT=.TRUE.
103 IF (MTRY.GE.1 .AND. MTRY.LE.45) GO TO 104
5104 FORMAT (/7H MTRY =,I4,24H IS NOT BETWEEN 1 AND 45)
WRITE (IWRITE,5104) MTRY
ABORT=.TRUE.
104 IF (.NOT.ABORT) GO TO 105
5105 FORMAT (/21H ERRONEOUS CARD IMAGE,10X,9A1,3I3)
WRITE (IWRITE,5105) LA,NLAMMX,IWT,MTRY
STOP
105 LAST=LA(1) .EQ. LATF(1)
REGINT=LA(2) .EQ. LATF(1)
NOBASE=LA(3) .EQ. LATF(1)
NONNEG=LA(4) .EQ. LATF(1)
PRY=LA(5) .EQ. LATF(1)
PRPREL=LA(6) .EQ. LATF(1)
PRFINL=LA(7) .EQ. LATF(1)
PLOTRS=LA(8) .EQ. LATF(1)
REPEAT=LA(9) .EQ. LATF(1)
5106 FORMAT (/7H LAST =,L2,1X,8HREGINT =,L2,2X,8HNOBASE =,L2,2X,
1 8HNONNEG =,L2,2X,5HPRY =,L2/9H PRPREL =,L2,2X,8HPRFINL =,L2,2X,
2 8HPLOTRS =,L2,2X,8HREPEAT =,L2/9H NLAMMX =,I2,1X,8HIWT =,
3 I2,2X,8HMTRY =,I3)
WRITE (IWRITE,5106) LAST,REGINT,NOBASE,NONNEG,PRY,PRPREL,PRFINL,
1 PLOTRS,REPEAT,NLAMMX,IWT,MTRY
IWTSGN=IWT
IWT=IABS(IWT)
IF (REGINT) GO TO 150
C-----------------------------------------------------------------------
C READS IN T VALUES (REGINT=.FALSE.).
C-----------------------------------------------------------------------
c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c READ (IREAD,5110) N
5110 FORMAT (I5)
c READ (IREAD,5120) (T(J),J=1,N)
5120 FORMAT (5E15.6)
n=ifreadi(in);
call freadra(in,t,n)
c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
NINT=1
TSTART(1)=0.
DELTAT(1)=0.
NT(1)=N
GO TO 200
C-----------------------------------------------------------------------
C CALCULATES T VALUES FROM TSTART, TEND, AND NT (REGINT=.TRUE.).
C-----------------------------------------------------------------------
c 150 READ (IREAD,5110) NINT
150 nint=ifreadi(in)
IF (NINT.LE.NINTMX .AND. NINT.GE.1) GO TO 160
5160 FORMAT (///7H NINT =,I5,30H IS NOT BETWEEN 1 AND NINTMX =,I2)
WRITE (IWRITE,5160) NINT,NINTMX
STOP
5165 FORMAT (///9X,6HTSTART,11X,4HTEND,3X,2HNT)
160 WRITE (IWRITE,5165)
N=0
DO 170 J=1,NINT
c READ (IREAD,5170) TSTART(J),TEND,NT(J)
tstart(j)=freadr(in); tend=freadr(in); nt(j)=ifreadi(in)
5170 FORMAT (2E15.6,I5)
WRITE (IWRITE,5170) TSTART(J),TEND,NT(J)
IF (NT(J) .GE. 2) GO TO 175
5175 FORMAT (/24H ABOVE NT IS LESS THAN 2)
WRITE (IWRITE,5175)
STOP
175 DUM=(TEND-TSTART(J))/FLOAT(NT(J)-1)
N=N+1
T(N)=TSTART(J)
L=NT(J)
DELTAT(J)=DUM
DO 180 K=2,L
N=N+1
180 T(N)=T(N-1)+DUM
170 CONTINUE
C-----------------------------------------------------------------------
C READS IN Y VALUES.
C-----------------------------------------------------------------------
200 IF (N.GT.2*NLAMMX+1 .AND. N.LE.NMAX) GO TO 205
5205 FORMAT (///4H N =,I5,37H IS NOT BETWEEN 2*NLAMMX+1 AND NMAX =,I5)
WRITE (IWRITE,5205) N,NMAX
STOP
205 IF (PRY) WRITE (IWRITE,5004)
5004 FORMAT (///1X)
c READ (IREAD,5120) (Y(J),J=1,N)
call freadra(in,y,n)
210 IF (IWT .EQ. 4) GO TO 220
IF (.NOT.PRY) GO TO 500
5210 FORMAT (2(16X,1HT,12X,1HY)/)
WRITE (IWRITE,5210)
5220 FORMAT (2(1PE17.5,E13.5))
WRITE (IWRITE,5220) (T(J),Y(J),J=1,N)
GO TO 500
C-----------------------------------------------------------------------
C READS IN SPECIAL WEIGHTS (IWT=4).
C-----------------------------------------------------------------------
c 220 READ (IREAD,5120) (SQRTW(J),J=1,N)
220 call freadra(in,sqrtw,n)
DO 227 J=1,N
227 SQRTW(J)=SQRT(SQRTW(J))
IF (.NOT.PRY) GO TO 500
5230 FORMAT (2(17X,1HT,12X,1HY,5X,8HSQRT(WT))/)
WRITE (IWRITE,5230)
5240 FORMAT (2(1PE18.5,2E13.5))
WRITE (IWRITE,5240) (T(J),Y(J),SQRTW(J),J=1,N)
500 RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE WEIGHT. DOES INITIAL LEAST SQUARES FITS TO RAW DATA TO
C GENERATE WEIGHTS, EVALUATES DATA TRANSFORMS, TRANSFORM AND
C EXPONENTIAL SUMS AT GRID POINTS FOR INTERPOLATION, AND OTHER
C QUANTITIES FOR LATER USE.
C IF FINAL=.TRUE. USES LAMF (VALUES FROM FANLYZ).
C IF FINAL=.FALSE., GENERATES ITS OWN CRUDE STARTING LAMBDAS.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - LSTSQR, ETHEOR
C WHICH IN TURN CALL - EVAR, PIVOT, PIVOT1, ANLERR
C-----------------------------------------------------------------------
SUBROUTINE WEIGHT (FINAL,Y,T,SQRTW,TIMU,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
EXTERNAL EVAR
LOGICAL FINAL, REGINT, NOBASE, WTED, FAILED, LDUM, LBLOKA, PRPREL,
1 LBLOKB, FFAIL
REAL LAMNMX, LAMST, LAMF, MU
LONG ALPH(3,2), LABEL
COMMON /BLOK/ TRBASE(11), GF(11,126), DGFL(11,126), DGRIDF,
1 GFSTL, FHAT(11), DELTA, VARMIN, FHATMX, SUMWTY,
2 SUMWT, ZERO(6), CONVRG, LAMST(5), LAMF(5,5,2), DGRIDE, GESTL,
3 VAR, LAMNMX(2,2), PALPHA(6), PLAM(6), PRECIS, RBLOKA(652)
COMMON /BLOK/ JGEADD, NBINOM(20,5), NVGRID(5),
1 NPERG2(5), MDIMA, IBASE, IWTSGN, MINTER, MLAMMX, NGRID2(5),
2 NLAMMX, IWT, LABEL(20), NINT, NF, N, IWRITE, IBLOKA(12), NSIGFG,
3 REGINT, FAILED, PRPREL, NOBASE, WTED, LBLOKA(6), FFAIL(5,2),
4 LBLOKB(14)
a lbloka 3->6, lblokb 12->14
DIMENSION E(NDIME,6), WTALP(6), WTLAM(6), MU(11)
DIMENSION Y(NMAX), T(NMAX), SQRTW(NMAX), TIMU(NMAX),
1 TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX), GSE(NDIMG,4),
2 TNLAMN(2), GAMMA(21)
C++++++++++
ALPH(1,1)='TRAN'; ALPH(2,1)='SFOR'; ALPH(3,2)='MS'
ALPH(1,2)='RAW' ; ALPH(2,2)='DATA'; ALPH(3,2)=' '
TNLAMN(1)=.2 ; TNLAMN(2)=.02
T1LAMX=2.08
NGRID1=240
NG2MAX=20
C++++++++++
IF (FINAL) GO TO 120
FHATMX=0.
PRECIS=EXP(-2.303*(FLOAT(NSIGFG)-2.))
PRECIS=AMAX1(1.E-7,PRECIS)
MDIMA=2*MLAMMX+1
K=MLAMMX+1
DO 105 J=1,K
105 ZERO(J)=0.
TN=-1.E+20
DDUM=-1.E+20
DO 106 K=1,5
DUM=1.E+20
DO 107 J=1,N
IF (K .EQ. 1) TN=AMAX1(TN,T(J))
107 IF (T(J).GT.DDUM .AND. T(J).LT.DUM) DUM=T(J)
DDUM=DUM+1.E-5*ABS(DUM)
IF (K .EQ. 1) T1=DUM
IF (K .EQ. 2) T2=DUM
106 CONTINUE
T5=DUM
TNL1=-1.E+20
DUM=TN-1.E-5*ABS(TN)
DO 108 J=1,N
108 IF (T(J) .LT. DUM) TNL1=AMAX1(TNL1,T(J))
DUM=.25*(T5-T1)
DELTA=DUM-T1
C-----------------------------------------------------------------------
C LAMNMX(J,K) = MIN. (J=1) AND MAX. (J=2) ALLOWED VALUES OF LAMBDA
C FOR LEAST SQUARES PARAMETERS (K=1) AND FOR STARTING VALUES (K=2).
C-----------------------------------------------------------------------
LAMNMX(2,1)=T1LAMX/(T1+DELTA)
LAMNMX(2,2)=AMIN1(LAMNMX(2,1),.693/(T1+DELTA))
LAMNMX(1,2)=TNLAMN(1)/(TN+DELTA)
IF (NOBASE) GO TO 110
NPARMX=2*NLAMMX+1
IBASE=1
LAMNMX(1,1)=LAMNMX(1,2)
GO TO 115
110 IBASE=0
NPARMX=2*NLAMMX
FHAT(NPARMX+1)=0.
LAMNMX(1,1)=TNLAMN(2)/(TN+DELTA)
C-----------------------------------------------------------------------
C IF NECESSARY, INCREASE DELTA SO THAT DISTANCE BETWEEN FIRST TWO
C POINTS ON X=ALOG(T+DELTA) AXIS IS LESS THAN A HALF-PERIOD OF
C THE HIGHEST FREQUENCY USED IN THE TRANSFORMS
C-----------------------------------------------------------------------
115 XRANGE=ALOG((TN+DELTA)**2/((T1+DELTA)*(TNL1+DELTA)))
IF (ALOG((T2+DELTA)/(T1+DELTA))*FLOAT(2*NLAMMX) .LT. XRANGE)
1 GO TO 150
DELTA=DELTA+DUM
GO TO 115
120 IWT=IABS(IWTSGN)
IF (IWT.EQ.1 .OR. IWT.EQ.4) RETURN
DO 130 J=1,N
130 Y(J)=Y(J)/SQRTW(J)
150 IF (IWT .EQ. 4) GO TO 400
SUMWT=FLOAT(N)
SUMWTY=0.
VARMIN=0.
DO 160 J=1,N
SUMWTY=SUMWTY+Y(J)
VARMIN=VARMIN+Y(J)*Y(J)
160 SQRTW(J)=1.
VARMIN=VARMIN*CONVRG**2
IF (IWT.EQ.1) GO TO 400
WTED=.FALSE.
VARBES=1.E+20
RATIOL=LAMNMX(2,2)/LAMNMX(1,2)
C-----------------------------------------------------------------------
C START OF MAIN LOOP FOR STEPWISE ANALYSIS
C-----------------------------------------------------------------------
DO 200 JLAM=1,NLAMMX
NF=2*JLAM+IBASE
5200 FORMAT (1H1,12A6//5X,
1 41HPRELIMINARY ANALYSIS TO DETERMINE WEIGHTS/5X,
2 17HANALYSIS ASSUMING,I2,11H COMPONENTS)
IF (PRPREL) WRITE (IWRITE,5200) LABEL,JLAM
IF (FINAL) GO TO 220
FFAIL(JLAM,1)=.TRUE.
DUM=RATIOL**(1./FLOAT(JLAM+1))
LAMST(1)=LAMNMX(1,2)*DUM
IF (JLAM .EQ. 1) GO TO 240
DO 210 J=2,JLAM
210 LAMST(J)=LAMST(J-1)*DUM
GO TO 240
220 DO 230 J=1,JLAM
230 LAMST(J)=LAMF(J,JLAM,1)
240 CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,GSE,NDIMG,
1 SQRTW,TIMU,E,NMAX,TSTART,DELTAT,NT,NINT,1,FFAIL(JLAM,1),NDIME,
2 .FALSE.,PRPREL)
IF (FAILED .OR. VAR.GE.VARBES) GO TO 300
VARBES=VAR
NLINWT=JLAM+IBASE
DO 260 J=1,NLINWT
WTALP(J)=PALPHA(J)
260 WTLAM(J)=PLAM(J)
IF (IWTSGN .LT. 0) NLINWT=JLAM
200 CONTINUE
300 IF (VARBES .LT. 1.E+20) GO TO 310
5300 FORMAT (////34H 1-COMPONENT ANALYSIS TO DETERMINE,
1 54H WEIGHTS SOMEHOW FAILED - UNIT WEIGHTS WILL BE USED. ,
2 11(2H**))
WRITE (IWRITE,5300)
IWT=1
GO TO 400
C-----------------------------------------------------------------------
C CALCULATE WEIGHTS INCLUDING ERRFIT (STANDARD ERROR IN REGION OF
C THE THEORETICAL CURVE WHERE IT IS CLOSEST TO Y=0., WHERE A
C WEIGHT COULD OTHERWISE BE DISASTROUSLY LARGE IF ERRFIT WERE
C NOT ACCOUNTED FOR).
C-----------------------------------------------------------------------
310 DUM=1.E+20
DO 320 K=1,N
S=0.
DDUM=T(K)
DO 330 J=1,NLINWT
330 S=S+WTALP(J)*EXP(-DDUM*WTLAM(J))
TIMU(K)=S
IF (ABS(S) .GE. DUM) GO TO 320
DUM=ABS(S)
L=K
320 CONTINUE
K=MIN0(N,L+5)
I=MAX0(1,K-9)
DUM=0.
DDUM=0.
IF (IWTSGN.LT.0 .AND. .NOT.NOBASE) DDUM=WTALP(NLINWT+1)
DO 340 J=I,K
340 DUM=DUM+(Y(J)-TIMU(J)-DDUM)**2
ERRFIT=SQRT(DUM/FLOAT(K-I+1))
DO 350 K=1,N
SQRTW(K)=1./(ABS(TIMU(K))+ERRFIT)
IF (IWT .EQ. 2) SQRTW(K)=SQRT(SQRTW(K))
350 CONTINUE
L=1
IF (FINAL) L=2
5350 FORMAT (/////41H PARAMETERS USED TO GENERATE WEIGHTS FOR ,
1 3A4//5X,6HLAMBDA,9X,5HALPHA,
2 /49HERRFIT (UNCERTAINTY TERM ADDED TO ABSOLUTE VALUES,
3 24H OF THEORETICAL CURVE) =,1PE9.2/(1X,1PE10.3,E14.3))
WRITE (IWRITE,5350) (ALPH(J,L),J=1,3),ERRFIT,(WTLAM(K),WTALP(K),
1 K=1,NLINWT)
400 WTED=IWT.NE.1
IF (FINAL) GO TO 600
C-----------------------------------------------------------------------
C INITIAL COMPUTATION OF QUANTITIES TO BE USED LATER
C-----------------------------------------------------------------------
DGRIDE=ALOG(LAMNMX(2,2)/LAMNMX(1,2))/FLOAT(NGRID1-1)
RGRIDE=EXP(DGRIDE)
JGEADD=MINTER/2
IF (NOBASE) JGEADD=JGEADD+(ALOG(LAMNMX(1,2)/LAMNMX(1,1))/
1 DGRIDE+.5)
GEST=LAMNMX(1,2)/RGRIDE**(JGEADD+1)
GESTL=ALOG(GEST)
NGE1=NGRID1+JGEADD+(ALOG(2.*LAMNMX(2,1)/LAMNMX(2,2))/DGRIDE+
1 .5)+MINTER/2
IF (NGE1 .LE. 500) GO TO 401
5401 FORMAT (////20H FIX LAMNMX. NGE1 =,I9,1X,51(2H**))
WRITE (IWRITE,5401) NGE1
STOP
401 JGEADD=JGEADD+1
IF (NLAMMX .EQ. 1) GO TO 403
DO 402 J=2,NLAMMX
402 NPERG2(J)=NGRID1/NGRID2(J)
403 GAMMA(1)=1.
DO 404 J=1,NG2MAX
404 GAMMA(J+1)=GAMMA(J)*FLOAT(J)
DO 405 J=1,NLAMMX
DO 406 K=J,NG2MAX
L=K-J+1
406 NBINOM(K,J)=(GAMMA(K)/(GAMMA(L)*GAMMA(J))+.5)
IF (J .EQ. 1) GO TO 405
L=NGRID2(J)+1
K=L-J
NVGRID(J)=(GAMMA(L)/(GAMMA(J+1)*GAMMA(K))+.5)
405 CONTINUE
NGF=NGRID1/2+2*(MINTER/2)
DGRIDF=ALOG(LAMNMX(2,1)/LAMNMX(1,1))/FLOAT(NGRID1/2-1)
RGRIDF=EXP(DGRIDF)
GFST=LAMNMX(1,1)/RGRIDF**(MINTER/2+1)
GFSTL=ALOG(GFST)
C-----------------------------------------------------------------------
C COMPUTE FHAT (TRANSFORMS OF DATA) AND GF AND DGFL (INTERPOLATION
C GRID POINTS FOR TRANSFORMS AND THEIR DERIVATIVES W.R.T.
C ALOG(LAMBDA))
C-----------------------------------------------------------------------
IF (.NOT.REGINT) GO TO 420
DO 410 J=1,NINT
410 TSTART(J)=TSTART(J)+DELTA
420 DO 430 J=1,N
430 T(J)=T(J)+DELTA
MU(1)=6.283185/XRANGE
DO 440 K=1,NPARMX
IF (K .GT. 1) GO TO 450
DO 445 J=1,N
FHATMX=FHATMX+SQRTW(J)*ABS(Y(J))
445 TIMU(J)=SQRTW(J)
GO TO 470
450 L=K/2
IF (L+L .LT. K) GO TO 460
DUM=MU(1)
IF (L .GT. 1) DUM=MU(L-1)+DUM
MU(L)=DUM
DO 455 J=1,N
455 TIMU(J)=COS(ALOG(T(J))*DUM)*SQRTW(J)
GO TO 470
460 DUM=MU(L)
DO 465 J=1,N
465 TIMU(J)=SIN(ALOG(T(J))*DUM)*SQRTW(J)
470 DUM=0.
DDUM=0.
DO 475 J=1,N
DUM=DUM+Y(J)*TIMU(J)
475 DDUM=DDUM+TIMU(J)
FHAT(K)=DUM
TRBASE(K)=DDUM
DUM=GFST
DO 480 L=1,NGF
DUM=DUM*RGRIDF
S=0.
SS=0.
IF (REGINT) GO TO 500
DO 490 J=1,N
DDUM=EXP(-DUM*T(J))*TIMU(J)
S=S+DDUM
490 SS=SS-DDUM*T(J)
GO TO 550
C-----------------------------------------------------------------------
C USE RECURSION RELATIONS IF REGINT=.TRUE.
C-----------------------------------------------------------------------
500 J=0
DO 510 JJ=1,NINT
DDUM=EXP(-DUM*DELTAT(JJ))
DDDUM=EXP(-DUM*TSTART(JJ))
II=NT(JJ)
DO 520 I=1,II
J=J+1
SSS=DDDUM*TIMU(J)
S=S+SSS
SS=SS-T(J)*SSS
DDDUM=DDDUM*DDUM
520 CONTINUE
510 CONTINUE
550 GF(K,L)=S
DGFL(K,L)=SS*DUM
480 CONTINUE
440 CONTINUE
IF (PRPREL) WRITE (IWRITE,5560) DELTA,FHAT(1),(MU(K),FHAT(2*K),
1 FHAT(2*K+1),K=1,NLAMMX)
5560 FORMAT (////6X,2HMU,13X,
1 38HREAL AND IMAGINARY PARTS OF TRANSFORMS,15X,7HDELTA =,1PE9.2
2 /4X,4H.000,1PE17.3/(1X,0PF7.3,1PE17.3,E14.3))
DO 560 J=1,N
560 T(J)=T(J)-DELTA
IF (.NOT.REGINT) GO TO 600
DO 570 J=1,NINT
570 TSTART(J)=TSTART(J)-DELTA
600 IF (.NOT.WTED) GO TO 620
VARMIN=0.
SUMWT=0.
SUMWTY=0.
DO 610 J=1,N
Y(J)=Y(J)*SQRTW(J)
VARMIN=VARMIN+Y(J)*Y(J)
IF (NOBASE) GO TO 610
SUMWT=SUMWT+SQRTW(J)*SQRTW(J)
SUMWTY=SUMWTY+SQRTW(J)*Y(J)
610 CONTINUE
VARMIN=VARMIN*CONVRG**2
620 IF (FINAL .OR. .NOT.REGINT) GO TO 900
C-----------------------------------------------------------------------
C COMPUTE GSE (INTERPOLATION GRID POINTS FOR EXPONENTIAL SUMS)
C-----------------------------------------------------------------------
DUM=GEST
L=NGE1-(.69315/DGRIDE)
DO 630 K=1,NGE1
DUM=DUM*RGRIDE
LDUM=K .LE. L
CALL ETHEOR (DUM,.TRUE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,NT,
1 NINT,S,SS,SSS,DDDUM,.TRUE.,LDUM)
GSE(K,1)=S
GSE(K,2)=SS
GSE(K,3)=SSS
GSE(K,4)=DDDUM
630 CONTINUE
900 RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE FANLYZ. DOES CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS
C OF RAW DATA (USING GRID SEARCH IF NECESSARY) TO GET STARTING
C VALUES FOR ANALYSIS OF TRANSFORMS, WHICH YEILD LAMF (STARTING
C VALUES FOR FINAL ANALYSIS OF RAW DATA IN YANLYZ).
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - LSTSQR, EVAR
C WHICH IN TURN CALL - PIVOT, PIVOT1, ETHEOR, VARF
C-----------------------------------------------------------------------
SUBROUTINE FANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
EXTERNAL VARF, EVAR
LOGICAL PRPREL, WTED, FAILED, SEARCH, LDUM, FFAIL, LBLOKA,
1 LBLOKB, LBLOKC, LBLOKD, SAVEPK
REAL LAMNMX, LAMST, LAMF
LONG LABEL
COMMON /BLOK/ RBLOKA(2808), LAMST(5), LAMF(5,5,2), DGRIDE, GESTL,
1 VAR, LAMNMX(2,2), PALPHA(6), PLAM(6), RBLOKB(61), VGRID2(252),
2 RBLOKC(340), JGEADD, IBLOKA(100), NVGRID(5), NPERG2(5), IBLOKB,
3 IBASE, IBLOKC(2), MLAMMX, NGRID2(5), NLAMMX, IWT,
4 LABEL(20), NINT, NF, N, IWRITE, IBLOKD(4), MTRY, IBLOKE(2),
5 NFTRY(5), IBLOKF, LBLOKA, FAILED, PRPREL, LBLOKB, WTED
COMMON /BLOK/ LBLOKC(6), FFAIL(5,2), LBLOKD(14)
a lblokc 3->6, lblokd 12->14
DIMENSION E(NDIME,6), LENDEQ(5), LEQ(5), JCENTR(20), LSTART(5,45)
DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4),
1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX), NT(NINTMX)
EQUIVALENCE (LENDEQ(1),LEND1), (LENDEQ(2),LEND2), (LENDEQ(3),
1 LEND3), (LENDEQ(4),LEND4), (LENDEQ(5),LEND5),
2 (LEQ(1),L1), (LEQ(2),L2), (LEQ(3),L3), (LEQ(4),L4), (LEQ(5),L5)
C++++++++
LLPEAK=1
C++++++++
WTED=IWT .NE. 1
LAMST(1)=SQRT(LAMNMX(2,2)*LAMNMX(1,2))
C-----------------------------------------------------------------------
C START OF MAIN LOOP FOR STEPWISE ANALYSIS
C-----------------------------------------------------------------------
DO 200 JLAM=1,NLAMMX
C-----------------------------------------------------------------------
C INITIALIZE FOR FITS ASSUMING JLAM COMPONENTS
C-----------------------------------------------------------------------
5150 FORMAT (1H1,5X,12A6/)
IF (PRPREL) WRITE (IWRITE,5150) LABEL
VARE=1.E+20
VARTRY=1.E+20
NF=2*JLAM+IBASE
IF (JLAM .EQ. 1) GO TO 220
LVGRID=NVGRID(JLAM)
DO 210 J=1,LVGRID
210 VGRID2(J)=0.
220 JLAML1=JLAM-1
JLAMP1=JLAM+1
IF (PRPREL) WRITE (IWRITE,5210) JLAM
5210 FORMAT (5X,45HINITIAL ANALYSIS OF THE RAW DATA AND THEN THE,
2 22H TRANSFORMS - ASSUMING,I2,11H COMPONENTS)
C-----------------------------------------------------------------------
C START OF LOOP FOR UP TO MTRY JLAM-COMPONENT FITS
C-----------------------------------------------------------------------
DO 300 JTRY=1,MTRY
5310 FORMAT (////5X,22HFIT USING STARTING SET,I3)
IF (PRPREL) WRITE (IWRITE,5310) JTRY
IF (JLAM.EQ.1) GO TO 600
IF (JTRY-2) 315,350,500
C-----------------------------------------------------------------------
C SELECT FIRST SET OF STARTING VALUES
C-----------------------------------------------------------------------
315 DO 320 J=1,JLAML1
320 LAMST(J)=LAMF(J,JLAML1,2)
IF (LPEAK.EQ.1 .OR. LPEAK.EQ.JLAM) GO TO 345
LAMST(JLAM)=SQRT(LAMF(LPEAK-1,JLAML1,2)*LAMF(LPEAK,JLAML1,2))
GO TO 600
345 IF (LPEAK .EQ. 1) LAMST(JLAM)=SQRT(LAMNMX(1,2)*LAMF(1,JLAML1,2))
IF (LPEAK .EQ. JLAM) LAMST(JLAM)=SQRT(LAMF(JLAML1,JLAML1,2)*
1 LAMNMX(2,2))
GO TO 600
C-----------------------------------------------------------------------
C DO COARSE GRID SEARCH
C-----------------------------------------------------------------------
350 L=NGRID2(JLAM)-JLAM
DO 360 J=1,JLAM
L=L+1
360 LENDEQ(J)=L
IF (JLAM .EQ. MLAMMX) GO TO 368
DO 365 J=JLAMP1,MLAMMX
365 LENDEQ(J)=LENDEQ(JLAM)
368 JCENTR(1)=JGEADD-1+(NPERG2(JLAM)+1)/2
L=NGRID2(JLAM)
DO 370 J=2,L
370 JCENTR(J)=JCENTR(J-1)+NPERG2(JLAM)
SEARCH=.TRUE.
C-----------------------------------------------------------------------
C ENTRY POINT TO FIND STARTING LAMBDA VALUES FROM GRID
C-----------------------------------------------------------------------
400 LOC=0
DO 410 L1=1,LEND1
LST2=L1+1
DO 420 L2=LST2,LEND2
LST3=L2+1
IF (JLAM .LT. 3) LST3=LEND3
DO 430 L3=LST3,LEND3
LST4=L3+1
IF (JLAM .LT. 4) LST4=LEND4
DO 440 L4=LST4,LEND4
LST5=L4+1
IF (JLAM .LT. 5) LST5=LEND5
DO 450 L5=LST5,LEND5
LOC=LOC+1
IF (SEARCH) GO TO 498
IF (VGRID2(LOC) .GT. VGMAX) GO TO 496
MINDIF=1000
IF (JTRY .LE. 2) GO TO 493
L=JTRY-1
DO 491 J=2,L
I=0
DO 492 K=1,JLAM
492 I=I+IABS(LEQ(K)-LSTART(K,J))
MINDIF=MIN0(MINDIF,I)
491 CONTINUE
IF (MINDIF-MAXDIF) 496,493,494
493 IF (VGRID2(LOC) .GE. VGBEST) GO TO 496
494 VGBEST=VGRID2(LOC)
MAXDIF=MINDIF
DO 495 J=1,JLAM
495 LSTART(J,JTRY)=LEQ(J)
496 IF (LOC .LT. LVGRID) GO TO 450
DO 497 J=1,JLAM
L=LSTART(J,JTRY)
497 LAMST(J)=EXP(GESTL+FLOAT(JCENTR(L))*DGRIDE)
GO TO 600
498 DO 499 J=1,JLAM
L=LEQ(J)
499 PLAM(J)=EXP(GESTL+FLOAT(JCENTR(L))*DGRIDE)
CALL EVAR (0.,JLAM,DUM,T,SQRTW,NMAX,TSTART,DELTAT,NT,
1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,JLAM+IBASE,1.E+20,0,NDIME,.TRUE.,
2 GSE,NDIMG,J)
450 CONTINUE
440 CONTINUE
430 CONTINUE
420 CONTINUE
410 CONTINUE
C-----------------------------------------------------------------------
C FIND LOWEST VARIANCE ON COARSE GRID
C-----------------------------------------------------------------------
500 VGMAX=1.E+20
DO 510 J=1,LVGRID
510 VGMAX=AMIN1(VGMAX,VGRID2(J))
IF (JTRY .GT. 2) VGMAX=10.*VGMAX
MAXDIF=-1
SEARCH=.FALSE.
VGBEST=1.E+20
GO TO 400
C-----------------------------------------------------------------------
C DO LEAST SQUARES FITS
C-----------------------------------------------------------------------
600 LDUM=JTRY.EQ.1
CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,
1 EVAR,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,
2 2, LDUM ,NDIME,.TRUE.,PRPREL)
FAILED=FAILED .OR. VAR.GT.VARE
IF (VAR .GT. VARE) GO TO 300
VARE=VAR
DO 605 J=1,JLAM
605 LAMF(J,JLAM,1)=PLAM(J)
SAVEPK=VARTRY .GE. 1.E+20
IF (FAILED) GO TO 610
5600 FORMAT (/51H END OF FIT TO RAW DATA, START OF FIT TO TRANSFORMS)
IF (PRPREL) WRITE (IWRITE,5600)
CALL LSTSQR (PLAM,JLAM,.FALSE.,NF,Y,T,JLAM+IBASE,
1 VARF,GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,
2 3,.FALSE.,NDIME,.TRUE.,PRPREL)
IF (FAILED .AND. VAR.GE.VARTRY) GO TO 300
SAVEPK=.TRUE.
VARTRY=VAR
DO 608 J=1,JLAM
608 LAMF(J,JLAM,2)=PLAM(J)
610 IF (.NOT.SAVEPK .AND. JLAM.NE.1) GO TO 300
IF (JLAM .EQ. NLAMMX) GO TO 670
C-----------------------------------------------------------------------
C FINDS MAX SUM OF ABSOLUTE ADJACENT AMPLITUDES FOR LATER DETERMINATION
C OF STARTING VALUE FOR LAMST(JLAMP1)
C-----------------------------------------------------------------------
DUM=ABS(PALPHA(1))+ABS(PALPHA(JLAMP1))
LLPEAK=1
IF (DUM .GT. ABS(PALPHA(JLAM))) GO TO 620
LLPEAK=JLAMP1
DUM=ABS(PALPHA(JLAM))
620 IF (JLAM .EQ. 1) GO TO 690
DO 630 J=2,JLAM
DDUM=ABS(PALPHA(J-1))+ABS(PALPHA(J))
IF (DDUM .LE. DUM) GO TO 630
DUM=DDUM
LLPEAK=J
630 CONTINUE
670 IF (.NOT.FAILED .OR. JLAM.EQ.1) GO TO 690
300 CONTINUE
JTRY=MTRY
690 LPEAK=LLPEAK
NFTRY(JLAM)=JTRY
FFAIL(JLAM,2)=FAILED
FFAIL(JLAM,1)=VARTRY .GE. 1.E+20
IF (.NOT.FFAIL(JLAM,1)) GO TO 200
DO 695 J=1,JLAM
695 LAMF(J,JLAM,2)=LAMF(J,JLAM,1)
200 CONTINUE
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE YANLYZ. FOR CONSTRAINED STEPWISE LEAST SQUARES ANALYSIS
C OF RAW DATA USING LAMF (STARTING VALUES OBTAINED FROM FANLYZ).
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - LSTSQR, FISHER, RESIDU
C WHICH IN TURN CALL - EVAR, ETHEOR, PIVOT, PIVOT1, PLPRES
C-----------------------------------------------------------------------
SUBROUTINE YANLYZ (T,Y,SQRTW,YLYFIT,E,NMAX,NDIME,
1 TSTART,DELTAT,NT,NINTMX,GSE,NDIMG)
EXTERNAL EVAR
LOGICAL NOBASE, PRFINL, PLOTRS, REPEAT, WTED, FAILED, YFAIL, DONE,
1 REGINT, FFAIL, LBLOKA, LBLOKB, LBLOKC
REAL LAMNMX, LAMST, LAMF
LONG NEXBES(2,5), AFFAIL(9,2), LABEL
COMMON /BLOK/ RBLOKA(2808), LAMST(5), LAMF(5,5,2), RBLOKB(2),
1 VAR, LAMNMX(2,2), PALPHA(6), PLAM(6), RBLOKC, SIGYY, ENPHI,
2 RBLOKD(36), SIGMAP(11), PCERR(11), RBLOKE(592), IBLOKA(112),
3 IBASE, IBLOKB(8), NLAMMX, IWT, LABEL(20), NINT, NF, N,
4 IWRITE, IBLOKC, LINEPG, IBLOKD(3), ITER, JLAMP1, NFTRY(5),
5 IBLOKE, REGINT, FAILED, LBLOKA, NOBASE, WTED, REPEAT, PLOTRS
COMMON /BLOK/ LBLOKB, FFAIL(5,2), LBLOKC(11), PRFINL
DIMENSION E(NDIME,6), YFAIL(5), ASAVE(6,5,6),
1 PROBSV(5), VARSV(5), SIGYSV(5), ITERSV(5),
2 ENPHSV(5), PROBLN(5), DONE(5)
DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), GSE(NDIMG,4),
1 YLYFIT(NMAX), TSTART(NINTMX), DELTAT(NINTMX),
2 NT(NINTMX), PUNCOR(6)
C++++++++
NEXBES(1,1)=' ' ; NEXBES(2,1)=' '
NEXBES(1,2)='SECO'; NEXBES(2,2)='ND'
NEXBES(1,3)=' THI'; NEXBES(2,3)='RD'
NEXBES(1,4)='FOUR'; NEXBES(2,4)='TH'
NEXBES(1,5)=' FIF'; NEXBES(2,5)='TH'
AFFAIL(1,1)=' '; AFFAIL(2,1)=' '; AFFAIL(3,1)=' '
AFFAIL(4,1)=' '; AFFAIL(5,1)=' '; AFFAIL(6,1)=' '
AFFAIL(7,1)=' '; AFFAIL(8,1)=' '; AFFAIL(9,1)=' '
AFFAIL(1,2)='(NO'; AFFAIL(2,2)='EXAC'; AFFAIL(3,2)='T FI'
AFFAIL(4,2)='T TO'; AFFAIL(5,2)=' TRA'; AFFAIL(6,2)='NSFO'
AFFAIL(7,2)='RMS'; AFFAIL(8,2)='FOUN' ; AFFAIL(9,2)='D'
C++++++++
WTED=IWT .NE. 1
JLAMBS=1
VARBES=1.E+20
VARMIN=1.E+20
ABSYMX=0.
DO 180 J=1,N
180 ABSYMX=AMAX1(ABSYMX,ABS(Y(J)))
C-----------------------------------------------------------------------
C START OF MAIN LOOP FOR STEPWISE ANALYSIS
C-----------------------------------------------------------------------
DO 200 JLAM=1,NLAMMX
JLAMP1=JLAM+1
NF=2*JLAM+IBASE
IF (.NOT.PRFINL) GO TO 210
5200 FORMAT (1H1,5X,12A6/)
WRITE (IWRITE,5200) LABEL
5240 FORMAT (5X,23HFINAL ANALYSIS ASSUMING,I2,11H COMPONENTS)
WRITE (IWRITE,5240) JLAM
210 DO 220 J=1,JLAM
220 LAMST(J)=LAMF(J,JLAM,2)
CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,
1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,
2 4,FFAIL(JLAM,2),NDIME,.FALSE.,PRFINL)
YFAIL(JLAM)=FAILED
IF (.NOT.FAILED) GO TO 225
IF (FFAIL(JLAM,1)) GO TO 250
5222 FORMAT (////43H SECOND ATTEMPT USING STARTING LAMBDAS FROM,
1 29H INITIAL ANALYSIS OF RAW DATA)
IF (PRFINL) WRITE (IWRITE,5222)
DO 222 J=1,JLAM
222 LAMST(J)=LAMF(J,JLAM,1)
CALL LSTSQR (LAMST,JLAM,.TRUE.,N,Y,T,JLAM+IBASE,EVAR,
1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,
2 4,.FALSE.,NDIME,.FALSE.,PRFINL)
YFAIL(JLAM)=FAILED
IF (FAILED) GO TO 250
C-----------------------------------------------------------------------
C SAVE VALUES FOR LATER USE AND FOR FINAL OUTPUT
C-----------------------------------------------------------------------
225 SIGYSV(JLAM)=SIGYY
VARSV(JLAM)=VAR
ITERSV(JLAM)=ITER
ENPHSV(JLAM)=ENPHI
I=JLAM
DO 230 K=1,JLAM
ASAVE(K,JLAM,1)=PALPHA(K)
ASAVE(K,JLAM,2)=PLAM(K)
ASAVE(K,JLAM,3)=SIGMAP(K)
I=I+1
ASAVE(K,JLAM,4)=SIGMAP(I)
ASAVE(K,JLAM,5)=PCERR(K)
ASAVE(K,JLAM,6)=PCERR(I)
230 CONTINUE
IF (NOBASE) GO TO 240
ASAVE(JLAMP1,JLAM,2)=0.
ASAVE(JLAMP1,JLAM,3)=0.
ASAVE(JLAMP1,JLAM,5)=0.
ASAVE(JLAMP1,JLAM,1)=PALPHA(JLAMP1)
ASAVE(JLAMP1,JLAM,4)=SIGMAP(NF)
ASAVE(JLAMP1,JLAM,6)=PCERR(NF)
C-----------------------------------------------------------------------
C TEST IF THIS SOLUTION IS THE BEST SO FAR
C-----------------------------------------------------------------------
240 IF (VARBES .GE. 1.E+20) GO TO 248
IF (VAR .GE. VARMIN) GO TO 200
K=N-NF
DUM=VARBES/VAR-1.
L=2*(JLAM-JLAMBS)
IF (FISHER(FLOAT(K)*DUM/(FLOAT(L)*(1.+ENPHI*FLOAT((NF+2)*N)/
1 FLOAT(NF*K))),L,K) .LE. .5) GO TO 250
248 JLAMBS=JLAM
VARBES=VAR
250 VARMIN=AMIN1(VARMIN,VAR)
200 CONTINUE
C-----------------------------------------------------------------------
C END OF MAIN LOOP FOR STEPWISE ANALYSIS
C-----------------------------------------------------------------------
IF (VARBES .LT. 1.E+20) GO TO 300
5300 FORMAT (///50H ALL ANALYSES OF THE RAW DATA FAILED - CHECK INPUT,
1 22H DATA FOR GROSS ERRORS )
WRITE (IWRITE,5300)
RETURN
C-----------------------------------------------------------------------
C CALCULATE PROBSV AND PROBLN (PROBABILITY AND UNCORRECTED PROBABILITY
C (I.E., USING ENPHI=0.) THAT THE BEST SOLUTION IS ACTUALLY
C BETTER THAN THE OTHERS)
C-----------------------------------------------------------------------
300 PNXBES=1.E+20
JLAMNX=JLAMBS
DO 310 J=1,NLAMMX
PROBSV(J)=2.
IF (J.EQ.JLAMBS .OR. YFAIL(J)) GO TO 310
I=2*MAX0(JLAMBS,J)+IBASE
K=N-I
L=2*(JLAMBS-J)
DUM=VARSV(J)/VARBES
IF (L .LT. 0) GO TO 315
DDUM=FLOAT(K)*(DUM-1.)/FLOAT(L)
DDDUM=1.+ENPHSV(JLAMBS)*FLOAT((I+2)*N)/FLOAT(I*K)
PROBSV(J)=FISHER(DDUM/DDDUM,L,K)
PROBLN(J)=FISHER(DDUM,L,K)
GO TO 325
315 IF (DUM .LT. 1.) GO TO 320
PROBSV(J)=1.
PROBLN(J)=1.
GO TO 325
320 DDUM=FLOAT(K)*(1.-1./DUM)/FLOAT(L)
DDDUM=1.+ENPHSV(J)*FLOAT((I+2)*N)/FLOAT(I*K)
PROBLN(J)=1.-FISHER(DDUM,-L,K)
PROBSV(J)=1.-FISHER(DDUM/DDDUM,-L,K)
325 IF (PROBSV(J) .GE. PNXBES) GO TO 310
PNXBES=PROBSV(J)
PNXLIN=PROBLN(J)
JLAMNX=J
310 CONTINUE
PROBSV(JLAMBS)=0.
IF (PNXBES .LT. 1.E+20) GO TO 330
PNXBES=1.
PNXLIN=1.
C-----------------------------------------------------------------------
C FINAL SUMMARY OF RESULTS - SUMMARIZES UP TO BEST 5 SOLUTIONS.
C-----------------------------------------------------------------------
330 CALL RESIDU (JLAMBS,PLOTRS,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)
5330 FORMAT(35H1DISCRETE - VERSION 1A (MARCH 1976),/12A6/)
340 WRITE (IWRITE,5330) LABEL
DO 350 J=1,NLAMMX
350 DONE(J)=YFAIL(J)
L=17+JLAMBS+IBASE
5350 FORMAT (5X,A4,A3,13HBEST SOLUTION,I13,11H COMPONENTS)
WRITE (IWRITE,5350) (NEXBES(J,1),J=1,2),JLAMBS
IF (PNXBES .GT. .95) GO TO 360
5360 FORMAT (48H ++++++ NOTE - THE APPROX. PROBABILITY THAT THIS,
1 34H SOLUTION IS ACTUALLY BEST IS ONLY,F6.3,1H.,
2 /35HLOOK AT SECOND BEST SOLUTION ALSO. ,6(1H+))
WRITE (IWRITE,5360) PNXBES
L=L-1
360 DO 370 J=1,NLAMMX
IF (J .GT. 1) GO TO 375
DONE(JLAMBS)=.TRUE.
I=JLAMBS
GO TO 400
375 DUM=1.E+20
DO 380 K=1,NLAMMX
IF (PROBSV(K).GT.DUM .OR. DONE(K)) GO TO 380
DUM=PROBSV(K)
I=K
380 CONTINUE
IF (DUM.GE.1.E+20 .OR. J.GT.5) GO TO 500
DONE(I)=.TRUE.
CALL RESIDU (I,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)
L=L+15+I+IBASE
IF (L-2 .LE. LINEPG) GO TO 385
L=14+I+IBASE
WRITE (IWRITE,5999)
5999 FORMAT (1H1)
GO TO 390
385 WRITE (IWRITE,5001)
WRITE (IWRITE,5001)
5001 FORMAT (1X)
390 IF (DUM .GT. .95) WRITE (IWRITE,5350) (NEXBES(K,J),K=1,2),I
IF (DUM .LE. .95) WRITE (IWRITE,5385) (NEXBES(K,J),K=1,2),I
5385 FORMAT (5X,A4,A3,13HBEST SOLUTION,I13,11H COMPONENTS,3X,
1 30H- A SIGNIFICANT POSSIBILITY ,25(1H+))
400 WRITE (IWRITE,5001)
WRITE (IWRITE,5400)
5400 FORMAT (1X,54(1H*))
WRITE (IWRITE,5402) NFTRY(I)
5402 FORMAT (2H *,4X,22HALPHA+-STD ERR PERCENT,3X,
1 27HLAMBDA+-STD ERR PERCENT* ,
2 38HSTART LAMBDA (FROM FIT TO TRANSFORMS -,
3 I4,7H TRIES))
5404 FORMAT (2H *,1PE9.3,2H+-,E7.1,0PF8.3,1PE9.3,2H+-,
1 E7.1,0PF8.3,1H*,1PE10.3,1X, 9A4)
JJ=1
IF (FFAIL(I,2)) JJ=2
WRITE (IWRITE,5404) (ASAVE(K,I,1),ASAVE(K,I,4),ASAVE(K,I,6),
1 ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5),LAMF(K,I,2),
2 (AFFAIL(KK,JJ),KK=1,9),K=1,I)
IF (NOBASE) GO TO 410
K=I+1
WRITE (IWRITE,5404) ASAVE(K,I,1),ASAVE(K,I,4),ASAVE(K,I,6),
1 ASAVE(K,I,2),ASAVE(K,I,3),ASAVE(K,I,5)
410 WRITE (IWRITE,5400)
IF (J.EQ.1 .AND. PNXBES.GT..95) WRITE (IWRITE,5410) JLAMNX,
1 JLAMBS, PNXBES
5410 FORMAT (43X,20(1H*)/34H APPROXIMATE PROBABILITY THAT THIS,
1 56H SOLUTION IS REALLY BETTER THAN THE SECOND BEST SOLUTION,
2 /9H = * PNG(,I1,1H/,I1,3H) =,F6.3,2H */43X,20(1H*))
5412 FORMAT (/5H PNG(,I1,1H/,I1,1H),12X,1H=,F21.3,13X,4HNPHI,13X,
1 1H=,1PE11.2,13X,25H(UNCORRECTED PNG WOULD BE,0PF8.3,1H))
IF (J .GT. 1) WRITE (IWRITE,5412) I,JLAMBS,PROBSV(I),ENPHSV(I),
1 PROBLN(I)
5414 FORMAT (22H ITERATIONS IN FIT =,I21,3X,
1 18HSTD. DEV. OF FIT =,1PE11.4,
2 /27HSIGNAL/NOISE RATIO OF FIT =,0PF7.0)
DDUM=ABSYMX/SIGYSV(I)
IF (WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I)
IF (.NOT.WTED) WRITE (IWRITE,5414) ITERSV(I),SIGYSV(I),DDUM
5416 FORMAT (21H LAMBDA HELD BETEWEEN,1PE9.2,4H AND,E9.2,3X,
1 4HNPHI,3X,1H=,E11.2)
IF (J .EQ. 1) WRITE (IWRITE,5416) LAMNMX(1,1),LAMNMX(2,1),
1 ENPHSV(JLAMBS)
5418 FORMAT (/4H LAG,15X,5(2X,3HK =,I2),3X,16(1H*)/
1 20H PROB. RESIDU UNCORR,5F9.3,4X,
2 22HWEIGHTED AV = * PUNC =,F6.3,2H *)
WRITE (IWRITE,5418) (K,K=1,5),PUNCOR
5422 FORMAT (1X,71(1H=))
IF (L-1 .LE. LINEPG) WRITE (IWRITE,5422)
IF (L .LE. LINEPG) WRITE (IWRITE,5422)
370 CONTINUE
500 IF (.NOT.REPEAT) RETURN
REPEAT=.FALSE.
CALL RESIDU (JLAMBS,.FALSE.,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)
GO TO 340
END
C-----------------------------------------------------------------------
C SUBROUTINE LSTSQR. CONSTRAINED LEAST SQUARES FIT USING STEPWISE
C REGRESSION APPROACH.
C N = NO. OF DATA POINTS TO BE FIT, I.E., =N (RAWDAT=.TRUE.) OR
C =2*JLAM+IBASE (RAWDAT=.FALSE.).
C Y,T = RAW DATA, REGARDLESS OF RAWDAT.
C ITYPE=1,2,3,4 FOR PRELIMINARY DETERMINATION OF WEIGHTS, FIT TO RAW
C DATA WITH INTERPOLATION, FIT TO TRANSFORMS, AND FINAL FIT TO
C RAW DATA, RESPECTIVELY.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - VARDUM (A DUMMY NAME), PIVOT, ANLERR
C WHICH IN TURN CALL - PIVOT1, ETHEOR
C-----------------------------------------------------------------------
SUBROUTINE LSTSQR (LAMST,JLAM,RAWDAT,N,Y,T,NLIN,VARDUM,
1 GSE,NDIMG,SQRTW,YLYFIT,E,NMAX,TSTART,DELTAT,NT,NINT,
2 ITYPE,SPLIT,NDIME,INTER,PRLSQ)
EXTERNAL VARDUM
REAL LAMNMX, LAMST
LONG HOLRTH(2), LABEL
LOGICAL RAWDAT, NOBASE, PRLSQ, FAILED, ALLPIV(2), PIVLAM,
1 LDUM, INTER, PIVALP, LBLOKA, LBLOKB, LBLOKC, LBLOKD, SPLIT
COMMON /BLOK/ RBLOKA(2796), DELTA, VARMIN, FHATMX, RBLOKB(8),
1 CONVRG, RBLOKC(57), VAR, LAMNMX(2,2), PALPHA(6), PLAM(6),
2 PRECIS, RBLOKD(38), SIGMAP(11), PCERR(11), RBLOKE(252),
3 PLMTRY(6), DELTAP(10), DFCAPL(6,5), SFDE(5), ADUB(11,11),
4 DDSE(5,5), DTOLER(11), F(11,6), DF(11,6), IBLOKA(111), MDIMA,
5 IBLOKB(52), NF, IBLOKC, IWRITE, IBLOKD(5), ITER, JLAMP1
a iblokb 32->52
COMMON /BLOK/ IBLOKE(6), LBLOKA, FAILED, LBLOKB, NOBASE,
1 LBLOKC(18), PIVALP(11), LBLOKD
a lblokc 14->18
DIMENSION PIVLAM(5), R(6), DALPL(6,5)
DIMENSION LAMST(JLAM), Y(NMAX), T(NMAX),
1 SQRTW(NMAX), YLYFIT(NMAX),
2 E(NDIME,NLIN), TSTART(NINT), QT(3), VT(3), GSE(NDIMG,4),
3 DELTAT(NINT), NT(NINT), MXITER(4), RLAMMN(2)
C++++++++
qmin=.001
holrth(1)=' '; holrth(2)='*'
mconv=3; nabort=4
mxiter(1)=30; mxiter(2)=30; mxiter(3)=40; mxiter(4)=40
rlammn(1)=1.1; rlammn(2)=2
dflat=-1.e-2
c++++++++
ISPLIT=1
IF (SPLIT) ISPLIT=2
JLAMP1=JLAM+1
NCONV=0
NVARUP=0
NFLAT=0
ALLPIV(2)=.TRUE.
IF (NOBASE) PALPHA(JLAMP1)=0.
PLAM(JLAMP1)=0.
PLMTRY(JLAMP1)=0.
PIVALP(JLAMP1)=.TRUE.
C-----------------------------------------------------------------------
C PUT PLAM IN ASCENDING ORDER AND SEPARATE THEM IF NECESSARY
C-----------------------------------------------------------------------
DO 150 J=1,JLAM
PIVLAM(J)=.FALSE.
150 DELTAP(J)=0.
DO 160 J=1,JLAM
DUM=1.E+20
DO 165 K=1,JLAM
IF (PIVLAM(K) .OR. LAMST(K).GT.DUM) GO TO 165
L=K
DUM=LAMST(K)
165 CONTINUE
PIVLAM(L)=.TRUE.
IF (J .EQ. 1) PLAM(J)=AMAX1(DUM,LAMNMX(1,ISPLIT))
IF (J .GT. 1) PLAM(J)=AMIN1(AMAX1(DUM,RLAMMN(ISPLIT)*PLAM(J-1)),
1 LAMNMX(2,1))
160 CONTINUE
VAROLD=1.E+20
VARG=1.E+20
ITER=-1
Q=0.
IF (PRLSQ) WRITE (IWRITE,5205)
5205 FORMAT (/4H ITR,1X,8HVARIANCE,2X,6HDAMP Q,1X,8HBASELINE,
1 5(4X,5HALPHA,2X,6HLAMBDA))
C-----------------------------------------------------------------------
C START OF MAIN LOOP FOR THE LEAST SQUARES FIT
C-----------------------------------------------------------------------
210 NTRY=1
215 ITER=ITER+1
220 CALL VARDUM (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,
1 NINT,.TRUE.,LDUM,E,Y,YLYFIT,NLIN,VARG,NTRY,NDIME,INTER,
2 GSE,NDIMG,IERROR)
IF (IERROR .EQ. 5) GO TO 800
ALLPIV(1)=LDUM
IF (VAR.LE.VARG .OR. NTRY.EQ.2) GO TO 240
NTRY=2
Q=QG
GO TO 220
240 DO 242 J=1,JLAM
R(J)=PALPHA(J)
IF (.NOT.RAWDAT) R(J)=PALPHA(J)*EXP(-PLMTRY(J)*DELTA)
242 PLAM(J)=PLMTRY(J)
DUM=VAR-VAROLD
IF (.NOT.PRLSQ) GO TO 250
DDUM=HOLRTH(1)
IF (DUM .GT. 0.) DDUM=HOLRTH(2)
DO 244 J=1,JLAM
PCERR(J)=HOLRTH(1)
IF (.NOT.PIVLAM(J)) PCERR(J)=HOLRTH(2)
244 CONTINUE
5240 FORMAT (1X,I2,1PE9.4,A1,E8.2,A1,E8.2,5(1P2E8.2,A1))
WRITE (IWRITE,5240) ITER,VAR,DDUM,Q,HOLRTH(NTRY),PALPHA(JLAMP1),
1 (R(J),PLAM(J),PCERR(J),J=1,JLAM)
C-----------------------------------------------------------------------
C TESTS FOR CONVERGENCE AND DIVERGENCE
C-----------------------------------------------------------------------
250 IF ((ABS(DUM).GT.AMAX1(VARMIN,VAROLD*CONVRG) .AND. RAWDAT) .OR.
1 (VAR.GT.FLOAT(NF)*(CONVRG*FHATMX)**2 .AND. .NOT.RAWDAT))
2 GO TO 270
NCONV=NCONV+1
IF (NCONV.GE.MCONV .OR. .NOT.RAWDAT) GO TO 800
NVARUP=0
GO TO 300
270 NCONV=0
IF (DUM .LE. 0.) GO TO 275
NVARUP=NVARUP+1
IF (NVARUP.LT.NABORT .OR. ITER.LT.JLAM) GO TO 280
IERROR=1
GO TO 800
275 NVARUP=0
280 IF (ITER.LT.MXITER(ITYPE)) GO TO 290
IERROR=2
GO TO 800
290 IF (DUM/VAROLD.LT.DFLAT .OR. RAWDAT) GO TO 295
NFLAT=NFLAT+1
IF (NFLAT .LT. NABORT) GO TO 300
IERROR=3
GO TO 800
295 NFLAT=0
C-----------------------------------------------------------------------
C FINISH COMPUTATION OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA)
C-----------------------------------------------------------------------
300 DO 310 K=1,JLAM
DO 320 J=1,NLIN
DFCAPL(J,K)=DFCAPL(J,K)*PALPHA(K)
320 CONTINUE
DFCAPL(K,K)=DFCAPL(K,K)+SFDE(K)
310 CONTINUE
C-----------------------------------------------------------------------
C CALCULATE DALPL (PARTIAL OF PALPHA W.R.T. LAMBDA) FROM MATRIX
C PRODUCT ADUB*DFCAPL
C-----------------------------------------------------------------------
DO 400 J=1,NLIN
DO 400 K=1,JLAM
400 DALPL(J,K)=0.
DO 410 K=1,JLAM
IF (.NOT.PIVALP(K)) GO TO 410
DO 420 J=1,NLIN
IF (.NOT.PIVALP(J)) GO TO 420
DUB=0.
DO 430 I=1,NLIN
IF (PIVALP(I)) DUB=DUB+ADUB(J,I)*DFCAPL(I,K)
430 CONTINUE
DALPL(J,K)=DUB
420 CONTINUE
410 CONTINUE
C-----------------------------------------------------------------------
C CALCULATE ADUB (NORMAL EQUATION MATRIX FOR NONLINEAR LEAST SQUARES)
C AND INVERT IT
C-----------------------------------------------------------------------
IF (.NOT.RAWDAT) GO TO 550
DO 510 K=1,JLAM
DO 520 J=1,K
DUM=0.
DO 530 L=1,NLIN
530 DUM=DUM-DALPL(L,J)*DFCAPL(L,K)
IF (ABS(PALPHA(K)) .GT. ABS(PALPHA(J))) GO TO 540
DUM=DUM+DALPL(J,K)*SFDE(J)
DUM=DUM+DALPL(K,J)*SFDE(K)
GO TO 545
540 DUM=DUM+DALPL(K,J)*SFDE(K)
DUM=DUM+DALPL(J,K)*SFDE(J)
545 DDUM=PALPHA(K)*PALPHA(J)*DDSE(J,K)
ADUB(J,K)=DUM+DDUM
520 CONTINUE
DTOLER(K)=DDUM
ADUB(K,JLAMP1)=PALPHA(K)*SFDE(K)
510 CONTINUE
GO TO 590
550 DO 560 J=1,JLAM
DTOLER(J)=0.
R(J)=0.
DO 565 K=J,JLAMP1
565 ADUB(J,K)=0.
560 CONTINUE
DO 570 NU=1,N
DO 575 K=1,JLAM
IF (.NOT.PIVALP(K)) GO TO 575
DUM=0.
DO 578 J=1,NLIN
578 DUM=DUM+F(NU,J)*DALPL(J,K)
DDUM=DF(NU,K)*PALPHA(K)
DTOLER(K)=DTOLER(K)+DDUM*DDUM
R(K)=DUM+DDUM
575 CONTINUE
DDUM=YLYFIT(NU)
DO 580 K=1,JLAM
IF (.NOT.PIVALP(K)) GO TO 580
DUM=R(K)
DO 585 J=1,K
585 ADUB(J,K)=ADUB(J,K)+R(J)*DUM
ADUB(K,JLAMP1)=ADUB(K,JLAMP1)+R(K)*DDUM
580 CONTINUE
570 CONTINUE
590 DO 595 J=1,JLAM
DTOLER(J)=PRECIS*AMAX1(DTOLER(J),ADUB(J,J))
595 R(J)=ADUB(J,JLAMP1)
ADUB(JLAMP1,JLAMP1)=0.
CALL PIVOT (ADUB,MDIMA,JLAM,.FALSE.,.FALSE.,.TRUE.,LAMNMX(1,1),
1 LAMNMX(2,1),PLAM,DTOLER,DELTAP,LDUM,QG,PIVLAM,JLAM,QMAX,IERROR)
IF (IERROR .EQ. 5) GO TO 800
ALLPIV(2)=LDUM
C-----------------------------------------------------------------------
C USE MODIFICATION-A OF G. E. P. BOX TO ESTIMATE BEST Q (FRACTIONAL
C STEP SIZE)
C-----------------------------------------------------------------------
VAROLD=VAR
DDUM=0.
DO 600 J=1,JLAM
600 DDUM=DDUM+DELTAP(J)*R(J)
IF (DDUM .GT. 0.) GO TO 610
5610 FORMAT (42H GAUSS VECTOR POINTING IN WRONG DIRECTION.,2X,
1 14(2H**))
IF (PRLSQ) WRITE (IWRITE,5610)
NTRY=2
Q=AMIN1(QMIN,QMAX)
GO TO 215
610 CALL VARDUM (QG,JLAM,VARG,T,SQRTW,NMAX,TSTART,DELTAT,NT,
1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,
2 GSE,NDIMG,IERROR)
IF (IERROR .EQ. 5) GO TO 800
ALLPIV(1)=LDUM
DUM=VARG-VAROLD+2.*DDUM*QG
IF (DUM .LE. 0.) DUM=.25*DDUM*QG*QG
Q=AMIN1(DDUM*QG*QG/DUM,4.,QMAX)
IF (VARG.LE.VAROLD .OR. (Q.GE..125*QG .AND. RAWDAT)) GO TO 210
C-----------------------------------------------------------------------
C PUT QUADRATIC PARABOLA THRU ACTUAL VARIANCE SURFACE WHEN A
C STRONG NONLINEARITY IS INDICATED BY A SMALL Q
C-----------------------------------------------------------------------
QT(1)=AMIN1(.125*QG,AMAX1(Q,.02))
IF (Q .GT. .125*QG) QT(1)=Q
CALL VARDUM (QT(1),JLAM,DUM,T,SQRTW,NMAX,TSTART,DELTAT,NT,
1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,
2 GSE,NDIMG,IERROR)
IF (IERROR .EQ. 5) GO TO 800
VT(1)=DUM
IF (VT(1) .GE. VAROLD) GO TO 615
QT(2)=4.*QT(1)
IF (QT(1) .LE. .125*QG) GO TO 617
QT(2)=0.
VT(2)=VAROLD
QT(3)=QG
VT(3)=VARG
GO TO 630
615 QT(2)=.5*QT(1)
QT(2)=AMIN1(QT(2),.05)
617 CALL VARDUM (QT(2),JLAM,DUM,T,SQRTW,NMAX,TSTART,DELTAT,NT,
1 NINT,.FALSE.,LDUM,E,Y,YLYFIT,NLIN,VAR,NTRY,NDIME,INTER,
2 GSE,NDIMG,IERROR)
IF (IERROR .EQ. 5) GO TO 800
VT(2)=DUM
IF (VT(1).GE.VAROLD .OR. VT(2).GE.VT(1)) GO TO 620
QT(3)=QG
VT(3)=VARG
GO TO 630
620 QT(3)=0.
VT(3)=VAROLD
630 VARG=AMIN1(VT(1),VT(2),VARG)
IF (VT(1) .LE. VARG) QG=QT(1)
IF (VT(2) .LE. VARG) QG=QT(2)
DUM=(VT(1)-VT(3))/(QT(1)-QT(3))
DDUM=(DUM-(VT(2)-VT(3))/(QT(2)-QT(3)))/(QT(1)-QT(2))
IF (DDUM .GT. 0.) GO TO 650
640 Q=AMIN1(Q,.125*QT(2),.01)
GO TO 210
650 DUM=.5*((QT(1)+QT(3))-DUM/DDUM)
IF (DUM .LE. 0.) GO TO 640
Q=AMIN1(DUM,8.,QMAX)
GO TO 210
C-----------------------------------------------------------------------
C CHECK FOR ABNORMAL EXIT AND CALL ANLERR
C-----------------------------------------------------------------------
800 FAILED=IERROR.NE.0 .OR. .NOT.ALLPIV(1)
IF (.NOT.RAWDAT) GO TO 890
FAILED=FAILED .OR. .NOT.ALLPIV(2)
IF (.NOT.FAILED .OR. ITYPE.EQ.4) CALL ANLERR (JLAM,T,SQRTW,NMAX,E,
1 NLIN,NDIME,TSTART,DELTAT,NT,NINT,N,ITYPE,INTER,PRLSQ)
890 IF (.NOT.(PRLSQ .AND. FAILED)) GO TO 950
5889 FORMAT (/47H VARIANCE DID NOT SIGNIFICANTLY DECREASE IN THE,
1 5H LAST,I2,7H TIMES.,5X,6(11H***TOO SLOW))
IF (IERROR .EQ. 3) WRITE (IWRITE,5889) NABORT
5890 FORMAT (/41H MAXIMUM ITERATIONS WITHOUT CONVERGENCE. ,
1 7(13H***ITERATIONS))
IF (IERROR .EQ. 2) WRITE (IWRITE,5890)
5892 FORMAT (/28H VARIANCE INCREASED THE LAST,I2,8H TIMES. ,
1 3X,7(13H***DIVERGENCE))
IF (IERROR .EQ. 1) WRITE (IWRITE,5892) NVARUP
5894 FORMAT (/42H NEARLY SINGULAR NORMAL EQUATIONS MATRIX. ,
1 6(15H****SINGULARITY))
IF (IERROR .EQ. 5) WRITE (IWRITE,5894)
C-----------------------------------------------------------------------
C PUT PLAM IN ASCENDING ORDER
C-----------------------------------------------------------------------
950 DO 960 J=1,JLAM
DUM=1.E+20
DO 985 K=J,JLAM
IF (PLAM(K) .GT. DUM) GO TO 985
L=K
DUM=PLAM(K)
985 CONTINUE
PLAM(L)=PLAM(J)
PLAM(J)=DUM
DUM=PALPHA(L)
PALPHA(L)=PALPHA(J)
PALPHA(J)=DUM
IF (.NOT.RAWDAT) GO TO 960
DUM=SIGMAP(L)
SIGMAP(L)=SIGMAP(J)
SIGMAP(J)=DUM
K=L+JLAM
JK=J+JLAM
DUM=SIGMAP(K)
SIGMAP(K)=SIGMAP(JK)
SIGMAP(JK)=DUM
PCERR(J)=100.*SIGMAP(J)/PLAM(J)
IF (ABS(PALPHA(J)) .GT. 0.) PCERR(JK)=100.*SIGMAP(JK)/ABS(PALPHA
1 (J))
960 CONTINUE
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE EVAR. FOR FRACTIONAL STEP SIZE Q, EVALUATES PALPHA
C (LINEAR LEAST SQUARES PARAMETERS), VAR (WEIGHTED VARIANCE) AND
C ACCUMULATES IT IN VGRID2 (COARSE GRID SUM), AND, IF
C INVERT=.TRUE., EVALUATES ARRAYS SFDE, DFCAPL, AND FULL UPPER
C TRIANGLES OF DDSE AND SE FOR LATER USE IN NONLINEAR LEAST
C SQUARES.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - ETHEOR, PIVOT
C WHICH IN TURN CALL - PIVOT1
C-----------------------------------------------------------------------
SUBROUTINE EVAR (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,NINT,
1 INVERT,ALLPIV,E,Y,YLYDUM,NLIN,VARG,NTRY,NDIME,INTER,GSE,NDIMG,
2 IERROR)
LOGICAL INVERT, ALLPIV, LDDUM, INTER, LBLOKA, LBLOKB, LBLOKC,
1 LBLOKD, NOTINV, NOBASE, WTED, REGINT, NONNEG, PIVALP
REAL LAMNMX
LONG LABEL
COMMON /BLOK/ RBLOKA(2799), SUMWTY, SUMWT, ZERO(6),
1 RBLOKB(56), DGRIDE, GESTL, RBLOKC, LAMNMX(2,2), PALPHA(6),
2 PLAM(6), PRECIS, RBLOKD(2), SE(6,6), RBLOKE(22), VGRID2(252),
3 PLMTRY(6), DELTAP(10), DFCAPL(6,5), SFDE(5), ADUB(11,11),
4 DDSE(5,5), DTOLER(11), RBLOKF(132), JGEADD, NBINOM(20,5),
5 IBLOKA(5), NPERG2(5), MDIMA, IBLOKB(4), NGRID2(5),IBLOKC(44), N
a iblokc 24->44
COMMON /BLOK/ IBLOKD(14), REGINT, LBLOKA(4), NOBASE, WTED,
1 LBLOKB(4), NONNEG, LBLOKC(10), PIVALP(11), LBLOKD
a lbloka 2->4, lblokb 2->4
DIMENSION ADUM(6), FACTOR(5), XEXP(5), JG2(5)
DIMENSION T(NMAX), SQRTW(NMAX), TSTART(NINT), DELTAT(NINT),
1 NT(NINT), E(NDIME,NLIN), Y(NMAX), YLYDUM(1), GSE(NDIMG,4)
DO 110 J=1,JLAM
110 PLMTRY(J)=PLAM(J)+Q*DELTAP(J)
NLINP=NLIN+1
DO 160 K=1,NLINP
160 ADUB(K,NLINP)=0.
IF (NOBASE) GO TO 170
PLMTRY(NLIN)=0.
ADUB(NLIN,NLIN)=SUMWT
ADUB(NLIN,NLINP)=SUMWTY
170 NOTINV=.NOT.INVERT
IF (REGINT) GO TO 300
C-----------------------------------------------------------------------
C ORDINARY COMPUTATION OF E (EXPONENTIALS), ADUB (NORMAL EQUATIONS
C MATRIX FOR LINEAR LEAST SQUARES) AND, IF INVERT=.TRUE., PART
C OF DFCAPL (PARTIAL OF FCAP W.R.T. LAMBDA) WHEN REGINT=.FALSE.
C-----------------------------------------------------------------------
DO 190 K=1,NLIN
DO 195 J=1,K
IF (J .GT. JLAM) GO TO 190
ADUB(J,K)=0.
DFCAPL(K,J)=0.
195 IF (K .LE. JLAM) DDSE(J,K)=0.
190 CONTINUE
DO 200 J=1,N
DUM=-T(J)
DO 210 K=1,JLAM
ADUM(K)=SQRTW(J)*EXP(PLMTRY(K)*DUM)
210 E(J,K)=ADUM(K)
IF (NOBASE) GO TO 215
ADUM(NLIN)=SQRTW(J)
E(J,NLIN)=ADUM(NLIN)
215 DO 220 K=1,NLIN
DDUM=ADUM(K)
LDDUM=K .LE. JLAM
DO 230 L=1,K
IF (L .GT. JLAM) GO TO 220
DDDUM=ADUM(L)*DDUM
ADUB(L,K)=ADUB(L,K)+DDDUM
IF (NOTINV) GO TO 230
DDDUM=DDDUM*DUM
DFCAPL(K,L)=DFCAPL(K,L)-DDDUM
IF (LDDUM) DDSE(L,K)=DDSE(L,K)+DDDUM*DUM
230 CONTINUE
ADUB(K,NLINP)=ADUB(K,NLINP)+DDUM*Y(J)
220 CONTINUE
200 CONTINUE
GO TO 500
C-----------------------------------------------------------------------
C EVALUATE ADUB AND, IF INVERT=.TRUE., PART OF DFCAPL AND DDSE
C (SECOND DERIVATIVES OF EXPONENTIAL SUMS W.R.T. LAMBDA)
C BY HERMITE AND LAGRANGE INTERPOLATION, IF REGINT=INTER=.TRUE.
C-----------------------------------------------------------------------
300 IF (.NOT.INTER) GO TO 400
DO 310 K=1,NLIN
LDDUM=K .LE. JLAM
DO 320 J=1,K
IF (J .GT. JLAM) GO TO 310
XLAM=PLMTRY(K)+PLMTRY(J)
DUM=(ALOG(XLAM)-GESTL)/DGRIDE
JZ=(DUM+.5)
HZ=DUM-FLOAT(JZ)
IF (ABS(HZ) .GT. 1.E-5) GO TO 330
ADUB(J,K)=GSE(JZ,1)
IF (NOTINV) GO TO 320
DFCAPL(K,J)=-GSE(JZ,2)/XLAM
IF (LDDUM) DDSE(J,K)=(GSE(JZ,3)-GSE(JZ,2))/(XLAM*XLAM)
GO TO 320
330 JM=(DUM)
JP=JM+1
H=DUM-FLOAT(JM)
HH=1.-H
DDUM=1.+3.*H
DDDUM=1.+3.*HH
DX=H*DGRIDE
DDX=-HH*DGRIDE
ADUB(J,K)=HH*HH*HH*((DDUM+6.*H*H)*GSE(JM,1)+DX*(DDUM*GSE(JM,2)+
1 .5*DX*GSE(JM,3)))+H*H*H*((DDDUM+6.*HH*HH)*GSE(JP,1)+
2 DDX*(DDDUM*GSE(JP,2)+.5*DDX*GSE(JP,3)))
IF (NOTINV) GO TO 320
HH=HZ*HZ
DX=(HH-1.)*HZ*.5
DUM=1./(HZ+1.)
DDUM=1./HZ
DDDUM=1./(HZ-1.)
H=DX*DX*(((3.*HZ+4.)*GSE(JZ-1,2)*DUM+DGRIDE*GSE(JZ-1,3))*DUM+
1 4.*DDUM*(DDUM*GSE(JZ,2)+DGRIDE*GSE(JZ,3))+DDDUM*((4.-3.*HZ)*
2 GSE(JZ+1,2)*DDDUM+DGRIDE*GSE(JZ+1,3)))
DFCAPL(K,J)=-H/XLAM
IF (.NOT.LDDUM) GO TO 320
DUM=DX*(HH-4.)*(HH-9.)*(GSE(JZ+3,3)/(HZ-3.)+GSE(JZ-3,3)/(HZ+3.)-
1 6.*(GSE(JZ+2,3)/(HZ-2.)+GSE(JZ-2,3)/(HZ+2.))+15.*
2 (GSE(JZ+1,3)*DDDUM+GSE(JZ-1,3)*DUM)-20.*GSE(JZ,3)*DDUM)/360.
DDSE(J,K)=(DUM-H)/(XLAM*XLAM)
320 CONTINUE
DUM=(ALOG(PLMTRY(K))-GESTL)/DGRIDE
JZ=(DUM+.5)
HZ=DUM-FLOAT(JZ)
IF (ABS(HZ) .GT. 1.E-5) GO TO 340
ADUB(K,NLINP)=GSE(JZ,4)
GO TO 310
340 HH=HZ*HZ
ADUB(K,NLINP)=HZ*(HH-1.)*(HH-4.)*(HH-9.)*(GSE(JZ+3,4)/(HZ-3.)+
1 GSE(JZ-3,4)/(HZ+3.)-6.*(GSE(JZ+2,4)/(HZ-2.)+GSE(JZ-2,4)/
2 (HZ+2.))+15.*(GSE(JZ+1,4)/(HZ-1.)+GSE(JZ-1,4)/(HZ+1.))
3 -20.*GSE(JZ,4)/HZ)/720.
310 CONTINUE
GO TO 500
C-----------------------------------------------------------------------
C EVALUATE ADUB, PART OF DFCAPL, AND DDSE BY CALLING ETHEOR
C IF REGINT=.TRUE. AND INTER=.FALSE.
C-----------------------------------------------------------------------
400 DO 410 K=1,NLIN
LDDUM=K .LE. JLAM
DO 420 J=1,K
IF (J .GT. JLAM) GO TO 410
XLAM=PLMTRY(K)+PLMTRY(J)
CALL ETHEOR (XLAM,INVERT,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,
1 NT,NINT,DUM,DDUM,DDDUM,S,.TRUE.,.FALSE.)
ADUB(J,K)=DUM
IF (NOTINV) GO TO 420
DFCAPL(K,J)=-DDUM/XLAM
IF (LDDUM) DDSE(J,K)=(DDDUM-DDUM)/(XLAM*XLAM)
420 CONTINUE
CALL ETHEOR (PLMTRY(K),.FALSE.,WTED,T,SQRTW,Y,NMAX,TSTART,DELTAT,
1 NT,NINT,DUM,DDUM,DDDUM,S,.FALSE.,.TRUE.)
ADUB(K,NLINP)=S
410 CONTINUE
500 DO 510 K=1,NLIN
DTOLER(K)=10.*PRECIS*ADUB(K,K)
DO 515 J=1,K
515 SE(J,K)=ADUB(J,K)
510 CONTINUE
C-----------------------------------------------------------------------
C SOLVE NORMAL EQUATIONS FOR PALPHA USING STEPWISE REGRESSION
C-----------------------------------------------------------------------
CALL PIVOT (ADUB,MDIMA,NLIN,.FALSE.,INVERT,NONNEG,0.,
1 1.E+20,ZERO,DTOLER,PALPHA,ALLPIV,DUM,PIVALP,JLAM,S,IERROR)
IF (IERROR .NE. 5) GO TO 518
VAR=1.E+20
RETURN
C-----------------------------------------------------------------------
C CALCULATE VAR AND, IF INVERT=.TRUE., SFDE (SUMS OF RESIDUALS
C TIMES DERIVATIVES OF NONLINEAR TERMS)
C-----------------------------------------------------------------------
518 VAR=0.
IF (NOTINV) GO TO 528
DO 520 K=1,JLAM
SFDE(K)=0.
DO 525 J=1,K
525 DFCAPL(J,K)=DFCAPL(K,J)
520 CONTINUE
528 IF (.NOT.REGINT) GO TO 600
K=0
DO 530 J=1,NINT
DO 540 L=1,JLAM
FACTOR(L)=EXP(-PLMTRY(L)*DELTAT(J))
540 XEXP(L)=EXP(-PLMTRY(L)*TSTART(J))
KK=NT(J)
DO 550 LL=1,KK
K=K+1
DUM=SQRTW(K)
YLYFIT=Y(K)
IF (.NOT.NOBASE) YLYFIT=YLYFIT-PALPHA(NLIN)*DUM
DO 560 L=1,JLAM
ADUM(L)=XEXP(L)*DUM
YLYFIT=YLYFIT-ADUM(L)*PALPHA(L)
XEXP(L)=XEXP(L)*FACTOR(L)
560 CONTINUE
VAR=VAR+YLYFIT*YLYFIT
IF (NOTINV) GO TO 550
DDUM=T(K)
DO 570 L=1,JLAM
570 SFDE(L)=SFDE(L)-DDUM*YLYFIT*ADUM(L)
550 CONTINUE
530 CONTINUE
GO TO 700
600 DO 610 J=1,N
YLYFIT=Y(J)
DO 620 K=1,NLIN
620 YLYFIT=YLYFIT-PALPHA(K)*E(J,K)
VAR=VAR+YLYFIT*YLYFIT
IF (NOTINV) GO TO 610
DUM=T(J)
DO 630 K=1,JLAM
630 SFDE(K)=SFDE(K)-DUM*YLYFIT*E(J,K)
610 CONTINUE
700 IF (JLAM.EQ.1 .OR. .NOT.INTER) GO TO 800
C-----------------------------------------------------------------------
C ADD VAR INTO VGRID2 (COARSE GRID SUM)
C-----------------------------------------------------------------------
J2MAX=0
DO 710 L=1,JLAM
IF (PLMTRY(L).LT.LAMNMX(1,2) .OR. PLMTRY(L).GT.LAMNMX(2,2))
1 GO TO 800
JG2(L)=(((ALOG(PLMTRY(L))-GESTL)/DGRIDE+.5)-JGEADD)/
1 NPERG2(JLAM)+1
IF (JG2(L) .LT. J2MAX) GO TO 710
J2MAX=JG2(L)
LMAX=L
710 CONTINUE
JG2(LMAX)=0
LOC=1
DO 720 J=1,JLAM
J2NEXT=0
DO 730 L=1,JLAM
IF (JG2(L) .LT. J2NEXT) GO TO 730
J2NEXT=JG2(L)
LNEXT=L
730 CONTINUE
L=J2NEXT+1
IF (J2MAX-L) 800,750,740
740 LL=J2MAX-1
DO 745 K=L,LL
KK=NGRID2(JLAM)-K+1
745 LOC=LOC+NBINOM(KK,J)
750 J2MAX=J2NEXT
JG2(LNEXT)=0
720 CONTINUE
DUM=VAR
LOC=MIN0(LOC,252)
IF (.NOT.INVERT .OR. (VAR.GT.VARG .AND. NTRY.EQ.1)) DUM=DUM*1.E-6
VGRID2(LOC)=VGRID2(LOC)+DUM
800 RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE VARF. SAME AS EVAR, EXCEPT FOR TRANSFORMS INSTEAD OF
C EXPONENTIALS AND DOES NO ACCUMULATION IN A COARSE GRID.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - PIVOT
C WHICH IN TURN CALL - PIVOT1
C-----------------------------------------------------------------------
SUBROUTINE VARF (Q,JLAM,VAR,T,SQRTW,NMAX,TSTART,DELTAT,NT,NINT,
1 INVERT,ALLPIV,E,Y,YLYFIT,NLIN,VARG,NTRY,NDIME,INTDUM,GSE,NDIMG,
2 IERROR)
LOGICAL INVERT, ALLPIV, INTDUM, LBLOKA, LBLOKB, LBLOKC,
1 LBLOKD, NOBASE, NONNEG, PIVALP, NOTINV
LONG LABEL
COMMON /BLOK/ TRBASE(11), GF(11,126), DGFL(11,126), DGRIDF,
1 GFSTL, FHAT(11), RBLOKA(5), ZERO(6), RBLOKB(63),
2 PALPHA(6), PLAM(6), PRECIS, RBLOKC(312), PLMTRY(6),
3 DELTAP(10), DFCAPL(6,5), RBLOKD(5), ADUB(11,11), RBLOKE(25),
4 DTOLER(11), F(11,6), DF(11,6), IBLOKA(111), MDIMA, IBLOKB(52),
5 NF, IBLOKC(15), LBLOKA(6), NOBASE, LBLOKB(6), NONNEG, LBLOKC(10)
a iblokb 32->52
COMMON /BLOK/ PIVALP(11), LBLOKD
a lbloka 3->6, lblokb 3->6
DIMENSION T(1), SQRTW(1), TSTART(1), DELTAT(1), NT(1), E(1,1),
1 Y(1), YLYFIT(NMAX), ADUM(9), GSE(1,1)
NLINP=NLIN+1
NOTINV=.NOT.INVERT
IF (NOBASE) GO TO 190
PLMTRY(NLIN)=0.
DO 130 J=1,NF
130 F(J,NLIN)=TRBASE(J)
C-----------------------------------------------------------------------
C COMPUTE F (NONLINEAR TERMS) AND, IF INVERT=.TRUE., DF (THEIR
C DERIVATIVES W.R.T. LAMBDA) BY HERMITE AND LAGRANGE
C INTERPOLATION.
C-----------------------------------------------------------------------
190 DO 200 J=1,JLAM
PLMTRY(J)=PLAM(J)+Q*DELTAP(J)
RLAM=1./PLMTRY(J)
DUM=(ALOG(PLMTRY(J))-GFSTL)/DGRIDF
JZ=(DUM+.5)
HZ=DUM-FLOAT(JZ)
IF (ABS(HZ) .GT. 1.E-5) GO TO 220
DO 210 K=1,NF
F(K,J)=GF(K,JZ)
210 IF (INVERT) DF(K,J)=DGFL(K,JZ)*RLAM
GO TO 200
220 ADUM(7)=1./HZ
ADUM(1)=4.*ADUM(7)*ADUM(7)
ADUM(2)=4.*DGRIDF*ADUM(7)
ADUM(8)=1./(HZ+1.)
ADUM(3)=(3.*HZ+4.)*ADUM(8)*ADUM(8)
ADUM(4)=DGRIDF*ADUM(8)
ADUM(9)=1./(HZ-1.)
ADUM(5)=(4.-3.*HZ)*ADUM(9)*ADUM(9)
ADUM(6)=DGRIDF*ADUM(9)
HH=HZ*HZ
DDUM=(HH-1.)*HZ*.5
DDDUM=DDUM*DDUM
DO 230 K=1,NF
230 F(K,J)=DDDUM*(ADUM(1)*GF(K,JZ)+ADUM(2)*DGFL(K,JZ)+ADUM(3)*
1 GF(K,JZ-1)+ADUM(4)*DGFL(K,JZ-1)+ADUM(5)*GF(K,JZ+1)+ADUM(6)*
2 DGFL(K,JZ+1))
IF (NOTINV) GO TO 200
ADUM(1)=-20.*ADUM(7)
ADUM(2)=15.*ADUM(8)
ADUM(3)=15.*ADUM(9)
ADUM(4)=-6./(HZ+2.)
ADUM(5)=-6./(HZ-2.)
ADUM(6)=1./(HZ+3.)
ADUM(7)=1./(HZ-3.)
DDDUM=DDUM*(HH-4.)*(HH-9.)*RLAM/360.
DO 240 K=1,NF
240 DF(K,J)=DDDUM*(ADUM(1)*DGFL(K,JZ)+ADUM(2)*DGFL(K,JZ-1)+ADUM(3)*
1 DGFL(K,JZ+1)+ADUM(4)*DGFL(K,JZ-2)+ADUM(5)*DGFL(K,JZ+2)+
2 ADUM(6)*DGFL(K,JZ-3)+ADUM(7)*DGFL(K,JZ+3))
200 CONTINUE
C-----------------------------------------------------------------------
C COMPUTE AND INVERT ADUB (NORMAL EQUATIONS MATRIX FOR LINEAR
C LEAST SQUARES
C-----------------------------------------------------------------------
ADUB(NLINP,NLINP)=0.
DO 300 K=1,NLIN
DO 310 J=1,K
DUM=0.
DO 320 L=1,NF
320 DUM=DUM+F(L,J)*F(L,K)
ADUB(J,K)=DUM
310 CONTINUE
DTOLER(K)=10.*PRECIS*ADUB(K,K)
DUM=0.
DO 330 L=1,NF
330 DUM=DUM+F(L,K)*FHAT(L)
ADUB(K,NLINP)=DUM
300 CONTINUE
CALL PIVOT (ADUB,MDIMA,NLIN,.FALSE.,INVERT,NONNEG,0.,
1 1.E+20,ZERO,DTOLER,PALPHA,ALLPIV,DUM,PIVALP,JLAM,DDUM,IERROR)
IF (IERROR .EQ. 5) RETURN
C-----------------------------------------------------------------------
C CALCULATE VAR (WEIGHTED VARIANCE OF FIT) AND YLYFIT (WEIGHTED
C RESIDUALS), AND,IF INVERT=.TRUE., PART OF DFCAPL FOR
C LATER USE IN NONLINEAR LEAST SQUARES
C-----------------------------------------------------------------------
VAR=0.
DO 410 L=1,NF
DUM=FHAT(L)
DO 420 J=1,NLIN
420 DUM=DUM-PALPHA(J)*F(L,J)
YLYFIT(L)=DUM
VAR=VAR+DUM*DUM
410 CONTINUE
IF (NOTINV .OR. (VAR.GT.VARG .AND. NTRY.EQ.1)) RETURN
DO 450 K=1,JLAM
DO 470 J=1,NLIN
DUM=0.
DO 480 L=1,NF
480 DUM=DUM-F(L,J)*DF(L,K)
DFCAPL(J,K)=DUM
470 CONTINUE
450 CONTINUE
RETURN
END
C-----------------------------------------------------------------------
C++++++++++++
FUNCTION EXPSMA(S)
EXPSMA=S*(720-S*(360-S*(120-S*(30-S*(6-S)))))/720
RETURN
END
C+++++++++++
C SUBROUTINE ETHEOR. IF DOE=.TRUE., PUTS EXPONENTIAL SUM AND, IF
C INVERT=.TRUE., ITS FIRST AND SECOND DERIVATIVES W.R.T.
C ALOG(LAMBDA) IN SE, DSE, DDSE. IF DOYE=.TRUE., PUTS SUM OF Y
C TIMES EXPONENTIALS IN SYE. DONE IN SINGLE PRECISION.
C CAN ONLY BE USED WHEN REGINT=.TRUE..
C IF WTED=.TRUE., USES RECURSION RELATIONS, AND IF WTED=.FALSE. USES
C GEOMETRIC SUM FORMULAS FOR SE, DSE, AND DDSE.
C NO COMMON BLOCKS NEEDED.
C-----------------------------------------------------------------------
C CALLS NO SUBPROGRAMS.
C-----------------------------------------------------------------------
SUBROUTINE ETHEOR (PLAMBD,INVERT,WTED,T,SQRTW,Y,NMAX,TSTART,
1 DELTAT,NT,NINT,SE,DSE,DDSE,SYE,DOE,DOYE)
REAL AN, B, C, DDDUM, DDSEI, DDUM,
1 DEDNL, DL, DNL, DSEI, DUM, EDL, ELTS, ENTS, RDEDL, SEI, S,
2 SS, SYEI, DDT, DTS, DPLAM
LOGICAL INVERT, NOTINV, WTED, DOE, DOYE, LDDUM
DIMENSION SQRTW(NMAX), TSTART(NINT), DELTAT(NINT), NT(NINT),
1 T(NMAX), Y(NMAX)
A EXPSMA(S)=S*(720.D0-S*(360.D0-S*(120.D0-S*(30.D0-S*(6.D0-S)))))/
A 1 720.D0
DPLAM=(PLAMBD)
IF (DOE) SEI=0.D0
IF (INVERT) DSEI=0.D0
IF (INVERT) DDSEI=0.D0
NOTINV=.NOT.INVERT
IF (.NOT.(DOYE .OR. WTED)) GO TO 200
IF (DOYE) SYEI=0.D0
LDDUM=.NOT.(DOE .AND. WTED)
J=0
DO 110 L=1,NINT
DDUM=EXP(-DPLAM*(DELTAT(L)))
DDDUM=EXP(-DPLAM*(TSTART(L)))
KK=NT(L)
DO 120 K=1,KK
J=J+1
DUM=DDDUM*SQRTW(J)
DDDUM=DDDUM*DDUM
IF (DOYE) SYEI=SYEI+DUM*Y(J)
IF (LDDUM) GO TO 120
DUM=DUM*SQRTW(J)
SEI=SEI+DUM
IF (NOTINV) GO TO 120
DUM=DUM*T(J)
DSEI=DSEI-DUM
DDSEI=DDSEI+DUM*T(J)
120 CONTINUE
110 CONTINUE
IF (DOYE) SYE=(SYEI)
IF (.NOT.DOE) RETURN
IF (WTED) GO TO 300
C-----------------------------------------------------------------------
C USE GEOMETRIC SUM FORMULAS WHEN WTED=.FALSE.
C-----------------------------------------------------------------------
200 DO 210 L=1,NINT
DDT=(DELTAT(L))
DTS=(TSTART(L))
DL=DPLAM*DDT
AN=(FLOAT(NT(L)))
DNL=AN*DL
IF (ABS(DL) .LE. 1.D-3) GO TO 214
RDEDL=1.D0/(1.D0-EXP(-DL))
DEDNL=1.D0-EXP(-DNL)
GO TO 218
214 RDEDL=1.D0/EXPSMA(DL)
IF (ABS(DNL) .LE. 1.D-3) DEDNL=EXPSMA(DNL)
IF (ABS(DNL) .GT. 1.D-3) DEDNL=1.D0-EXP(-DNL)
218 ELTS=EXP(-DPLAM*DTS)
S=ELTS*DEDNL*RDEDL
SEI=SEI+S
IF (NOTINV) GO TO 210
B=DDT*RDEDL
EDL=EXP(-DL)
ENTS=ELTS*EXP(-DNL)
C=AN*ENTS-S*EDL
SS=B*C-DTS*S
DSEI=DSEI+SS
DDSEI=DDSEI-DTS*SS-B*(EDL*(B*C+SS-S*DDT)+AN*(AN*DDT+DTS)*ENTS)
210 CONTINUE
300 SE=(SEI)
IF (NOTINV) RETURN
DUM=DPLAM*DSEI
DSE=(DUM)
DDSE=(DUM+DDSEI*DPLAM*DPLAM)
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE PIVOT. DOES GAUSS-JORDAN PIVOTS ON ALL NP DIAGONAL
C ELEMENTS OF AUGMENTED MATRIX ADUB EXCEPT THOSE THAT HAVE A
C TOLERANCE LESS THAN DTOLER OR THAT CAUSE A VIOLATION OF THE
C CONSTRAINT PMIN .LE. P .LE. PMAX.
C ONLY NEEDS FULL UPPER TRIANGLE OF MATRIX ON INPUT, AND OVERWRITES IT
C WITH INVERSE AND SOLUTION.
C NEEDS NO COMMON BLOCKS.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - PIVOT1
C-----------------------------------------------------------------------
SUBROUTINE PIVOT (ADUB,MDIMA,NP,ONLYIN,INVERT,CONSTR,PMIN,
1 PMAX,P,DTOLER,DELTAP,ALLPIV,QG,PIV,JLAM,QMAX,IERROR)
LOGICAL ONLYIN, ALLPIV, CONSTR, INVERT, PIV(NP)
DIMENSION DELTAP(NP), P(NP), ADUB(MDIMA,MDIMA), DTOLER(NP)
IERROR=0
QG=1.
QMAX=999.
NPP=NP+1
IF (ONLYIN) NPP=NP
DO 110 J=1,NP
110 PIV(J)=.FALSE.
IF (NPP .EQ. 1) GO TO 190
DO 120 J=2,NPP
L=J-1
DO 120 K=1,L
ADUB(J,K)=ADUB(K,J)
120 CONTINUE
C-----------------------------------------------------------------------
C MAIN LOOP FOR COMPLETE PIVOTING
C-----------------------------------------------------------------------
190 DO 200 NPIV=1,NP
C-----------------------------------------------------------------------
C FIND NEXT PIVOT ELEMENT
C-----------------------------------------------------------------------
DUM=0.
DO 210 J=1,NP
IF (PIV(J) .OR. ADUB(J,J).LE.DTOLER(J)) GO TO 210
DDUM=ADUB(J,NPP)*ADUB(J,NPP)/ADUB(J,J)
IF (DDUM .LE. DUM) GO TO 210
DUM=DDUM
L=J
210 CONTINUE
IF (DUM .LE. 0.) GO TO 300
PIV(L)=.TRUE.
CALL PIVOT1 (ADUB,MDIMA,NPP,L)
200 CONTINUE
NPIV=NP+1
300 NPIV=NPIV-1
IF (ONLYIN .OR. .NOT.CONSTR) GO TO 400
C-----------------------------------------------------------------------
C ENFORCE CONSTRAINTS BY MAKING QG.LT.1. OR BY UNPIVOTING
C-----------------------------------------------------------------------
DO 330 J=1,JLAM
DUM=1.E+20
DO 340 K=1,JLAM
IF (.NOT.PIV(K)) GO TO 340
RDELTA=1./ADUB(K,NPP)
DDUM=AMAX1((PMIN-P(K))*RDELTA,(PMAX-P(K))*RDELTA)
IF (DDUM .GE. DUM) GO TO 340
DUM=DDUM
L=K
340 CONTINUE
IF (DUM .GT. 1.E-5) GO TO 350
NPIV=NPIV-1
PIV(L)=.FALSE.
CALL PIVOT1 (ADUB,MDIMA,NPP,L)
330 CONTINUE
C-----------------------------------------------------------------------
C IF ALL ELEMENTS ARE SOMEHOW UNPIVOTED, TAKE ERROR EXIT
C-----------------------------------------------------------------------
350 IF (NPIV .GT. 0) GO TO 390
IERROR=5
RETURN
390 QG=AMIN1(DUM,1.)
QMAX=DUM
400 ALLPIV=NPIV .EQ. NP
IF (ONLYIN) RETURN
DO 430 J=1,NP
DELTAP(J)=0.
430 IF (PIV(J)) DELTAP(J)=ADUB(J,NPP)
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE PIVOT1. DOES A GAUSS-JORDAN PIVOT ON DIAGONAL ELEMENT LPIV
C NEEDS NO COMMON BLOCKS.
C-----------------------------------------------------------------------
C CALLS NO SUBPROGRAMS
C-----------------------------------------------------------------------
SUBROUTINE PIVOT1 (ADUB,MDIMA,NPP,LPIV)
DIMENSION ADUB(MDIMA,MDIMA)
DUB=1./ADUB(LPIV,LPIV)
DO 110 J=1,NPP
IF (J .EQ. LPIV) GO TO 110
DDUB=ADUB(J,LPIV)*DUB
DO 120 K=1,NPP
120 IF (K .NE. LPIV) ADUB(J,K)=ADUB(J,K)-ADUB(LPIV,K)*DDUB
110 CONTINUE
DO 130 J=1,NPP
ADUB(J,LPIV)=-ADUB(J,LPIV)*DUB
130 ADUB(LPIV,J)=ADUB(LPIV,J)*DUB
ADUB(LPIV,LPIV)=DUB
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE ANLERR. DOES FINAL ERROR ANALYSIS OF FIT AND CALCULATES
C NONLINEARITY TERM ENPHI OF BEALE FOR LATER USE.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - PIVOT
C WHICH IN TURN CALL - PIVOT1
C-----------------------------------------------------------------------
SUBROUTINE ANLERR (JLAM,T,SQRTW,NMAX,E,NLIN,NDIME,TSTART,DELTAT,
1 NT,NINT,N,ITYPE,INTER,PRLSQ)
LOGICAL LDUM, NOBASE, PIVALP, FAILED, PRLSQ, REGINT,
1 LBLOKA, LBLOKB, LBLOKC, INTER
LONG HOLRTH, ALPLAM(5,2), LABEL
COMMON /BLOK/ RBLOKA(2800), SUMWT, RBLOKB(64), VAR, RBLOKC(4),
1 PALPHA(6), PLAM(6), PRECIS, SIGYY, ENPHI, SE(6,6), SIGMAP(11),
2 PCERR(11), RBLOKD(268), DFCAPL(6,5), RBLOKE(5), ADUB(11,11),
3 DDSE(5,5), DTOLER(11), RBLOKF(132), IBLOKA(111), MDIMA, IBLOKB,
4 IWTSGN, IBLOKC(50), NF, IBLOKD, IWRITE, IBLOKE(6), JLAMP1,
5 IBLOKF(6), REGINT, FAILED, LBLOKA, NOBASE, LBLOKB(18), PIVALP(11)
a iblokc 30->50
a lblokb 14->18
COMMON /BLOK/ LBLOKC
DIMENSION CUI(11), CSTAR(11), DELT(10,10),
1 CCSUM(11,10), ADUM(6)
DIMENSION SUM(4), T(NMAX), E(NDIME,NLIN), SQRTW(NMAX),
1 TSTART(NINT), DELTAT(NINT), NT(NINT)
c+++++++++
holrth='base'
alplam(1,1)='lam1'; alplam(2,1)='lam2'; alplam(3,1)='lam3'
alplam(4,1)='lam4'; alplam(5,1)='lam5'
alplam(1,2)='alp1'; alplam(2,2)='alp2'; alplam(3,2)='alp3'
alplam(4,2)='alp4'; alplam(5,2)='alp5'
c++++++++
C-----------------------------------------------------------------------
C COMPUTE CORRELATION COEFFICIENTS AND STANDARD ERROR ESTIMATES
C USING FULL NORMAL EQUATIONS MATRIX (ADUB) WITH ALL PARAMETERS
C ASSUMED INDEPENDENT
C-----------------------------------------------------------------------
DO 110 K=1,JLAM
KK=K+JLAM
DO 120 J=1,K
JJ=J+JLAM
ADUB(J,K)=PALPHA(J)*PALPHA(K)*DDSE(J,K)
ADUB(J,KK)=-DFCAPL(K,J)*PALPHA(J)
ADUB(K,JJ)=-DFCAPL(J,K)*PALPHA(K)
ADUB(JJ,KK)=SE(J,K)
120 CONTINUE
DTOLER(K)=ADUB(K,K)*PRECIS
DTOLER(KK)=ADUB(KK,KK)*PRECIS
IF (NOBASE) GO TO 110
ADUB(K,NF)=-DFCAPL(NLIN,K)*PALPHA(K)
ADUB(KK,NF)=SE(K,NLIN)
110 CONTINUE
IF (NOBASE) GO TO 150
ADUB(NF,NF)=SUMWT
DTOLER(NF)=SUMWT*PRECIS
150 CALL PIVOT (ADUB,MDIMA,NF,.TRUE.,.TRUE.,.FALSE.,0.,
1 0.,CUI,DTOLER,CSTAR,LDUM,DUM,PIVALP,JLAM,DDUM,J)
SIGYY=SQRT(VAR/FLOAT(N-NF))
IF (.NOT.LDUM .OR. J.EQ.5) GO TO 172
DO 170 J=1,NF
IF (ADUB(J,J) .LE. 0.) GO TO 172
CSTAR(J)=1./SQRT(ADUB(J,J))
SIGMAP(J)=SIGYY/CSTAR(J)
170 CONTINUE
GO TO 174
172 FAILED=.TRUE.
5172 FORMAT (//51H0SINGULARITY IN INVERTING FULL LEAST SQUARES MATRIX,
1 30H - NO CORRELATIONS CALCULATED.,3X,24(2H**))
IF (PRLSQ) WRITE (IWRITE,5172)
RETURN
174 IF (INTER) RETURN
K=JLAM
DO 176 J=1,JLAM
PCERR(J)=100.*SIGMAP(J)/PLAM(J)
CUI(J)=ALPLAM(J,1)
K=K+1
PCERR(K)=100.*SIGMAP(K)/ABS(PALPHA(J))
CUI(K)=ALPLAM(J,2)
176 CONTINUE
IF (NOBASE) GO TO 178
CUI(NF)=HOLRTH
PCERR(NF)=100.*SIGMAP(NF)/ABS(PALPHA(JLAMP1))
5176 FORMAT (//5X,24HCORRELATION COEFFICIENTS/7X,10(3X,A4))
178 IF (.NOT.PRLSQ) GO TO 190
K=NF-1
WRITE (IWRITE,5176) (CUI(J),J=1,K)
DO 179 J=2,NF
DUM=CSTAR(J)
JK=J-1
DO 180 K=1,JK
180 DTOLER(K)=ADUB(J,K)*DUM*CSTAR(K)
5180 FORMAT (1X,A4,2X,10F7.3)
WRITE (IWRITE,5180) CUI(J),(DTOLER(K),K=1,JK)
179 CONTINUE
C-----------------------------------------------------------------------
C TO PREVENT A VERY SMALL COMPONENT (EFFECTIVELY A SECOND BASELINE)
C THAT IS HIGHLY CORRELATED WITH THE BASELINE FROM INTERFERING
C WITH THE CALCULATION OF THE WEIGHTS
C-----------------------------------------------------------------------
190 IF (IWTSGN.GT.0 .OR. ITYPE.NE.1) GO TO 200
DUM=1.E+20
DO 195 J=1,JLAM
IF (PLAM(J) .GT. DUM) GO TO 195
DUM=PLAM(J)
L=J
195 CONTINUE
IF (PCERR(L) .LT. 1000.) GO TO 200
FAILED=.TRUE.
5195 FORMAT (//47H SMALLEST LAMBDA HAS TOO HIGH A STANDARD ERROR.,
1 60H THIS COULD LEAD TO INACCURATE WEIGHTS - SOLUTION REJECTED.,
2 3X,11(2H**))
IF (PRLSQ) WRITE (IWRITE,5195)
200 IF (FAILED .OR. ITYPE.NE.4) GO TO 282
C-----------------------------------------------------------------------
C CALCULATION OF ENPHI (EQ. 3.25 OF E. M. L. BEALE FOR F RATIO
C CORRECTION FACTOR DUE TO NONLINEARITY)
C-----------------------------------------------------------------------
NUJ=JLAM
DO 210 J=1,JLAM
DUM=ADUB(J,J)
NUJ=NUJ+1
NUK=JLAM
DO 220 K=1,JLAM
NUK=NUK+1
DELT(J,K)=DUM*ADUB(K,K)+2.*ADUB(J,K)*ADUB(J,K)
DELT(J,NUK)=DUM*ADUB(K,NUK)+2.*ADUB(J,K)*ADUB(J,NUK)
DELT(NUJ,K)=ADUB(J,NUJ)*ADUB(K,K)+2.*ADUB(J,K)*ADUB(NUJ,K)
DELT(NUJ,NUK)=ADUB(J,NUJ)*ADUB(K,NUK)+ADUB(J,K)*ADUB(NUJ,NUK)
1 +ADUB(J,NUK)*ADUB(NUJ,K)
DO 230 L=1,4
230 SUM(L)=0.
NU=0
DO 240 L=1,NINT
CUI(1)=EXP(-PLAM(J)*DELTAT(L))
CUI(2)=EXP(-PLAM(K)*DELTAT(L))
CUI(3)=EXP(-PLAM(J)*TSTART(L))
CUI(4)=EXP(-PLAM(K)*TSTART(L))
II=NT(L)
DO 242 I=1,II
NU=NU+1
TNU=T(NU)
IF (REGINT) GO TO 244
SS=E(NU,J)
DDDUM=TNU*E(NU,K)
GO TO 246
244 SS=SQRTW(NU)*CUI(3)
DDDUM=TNU*SQRTW(NU)*CUI(4)
CUI(3)=CUI(3)*CUI(1)
CUI(4)=CUI(4)*CUI(2)
246 S=PALPHA(J)*TNU*SS
DDUM=PALPHA(K)*TNU*DDDUM
SUM(1)=SUM(1)-S*DDUM
SUM(2)=SUM(2)+S*DDDUM
SUM(3)=SUM(3)+SS*DDUM
SUM(4)=SUM(4)-SS*DDDUM
242 CONTINUE
240 CONTINUE
CCSUM(J,K)=SUM(1)
CCSUM(J,NUK)=SUM(2)
CCSUM(NUJ,K)=SUM(3)
CCSUM(NUJ,NUK)=SUM(4)
220 CONTINUE
IF (NOBASE) GO TO 210
S=0.
SS=0.
NU=0
DO 250 L=1,NINT
DUM=EXP(-PLAM(J)*DELTAT(L))
DDUM=EXP(-PLAM(J)*TSTART(L))
II=NT(L)
DO 252 I=1,II
NU=NU+1
IF (REGINT) GO TO 254
SSS=E(NU,J)
GO TO 256
254 SSS=DDUM*SQRTW(NU)
DDUM=DDUM*DUM
256 DDDUM=T(NU)*SSS*SQRTW(NU)
SS=SS-DDDUM
S=S+PALPHA(J)*T(NU)*DDDUM
252 CONTINUE
250 CONTINUE
CCSUM(NF,J)=S
CCSUM(NF,NUJ)=SS
210 CONTINUE
LDUM=.NOT.NOBASE
ENPHI=0.
DO 260 NU=1,N
TNU=T(NU)
DO 261 J=1,NLIN
261 ADUM(J)=SQRTW(NU)*EXP(-PLAM(J)*TNU)
NUJ=JLAM
DO 262 J=1,JLAM
NUJ=NUJ+1
SUM(1)=0.
SUM(2)=0.
NUK=JLAM
DO 264 K=1,JLAM
S=0.
IF (LDUM) S=SQRTW(NU)*ADUB(K,NF)
NUI=JLAM
DO 266 I=1,JLAM
NUI=NUI+1
CUI(NUI)=ADUM(I)
CUI(I)=-PALPHA(I)*TNU*CUI(NUI)
S=S+ADUB(K,I)*CUI(I)+ADUB(K,NUI)*CUI(NUI)
266 CONTINUE
SUM(1)=SUM(1)+CCSUM(K,J)*S
SUM(2)=SUM(2)+CCSUM(K,NUJ)*S
NUK=NUK+1
S=0.
IF (LDUM) S=SQRTW(NU)*ADUB(NUK,NF)
NUI=JLAM
DO 268 I=1,JLAM
NUI=NUI+1
S=S+ADUB(NUK,I)*CUI(I)+ADUB(NUK,NUI)*CUI(NUI)
268 CONTINUE
SUM(1)=SUM(1)+CCSUM(NUK,J)*S
SUM(2)=SUM(2)+CCSUM(NUK,NUJ)*S
264 CONTINUE
IF (NOBASE) GO TO 272
S=0.
IF (LDUM) S=SQRTW(NU)*ADUB(NF,NF)
NUI=JLAM
DO 270 I=1,JLAM
NUI=NUI+1
S=S+ADUB(NF,I)*CUI(I)+ADUB(NF,NUI)*CUI(NUI)
270 CONTINUE
SUM(1)=SUM(1)+CCSUM(NF,J)*S
SUM(2)=SUM(2)+CCSUM(NF,NUJ)*S
272 DUM=TNU*ADUM(J)
CSTAR(J)=PALPHA(J)*TNU*DUM-SUM(1)
CSTAR(NUJ)=-2.*(DUM+SUM(2))
262 CONTINUE
SUM(1)=0.
NUJ=JLAM
DO 280 J=1,JLAM
NUJ=NUJ+1
S=0.
SS=0.
NUK=JLAM
DO 281 K=1,JLAM
NUK=NUK+1
S=S+CSTAR(K)*DELT(J,K)+CSTAR(NUK)*DELT(J,NUK)
SS=SS+CSTAR(K)*DELT(NUJ,K)+CSTAR(NUK)*DELT(NUJ,NUK)
281 CONTINUE
SUM(1)=SUM(1)+CSTAR(J)*S+CSTAR(NUJ)*SS
280 CONTINUE
ENPHI=ENPHI+SUM(1)
260 CONTINUE
ENPHI=ENPHI*SIGYY**2/FLOAT(NF+2)
5282 FORMAT (36H STANDARD DEVIATION OF FIT = SIGYY =,1PE12.5//
1 6X,27HLAMBDA+-STD. ERROR PERCENT,7X,
2 26HALPHA+-STD. ERROR PERCENT)
5002 FORMAT (/1X)
282 IF (.NOT.PRLSQ) RETURN
WRITE (IWRITE,5002)
IF (.NOT. (FAILED .OR. ITYPE.NE.4)) WRITE (IWRITE,5283) ENPHI
5283 FORMAT ( 7H NPHI =,1PE10.3)
WRITE (IWRITE,5282) SIGYY
K=JLAM
DO 284 J=1,JLAM
K=K+1
5284 FORMAT (2(1PE14.4,2H+-,E10.3,0PF9.3))
WRITE (IWRITE,5284) PLAM(J),SIGMAP(J),PCERR(J),PALPHA(J),
1 SIGMAP(K),PCERR(K)
284 CONTINUE
5286 FORMAT (4X,8HBASELINE,19X,1PE14.4,2H+-,E10.3,0PF9.3)
IF (.NOT.NOBASE) WRITE(IWRITE,5286) PALPHA(JLAMP1),SIGMAP(NF),
1 PCERR(NF)
RETURN
END
C-------------------------- VERSION 1A OR 1B ---------------------------
C FUNCTION FISHER. CALCULATES PROBABILITY THAT FS(NU1,NU2) IS LESS
C THAN F FOR FISHER F-TEST.
C NU2 MUST BE .GE. 3
C SUMS UNTIL A TERM IS FOUND LESS THAN TOL*(CURRENT VALUE OF SUM).
C EQUATION NUMBERS REFER TO ABRAMOWITZ AND STEGUN, P. 946.
C NO COMMON BLOCKS NEEDED.
C-----------------------------------------------------------------------
C CALLS NO SUBPROGRAMS.
C-----------------------------------------------------------------------
FUNCTION FISHER (FS,NU1,NU2)
TOL=1.E-9
RNU1=FLOAT(NU1)
RNU2=FLOAT(NU2)
X=RNU2/(RNU2+RNU1*FS)
IF (X.LT.1. .OR. FS.LT.0.) GO TO 105
FISHER=0.
RETURN
105 IF (MOD(NU1,2) .NE. 0) GO TO 200
C-----------------------------------------------------------------------
C EVALUATES EQ. 26.6.4
C-----------------------------------------------------------------------
Q=1.
IF (NU1 .EQ. 2) GO TO 150
TERM=1.
TX=1.-X
TN=RNU2-2.
K=NU1-2
DO 110 J=2,K,2
TN=TN+2.
TERM=TERM*TN*TX/FLOAT(J)
Q=Q+TERM
IF (TERM .LE. TOL*Q) GO TO 150
110 CONTINUE
150 FISHER=1.-X**(.5*RNU2)*Q
RETURN
200 IF (MOD(NU2,2) .EQ. 0) GO TO 300
C-----------------------------------------------------------------------
C EVALUATES EQ. 26.6.8
C-----------------------------------------------------------------------
TH=ATAN(SQRT(RNU1*FS/RNU2))
CTH=COS(TH)
CC=CTH*CTH
STH=SIN(TH)
SS=STH*STH
TERM=1.
Q=1.
K=NU2-2
PR=1.
TN=0.
IF (K .LT. 3) GO TO 220
DO 210 J=3,K,2
TN=TN+2.
DUM=TN/FLOAT(J)
PR=PR*DUM
TERM=TERM*CC*DUM
Q=Q+TERM
IF (TERM .LE. TOL*Q) GO TO 220
210 CONTINUE
J=K
220 FISHER=TH+STH*CTH*Q
Q=0.
IF (NU1 .EQ. 1) GO TO 290
TN=TN+2.
PR=PR*TN
IF (J .EQ. K) GO TO 240
L=J+2
DO 230 J=L,K,2
TN=TN+2.
230 PR=PR*TN/FLOAT(J)
240 TN=RNU2-1.
Q=1.
TERM=1.
K=NU1-2
IF (K .LT. 3) GO TO 260
DO 250 J=3,K,2
TN=TN+2.
TERM=TERM*TN*SS/FLOAT(J)
Q=Q+TERM
IF (TERM .LE. TOL*Q) GO TO 260
250 CONTINUE
260 FISHER=FISHER-Q*PR*STH*CTH**NU2
290 FISHER=FISHER*.6366198
RETURN
C-----------------------------------------------------------------------
C EVALUATES EQ. 26.6.5
C-----------------------------------------------------------------------
300 TERM=1.
Q=1.
K=NU2-2
TN=RNU1-2.
DO 310 J=2,K,2
TN=TN+2.
TERM=TERM*TN*X/FLOAT(J)
Q=Q+TERM
IF (TERM .LE. TOL*Q) GO TO 320
310 CONTINUE
320 FISHER=SQRT(1.-X)*Q
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE RESIDU. CALCULATES WEIGHTED RESIDUALS FROM ASAVE(J,JLAM,1)
C AND ASAVE(J,JLAM,2), WHICH CORRESPOND TO ALPHA AND LAMBDA,
C RETURNS PUNCOR (APPROX. PROBABILITY THAT THE RESIDUALS ARE
C UNCORRELATED FOR LAGS 1 THRU 5) AND WEIGHTED AVERAGE IN
C PUNCOR(6), AND, IF PLOT=.TRUE., CALLS PLPRES TO PLOT WEIGHTED
C RESIDUALS ON PRINTER.
C-----------------------------------------------------------------------
C CALLS SUBPROGRAMS - FISHER, PLPRES
C-----------------------------------------------------------------------
SUBROUTINE RESIDU (JLAM,PLOT,T,Y,SQRTW,YLYFIT,NMAX,PUNCOR,ASAVE)
LOGICAL PLOT, NOBASE, WTED, LDUM, LBLOKA, LBLOKB
COMMON /BLOK/ RBLOKA(3535), IBLOKA(165), N, IWRITE, IBLOKB,
1 LINEPG, IBLOKC(11), LBLOKA(6), NOBASE, WTED, LBLOKB(30)
a ibloka 145->165
a lbloka 3->6, lblokb 25->30
DIMENSION ASAVE(6,5,6)
DIMENSION T(NMAX), Y(NMAX), SQRTW(NMAX), YLYFIT(NMAX), PUNCOR(6)
JLAMP1=JLAM+1
LDUM=.NOT.NOBASE
IF (LDUM) DUM=ASAVE(JLAMP1,JLAM,1)
DO 110 K=1,N
YLYFIT(K)=Y(K)
110 IF (LDUM) YLYFIT(K)=YLYFIT(K)-DUM*SQRTW(K)
DO 120 J=1,JLAM
DUM=ASAVE(J,JLAM,1)
DDUM=ASAVE(J,JLAM,2)
DO 125 K=1,N
125 YLYFIT(K)=YLYFIT(K)-DUM*EXP(-DDUM*T(K))*SQRTW(K)
120 CONTINUE
RSUMST=0.
IF (LDUM) GO TO 140
DO 130 K=1,N
130 RSUMST=RSUMST+YLYFIT(K)
140 RSUMEN=RSUMST
C-----------------------------------------------------------------------
C CALCULATE AUTOCOVARIANCES AND PROBABILITIES THAT RESIDUALS ARE
C UNCORRELATED (HAMILTON, P. 183)
C-----------------------------------------------------------------------
PUNCOR(6)=0.
DO 150 L=1,5
NEND=N-L
RSUMST=RSUMST-YLYFIT(NEND+1)
RSUMEN=RSUMEN-YLYFIT(L)
STBAR=RSUMST/FLOAT(NEND)
ENBAR=RSUMEN/FLOAT(NEND)
SS=0.
SE=0.
EE=0.
J=L
DO 152 K=1,NEND
S=YLYFIT(K)-STBAR
J=J+1
E=YLYFIT(J)-ENBAR
SS=SS+S*S
SE=SE+S*E
EE=EE+E*E
152 CONTINUE
RSQ=SE*SE/(SS*EE)
PUNCOR(L)=1.-FISHER(FLOAT(NEND-2)*RSQ/(1.-RSQ),1,NEND-2)
PUNCOR(6)=PUNCOR(6)+PUNCOR(L)/FLOAT(L)
150 CONTINUE
PUNCOR(6)=PUNCOR(6)/2.283333
IF (PLOT) CALL PLPRES (YLYFIT,NMAX,N,WTED,IWRITE,LINEPG)
RETURN
END
C-----------------------------------------------------------------------
C SUBROUTINE PLPRES. PLOTS RESIDUALS ON PRINTER WITH X-AXIS HORIZONTAL.
C NEEDS NO COMMON BLOCKS
C-----------------------------------------------------------------------
C CALLS NO EXTERNAL SUBPROGRAMS.
C-----------------------------------------------------------------------
SUBROUTINE PLPRES (YLYFIT,NMAX,N,WTED,IWRITE,LINEPG)
LOGICAL WTED
DIMENSION YLYFIT(NMAX), JCHAR(8), LINE(131), LABEL(6), BOUND(21),
1 LCHARJ(20), LINE1(20)
c++++++++
jchar(1)=42; jchar(2)=45;
jchar(3)=111 .shift. 8 + 117
jchar(4)= 32 .shift. 8 + 108
jchar(5)= 32 .shift. 8 + 32
jchar(6)= 32 .shift. 8 + 48
jchar(7)= 61 ; jchar(8)= 43
c++++++++
IF (LINEPG .GE. 17) GO TO 100
5100 FORMAT (/////47H NO PLOT OF RESIDUALS POSSIBLE BECAUSE LINEPG =,
1 I2,2X,40(2H**)//)
WRITE (IWRITE,5100) LINEPG
RETURN
100 MPLOT=(LINEPG-1)/16
MLINE=MIN0(20,(LINEPG-1)/MPLOT-3)
RMIN=1.E+20
RMAX=-1.E+20
DO 110 J=1,N
RMIN=AMIN1(RMIN,YLYFIT(J))
110 RMAX=AMAX1(RMAX,YLYFIT(J))
5200 FORMAT(55H1PLOT OF RESIDUALS MULTIPLIED BY SQUARE-ROOT OF WEIGHTS
1 ,21H - FOR BEST SOLUTION.,4X,15HMAX RESIDUAL=U=,1PE8.1,4X,
2 15HMIN RESIDUAL=L=,1PE8.1)
IF (WTED) WRITE (IWRITE,5200) RMAX,RMIN
5202 FORMAT (18H1PLOT OF RESIDUALS
1 ,21H - FOR BEST SOLUTION.,4X,15HMAX RESIDUAL=U=,1PE8.1,4X,
2 15HMIN RESIDUAL=L=,1PE8.1)
IF (.NOT.WTED) WRITE (IWRITE,5202) RMAX,RMIN
DELTA=(RMAX-RMIN)/FLOAT(MLINE-1)
BOUND(1)=RMAX+.5*DELTA
K=MLINE+1
DO 120 J=2,K
BOUND(J)=BOUND(J-1)-DELTA
120 IF (BOUND(J)*BOUND(J-1) .LT. 0.) JAXIS=J-1
LABEL(1)=-110
DO 130 J=2,6
130 LABEL(J)=LABEL(J-1)+20
K=MLINE-1
DO 140 J=2,K
LINE1(J)=JCHAR(5)
140 LCHARJ(J)=5
LINE1(1)=JCHAR(3)
LINE1(JAXIS)=JCHAR(6)
LINE1(MLINE)=JCHAR(4)
LCHARJ(1)=7
LCHARJ(JAXIS)=2
LCHARJ(MLINE)=7
NPOINT=0
DO 200 NPAGE=1,30
IF (NPAGE .GT. 1) WRITE(IWRITE,5999)
5999 FORMAT (1H1)
DO 210 NPLOT=1,MPLOT
NST=NPOINT+1
NEND=NPOINT+130
NPOINT=NEND
NLIM=MIN0(NEND,N)
DO 220 NLINE=1,MLINE
LCHAR=LCHARJ(NLINE)
LINE(1)=LINE1(NLINE)
BMAX=BOUND(NLINE)
BMIN=BOUND(NLINE+1)
K=1
DO 230 J=NST,NLIM
K=K+1
LINE(K)=JCHAR(LCHAR)
IF (YLYFIT(J).LT.BMAX .AND. YLYFIT(J).GE.BMIN) LINE(K)=JCHAR(1)
230 CONTINUE
K=NLIM-NST+2
IF (NLINE .NE. MLINE) GO TO 235
DO 232 J=21,K,20
232 IF (LINE(J) .NE. JCHAR(1)) LINE(J)=JCHAR(8)
5230 FORMAT (A2,130A1)
235 WRITE (IWRITE,5230) (LINE(J),J=1,K)
220 CONTINUE
DO 240 J=1,6
240 LABEL(J)=LABEL(J)+130
5240 FORMAT (3X,6(16X,I4)/)
WRITE (IWRITE,5240) LABEL
IF (NLIM .EQ. N) RETURN
210 CONTINUE
200 CONTINUE
5250 FORMAT (49H NO MORE THAN 30 PAGES OF RESIDUAL PLOTS ALLOWED.,
1 2X,38(2H**))
WRITE (IWRITE,5250)
RETURN
END
▶EOF◀