C monanova.f
C This software implements the method called MONANOVA.
C The author of this software is Joseph B Kruskal

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C For explanation of the method this software implements, see
C "Analysis of Factorial Experiments by Estimating Monotone
C Transformations of the Data" (1965) by Joseph B. Kruskal, in Journal
C of the Royal Statistical Society, Series B, 27(2), pp 251-263

C "MONANOVA, a Fortran-IV Program for Monotone Analysis of Variance
C (Non-Metric Analysis of Factorial Experiments)" (1969) by Joseph B.
C Kruskal and Frank J. Carmone in
C Behavioral Science 14, pp 165-166

C "MONANOVA, A Fortran IV Program for Monotone Analysis of Variance"
C (1969), Joseph B. Kruskal and Frank J. Carmone, in Journal of
C Marketing Research 6, p 497
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

C                                                                       MONA   1
C            M O N A N O V A   FOUR(4)  CHARACTER VERSION               MONA   2
C               NONMETRIC ANALYSIS OF FACTORIAL DESIGN                  MONA   3
C                                                                       MONA  A3
C                  DISTRIBUTED BY                                       MONA  B3
C             MARKETING SCIENCE INSTITUTE                               MONA  C3
C             PHILADELPHIA,PENNA.  19104                                MONA  D3
C                                                                       MONA   4
C        THIS PROGRAM WAS WRITTEN BY                                    MONA   5
C                                                                       MONA   6
C             DR. J. B. KRUSKAL                                         MONA   7
C             BELL TELEPHONE LABORATORIES                               MONA   8
C             MURRAY HILL, N. J.                                        MONA   9
C                                                                       MONA  10
C                                                                       MONA  11
C        ADAPTED FOR IBM S/360 BY                                       MONA  12
C                                                                       MONA  13
C             F. J. CARMONE, JR.                                        MONA  14
C             MARKETING SCIENCE INSTITUTE                               MONA  15
C             PHILADELPHIA, PA.                                         MONA  16
C          NEW ADDRESS
C              DEPARTMENT OF ECONOMICS
C              UNIV OF WATERLOO
C              WATERLOO, ONTARIO   CANADA
C                                                                       MONA  17
C        JULY  1968                                                     MONA  18
C                                                                       MONA  19
C                                                                       MONA  20
C        REFERENCE                                                      MONA  21
C        KRUSKAL, J. B. "ANALYSIS OF FACTORIAL EXPERIMENTS BY           MONA  22
C          ESTIMATING MONOTONE TRANSFORMATIONS OF THE DATA",            MONA  23
C          JOURNAL OF ROYAL STATISTICAL SOCIETY,SERIES B, VOL. 27,      MONA  24
C          NO. 2 (1965), PP. 251-63.                                    MONA  25
C                                                                       MONA  26
C                                                                       MONA  27
C                                                                       MONA  28
C      GENERAL REMARKS                                                  MONA  29
C                                                                       MONA  30
C        THIS PROGRAM CONTAINS THE FOLLOWING MAJOR ROUTINES             MONA  31
C                                                                       MONA  32
C        OPTEX   MAIN ROUTINE                                           MONA  33
C        OPSOL   CALCULATES STEP SIZE AND OTHER VARIABLES FOR TERMINATIOMONA  34
C        CALC    CALCULATES STRESS AND READS IN DATA                    MONA  35
C        MFIT    PERFORMS WEIGHTED LEAST SQUARES REGRESSION             MONA  36
C        CCACT   READS AND INTERPERTS CONTROL CARDS                     MONA  37
C        SORT    SORTS ARRAYS BY USING POINTERS; SIMPLE, BUT VERY FAST  MONA  38
C        PLOTR   PLOTS TWO MODELS VERSUS ORIGINAL DATA                  MONA  39
C                                                                       MONA  40
C                                                                       MONA  41
C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE.                   MONA  42
C                                                                       MONA  43
C ALL NORMAL INPUT-OUTPUT HANDLED BY   CALC   AND  CCACT.               MONA  44
C                                                                       MONA  45
C MFIT  HAS EMERGENCY DIAGNOSTIC OUTPUT.                                MONA  46
C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS CONTROLLED BY      MONA  47
C THESE VARIABLES.  LREAD, LPRINT, LPUNCH, LSCRAT.                      MONA  48
C UNIT NUMBERS 5, 6, 7, 2 HAVE BEEN USED RESPECTIVELY.                  MONA  49
C TO CHANGE THESE ASSIGNMENTS, IT SUFFICES TO CHANGE THE VALUES         MONA  50
C FOR THE FOUR VARIABLES SET IN THE BLOCK DATA SUBPROGRAM.              MONA  51
C                                                                       MONA  52
C NOTE THAT THE SCRATCH UNIT IS USED IN A VERY MINOR WAY.               MONA  53
C IT IS USED ONLY BY  CCACT                                             MONA  54
C MANY INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THINGMONA  55
C                                                                       MONA  56
       CALL OPTEX                                                       MONA  57
       CALL EXIT                                                        MONA  58
       END                                                              MONA  59
C                                                                       OPTS   1
       SUBROUTINE OPTSOL                                                OPTS   2
     1           ( SRATST, STRMIN, NOIT, SRTAVW, COSAVW, ACSAVW,        OPTS   3
     2 LCONSW, STRESS )                                                 OPTS   4
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    OPTS   5
C              INITIALIZE UNIVERSAL VARIABLES                           OPTS   6
  100  ITNO=0                                                           OPTS   7
       GRMULT=0.0                                                       OPTS   8
       GRMAG=1.0                                                        OPTS   9
       SRATAV=1.2                                                       OPTS  10
       COSAV=0.0                                                        OPTS  11
       ACSAV=0.2                                                        OPTS  12
       LXMIN=1                                                          OPTS  13
C              INITIALIZE CONFIGURATION AND GRADIENTS                   OPTS  14
  200  GO TO (210,220), LCONSW                                          OPTS  15
C              CREATE INITIAL CONFIGURATION                             OPTS  16
210    CALL CALC(  2,  0.0,  0.0,  0.0 )                                OPTS  17
C              NORMALIZE CONFIGURATION AND PRINT INPUT                  OPTS  18
220    CALL CALC(  6,  0.0,  0.0,  0.0 )                                OPTS  19
C              PERFORM ITERATION                                        OPTS  20
  300  STRLST=STRESS                                                    OPTS  21
       CALL CALC(  4,  STRESS, GRGR, GRGL )                             OPTS  22
305    GLMAG=GRMAG                                                      OPTS  23
       GRMAG = SQRT ( GRGR )                                            OPTS  24
       IF(STRESS) 9999, 320, 310                                        OPTS  25
  310  IF(GRGR ) 9999, 320, 330                                         OPTS  26
320    LXMIN=2                                                          OPTS  27
       GO TO 1000                                                       OPTS  28
  330  IF(ITNO) 9999, 340, 350                                          OPTS  29
340    SRAT=1.2                                                         OPTS  30
       CAGRGL=0.0                                                       OPTS  31
       GO TO 360                                                        OPTS  32
  350  SRAT=STRLST/STRESS                                               OPTS  33
       CAGRGL = GRGL / (GRMAG*GLMAG)                                    OPTS  34
C              GO TO UPDATING SECTION                                   OPTS  35
360    GO TO 1000                                                       OPTS  36
380    CALL CALC(  5,  GRMULT,  0.0,  0.0 )                             OPTS  37
       ITNO=ITNO+1                                                      OPTS  38
       GO TO 300                                                        OPTS  39
C              TERMINATION HAS OCCURRED. PERFORM TERMINAL ACTIONS.      OPTS  40
C              NORMALIZE CONFIGURATION.                                 OPTS  41
400    CALL CALC(  6,  0.0,  0.0,-50.0 )                                OPTS  42
C              AGAIN CALCULATE STRESS TO GUARANTEE AGAINST ROUNDING     OPTS  43
C              EFFECTS DURING NORMALIZATION                             OPTS  44
       CALL CALC(  4,  STRESS, GRGR, GRGL )                             OPTS  45
       RETURN                                                           OPTS  46
