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