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  ***************************************************************
C     OPTIONS /NOF77
C
C
C
      BLOCK DATA
      COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
      COMMON/OPTEX2/SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,LCNTRL,
     +CUTOFF,MFITCO,LPUNSW,LPRINS,LCONSW
     
      COMMON/OPTEX3/MTAB,TAB
      INTEGER LREAD,LPRINT,LPUNCH,LSCRAT
      DATA LREAD,LPRINT,LPUNCH,LSCRAT/5,6,7,8/
      INTEGER MTAB(4)/1,18,20,22/
      INTEGER*4 TAB(4,22)/'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/
       END

      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
      INTEGER LREAD,LPRINT,LPUNCH,LSCRAT
      COMMON/OPTEX2/LPAR,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,LCNTRL,
     +CUTOFF,MFITCO,LPUNSW,LPRINS,LCONSW

      COMMON/OPTEX3/MTAB,TAB
      EQUIVALENCE(LPAR,PAR)
      EQUIVALENCE(TAB,ATAB)
      OPEN (5,FILE='MONANOVA.DAT', STATUS='old')
      OPEN (6,FILE='monanova.out',STATUS='UNKNOWN')
      OPEN (7,FILE='MONANOVA.XHA',STATUS='UNKNOWN')
      OPEN (8,FILE='OUTT8',STATUS='UNKNOWN')
      CALL OPTEX
      STOP
      END
C     
      SUBROUTINE OPTSOL(SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,
     +LCONSW,STRESS)
      COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
      INTEGER 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
      INTEGER 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)
      CHARACTER*1 FMTCON(7),FMTDAT(7)
      COMMON/OPTEX1/LREAD,LPRINT,LPUNCH,LSCRAT
      INTEGER 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
C     PRINT *,LLL
      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,801)(FMTCON(I),I=1,7)
  801 FORMAT(7A1)
C     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,801)FMTDAT
C     PRINT *,"FMTDAT",FMTDAT
C     PRINT *,"THIS IS THE END"
      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
      INTEGER 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,END=1001)(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,END=1001)TYPE
      IF(TYPE.NE.XTYPE)GOTO999
      GOTO (400,410,420),XTYPE
 400  READ(LSCRAT,3,END=1001)(WORD(L),L=1,3)
      GOTO 510
 410  READ(LSCRAT,4,END=1001)INTPAR
      LPAR(PARNO)=INTPAR
      GOTO 430
 420  READ(LSCRAT,5,END=1001)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/' '/
      INTEGER LREAD,LPRINT,LPUNCH,LSCRAT
C     LREAD=5
C     LPRINT=6
C     LPUNCH=7
C     LSCRAT=8
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
      INTEGER 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
      INTEGER 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
      INTEGER 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