9999   CALL OPTERR(  6HOPTSOL  )                                        OPTS  47
      CALL EXIT                                                         OPTS  48
C      UPDATE VARIOUS VARIABLES. DECIDE WHETHER TO TERMINATE.           OPTS  49
1000   IF(ITNO) 9999, 1010, 1100                                        OPTS  50
1010   WRITE(LPRINT,1)                                                  OPTS  51
    1  FORMAT(     '0HISTORY OF COMPUTATION.'/1X)                       OPTS  52
       WRITE(LPRINT,2)                                                  OPTS  53
    2  FORMAT(    ' ITERATION STRESS    SRAT  SRTAVG CAGRGL  COSAV',    OPTS  54
     1 '  ACSAV    GRMAG   GRMULT     STEP'   )                         OPTS  55
       LFOOT=1                                                          OPTS  56
1100   GO TO (1110, 1300 ),LXMIN                                        OPTS  57
1110   COSAV = CAGRGL*COSAVW  +  COSAV*(1.0-COSAVW)                     OPTS  58
       ACSAV = ABS (CAGRGL)*ACSAVW  +  ACSAV*(1.0-ACSAVW)               OPTS  59
       SRATAV = (SRAT**SRTAVW) * (SRATAV**(1.0-SRTAVW))                 OPTS  60
       IF(ITNO) 9999, 1200, 1250                                        OPTS  61
1200   GRMULT=25.0*STRESS                                               OPTS  62
       STEP=GRMULT*GRMAG                                                OPTS  63
       GO TO 1300                                                       OPTS  64
1250   ANG=4.0**COSAV                                                   OPTS  65
      TEMP1 = 1.0 + (AMIN1 (1.0,1.0/SRATAV) ) ** 5                      OPTS  66
       TEMP3=ABS (COSAV)                                                OPTS  67
       TEMP2= 1.0 + ( (ACSAV-TEMP3) / (1.0 - TEMP3 ) ) ** 3             OPTS  68
       BIAS    = 1.6 / (TEMP1 * TEMP2 )                                 OPTS  69
      GOODLK=SQRT(AMIN1(1.0,SRAT))                                      OPTS  70
       STEP = STEP * ANG * BIAS * GOODLK                                OPTS  71
       GRMULT=STEP/GRMAG                                                OPTS  72
 1300  WRITE(LPRINT,3 )  ITNO,STRESS,SRAT,SRATAV,CAGRGL,                OPTS  73
     1COSAV,ACSAV,GRMAG,GRMULT,STEP                                     OPTS  74
    3  FORMAT(I10,F7.3,2F8.4,3F7.3,3F9.5)                               OPTS  75
 2000  IF(STRESS) 9999, 2010, 2020                                      OPTS  76
2010   LSTOP=1                                                          OPTS  77
       GO TO 3000                                                       OPTS  78
 2020  IF(GRMAG) 9999, 2030, 2040                                       OPTS  79
2030   LSTOP=2                                                          OPTS  80
       GO TO 3000                                                       OPTS  81
 2040  TEMP1 = 0.5 * (1.0+SRATST)                                       OPTS  82
       TEMP2 = ABS (TEMP1 - 1.0)                                        OPTS  83
       IF( ABS (SRAT-TEMP1) - TEMP2 ) 2050, 2050, 2070                  OPTS  84
 2050  IF( ABS (SRATAV-TEMP1) - TEMP2 ) 2060, 2060, 2070                OPTS  85
2060   LSTOP=3                                                          OPTS  86
       GO TO 3000                                                       OPTS  87
 2070  IF(STRESS-STRMIN) 2080, 2080, 2090                               OPTS  88
2080   LSTOP=4                                                          OPTS  89
       GO TO 3000                                                       OPTS  90
 2090  IF(NOIT-ITNO) 9999, 2100, 2110                                   OPTS  91
2100   LSTOP=5                                                          OPTS  92
       GO TO 3000                                                       OPTS  93
2110   GO TO (380, 400), LFOOT                                          OPTS  94
 3000  GO TO ( 3010, 3020 ), LFOOT                                      OPTS  95
 3010  LFOOT=2                                                          OPTS  96
       WRITE(LPRINT, 2)                                                 OPTS  97
       WRITE(LPRINT, 30)                                                OPTS  98
   30  FORMAT(1X)                                                       OPTS  99
 3020  GO TO (       3030, 3040, 3050, 3060, 3070), LSTOP               OPTS 100
 3030  WRITE(LPRINT, 31)                                                OPTS 101
   31  FORMAT(24H ZERO STRESS WAS REACHED )                             OPTS 102
       GO TO 2020                                                       OPTS 103
 3040  WRITE(LPRINT,32)                                                 OPTS 104
 32    FORMAT(21H0MINIMUM WAS ACHIEVED )                                OPTS 105
       GO TO 2040                                                       OPTS 106
 3050  WRITE(LPRINT,33)                                                 OPTS 107
   33  FORMAT(    ' APPROXIMATE MINIMUM WAS REACHED.')                  OPTS 108
       GO TO 2070                                                       OPTS 109
 3060  WRITE(LPRINT,34)                                                 OPTS 110
   34  FORMAT(32H SATISFACTORY STRESS WAS REACHED )                     OPTS 111
       GO TO 2090                                                       OPTS 112
 3070  WRITE(LPRINT,35)                                                 OPTS 113
   35  FORMAT(39H MAXIMUM NUMBER OF ITERATIONS WERE USED )              OPTS 114
       GO TO 2110                                                       OPTS 115
      END                                                               OPTS 116
C                                                                       CCAC   1
       SUBROUTINE CCACT                                                 CCAC   2
CCCACT             CCACT  FOR MDSCAL                                    CCAC   3
      REAL*4  ATAB(4,22)                                                CCAC   4
      INTEGER*4  TAB(4,22),MTAB(4)                                      CCAC   5
       EQUIVALENCE                         (TAB,ATAB)                   CCAC   6
       INTEGER                         LPAR(   25 )                     CCAC   7
       REAL                             PAR(   25 )                     CCAC   8
       EQUIVALENCE           (LPAR,PAR)                                 CCAC   9
       INTEGER  IPARAM(5)                                               CCAC  10
       REAL    C(73),WORD(3)                                            CCAC  11
       REAL CHTAB(13)                                                   CCAC  12
       REAL BLANK, DOT                                                  CCAC  13
       INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO           CCAC  14
       COMMON  /OPTEX1/  LREAD, LPRINT, LPUNCH, LSCRAT                  CCAC  15
       COMMON  /OPTEX2/  LPAR                                           CCAC  16
       COMMON  /OPTEX3/  MTAB,  TAB                                     CCAC  17
       DATA    BLANK, DOT, EQUALS, COMMA, C(73) /1H ,1H.,1H=,1H,,1H  /  CCAC  18
       DATA CHTAB/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H.,1H+,1H-/  CCAC  19
       DATA ASTER/1H*/                                                  CCAC  20
 1     FORMAT(72A1)                                                     CCAC  21
 11    FORMAT(1H0,72A1)                                                 CCAC  22
 2     FORMAT(1X,I1)                                                    CCAC  23
 3     FORMAT(1X,18A4)                                                  CCAC  24
 4     FORMAT(1X,I18)                                                   CCAC  25
 5     FORMAT(1X,F18.3)                                                 CCAC  26
 6     FORMAT(1X,72A1)                                                  CCAC  27
 9     FORMAT(51H MDSCAL DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.)   CCAC  28
  99    FORMAT( 12H ITEM NUMBER, I5, 20H   HAS EXPECTED TYPE ,  I5,     CCAC  29
     1  18H   AND ACTUAL TYPE,  I5, 16H . THIS ITEM IS', 3A4, 1H' )     CCAC  30
C      READ AND PRINT CONTROL CARD                                      CCAC  31
 100   READ  (LREAD,1)  (C(I),I=1,72)                                   CCAC  32
       WRITE(LPRINT,11) (C(I),I=1,72)                                   CCAC  33
