Google
 

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)