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