C      ANALYZE CONTROL CARD INTO TOKENS                                 CCAC  34
C              EACH 'TOKEN' IS DELIMITED BY BLANKS                      CCAC  35
C              THERE ARE THREE TYPES OF TOKENS.                         CCAC  36
C              ALPHABETIC, INTEGER, DECIMAL                             CCAC  37
C              ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT  CCAC  38
C              DECIMALS       DISTINGUISHED BY DECIMAL POINT            CCAC  39
       REWIND LSCRAT                                                    CCAC  40
       BLSW=1                                                           CCAC  41
       NUMSW=1                                                          CCAC  42
       DECSW=1                                                          CCAC  43
       K=0                                                              CCAC  44
       DO 300 I=1,73                                                    CCAC  45
       X=C(I)                                                           CCAC  46
       GO TO (110,120), BLSW                                            CCAC  47
 110   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 300CCAC  48
      IF(X.EQ.ASTER) GO TO 310                                          CCAC  49
       BLSW=2                                                           CCAC  50
       JA=I                                                             CCAC  51
       DO 115 L=1,13                                                    CCAC  52
 115   IF(X.EQ.CHTAB(L)) NUMSW=2                                        CCAC  53
       IF(X.EQ.DOT) DECSW=2                                             CCAC  54
       GO TO 300                                                        CCAC  55
 120   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 130CCAC  56
      IF(X.EQ.ASTER) GO TO 130                                          CCAC  57
       IF(X.EQ.DOT) DECSW=2                                             CCAC  58
       GO TO 300                                                        CCAC  59
 130   K=K+1                                                            CCAC  60
       JB = MIN0 (I-1,JA+16)                                            CCAC  61
       JC=18-(JB-JA+1)                                                  CCAC  62
       TYPE=1                                                           CCAC  63
       IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1                                CCAC  64
       WRITE (LSCRAT,2) TYPE                                            CCAC  65
       GO TO (140,150),NUMSW                                            CCAC  66
 140   WRITE (LSCRAT,6) (C(J),J=JA,JB), (BLANK,J=1,JC)                  CCAC  67
       GO TO 160                                                        CCAC  68
 150   WRITE (LSCRAT,6) (BLANK,J=1,JC), (C(J),J=JA,JB)                  CCAC  69
       GO TO 160                                                        CCAC  70
 160   BLSW=1                                                           CCAC  71
       NUMSW=1                                                          CCAC  72
       DECSW=1                                                          CCAC  73
      IF(X.EQ.ASTER) GO TO 310                                          CCAC  74
       GO TO 300                                                        CCAC  75
 300   CONTINUE                                                         CCAC  76
C      ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY              CCAC  77
  310  KB=K                                                             CCAC  78
       IF(KB.EQ.0) RETURN                                               CCAC  79
       REWIND LSCRAT                                                    CCAC  80
       XTYPE=1                                                          CCAC  81
       IPARAM(1) = 1                                                    CCAC  82
       DO 1000 J=1,KB                                                   CCAC  83
       READ  (LSCRAT,2) TYPE                                            CCAC  84
       IF(TYPE.NE.XTYPE) GO TO 999                                      CCAC  85
       GO TO (400,410,420), XTYPE                                       CCAC  86
 400   READ  (LSCRAT,3) (WORD(L),L=1,3)                                 CCAC  87
       GO TO 510                                                        CCAC  88
 410   READ  (LSCRAT,4) INTPAR                                          CCAC  89
       LPAR(PARNO)=INTPAR                                               CCAC  90
       GO TO 430                                                        CCAC  91
 420   READ  (LSCRAT,5) DECPAR                                          CCAC  92
       PAR(PARNO)=DECPAR                                                CCAC  93
 430   XTYPE=1                                                          CCAC  94
       GO TO 1000                                                       CCAC  95
 510   TABNO = IPARAM(1)                                                CCAC  96
       MA = MTAB(TABNO)                                                 CCAC  97
       MB = MTAB(TABNO+1) - 1                                           CCAC  98
       LMISSW = 0                                                       CCAC  99
       DO 700 M=MA,MB                                                   CCAC 100
       IF(WORD(1).NE.ATAB(1,M)) GO TO 700                               CCAC 101
       LMISSW = 1                                                       CCAC 102
 600   XTYPE=1                                                          CCAC 103
       IPARAM(1) = 1                                                    CCAC 104
       PARNO=TAB(3,M)                                                   CCAC 105
      LTEMP = TAB(2,M)                                                  CCAC 106
       GO TO (610,620,630,640,650),  LTEMP                              CCAC 107
C      NAME OF INTEGER PARAMETER                                        CCAC 108
 610   XTYPE=2                                                          CCAC 109
       GO TO 700                                                        CCAC 110
C      NAME OF DECIMAL PARAMETER                                        CCAC 111
 620   XTYPE=3                                                          CCAC 112
       GO TO 700                                                        CCAC 113
C      IMPLICITLY SPECIFIED INTEGER PARAMETER                           CCAC 114
 630   LPAR(PARNO) = TAB(4,M)                                           CCAC 115
       GO TO 700                                                        CCAC 116
C      IMPLICITLY SPECIFIED DECIMAL PARAMETER                           CCAC 117
 640   PAR(PARNO) = ATAB(4,M)                                           CCAC 118
       GO TO 700                                                        CCAC 119
C      INTERNAL PARAMETER OF CCACT PROGRAM                              CCAC 120
 650   IPARAM(PARNO) = TAB(4,M)                                         CCAC 121
       GO TO 700                                                        CCAC 122
 700   CONTINUE                                                         CCAC 123
       IF(LMISSW.EQ.0)  GO TO 999                                       CCAC 124
 1000  CONTINUE                                                         CCAC 125
 1001  RETURN                                                           CCAC 126
 999   WRITE (LPRINT,9)                                                 CCAC 127
       WRITE  (LPRINT,99)  K,  XTYPE, TYPE, (WORD(L),L=1,3)             CCAC 128
      CALL EXIT                                                         CCAC 129
       END                                                              CCAC 130
       SUBROUTINE OPTERR( NAME )                                        OPTR   1
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    OPTR   2
      REAL*8 NAME                                                       OPTR   3
      WRITE(LPRINT,1) NAME                                              OPTR   4
1      FORMAT(     '0MONANOVA ERROR STOP FROM ', A6 )                   OPTR   5
      CALL EXIT                                                         OPTR   6
       END                                                              OPTR   7
C                                                                       OPTX   1
       SUBROUTINE OPTEX                                                 OPTX   2
       DIMENSION CC(12), TITLE(20)                                      OPTX   3
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    OPTX   4
      COMMON /OPTEX2/ SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW,          OPTX   5
     1 LCNTRL ,CUTOFF,MFITCO,LPUNSW,LPRINS,LCONSW                       OPTX   6
      COMMON /TIL/ TITLE                                                OPTX   7
       DATA  BLANK/1H /                                                 OPTX  A7
C              INITIALIZE OPTEX                                         OPTX   8
   10 CONTINUE                                                          OPTX   9
       WRITE(LPRINT, 1)                                                 OPTX  10
1      FORMAT(1H1)                                                      OPTX  11
       NOIT = 50                                                        OPTX  12
       STRMIN = 0.01                                                    OPTX  13
       SRATST=1.001                                                     OPTX  14
       SRTAVW=0.33334                                                   OPTX  15
       COSAVW=0.6666                                                    OPTX  16
       ACSAVW=0.6666                                                    OPTX  17
       LCONSW=1                                                         OPTX  18
       DO  11 K = 1,20                                                  OPTX A18
   11  TITLE(K) = BLANK                                                 OPTX B18
C              PERFORM INITIAL ACTIONS                                  OPTX  19
       CALL CALC( 9,  0.0,  0.0,  0.0 )                                 OPTX  20
12     FORMAT(20A4)                                                     OPTX  21
13     FORMAT(1X,20A4)                                                  OPTX  22
  100 LCNTRL=1                                                          OPTX  23
      CALL CCACT                                                        OPTX  24
       GO TO (100, 200, 300, 400, 500, 600,700), LCNTRL                 OPTX  25
