SUBROUTINE OBSWT(ND,ALPHA) C J.R.M. 6/27/91 C THIS IS A ROUTINE FOR DATA ENTRY C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*78 ALPHA COMMON /DATA/ YKP(151),XKP(151),SXKP(151),SYKP(151) DIMENSION XQ(200),YQ(200) C C NOTE: IF YOUR FORTRAN COMPILER CAN'T HANDLE AN "OPEN" SUCH AS C THAT BELOW, CHANGE: FILE=' ' TO: FILE=GENIN C AND USE BATCH FILE: RNGL.BAT FOR INPUT: i.e. PUT C RNGL YOURFILENAME on the command line and push ENTER. C OPEN(7,FILE=' ') READ(7,31) ALPHA WRITE(*,33) ALPHA 31 FORMAT(A78) 33 FORMAT(/2X,A78/) C C ## MENU 0: ENTER DATA FORMAT TYPE ## C WRITE(*,115) 115 FORMAT(5X,' ## MENU 0: DATA FORMAT TYPE ##'/) WRITE(*,121) 121 FORMAT(3X,'ENTER NUMBER TO MATCH INPUT DATA TYPE'/ 1 5X,'1: 2 COLUMNS: X-COLUMN Y-COLUMN (ALL STANDARD DEVIATIONS W 2ILL BE SET TO 1)' 3 4X,'2: 4 COLUMNS: X-COL SIG(X) Y-COL SIG(Y)'/ 4 5X,'3: 2 COLUMNS: X FOLLOWED BY SIG(X); Y FOLLOWED BY SIG(Y)'/ 5 5X,'4: 2 COLUMNS: X FOLLOWED BY WT(X); Y FOLLOWED BY WT(Y)'/) READ(*,*) ITYP IND = 0 C IF(ITYP.EQ.1) THEN DO 80 I = 1,200 READ(7,*,END=135) XKP(I),YKP(I) SXKP(I) = 1.D0 SYKP(I) = 1.D0 IND = IND + 1 80 CONTINUE 135 ND = IND C ELSEIF(ITYP.EQ.2) THEN DO 90 I = 1,200 READ(7,*,END=140) XKP(I),SXKP(I),YKP(I),SYKP(I) SXKP(I) = SXKP(I) + 1.D-36 SYKP(I) = SYKP(I) + 1.D-36 IND = IND + 1 90 CONTINUE 140 ND = IND C ELSEIF(ITYP.EQ.3.OR.ITYP.EQ.4) THEN DO 50 I = 1,400 READ(7,*,END=150) XQ(I),YQ(I) IND = IND + 1 50 CONTINUE 150 CONTINUE ND = IND/2 DO 60 I = 1,ND XKP(I) = XQ(I) YKP(I) = YQ(I) SXKP(I) = XQ(I+ND) + 1.D-36 SYKP(I) = YQ(I+ND) + 1.D-36 IF(ITYP.EQ.4) THEN SXKP(I) = 1.D0/DSQRT(SXKP(I)) SYKP(I) = 1.D0/DSQRT(SYKP(I)) ENDIF 60 CONTINUE ENDIF CLOSE(7) C RETURN END C SUBROUTINE PRINTOUT(SUMSXR,SUMSYR,S2,A,AA,SAK,SAKR,RC, 1 SDER,IMEND,NP) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(10),AA(10),SAK(10),RC(10,10),SAKR(10),SDER(10) C SAVE IF(IMEND.GT.1) GOTO 105 PRINT 400, (SDER(K),K =1,NP) PRINT 590, SUMSXR,SUMSYR,S2 PRINT 650, (A(K),K=1,NP) PRINT 670, (AA(K),K=1,NP) PRINT 680 IF(IMEND.EQ.0) GOTO 120 PRINT 910, (SAK(I),I=1,NP) PRINT 920, (SAKR(I),I=1,NP) PRINT 930 DO 940 I = 1,NP 940 PRINT 950, (RC(I,J),J=1,NP) PRINT 680 IF(IMEND.EQ.1) GOTO 120 C 105 WRITE(10,400) (SDER(K),K =1,NP) WRITE(10,590) SUMSXR,SUMSYR,S2 WRITE(10,650) (A(K),K=1,NP) WRITE(10,670) (AA(K),K=1,NP) WRITE(10,910) (SAK(I),I=1,NP) WRITE(10,920) (SAKR(I),I=1,NP) WRITE(10,930) DO 945 I = 1,NP 945 WRITE(10,950) (RC(I,J),J=1,NP) WRITE(10,680) C 400 FORMAT(4X,'(A/S)*(DS/DA):' 1P5E11.2) 590 FORMAT(/5X,6HS(X) = 1PD17.7/ 1 5X,6HS(Y) = D17.7/5X,14HS(X) + S(Y) = D20.11/) 650 FORMAT(/5X,16HPARAMETER VALUES /(2X,1P4D17.9)) 670 FORMAT(5X,10HINCREMENTS /(2X,1P4D17.5)) 680 FORMAT(/4X,19H *-*-*-*-*-*-*-*-* ) 910 FORMAT(//5X,33HSTANDARD DEVIATIONS OF PARAMETERS /(2X,1P4E17.5)) 920 FORMAT(//5X,42HRELATIVE STANDARD DEVIATIONS OF PARAMETERS / 1 (2X,1P4E17.5)) 930 FORMAT(//5X,28HPARAMETER CORRELATION MATRIX) 950 FORMAT (5X,1P5E13.4) C 120 RETURN END C SUBROUTINE PWGM(ALAM,DIVI,WX,WY,WXI,WYI,XKEEP,YKEEP,A,N,KM,IFPM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALAM(2),XKEEP(151),YKEEP(151),WX(151),WY(151),AP(2,2), 1 SDP(2,2),RSD(2,2),SP(2),SWM(2),SDX(151),SDY(151),WW(151),A(10), 1 WXI(151),WYI(151) C SAVE IF(IFPM.EQ.0) THEN SWW = 0.D0 SWX = 0.D0 SWY = 0.D0 SWXX = 0.D0 SWYY = 0.D0 SWXY = 0.D0 SWM(1) = 0.D0 SWM(2) = 0.D0 SWX2 = 0.D0 DO 737 I = 1,N IF(ALAM(1).NE.0.D0) THEN WW(I) = WYI(I) ELSE WW(I) = WXI(I) ENDIF SWW = SWW + WW(I) SWX = SWX + WW(I)*XKEEP(I) SWY = SWY + WW(I)*YKEEP(I) 737 CONTINUE SWWI = 1.D0/SWW XMN = SWX*SWWI YMN = SWY*SWWI DO 347 I = 1,N SDX(I) = XKEEP(I) - XMN SDY(I) = YKEEP(I) - YMN SWXX = SWXX + WW(I)*SDX(I)*SDX(I) SWYY = SWYY + WW(I)*SDY(I)*SDY(I) SWXY = SWXY + WW(I)*SDX(I)*SDY(I) SWX2 = SWX2 + WW(I)*XKEEP(I)**2 347 CONTINUE ALAM(2) = (SWYY/SWXX) SGG = DSIGN(1.D0,SWXY) GAM = (ALAM(1)*SWXX - SWYY)/(2.D0*SWXY) SURD = SGG*DSQRT(GAM*GAM + ALAM(1)) C C AP(K,J) K=1: PW; K=2: GM. J=1: INTERCEPT; J=2: SLOPE C IF((GAM.GE.0.AND.SURD.GE.0).OR.(GAM.LT.0.AND.SURD.LT.0)) THEN AP(1,2) = ALAM(1)/(GAM + SURD) ELSE AP(1,2) = -GAM + SURD ENDIF AP(2,2) = SGG*DSQRT(ALAM(2)) DO 356 KK = 1,KM AP(KK,1) = YMN - AP(KK,2)*XMN DO 367 JJ = 1,N SWM(KK) = SWM(KK) + WW(JJ)*(AP(KK,2)*SDX(JJ) - SDY(JJ))**2 367 CONTINUE SDP(KK,2) = DSQRT(DIVI*SWM(KK)/SWXX) SDP(KK,1) = SDP(KK,2)*DSQRT(SWX2*SWWI) RSD(KK,1) = SDP(KK,1)/DABS(AP(KK,1)) RSD(KK,2) = SDP(KK,2)/DABS(AP(KK,2)) C C CALCULATE CONVERGED OBJECTIVE FUNCTION VALUES C SP(KK) = 0.D0 DO 455 I = 1,N WF = 1.D0/(WY(I) + (AP(KK,2)*AP(KK,2))*WX(I)) SP(KK) = SP(KK) + WF*(YKEEP(I) - AP(KK,1) - AP(KK,2)*XKEEP(I))**2 455 CONTINUE 356 CONTINUE A(1) = AP(1,1) A(2) = AP(1,2) C DO 713 KK = 1,KM IF(KK.EQ.1) THEN PRINT 835 835 FORMAT(/10X,'******* PW SOLUTION *******'/) ELSE PRINT 837 837 FORMAT(/10X,'******* GM SOLUTION *******'/) ENDIF WRITE(*,641) ALAM(KK),SP(KK) 641 FORMAT(2X,'LAMBDA =',1PD14.5,5X,'SUM OF SQUARES =',1PD18.9) PRINT 656, AP(KK,1),AP(KK,2) 656 FORMAT(/5X,16HPARAMETER VALUES /(2X,1P2D18.9)) PRINT 912, SDP(KK,1),SDP(KK,2) 912 FORMAT(/5X,33HSTANDARD DEVIATIONS OF PARAMETERS / 1 (2X,1P2D14.5)) PRINT 922, RSD(KK,1),RSD(KK,2) 922 FORMAT(/5X,42HRELATIVE STANDARD DEVIATIONS OF PARAMETERS / 1 (2X,1P2D14.5)) PRINT 680 680 FORMAT(/4X,19H *-*-*-*-*-*-*-*-* ) C 713 CONTINUE C ELSEIF(IFPM.EQ.1) THEN DO 714 KK = 1,KM IF(KK.EQ.1) THEN WRITE(10,835) ELSE WRITE(10,837) ENDIF WRITE(10,641) ALAM(KK),SP(KK) WRITE(10,656) AP(KK,1),AP(KK,2) WRITE(10,912) SDP(KK,1),SDP(KK,2) WRITE(10,922) RSD(KK,1),RSD(KK,2) WRITE(10,680) 714 CONTINUE C ENDIF C RETURN END C SUBROUTINE DERXE(X,DX,I) C C THIS ROUTINE CALCULATES THE X DERIVATIVES C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*78 ALPHA DIMENSION P(2),X(151) COMMON /SECDER/ P,DEL,XX COMMON /IN/ ND,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA COMMON /EXTRA/ FAC C DO 40 K = 1,2 IDEL = (-1)**K XFAC = 0.D0 IF(XX.NE.0.D0) GOTO 30 I1 = I + 1 DO 10 J = I1,ND XX1 = X(J) IF(XX1.EQ.0.D0) GOTO 10 XFAC = DABS(XX1*FAC) GOTO 30 10 CONTINUE PRINT 20 20 FORMAT(/10X,14HERROR IN DERXE) STOP 30 CONTINUE DEL = (FAC*DABS(XX) + XFAC)*IDEL HOLD = XX XX = XX + DEL P(K) = SFUNC(ISD) 40 XX = HOLD DX = (P(2) - P(1))/(2.D0*DEL) RETURN END C SUBROUTINE DERAE(FA) C C THIS GETS THE PARTIALS OF F W.R.T. A C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*78 ALPHA DIMENSION P(2),FA(10) COMMON /IN/ ND,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA COMMON /EXTRA/ FAC C DO 40 J = 1,NP IF(NOUT.EQ.0) GOTO 20 DO 10 L = 1,NOUT IF(J.EQ.NS(L)) GOTO 40 10 CONTINUE 20 CONTINUE DO 30 K = 1,2 IDEL = (-1)**K FACS = IDEL *FAC IF(A(J).EQ.0) THEN DEL = FACS ELSE DEL = FACS*A(J) ENDIF HOLD = A(J) A(J) = A(J) + DEL P(K) = SFUNC(ISD) 30 A(J) = HOLD FA(J) = (P(2) - P(1))/(2.D0*DEL) 40 CONTINUE RETURN END C SUBROUTINE SCALC(IA,IB) C C PARAMETER ADJUSTMENT FOR DERIVATIVES C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*78 ALPHA COMMON /IN/ ND,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA COMMON /DATA/ Y(151),X(151),SX(151),SY(151) COMMON /EXTRY/ DEL1,DEL2,J,K,XPART COMMON /DON/ RX(151),RY(151),XKEEP(151),YPART,XPORT,XHOLD(151) COMMON /WTI/ WXI(151),WYI(151) C A1 = A(J) A2 = A(K) A(J) = A(J) + IA*DEL1 A(K) = A(K) + IB*DEL2 IF(NOUT.EQ.0) GOTO 20 DO 10 I = 1,NOUT IF(J.EQ.NS(I)) A(J) = A(J) - IA*DEL1 IF(K.EQ.NS(I)) A(K) = A(K) - IB*DEL2 10 CONTINUE 20 CONTINUE CALL XCHNGE(WXI,WYI) XPART = XPORT A(J) = A1 A(K) = A2 RETURN END C SUBROUTINE XCHNGE(WXI,WYI) C C THIS SUBROUTINE ADJUSTS THE X VECTOR C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*78 ALPHA COMMON /IN/ ND,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA COMMON /DATA/ Y(151),X(151),SX(151),SY(151) COMMON /SECDER/ P1,P2,DEL,XX COMMON /DON/ RX(151),RY(151),XKEEP(151),YPART,XPART,XHOLD(151) DIMENSION WXI(151),WYI(151) C YPART = 0.D0 XPART = 0.D0 TDLS = 2.D0/(DEL*DEL) DO 20 I = 1,ND XX = X(I) YY = Y(I) CALL DERXE(X,DX,I) DDX = DX SSC = SFUNC(ISD) SSL = P1 SECD = (DEL*DDX + SSL - SSC)*TDLS FNUM = (YY - SSC)*DDX*WYI(I) FDEN = (DDX**2 - (YY -SSC)*SECD)*WYI(I) + WXI(I) SUB = (XKEEP(I) - XX)*WXI(I) DELX = (FNUM + SUB)/FDEN XHOLD(I) = XX + DELX DELXX = XKEEP(I) - XHOLD(I) RX(I) = DELXX XX = XHOLD(I) TEMP = YY - SFUNC(ISD) RY(I) = TEMP IF(WYI(I)/WXI(I).LT.1.D20) GOTO 10 RY(I) = 0.D0 TEMP = 0.D0 10 CONTINUE YPART = YPART + TEMP**2*WYI(I) XPART = XPART + DELXX**2*WXI(I) 20 CONTINUE C RETURN END C SUBROUTINE PEARSN(X,Y,N,R,PROB,Z) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (TINY=1.D-20) DIMENSION X(N),Y(N) AX=0.D0 AY=0.D0 DO 11 J=1,N AX=AX+X(J) AY=AY+Y(J) 11 CONTINUE AX=AX/N AY=AY/N SXX=0.D0 SYY=0.D0 SXY=0.D0 DO 12 J=1,N XT=X(J)-AX YT=Y(J)-AY SXX=SXX+XT**2 SYY=SYY+YT**2 SXY=SXY+XT*YT 12 CONTINUE R=SXY/DSQRT(SXX*SYY) C Z=0.5*DLOG(((1.+R)+TINY)/((1.-R)+TINY)) C DF=N-2 C T=R*DSQRT(DF/(((1.-R)+TINY)*((1.+R)+TINY))) C PROB=BETAI(0.5*DF,0.5,DF/(DF+T**2)) C PROB=ERFCC(ABS(Z*DSQRT(N-1.D0))/1.414213562) RETURN END C SUBROUTINE DMATIN(A,MP,CINV,ISD) C C CARRY OUT DOUBLE PRECISION MATRIX INVERSION C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(10,10),CINV(10,10),INDX(10) C DO 121 I = 1,MP DO 111 J = 1,MP CINV(I,J) = 0.D0 111 CONTINUE CINV(I,I) = 1.D0 121 CONTINUE CALL LUDCMP(A,MP,10,INDX,D,ISD) DO 131 J = 1,MP CALL LUBKSB(A,MP,10,INDX,CINV(1,J)) 131 CONTINUE C IF(ISD.EQ.1) WRITE(*,*) '!!!** SINGULAR MATRIX - Beware *!!!' RETURN END C SUBROUTINE LUDCMP(A,N,NP,INDX,D,ISD) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NMAX=18,TINY=1.0D-20) DIMENSION A(NP,NP),INDX(*),VV(NMAX) C D=1.D0 ISD = 0 DO 12 I=1,N AAMAX=0.D0 DO 11 J=1,N IF (DABS(A(I,J)).GT.AAMAX) AAMAX=DABS(A(I,J)) 11 CONTINUE IF (AAMAX.EQ.0.D0) THEN write(*,*) 'Singular matrix.' AAMAX = 1.D0 ISD = 1 ENDIF VV(I)=1.D0/AAMAX 12 CONTINUE DO 19 J=1,N IF (J.GT.1) THEN DO 14 I=1,J-1 SUM=A(I,J) IF (I.GT.1)THEN DO 13 K=1,I-1 SUM=SUM-A(I,K)*A(K,J) 13 CONTINUE A(I,J)=SUM ENDIF 14 CONTINUE ENDIF AAMAX=0.D0 DO 16 I=J,N SUM=A(I,J) IF (J.GT.1)THEN DO 15 K=1,J-1 SUM=SUM-A(I,K)*A(K,J) 15 CONTINUE A(I,J)=SUM ENDIF DUM=VV(I)*DABS(SUM) IF (DUM.GE.AAMAX) THEN IMAX=I AAMAX=DUM ENDIF 16 CONTINUE IF (J.NE.IMAX)THEN DO 17 K=1,N DUM=A(IMAX,K) A(IMAX,K)=A(J,K) A(J,K)=DUM 17 CONTINUE D=-D VV(IMAX)=VV(J) ENDIF INDX(J)=IMAX IF(J.NE.N)THEN IF(A(J,J).EQ.0.D0) A(J,J)=TINY DUM=1.D0/A(J,J) DO 18 I=J+1,N A(I,J)=A(I,J)*DUM 18 CONTINUE ENDIF 19 CONTINUE IF(A(N,N).EQ.0.D0) A(N,N)=TINY RETURN END C SUBROUTINE LUBKSB(A,N,NP,INDX,B) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(NP,NP),INDX(*),B(NP) C II=0 DO 12 I=1,N LL=INDX(I) SUM=B(LL) B(LL)=B(I) IF (II.NE.0) THEN DO 11 J=II,I-1 SUM=SUM-A(I,J)*B(J) 11 CONTINUE ELSE IF (DABS(SUM).GT.1.D-200) THEN II=I ENDIF B(I)=SUM 12 CONTINUE DO 14 I=N,1,-1 SUM=B(I) IF(I.LT.N)THEN DO 13 J=I+1,N SUM=SUM-A(I,J)*B(J) 13 CONTINUE ENDIF B(I)=SUM/A(I,I) 14 CONTINUE RETURN END