Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmdx85.for
There is 1 other file named bmdx85.for in the archive. Click here to see a list.
C NON-LINEAR LEAST SQUARES - MAIN PROGRAM APRIL 4, 1966
DIMENSION DATE(4),MT(19)
DIMENSION PLOW(20)
DIMENSION A(10000),P(100),P0(100),P1(100),F(180)
DOUBLE PRECISION GUESS,PROBLM,FINISH,AMON,AMOX,PP,PC
DOUBLE PRECISION TO,EP
C
C
DATA DATE/'MAY ','14, ','1969',' '/
DATA YES/'YES'/,ONO/'NO'/,GUESS/'PARAM '/,PROBLM/'PROBLM'/,
1FINISH/'FINISH'/,AMON/'MINMUM'/,AMOX/'MAXMUM'/
C
LLK=0
CALL USAGEB('BMDX85')
15 READ(5,1)PP,PC,NV,ND,NC,NP,TO,EP,MXIT,NF,IT,RE,IW
16 IF(NP.LE.0.OR.NP.GT.100)GO TO 700
IF(ND.LE.0.OR.ND.GT.NV)GO TO 800
NN=MIN0(NP,11)
LLK=LLK+1
1 FORMAT(2A6,2I4,I6,I4,2A6,3I4,A3,I4)
TOL=0.00001
CALL ATOF(TO,6,TOL)
EPS=0.00001
CALL ATOF(EP,6,EPS)
IF (MXIT.LT.1) MXIT=100
IF(IW.LT.0.OR.IW.GT.NV)IW=0
TOL=ABS(TOL)
IF(RE.NE.YES) RE=ONO
IF(IT.EQ.0) IT=5
IF(RE.EQ.YES) REWIND IT
NF=MAX0(1,NF)
NF=MIN0(NF,9)
WRITE(6,1000)DATE,PC,NV,ND,IW,NC,NP,TOL
WRITE (6,1001) EPS,MXIT,NF,IT,RE
IF(NF.GE.1 .AND. NF.LE.10) GO TO 11
NF = 1
WRITE(6,121)
121 FORMAT('0THE NUMBER OF VARIABLE FORMAT CARDS CANNOT BE .LT. 1 OR
1 .GT. 10. NF IS SET EQUAL TO 1.')
11 NF = NF * 18
READ (5,2) (F(I),I=1,NF)
2 FORMAT (18A4)
WRITE (6,1002) (F(I),I=1,NF)
1000 FORMAT('1BMDX85 - NON LINEAR LEAST SQUARES - REVISED ',4A4/
X/41H0HEALTH SCIENCES COMPUTING FACILITY, UCLA
X/32H0PROBLM CODE A6
X/32H0NUMBER OF VARIABLES I6
X/32H0INDEX OF THE DEPENDENT VARIABLE I6
X/'0INDEX OF THE WEIGHTING VARIABLE',I6
X/32H0NUMBER OF CASES I6
X/32H0NUMBER OF PARAMETERS I6
X/32H0TOLERANCE F10.6)
1001 FORMAT(32H0EPSILON F10.6
X/32H0MAXIMUM NUMBER OF ITERATIONS I6
X/32H0NUMBER OF VARIABLE FORMAT CARDS I6
X/32H0ALTERNATE INPUT TAPE NUMBER I6
X/32H0REWIND OPTION 3X,A3)
1002 FORMAT(16H0VARIABLE FORMAT 16X,24A4/(32X,24A4))
J=1
DO 3 L=1,NC
READ (IT,F) (P(I),I=1,NV)
C CALL TRANS(P,L,NV,SELECT)
DO 3 I=1,NV
A(J)=P(I)
3 J=J+1
NP1=NP+1
L1=1+NV*NC
L2=L1+NP1*NP1
C
C
IF(L2.GT.10000)GO TO 600
LLJ=0
DO 7 I=1,NP
P1(I)=1.E20
7 P0(I)=-1.E20
81 READ(5,25)PP,MT
25 FORMAT(A6,A2,18A4)
REWIND 1
WRITE(1,25)PP,MT
REWIND 1
IF(PP.NE.AMON.AND.PP.NE.AMOX)GO TO 40
IF(PP.EQ.AMON)READ(1,12)(P0(I),I=1,NN)
IF(PP.EQ.AMOX)READ(1,12)(P1(I),I=1,NN)
NQ=11
42 IF(NQ.GE.NP)GO TO 81
NR=NQ+1
NQ=MIN0(NQ+11,NP)
IF(PP.EQ.AMON)READ(5,12)(P0(I),I=NR,NQ)
IF(PP.EQ.AMOX)READ(5,12)(P1(I),I=NR,NQ)
GO TO 42
40 IF(PP.EQ.PROBLM)GO TO 22
IF(PP.EQ.FINISH) GO TO 24
IF(PP.EQ.GUESS) GO TO 18
WRITE(6,23)PP,MT
C
23 FORMAT(' NON CONTROL CARD (',A6,A2,18A4,')')
GO TO 81
22 READ(1,1)PP,PC,NV,ND,NC,NP,TO,EP,MXIT,NF,IT,RE,IW
GO TO 16
C
24 STOP
18 READ(1,12)(P(I),I=1,NN)
12 FORMAT(6X,11F6.0)
NQ=11
47 IF(NQ.GE.NP)GO TO 45
NR=NQ+1
NQ=MIN0(NQ+11,NP)
READ(5,12)(P(I),I=NR,NQ)
GO TO 47
45 CALL MINIZ(A,A(L1),P,NV,NC,NP1,TOL,EPS,EE,ND,MXIT,LLK,P0,P1,IW)
GO TO 81
600 WRITE (6,601)
601 FORMAT (' TOTAL AMOUNT OF DATA, N*(T+1) EXCEEDS 1500 VALUES, I.E.
1THE PROBLEM IS TOO LARGE'//)
C
602 GO TO 24
700 WRITE(6,701)
701 FORMAT(' NUMBER OF PARAMETERS EXCEEDS 100')
GO TO 602
800 WRITE(6,801)
801 FORMAT(' INDEX OF DEPENDENT VARIABLE EXCEEDS NUMBER OF VARIABLES')
GO TO 602
END
SUBROUTINE MINIZ(X,A,P,NV,NC,NP1,TOL,EPS,EE,ND,MXIT,LLK,P0,P1,IW)
C SUBROUTINE MINIZ FOR BMDX85 APRIL 4, 1966
DIMENSION A(NP1,NP1),X(NV,NC),P(100),FP(100),D(100),P0(100),P1(100
1),B(10)
DIMENSION PP0(100)
C
C
C
IF(IW.EQ.0)C=1.0
NP=NP1-1
MNAIL = 10
NAIL = 0
DO 557 I=1,NP
557 PP0(I)=P(I)
C
EPS1=5.*EPS+1.
WRITE(6,23)(P0(I),I=1,NP)
23 FORMAT('1MINIMA',15X,1P9E12.4/(22X,1P9E12.4))
WRITE(6,22)(P1(I),I=1,NP)
22 FORMAT('0MAXIMA',15X,1P9E12.4/(22X,1P9E12.4))
WRITE(6,195)
195 FORMAT('0ITERATION ERROR PARAMETERS'/14X,'MEAN'/
1' NO SHF SQUARE')
194 LLL=5
LIT=0
HU=NC
H=0.
LLLL=5
E0=1.0E30
DO 556 I=1,NP
556 P(I)=PP0(I)
FNC=NC-NP
29 DO 85 I=1,NP1
DO 85 J=1,I
85 A(I,J)=0.
DO 2 L=1,NC
C CALL FUN(F(L),FP,X(1,L),P,LLK)
CALL FUN(Q,FP,P,X(1,L))
FP(NP1)=X(ND,L)-Q
IF(IW.NE.0)C=X(IW,L)
IF(C.LT.(1.E-20)) GO TO 2
DO 1 I=1,NP1
FPIC=FP(I)*C
DO 1 J=1,I
1 A(I,J)=A(I,J)+FPIC*FP(J)
2 CONTINUE
EE=A(NP1,NP1)/FNC
6 IF(EE.LE.E0) GO TO 871
NAIL=NAIL+1
IF(NAIL.GT.MNAIL) GO TO 871
DO 872 I=1,NP
P(I)=PP0(I)+D(I)/2.
872 D(I)=D(I)/2.
GO TO 29
C
C
871 WRITE(6,198)LIT,NAIL,EE,(P(I),I=1,NP)
198 FORMAT(I6,1X,I2,1X,1P10E12.4/(22X,1P9E12.4))
NAIL=0
8191 LIT=LIT+1
CALL STEP(A,D,P,P0,P1,NP1,TOL)
IF(LIT.GT.MXIT) GO TO 200
IF(EE.EQ.0.) GO TO 470
PC=ABS((E0-EE)/EE)
E0=EE
IF(PC.LT.EPS) GO TO 100
LLL=LLLL
121 DO 28 I=1,NP
PP0(I)=P(I)
28 P(I)=P(I)+D(I)
GO TO 29
200 WRITE (6,201)
201 FORMAT(30H0THE PROCESS IS NOT CONVERGING)
GO TO 470
100 LLL=LLL-1
IF(LLL.NE.0) GO TO 121
470 DO 778 I=1,NP
778 FP(I)=SQRT(ABS(A(I,I)*EE))
WRITE (6,779) (FP(I),I=1,NP)
779 FORMAT('0ASYMPTOTIC STANDARD DEVIATIONS OF THE PARAMETERS'//
X(10X,1P10E12.4))
DO 841 I=1,NP
FP(I)=SQRT(ABS(A(I,I)))
DO 841 J=1,I
841 A(J,I)=A(I,J)
WRITE(6,774)
774 FORMAT('0ASYMPTOTIC CORRELATION MATRIX OF THE PARAMETERS')
L1=0
371 L0=L1+1
L1=MIN0(NP,L1+10)
WRITE(6,372)(L,L=L0,L1)
372 FORMAT('0',10I12)
DO 842 I=1,NP
K=0
DO 373 J=L0,L1
K=K+1
IF(FP(I).EQ.0.0.OR.FP(J).EQ.0.0.OR.A(I,J).EQ.0.0)GO TO 843
B(K)=A(I,J)/(FP(I)*FP(J))
GO TO 373
843 B(K)=0.0
373 CONTINUE
842 WRITE(6,374)I,(B(J),J=1,K)
374 FORMAT(I4,2X,10F12.5)
IF(L1.NE.NP)GO TO 371
1515 WRITE(6,553)
553 FORMAT('1CASE F Y-F STANDARD VARIABLES
1'/30X,'DEVIATION'/30X,'OF ESTIMATE')
DO 551 L=1,NC
SDEST=0.0
CALL FUN(Q,FP,P,X(1,L))
DO 24 I=1,NP
DO 24 J=1,NP
24 SDEST=SDEST+FP(I)*FP(J)*A(I,J)
SDEST=SQRT(SDEST*EE)
TT=X(ND,L)-Q
551 WRITE(6,552)L,Q,TT,SDEST,(X(I,L),I=1,NV)
552 FORMAT(1X,I4,10F12.5/(41X,7F12.5))
RETURN
END
SUBROUTINE STEP(A,D,P,P0,P1,NP1,TOL)
C SUBROUTINE STEP FOR BMDX85 APRIL 4, 1966
DIMENSION A(NP1,NP1),V(100),D(100),IN(100),P0(100),P1(100),P(100)
C
NP=NP1-1
DO 20 I=1,NP1
V(I)=SQRT(A(I,I))
DO 20 J=1,I
IF (V(I)*V(J).EQ.0.) GO TO 212
A(I,J)=A(I,J)/(V(I)*V(J))
GO TO 20
212 A(I,J)=0.
20 CONTINUE
DO 2 I=1,NP
2 IN(I)=0
8 K=0
R=0.
DO 1 I=1,NP
IF(IN(I).NE.0 .OR. A(I,I).LE.TOL) GO TO 1
Q=A(NP1,I)*A(NP1,I)/A(I,I)
IF(Q.LE.R) GO TO 1
K=I
R=Q
1 CONTINUE
IF(K.EQ.0) GO TO 3
C=-1.
4 DO 5 I=1,K
D(I)=A(K,I)
5 A(K,I)=0.
PP=D(K)
DO 6 I=K,NP1
D(I)=A(I,K)
6 A(I,K)=0.
D(K)=C
IN(K)=IN(K)+1
DO 7 I=1,NP1
IF (PP.EQ.0..OR.D(I).EQ.0.) GO TO 24
Y=D(I)/PP
GO TO 25
24 Y=0.
25 CONTINUE
DO 7 J=1,I
7 A(I,J)=A(I,J)-D(J)*Y
GO TO 8
3 R=1.
DO 9 I=1,NP
IF(IN(I).NE.1) GO TO 9
IF (V(I).EQ.0.) GO TO 15
D(I)=A(NP1,I)*V(NP1)/V(I)
C
C
C
16 IF (D(I).EQ.0.) GO TO 17
H=AMAX1((P0(I)-P(I))/D(I),(P1(I)-P(I))/D(I))
GO TO 18
15 D(I)=0.0
17 H=0.
18 CONTINUE
IF(H.GE.R) GO TO 9
R=H
K=I
9 CONTINUE
C=1.
IF(R.LE.TOL) GO TO 4
DO 10 I=1,NP
IF(IN(I).NE.1) GO TO 11
D(I)=D(I)*R
DO 21 J=1,I
IF ((V(I)*V(J)).EQ.0.0) GO TO 30
A(I,J)=-A(I,J)/(V(I)*V(J))
GO TO 21
30 A(I,J)=0.
21 A(J,I)=A(I,J)
GO TO 10
11 D(I)=0.
DO 22 J=1,NP
A(I,J)=0.
22 A(J,I)=0.
10 CONTINUE
RETURN
END
SUBROUTINE ATOF(A,N,F)
DIMENSION A(1)
LOGICAL BLANK
BLANK=.TRUE.
S=1.0
NUMB=0
TEN=1.0
DIV=1.0
DO 10 I=1,N
L=INTCHR(A,I)
IF(L.EQ.36) GO TO 10
BLANK=.FALSE.
IF(L.NE.38) GO TO 2
S=-1.0
GO TO 10
2 IF(L.NE.44) GO TO 4
TEN=10.0
GO TO 10
4 IF(L.GT.9) GO TO 9
NUMB=NUMB*10+L
DIV=DIV*TEN
9 CONTINUE
10 CONTINUE
IF(BLANK)RETURN
F=S*FLOAT(NUMB)/DIV
RETURN
END
FUNCTION INTCHR(STRING,N)
DIMENSION SEQ(50),STRING(1),EBCD(5)
DATA SEQ/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
X 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,
X 1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,
X 1HU,1HV,1HW,1HX,1HY,1HZ,1H ,1H+,1H-,1H*,
X 1H/,1H(,1H),1H,,1H.,1H',1H=,1H$,1H ,1H /
DATA EBCD/1H+,1H(,1H),1H',1H=/
CALL GETCHR(STRING,N,CHR)
IF (CHR.NE.SEQ(37)) GO TO 2
INTCHR = 36
GO TO 10
2 DO 1 I=1,48
IF(SEQ(I).EQ.CHR) GO TO 9
1 CONTINUE
I=51
IF(EBCD(1).EQ.CHR) I=38
IF(EBCD(2).EQ.CHR) I=42
IF(EBCD(3).EQ.CHR) I=43
IF(EBCD(4).EQ.CHR) I=46
IF(EBCD(5).EQ.CHR) I=47
9 INTCHR=I-1
10 RETURN
END
SUBROUTINE FUN(F,D,P,X)
DIMENSION D(1),P(1),X(1)