C              READ DATA AND PERFORM INITIAL MANIPULATIONS              OPTX  26
200    CALL CALC(  3,  0.0,  0.0,  0.0 )                                OPTX  27
       WRITE(LPRINT,2)                                                  OPTX  28
       GO TO 100                                                        OPTX  29
C              READ CONFIGURATION                                       OPTX  30
300    CALL CALC(  1,  0.0,  0.0,  0.0 )                                OPTX  31
       LCONSW=2                                                         OPTX  32
       WRITE(LPRINT,2)                                                  OPTX  33
       GO TO 100                                                        OPTX  34
C              READ PARAMETERS                                          OPTX  35
400    CALL CALC(  7,  0.0,  0.0,  0.0 )                                OPTX  36
       GO TO 100                                                        OPTX  37
C              SOLVE THE SYSTEM                                         OPTX  38
500    CALL OPTSOL(SRATST, STRMIN, NOIT, SRTAVW, COSAVW, ACSAVW,        OPTX  39
     1 LCONSW, STRESS )                                                 OPTX  40
C              CREATE OUTPUT                                            OPTX  41
       WRITE(LPRINT,50)  STRESS                                         OPTX  42
50     FORMAT(    '0FINAL CONFIGURATION HAS STRESS OF', 2PF7.1,         OPTX  43
     1 ' PERCENT.'  )                                                   OPTX  44
       CALL CALC(  8,  0.0,  0.0,  0.0 )                                OPTX  45
       GO TO 10                                                         OPTX  46
C              PERFORM TERMINAL FUNCTIONS                               OPTX  47
600    CALL CALC( 10,  0.0,  0.0,  0.0 )                                OPTX  48
C      CONTROL MAY OR MAY NOT RETURN HERE.                              OPTX  49
       RETURN                                                           OPTX  50
2      FORMAT(1H    )                                                   OPTX  51
  700 READ(LREAD,12)  TITLE                                             OPTX  52
      WRITE(LPRINT,13)  TITLE                                           OPTX  53
      GO TO 100                                                         OPTX  54
       END                                                              OPTX  55
C                                                                       MFIT   1
      SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,IPOINT)                MFIT   2
C      FIT     PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION      MFIT   3
C      THIS ROUTINE FINDS THE VALUES FOR  DHAT  WHICH ARE MONOTONIC     MFIT   4
C      AND WHICH BEST APPROXIMATE  THE VALUES OF  DIST,                 MFIT   5
C      IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS,             MFIT   6
C      EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN  WW ,            MFIT   7
C      IS A MINIMUM.                                                    MFIT   8
C      IT DIRECTLY USES  SORT.                                          MFIT   9
      DIMENSION    IPOINT(1),   DIST(1),DHAT(1),LBLOCK(1)               MFIT  10
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    MFIT A10
C      FORM FIRST APPROXIMATION TO CORRECT PARTITON                     MFIT  11
C     IF LFITSW=0, DO PLAIN REGRESSION                                  MFIT  12
C     IF LFITSW GT 0, DO BLOCK REGRESSION                               MFIT  13
C              IF LFITSW=1,  USE PRIMARY APPROACH.  THUS WE SORT EACH   MFIT  14
C              BLOCK INDICATED BY LBLOCK ACCORDING TO DIST VALUES.      MFIT  15
C              BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES.  MFIT  16
C              THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START.  MFIT  17
C              IF LFITSW=2,   USE SECONDARY APPROACH. THUS WE START WITHMFIT  18
C              PARTITION INTO BLOCKS INDICATED BY LBLOCK.               MFIT  19
C              WITHIN EACH BLOCK  OF WHATEVER SIZE, THE FIRST DHAT VALUEMFIT  20
C      GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE MFIT  21
C              THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THATMFIT  22
C              BLOCK.  SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK   MFIT  23
C              VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEMFIT  24
C              BLOCK.                                                   MFIT  25
C      BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME MFIT  26
C              EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE  MFIT  27
C              LAST DHAT VALUE IN THE BLOCK.                            MFIT  28
      IF(LFITSW.GT.0)GOTO100                                            MFIT  29
      DO90M=1,MM                                                        MFIT  30
      DHAT(M)=DIST(M)                                                   MFIT  31
 90   LBLOCK(M)=1                                                       MFIT  32
      GOTO500                                                           MFIT  33
  100  MA=1                                                             MFIT  34
  110 MAP=IPOINT(MA)                                                    MFIT  35
      K=LBLOCK(MAP)                                                     MFIT  36
       MB=MA+K-1                                                        MFIT  37
       GO TO ( 200, 300 ), LFITSW                                       MFIT  38
C              PRIMARY APPROACH                                         MFIT  39
  200 CONTINUE                                                          MFIT  40
       IF(K-1) 9999, 220, 210                                           MFIT  41
C              IF K=1, SAVE SORTING TIME                                MFIT  42
  210 CALL SORT(IPOINT(MA),  K,DIST( 1),+1)                             MFIT  43
C              'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC   MFIT  44
C              ORDER.                                                   MFIT  45
  220  DO 230 M=MA,MB                                                   MFIT  46
      MP=IPOINT(M)                                                      MFIT  47
      DHAT(MP)=DIST(MP)                                                 MFIT  48
      LBLOCK(MP)=1                                                      MFIT  49
 230  CONTINUE                                                          MFIT  50
       GO TO 400                                                        MFIT  51
C              SECONDARY APPROACH                                       MFIT  52
 300  CONTINUE                                                          MFIT  53
       TEMP1=0.0                                                        MFIT  54
       DO 310 M=MA,MB                                                   MFIT  55
      MP=IPOINT(M)                                                      MFIT  56
      TEMP1=TEMP1+DIST(MP)                                              MFIT  57
 310   CONTINUE                                                         MFIT  58
      MAP=IPOINT(MA)                                                    MFIT  59
      DHAT(MAP)=TEMP1                                                   MFIT  60
C              PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST              MFIT  61
  400  MA=MA+K                                                          MFIT  62
       IF(MA-MM-1) 110, 500, 9999                                       MFIT  63
C      START MAIN COMPUTATION      *************************************MFIT  64
  500  MA=1                                                             MFIT  65
  510  LUD=2                                                            MFIT  66
      MAP=IPOINT(MA)                                                    MFIT  67
       NSATIS=0                                                         MFIT  68
  520  K=LBLOCK(MAP)                                                    MFIT  69
       MB=MA+K-1                                                        MFIT  70
      MBP=IPOINT(MB)                                                    MFIT  71
      WT=FLOAT(K)                                                       MFIT  72
 550   DAV = DHAT(MAP)/ WT                                              MFIT  73
       GO TO ( 600, 700 ), LUD                                          MFIT  74
C              IS BLOCK DOWN-SATISFIED. IF NOT, MERGE                   MFIT  75
  600 CONTINUE                                                          MFIT  76
       IF(MA-1) 9999, 630, 610                                          MFIT  77
  610  MBD=MA-1                                                         MFIT  78
      MBDP=IPOINT(MBD)                                                  MFIT  79
       KD=LBLOCK(MBDP)                                                  MFIT  80
       MAD=MBD-KD+1                                                     MFIT  81
      MADP=IPOINT(MAD)                                                  MFIT  82
      WTD=FLOAT(KD)                                                     MFIT  83
 617   DAVD = DHAT(MADP)/ WTD                                           MFIT  84
       IF( DAVD-DAV ) 630, 620, 620                                     MFIT  85
  620  KNEW=K+KD                                                        MFIT  86
       LBLOCK(MADP)=KNEW                                                MFIT  87
       LBLOCK(MBP)=KNEW                                                 MFIT  88
       DTONEW = DHAT(MADP ) + DHAT(MAP)                                 MFIT  89
       DHAT(MADP)= DTONEW                                               MFIT  90
       NSATIS=0                                                         MFIT  91
       MA=MAD                                                           MFIT  92
      MAP=MADP                                                          MFIT  93
      MAP=MADP                                                          MFIT  94
       GO TO 800                                                        MFIT  95
  630  NSATIS=NSATIS+1                                                  MFIT  96
       GO TO 800                                                        MFIT  97
