C DISTRIBUTED BY
C MARKETING SCIENCE INSTITUTE
C PHILADELPHIA,PENNA. 19104
C
C THIS PROGRAM WAS WRITTEN BY:
C
C DR. J. B. KRUSKAL
C BELL TELEPHONE LABORATORIES
C MURRAY HILL, N. J.
C
C
C ADAPTED FOR IBM S/360 BY:
C
C F. J. CARMONE, JR.
C MARKETING SCIENCE INSTITUTE
C PHILADELPHIA, PA.
C NEW ADDRESS:
C DEPARTMENT OF ECONOMICS
C UNIV OF WATERLOO
C WATERLOO, ONTARIO CANADA
C
C JULY 1968
C
C This Version contains a minor change by M. Birnbaum 1970
C to allow the program to create a file ("punch") of XHATS for
C analysis of variance. An example application of this aspect
C of the program is given in Birnbaum & Veit (1974).
C
C REFERENCES:
C Kruskal, J. B., & Carmone, F. J. (1969). MONANOVA: A FORTRAN IV
C program for monotone analysis of variance. Behavioral Science,
C 14, 165-166.
C
C KRUSKAL, J. B. 'ANALYSIS OF FACTORIAL EXPERIMENTS BY
C ESTIMATING MONOTONE TRANSFORMATIONS OF THE DATA',
C JOURNAL OF ROYAL STATISTICAL SOCIETY,SERIES B, VOL. 27,
C NO. 2 (1965), PP. 251-63.
C
C Birnbaum, M. H., & Veit, C. T. (1974). Scale-free tests
C of an averaging model for the size-weight illusion.
C Perception & Psychophysics, 16, 276-282.
C
C
C
C GENERAL REMARKS:
C
C THIS PROGRAM CONTAINS THE FOLLOWING MAJOR ROUTINES:
C
C OPTEX MAIN ROUTINE
C OPSOL CALCULATES STEP SIZE AND OTHER VARIABLES FOR TERMINATIO
C CALC CALCULATES STRESS AND READS IN DATA
C MFIT PERFORMS WEIGHTED LEAST SQUARES REGRESSION
C CCACT READS AND INTERPERTS CONTROL CARDS
C SORT SORTS ARRAYS BY USING POINTERS; SIMPLE, BUT VERY FAST
C PLOTR PLOTS TWO MODELS VERSUS ORIGINAL DATA
C
C
C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE.
C
C ALL NORMAL INPUT-OUTPUT HANDLED BY CALC AND CCACT.
C
C MFIT HAS EMERGENCY DIAGNOSTIC OUTPUT.
C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS CONTROLLED BY
C THESE VARIABLES. LREAD, LPRINT, LPUNCH, LSCRAT.
C UNIT NUMBERS 5, 6, 7, 8 HAVE BEEN USED RESPECTIVELY.
C TO CHANGE THESE ASSIGNMENTS, IT SUFFICES TO CHANGE THE VALUES
C FOR THE FOUR VARIABLES SET IN THE MAIN PROGRAM.
C
C NOTE THAT THE SCRATCH UNIT IS USED IN A VERY MINOR WAY.
C IT IS USED ONLY BY CCACT
C MANY INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THING
C
C
C ****************************************************************
C The line below instructs the LS FORTRAN compiler to NOT use Fortran 77
C but instead to use an earlier FORTRAN. OPTIONS /NOF77 means use
C old FORTRAN, F66. You may need to make a similar modification,
C or you may need to delete the line.
C ***************************************************************
OPTIONS /NOF77
C
C
C
Program MONANOV
INTEGER LPAR(25)
REAL PAR(25)
C REAL*4 ATAB(4,22)
REAL ATAB(4,22)
C INTEGER*4 TAB(4,22),MTAB(4)
INTEGER TAB(4,22),MTAB(4)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
COMMON/OPTEX2/LPAR
COMMON/OPTEX3/MTAB,TAB
EQUIVALENCE(LPAR,PAR)
EQUIVALENCE(TAB,ATAB)
OPEN (5,FILE='MONANOVA.DAT', STATUS='old')
OPEN (6,FILE='monanova.out',STATUS='new')
OPEN (7,FILE='MONANOVA.XHATS',STATUS='NEW')
OPEN (8,FILE='OUTT8',STATUS='UNKNOWN')
DATA LREAD,LPRINT,LPUNCH,LSCRAT/5,6,7,8/
DATA MTAB/1,18,20,22/
DATA TAB/'STRM',2,2,0,'SRAT',2,1,0,'ITER',1,3,0,'COSA',2,5,0,
+'ACSA',2,6,0,'SRTA',2,4,0,'DATA',3,7,2,'CONF',3,7,3,'PARA',3,7,
+4,'COMP',3,7,5,'STOP',3,7,6,'TITL',3,7,7,'CUTO',2,8,0,'TIES',5,
+1,2,'CARD',3,10,1,'NOCA',3,10,2,'PRIN',5,1,3,'PRIM',3,9,1,
+'SECO',3,9,2,'NOIN',3,11,1,'INPU',3,11,2,0,0,0,0/
CALL OPTEX
STOP
END
C
SUBROUTINE OPTSOL(SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,
+LCONSW,STRESS)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
C INITIALIZE UNIVERSAL VARIABLES
100 ITNO=0
GRMULT=0.0
GRMAG=1.0
SRATAV=1.2
COSAV=0.0
ACSAV=0.2
LXMIN=1
C INITIALIZE CONFIGURATION AND GRADIENTS
200 GOTO (210,220),LCONSW
C CREATE INITIAL CONFIGURATION
210 CALL CALC(2,0.0,0.0,0.0)
C NORMALIZE CONFIGURATION AND PRINT INPUT
220 CALL CALC(6,0.0,0.0,0.0)
C PERFORM ITERATION
C
C
C STRESS=1.0
C
C
300 STRLST=STRESS
CALL CALC(4,STRESS,GRGR,GRGL)
305 GLMAG=GRMAG
GRMAG=SQRT(GRGR)
IF(STRESS)9999,320,310
310 IF(GRGR)9999,320,330
320 LXMIN=2
GOTO 1000
330 IF(ITNO)9999,340,350
340 SRAT=1.2
CAGRGL=0.0
GOTO 360
350 SRAT=STRLST/STRESS
CAGRGL=GRGL/(GRMAG*GLMAG)
C GO TO UPDATING SECTION
360 GOTO 1000
380 CALL CALC(5,GRMULT,0.0,0.0)
ITNO=ITNO+1
GOTO 300
C TERMINATION HAS OCCURRED. PERFORM TERMINAL ACTIONS.
C NORMALIZE CONFIGURATION.
400 CALL CALC(6,0.0,0.0,-50.0)
C AGAIN CALCULATE STRESS TO GUARANTEE AGAINST ROUNDING
C EFFECTS DURING NORMALIZATION
CALL CALC(4,STRESS,GRGR,GRGL)
RETURN
9999 CALL OPTERR('OPTSOL')
STOP
C UPDATE VARIOUS VARIABLES. DECIDE WHETHER TO TERMINATE.
1000 IF(ITNO)9999,1010,1100
1010 WRITE(LPRINT,1)
1 FORMAT(' HISTORY OF COMPUTATION.'/1X)
WRITE(LPRINT,2)
2 FORMAT(' ITERATION STRESS SRAT SRTAVG CAGRGL COSAV',
+' ACSAV GRMAG GRMULT STEP')
LFOOT=1
1100 GOTO (1110,1300),LXMIN
1110 COSAV=CAGRGL*COSAVW+COSAV*(1.0-COSAVW)
ACSAV=ABS(CAGRGL)*ACSAVW+ACSAV*(1.0-ACSAVW)
SRATAV=(SRAT**SRTAVW)*(SRATAV**(1.0-SRTAVW))
IF(ITNO)9999,1200,1250
1200 GRMULT=25.0*STRESS
STEP=GRMULT*GRMAG
GOTO 1300
1250 ANG=4.0**COSAV
TEMP1=1.0+(AMIN1(1.0,1.0/SRATAV))**5
TEMP3=ABS(COSAV)
TEMP2=1.0+((ACSAV-TEMP3)/(1.0-TEMP3))**3
BIAS=1.6/(TEMP1*TEMP2)
GOODLK=SQRT(AMIN1(1.0,SRAT))
STEP=STEP*ANG*BIAS*GOODLK
GRMULT=STEP/GRMAG
1300 WRITE(LPRINT,3)ITNO,STRESS,SRAT,SRATAV,CAGRGL,COSAV,ACSAV,GRMAG,
+GRMULT,STEP
3 FORMAT(I10,F7.3,2F8.4,3F7.3,3F9.5)
2000 IF(STRESS)9999,2010,2020
2010 LSTOP=1
GOTO 3000
2020 IF(GRMAG)9999,2030,2040
2030 LSTOP=2
GOTO 3000
2040 TEMP1=0.5*(1.0+SRATST)
TEMP2=ABS(TEMP1-1.0)
IF(ABS(SRAT-TEMP1)-TEMP2)2050,2050,2070
2050 IF(ABS(SRATAV-TEMP1)-TEMP2)2060,2060,2070
2060 LSTOP=3
GOTO 3000
2070 IF(STRESS-STRMIN)2080,2080,2090
2080 LSTOP=4
GOTO 3000
2090 IF(NOIT-ITNO)9999,2100,2110
2100 LSTOP=5
GOTO 3000
2110 GOTO (380,400),LFOOT
3000 GOTO (3010,3020),LFOOT
3010 LFOOT=2
WRITE(LPRINT,2)
WRITE(LPRINT,30)
30 FORMAT(1X)
3020 GOTO (3030,3040,3050,3060,3070),LSTOP
3030 WRITE(LPRINT,31)
31 FORMAT(' ZERO STRESS WAS REACHED')
GOTO 2020
3040 WRITE(LPRINT,32)
32 FORMAT(' MINIMUM WAS ACHIEVED')
GOTO 2040
3050 WRITE(LPRINT,33)
33 FORMAT(' APPROXIMATE MINIMUM WAS REACHED.')
GOTO 2070
3060 WRITE(LPRINT,34)
34 FORMAT(' SATISFACTORY STRESS WAS REACHED')
GOTO 2090
3070 WRITE(LPRINT,35)
35 FORMAT(' MAXIMUM NUMBER OF ITERATIONS WERE USED')
GOTO 2110
END
SUBROUTINE OPTERR(NAME)
C REAL*8 NAME
DOUBLE PRECISION NAME
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
WRITE(LPRINT,1)NAME
1 FORMAT(' MONANOVA ERROR STOP FROM ',A6)
STOP
END
C
SUBROUTINE CALC(LLL,AAA,BBB,CCC)
C FOR CALC NONMETRIC ANALYSIS OF FACTORIAL DESIG
DIMENSION DATA(501),INDEX(500,4),X(500),XHAT(500)
DIMENSION IPOINT(501),LTIES(501)
DIMENSION CON(4,100),GR(4,100),GL(4,100),GR2(4,100)
DIMENSION II(4),IT(4)
DIMENSION LPAR(25),FMTCON(20),FMTDAT(20)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
COMMON/OPTEX2/LPAR
EQUIVALENCE(PARVAL,LPARVA)
EQUIVALENCE(CUTOFF,LPAR(8)),(MFITCO,LPAR(9)),(LPUNSW,LPAR(10))
EQUIVALENCE(LPRINS,LPAR(11)),(LCONSW,LPAR(12))
DATA IXSBJ/0/
C ********************************************************************
C READ CONFIGURATION
GOTO (1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,11000),
+LLL
1000 READ(LREAD,100)KK,(II(K),K=1,KK)
WRITE(LPRINT,100)KK,(II(K),K=1,KK)
100 FORMAT(24I3)
READ(LREAD,101)FMTCON
WRITE(LPRINT,105)FMTCON
101 FORMAT(20A4)
105 FORMAT(2X,20A4)
DO 1010K=1,KK
IIK=II(K)
1010 READ(LREAD,FMTCON)(CON(K,I),I=1,IIK)
RETURN
C ********************************************************************
C CREATE INITIAL CONFIGURATION
2000 DO 2005K=1,KK
IIK=II(K)
DO 2005I=1,IIK
CON(K,I)=0
GR(K,I)=0
GL(K,I)=0
GR2(K,I)=0
2005 CONTINUE
DO 2010L=1,LDATA
M=IPOINT(L)
TEMP=DATA(M)
DO 2010K=1,KK
I=INDEX(M,K)
CON(K,I)=CON(K,I)+TEMP
2010 GR(K,I)=GR(K,I)+1.0
DO 2020K=1,KK
IIK=II(K)
DO 2020I=1,IIK
IF(GR(K,I).EQ.0.0)GOTO2020
CON(K,I)=CON(K,I)/GR(K,I)
GR(K,I)=0
2020 CONTINUE
RETURN
C ********************************************************************
C READ DATA
3000 CONTINUE
IF(IXSBJ.GT.0)GOTO3002
C IPUNQ NOT EQUAL TO ZERO MEANS TO PUNCH XHAT
READ(LREAD,100)KK,(II(K),K=1,KK),LREPLI,IPUNQ
C TO PUNCH XHAT SIMPLY TYPE THE NUMBER 1 ON THE PROBLEM CARD AFTER REPLI
WRITE(LPRINT,100)KK,(II(K),K=1,KK),LREPLI
READ(LREAD,101)FMTDAT
WRITE(LPRINT,105)FMTDAT
3002 IIS=0
IIP=1
DO 3010K=1,KK
IIS=IIS+II(K)
IIP=IIP*II(K)
3010 IT(K)=1
NOD=IIP*LREPLI
3020 READ(LREAD,FMTDAT)(DATA(L),L=1,NOD)
L=0
LR=0
DATMAX=-1.0E+20
DATMIN=-DATMAX
DO 3090LA=1,NOD
TDATA=DATA(LA)
IF(TDATA-CUTOFF)3045,3045,3030
3030 L=L+1
DATMAX=AMAX1(DATMAX,TDATA)
DATMIN=AMIN1(DATMIN,TDATA)
DATA(L)=TDATA
DO 3040K=1,KK
3040 INDEX(L,K)=IT(K)
3045 LR=LR+1
IF(LR-LREPLI)3090,3050,9993
3050 LR=0
IT(KK)=IT(KK)+1
DO 3080KC=1,KK
K=KK+1-KC
IF(IT(K)-II(K)-1)3090,3060,9993
3060 IT(K)=1
IF(K-1)9993,3100,3070
3070 IT(K-1)=IT(K-1)+1
3080 CONTINUE
3090 CONTINUE
3100 LDATA=L
LDATAP=LDATA+1
DATA(LDATAP)=1.0E+20
DO 3110L=1,LDATAP
3110 IPOINT(L)=L
CALL SORT(IPOINT(1),LDATA,DATA(1),+1)
RETURN
9993 CALL OPTERR('MOD1RD')
STOP
C ********************************************************************
C CALCULATE STRESS AND GRADIENT
4000 XBAR=0.0
DO 4010L=1,LDATA
TEMP=0.0
M=IPOINT(L)
DO 4005K=1,KK
I=INDEX(M,K)
4005 TEMP=TEMP+CON(K,I)
X(M)=TEMP
4010 XBAR=XBAR+TEMP
XBAR=XBAR/FLOAT(LDATA)
K=1
L1=1
M=IPOINT(1)
TEMP2=DATA(M)
DO 50L=1,LDATA
TEMP1=TEMP2
M=IPOINT(L+1)
TEMP2=DATA(M)
IF(TEMP1-TEMP2)30,20,30
20 K=K+1
GOTO 50
30 L2=L1+K-1
DO 40LA=L1,L2
M=IPOINT(LA)
LTIES(M)=K
40 CONTINUE
L1=L2+1
K=1
50 CONTINUE
CALL MFIT(X(1),XHAT(1),LTIES(1),LDATA,MFITCO,IPOINT)
U=0
T=0
DO 4030K=1,KK
IIK=II(K)
DO 4030I=1,IIK
GR(K,I)=0
4030 GR2(K,I)=0
DO 4050L=1,LDATA
M=IPOINT(L)
TEMP1=X(M)-XHAT(M)
TEMP2=X(M)-XBAR
U=U+TEMP1**2
T=T+TEMP2**2
DO 4040K=1,KK
I=INDEX(M,K)
GR(K,I)=GR(K,I)+TEMP1
4040 GR2(K,I)=GR2(K,I)+TEMP2
4050 CONTINUE
U=SQRT(U)
T=SQRT(T)
STRESS=U/T
GRGR=0
GRGL=0
IF(U)9994,4080,4060
4060 COEF1=1.0/(U*T)
COEF2=-U/T**3
DO 4070K=1,KK
IIK=II(K)
DO 4070I=1,IIK
GR(K,I)=-COEF1*GR(K,I)-COEF2*GR2(K,I)
GRGR=GRGR+GR(K,I)**2
4070 GRGL=GRGL+GR(K,I)*GL(K,I)
4080 CONTINUE
AAA=STRESS
BBB=GRGR
CCC=GRGL
RETURN
9994 CALL OPTERR('MOD1CA')
RETURN
C ********************************************************************
C MOVE TO NEW CONFIGURATION
5000 GRMULT=AAA
DO 5010K=1,KK
IIK=II(K)
DO 5010I=1,IIK
CON(K,I)=CON(K,I)+GRMULT*GR(K,I)
GL(K,I)=GR(K,I)
GR(K,I)=0
5010 GR2(K,I)=0
RETURN
C ********************************************************************
C NORMALIZE CONFIGURATION
6000 SMULT=0
C GO TO 6040
IF(LPRINS.NE.2)GOTO6040
IF(CCC.LT.0.0)GOTO6040
WRITE(LPRINT,602)
602 FORMAT('1',' SEQ. NO. DATA SUBSCRIPTS'/)
DO 6050L=1,LDATA
WRITE(LPRINT,601)L,DATA(L),(INDEX(L,K),K=1,KK)
601 FORMAT(I6,4X,F10.5,4I5)
6050 CONTINUE
IF(LCONSW.EQ.1)GOTO6040
WRITE(LPRINT,603)
603 FORMAT(///5X,' CONFIGURATION'/)
DO 8011K=1,KK
IIK=II(K)
WRITE(LPRINT,8005)IIK,(CON(K,I),I=1,IIK)
8011 CONTINUE
6040 DO 6020K=1,KK
AV=0
IIK=II(K)
DO 6005I=1,IIK
GL(K,I)=0
6005 AV=AV+CON(K,I)
AV=AV/FLOAT(IIK)
DO 6010I=1,IIK
CON(K,I)=CON(K,I)-AV
6010 SMULT=SMULT+CON(K,I)**2
6020 CONTINUE
SMULT=SQRT(FLOAT(IIS)/SMULT)
DO 6030K=1,KK
IIK=II(K)
DO 6030I=1,IIK
6030 CON(K,I)=SMULT*CON(K,I)
RETURN
C ********************************************************************
C READ PARAMETERS
7000 READ(LREAD,103)CUTOFF,MFITCO
103 FORMAT(2F10.4)
RETURN
C ********************************************************************
C PRINT, PUNCH, AND DRAW OUTPUT
8000 CALL OPTPRI
DO 8010K=1,KK
IIK=II(K)
WRITE(LPRINT,8005)IIK,(CON(K,I),I=1,IIK)
8005 FORMAT(1X,I3,10F7.3/(4X,10F7.3))
8010 CONTINUE
WRITE(LPRINT,8015)
8015 FORMAT('-SEQ.NO.',6X,'DATA',9X,'X',10X,'XHAT',5X,
+'SUBSCRIPTS...'/)
C *****************************************************************
DO 8016L=1,LDATA
8016 WRITE(LPRINT,8017)L,DATA(L),X(L),XHAT(L),(INDEX(L,K),K=1,KK)
8017 FORMAT(1X,I5,3X,3F12.4,4I6)
C*******************************************************************
C ***************** MODIFIED HERE BY M BIRNBAUM TO GIVE XHATS ******
IF(IPUNQ.EQ.0)GOTO8999
KLIM=KK-1
IPROD=1
DO 8887I=1,KLIM
8887 IPROD=IPROD*II(I)
DO 8889J=1,IPROD
IL=II(KK)*(J-1)+1
IU=IL+II(KK)-1
WRITE(LPRINT,8888)J,(XHAT(L),L=IL,IU)
8889 WRITE(LPUNCH,8888)J,(XHAT(L),L=IL,IU)
8888 FORMAT(1X,I4,1X,8F8.4/8F8.4)
C ********************************************************************
8999 WRITE(LPRINT,8001)
8001 FORMAT('1')
CALL PLOTR(DATA,XHAT,X,LDATA,6,+1)
IF(LPUNSW.EQ.2)RETURN
CALL OPTPUN
WRITE(LPUNCH,106)
106 FORMAT(2X,' (2X,10F7.3) ')
WRITE(LPUNCH,8020)KK,(II(K),K=1,KK)
8020 FORMAT(20I3)
DO 8040K=1,KK
IIK=II(K)
WRITE(LPUNCH,8030)(CON(K,I),I=1,IIK)
8030 FORMAT(2X,10F7.3)
8040 CONTINUE
RETURN
C ********************************************************************
C INITIAL ACTIONS
9000 CUTOFF=-1.23E+20
MFITCO=1
LPUNSW=2
RETURN
C ********************************************************************
10000 CONTINUE
RETURN
C ********************************************************************
C RECEIVE LOCAL PARAMETER FROM CONTROL CARD
11000 CONTINUE
C CALL OPTPAR(LPARAM,PARVAL)
LPARAM=3
GOTO (11010,11020,11030),LPARAM
11010 CUTOFF=PARVAL
RETURN
11020 MFITCO=LPARVA
RETURN
11030 GOTO 9981
9981 CALL OPTERR('MOD1SP')
STOP
END
C
SUBROUTINE CCACT
CCCACT CCACT FOR MDSCAL
C REAL*4 ATAB(4,22)
REAL ATAB(4,22)
C INTEGER*4 TAB(4,22),MTAB(4)
INTEGER TAB(4,22),MTAB(4)
INTEGER LPAR(25)
REAL PAR(25)
INTEGER IPARAM(5)
REAL C(73),WORD(3)
REAL CHTAB(13)
REAL BLANK,DOT
INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA,PARNO,TABNO
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
COMMON/OPTEX2/LPAR
COMMON/OPTEX3/MTAB,TAB
EQUIVALENCE(TAB,ATAB)
EQUIVALENCE(LPAR,PAR)
DATA BLANK,DOT,EQUALS,COMMA,C(73)/' ','.','=',',',' '/
DATA CHTAB/'0','1','2','3','4','5','6','7','8','9','.','+','-'/
DATA ASTER/'*'/
1 FORMAT(72A1)
11 FORMAT(' ',72A1)
2 FORMAT(1X,I1)
3 FORMAT(1X,18A4)
4 FORMAT(1X,I18)
5 FORMAT(1X,F18.3)
6 FORMAT(1X,72A1)
9 FORMAT(' MDSCAL DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.')
99 FORMAT(' ITEM NUMBER',I5,' HAS EXPECTED TYPE',I5,
+' AND ACTUAL TYPE',I5,' . THIS ITEM IS'',3A4,''')
C READ AND PRINT CONTROL CARD
100 READ(LREAD,1)(C(I),I=1,72)
WRITE(LPRINT,11)(C(I),I=1,72)
C ANALYZE CONTROL CARD INTO TOKENS
C EACH 'TOKEN' IS DELIMITED BY BLANKS
C THERE ARE THREE TYPES OF TOKENS.
C ALPHABETIC, INTEGER, DECIMAL
C ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT
C DECIMALS DISTINGUISHED BY DECIMAL POINT
REWIND LSCRAT
BLSW=1
NUMSW=1
DECSW=1
K=0
DO 300I=1,73
X=C(I)
GOTO (110,120),BLSW
110 IF(X.EQ.BLANK.OR.X.EQ.EQUALS.OR.X.EQ.COMMA)GOTO300
IF(X.EQ.ASTER)GOTO310
BLSW=2
JA=I
DO 115L=1,13
115 IF(X.EQ.CHTAB(L))NUMSW=2
IF(X.EQ.DOT)DECSW=2
GOTO 300
120 IF(X.EQ.BLANK.OR.X.EQ.EQUALS.OR.X.EQ.COMMA)GOTO130
IF(X.EQ.ASTER)GOTO130
IF(X.EQ.DOT)DECSW=2
GOTO 300
130 K=K+1
JB=MIN0(I-1,JA+16)
JC=18-(JB-JA+1)
TYPE=1
IF(NUMSW.EQ.2)TYPE=NUMSW+DECSW-1
WRITE(LSCRAT,2)TYPE
GOTO (140,150),NUMSW
140 WRITE(LSCRAT,6)(C(J),J=JA,JB),(BLANK,J=1,JC)
GOTO 160
150 WRITE(LSCRAT,6)(BLANK,J=1,JC),(C(J),J=JA,JB)
GOTO 160
160 BLSW=1
NUMSW=1
DECSW=1
IF(X.EQ.ASTER)GOTO310
GOTO 300
300 CONTINUE
C ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY
310 KB=K
IF(KB.EQ.0)RETURN
REWIND LSCRAT
XTYPE=1
IPARAM(1)=1
DO 1000J=1,KB
READ(LSCRAT,2)TYPE
IF(TYPE.NE.XTYPE)GOTO999
GOTO (400,410,420),XTYPE
400 READ(LSCRAT,3)(WORD(L),L=1,3)
GOTO 510
410 READ(LSCRAT,4)INTPAR
LPAR(PARNO)=INTPAR
GOTO 430
420 READ(LSCRAT,5)DECPAR
PAR(PARNO)=DECPAR
430 XTYPE=1
GOTO 1000
510 TABNO=IPARAM(1)
MA=MTAB(TABNO)
MB=MTAB(TABNO+1)-1
LMISSW=0
DO 700M=MA,MB
IF(WORD(1).NE.ATAB(1,M))GOTO700
LMISSW=1
600 XTYPE=1
IPARAM(1)=1
PARNO=TAB(3,M)
LTEMP=TAB(2,M)
GOTO (610,620,630,640,650),LTEMP
C NAME OF INTEGER PARAMETER
610 XTYPE=2
GOTO 700
C NAME OF DECIMAL PARAMETER
620 XTYPE=3
GOTO 700
C IMPLICITLY SPECIFIED INTEGER PARAMETER
630 LPAR(PARNO)=TAB(4,M)
GOTO 700
C IMPLICITLY SPECIFIED DECIMAL PARAMETER
640 PAR(PARNO)=ATAB(4,M)
GOTO 700
C INTERNAL PARAMETER OF CCACT PROGRAM
650 IPARAM(PARNO)=TAB(4,M)
GOTO 700
700 CONTINUE
IF(LMISSW.EQ.0)GOTO999
1000 CONTINUE
1001 RETURN
999 WRITE(LPRINT,9)
WRITE(LPRINT,99)K,XTYPE,TYPE,(WORD(L),L=1,3)
STOP
END
C
SUBROUTINE OPTEX
DIMENSION CC(12),TITLE(20)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
COMMON/OPTEX2/SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,LCNTRL,
+CUTOFF,MFITCO,LPUNSW,LPRINS,LCONSW
COMMON/TIL/TITLE
DATA BLANK/' '/
C INITIALIZE OPTEX
10 CONTINUE
WRITE(LPRINT,1)
1 FORMAT('1')
NOIT=50
STRMIN=0.01
SRATST=1.001
SRTAVW=0.33334
COSAVW=0.6666
ACSAVW=0.6666
LCONSW=1
DO 11K=1,20
11 TITLE(K)=BLANK
C PERFORM INITIAL ACTIONS
CALL CALC(9,0.0,0.0,0.0)
12 FORMAT(20A4)
13 FORMAT(1X,20A4)
100 LCNTRL=1
CALL CCACT
GOTO (100,200,300,400,500,600,700),LCNTRL
C READ DATA AND PERFORM INITIAL MANIPULATIONS
200 CALL CALC(3,0.0,0.0,0.0)
WRITE(LPRINT,2)
GOTO 100
C READ CONFIGURATION
300 CALL CALC(1,0.0,0.0,0.0)
LCONSW=2
WRITE(LPRINT,2)
GOTO 100
C READ PARAMETERS
400 CALL CALC(7,0.0,0.0,0.0)
GOTO 100
C SOLVE THE SYSTEM
500 CALL OPTSOL(SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,LCONSW,
+STRESS)
C CREATE OUTPUT
WRITE(LPRINT,50)STRESS
50 FORMAT(' FINAL CONFIGURATION HAS STRESS OF',2PF7.1,' PERCENT.')
CALL CALC(8,0.0,0.0,0.0)
CALL CALC(9,0.0,0.0,0.0)
GOTO 100
C PERFORM TERMINAL FUNCTIONS
600 CALL CALC(10,0.0,0.0,0.0)
C CONTROL MAY OR MAY NOT RETURN HERE.
RETURN
2 FORMAT(' ')
700 READ(LREAD,12)TITLE
WRITE(LPRINT,13)TITLE
GOTO 100
END
C
SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,IPOINT)
C FIT PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION
C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC
C AND WHICH BEST APPROXIMATE THE VALUES OF DIST,
C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS,
C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW ,
C IS A MINIMUM.
C IT DIRECTLY USES SORT.
DIMENSION IPOINT(501),DIST(500),DHAT(500),LBLOCK(501)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
C FORM FIRST APPROXIMATION TO CORRECT PARTITON
C IF LFITSW=0, DO PLAIN REGRESSION
C IF LFITSW GT 0, DO BLOCK REGRESSION
C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH
C BLOCK INDICATED BY LBLOCK ACCORDING TO DIST VALUES.
C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES.
C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START.
C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITH
C PARTITION INTO BLOCKS INDICATED BY LBLOCK.
C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUE
C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE
C THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THAT
C BLOCK. SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK
C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THE
C BLOCK.
C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME
C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE
C LAST DHAT VALUE IN THE BLOCK.
IF(LFITSW.GT.0)GOTO100
DO 90M=1,MM
DHAT(M)=DIST(M)
90 LBLOCK(M)=1
GOTO 500
100 MA=1
110 MAP=IPOINT(MA)
K=LBLOCK(MAP)
MB=MA+K-1
GOTO (200,300),LFITSW
C PRIMARY APPROACH
200 CONTINUE
IF(K-1)9999,220,210
C IF K=1, SAVE SORTING TIME
210 CALL SORT(IPOINT(MA),K,DIST(1),+1)
C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC
C ORDER.
220 DO 230M=MA,MB
MP=IPOINT(M)
DHAT(MP)=DIST(MP)
LBLOCK(MP)=1
230 CONTINUE
GOTO 400
C SECONDARY APPROACH
300 CONTINUE
TEMP1=0.0
DO 310M=MA,MB
MP=IPOINT(M)
TEMP1=TEMP1+DIST(MP)
310 CONTINUE
MAP=IPOINT(MA)
DHAT(MAP)=TEMP1
C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST
400 MA=MA+K
IF(MA-MM-1)110,500,9999
C START MAIN COMPUTATION *************************************
500 MA=1
510 LUD=2
MAP=IPOINT(MA)
NSATIS=0
520 K=LBLOCK(MAP)
MB=MA+K-1
MBP=IPOINT(MB)
WT=FLOAT(K)
550 DAV=DHAT(MAP)/WT
GOTO (600,700),LUD
C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE
600 CONTINUE
IF(MA-1)9999,630,610
610 MBD=MA-1
MBDP=IPOINT(MBD)
KD=LBLOCK(MBDP)
MAD=MBD-KD+1
MADP=IPOINT(MAD)
WTD=FLOAT(KD)
617 DAVD=DHAT(MADP)/WTD
IF(DAVD-DAV)630,620,620
620 KNEW=K+KD
LBLOCK(MADP)=KNEW
LBLOCK(MBP)=KNEW
DTONEW=DHAT(MADP)+DHAT(MAP)
DHAT(MADP)=DTONEW
NSATIS=0
MA=MAD
MAP=MADP
MAP=MADP
GOTO 800
630 NSATIS=NSATIS+1
GOTO 800
C IS BLOCK UP-SATISFIED. IF NOT, MERGE
700 CONTINUE
IF(MB-MM)710,730,9999
710 MAU=MB+1
MAUP=IPOINT(MAU)
KU=LBLOCK(MAUP)
MBU=MAU+KU-1
MBUP=IPOINT(MBU)
WTU=FLOAT(KU)
717 DAVU=DHAT(MAUP)/WTU
IF(DAV-DAVU)730,720,720
720 KNEW=K+KU
LBLOCK(MAP)=KNEW
LBLOCK(MBUP)=KNEW
DTONEW=DHAT(MAP)+DHAT(MAUP)
DHAT(MAP)=DTONEW
NSATIS=0
GOTO 800
730 NSATIS=NSATIS+1
GOTO 800
C PROCEED TO NEXT BLOCK IF READY.
800 LUD=3-LUD
C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETUR
IF(NSATIS-1)520,520,810
C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK.
810 CONTINUE
IF(MB-MM)820,900,9999
820 MA=MB+1
GOTO 510
C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT.
900 MA=1
910 MAP=IPOINT(MA)
K=LBLOCK(MAP)
MB=MA+K-1
IF(K-1)9999,940,920
920 TEMP1=DHAT(MAP)/FLOAT(K)
DO 930M=MA,MB
MP=IPOINT(M)
930 DHAT(MP)=TEMP1
GOTO 945
940 DHAT(MAP)=DIST(MAP)
945 MA=MB+1
LKK=7
IF(MA-MM-1)910,950,9999
950 RETURN
C TROUBLE EXIT
9999 WRITE(LPRINT,99)
99 FORMAT(
+'0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENT.')
RETURN
END
SUBROUTINE OPTPUN
DIMENSION TITLE(20)
COMMON/TIL/TITLE
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
WRITE(LPUNCH,19)TITLE
19 FORMAT(20A4/'CONFIGURATION')
RETURN
END
SUBROUTINE OPTPRI
DIMENSION TITLE(20)
COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
COMMON/TIL/TITLE
WRITE(LPRINT,13)TITLE
13 FORMAT(1X,20A4)
RETURN
END
C
SUBROUTINE PLOTR(X,Y,YP,NPOI,OUT,ID)
C
C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY -X- VS. -Y-.
C
C THE PARAMETERS -XMAX- -XMIN- YMAX- - YMIN- INDICATE THE UPPER
C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES
C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN.
C
C IT IS ASSUMED THAT THE ENTRIES HAVE -NPOI- ENTRIES
C
C THE PLOTTING IS DONE ON TAPE -OUT-.
C
C IF -ID- IS NEGATIVE, AXES WILL BE INCLUDED ON THE PLOT.
C IF -ID- IS POSITIVE, NO AXES WILL APPEAR.
C
C -LONG- IS THE LENGTH OF THE ARRAYS -X- AND -Y- IN THE CALLING PROGRAM.
C
C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965
C VERSION 3, APRIL 1967
C
C--------ADAPTED BY F. J. CARMONE,JR. 5/1/67
C
C
DIMENSION X(NPOI),Y(NPOI),ITEM(55,101),SMALL(21),HOLL(11),YP(
+NPOI)
REAL LABELS(55)
REAL ITEM
INTEGER OUT
DATA HOLL,DD,PLUS,BLANK/' ','X','2','3','4','5','6','7','8','9',
+'M','.','*',' '/,AIE,AMINUS/'I','-'/
DATA LABELS/15*' ','T','H','E',' ','X',' ','A','R','E',' ','T',
+'H','E',' ','L','I','N','E','A','R',' ','M','O','D','E','L',14*
+' '/
DATA AAB/'X'/,DOO/'0'/,PL/'+'/
3001 FORMAT(14X,
+'*.****.****.****.****.****.****.****.****.****.****.****.****.***
+*.****.****.****.****.****.****.****.*')
3300 FORMAT(' ',A1,F9.2,1X,105A1,F11.2)
3301 FORMAT(15X,A1,10(F9.4,A1))
3302 FORMAT(10X,11F10.4)
DO 115I=1,55
DO 115J=1,101
115 ITEM(I,J)=BLANK
116 CONTINUE
IF(ID.GT.0)GOTO119
DO 117I=1,55
117 ITEM(I,51)=AIE
DO 118I=1,101
118 ITEM(28,I)=AMINUS
119 CONTINUE
XMAX=X(1)
XMIN=X(2)
DO 121I=1,NPOI
IF(X(I).GT.XMAX)XMAX=X(I)
IF(X(I).LT.XMIN)XMIN=X(I)
121 CONTINUE
YMAX=Y(1)
YMIN=Y(2)
DO 123I=1,NPOI
IF(Y(I).GT.YMAX)YMAX=Y(I)
IF(YP(I).GT.YMAX)YMAX=YP(I)
IF(Y(I).LT.YMIN)YMIN=Y(I)
IF(YP(I).LT.YMIN)YMIN=YP(I)
123 CONTINUE
SPANX=XMAX-XMIN
SPANY=YMAX-YMIN
XMAX=XMAX+.05*SPANX
YMAX=YMAX+.05*SPANY
XMIN=XMIN-.05*SPANX
YMIN=YMIN-.05*SPANY
SPANX=XMAX-XMIN
SPANY=YMAX-YMIN
DELX=SPANX/100.0
DELY=SPANY/54.0
VALUE=YMAX+DELY
SMALL(1)=XMIN
DO 120I=1,20
120 SMALL(I+1)=SMALL(I)+5.0*DELX
DO 135II=1,NPOI
I=(YMAX-Y(II))/DELY+1.5
IP=(YMAX-YP(II))/DELY+1.5
J=(X(II)-XMIN)/DELX+1.5
IF(I.GT.56.OR.I.LT.1.OR.J.LT.1.OR.J.GT.101)GOTO135
IF(IP.GT.56.OR.IP.LT.1)GOTO135
IF(ITEM(I,J).NE.BLANK)GOTO200
ITEM(I,J)=DOO
GOTO 201
200 ITEM(I,J)=PL
201 IF(ITEM(IP,J).NE.BLANK)GOTO202
ITEM(IP,J)=AAB
GOTO 135
202 ITEM(IP,J)=PL
135 CONTINUE
WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)
WRITE(OUT,3302)(SMALL(I),I=1,21,2)
WRITE(OUT,3001)
DO 330I=1,55
VALUE=VALUE-DELY
A=BLANK
B=PLUS
L=I+2
IF(L/3*3-L)330,310,330
310 B=DD
IF(L/2*2-L)330,320,330
320 A=DD
330 WRITE(OUT,3300)LABELS(I),VALUE,A,B,(ITEM(I,J),J=1,101),B,A,VALUE
WRITE(OUT,3001)
WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)
WRITE(OUT,3302)(SMALL(I),I=1,21,2)
CONTINUE
RETURN
END
SUBROUTINE SORT(A,N,B,SWITCH)
C
C SORT ROUTINE MANIPULATES POINTERS IN 'A' WITHOUT MOVING
C VALUES IN 'B'
C THE N VALUES FROM A(1) TO A(N) ARE SORTED IN ALGEBRAIC ORDER.
C IF SWITCH IS +, ASCENDING ORDER.
C IF SWITCH IS -, DESCENDING ORDER.
C **** NOTE: THIS IS NOT THE SLOW SORT REFERRED TO IN THE WRITE UP BY
C KRUSKAL AND CARMONE,BUT A MODIFIED VERSION OF THE
C D. L. SHELL ALGORITHM.
C
DIMENSION A(500),B(501)
INTEGER SWITCH,A,TEMP
IF(N.LE.1)GOTO999
M=1
106 M=M+M
IF(M.LT.N)GOTO106
M=M-1
994 M=M/2
IF(M.EQ.0)GOTO999
KK=N-M
J=1
992 I=J
II=A(I)
996 IM=I+M
IMM=A(IM)
IF(SWITCH)810,810,800
800 IF(B(II).GT.B(IMM))GOTO110
GOTO 995
810 IF(B(II).LT.B(IMM))GOTO110
995 J=J+1
IF(J.GT.KK)GOTO994
GOTO 992
110 TEMP=A(I)
A(I)=A(IM)
A(IM)=TEMP
140 I=I-M
IF(I.LT.1)GOTO995
II=A(I)
GOTO 996
999 CONTINUE
RETURN
END
C
C