C              IS BLOCK UP-SATISFIED. IF NOT, MERGE                     MFIT  98
  700 CONTINUE                                                          MFIT  99
       IF(MB-MM) 710, 730, 9999                                         MFIT 100
  710  MAU=MB+1                                                         MFIT 101
      MAUP=IPOINT(MAU)                                                  MFIT 102
      KU=LBLOCK(MAUP)                                                   MFIT 103
       MBU=MAU+KU-1                                                     MFIT 104
      MBUP=IPOINT(MBU)                                                  MFIT 105
      WTU=FLOAT(KU)                                                     MFIT 106
 717   DAVU = DHAT(MAUP)/ WTU                                           MFIT 107
       IF( DAV-DAVU ) 730, 720, 720                                     MFIT 108
  720  KNEW=K+KU                                                        MFIT 109
       LBLOCK(MAP)=KNEW                                                 MFIT 110
       LBLOCK(MBUP)=KNEW                                                MFIT 111
       DTONEW = DHAT(MAP)+ DHAT(MAUP)                                   MFIT 112
       DHAT(MAP)= DTONEW                                                MFIT 113
       NSATIS=0                                                         MFIT 114
       GO TO 800                                                        MFIT 115
  730  NSATIS=NSATIS+1                                                  MFIT 116
       GO TO 800                                                        MFIT 117
C              PROCEED TO NEXT BLOCK IF READY.                          MFIT 118
  800  LUD = 3-LUD                                                      MFIT 119
C              QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURMFIT 120
       IF(NSATIS-1) 520, 520, 810                                       MFIT 121
C              QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK.  MFIT 122
  810 CONTINUE                                                          MFIT 123
       IF(MB-MM) 820, 900, 9999                                         MFIT 124
  820  MA=MB+1                                                          MFIT 125
       GO TO 510                                                        MFIT 126
C      MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT.                MFIT 127
  900  MA=1                                                             MFIT 128
  910 MAP  =IPOINT(MA)                                                  MFIT 129
      K=LBLOCK(MAP)                                                     MFIT 130
       MB=MA+K-1                                                        MFIT 131
       IF(K-1) 9999, 940, 920                                           MFIT 132
 920  TEMP1=DHAT(MAP)/FLOAT(K)                                          MFIT 133
       DO 930 M=MA,MB                                                   MFIT 134
      MP=IPOINT(M)                                                      MFIT 135
  930  DHAT(MP)=TEMP1                                                   MFIT 136
       GO TO 945                                                        MFIT 137
 940   DHAT(MAP) = DIST(MAP)                                            MFIT 138
 945   MA = MB + 1                                                      MFIT 139
      LKK=7                                                             MFIT 140
       IF(MA-MM-1) 910, 950, 9999                                       MFIT 141
  950  RETURN                                                           MFIT 142
C      TROUBLE EXIT                                                     MFIT 143
9999   WRITE(LPRINT,99)                                                 MFIT 144
 99   FORMAT(58H0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENMFIT 145
     1T.)                                                               MFIT 146
      RETURN                                                            MFIT 147
      END                                                               MFIT 148
       SUBROUTINE  OPTPUN                                               OPUN   1
      DIMENSION  TITLE(20)                                              OPUN   2
      COMMON /TIL/ TITLE                                                OPUN   3
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    OPUN   4
      WRITE(LPUNCH,19)  TITLE                                           OPUN   5
   19 FORMAT(20A4/'CONFIGURATION')                                      OPUN   6
      RETURN                                                            OPUN   7
      END                                                               OPUN   8
       SUBROUTINE  OPTPRI                                               OPRI   1
      DIMENSION  TITLE(20)                                              OPRI   2
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    OPRI   3
      COMMON /TIL/ TITLE                                                OPRI   4
      WRITE(LPRINT,13)  TITLE                                           OPRI   5
   13 FORMAT(1X,20A4)                                                   OPRI   6
      RETURN                                                            OPRI   7
      END                                                               OPRI   8
       SUBROUTINE  SORT(A,N,B,SWITCH)                                   SORT   1
C                                                                       SORT   2
C       SORT ROUTINE MANIPULATES POINTERS IN "A" WITHOUT MOVING         SORT   3
C        VALUES IN "B"                                                  SORT   4
C      THE N VALUES FROM A(1) TO A(N) ARE SORTED IN ALGEBRAIC ORDER.    SORT   5
C      IF SWITCH IS +, ASCENDING ORDER.                                 SORT   6
C      IF SWITCH IS -, DESCENDING ORDER.                                SORT   7
C   **** NOTE   THIS IS NOT THE SLOW SORT REFERRED TO IN THE WRITE UP BYSORT   8
C                  KRUSKAL AND CARMONE,BUT A MODIFIED VERSION OF THE    SORT   9
C                  D. L. SHELL ALGORITHM.                               SORT  10
C                                                                       SORT  11
       DIMENSION  A(1),B(1)                                             SORT  12
       INTEGER  SWITCH,A,TEMP                                           SORT  13
       IF(N.LE.1) GO TO 999                                             SORT  14
      M=1                                                               SORT  15
  106 M=M+M                                                             SORT  16
      IF(M.LE.N) GO TO 106                                              SORT  17
      M=M-1                                                             SORT  18
 994       M = M/2                                                      SORT  19
       IF(M.EQ.0) GO TO 999                                             SORT  20
       KK = N-M                                                         SORT  21
       J = 1                                                            SORT  22
992     I = J                                                           SORT  23
      II=A(I)                                                           SORT  24
996     IM = I + M                                                      SORT  25
      IMM=A(IM)                                                         SORT  26
       IF(SWITCH)  810,810,800                                          SORT  27
800     IF(B(II).GT.B(IMM)) GO TO 110                                   SORT  28
       GO TO 995                                                        SORT  29
810    IF(B(II).LT.B(IMM))  GO TO 110                                   SORT  30
995     J = J+1                                                         SORT  31
       IF(J.GT.KK)  GO TO 994                                           SORT  32
       GO TO 992                                                        SORT  33
 110   TEMP=A(I )                                                       SORT  34
       A(I )=A(IM )                                                     SORT  35
       A(IM )=TEMP                                                      SORT  36
140       I = I-M                                                       SORT  37
       IF(I.LT.1)  GO TO 995                                            SORT  38
      II=A(I)                                                           SORT  39
        GO TO 996                                                       SORT  40
 999   CONTINUE                                                         SORT  41
       RETURN                                                           SORT  42
      END                                                               SORT  43
C                                                                       BDAT   1
       BLOCK   DATA                                                     BDAT   2
       INTEGER                         LPAR(   25 )                     BDAT   3
       REAL                             PAR(   25 )                     BDAT   4
       EQUIVALENCE           (LPAR,PAR)                                 BDAT   5
      REAL*4  ATAB(4,22)                                                BDAT   6
      INTEGER*4  TAB(4,22),MTAB(4)                                      BDAT   7
      EQUIVALENCE  (TAB,ATAB)                                           BDAT   8
       COMMON  /OPTEX1/  LREAD, LPRINT, LPUNCH, LSCRAT                  BDAT   9
       COMMON  /OPTEX2/  LPAR                                           BDAT  10
       COMMON  /OPTEX3/  MTAB,  TAB                                     BDAT  11
       DATA    LREAD, LPRINT, LPUNCH, LSCRAT     /5, 6,7,2/             BDAT  12
      DATA  MTAB /1,18,20,22/                                           BDAT  13
      DATA  TAB/4HSTRM,2,2,0,'SRAT',2,1,0,'ITER',1,3,0,'COSA',2,5,0,    BDAT  14
     1'ACSA',2,6,0, 'SRTA',2,4,0,'DATA',3,7,2,'CONF',3,7,3,'PARA',3,7,4,BDAT  15
     2 'COMP',3,7,5,'STOP',3,7,6, 'TITL',3,7,7, 'CUTO',2,8,0,'TIES',    BDAT  16
     35,1,2,'CARD',3,10,1,'NOCA',3,10,2,'PRIN',5,1,3,                   BDAT  17
     4    'PRIM',3,9,1,'SECO',3,9,2,'NOIN',3,11,1,'INPU',3,11,2,0,0,0,0/BDAT  18
      END                                                               BDAT  19
C                                                                       CALC   1
       SUBROUTINE CALC( LLL, AAA, BBB, CCC )                            CALC   2
C      FOR         CALC            NONMETRIC ANALYSIS OF FACTORIAL DESIGCALC   3
C      DIMENSION DATA(501), INDEX(500, 4), X(500), XHAT(500)            CALC   4
       DIMENSION DATA(2001), INDEX(2000,5), X(2000), XHAT(2000)
C      DIMENSION IPOINT(501),LTIES(501)                                 CALC   5
       DIMENSION IPOINT(2001), LTIES(2001)
C      DIMENSION CON(4,100), GR(4,100), GL(4,100), GR2(4,100)           CALC   6
       DIMENSION  CON(20,10), GR(20,10), GL(20,10), GR2(20,10)
C      DIMENSION II(4), IT(4)                                           CALC   7
       DIMENSION II(20),  IT(20)
      DIMENSION  LPAR(25),FMTCON(20),FMTDAT(20)                         CALC   8
       COMMON  /OPTEX1/   LREAD,LPRINT,LPUNCH,LSCRAT                    CALC   9
      COMMON  /OPTEX2/  LPAR                                            CALC  10
       EQUIVALENCE(PARVAL,LPARVA)                                       CALC  11
      EQUIVALENCE (CUTOFF,LPAR(8)),(MFITCO,LPAR(9)),(LPUNSW,LPAR(10))   CALC  12
      EQUIVALENCE  (LPRINS,LPAR(11)),(LCONSW,LPAR(12))                  CALC  13
       GO TO ( 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000,    CALC  14
     1 1, 11000 ), LLL                                              CALC  150000
C   ********************************************************************CALC  16
C              READ CONFIGURATION                                       CALC  17
1000   READ(LREAD,100) KK,(II(K),K=1,KK)                                CALC  18
       WRITE(LPRINT,100) KK,(II(K),K=1,KK)                              CALC  19
  100 FORMAT(24I3)                                                      CALC  20
      READ(LREAD,101)  FMTCON                                           CALC  21
      WRITE(LPRINT,105)  FMTCON                                         CALC  22
  101 FORMAT(20A4)                                                      CALC  23
  105  FORMAT(2X,20A4)                                                  CALC  24
       DO 1010 K=1,KK                                                   CALC  25
       IIK=II(K)                                                        CALC  26
1010   READ(LREAD,FMTCON)(CON(K,I),I=1,IIK)                             CALC  27
       RETURN                                                           CALC  28
C   ********************************************************************CALC  29
C              CREATE INITIAL CONFIGURATION                             CALC  30
2000  DO 2005 K=1,KK                                                    CALC  31
      IIK=II(K)                                                         CALC  32
      DO 2005 I=1,IIK                                                   CALC  33
      CON(K,I)=0                                                        CALC  34
      GR (K,I)=0                                                        CALC  35
       GL(K,I)=0                                                        CALC  36
      GR2(K,I)=0                                                        CALC  37
 2005 CONTINUE                                                          CALC  38
      DO 2010 L=1,LDATA                                                 CALC  39
       M=IPOINT(L)                                                      CALC  40
       TEMP=DATA(M)                                                     CALC  41
      DO 2010 K=1,KK                                                    CALC  42
       I=INDEX(M,K)                                                     CALC  43
      CON(K,I)=CON(K,I)+TEMP                                            CALC  44
 2010 GR(K,I)=GR(K,I)+1.0                                               CALC  45
      DO 2020 K=1,KK                                                    CALC  46
      IIK=II(K)                                                         CALC  47
      DO 2020 I=1,IIK                                                   CALC  48
      CON(K,I)=CON(K,I)/GR(K,I)                                         CALC  49
      GR (K,I)=0                                                        CALC  50
 2020 CONTINUE                                                          CALC  51
      RETURN                                                            CALC  52
C   ********************************************************************CALC  53
C              READ DATA                                                CALC  54
 3000  READ(LREAD,100) KK,(II(K),K=1,KK),LREPLI                         CALC  55
      WRITE(LPRINT,100) KK,(II(K),K=1,KK),LREPLI                        CALC  56
      READ(LREAD,101)  FMTDAT                                           CALC  57
      WRITE(LPRINT,105)  FMTDAT                                         CALC  58
       IIS=0                                                            CALC  59
       IIP=1                                                            CALC  60
       DO 3010 K=1,KK                                                   CALC  61
       IIS=IIS+II(K)                                                    CALC  62
       IIP=IIP*II(K)                                                    CALC  63
 3010  IT(K)=1                                                          CALC  64
       NOD=IIP*LREPLI                                                   CALC  65
 3020  READ(LREAD,FMTDAT)  (DATA(L), L=1,NOD )                          CALC  66
      L=0                                                               CALC  67
      LR=0                                                              CALC  68
       DATMAX=-1.0E+20                                                  CALC  69
       DATMIN= -DATMAX                                                  CALC  70
       DO 3090 LA=1,NOD                                                 CALC  71
       TDATA= DATA(LA)                                                  CALC  72
       IF( TDATA - CUTOFF )  3045, 3045, 3030                           CALC  73
 3030  L=L+1                                                            CALC  74
       DATMAX =AMAX1 ( DATMAX, TDATA )                                  CALC  75
       DATMIN =AMIN1  (DATMIN, TDATA )                                  CALC  76
       DATA(L) = TDATA                                                  CALC  77
       DO 3040 K=1,KK                                                   CALC  78
 3040 INDEX(L,K)=IT(K)                                                  CALC  79
3045   LR=LR+1                                                          CALC  80
       IF( LR-LREPLI ) 3090, 3050, 9993                                 CALC  81
 3050 LR=0                                                              CALC  82
       IT(KK)=IT(KK)+1                                                  CALC  83
       DO 3080 KC=1,KK                                                  CALC  84
       K=KK+1-KC                                                        CALC  85
       IF( IT(K)-II(K)-1 ) 3090, 3060, 9993                             CALC  86
 3060  IT(K)=1                                                          CALC  87
       IF(K-1) 9993, 3100, 3070                                         CALC  88
 3070  IT(K-1)=IT(K-1)+1                                                CALC  89
 3080  CONTINUE                                                         CALC  90
 3090  CONTINUE                                                         CALC  91
 3100  LDATA=L                                                          CALC  92
       LDATAP = LDATA + 1                                               CALC  93
       DATA(LDATAP) = 1.0E+20                                           CALC  94
       DO 3110 L=1,LDATAP                                               CALC  95
 3110  IPOINT(L) = L                                                    CALC  96
      CALL SORT(IPOINT(1),LDATA,DATA(1),+1)                             CALC  97
       RETURN                                                           CALC  98
9993   CALL OPTERR(  6HMOD1RD    )                                      CALC  99
      CALL EXIT                                                         CALC 100
C   ********************************************************************CALC 101
C              CALCULATE STRESS AND GRADIENT                            CALC 102
4000  XBAR=0.0                                                          CALC 103
      DO 4010 L=1,LDATA                                                 CALC 104
      TEMP=0.0                                                          CALC 105
       M=IPOINT(L)                                                      CALC 106
       DO 4005 K=1,KK                                                   CALC 107
       I=INDEX(M,K)                                                     CALC 108
4005  TEMP=TEMP+CON(K,I)                                                CALC 109
      X(M)=TEMP                                                         CALC 110
 4010 XBAR=XBAR+TEMP                                                    CALC 111
      XBAR=XBAR/FLOAT (LDATA)                                           CALC 112
      K=1                                                               CALC 113
      L1=1                                                              CALC 114
      M=IPOINT(1)                                                       CALC 115
       TEMP2=DATA(M)                                                    CALC 116
      DO   50 L=1,LDATA                                                 CALC 117
      TEMP1=TEMP2                                                       CALC 118
      M=IPOINT(L+1)                                                     CALC 119
       TEMP2=DATA(M)                                                    CALC 120
      IF(TEMP1-TEMP2)   30,   20,   30                                  CALC 121
   20 K=K+1                                                             CALC 122
      GO TO   50                                                        CALC 123
   30 L2=L1+K-1                                                         CALC 124
      DO   40 LA=L1,L2                                                  CALC 125
      M=IPOINT(LA)                                                      CALC 126
      LTIES(M)=K                                                        CALC 127
   40 CONTINUE                                                          CALC 128
      L1=L2+1                                                           CALC 129
       K=1                                                              CALC 130
   50 CONTINUE                                                          CALC 131
      CALL MFIT(X(1),XHAT(1),LTIES(1),LDATA,MFITCO,IPOINT)              CALC 132
      U=0                                                               CALC 133
      T=0                                                               CALC 134
      DO 4030 K=1,KK                                                    CALC 135
      IIK=II(K)                                                         CALC 136
      DO 4030 I=1,IIK                                                   CALC 137
       GR(K,I)=0                                                        CALC 138
 4030 GR2(K,I)=0                                                        CALC 139
      DO 4050 L=1,LDATA                                                 CALC 140
       M=IPOINT(L)                                                      CALC 141
      TEMP1=X(M)-XHAT(M)                                                CALC 142
      TEMP2=X(M)-XBAR                                                   CALC 143
      U=U+TEMP1**2                                                      CALC 144
      T=T+TEMP2**2                                                      CALC 145
      DO 4040 K=1,KK                                                    CALC 146
       I=INDEX(M,K)                                                     CALC 147
      GR(K,I)=GR(K,I)+TEMP1                                             CALC 148
 4040 GR2(K,I)=GR2(K,I)+TEMP2                                           CALC 149
 4050 CONTINUE                                                          CALC 150
      U=SQRT (U)                                                        CALC 151
      T=SQRT (T)                                                        CALC 152
      STRESS=U/T                                                        CALC 153
      GRGR=0                                                            CALC 154
      GRGL=0                                                            CALC 155
       IF(U) 9994, 4080, 4060                                           CALC 156
 4060 COEF1=1.0/(U*T)                                                   CALC 157
      COEF2=-U/T**3                                                     CALC 158
      DO 4070 K=1,KK                                                    CALC 159
      IIK=II(K)                                                         CALC 160
      DO 4070 I=1,IIK                                                   CALC 161
      GR(K,I)=-COEF1*GR(K,I) - COEF2*GR2(K,I)                           CALC 162
      GRGR=GRGR+GR(K,I)**2                                              CALC 163
 4070 GRGL=GRGL+GR(K,I)*GL(K,I)                                         CALC 164
 4080  CONTINUE                                                         CALC 165
       AAA=STRESS                                                       CALC 166
       BBB=GRGR                                                         CALC 167
       CCC=GRGL                                                         CALC 168
       RETURN                                                           CALC 169
9994   CALL OPTERR(  6HMOD1CA  )                                        CALC 170
       RETURN                                                           CALC 171
C   ********************************************************************CALC 172
C              MOVE TO NEW CONFIGURATION                                CALC 173
5000   GRMULT=AAA                                                       CALC 174
      DO 5010 K=1,KK                                                    CALC 175
      IIK=II(K)                                                         CALC 176
      DO 5010 I=1,IIK                                                   CALC 177
      CON(K,I)=CON(K,I)+GRMULT*GR(K,I)                                  CALC 178
      GL(K,I)=GR(K,I)                                                   CALC 179
      GR (K,I)=0                                                        CALC 180
 5010 GR2(K,I)=0                                                        CALC 181
      RETURN                                                            CALC 182
C   ********************************************************************CALC 183
C              NORMALIZE CONFIGURATION                                  CALC 184
6000   SMULT=0                                                          CALC 185
      IF(LPRINS.NE.2)  GO TO 6040                                       CALC 186
      IF(CCC.LT.0.0) GO TO 6040                                         CALC 187
      WRITE(LPRINT,602)                                                 CALC 188
  602 FORMAT(1H1,' SEQ. NO.   DATA     SUBSCRIPTS'/)                    CALC 189
      DO 6050 L=1,LDATA                                                 CALC 190
      WRITE(LPRINT,601)  L,DATA(L),(INDEX(L,K),K=1,KK)                  CALC 191
  601 FORMAT(I6,4X,F10.5,20I5)
 6050 CONTINUE                                                          CALC 193
      IF(LCONSW.EQ.1) GO TO 6040                                        CALC 194
      WRITE(LPRINT,603)                                                 CALC 195
  603 FORMAT(///5X,' CONFIGURATION'/)                                   CALC 196
       DO 8011 K=1,KK                                                   CALC 197
       IIK=II(K)                                                        CALC 198
       WRITE(LPRINT,8005)  IIK,(CON(K,I),I=1,IIK)                       CALC 199
 8011 CONTINUE                                                          CALC 200
 6040  DO 6020 K=1,KK                                                   CALC 201
      AV=0                                                              CALC 202
       IIK=II(K)                                                        CALC 203
       DO 6005 I=1,IIK                                                  CALC 204
       GL(K,I)=0                                                        CALC 205
6005   AV=AV+CON(K,I)                                                   CALC 206
       AV=AV/FLOAT (IIK)                                                CALC 207
       DO 6010 I=1,IIK                                                  CALC 208
       CON(K,I)=CON(K,I)-AV                                             CALC 209
 6010  SMULT=SMULT+CON(K,I)**2                                          CALC 210
 6020  CONTINUE                                                         CALC 211
       SMULT=SQRT ( FLOAT (IIS)/SMULT)                                  CALC 212
       DO 6030 K=1,KK                                                   CALC 213
       IIK=II(K)                                                        CALC 214
       DO 6030 I=1,IIK                                                  CALC 215
 6030  CON(K,I)=SMULT*CON(K,I)                                          CALC 216
       RETURN                                                           CALC 217
C   ********************************************************************CALC 218
C              READ PARAMETERS                                          CALC 219
7000   READ(LREAD,103) CUTOFF,MFITCO                                    CALC 220
  103 FORMAT(2F10.4)                                                    CALC 221
       RETURN                                                           CALC 222
C   ********************************************************************CALC 223
C              PRINT, PUNCH, AND DRAW OUTPUT                            CALC 224
8000   CALL OPTPRI                                                      CALC 225
       DO 8010 K=1,KK                                                   CALC 226
       IIK=II(K)                                                        CALC 227
       WRITE(LPRINT,  8005)  IIK,(CON(K,I),I=1,IIK)                     CALC 228
 8005  FORMAT(1X,I3,10F7.3/  (4X,10F7.3) )                              CALC 229
 8010  CONTINUE                                                         CALC 230
      WRITE(LPRINT,8001)                                                CALC 231
 8001 FORMAT(1H1)                                                       CALC 232
      CALL PLOTR(DATA,XHAT,X,LDATA,6,+1)                                CALC 233
      IF(LPUNSW.EQ.2) RETURN                                            CALC 234
       CALL OPTPUN                                                      CALC 235
       WRITE(LPUNCH, 106 )                                              CALC 236
  106 FORMAT(2X, ' (2X,10F7.3) ')                                       CALC 237
       WRITE(LPUNCH,8020) KK,(II(K),K=1,KK)                             CALC 238
 8020  FORMAT(20I3)                                                     CALC 239
       DO 8040 K=1,KK                                                   CALC 240
       IIK=II(K)                                                        CALC 241
       WRITE(LPUNCH,8030) (CON(K,I),I=1,IIK)                            CALC 242
 8030  FORMAT(2X,10F7.3)                                                CALC 243
 8040  CONTINUE                                                         CALC 244
       RETURN                                                           CALC 245
C   ********************************************************************CALC 246
C              INITIAL ACTIONS                                          CALC 247
9000   CUTOFF=  -1.23E+20                                               CALC 248
       MFITCO=1                                                         CALC 249
      LPUNSW=2                                                          CALC 250
       RETURN                                                           CALC 251
C   ********************************************************************CALC 252
1 CONTINUE                                                          CALC 2530000
      RETURN                                                            CALC 254
C   ********************************************************************CALC 255
C              RECEIVE LOCAL PARAMETER FROM CONTROL CARD                CALC 256
11000 CONTINUE                                                          CALC 257
C      CALL OPTPAR(LPARAM,PARVAL)                                       CALC 258
       LPARAM = 3                                                       CALCA259
       GO TO (11010, 11020, 11030), LPARAM                              CALC 259
11010  CUTOFF=PARVAL                                                    CALC 260
       RETURN                                                           CALC 261
11020  MFITCO=LPARVA                                                    CALC 262
       RETURN                                                           CALC 263
11030  GO TO 9981                                                       CALC 264
9981   CALL OPTERR(  6HMOD1SP )                                         CALC 265
      CALL EXIT                                                         CALC 266
       END                                                              CALC 267
C                                                                       PLOT   1
      SUBROUTINE PLOTR(X,Y,YP,NPOI,OUT,ID)                              PLOT   2
C                                                                       PLOT   3
C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY -X- VS. -Y-.             PLOT   4
C                                                                       PLOT   5
C THE PARAMETERS -XMAX- -XMIN- YMAX- - YMIN- INDICATE THE UPPER         PLOT   6
C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES    PLOT   7
C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN.            PLOT   8
C                                                                       PLOT   9
C IT IS ASSUMED THAT THE ENTRIES HAVE -NPOI- ENTRIES                    PLOT  10
C                                                                       PLOT  11
C THE PLOTTING IS DONE ON TAPE -OUT-.                                   PLOT  12
C                                                                       PLOT  13
C IF -ID- IS NEGATIVE, AXES WILL BE INCLUDED ON THE PLOT.               PLOT  14
C IF -ID- IS POSITIVE, NO AXES WILL APPEAR.                             PLOT  15
C                                                                       PLOT  16
C -LONG- IS THE LENGTH OF THE ARRAYS -X- AND -Y- IN THE CALLING PROGRAM.PLOT  17
C                                                                       PLOT  18
C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965                           PLOT  19
C      VERSION 3, APRIL 1967                                            PLOT  20
C                                                                       PLOT  21
C--------ADAPTED BY F. J. CARMONE,JR.  5/1/67                           PLOT  22
C                                                                       PLOT  23
C                                                                       PLOT  24
      DIMENSION X(NPOI),Y(NPOI),ITEM(55,101),SMALL(21),HOLL(11),YP(NPOI)PLOT  25
      REAL LABELS(55)                                                   PLOT  26
      REAL ITEM                                                         PLOT  27
      INTEGER OUT                                                       PLOT  28
      DATA HOLL,DD,PLUS,BLANK/1H ,1HX,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,  PLOT  29
     11HM,1H.,1H*,1H /,AIE,AMINUS/1HI,1H-/                              PLOT  30
      DATA LABELS/15*1H ,'T','H','E',' ','X',' ','A','R','E',' ',       PLOT  31
     1'T','H','E',' ','L','I','N','E','A','R',' ','M','O','D','E','L',  PLOT  32
     1 14*1H /                                                          PLOT  33
      DATA AAB/1HX/,DOO/'0'/,PL/'+'/                                    PLOT  34
 3001 FORMAT(14X,103H*.****.****.****.****.****.****.****.****.****.****PLOT  35
     1.****.****.****.****.****.****.****.****.****.****.*)             PLOT  36
 3300 FORMAT('  ',A1,F9.2,1X,105A1,F11.2)                               PLOT  37
 3301 FORMAT(15X,A1,10(F9.4,A1))                                        PLOT  38
 3302 FORMAT(10X,11F10.4)                                               PLOT  39
      DO 115 I=1,55                                                     PLOT  40
      DO 115 J=1,101                                                    PLOT  41
  115 ITEM(I,J)=BLANK                                                   PLOT  42
  116 CONTINUE                                                          PLOT  43
      IF(ID.GT.0)GO TO 119                                              PLOT  44
      DO 117 I=1,55                                                     PLOT  45
  117 ITEM(I,51)=AIE                                                    PLOT  46
      DO 118 I=1,101                                                    PLOT  47
  118 ITEM(28,I)=AMINUS                                                 PLOT  48
  119 CONTINUE                                                          PLOT  49
      XMAX=X(1)                                                         PLOT  50
      XMIN =X(2)                                                        PLOT  51
      DO 121 I=1,NPOI                                                   PLOT  52
      IF(X(I).GT.XMAX)XMAX=X(I)                                         PLOT  53
      IF(X(I).LT.XMIN)XMIN=X(I)                                         PLOT  54
  121 CONTINUE                                                          PLOT  55
      YMAX=Y(1)                                                         PLOT  56
      YMIN=Y(2)                                                         PLOT  57
      DO 123 I=1,NPOI                                                   PLOT  58
      IF(Y(I).GT.YMAX)YMAX=Y(I)                                         PLOT  59
      IF(YP(I).GT.YMAX)YMAX=YP(I)                                       PLOT  60
      IF(Y(I).LT.YMIN)YMIN=Y(I)                                         PLOT  61
      IF(YP(I).LT.YMIN)YMIN=YP(I)                                       PLOT  62
  123 CONTINUE                                                          PLOT  63
      SPANX=XMAX-XMIN                                                   PLOT  64
      SPANY=YMAX-YMIN                                                   PLOT  65
      XMAX=XMAX+.05*SPANX                                               PLOT  66
      YMAX=YMAX+.05*SPANY                                               PLOT  67
      XMIN=XMIN-.05*SPANX                                               PLOT A67
      YMIN=YMIN-.05*SPANY                                               PLOT B67
      SPANX=XMAX-XMIN                                                   PLOT  68
      SPANY=YMAX-YMIN                                                   PLOT  69
      DELX=SPANX/100.0                                                  PLOT  70
      DELY=SPANY/54.0                                                   PLOT  71
      VALUE=YMAX+DELY                                                   PLOT  72
      SMALL(1)=XMIN                                                     PLOT  73
      DO 120 I=1,20                                                     PLOT  74
  120 SMALL(I+1)=SMALL(I)+5.0*DELX                                      PLOT  75
      DO 135 II=1,NPOI                                                  PLOT  76
      I=(YMAX-Y(II))/DELY+1.5                                           PLOT  77
      IP=(YMAX-YP(II))/DELY+1.5                                         PLOT  78
      J=(X(II)-XMIN)/DELX+1.5                                           PLOT  79
      IF(I.GT.56.OR.I.LT.1.OR.J.LT.1.OR.J.GT.101)GO TO 135              PLOT  80
      IF(IP.GT.56.OR.IP.LT.1)GO TO 135                                  PLOT  81
  201 IF(ITEM(IP,J).NE.BLANK) GO TO 202                                 PLOT  86
      ITEM(IP,J)=AAB                                                    PLOT  87
       GO TO 203
  202 ITEM(IP,J)=PL                                                     PLOT  89
  203  CONTINUE
      ITEM( I,J)=DOO                                                    PLOT  83
  135 CONTINUE                                                          PLOT  90
      WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)                          PLOT  91
      WRITE(OUT,3302)(SMALL(I),I=1,21,2)                                PLOT  92
      WRITE(OUT,3001)                                                   PLOT  93
      DO 330 I=1,55                                                     PLOT  94
      VALUE=VALUE-DELY                                                  PLOT  95
      A=BLANK                                                           PLOT  96
      B=PLUS                                                            PLOT  97
      L=I+2                                                             PLOT  98
      IF(L/3*3-L)330,310,330                                            PLOT  99
  310 B=DD                                                              PLOT 100
      IF(L/2*2-L)330,320,330                                            PLOT 101
  320 A=DD                                                              PLOT 102
  330 WRITE(OUT,3300)LABELS(I),VALUE,A,B,(ITEM(I,J),J=1,101),B,A,VALUE  PLOT 103
      WRITE(OUT,3001)                                                   PLOT 104
      WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)                          PLOT 105
      WRITE(OUT,3302)(SMALL(I),I=1,21,2)                                PLOT 106
      CONTINUE                                                          PLOT 107
      RETURN                                                            PLOT 108
      END                                                               PLOT 109