C-PC--------------------------------------------------------------------
C     A L S C A L  PC  Version                     Last change: Febr. 90
C-----------------------------------------------------------------------

C     A L S C A L 8 4                                 VERSION 84.1
C     ===============                                 JUNE 1983

C COPYRIGHT 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 01JUN83

C     ALSCAL PERFORMS METRIC OR NONMETRIC:
C       MULTIDIMENSIONAL SCALING (SIMPLE, REPLICATED, WEIGHTED,
C         GENERALIZED, ASYMMETRIC, OR WEIGHTED ASYMMETRIC MODELS).
C       MULTIDIMENSIONAL UNFOLDING (SIMPLE, REPLICATED OR WEIGHTED)

C     ALSCAL ANALYZES DATA WHICH ARE EITHER
C       TWO OR THREE WAY
C       SYMMETRIC, ASYMMETRIC OR RECTANGULAR
C       ROW OR MATRIX CONDITIONAL, OR UNCONDITIONAL
C       REPLICATED OR UNREPLICATED
C       BINARY, NOMINAL, ORDINAL, INTERVAL OR RATIO
C       DISCRETE OR CONTINUOUS
C       WITH OR WITHOUT MISSING DATA (ANY PATTERN)

C-----------------------------------------------------------------------

C     THE ALTERNATING LEAST SQUARES ALGORITHM IS DESCRIBED IN:

C     1) TAKANE, Y., YOUNG, F.W., AND DELEEUW, J. NONMETRIC INDIVIDUAL
C        DIFFERENCES MULTIDIMENSIONAL SCALING. PSYCHOMETRICA,
C        VOL. 42, NO. 1, MARCH 1977, PP. 7-67.
C     2) YOUNG, F.W., TAKANE, Y. AND LEWYCKYJ, R. THREE NOTES ON ALSCAL
C        PSYCHOMETRIKA, VOL. 43, 1978, P 433-435.
C     3) YOUNG, F.W. AN ASYMMETRIC EUCLIDEAN MODEL FOR MULTIPROCESS
C        ASYMMETRIC DATA.  US-JAPAN MDS SEMINAR PROCEEDINGS, 1975
C     4) YOUNG, F.W. ENHANCEMENTDS IN ALSCAL-82.  SUGI PROCEEDINGS,
C        1982, 7, 633-642.
C     5) YOUNG, F.W. THE GENERAL EUCLIDEAN MODEL. IN: LAW, H.G., ET.AL.
C        (EDS.), THREE-MODE MODELS FOR DATA ANALYSIS. PRAEGER 1984



C THIS IS THE MAIN ROUTINE FOR THE FIXED CORE VERSION OF ALSCAL 84.1

C     THIS ROUTINE CALLS STEP0 WHICH CALCULATES THE NUMBER OF WORDS
C     OF STORAGE REQUIRED BY THE PROBLEM.  IT THEN CHECKS TO SEE IF
C     THE PROGRAM HAS SUFFICIENT MEMORY ALLOCATED.  IF NOT THE PROGRAM
C     TERMINATES WITH A MESSAGE INDICATING THE AMOUNT OF MEMORY NEEDED.
C     IF SUFFICIENT MEMORY IS ALLOCATED THIS ROUTINE CALLS DRIVER
C     WHICH PROCEEDS WITH THE ANALYSIS OF THE DATA.

C-----------------------------------------------------------------------

C PROGRAM LIMITS

C     A) THE NUMBER OF DIMENSIONS OF THE MULTIDIMENSIONAL SOLUTION MAY
C        NEVER EXCEED SIX.
C     B) ALL NOMINAL OR ORDINAL OBSERVATIONS MUST BE NOT GREATER THAN
C        9.0E20.
C     C) THE MAXIMUM TOTAL PROBLEM SIZE MAY NOT EXCEED THE VALUE SET
C        FOR VARIABLE MAXSIZ BELOW.

C-----------------------------------------------------------------------

C I/O DEVICE ASSIGNMENTS

C     ALL I/O DEVICES ARE SPECIFIED IN THE MAIN ROUTINE.
C     THE PROGRAM USES 3 SCRATCH FILES IN ADDITION TO THE 3 STANDARD
C     UNITS.   THE STANDARD UNITS ARE ASSUMED TO BE NUMBERED 1, 2, AND
C     3, AND THE SCRATCH FILES 11, 12, 13.  THESE SPECIFICATIONS
C     MAY BE CHANGED IN THE MAIN ROUTINE.

C-----------------------------------------------------------------------

C VARIABLES USED IN THIS ROUTINE

C     AREA     THIS ARRAY CONTAINS THE WORK AREA FOR ALL DIMENSIONED
C                  VARIABLES IN THE REMAINING ROUTINES.  THE LENGTH
C                  OF THIS ARRAY DETERMINES THE MAXIMUM PROBLEM SIZE.
C     MAXSIZ   THE MAXIMUM PROBLEM SIZE IN WORDS
C     PRSIZE   THE SIZE OF THE PROBLEM BEING ANALYZED IN WORDS
C     EOJ      END OF JOB INDICATOR
C     IONUMS   COMMON AREA WITH INPUT, OUTPUT AND SCRATCH UNIT NUMBERS

C-----------------------------------------------------------------------

C INSTRUCTIONS FOR SETTING MAXIMUM PROBLEM SIZE

C     A) CHOOSE AN ESTIMATED MAXIMUM PROBLEM SIZE, IN WORDS.
C        CONSULTATION WITH YOUR COMPUTER CENTER, PLUS KNOWLEDGE OF
C        THE NATURE OF YOUR DATA, PLUS CALCULATING THE VALUE OF THE
C        APPROXIMATE PROBLEM SIZE FORMULA GIVEN IN THE USERS MANUAL
C        WILL HELP IN DETERMING THIS VALUE.
C        NOTE THAT THIS VALUE DETERMINE THE AMOUNT OF CORE, IN WORDS,
C        NEEDED BY ALSCAL.
C     B) SET THE LENGTH OF ARRAY AREA EQUAL TO YOUR SELECTED MAXIMUM
C        PROBLEM SIZE.
C     C) SET VARIABLE MAXSIZ EQUAL TO THE VALUE YOU HAVE SELECTED.
C     D) RECOMPILE THE MAIN ROUTINE, AND RELINK ALSCAL.
C     E) ANALYZE YOUR DATA.  IF THE MAXIMUM PROBLEM SIZE IS TOO SMALL,
C        SEE BELOW.

C-----------------------------------------------------------------------

C INSTRUCTIONS FOR CHANGING MAXIMUM PROBLEM SIZE

C     A) DETERMINE MAXIMUM PROBLEM SIZE IN WORDS.  THIS IS EQUAL TO THE
C        ESTIMATED PROBLEM SIZE PRINTED BY ALSCAL WHEN THE MAXIMUM
C        PROBLEM SIZE IS TOO SMALL TO PERMIT A PROBLEM TO RUN.
C     B) CHANGE THE LENGTH OF ARRAY AREA TO BE AT LEAST AS LONG AS THE
C        ESTIMATED PROBLEM SIZE.  TO BE ON THE SAFE SIDE, ADD
C        AT LEAST 10 PERCENT TO THE SIZE ESTIMATED BY ALSCAL.
C     C) CHANGE THE VALUE OF MAXSIZ TO EQUAL THE LENGTH OF AREA.
C     D) RECOMPILE MAIN AND RELINK ALSCAL.
C     E) YOU ARE NOW READY TO TRY TO REANNALYZE YOUR DATA.



      SUBROUTINE MDS(NROW_MY,NCOL_MY,XX_MY,NDTYP_MY,NDIM_MY,NDMN_MY,
     *X1_MY,X2_MY,X3_MY,X4_MY,X5_MY,X6_MY,STRSS1)
C     ==============

      INTEGER PRSIZE,EOJ,AREA
      CHARACTER COPR*72,BLANK*1
C???      DIMENSION AREA(38000),XX_MY(NROW_MY,NCOL_MY)
      DIMENSION AREA(10000000),XX_MY(NROW_MY,NCOL_MY)
C-----FOLLOWING LINE ADDED 12/1/88 FOR COPYRIGHT COMMON-BLOCK
      COMMON /CCOPR/COPR
     +       /IONUMS/NDP,NDQ,NDR,NDPP
      COPR='COPYRIGHT 1977 BY  F. W. YOUNG, Y. TAKANE & R. J. LEWYCKYJ'
      DATA BLANK/' '/
C???      DATA MAXSIZ/38000/
      DATA MAXSIZ/10000000/

C----------------------------------------------------------------------

C      I/O DEVICE SPECIFICATION
C      ------------------------

C  IN=1 IS THE READER
C  NDP=11 IS A SCRATCH UNIT
C  NDQ=12 IS A SCRATCH UNIT
C  NDPP=13 IS A SCRATCH UNIT

      IN=1
      NDP=11
      NDQ=12
      NDPP=13

C-----------------------------------------------------------------------

      OPEN(NDP,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN(NDQ,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN(NDPP,STATUS='SCRATCH',FORM='UNFORMATTED')

C----------------------------------------------------------------------

    1 EOJ=-1
      CALL STEP0(PRSIZE,EOJ,NROW_MY,NCOL_MY,NDTYP_MY,NDIM_MY,NDMN_MY)
C???      CALL STEP0(PRSIZE,EOJ)
      IF(EOJ.NE.0) GOTO 99
      IF(PRSIZE.GT.MAXSIZ) GOTO 999
      DO 10 I=1,PRSIZE
   10 AREA(I)=0
C???      WRITE(*,20) MAXSIZ,PRSIZE
C???   20 FORMAT(/' MAXIMUM PROBLEM SIZE :',I8,' WORDS',
C???     +/       ' REQUIRED PROBLEM SIZE:',I8)
      CALL DRIVER(AREA,EOJ,NROW_MY,NCOL_MY,XX_MY,X1_MY,X2_MY,X3_MY,
     *X4_MY,X5_MY,X6_MY,STRSS1)
C???      CALL DRIVER(AREA,EOJ)
C???      IF(EOJ.EQ.0) GOTO 1
      GOTO 99
  999 WRITE(*,9001)PRSIZE,PRSIZE
 9001 FORMAT(/' ALSCAL FATAL ERROR:  NOT ENOUGH MEMORY AVAILABLE.'/
     1' *****RECOMPILE PROGRAM SO THAT THE MAIN ROUTINE HAS ITS DIMENSIO
     2N AND DATA STATEMENTS CHANGED TO READ:'/
     3/'          DIMENSION AREA(',I7,')'/
     4 '          DATA MAXSIZ/',I7,'/'/
     5/' *****RUN ABORTED.')


   99 CLOSE(NDP)
      CLOSE(NDQ)
      CLOSE(NDPP)

      RETURN
      END


C***********************************************************************

      SUBROUTINE STEP0(NWORDS,EOJ,NROW_MY,NCOL_MY,NDTYP_MY,NDIM_MY,
     *NDMN_MY)
C???      SUBROUTINE STEP0(NWORDS,EOJ)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 01JUN83

C  THIS ROUTINE READS AND CHECKS THE PROBLEM PARAMETERS.
C  IT THEN CALCULATES THE STORAGE EXTENT AND STRUCTURE FOR THE PROBLEM.
C      THE COMMON /STARTS/ HOLDS THE STARTING POSITIONS OF VARIOUS
C  ARRAYS USED BY THE PROBLEM WITHIN A LARGE BLOCK OF STORAGE
C  ALLOCATED BY THE MAIN ROUTINE. THE NAMING CONVENTION USED IS
C  TO PREFIX THE ARRAY NAME BY THE LETTER J. IE. JXX GIVES THE
C  STARTING POSITION OF THE ARRAY XX  ETC.

C----------------------------------------------------------------------

      CHARACTER*80,FMT,COPR*72
      CHARACTER*8,CCC*40
C-PC  CHARACTER*8,FIX,DYN,TYPE,D*1,F*1,VER*1
      COMMON /CCOPR/COPR
C-PC---------------------
     +       /CPPLOT/NROW
C------------------------

C-----FOLLOWING LINE ADDED 8/5/82 FOR GEMSCAL
     +       /PDCOM/IFLPDS,NDIR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
     +       /STARTS/JX,JWA,JWD,JXX,JDS,JCFR,JCFL,JW,JWS,JZZ,JTR,JFK,
     + JCV,JCW,JU11,JU12,JU22,JUB1,JUB2,
     + JXN,JPHSUB,JPHSTI,JIX,JIY,JIZ,JNDSR,JNAD,
     + JFTLN,JDIST,JPIJP,JXK
C-----JFTLN AND JDIST ADDED TO THE ABOVE STATEMENT 7/8/82
C-----JPIJP AND JXK ADDED TO THE ABOVE STATEMENT FOR GEMSCAL 8/5/82
     +      /INICON/INITX,INITW,INITWS,INITXC
      INTEGER EOJ
C-----THE FOLLOWING THREE STATEMENTS ADDED 7/7/82
      DOUBLE PRECISION CUT,STMIN,STMIND
      INTEGER DEBUG,ICNSTR,NOULB
      COMMON /PRMBLK/CUT,STMIN,DEBUG,ICNSTR,NOULB

C----------------------------------------------------------------------

C      FIXED CONSTANTS
C      ---------------

C  MAXDIM=6 DEFINES THE MAXIMUM NUMBER OF DIMENSIONS.
C      THIS VALUE MAY NOT BE CHANGED.
C  BIG=9.0E20 DEFINES THE LARGEST NOMINAL OR ORDINAL DATUM ALLOWED.
C      THIS MAY BE CHANGED.
C  EPSID=.001 DEFINES THE DEFAULT ITERATION CONVERGENCE CRITERION.
C      THIS MAY BE CHANGED.
C  MAXITD=30 DEFINES THE DEFAULT MAXIMUM NUMBER OF ITERATIONS.
C      THIS MAY BE CHANGED.
C  STMIND=.005 DEFINES THE DEFAULT MINIMUN S-STRESS.
C      THIS MAY BE CHANGED.

C----------------------------------------------------------------------

      PARAMETER(MAXDIM=6,EPSID=.001,MAXITD=30,STMIND=.005)
C-PC  DATA D,F,DYN,FIX/'D','F','DYNAMIC ','  FIXED '/
      BIG=9.0E20

C-PC  VER=D
C-PC  TYPE=DYN
      IF(EOJ.EQ.0)GOTO 2
C-PC  VER=F
C-PC  TYPE=FIX
      EOJ=0
    2 NADCT=0

C-------------------------------------------------


C     DATA SPECIFICATIONS
C     -------------------

C     NROW     NUMBER OF ROW STIMULI (NO MAXIMUM)
C     NCOL     NUMBER OF COLUMN STIMULI (NO MAXIMUM)
C     NS       NUMBER OF MATRICES (NO MAXIMUM)
C     NDTYP    MEASUREMENT LEVEL
C                   =1   RATIO (POLYNOMIAL W/O INTERCEPT)
C                   =2   INTERVAL (POLYNOMIAL WITH INTERCEPT)
C                   =3   ORDINAL(DEFAULT)
C                   =4   NOMINAL
C     NSIM     DATA FORM
C                   =0   SYMMETRIC-DISSIMILARITY(DEFAULT)
C                   =1   SYMMETRIC-SIMILARITY
C                   =2   ASYMMETRIC-DISSIMILARTY
C                   =3   ASYMMETRIC-SIMILARITY
C                   =4   RECTANGULAR-DISSIMILARITY
C                   =5   RECTANGULAR-SIMILARITY
C     NPS      MEASUREMENT PROCESS (ONLY WHEN NDTYP=3)
C                   =1   DISCRETE(DEFAULT)
C                   =2   CONTINUOUS
C     NWC      MEASUREMENT CONDITIONALITY
C                   =1   MATRIX CONDITIONAL(DEFAULT)
C                   =2   ROW CONDITIONAL
C                   =3   UNCONDITIONAL
C     NDEG     DEGREES IN POLYNOMIAL (WHEN NDTYP=1 OR 2)
C     NDMX     MAXIMUM NUMBER OF ORDINAL OBSERVATION CATEGORIES.
C                   (DEFAULT= TOTAL NUMBER OF CELLS, OR 1000,
C                             WHICHEVER IS SMALLER)
C     CUT      CUTOFF FOR MISSING DATA (DEFAULT 0.0)


C---CUT ADDED TO THE FOLLOWING STATEMENT 7/7/82
C???      READ(IN,101,ERR=905,END=905)
C     +NROW,NCOL,NS,NDTYP,NSIM,NPS,NWC,NDEG,NDMX,CUT
C  101 FORMAT(9I4,F8.4)

      NROW=NROW_MY
      NCOL=NCOL_MY
      NS=1
      NDTYP=NDTYP_MY
      NSIM=0
      NPS=1
      NWC=1
      NDEG=0
      NDMX=100000
      CUT=0

      IF(NROW.LT.3)GOTO 903
      IF(NS.LE.0)GOTO 904


      NB=NROW
      IF(NSIM.GT.3)NB=NB+NCOL
      NBS=NB**2
      NBNBNS=NB**2*NS
      NC=NB*(NB-1)
      NC2=NC/2
      NDX=MIN0(MAXDIM,NB-2)

      NCST=1

      IF(NDTYP.EQ.1) THEN
	NCST=0
      ELSE
         IF(NDTYP.NE.2.AND.NDTYP.NE.4) NDTYP=3
      ENDIF

      IF(NSIM.LT.0.OR.NSIM.GT.5) NSIM=0
      IF(NSIM.GT.1) NC2=NBS
      NT=NC2*NS

      IF(NPS.NE.2) NPS=1
      IF(NPS.EQ.2.AND.NDTYP.EQ.4) THEN
	WRITE(*,9910)
	WRITE(*,228)
	NPS=1
      ENDIF
  228 FORMAT(8X,'THE PROGRAM WILL CONTINUE WITH THE DISCRETE'
     1,' MEASUREMENT PROCESS (NPS=1)')

      IF(NWC.NE.2) THEN
         IF(NWC.EQ.3) THEN
            NWC=0
         ELSE
            NWC=1
         ENDIF
      ENDIF


      IF(NDEG.LT.1) NDEG=1
      IF(NDTYP.LE.2) THEN
	IF(NDEG.LE.1) GOTO 813
	IF(NDEG.NE.2.AND.NDEG.NE.3) NDEG=4
      ENDIF
  813 IF(NROW.NE.NCOL.AND.NSIM.LT.4)GOTO 902


C     ANALYSIS SPECIFICATIONS
C     -----------------------

C     NWE      MODEL TYPE
C                   =0   SIMPLE EUCLIDIAN MODEL (DEFAULT)
C                   =1   INDIVIDUAL DIFFERENCES MODEL
C                   =2   MULTIPROCESS ASYMMETRIC MODEL
C                   =3   MULTIPROCESS ASYMMETRIC INDIVIDUAL DIFFERENCES
C                        MODEL
C                   =4   GENERALIZED EUCLIDEAN MODEL
C     NDIM     NUMBER OF DIMENSIONS IN THE SOLUTION  (MAXIMUM)
C     NDMN     NUMBER OF DIMENSIONS IN THE SOLUTION  (MINIMUM)
C     NNC=1    IF WEIGHTS NOT CONSTRAINED TO BE POSITIVE
C              (DEFAULT IS NONNEGATIVITY RESTRICTIONS)
C     MAXIT    MAXIMUM NUMBER OF ITERATIONS (DEFAULT IS 30)
C     EPSI     CONVERGENCE CRITERION (DEFAULT IS 0.001)
C     STMIN    MINIMUM STRESS (DEFAULT IS .005)
C     NDIR     NUMBER OF GEMSCAL DIRECTIONS


C-----STMIN ADDED TO FOLLOWING STATEMENT 7/7/82
C-----NDIR ADDED TO FOLLOWING STATEMENT FOR GEMSCAL 8/5/82
C???      READ(IN,102)NWE,NDIM,NDMN,NNC,MAXIT,EPSI,STMIN,NDIR
C  102 FORMAT(5I4,2F8.0,I4)

      NWE=0
C???      NDIM=2
C???      NDMN=2
      NDIM=NDIM_MY
      NDMN=NDMN_MY
      NNC=0
      MAXIT=100
      EPSI=0.001
      STMIN=0.005
      NDIR=0

      IF(NWE.NE.1.AND.NWE.NE.2.AND.NWE.NE.2.AND.NWE.NE.4) NWE=0

C-----FOLLOWING STATMENTS ADDED FOR GEMSCAL 8/5/82
      IFLPDS=0
      IF(NWE.EQ.4) THEN
	IFLPDS=1
	NWE=1
      ENDIF
      IF(NWE.EQ.1.AND.NS.EQ.1) THEN
	WRITE(*,9910)
	WRITE(*,333)
      ENDIF
C-----FOLLOWING STATMENT ADDED FOR GEMSCAL 8/5/82
  333 FORMAT(8X,'THE INDIVIDUAL DIFFERENCES OR GENERALIZED EUCLIDEAN',
     1' MODEL',8X,'CAN NOT BE USED WITH ONLY ONE MATRIX.'/
     2 8X,'ANALYSIS CONTINUES WITH EUCLIDEAN MODEL (NWE=0)')
      IF(NWE.EQ.1.AND.NS.EQ.1) NWE=0
      IF(NWE.EQ.3.AND.NS.EQ.1) THEN
	WRITE(*,9910)
	WRITE(*,9238)
      END IF
 9238 FORMAT(8X,'THE INDIVIDUAL DIFFERENCES ASYMMETRIC MODEL CANNOT',
     1' BE USED WITH ONLY ONE MATRIX.'/
     2 8X,'ANALYSIS CONTINUES WITH THE ASYMMETRIC MODEL (NWE=2)')
      IF(NWE.EQ.3.AND.NS.EQ.1) NWE=2
      IF(NWC.NE.0.OR.(NWE.NE.0.AND.NWE.NE.2).OR.NS.NE.1)GOTO 17
      NWC=1
      WRITE(*,9910)
      WRITE(*,219)
  219 FORMAT(8X,'PROGRAM WILL CONTINUE TREATING THE DATA AS MATRIX',
     1' CONDITIONAL (NWC=1)')
   17 CONTINUE
      IF(NWE.LE.1.OR.NSIM.GT.1)GOTO 20
      WRITE(*,9910)
      WRITE(*,21)
   21 FORMAT(8X,'THE ASYMMETRIC MODEL REQUIRES ASYMMETRIC DATA.',
     1/8X,' THE MODEL HAS BEEN CHANGED AS FOLLOWS')
      NWE=NWE-2

   20 IF(NWC.NE.2.OR.NSIM.GE.2)GOTO 22
      WRITE(*,9910)
      WRITE(*,23)
   23 FORMAT(8X,'ROW CONDITIONAL DATA ARE NOT PERMITTED WITH SYMMETRIC',
     1' DATA.',
     2/8X,'THE DATA WILL BE TREATED AS MATRIX CONDITIONAL. (NWC=1)')
      NWC=1

   22 IF(NWE.GT.1.AND.NWC.EQ.2) THEN
      WRITE(*,9910)
      WRITE(*,9901)
      NWC=1
      END IF
 9901 FORMAT(8X,'ROW CONDITIONAL DATA ARE NOT PERMITTED WITH THE',
     1' ASYMMETRIC MODELS.'/
     28X,'ANALYSIS CONTINUES WITH MATRIX  CONDITIONAL DATA (NWC=1)')
      IF(NWE.GT.1.AND.NDIM.GT.5) NDIM=5
      IF(NDIM.EQ.0) NDIM=2
      IF(NDIM.GT.NDX) THEN
	WRITE(*,9910)
	WRITE(*,4003) NDX
      END IF
 4003 FORMAT(8X,'THE MAXIMUM DIMENSIONALITY MAY NOT EXCEED',I2)
      IF(NDIM.GT.NDX) NDIM=NDX
      NDX=NDIM
      NDXS=NDX**2
      NDXP=NDX+1
      IF(NDMN.LE.0.OR.NDMN.GT.NDIM)NDMN=NDIM
      IF(NDMN.EQ.1.AND.NWE.GT.0) THEN
	WRITE(*,9910)
	WRITE(*,9902)
      END IF
      IF(NDMN.EQ.1.AND.NWE.GT.0.AND.NDIM.EQ.1)NDIM=2
      IF(NDMN.EQ.1.AND.NWE.GT.0)NDMN=2
 9902 FORMAT(8X,'ONE-DIMENSIONAL WEIGHTED MODELS NOT PERMITTED.'/
     18X,'ANALYSIS CONTINUES WITHOUT A ONE-DIMENSIONAL SOLUTION.')
      ND=NDIM-NDMN+1
      NSC=1
      IF(NWC.EQ.0)NSC=0

      IF(NNC.NE.1) NNC=0


C     I/O OPTIONS
C     -----------

C     INITX    INITIAL STIMULUS COORDINATES
C                  =0   COMPUTE
C                  =1   COMPUTE AND PRINT
C                  =2   READ AND PRINT
C                  =3   READ, PRINT AND FIX
C     INITXC   INITIAL COLUMN STIMULUS COORDINATES
C                  =0   COMPUTE
C                  =1   COMPUTE AND PRINT
C                  =2   READ AND PRINT
C                  =3   READ, PRINT AND FIX
C     INITW    INITIAL SUBJECT WEIGHTS
C                  =0   COMPUTE
C                  =1   COMPUTE AND PRINT
C                  =2   READ AND PRINT
C                  =3   READ, PRINT AND FIX
C     INITWS   INITIAL STIMULUS WEIGHTS
C                  =0   COMPUTE
C                  =1   COMPUTE AND PRINT
C                  =2   READ AND PRINT
C                  =3   READ, PRINT AND FIX
C     NOULB    UPPER AND LOWER BOUNDS TO ESTIMATE MISSING DATA
C                  =0   YES (DEFAULT)
C                  =1   NO  (ALSCAL 4 METHOD)
C     ICNSTR   CONSTRAIN MISSING DATA (NOT FULLY IMPLEMENTED)
C     DEBUG    FULL DEBUGGING OUTPUT



C-----NOULB, ICNSTR, AND DEBUG ADDED TO FOLLOWING STATEMENT 7/7/82

      INITX=0
      INITXC=0
      INITW=0
      INITWS=0
      NOULB=0
      ICNSTR=0
      DEBUG=0

      ICNSTR=0


C-----ALGORITHMIC OPTIONS SECTION COMPRISED OF MOVED AND NEW STATEMENTS
      IF(MAXIT.EQ.0)MAXIT=MAXITD
      IF(EPSI.LE.0)EPSI=EPSID
      IF(STMIN.LE.0.0)STMIN=STMIND
      IF(NDMX.LE.0)NDMX=MIN0(NT,1000)
      IF(NDTYP.NE.3)NDMX=0

C----------------------------------------------------------------------

C     CALCULATION OF STORAGE REQUIREMENTS FOR ARRAYS
C     ----------------------------------------------

C    IN THIS SECTION WE CALCULATE THE NUMBER OF WORDS OF STORAGE
C  REQUIRED FOR VARIOUS ARRAYS USED BY THE PROBLEM. THIS IS DONE
C  BY ARRANGING THE ARRAYS END TO END AND ADDING THEIR LENGTHS
C  TOGETHER. IN THIS PROCESS AT EACH STEP THE STARTING PLACE OF
C  THE NEXT ARRAY IS CALCULATED BY ADDING THE LENGTH OF THE CURRENT
C  ARRAY TO ITS STARTING POINT. THUS THE STARTING POINT OF THE LAST
C  ARRAY PLUS ITS LENGTH GIVES THE SIZE OF THE STORAGE BLOCK NEEDED.
C     AS A CONVENTION THE STARTING POINT OF AN ARRAY IS GIVEN BY
C  THE NAME OF THE ARRAY PREFIXED BY J . THUS
C      JNEXT=JCURR + LCURR
C  WHERE LCURR IS USUALLY AN EXPRESSION IN THE PARAMETERS OF THE PRO-
C  BLEM. IF MORE THAN ONE ARRAY IS TO SHARE THE SAME STORAGE THEN
C      JNEXT=JCURR + (MAXIMUM OF THE LENGTHS OF ARRAYS SHARING STO-
C   RAGE WITH CURR) SO AS NOT TO OVERWRITE THE NEXT ARRAY.
C     THE STARTING POINTS OF THE ARRAYS ARE STORED IN THE COMMON
C  BLOCK /STARTS/.
C     THE FINAL CALCULATED NUMBER OF WORDS OF STORAGE REQUIRED IS SAVED
C  IN NWORDS.

C----------------------------------------------------------------------

      JX=1
C   X IS USED AS A 2*NB*NS ARRAY IN INNER
      JWA=JX+MAX0(NT,2*NB*NS)
      JWD=JWA+NBNBNS
      JXX=JWD+MAX0(NT,NDXS*NB,NBS*NDX,2778)
C   XX IS PASSED TO AM IN INIT AND USED THERE AS AN NS*NS ARRAY
      JDS=JXX+MAX0(NBS,NS*NS)
      JCFR=JDS+NBS
      LCFR=NB*NDXP
      JCFL=JCFR+LCFR
C  CFL GETS PASSED TO XEQ IN INIT AND FROM THERE IN TWO PIECES
C  TO U AND V IN CJEIG. FOR THIS CALL THE LENGTH OF U+V IS 4*NS.
C  OTHERWISE CFL IS USED AS AN NB*NDXP ARRAY
      JW=JCFL+MAX0(4*NS,LCFR)
      JWS=JW+NS*NDXP
      JZZ=JWS+LCFR
C  ZZ GETS PASSED FROM INSWM TO CJEIG AS SCRATCH SPACE OF LENGTH 2*NB.
C  OTHERWISE ZZ IS USED AS A SQUARE (NDIM*(NDIM+1))/2 ARRAY.
      JU11=JZZ+MAX0(2*NB,((NDIM+1)*NDIM/2)**2)
      JU12=JU11+NDXS
      JU22=JU12+NDXS
      JUB1=JU22+NDXS
      JUB2=JUB1+NDIM
      JXN=JUB2+NDIM
      JTR=JXN+NDIM
C  TR IS USED AS SCRATCH SPACE IN SEVERAL SUBROUTINES
      JFK=JTR+MAX0(NB,NS)
      JCV=JFK+NDXP**2
      LCV=((NDIM+1)*NDIM)/2
      JCW=JCV+LCV
C       LCW=LCV
      JPHSUB=JCW+LCV
      JPHSTI=JPHSUB+NS
      JNDSR=JPHSTI+NB
      JNAD=JNDSR
      IF(NWC.EQ.1.AND.NDTYP.EQ.4)JNAD=JNDSR+NS
      IF(NWC.EQ.2.AND.NDTYP.EQ.4)JNAD=JNDSR+NS*NB
      JIX=JNAD
      IF(NWC.EQ.1)JIX=JNAD+NS
      IF(NWC.EQ.2)JIX=JNAD+NS*NB
C-----IX, IY, AND IZ CHANGED TO INTEGER 7/8/82
C     THUS, THEIR LENGTHS ARE NO LONGER HALF OF THEIR DIMENSIONS.
      JIY=JIX+(NT+1)
      JIZ=JIY+(NDMX+1)
C-----NEXT STATEMENTS ADDED FOR EQUIVALENCE TO SAS VERSION 7/8/82
      JFTLN=JIZ+1
      IF(NDTYP.NE.3)GOTO 5100
      IF(NWC.EQ.1)JFTLN=JIZ+(NS+1)
      IF(NWC.EQ.2)JFTLN=JIZ+(NB*NS+1)
C-----SPACE FOR FTLN NOT BEING USED IN VERSION 4.04 AND UP
 5100 JDIST=JFTLN+2
      IF(NDTYP.NE.2.OR.NDEG.NE.1.OR.NWC.EQ.0)GOTO 510
      IF(NWC.EQ.1)JDIST=JFTLN+NS*2
      IF(NWC.EQ.2)JDIST=JFTLN+NS*NB*2
C-----FOLLOWING SEVEN STATEMENTS ADDED FOR GEMSCAL 8/5/82
  510 IF(IFLPDS.NE.1)GOTO 7381
      JPIJP=JDIST+NBS
      JXK=JPIJP+NBS
      NWORDS=JXK+NB*6
      RETURN
 7381 JPIJP=JDIST
      JXK=JDIST
      NWORDS=JDIST
      RETURN

C----------------------------------------------------------------------

  902 WRITE(*,9900)
      WRITE(*,9918)
 9918 FORMAT (8X,'NUMBER OF ROWS MUST EQUAL NUMBER OF COLUMNS FOR NON-RE
     1CTANGULAR DATA.')
      EOJ=1
      RETURN
  903 WRITE(*,9900)
      WRITE(*,9926)
 9926 FORMAT(8X,' NUMBER OF STIMULI LESS THAN 3'/)
      EOJ=1
      RETURN
  904 WRITE(*,9900)
      WRITE(*,9927)
 9927 FORMAT(8X,' TOTAL NUMBER OF MATRICES LESS THAN 1'/)
      EOJ=1
      RETURN
  905 WRITE(*,9925)
 9925 FORMAT(/' ALSCAL MESSAGE: ALL CARDS AND PROBLEMS READ.  NORMAL END
     1 OF JOB.')
       EOJ=1
       RETURN
 9900 FORMAT(/' ALSCAL FATAL ERROR:  COMPUTATIONS TERMINATED.')
 9910 FORMAT(/' ALSCAL WARNING:  INCONSISTENT CONTROL PARAMETERS.')
      END




C**********************************************************************

      SUBROUTINE DRIVER(AREA,EOJ,NROW_MY,NCOL_MY,XX_MY,X1_MY,X2_MY,
     *X3_MY,X4_MY,X5_MY,X6_MY,STRSS1)
C???      SUBROUTINE DRIVER(AREA,EOJ)
C     =================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/08/82

C     THIS SUBROUTINE CONTROLS THE FLOW OF THE PROGRAM

C   THE ARRAY "AREA" IS A BLOCK OF MEMORY WHICH IS ALLOCATED BY THE
C   MAIN ROUTINE, AND IS SUBDIVIDED HERE FOR USE IN THE SUBROUTINES

C----------------------------------------------------------------------

      INTEGER EOJ
      CHARACTER COPR*72
      DIMENSION AREA(1),XX_MY(NROW_MY,NCOL_MY)

C----------------------------------------------------------------------
C     INTEGER*2 IX(8000),IY(3000),IZ(525),ITEM(55,101)
C     DIMENSION FMT(20)
C     DIMENSION X(8000),WA(8000),WD(8000)
C     DIMENSION UA(875),U11(36),U12(36),U22(36),R(36),UB1(6),
C    1UB2(6),BK(6),WK(6),WK2(6),WK3(36),WK4(6),XN(6)
C     DIMENSION XX(1600),CFR(280),CFL(280),W(280),TR(40),FK(49)
C     DIMENSION WS(280),DS(1600),ZZ(225),CV(15),CW(15),ZX(280)
C     DIMENSION XEQ(160),NDSB(35),NDSR(525),NAD(525),FMRR(40)
C     EQUIVALENCE (UA(1),U11(1)),(UA(37),U12(1)),
C    1(UA(73),U22(1)),(UA(109),R(1)),(UA(145),UB1(1)),(UA(151),UB2(1)),
C    2(UA(157),BK(1)),(UA(163),WK(1)),(UA(169),WK2(1)),(UA(175),WK3(1)),
C    3(UA(211),WK4(1)),(UA(217),XN(1))
C     EQUIVALENCE (WD(1),ITEM(1,1)),(CFL(1),XEQ(1))
C     EQUIVALENCE (TR(1),FMRR(1)),(ZZ(1),ZX(1))
C----------------------------------------------------------------------

      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
C    +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
     +       /STARTS/JX,JWA,JWD,JXX,JDS,JCFR,JCFL,JW,JWS,JZZ,JTR,JFK,
     + JCV,JCW,JU11,JU12,JU22,JUB1,JUB2,
     + JXN,JPHSUB,JPHSTI,JIX,JIY,JIZ,JNDSR,JNAD,
     + JFTLN,JDIST,JPIJP,JXK
C-----JFTLN AND JDIST ADDED TO THE ABOVE STATEMENT 7/8/82
C-----JPIJP AND JXK ADDED 5MAY83 FWY

      CALL STEP1(AREA(JIX),AREA(JX),AREA(JWA),AREA(JXX),AREA(JDS),
     2 AREA(JNAD),IRET,XX_MY)
      IF(IRET.EQ.1) THEN
	EOJ=1
	RETURN
      ENDIF
      NPARM=NB
      IF(NWE.EQ.1.OR.NWE.EQ.3)NPARM=NPARM+NS
      IF(NWE.GT.1)NPARM=NPARM+NB

      DO 10 IJKL=1,ND
      NP=NPARM*(NDXP-IJKL)
      IF(NP.GT.NT)GOTO 7
      IF(NT.GE.2.5*NP)GOTO 9
      IF(NDTYP.LT.3)GOTO 9
C-PC-----------------------
C???      WRITE(*,9910)NP,NT
C--------------------------
C     WRITE(*,9910)NP
      GOTO 9
    7 WRITE(*,9920)NP,NT
C-----------------------
      STOP

    9 CALL STEP2(AREA(JIX),AREA(JIY),AREA(JIZ),AREA(JX),AREA(JWA),
     2 AREA(JWA),AREA(JWD),AREA(JCFR),AREA(JCFL),AREA(JW),AREA(JTR),
     3 AREA(JFK),AREA(JXX),AREA(JWS),AREA(JDS),AREA(JZZ),
     4 AREA(JCV),AREA(JCW),AREA(JNDSR),AREA(JNAD),IJKL)

      CALL STEP3(AREA(JIX),AREA(JIY),AREA(JIZ),AREA(JX),AREA(JWD),
     2 AREA(JWA),AREA(JWA),AREA(JWA),AREA(JU11),AREA(JU12),AREA(JU22),
     3 AREA(JUB1),AREA(JUB2),AREA(JXN),AREA(JXX),AREA(JCFL),AREA(JW),
     4 AREA(JTR),AREA(JFK),AREA(JWS),AREA(JNDSR),AREA(JNAD),IJKL,IRE)
      IF(IRE.EQ.1) GOTO 10

      CALL STEP3A(AREA(JIX),AREA(JIY),AREA(JIZ),AREA(JX),AREA(JWA),
     2 AREA(JWD),AREA(JXX),AREA(JCFL),AREA(JW),
     3 AREA(JWS),AREA(JNDSR),AREA(JNAD),AREA(JPHSUB),AREA(JPHSTI),
     4 STRSS1)

      CALL STEP4(AREA(JX),AREA(JWA),AREA(JXX),AREA(JCFL),
     2 AREA(JW),AREA(JTR),AREA(JWS),
     3 AREA(JDIST),AREA(JDS),AREA(JPIJP),AREA(JXK),AREA(JTR),IJKL,
     4 X1_MY,X2_MY,X3_MY,X4_MY,X5_MY,X6_MY)
C-----PRECEEDING LINE ADDED 5MAY83 FWY

   10 CONTINUE
      RETURN

C-PC-------------------------------------------------------------------
C  WARNING 9910 HAS BEEN CHANGED TO AVOID MISUNDERSTANDING.
C  THE WARNING IS PRINTED, IF
C  - THE MEASUREMENT LEVEL IS ORDINAL OR NOMINAL  .AND.
C  - THE NUMBER OF PARAMETERS IS < 2.5*NUMBER OF INPUT DATA.
 9910 FORMAT(//' ALSCAL WARNING:',
     1 I5,' PARAMETERS,',I5,' OBSERVATIONS',
     2/8X,'- THE NUMBER OF PARAMETERS BEING COMPUTED MAY BE TOO LARGE',
     3/8X,'- OR THE NUMBER OF OBSERVATIONS MAY BE TOO SMALL',
     4/8X,'FOR RELIABLE RESULTS.')
C----------------------------------------------------------------------
C9910 FORMAT(//' ALSCAL WARNING:  INCONSISTENT CONTROL PARAMETERS'/
C    1 8X,'THE NUMBER OF PARAMETERS BEING COMPUTED (',I5,
C    2 ') MAY BE TOO LARGE',/8X,'FOR RELIABLE RESULTS')

 9920 FORMAT(//' ALSCAL FATAL ERROR:  COMPUTATIONS TERMINATED'/
     1 8X,'THE NUMBER OF PARAMETERS (',I5,') EXCEEDS THE NUMBER',
     2' OF OBSERVATIONS(',I5,').')
      END


C**********************************************************************

      SUBROUTINE STEP1(IX,X,WA,XX,DS,NAD,IER,XX_MY)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 01JUN83

C     DATA INPUT        ROUTINE

C----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 7/8/82
      INTEGER IX(NT)
      DIMENSION X(NT),WA(NB,NB,NS),XX(NB,NB),DS(NB,NB),
     2 NAD(NS,NB),XX_MY(NB,NB)
      CHARACTER*80,FMT,COPR*72
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
C-----THE FOLLOWING THREE STATEMENTS ADDED 7/9/82
      DOUBLE PRECISION CUT,STMIN,AX
      INTEGER DEBUG,ICNSTR,NOULB
      COMMON /PRMBLK/CUT,STMIN,DEBUG,ICNSTR,NOULB

C----------------------------------------------------------------------

C  X CONTAINS THE ACTIVE DATA ( WITH THE MISSING DATA FLAGGED)
C  IN STANDARD ORDER WITH SIMILARITIES CONVERTED TO DISSIMILARITIES
C  AND SCALED TO KNOWN UNITS.
C  WA IS THE SAME AS X EXCEPT
C    1) HAS MEANS IN PLACE OF MISSING DATA
C    2) HAS ADDITIVE CONSTANT ESTIMATED
C    3) IS ALWAYS SYMMETRIC
C  DS CONTAINS THE SUM,OVER SUBJECTS, OF INFORMATION IN X
C  IX CONTAINS THE MISSING DATA PATTERN 1=MISSING, 0=ACTIVE
C  XX IS USED AS TEMPORARY SPACE TO CONSTRUCT EACH MATRIX OF
C        DATA TO BE USED BY THE INITIALIZATION ROUTINES.
C        XX IS PROCESSED SO THAT:
C           1) MISSING DATA ESTIMATED BY MEAN UPPER-LOWER BOUND
C           2) ASYMMETRIC DATA ARE SYMMETRIZED
C           3) SIMILARITIES ARE MADE DISSIMILARITIES
C           4) CONDITIONAL DATA ARE NORMALIZED
C  WA CONTAINS THE SAME INFORMATION AS XX, EXCEPT
C           1) FOR ALL MATRICES
C           2) ADDITIVE CONSTANT IS ESTIMATED
C           3) RESULT IS SQUARED

C  CERTAIN IMPORTANT DATA FLAGS AND THEIR VALUES ARE:
C      NSIM: 0,1=SYM 2,3=ASYM 4,5=RECT EVEN=DISSIM ODD=SIMILARITIES
C      NWC:  0=UNCONDITIONAL 1=MATCON 2=ROWCON

C----------------------------------------------------------------------

      NAB=0
      NX=2
      IF(NSIM.GT.1) NX=1
      I1=NB
      DO 4001 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 4001 J=1,I1
 4001 DS(I,J)=0.0
      NN=0
      IER=0
      REWIND NDPP
      DO 10 L=1,NS
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/8/82
      IF(DEBUG.NE.0)WRITE(*,4040)
 4040 FORMAT(/' RAW DATA AS READ FROM DATA SET')
      NTEMP2=0
      XX(1,1)=0.0
      DO 11 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      IF(NSIM.GT.3.AND.I.LE.NCOL)GOTO 7
      J1=I1
      IF(NSIM.GT.3)J1=NCOL
      DO 77,J=1,J1
 77      XX(I,J)=XX_MY(I,J)
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/8/82
      IF(DEBUG.NE.0)WRITE(*,40)(XX(I,J),J=1,J1)
   40 FORMAT(6X,10G12.4)
      IF(NSIM.LT.4)GOTO 9
      J1=J1+1
      DO 6 J=J1,I1
C-----XX(I,J) CHANGED FROM 0.0 7/8/82
    6 XX(I,J)=CUT-1.0
      GOTO 9
    7 DO 8 J=1,NB
C-----XX(I,J) CHANGED FROM 0.0 7/8/82
    8 XX(I,J)=CUT-1.0
    9 WRITE(NDPP)(XX(I,J),J=1,I1)
      IF(NSIM.GT.1)GOTO 11
      DO 12 J=1,I1
   12 XX(J,I)=XX(I,J)
   11 CONTINUE
C--------------------------------------------------------------------
C-----THE FOLLOWING TEN LINES ADDED 7/9/82
C*****THE SAS VERSION READS IN A TRANSPOSED MATRIX AND THUS MUST
C     RETRANSPOSE IT.  THIS IS NOT THE CASE WITH THE STAND ALONE
C     VERSION, SO THIS SECTION OF CODE IS NOT NEEDED.*****

C     IF THE MATRIX IS NOT SYMMETRIC, TRANSPOSE IT

C     IF(NSIM.LT.2)GOTO 21
C     DO 20 I=2,NB
C     IM1=I-1
C     DO 20 J=1,IM1
C     TEMP=XX(I,J)
C     XX(I,J)=XX(J,I)
C  20 XX(J,I)=TEMP
C--------------------------------------------------------------------

C     CHECK FOR MISSING DATA AND REPLACE BY MEANS

 21   DO 5 I=1,NB
    5 XX(I,I)=0.0
      FMR=0.0
      N=NN
      IFLXER=0
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/9/82
      IF(DEBUG.GT.0)WRITE(*,4848)
 4848 FORMAT(/' MISSING DATA PATTERN')
      DO 2337 I=NX,NB
      NTEMP1=0
      T=0.0
      K=N
      IF(NSIM.LE.1)I1=I-1
      DO 337 J=1,I1
      N=N+1
C---THE FOLLOWING STATEMENT IS REPLACED BY THE NEXT TO FIX DATA CUTOFF
C     IF(XX(I,J).GT.0.0.OR.I.EQ.J)GOTO 338
      IF(XX(I,J).GE.CUT.OR.I.EQ.J)GOTO 338
      IX(N)=1
C-----THE FOLLOWING THREE STATEMENTS ADDED 7/9/82
      IF(ICNSTR.NE.1.OR.NSIM.LT.4)GOTO 337
      IF(I.GT.NCOL.AND.J.LE.NCOL)GOTO 337
      IX(N)=0
      GOTO 337
  338 NTEMP1=NTEMP1+1
      T=T+XX(I,J)
      IX(N)=0
  337 CONTINUE
      IF(NTEMP1.EQ.1.AND.NWC.EQ.2) IFLXER=1
       IF(NTEMP1.EQ.0.AND.NWC.EQ.2)WRITE(*,9001)L,I
 9001 FORMAT(/' ALSCAL WARNING: ALL DATA MISSING FOR MATRIX',I4,
     *'STIMULUS',I4)
      NTEMP2=NTEMP2+NTEMP1
      FMR=FMR+T
C-----THE FOLLOWING TWO STATEMENTS REPLACED BY THE NEXT THREE 7/9/82
C     IF(NWC.NE.2.OR.IFLXER.EQ.1)GOTO 2337
C     NAD(L,I)=NTEMP1
      IF(NWC.EQ.2)NAD(L,I)=NTEMP1
      IF(NWC.NE.2)GOTO 2337
      IF(NWC.EQ.2.AND.NTEMP1.EQ.1)GOTO 2337
      T=T/(NTEMP1-1)
      DO 340 J=1,I1
      K=K+1
  340 IF(IX(K).NE.0)XX(I,J)=T
 2337 CONTINUE
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/9/82
      IF(DEBUG.GT.0)WRITE(*,4789)(IX(IIIIII),IIIIII=1,N)
 4789 FORMAT(10(1X,10I1))
      NMISS=NC2-NTEMP2
      IF(NSIM.GT.3)NMISS=NMISS-NCOL*NB-(NB-NCOL)**2+NB
      IF(NMISS.GT.0) THEN
      WRITE(*,5318)L,NMISS
 5318 FORMAT(/' ALSCAL MESSAGE:  MATRIX', I4,' HAS',I5,' MISSING OBSERVA
     1TIONS.')
      END IF
      IF(NTEMP2.EQ.0)GOTO 901
      IF(IER.NE.0)GOTO 10
      NAB=NAB+NTEMP2
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/9/82
      IF(NSIM.LE.1)FMR=FMR/NTEMP2
      IF(NSIM.GT.1)FMR=FMR/(NTEMP2-NB)
      IF(NWC.EQ.2.AND.IFLXER.EQ.0)GOTO 341
C-----THE FOLLOWING STATEMENT ADDED 7/9/82
      IF(NWC.EQ.2)GOTO 342
      IF(NWC.EQ.1)NAD(L,1)=NTEMP2
      IF(NMISS.EQ.0.AND.NSIM.LT.4)GOTO 341
      N=NN
      DO 339 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 339 J=1,I1
      N=N+1
      IF(IX(N).NE.0)XX(I,J)=FMR
C-----THE FOLLOWING STATEMENT ADDED 7/9/82
      IF(IX(N).NE.0.AND.NSIM.LE.1)XX(J,I)=FMR
  339 CONTINUE
C-----THE FOLLOWING SEVEN STATEMENTS ADDED 7/9/82
      GOTO 341
  342 N=NN
      DO 343 I=1,NB
      DO 343 J=1,NB
      N=N+1
      IF(NAD(L,I).GT.1)GOTO 343
      IF(IX(N).NE.0)XX(I,J)=FMR
  343 CONTINUE

C     IF DATA ARE SIMILARITY (NSIM=1 OR 3) CONVERT THEM
C     INTO DISSIMILARITY

  341 IF(NSIM/2*2.EQ.NSIM) GOTO 114
      IF(NWC.EQ.2) GOTO 2005
      AX=XX(NB,1)
      DO 115 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 115 J=1,I1
  115 IF(XX(I,J).GT.AX) AX=XX(I,J)
      DO 117 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 116 J=1,I1
  116 XX(I,J)=AX-XX(I,J)
  117 XX(I,I)=0.0
      XX(1,1)=0.0
      GOTO 114
 2005 DO 2006 I=1,NB
      AX=XX(I,1)
      DO 2007 J=1,NB
 2007 IF(XX(I,J).GT.AX) AX=XX(I,J)
      DO 2008 J=1,NB
 2008 XX(I,J)=AX-XX(I,J)
      XX(I,I)=0.0
 2006 CONTINUE

C     INITIAL SCALING WITHIN SUBJECT

C-----THE FOLLOWING FIVE STATMENTS WERE ADDED 7/8/82
  114 IF(DEBUG.EQ.0) GOTO 3200
      WRITE(*,3399)
 3399 FORMAT(/' DATA WITH MISSING ESTIMATED BY ROW OR MATRIX MEAN')
      DO 3388 I=1,NB
 3388 WRITE(*,6789)(XX(I,J),J=1,NB)
 3200 AX=1.0
      IF(NSC.EQ.0) GOTO 5091
      AX=0.0
      DO 5092 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 5092 J=1,I1
 5092 AX=AX+XX(I,J)**2
      AX=DSQRT(AX/NC2)
      DO 5093 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 5093 J=1,I1
 5093 XX(I,J)=XX(I,J)/AX
C-----THE FOLLOWING FIVE STATEMENTS ADDED 7/9/82
      IF(NSIM.GT.1)GOTO 5091
      DO 5094 I=NX,NB
      I1=I-1
      DO 5094 J=1,I1
 5094 XX(J,I)=XX(I,J)

C     STORE ORIGINAL DATA

 5091 DO 4002 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 4002 J=1,I1
 4002 DS(I,J)=DS(I,J)+XX(I,J)
C-----THE FOLLOWING FIVE STATEMENTS WERE ADDED 7/8/82
      IF(DEBUG.EQ.0) GOTO 3210
      WRITE(*,3377)
 3377 FORMAT(/' NORMALIZED DATA AND MEANS')
      DO 3366 I=1,NB
 3366 WRITE(*,6789)(XX(I,J),J=1,NB)
 3210 N=NN
      DO 2009 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 2009 J=1,I1
      N=N+1
      X(N)=XX(I,J)
      IF(IX(N).NE.0) X(N)=BIG
 2009 CONTINUE
C-----MOST OF THE STATEMENTS IN THE FOLLOWING SECTION WERE ADDED 7/8/82

C     CHECK TO SEE IF MISSING DATA EXIST OR IT THE
C     DATA ARE RECTANGULAR (IMPLIED MISSING DATA).  IF SO
C     ESTIMATE BY MEAN OF UPPER AND LOWER BOUNDARIES
C     GENERATED BY LINE OF SIGHT METHOD.

      IF(NSIM.GT.3.OR.NMISS.GT.0) GOTO 7200
      IF(NSIM.LE.1) GOTO 3245

C     COME HERE WHEN THE DATA ARE ASYMMETRIC AND COMPLETE
C     TO GENERATE SYMMETRIC INICON DATA.

      DO 7150 I=2,NB
      IM1=I-1
      DO 7150 J=1,IM1
      XX(I,J)=(XX(I,J)+XX(J,I))/2.0
 7150 XX(J,I)=XX(I,J)
      GOTO 3245

C     COME HERE WHEN THERE IS MISSING DATA OR WHEN THE
C     DATA ARE RECTANGULR TO FLAG THE MISSING CELLS.

 7200 N=NN
      DO 135 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 135 J=1,I1
      N=N+1
      IF(IX(N).EQ.0)GOTO 135
      XX(I,J)=CUT-XX(I,J)
      IF(NSIM.LE.1)XX(J,I)=XX(I,J)
  135 CONTINUE
      IF(NSIM.LT.4)GOTO 130
      NCOLP1=NCOL+1
      DO 120 I=NCOLP1,NB
      DO 120 J=1,NCOL
  120 XX(J,I)=XX(I,J)

C     SYMMETRIZE THE ASYMMETRIC MATRIX

  130 IF(DEBUG.EQ.0)GOTO 3211
      WRITE(*,3370)
 3370 FORMAT(/' DATA AFTER MISSING FLAGGED')
      DO 3360 I=1,NB
 3360 WRITE(*,6789)(XX(I,J),J=1,NB)
 3211 IF(NSIM.LE.1) GOTO 2010
      NA=2
      IF(NSIM.GE.4)NA=NCOL+2
      NC=NA-1
      DO 2011 I=NA,NB
      IMI=I-1
      DO 2011 J=NC,IMI
      IF(XX(I,J).LT.CUT.OR.XX(J,I).LT.CUT)GOTO 2002
 2001 XX(I,J)=(XX(I,J)+XX(J,I))*.5
      XX(J,I)=XX(I,J)
      GOTO 2011
 2002 IF(XX(I,J).LT.CUT.AND.XX(J,I).LT.CUT) GOTO 2001
      IF(XX(I,J).LT.CUT)XX(I,J)=XX(J,I)
      IF(XX(J,I).LT.CUT)XX(J,I)=XX(I,J)
 2011 CONTINUE
 3821 IF(DEBUG.EQ.0)GOTO 2010
      WRITE(*,1355)
 1355 FORMAT(/' SYMMETRIC NORMALIZED DATA')
      DO 232 I=1,NB
  232 WRITE(*,6789)(XX(I,J),J=1,NB)
 6789 FORMAT(1X,15F8.2)

C     UNLESS THE USER HAS REQUESTED, VIA THE NOULB PARAMETER,
C     REPLACE MISSING DATA ESITMATES WITH MEAN OF UPPER AND
C     LOWER BOUNDS ON DISTANCES WHEREVER POSSIBLE.
C     IF NOT POSSIBLE OR IF REQUESTED USE MEAN CALCULATED ABOVE.

 2010 IF(NOULB.NE.0)GOTO 1492
      DO 150 I=2,NB
      IM1=I-1
      DO 150 J=1,IM1
      IF(XX(J,I).GT.CUT)GOTO 150
      MISALL=1
      BIGDIF=0.0
      SMLSUM=10.0E10
      DO 330 II=1,NB
      IF(II.EQ.I.OR.II.EQ.J)GOTO 330
      IA=I
      IB=J
      IC=II
      ID=II
      IF(IC.GE.IA)GOTO 320
      IA=II
      IC=I
  320 IF(ID.GE.IB)GOTO 325
      IB=II
      ID=J
  325 IF(XX(IC,IA).LE.CUT.OR.XX(ID,IB).LE.CUT)GOTO 330
      MISALL=0
      SUM=    XX(IC,IA)+XX(ID,IB)
      DIF=ABS(XX(IC,IA)-XX(ID,IB))
      IF(SMLSUM.GT.SUM)SMLSUM=SUM
      IF(BIGDIF.LT.DIF)BIGDIF=DIF
  330 CONTINUE
      IF(MISALL.EQ.0)XX(J,I)=(SMLSUM+BIGDIF)*.5
  150 CONTINUE
 1492 DO 160 I=2,NB
      IM1=I-1
      DO 160 J=1,IM1
      IF(XX(J,I).LT.CUT)XX(J,I)=CUT-XX(J,I)
  160 XX(I,J)=XX(J,I)
      IF(DEBUG.EQ.0)GOTO 3245
      WRITE(*,1354)
 1354 FORMAT(/' MISSING DATA NOW ESTIMATED BY MEAN U/L BOUNDS')
      DO 233 I=1,NB
  233 WRITE(*,6789)(XX(I,J),J=1,NB)
C-----END OF ADDED SECTION 7/8/82

C     ESTIMATION OF ADDITIVE CONSTANT

 3245 ALPH1=0.0
      IF(NDTYP.EQ.1)GOTO 19
      ALPH1=XX(2,1)
      DO 14 I=2,NB
      IMI=I-1
      DO 14 J=1,IMI
   14 IF(XX(I,J).LT.ALPH1) ALPH1=XX(I,J)
      DO 15 I=1,NB
      DO 15 J=2,NB
      IF(J.EQ.I) GOTO 15
      J1=J-1
      DO 16 K=1,J1
      IF(K.EQ.I) GOTO 16
      TRAN=XX(I,J)+XX(I,K)-XX(J,K)
      IF(TRAN.LT.ALPH1)ALPH1=TRAN
   16 CONTINUE
   15 CONTINUE


   19 WA(1,1,L)=0.0
      DO 18 I=2,NB
      WA(I,I,L)=0.0
      IMI=I-1
      DO 18 J=1,IMI
      WA(I,J,L)=(XX(I,J)-ALPH1)**2
   18 WA(J,I,L)=WA(I,J,L)
C-----THE FOLLOWING SEVEN LINES ADDED 7/8/82
      IF(DEBUG.EQ.0)GOTO 10
      WRITE(*,1382)
 1382 FORMAT(/' DATA SQUARED WITH ADDITIVE CONSTANT ESTIMATED',
     1' (INICON DATA)')
      DO 288 I=1,NB
  288 WRITE(*,6789)(WA(I,J,L),J=1,NB)
      GOTO 10
  901 IF(IER.EQ.0)WRITE(*,9900)
      WRITE(*,9002)L
9002  FORMAT(8X,' ALL DATA MISSING FOR MATRIX',I4)
 9900 FORMAT(/' ALSCAL FATAL ERROR:  COMPUTATIONS TERMINATED.')
      WRITE(*,9001)L,I
      IER=1
   10 NN=NN+NC2

      IF(IER.NE.0)RETURN
      REWIND NDP
      WRITE(NDP)X
      WRITE(NDP) WA
      REWIND NDP
C-----THE FOLLOWING STATEMENT ADDED 7/8/82
      IF(ICNSTR.EQ.1)NSIM=NSIM-2
      RETURN
      END



C*********************************************************************

      SUBROUTINE STEP2(IX,IY,IZ,X,WA,WC,WD,CFR,CFL,W,TR,FK,
C     ================

     1 XX,WS,DS,ZZ,CV,CW,NDSR,NAD,IJKL)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/21/82

C  THIS ROUTINE OBTAINS STARTING CONFIGURATIONS AND INITIAL OPTIMAL
C  SCALING TO INITIATE THE ITERATIVE PROCESS

C    WA & WC  REFER TO THE SAME ARRAY IN THE CALLING ROUTINE.
C    THIS  ARRAY  IS  PASSED IN AS DIFFERENT PARAMETERS TO PERMIT
C    REFERENCING THE ARRAYS USING DIFFERENT SHAPE PARAMETERS.

C-----------------------------------------------------------------------

      DOUBLE PRECISION T,TT,TRT,RNB,RNS,RNBS
C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 7/8/82
      INTEGER IX(NT),IY(1),IZ(1)
      DIMENSION X(NT),WA(NB,NB,NS),WC(1),WD(1),XX(NB,NB)
      DIMENSION CFR(NB,1),CFL(NB,1),W(NS,1),TR(1),FK(NDXP,1)
      DIMENSION WS(NB,1),DS(NB,NB),ZZ(1,1),CV(1),CW(1)
C   THE ARRAY ZZ IS NOT USED IN STEP2 BUT PASSED THROUGH TO INSWM
C   AND ITS DIMENSIONS HERE ARE DUMMY.
      DIMENSION NDSR(NS,NB),NAD(NS,NB)
      CHARACTER*80,FMT,COPR*72
      COMMON /CCOPR/COPR
C-----THE FOLLOWING THREE STATEMENTS ADDED 7/14/82
      DOUBLE PRECISION CUT,STMIN
      INTEGER DEBUG,ICNSTR,NOULB
      COMMON /PRMBLK/CUT,STMIN,DEBUG,ICNSTR,NOULB
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
     +       /INICON/INITX,INITW,INITWS,INITXC

C----------------------------------------------------------------------

C     COMPUTE INITIAL COORDINATES AND WEIGHTS
C     ---------------------------------------

      ISTFOR=1
      IF(NSIM.GT.3)ISTFOR=2
      RNS=1./NS
      RNB=1./NB
      RNBS=1./NBS
      STRSO=BIG
      NDIM=NDXP-IJKL

C     IF DESIRED, READ INITIAL STIMULUS COORDINATES

      LLL=0
      IF(NSIM.LT.4)GOTO 3200
      IF(INITX.GE.2.AND.INITXC.GE.2)GOTO 3201
      GOTO 446
 3200 IF(INITX.LT.2)GOTO 446
 3201 READ(IN,447)FMT
  447 FORMAT(A80)
      K=1
      IF(NSIM.GT.3)K=NCOL+1
      DO 448 I=K,NB
  448 READ(IN,FMT)(CFL(I,J),J=1,NDIM)
      GOTO 3026

C     OTHERWISE COMPUTE  PRODUCT MOMENT MATRICES

  446 LLL=LLL+1

      DO 24 J=1,NB
      DO 24 I=J,NB
   24 XX(I,J)=0.0
      TT=0.0
      DO 19 L=1,NS
      TRT=0.0
      DO 20 I=1,NB
      T=0.0
      DO 21 J=1,NB
   21 T=T+WA(J,I,L)
      TRT=TRT+T
   20 TR(I)=T*RNB
      TRT=TRT*RNBS
      T=0.0
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/12/82
      IF(DEBUG.GT.0)WRITE(*,2144)L
 2144 FORMAT(' SCALAR PRODUCTS FOR SUBJECT',I4)
      DO 222 I=2,NB
	DO 22 J=1,I-1
	  WAJIL=(WA(J,I,L)-TR(I)-TR(J)+TRT)*(-.5)
	  WA(J,I,L)=WAJIL
	  WA(I,J,L)=WAJIL
	  T=T+WAJIL*WAJIL*2
   22   CONTINUE
	WAIIL=TR(I)-.5*(WA(I,I,L)+TRT)
	WA(I,I,L)=WAIIL
	T=T+WAIIL*WAIIL
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/12/82
	IF(DEBUG.GT.0)WRITE(*,2145)(WA(I,J,L),J=1,I)
 2145   FORMAT(15F8.2)
  222 CONTINUE
      WA11L=TR(1)-.5*(WA(1,1,L)+TRT)
      WA(1,1,L)=WA11L
      T=T+WA11L*WA11L
      TT=TT+T

C     NORMALIZATION (WITHIN SUBJECT)

      IF(NSC.EQ.0) GOTO 19
C-----THE FOLLOWING TWO STATEMENTS ADDED 7/12/82
      IF(DEBUG.GT.0)WRITE(*,2143)T
 2143 FORMAT(' IN STEP2: T = ',F20.7)
      T=DSQRT(NBS/T)
      DO 122 J=1,NB
      DO 122 I=J,NB
      WA(I,J,L)=WA(I,J,L)*T
      WA(J,I,L)=WA(I,J,L)
  122 XX(I,J)=XX(I,J)+WA(I,J,L)
   19 CONTINUE

C     NORMALIZATION (ACROSS SUBJECTS)

      IF(NSC.NE.0)GOTO 126
      TT=DSQRT(NBNBNS/TT)
      DO 123 L=1,NS
      DO 123 J=1,NB
      DO 123 I=J,NB
      WA(I,J,L)=WA(I,J,L)*TT
      WA(J,I,L)=WA(I,J,L)
  123 XX(I,J)=XX(I,J)+WA(I,J,L)
  126 DO 23 J=1,NB
      DO 23 I=J,NB
      XX(I,J)=XX(I,J)*RNS
   23 XX(J,I)=XX(I,J)
C---NWT AND NWW ARE NOT NEEDED; A LATER IF CHECK IS CHANGED
C     NWT=0
C     NWW=0
      FMX=0.0
      IF(NWE/2*2.EQ.NWE) GOTO 25

C     COMPUTE COORDINATES AND WEIGHTS FOR THE WEIGHTED EUCLIDIAN MODEL
C     BY USING THE SCHONEMANN-DE LEEUW METHOD

      CALL INIT(W,CFL,CFR,NB,NS,NDIM,XX,WA,TR,FK,WD,CFL,XX,WS,NDX,NDXP)

      DO 855 J=1,NDIM
      DO 855 I=1,NS
  855 IF(W(I,J).LT.FMX) FMX=W(I,J)
      IF(FMX.EQ.0.0)GOTO 867
      DO 866 J=1,NDIM
      DO 866 I=1,NS
  866 W(I,J)=W(I,J)-FMX
C---NWT IS NOT NEEDED
C     NWT=1
  867 T=0.0
      DO 517 J=1,NDIM
      DO 517 I=1,NS
  517 T=T+W(I,J)**2
      T=DSQRT((NDIM*NS)/T)
      DO 518 J=1,NDIM
      DO 518 I=1,NS
  518 W(I,J)=W(I,J)*T
      IF(NWE.EQ.3) GOTO 4004
      DO 4005 J=1,NDIM
      DO 4005 I=1,NB
 4005 WS(I,J)=1.0
      GOTO 26

C     COMPUTE COORDINATES FOR THE UNWEIGHTED EUCLIDIAN MODEL
C     BY USING TORGERSON'S METHOD

   25 CALL CJEIG(XX,CFL,CFR,NB,NDIM+1,FK,WS ,TR,1,NDXP)
C   HERE WS IS BEING PASSED TO BE USED AS SCRATCH SPACE IN CJEIG

      DO 41 J=1,NDIM
      TR(J)=SQRT(TR(J))
      DO 4080 I=1,NB
      CFL(I,J)=CFL(I,J)*TR(J)
 4080 WS(I,J)=1.0
      DO 41 I=1,NS
   41 W(I,J)=1.0
      IF(NWE.EQ.0) GOTO 26

C     COMPUTE COORDINATES AND WEIGHTS FOR THE ASYMMETRIC
C     EUCLIDEAN MODEL

      DO 4010 J=1,NDIM
 4010 CW(J)=1.0
      GOTO 4009
 4004 DO 4007 J=1,NDIM
      T=0.0
      DO 4008 I=1,NS
 4008 T=T+W(I,J)
 4007 CW(J)=T*RNS
 4009 DO 4053 I=1,NB
      DS(I,I)=0.0
      IF(NSIM.GT.1) GOTO 4053
      DO 4052 J=1,I
 4052 DS(J,I)=DS(I,J)
 4053 CONTINUE
      CALL INSWM(DS,CFL,CFR,WS,XX,TR,CV,CW,FK,ZZ,WD,NADCT,NB,NDIM,
     1 (NDIM*(NDIM+1))/2,NDX,NDXP,NS)
      DO 4013 J=1,NDIM
      DO 4013 I=1,NB
 4013 IF(WS(I,J).LT.FMX) FMX=WS(I,J)
      IF(FMX.GE.0.0)GOTO 4011
      DO 4014 J=1,NDIM
      DO 4014 I=1,NB
 4014 WS(I,J)=WS(I,J)-FMX
C---NWW IS NOT NEEDED
C     NWW=1
 4011 T=0.0
      DO 4015 J=1,NDIM
      DO 4015 I=1,NB
 4015 T=T+WS(I,J)**2
      T=DSQRT((NDIM*NB)/T)
      DO 4016 J=1,NDIM
      DO 4016 I=1,NB
 4016 WS(I,J)=WS(I,J)*T

C  ALL PATHS OF THE PROGRAM FROM ABOVE COME TOGETHER TO STATEMENTS 26
C     OR 3026

   26 IF(LLL.EQ.2)RETURN

C     READ OR PRINT INITIAL CONFIGURATION AND/OR INITIALIZE
C     INITIAL WEIGHTS WHEN INITIAL COORDINATES HAVE BEEN READ

      IF(INITX.GE.2)GOTO 3201
 3026 IF(LLL.EQ.2)RETURN
      IF(INITXC.LT.2.OR.NSIM.LT.3)GOTO 3126
      READ(IN,447)FMT
      DO 3030 I=1,NCOL
 3030 READ(IN,FMT)(CFL(I,J),J=1,NDIM)
 3126 IF(INITW.LT.2)GOTO 450
      READ(IN,447)FMT
      DO 449 I=1,NS
  449 READ(IN,FMT)(W(I,J),J=1,NDIM)
      GOTO 452
  450 IF(INITX.LT.2)GOTO 452
      DO 451 J=1,NDIM
      DO 451 I=1,NS
  451 W(I,J)=1.0
  452 IF(INITWS.LT.2)GOTO 454
      READ(IN,447)FMT
      DO 453 I=1,NB
  453 READ(IN,FMT)(WS(I,J),J=1,NDIM)
      GOTO 456
  454 IF(INITX.LT.2)GOTO 456
      DO 455 J=1,NDIM
      DO 455 I=1,NB
  455 WS(I,J)=1.0
C---INITXC ADDED TO THE FOLLOWING STATEMENT 7/6/82
  456 IF(INITX.GT.0.OR.INITW.GT.0.OR.INITWS.GT.0.OR.INITXC.GT.0) THEN
	WRITE(*,457)
  457   FORMAT(' INITIAL CONFIGURATION'//)
      ENDIF
C---INITXC ADDED TO THE FOLLOWING STATEMENT 7/6/82
      IF(INITX.EQ.0.AND.INITXC.EQ.0)GOTO 461
      WRITE(*,458)(I,I=1,NDIM)
  458 FORMAT(' INITIAL STIMULUS SPACE'/
     1T29,'DIMENSION'/5X,'STIMULUS',6I12)
      IF(NSIM.GT.3)WRITE(*,470)
  470 FORMAT(/'      COLUMN')
      DO 459 I=1,NB
      M=I
      IF(NSIM.GT.3.AND.I.GT.NCOL)M=I-NCOL
      NCOL1=NCOL+1
      IF(NSIM.GT.3.AND.I.EQ.NCOL1)WRITE(*,471)
  471 FORMAT(/' ',8X,'ROW')
  459 WRITE(*,460)M,(CFL(I,J),J=1,NDIM)
  460 FORMAT(I11,5X,6F12.4)
  461 IF(INITW.EQ.0)GOTO 464
      WRITE(*,462)(I,I=1,NDIM)
  462 FORMAT(/' INITIAL SUBJECT WEIGHTS'/
     1T29,'DIMENSION'/5X,'SUBJECT',6I12)
      DO 463 I=1,NS
  463 WRITE(*,460)I,(W(I,J),J=1,NDIM)
  464 IF(INITWS.EQ.0)GOTO 467
      WRITE(*,465)(I,I=1,NDIM)
  465 FORMAT(/' INITIAL STIMULUS WEIGHTS'/
     1T29,'DIMENSION'/5X,'STIMULUS',6I12)
      DO 466 I=1,NB
  466 WRITE(*,460)I,(WS(I,J),J=1,NDIM)


C---------------------------------------------------------------------

C     PREPARE QUALITATIVE DATA FOR OPTIMAL SCALING
C     --------------------------------------------


 467  IF(IJKL.GT.1)GOTO 30
      IZ(1)=1
      N2=0
      IF(NDTYP-3) 30,31,144
   31 IF(NWC.EQ.1) GOTO 431
      IF(NWC.EQ.2) GOTO 2012

C     UNCONDITIONAL ORDINAL DATA

      DO 33 K=1,NT
      IX(K)=K
   33 WC(K)=X(K)
      CALL SHEL9(WC,IX,NT)
      CALL BLOC2(WC,IY,IZ,NT)
      GOTO 30

C     MATRIX CONDITIONAL ORDINAL DATA

  431 DO 432 L=1,NS
      DO 433 K=1,NC2
      M=N2+K
      IX(M)=K
  433 WC(M)=X(M)
      N=N2+1
      CALL SHEL9(WC(N),IX(N),NC2)
      IF(L.EQ.1)GOTO 524
      IZ(L)=IZ(L-1)
  524 CALL BLOC2(WC(N),IY,IZ(L),NC2)
  432 N2=N2+NC2
      GOTO 30

C     ROW CONDITIONAL ORDINAL DATA

 2012 N1=0
      DO 2013 L=1,NS
      DO 2014 I=1,NB
      DO 2015 J=1,NB
      M=N2+J
      IX(M)=J
 2015 WC(M)=X(M)
      N=N2+1
      CALL SHEL9(WC(N),IX(N),NB)
      N1=N1+1
      IF(N1.EQ.1) GOTO 2016
      IZ(N1)=IZ(N1-1)
 2016 CALL BLOC2(WC(N),IY,IZ(N1),NB)
 2014 N2=N2+NB
 2013 CONTINUE
      GOTO 30

  144 IF(NWC.EQ.1) GOTO 573
      IF(NWC.EQ.2) GOTO 2017

C     UNCONDITIONAL NOMINAL DATA

      NDCT=0
      DO 145 L=1,NT
      IF(IX(L).EQ.1) GOTO 1145
      IF(NDCT.EQ.0) GOTO 1146
      DO 146 K=1,NDCT
      IF(WC(K).EQ.X(L)) GOTO 147
  146 CONTINUE
 1146 NDCT=NDCT+1
      WC(NDCT)=X(L)
      IX(L)=NDCT
      GOTO 145
  147 IX(L)=K
      GOTO 145
 1145 IX(L)=-1
  145 CONTINUE
      IF(NAB.EQ.NT) GOTO 30
      NDCT=NDCT+1
      DO 2310 L=1,NT
 2310 IF(IX(L).EQ.-1) IX(L)=NDCT
      GOTO 30

C     MATRIX CONDITIONAL NOMINAL DATA

  573 DO 574 L=1,NS
      NDSBB=0
      N1=N2+1
      N2=N2+NC2
      DO 575 J=N1,N2
      IF(IX(J).EQ.1) GOTO 1575
      IF(NDSBB.EQ.0) GOTO 1577
      DO 576 K=1,NDSBB
      IF(WC(K).EQ.X(J)) GOTO 577
  576 CONTINUE
 1577 NDSBB=NDSBB+1
      WC(NDSBB)=X(J)
      IX(J)=NDSBB
      GOTO 575
  577 IX(J)=K
      GOTO 575
 1575 IX(J)=-1
  575 CONTINUE
      IF(NAD(L,1).EQ.NC2) GOTO 574
      NDSBB=NDSBB+1
      DO 2311 J=N1,N2
 2311 IF(IX(J).EQ.-1) IX(J)=NDSBB
  574 NDSR(L,1)=NDSBB
      GOTO 30

C     ROW CONDITIONAL NOMINAL DATA

 2017 DO 2018 L=1,NS
      DO 2019 I=1,NB
      NDSRR=0
      N1=N2+1
      N2=N2+NB
      DO 2020 J=N1,N2
      IF(IX(J).EQ.1) GOTO 2021
      IF(NDSRR.EQ.0) GOTO 2022
      DO 2023 K=1,NDSRR
      IF(WC(K).EQ.X(J)) GOTO 2024
 2023 CONTINUE
 2022 NDSRR=NDSRR+1
      WC(NDSRR)=X(J)
      IX(J)=NDSRR
      GOTO 2020
 2024 IX(J)=K
      GOTO 2020
 2021 IX(J)=-1
 2020 CONTINUE
      IF(NAD(L,I).EQ.NB) GOTO 2019
      NDSRR=NDSRR+1
      DO 2025 J=N1,N2
 2025 IF(IX(J).EQ.-1) IX(J)=NDSRR
 2019 NDSR(L,I)=NDSRR
 2018 CONTINUE
   30 CONTINUE

C-----------------------------------------------------------------------

C     PERFORM INITIAL OPTIMAL SCALING
C     -------------------------------

C---THE FOLLOWING TWO STATEMENTS ARE REPLACED BY THE NEXT STATEMENT
C     IF(NWE.EQ.0.AND.NDTYP.NE.1) GOTO 8888
C     IF(NWT.EQ.0.AND.NWW.EQ.0)RETURN
      IF(NWE.EQ.0.AND.NDTYP.EQ.1)RETURN
 8888 CALL DISTP(W,CFL,X,WC,WD,IX,IY,IZ,XX,NDSR,WS,NAD)
      IF(ISTFOR.EQ.2)RETURN
      IF(NWE.EQ.0.AND.STRSS.LT.0.5)RETURN
      IF(INITX.GT.1.OR.INITW.GT.1.OR.INITWS.GT.1.OR.INITXC.GT.1)RETURN
C-PC--------------------
      WRITE(*,8889)STRSS
C-----------------------
 8889 FORMAT(10X,'0',2X,F15.5)
      N=0
      DO 442 L=1,NS
      IF(NSIM.GT.1) GOTO 2026
      DO 444 I=1,NB
  444 WA(I,I,L)=0.0
      DO 443 I=2,NB
      I1=I-1
      DO 443 J=1,I1
      N=N+1
      WA(I,J,L)=X(N)
      WA(J,I,L)=X(N)
  443 CONTINUE
      GOTO 442
 2026 DO 2027 I=1,NB
      DO 2027 J=1,NB
      N=N+1
 2027 WA(I,J,L)=X(N)
      DO 2028 I=1,NB
      DO 2028 J=1,I
      WA(I,J,L)=(WA(I,J,L)+WA(J,I,L))*.5
      WA(J,I,L)=WA(I,J,L)
 2028 CONTINUE
  442 CONTINUE
      GOTO 446

      END


C*********************************************************************

      SUBROUTINE STEP3(IX,IY,IZ,X,G,WA,WB,WC,U11,U12,U22,UB1,UB2,
C     ================

     1 XN,XX,CFL,W,TR,FK,WS,NDSR,NAD,IJKL,IRET)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/08/82

C  THIS IS THE MAJOR COMPUTATIONAL ROUTINE. IT CONTROLS THE FLOW
C  OF THE ITERATIVE MINIMIZATION PROCESS, PERFORMS REGRESSIONS
C  TO CALCULATE WEIGHTS, AND OUTPUTS THE RESULTS

C    WA,WB & WC REFER TO THE SAME ARRAY IN THE CALLING ROUTINE
C    THESE ARRAYS ARE PASSED IN AS DIFFERENT PARAMETERS TO PERMIT
C    REFERENCING THE ARRAYS USING DIFFERENT NUMBERS OF DIMENSIONS

C-----------------------------------------------------------------------

      DOUBLE PRECISION T,TT,A,B,C
C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 7/8/82
      INTEGER IX(NT),IY(1),IZ(1)
      DIMENSION X(NT),G(NBS,1),WA(NB,NB,NS),WB(NBS,NS),WC(1),
     1 U11(NDX,1),U12(NDX,1),U22(NDX,1),UB1(1),UB2(1),XN(1)
      DIMENSION XX(NB,NB),CFL(NB,1),W(NS,1),TR(1),FK(NDXP,1)
      DIMENSION WS(NB,1),NDSR(NS,NB),NAD(NS,NB)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
C----- THE FOLLOWING THREE STATEMENTS ADDED 7/7/82
      DOUBLE PRECISION CUT,STMIN
      INTEGER DEBUG,ICNSTR,NOULB
      COMMON /PRMBLK/CUT,STMIN,DEBUG,ICNSTR,NOULB
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
C    +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
     +       /INICON/INITX,INITW,INITWS,INITXC


C     ITERATION BEGINS HERE
C     ---------------------

      IRET=0
      DO 2 LL=1,MAXIT

C     COMPUTE SQUARED DISTANCES,DISPARITIES AND STRESS

      CALL DISTP(W,CFL,X,WC,G,IX,IY,IZ,XX,NDSR,WS,NAD)

      DIF=STRSO-STRSS
C     TEST FOR EXIT
      IF(DIF.LT.0.0) GOTO 145
C----- FOLLOWING STATEMENT ADDED 7/7/82 (SAS 10/23/81)
      IF(STRSS.LT.STMIN) GOTO 141
      IF(DIF.LT.EPSI) GOTO 143
      IF(LL.GE.MAXIT) GOTO 3
      STRSO=STRSS

C     SET-UP FOR THE NEXT ITERATION

      A=1./(1.0-STRSS*STRSS)
      IF(NSIM.GT.1) GOTO 2029
C   NSIM <= 1 IMPLIES NC2=NB*(NB-1)/2 THEREFORE
      N=0
      DO 43 L=1,NS
      WA(1,1,L)=0.0
      DO 44 I=2,NB
      I1=I-1
      WA(I,I,L)=0.0
      DO 44 J=1,I1
      N=N+1
      WA(I,J,L)=X(N)*A
   44 WA(J,I,L)=WA(I,J,L)
   43 CONTINUE
      GOTO 2030
 2029 DO 2031 L=1,NBNBNS
 2031 WC(L)=X(L)*A
 2030 CONTINUE
      IF(NWE.EQ.0) GOTO 4017
      L=0
      DO 47 I=1,NB
      DO 47 J=1,NB
      L=L+1
      DO 47 K=1,NDIM
   47 G(L ,K)=(CFL(I,K)-CFL(J,K))**2
      IF(NWE.EQ.2) GOTO 52
      IF(INITW.EQ.3)GOTO 4533

C     (1) ESTIMATE SUBJECT WEIGHTS

      DO 48 I=1,NDIM
      DO 48 J=1,I
      T=0.0
      L=0
      DO 49 KK=1,NB
      TT=WS(KK,I)*WS(KK,J)
      DO 49 K=1,NB
      L=L+1
   49 T=T+G(L,I)*G(L,J)*TT
      FK(I,J)=T
   48 FK(J,I)=T
      CALL MINV(FK,NDIM,NDXP)
      DO 51 J=1,NDIM
      DO 51 I=1,NS
      T=0.0
      DO 50 K=1,NDIM
      L=0
      DO 50 M=1,NB
      TT=WS(M,K)*FK(J,K)
      DO 50 N=1,NB
      L=L+1
   50 T=T+WB(L,I)*G(L,K)*TT
   51 W(I,J)=T

C     CHECK FOR NEGATIVE SUBJECT WEIGHTS

      IF(NNC.EQ.1) GOTO 4533
      DO 471 K=1,NS
      DO 476 J=1,NDIM
  476 TR(J)=W(K,J)
      DO 472 J=1,NDIM
      IF(TR(J).LT.0.0) GOTO 473
  472 CONTINUE
      GOTO 471
  473 TR(J)=0.0
      IF(NDIM.EQ.1) GOTO 4070
      J1=J+1
      IF(J.EQ.NDIM) J1=1
      DO 482 LX=1,20
      DO 478 I=J1,NDIM
      A=0.0
      B=0.0
      DO 477 I1=1,NBS
      KK=(I1-1)/NB+1
      C=WB(I1,K)
      DO 474 I2=1,NDIM
  474 IF(I2.NE.I)C=C-G(I1,I2)*WS(KK,I2)*TR(I2)
      A=A+C*G(I1,I)*WS(KK,I)
  477 B=B+G(I1,I)**2*WS(KK,I)**2
      TR(I)=A/B
  478 IF(TR(I).LT.0.0) TR(I)=0.0
      DO 479 J=1,NDIM
      IF(ABS(TR(J)-W(K,J)).GT.0.0005) GOTO 483
  479 CONTINUE
      GOTO 4070
  483 DO 481 J=1,NDIM
  481 W(K,J)=TR(J)
      J1=1
  482 CONTINUE
      GOTO 471
 4070 DO 484 J=1,NDIM
  484 W(K,J)=TR(J)
  471 CONTINUE
      DO 355 J=1,NDIM
      DO 356 I=1,NS
      IF(W(I,J).NE.0.0) GOTO 355
  356 CONTINUE
      GOTO 4
  355 CONTINUE
 4533 IF(NWE.EQ.1) GOTO 4017
   52 IF(INITWS.EQ.3)GOTO 4017

C     (2) ESTIMATE STIMULUS WEIGHTS

      DO 4022 K=1,NB
      DO 4018 I=1,NDIM
      DO 4018 J=1,I
      T=0.0
      N=(K-1)*NB
      DO 4019 L=1,NB
      N=N+1
      TT=G(N,I)*G(N,J)
      DO 4019 KK=1,NS
 4019 T=T+TT*W(KK,I)*W(KK,J)
      FK(J,I)=T
 4018 FK(I,J)=T
      CALL MINV(FK,NDIM,NDXP)
      DO 4021 J=1,NDIM
      T=0.0
      DO 4020 I=1,NDIM
      KK=(K-1)*NB
      DO 4020 N=1,NB
      KK=KK+1
      TT=G(KK,I)*FK(I,J)
      DO 4020 L=1,NS
 4020 T=T+TT*WB(KK,L)*W(L,I)
 4021 WS(K,J)=T

C     CHECK FOR NEGATIVE STIMULUS WIEGHTS

      IF(NNC.EQ.1) GOTO 4022
      DO 4026 J=1,NDIM
 4026 TR(J)=WS(K,J)
      DO 4023 J=1,NDIM
      IF(TR(J).LT.0.0) GOTO 4024
 4023 CONTINUE
      GOTO 4022
 4024 TR(J)=0.0
      IF(NDIM.EQ.1) GOTO 4071
      J1=J+1
      IF(J.EQ.NDIM) J1=1
      DO 4025 LX=1,20
      DO 4027 I=J1,NDIM
      A=0.0
      B=0.0
      DO 5029 I1=1,NB
      II1=(K-1)*NB+I1
      DO 5029 I11=1,NS
      C=WB(II1,I11)
      DO 4028 I2=1,NDIM
 4028 IF(I2.NE.I)C=C-G(II1,I2)*TR(I2)*W(I11,I2)
      A=A+C*G(II1,I)*W(I11,I)
 5029 B=B+G(II1,I)**2*W(I11,I)**2
      TR(I)=A/B
 4027 IF(TR(I).LT.0.0) TR(I)=0.0
      DO 4030 J=1,NDIM
      IF(ABS(TR(J)-WS(K,J)).GT.0.0005) GOTO 4031
 4030 CONTINUE
      GOTO 4071
 4031 DO 4033 J=1,NDIM
 4033 WS(K,J)=TR(J)
      J1=1
 4025 CONTINUE
      GOTO 4022
 4071 DO 4032 J=1,NDIM
 4032 WS(K,J)=TR(J)
 4022 CONTINUE
      DO 4034 J=1,NDIM
      DO 4035 I=1,NB
      IF(WS(I,J).NE.0.0) GOTO 4034
 4035 CONTINUE
      GOTO 4
 4034 CONTINUE

C     (3) ESTIMATE CONFIGURATION

 4017 IF(INITX.EQ.3.AND.NSIM.LT.4)GOTO 2
      IF(INITX.EQ.3.AND.INITXC.EQ.3.AND.NSIM.GE.4)GOTO 2
      CALL INNER(CFL,W,X,WB,U11,U12,U22,UB1,UB2,XN,NB,NDIM,NS,NDX,NBS,
     1 WS)
    2 CONTINUE

C----------------------------------------------------------------------
  
C???    3 WRITE(*,9901)
C??? 9901 FORMAT(/' ALSCAL MESSAGE:  ITERATIONS STOPPED BECAUSE ')
C???      WRITE(*,299)MAXIT
C???  299 FORMAT(T19,'THIS IS ITERATION',I4//)
 3    STRSS=STRSO
      GOTO 9999
C----- THE FOLLOWING FOUR STATEMENTS ADDED 7/7/82 (SAS 10/23/81) 
C???  141 WRITE(*,9901)
C???      WRITE(*,142)STMIN
C???  142 FORMAT(T19,'S-STRESS LESS THAN',F9.6//)
 141  GOTO 9999  
C???  143 WRITE(*,9901)
C???      WRITE(*,144)EPSI
C???  144 FORMAT(T19,'SSTRESS IMPROVEMENT LESS THAN',F9.6//)
 143  GOTO 9999
  145 WRITE(*,9902)
 9902 FORMAT(/' ALSCAL WARNING:  ITERATIONS TERMINATED BECAUSE OF NEGATI
     *VE SSTRESS IMPROVEMENT.'/18X,'DATA MAY BE ILL-CONDITIONED OR THERE
     * MAY BE A PROGRAM ERROR.'//)
      GOTO 9999
    4 WRITE(*,9900)
 9900 FORMAT(/' ALSCAL WARNING:    COMPUTATIONS INTERRUPTED.')
      WRITE(*,298)
  298 FORMAT(8X,'A DIMENSION HAS ONLY ZERO WEIGHTS. SOLUTION SKIPPED')
      IF(IJKL.EQ.ND) GOTO 4099
      REWIND NDP
      READ(NDP) X
 2899 READ(NDP) WA
 4099 REWIND NDP
      IRET=1
      RETURN

C-PC-----------
C??? 9999 CALL HITR
C--------------
 9999 RETURN
      END


C***********************************************************************

      SUBROUTINE STEP3A(IX,IY,IZ,X,WC,WD,XX,CFL,W,WS,NDSR,NAD,
C     =================

     1 PHISUB,PHISTI,STRSS1)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 23MAY83

C     SUBROUTINE TO PRINT STRESS AND CORRELATIONS IN DISTANCES
C        NORMALIZES WEIGHTS ACCORDING TO CORRELATION

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 7/8/82
      INTEGER IX(NT),IY(1),IZ(1)
      DIMENSION X(NT),WC(1),WD(1),XX(NB,NB),CFL(NB,1)
      DIMENSION W(NS,1),WS(NB,1),PHISUB(NS),PHISTI(NB)
      DIMENSION NDSR(NS,NB),NAD(NS,NB)
      CHARACTER*72,COPR,LINE*120
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
      COMMON /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /IONUMS/NDP,NDQ,NDR,NDPP

      LINE=' '
      ISTFOR=1
      IF(NSIM.GT.3)ISTFOR=2

C     STRESS IN DISTANCES

      DO 80 I=1,NT
      IF(X(I).LT.0.0) GOTO 4081
      X(I)=SQRT(X(I))
      GOTO 4082
 4081 X(I)=-SQRT(-X(I))
 4082 IF(WC(I).LT.0.0) GOTO 4083
      WC(I)=SQRT(WC(I))
      GOTO 80
 4083 WC(I)=-SQRT(-WC(I))
   80 CONTINUE
C???      WRITE(*,8451)
C??? 8451 FORMAT(////' STRESS AND SQUARED CORRELATION (RSQ) IN DISTANCES'/
C???     */' RSQ VALUES ARE THE PROPORTION OF VARIANCE OF THE SCALED DATA ',
C???     *'(DISPARITIES) IN THE PARTITION'/' (ROW, MATRIX, OR ENTIRE DATA)',
C???     *' WHICH IS ACCOUNTED FOR BY THEIR CORRESPONDING DISTANCES.')
C???      WRITE(*,7483)ISTFOR
C-PC----------------------------------------------------------------
C???      CALL CLEAR
C???      WRITE (*,8452)
C??? 8452 FORMAT(/' STRESS AND SQUARED CORRELATION (RSQ) IN DISTANCES'/)
C-------------------------------------------------------------------

C??? 7483 FORMAT(/' KRUSKAL''S STRESS FORMULA',I2,' IS USED.'/)
      IF(NWC-1)8460,8470,8480

C UNCONDITIONAL STRESS OUTPUT

 8460 CALL MSTRS(NT,X,WC,STRSS1,STRSS2)
      PHI=1.0-STRSS2
      STRSS1=SQRT(STRSS1)
C???      WRITE(*,8461)STRSS1,PHI
C??? 8461 FORMAT(' STRESS = '    ,F5.3,'      RSQ = ',F5.3)
      GOTO 8490

C MATRIX CONDITIONAL STRESS OUTPUT

 8470 STRES1=0.0
      STRES2=0.0
C???      IF(NS.NE.1)WRITE(*,8473)
C??? 8473 FORMAT(4('   MATRIX     STRESS      RSQ '))
      N=1
      DO 8479 L=1,NS
      CALL MSTRS(NC2,X(N),WC(N),STRSS1,STRSS2)
      STRES1=STRES1+STRSS1
      STRES2=STRES2+STRSS2
      PHISUB(L)=1.0-STRSS2
      STRSS1=SQRT(STRSS1)
      LMOD4=MOD(L-1,4)
      LSTART=5+LMOD4*30
      LEND=LSTART+24
C???      IF(NS.NE.1)WRITE(LINE(LSTART:LEND),8471)L,STRSS1,PHISUB(L)
C??? 8471 FORMAT(I3,F13.3,F9.3)
C???      IF(LMOD4.EQ.3)WRITE(*,8475)LINE
C??? 8475 FORMAT(A)
      N=N+NC2
 8479 CONTINUE
C???      IF(LMOD4.NE.3)WRITE(*,8475)LINE(1:29+LMOD4*30)
      LINE=' '
      PHI=1.0-STRES2/NS
      STRES1=SQRT(STRES1/NS)
C???      IF(NS.NE.1)WRITE(*,8472)STRES1,PHI
C???      IF(NS.NE.1)WRITE(*,9472)STRES1,PHI
C??? 8472 FORMAT(/'   OVERALL',F10.3, F9.3/)
C???      IF(NS.EQ.1)WRITE(*,8461)STRES1,PHI
      GOTO 8490

C ROW CONDITIONAL STRESS OUTPUT

 8480 STRES1=0.0
      STRES2=0.0
      DO 8485 J=1,NB
 8485 PHISTI(J)=0.0
      N1=1

C----- THE FOLLOWING TWO STATEMENTS WERE ADDED IN CONJUNCTION
C----- WITH THE ERROR NOTED BELOW. 7/6/82 (SAS CHANGE 12/1/80)

      NBORNR=NB
      IF(NSIM.GT.3)NBORNR=NB-NCOL

      DO 8489 I=1,NS
C???      WRITE(*,8481) I
C??? 8481 FORMAT(' MATRIX',I4/)
C???      IF(NSIM.GT.3)WRITE(*,2000)
C??? 2000 FORMAT(' ROW STIMULI')
C???      WRITE(*,2002)
C??? 2002 FORMAT(4('  STIMULUS    STRESS      RSQ '))
      STR1=0.0
      STR2=0.0
      N=N1
      DO 8488 J=1,NB

C----- THE FOLLOWING STATEMENT WAS MOVED TO CORRECT AND ERROR IN
C----- CALCULATING AVERAGE STRESS AND RSQ INDICES FOR RECTANGULAR
C----- ROW CONDITIONAL DATA. 7/6/82 (SAS CHANGE 12/1/80)
C----- THE ERROR WAS THAT THE PERFECT FIT INDICES FOR COLUMNS
C----- WERE BEING INCLUDED IN THE SUMMATION PROCESS.

      IF(NSIM.GT.3.AND.J.LE.NCOL)GOTO 2001
      CALL MSTRS(NB,X(N), WC(N), STRSS1,STRSS2)
      STR1=STR1+STRSS1
      STR2=STR2+STRSS2
      PHIROW=1.0-STRSS2
      PHISTI(J)=PHISTI(J)+PHIROW
      STRSS1=SQRT(STRSS1)
C     IF(NSIM.GT.3.AND.J.LE.NCOL)GOTO 2001
      K=J
      IF(NSIM.GT.3)K=J-NCOL
      LMOD4=MOD(K-1,4)
      LSTART=5+LMOD4*30
      LEND=LSTART+24
C???      WRITE(LINE(LSTART:LEND),8471)K,STRSS1,PHIROW
C???      IF(LMOD4.EQ.3)WRITE(*,8475)LINE
 2001 N=N+NB
 8488 CONTINUE
C???      IF(LMOD4.NE.3)WRITE(*,8475)LINE(1:29+LMOD4*30)
      LINE=' '
      STRES1=STRES1+STR1
      STRES2=STRES2+STR2

C----- THE FOLLOWING TWO STATEMENTS WERE CHANGED TO COMPUTE
C----- AVERAGE STRESS AND RSQ INDICES CORRECTLY FOR RECTANGULAR
C----- ROW CONDITIONAL DATA.  THE ERROR WAS THAT THE DIVISION
C----- WAS BY THE TOTAL NUMBER OF SIMULI, NOT THE NUMBER OF ROWS.
C----- 7/6/82 (SAS CHANGE 12/1/80)
C     PHISUB(I)=1.0-STR2/NB
C     STR1=SQRT(STR1/NB)

      PHISUB(I)=1.0-STR2/NBORNR
      STR1=SQRT(STR1/NBORNR)
C???      WRITE(*,8472)STR1,PHISUB(I)
      N1=N1+NC2
 8489 CONTINUE
      IF (NS.EQ.1)GOTO 8490
C???      WRITE(*,8491)
C??? 8491 FORMAT(/' AVERAGED OVER MATRICES')
C???      IF(NSIM.LT.4)WRITE(*,8492)
C???      IF(NSIM.GT.3)WRITE(*,8493)
C??? 8492 FORMAT(/'  STIMULUS',6X,'RSQ')
C??? 8493 FORMAT(/'  ROW STIMULUS  RSQ')
      DO 8486 I=1,NB
      PHISTI(I)=PHISTI(I)/NS
      K=I
      IF(NSIM.GT.3.AND.K.LE.NCOL)GOTO 8486
      IF(NSIM.GT.3)K=I-NCOL
C???      WRITE(*,8487)K,PHISTI(I)
 8486 CONTINUE
C??? 8487 FORMAT(T5,I3,F13.3,F9.3)
      NA=NB
      IF(NSIM.GT.3)NA=NB-NCOL
      PHI=1.0-STRES2/(NA*NS)
      STRES1=SQRT(STRES1/(NA*NS))
C???      WRITE(*,9472)STRES1,PHI
C??? 9472 FORMAT(/' OVERALL STRESS =',F6.3,' AND RSQ =',F6.3)

C     NORMALIZE DERIVED SOLUTION

 8490 CALL NORMX(CFL,W,NB,NS,NDIM,NWE)

      ICONFL=0
      IF(NWC.EQ.0)ICONFL=1
      IF(NWE.EQ.1.OR.NWE.EQ.3)CALL NORMW(W,NS,NDIM,PHI,PHISUB,ICONFL)
      ICONFL=1
      IF(NWC.EQ.2)ICONFL=0
      IF (NWE.GE.2) CALL NORMW(WS,NB,NDIM,PHI,PHISTI,ICONFL)
      NAD(1,1)=-NAD(1,1)
      CALL DISTP(W,CFL,X,WC,WD,IX,IY,IZ,XX,NDSR,WS,NAD)
      RETURN
      END


C*********************************************************************

      SUBROUTINE STEP4(X,WC,DISP,CFL,W,TR,WS,
C     ================

     *DIST,DUMMY,PIJP,XK,ROW,IJKL,X1_MY,X2_MY,X3_MY,X4_MY,X5_MY,X6_MY)
C-----LINE ABOVE ADDED TO SUBROUTINE STATEMENT FOR GEMSCAL 8/9/82
C-----IJKL MOVED TO PROPER POSITON IN LINE ABOVE 5MAY83 FWY
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 23MAY83 FWY

C                   THIS SUBROUTINE OUTPUTS THE RESULTS

C-----------------------------------------------------------------------

C-----XX CHANGED TO DISP 8/11/82
      DIMENSION X(NT),WC(NBNBNS),DISP(NB,NB),CFL(NB,1)
      DIMENSION W(NS,1),TR(1),WS(NB,1)
      DIMENSION WL(6)
      DIMENSION X1_MY(NB),X2_MY(NB),X3_MY(NB)
      DIMENSION X4_MY(NB),X5_MY(NB),X6_MY(NB)
      CHARACTER*80,FMT,COPR*72
C-PC------------------------------------------------------------------
      CHARACTER*1 ALFA(52),ALFAL(26),ALFAS(26),ALFA1*52
      EQUIVALENCE(ALFA,ALFA1),(ALFA1(1:26),ALFAL),(ALFA1(27:52),ALFAS)
C---------------------------------------------------------------------

      COMMON /CCOPR/COPR
C-----FOLLOWING THREE STATEMENTS ADDED FOR GEMSCAL 8/9/82
      DIMENSION DIST(NB,NB),DUMMY(NB,NB),PIJP(NB,NB),XK(NB,6),ROW(NB)
      COMMON /PDCOM/IFLPDS,NDIR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP
C-PC--------------------------------------------------------------------
      DATA ALFA1/'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxy'/
C-----------------------------------------------------------------------

C-----FOLLOWING FOUR STATEMENTS ADDED 8/9/82
      IF(NDIR.EQ.0)NDIR=NDIM
C-----FOLLOWING STATEMENT ADDED 19/10/82 FOR GEMSCAL
      DO 80 I=1,NT
      IF(X(I).LT.0.0) GOTO 4081
      X(I)=SQRT(X(I))
      GOTO 4082
 4081 X(I)=-SQRT(-X(I))
 4082 IF(WC(I).LT.0.0) GOTO 4083
      WC(I)=SQRT(WC(I))
      GOTO 80
 4083 WC(I)=-SQRT(-WC(I))
   80 CONTINUE

C      PRINT RESULTS

      IF(NDIM.GT.1)CALL ARNGE(W,CFL,WS,NB,NS,NDIM,TR)
      DO 64 I=1,NB
         X1_MY(I)=CFL(I,1)
         IF(NDIM.GT.1) THEN
            X2_MY(I)=CFL(I,2)
            IF(NDIM.GT.2) THEN
               X3_MY(I)=CFL(I,3)
               IF(NDIM.GT.3) THEN
                  X4_MY(I)=CFL(I,4)
                  IF(NDIM.GT.4) THEN
                     X5_MY(I)=CFL(I,5)
                     IF(NDIM.GT.5) THEN
                        X6_MY(I)=CFL(I,6)
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
 64      KK=MOD(I-1,52)+1
C???  64         WRITE(*,521) I,ALFA(KK),(CFL(I,J),J=1,NDIM)
C---------------------------------------------------------

 9999 RETURN
      END


C***********************************************************************

      SUBROUTINE INIT(W,CFL,CFR,NB,NS,ND,XX,WA,TR,FK,C,XEQ,AM,WS,
C     ===============

     * NDX,NDXP)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 12/14/77
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     INITIAL CONFIGURATION ESTIMATION IN INDSCAL
C     BY SCHONEMANN-DE LEEUW METHOD

C     YOSHIO TAKANE JULY 1974

C----------------------------------------------------------------------

      DIMENSION WA(NB,NB,NS),WS(NB,1),CFL(NB,1),CFR(NB,1)
      DIMENSION C(NDX,NDX,1),XX(NB,NB),TR(1),FK(NDXP,1)
      DIMENSION AM(NS,1),XEQ(NS,1),W(NS,1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION T1,T2,T3,XM,XS


      CALL CJEIG(XX,CFR,CFL,NB,ND+1,FK,WS,TR,1, NDXP)
C  WS IS BEING PASSED TO CJEIG FOR USE AS SCRATCH SPACE


      DO 12 K=1,NS
      DO 12 K1=1,ND
      DO 12 K2=1,K1
      T1=0.0
      T2=TR(K1)*TR(K2)
      T2=1.0D0/DSQRT(T2)
      DO 13 J=1,NB
      T3=0.0
      DO 10 I=1,NB
   10 T3=T3+WA(I,J,K)*CFR(I,K1)
   13 T1=T1+T3*CFR(J,K2)
      T1=T1*T2
      C(K1,K2,K)=T1
   12 C(K2,K1,K)=T1
      DO 50 J=1,ND
      DQ=SQRT(TR(J))
      DO 50 I=1,NB
   50 CFR(I,J)=CFR(I,J)*DQ

C     FIND FEASIBLE WEIGHTS
C          FORMING A MATRIX

      IF(ND.LE.1) GOTO 301
      DO 14 K=1,NS
      T1=0.0
      DO 17 J=1,ND
      DO 17 I=1,ND
   17 T1=T1+C(I,J,K)
   14 TR(K)=T1/NB
      DO 15 I=1,NS
      DO 15 J=1,I
      T1=0.0
      DO 16 K=1,ND
      DO 16 L=1,ND
   16 T1=T1+C(K,L,I)*C(L,K,J)
      AM(I,J)=T1/NB-TR(I)*TR(J)
   15 AM(J,I)=AM(I,J)


      CALL CJEIG(AM,XEQ(1,1),XEQ(1,3),NS,2,FK,W,TR,1,NDXP)


      DO 18 I=1,ND
      DO 18 J=1,I
      T1=0.0
      DO 19 K=1,NS
   19 T1=T1+C(I,J,K)*XEQ(K,1)
      AM(I,J)=T1
   18 AM(J,I)=T1

C     FIND THE TRANSFORMATION MATRIX


      CALL EIGK(AM,TR,ND,NS)


C     FIND THE INTITIAL CONFIGURATION

      DO 20 I=1,NB
      DO 20 J=1,ND
      T1=0.0
      DO 27 K=1,ND
   27 T1=T1+CFR(I,K)*AM(K,J)
   20 CFL(I,J)=T1

C     FIND THE INITIAL WEIGHTS

      DO 22 I=1,ND
      DO 22 J=1,I
      T1=0.0
      DO 23 K=1,NB
   23 T1=T1+CFL(K,I)*CFL(K,J)
      FK(I,J)=T1
   22 FK(J,I)=T1


      CALL MINV(FK,ND,NDXP)


      DO 21 L=1,NS
      DO 24 I=1,ND
      DO 24 J=1,I
      T1=0.0
      DO 251 K2=1,NB
      T2=0.0
      DO 25 K1=1,NB
   25 T2=T2+WA(K1,K2,L)*CFL(K1,I)
  251 T1=T1+T2*CFL(K2,J)
      AM(I,J)=T1
   24 AM(J,I)=T1
      DO 28 I=1,ND
      T1=0.0
      DO 26 K1=1,ND
      DO 26 K2=1,ND
   26 T1=T1+FK(I,K1)*AM(K1,K2)*FK(K2,I)
   28 W(L,I)=T1
   21 CONTINUE

C     NORMALIZATION

      DO 57 J=1,ND
      XM=0.0
      XS=0.0
      DO 58 I=1,NB
      XM=XM+CFL(I,J)
   58 XS=XS+CFL(I,J)**2
      XM=XM/NB
      XS=XS/NB-XM**2
      DO 59 I=1,NS
   59 W(I,J)=W(I,J)*XS
      XS=DSQRT(XS)
      DO 60 I=1,NB
   60 CFL(I,J)=(CFL(I,J)-XM)/XS
   57 CONTINUE
      RETURN


  301 DO 302 L=1,NS
  302 W(L,1)=C(1,1,L)
      RETURN
      END


C***********************************************************************

      SUBROUTINE INSWM(DS,CFL,CFR,WS,XX,TR,CV,CW,FK,ZZ,C,NADCT,
C     ================

     1NB,NDIM,ND,NDX,NDXP,NS)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 02/21/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     INITIAL CONFIGURATIONS FOR THE STIMULUS WEIGHTED MODEL

C----------------------------------------------------------------------

      DIMENSION CFL(NB,1),CFR(NB,1),WS(NB,1),XX(NB,NB)
      DIMENSION DS(NB,NB),TR(1),CV(1),CW(1)
      DIMENSION FK(NDXP,1),ZZ(ND,1),C(NDX,NDX,1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION T,TT

      IF(NADCT.EQ.1) GOTO 403
      T=0.0
      DO 17 J=1,NB
      DO 17 I=1,NB
   17 T=T+DS(I,J)**2
      T=T/((NB*NS)**2)
      DO 57 J=1,NB
      DO 57 I=1,NB
   57 DS(I,J)=DS(I,J)/T
      NADCT=1
  403 CONTINUE
      DO 65 J=1,NDIM
   65 CFR(1,J)=CW(J)
      DO 20 I=1,NB
      DO 21 J=1,NB
      DO 22 K=1,NDIM
   22 XX(J,K)=(CFL(I,K)-CFL(J,K))**2*CFR(1,K)
      IF(NDIM.EQ.1) GOTO 21
      N=NDIM
      DO 23 K=2,NDIM
      NM2=K-1
      DO 23 KK=1,NM2
      N=N+1
   23 XX(J,N)=(CFL(I,KK)-CFL(J,KK))*(CFL(I,K)-CFL(J,K))*
     1 2.0*SQRT(CFR(1,K)*CFR(1,KK))
   21 CONTINUE
      DO 24 K=1,ND
      DO 24 J=1,ND
      ZZ(J,K)=0.0
      DO 24 KK=1,NB
   24 ZZ(J,K)=ZZ(J,K)+XX(KK,J)*XX(KK,K)

      CALL MINV(ZZ,ND,ND)

      DO 25 J=1,ND
      CW(J)=0.0
      DO 25 K=1,NB
   25 CW(J)=CW(J)+XX(K,J)*DS(I,K)
      DO 26 J=1,ND
      CV(J)=0.0
      DO 26 K=1,ND
   26 CV(J)=CV(J)+ZZ(J,K)*CW(K)
      DO 27 J=1,NDIM
   27 C(J,J,I)=CV(J)
      IF(NDIM.EQ.1) GOTO 151
      N=NDIM
      DO 28 K=2,NDIM
      NM2=K-1
      DO 28 KK=1,NM2
      N=N+1
      C(K,KK,I)=CV(N)
   28 C(KK,K,I)=CV(N)
  151 CONTINUE
   20 CONTINUE
      DO 29 K=1,NB
      TR(K)=0.0
      DO 30 J=1,NDIM
      DO 30 I=1,NDIM
   30 TR(K)=TR(K)+C(I,J,K)
   29 TR(K)=TR(K)/NB
      DO 31 J=1,NB
      DO 31 I=1,NB
      T=0.0
      DO 32 K=1,NDIM
      DO 32 L=1,NDIM
   32 T=T+C(K,L,I)*C(L,K,J)
   31 XX(I,J)=T/NB-TR(I)*TR(J)

      IF(NDIM.EQ.1) GOTO 130
      CALL CJEIG(XX,WS,CFR,NB,2,FK,ZZ,TR,1,NDXP)
      GOTO 131
  130 DO 132 K=1,NB
  132 WS(K,1)=TR(K)
  131 CONTINUE

      DO 33 J=1,NDIM
      DO 33 I=1,NDIM
      XX(I,J)=0.0
      DO 133 K=1,NB
  133 XX(I,J)=XX(I,J)+C(I,J,K)*WS(K,1)
   33 CONTINUE

      CALL EIGK(XX,TR,NDIM,NB)

      DO 60 J=1,NDIM
      DO 64 I=1,NB
      T=0.0
      DO 63 K=1,NDIM
   63 T=T+CFL(I,K)*XX(K,J)
   64 CFR(I,J)=T
   60 CONTINUE
      DO 35 J=1,NDIM
      DO 35 I=1,NB
      CFL(I,J)=CFR(I,J)
      T=0.0
      DO 36 KK=1,NDIM
      TT=0.0
      DO 37 K=1,NDIM
   37 TT=TT+XX(K,J)*C(K,KK,I)
   36 T=T+TT*XX(KK,J)
   35 WS(I,J)=T
      RETURN
      END


C**********************************************************************

      SUBROUTINE DISTP(W,CFL,X,WC,WD,IX,IY,IZ,XX,NDSR,WS,NAD)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/08/82

C  THIS ROUTINE COMPUTES DISTANCES AND DISPARITIES (OPTIMALLY SCALED
C  DATA) FOR A GIVEN SET OF COORDINATES AND WEIGHTS BY THE YOUNG,
C  DE LEEUW AND TAKANE METHOD

C----------------------------------------------------------------------

      DOUBLE PRECISION S
C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 7/8/82
      INTEGER IX,IY,IZ
      DIMENSION W(NS,1),CFL(NB,1),X(1),WC(1),WD(NT),IX(NT),IY(1),IZ(1),
     1 XX(1),WS(NB,1),NDSR(NS,NB),NAD(NS,NB)
      DIMENSION R(5,5),ALPH(5)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
C    +       /BLOCK3/FMT
     +       /IONUMS/NDP,NDQ,NDR,NDPP

C ---------------------------------------------------------------------

C     COMPUTE SQUARED DISTANCES
C     -------------------------

      LAST=0
      IF(NAD(1,1).LT.0)LAST=1
      IF(NAD(1,1).LT.0)NAD(1,1)=-NAD(1,1)
      OSEST=BIG
      STRSS=0.0
      LIN=0
      IF(NDEG.EQ.1.AND.NCST.EQ.1)LIN=1
      N=0
      NX=2
      IF(NSIM.GT.1) NX=1
      I1=NB
      DO 37 L=1,NS
      DO 35 I=NX,NB
      IF(NSIM.LE.1)I1=I-1
      DO 35 J=1,I1
      N=N+1
      S=0.0
      IF(NWE-1)25,30,38
   25 DO 26 K=1,NDIM
   26 S=S+(CFL(I,K)-CFL(J,K))**2
      GOTO 34
   30 DO 33 K=1,NDIM
   33 S=S+W(L,K)*(CFL(I,K)-CFL(J,K))**2
      GOTO 34
   38 DO 36 K=1,NDIM
   36 S=S+W(L,K)*WS(I,K)*(CFL(I,K)-CFL(J,K))**2
   34 WC(N)=S
      WD(N)=S
   35 CONTINUE
   37 CONTINUE

C---------------------------------------------------------------------

C     PERFORM OPTIMAL SCALING (COMPUTE DISPARITIES)
C     ---------------------------------------------

      N=1
      IF(NDTYP-3) 134,32,42

C -----  ORDINAL DATA  -----

   32 IF(NWC.EQ.1) GOTO 434
      IF(NWC.EQ.2) GOTO 506

C     UNCONDITIONAL ORDINAL DATA

      IF(NPS.NE.1) CALL PRS(IX,IY,WC,NAB,XX,IZ,1)
      IF(NPS.EQ.1) CALL SES(IX,IY,WD,NAB,IZ,1)
      CALL TRS(WD,X,IX,NAB)
      IF(NAB.EQ.NT)GOTO 150
      NAB1=NAB+1
      DO 151 I=NAB1,NT
      J=IX(I)
      IF(LAST.EQ.1)OSEST=WC(J)
  151 X(J)=OSEST
  150 IF(LAST.EQ.1)RETURN
      CALL LENGTH(WC,X,NT)
      CALL MSTRS(NT,X,WC,STRSS,STRSS2)
      STRSS=SQRT(STRSS)
      RETURN

C     MATRIX CONDITIONAL ORDINAL DATA

  434 DO 435 L=1,NS
      NC3=NAD(L,1)
      IF(NPS.NE.1) CALL PRS(IX(N),IY,WC(N),NC3,XX,IZ,L)
      IF(NPS.EQ.1) CALL SES(IX(N),IY,WD(N),NC3,IZ,L)
      CALL TRS(WD(N),X(N),IX(N),NC3)
      IF(NC3.GE.NC2) GOTO 453
      NI=N+NC3
      NE=N+NC2-1
      DO 152 I=NI,NE
      J=IX(I)+N-1
      IF(LAST.EQ.1)OSEST=WC(J)
  152 X(J)=OSEST
  453 IF(LAST.EQ.1)GOTO 1436
      CALL LENGTH(WC(N),X(N),NC2)
      CALL MSTRS(NC2,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 1436 N=N+NC2
  435 CONTINUE
      STRSS=SQRT(STRSS/NS)
      RETURN

C     ROW CONDITIONAL ORDINAL DATA

  506 N1=0
      DO 507 L=1,NS
      DO 507 I=1,NB
      NC3=NAD(L,I)
      N1=N1+1
      IF(NSIM.GT.3.AND.I.LE.NCOL)GOTO 504
      IF(NPS.NE.1) CALL PRS(IX(N),IY,WC(N),NC3,XX,IZ,N1)
      IF(NPS.EQ.1) CALL SES(IX(N),IY,WD(N),NC3,IZ,N1)
      CALL TRS(WD(N),X(N),IX(N),NC3)
  504 IF(NC3.GE.NB)GOTO 505
      IF(NSIM.GT.3.AND.I.LE.NCOL)X(N+I-1)=0.0
      NI=N+NC3
      NE=N+NB-1
      DO 508 J=NI,NE
      K=IX(J)+N-1
      IF(LAST.EQ.1)OSEST=WC(K)
  508 X(K)=OSEST
  505 IF(LAST.EQ.1)GOTO 1506
      CALL LENGTH(WC(N),X(N),NB)
      CALL MSTRS(NB,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 1506 N=N+NB
  507 CONTINUE
      NA=NB
      IF(NSIM.GT.3)NA=NB-NCOL
      STRSS=SQRT(STRSS/(NA*NS))
      RETURN

C  -----  NOMINAL DATA  -----

   42 IF(NWC.EQ.1) GOTO 436
      IF(NWC.EQ.2) GOTO 509

C     UNCONDITIONAL NOMINAL DATA

      CALL CATSES(WC,X,IX,NT,WD,NDCT)
      IF(NAB.EQ.NT.OR.NPS.NE.1)GOTO 43
      DO 201 I=1,NT
      IF(LAST.EQ.1)OSEST=WC(I)
  201 IF(IX(I).EQ.NDCT) X(I)=OSEST
   43 IF(LAST.EQ.1)RETURN
      CALL LENGTH(WC,X,NT)
      CALL MSTRS(NT,X,WC,STRSS,STRSS2)
      STRSS=SQRT(STRSS)
      RETURN

C     MATRIX CONDITIONAL NOMINAL DATA

  436 DO 437 L=1,NS
      CALL CATSES(WC(N),X(N),IX(N),NC2,WD,NDSR(L,1))
      IF(NAD(L,1).EQ.NC2.OR.NPS.NE.1) GOTO 637
      NE=N+NC2-1
      DO 202 I=N,NE
      IF(LAST.EQ.1)OSEST=WC(I)
  202 IF(IX(I).EQ.NDSR(L,1)) X(I)=OSEST
  637 IF(LAST.EQ.1)GOTO 2436
      CALL LENGTH(WC(N),X(N),NC2)
      CALL MSTRS(NC2,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 2436 N=N+NC2
  437 CONTINUE
      STRSS=SQRT(STRSS/NS)
      RETURN

C     ROW CONDITIONAL NOMINAL DATA

  509 DO 510 L=1,NS
      DO 510 I=1,NB
      IF(NSIM.GT.3.AND.I.LE.NCOL)GOTO 503
      CALL CATSES(WC(N),X(N),IX(N),NB,WD,NDSR(L,I))
  503 IF(NAD(L,I).EQ.NB.OR.NPS.NE.1) GOTO 610
      IF(NSIM.GT.3.AND.I.LE.NCOL)X(N+I-1)=0.0
      NE=N+NB-1
      DO 511 J=N,NE
      IF(LAST.EQ.1)OSEST=WC(J)
  511 IF(IX(J).EQ.NDSR(L,I)) X(J)=OSEST
  610 IF(LAST.EQ.1)GOTO 1510
      CALL LENGTH(WC(N),X(N),NB)
      CALL MSTRS(NB,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 1510 N=N+NB
  510 CONTINUE
      NA=NB
      IF(NSIM.GT.3)NA=NB-NCOL
      STRSS=SQRT(STRSS/(NA*NS))
      RETURN

C     NUMERICAL DATA

  134 REWIND NDP
      READ(NDP)WD
      REWIND NDP
      IF(LIN.EQ.1)GOTO 411
      DO 400 I=1,NT
	WDI=WD(I)
	IF(WDI.NE.BIG) WD(I)=WDI*WDI
  400 CONTINUE
  411 IF(NWC.EQ.1) GOTO 310
      IF(NWC.EQ.2) GOTO 512

C     UNCONDITIONAL NUMERICAL DATA

      CALL POLYF(X,WC,NT,WD,R,ALPH,NDEG,NCST,NAB)
      IF(LIN.EQ.1)CALL LINT(X,WC,NT,WD)
      IF(LAST.NE.1)GOTO 420
      DO 415 I=1,NT
  415 IF(X(I).EQ.BIG)X(I)=WC(I)
      RETURN
  420 CALL LENGTH(WC,X,NT)
      CALL MSTRS(NT,X,WC,STRSS,STRSS2)
      STRSS=SQRT(STRSS)
      RETURN

C     MATRIX CONDITIONAL NUMERICAL DATA

  310 DO 311 L=1,NS
      CALL POLYF(X(N),WC(N),NC2,WD(N),R,ALPH,NDEG,NCST,NAD(L,1))
      IF(LIN.EQ.1)CALL LINT(X(N),WC(N),NC2,WD(N))
      IF(LAST.NE.1)GOTO 300
      DO 305 I=1,NC2
      NI=N+I-1
  305 IF(X(NI).EQ.BIG)X(NI)=WC(NI)
      GOTO 1311
  300 CALL LENGTH(WC(N),X(N),NC2)
      CALL MSTRS(NC2,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 1311 N=N+NC2
  311 CONTINUE
      STRSS=SQRT(STRSS/NS)
      RETURN

C     ROW CONDITIONAL NUMERICAL DATA

  512 DO 513 L=1,NS
      DO 513 I=1,NB
      IF(NSIM.GT.3.AND.I.LE.NCOL)GOTO 514
      IF(NAD(L,I).LE.1)GOTO 514
      CALL POLYF(X(N),WC(N),NB,WD(N),R,ALPH,NDEG,NCST,NAD(L,I))
      IF(LIN.EQ.1)CALL LINT(X(N),WC(N),NB,WD(N))
      GOTO 516
  514 DO 515 J=1,NB
      K=J+N-1
  515 X(K)=WC(K)
  516 IF(LAST.NE.1)GOTO 500
      DO 501 J=1,NB
      NI=N+J-1
  501 IF(X(NI).EQ.BIG)X(NI)=WC(NI)
      GOTO 1513
  500 CALL LENGTH(WC(N),X(N),NB)
      CALL MSTRS(NB,X(N),WC(N),STRSS1,STRSS2)
      STRSS=STRSS+STRSS1
 1513 N=N+NB
  513 CONTINUE
      NA=NB
      IF(NSIM.GT.3)NA=NB-NCOL
      STRSS=SQRT(STRSS/(NA*NS))
      RETURN
      END


C***********************************************************************

      SUBROUTINE PRS(IR,IBK,D,N,W,IZ,L)
C     ==============

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 08/12/82

C     PREPARES TIES FOR CONTINUOUS ORDINAL TRANSFORMATION
C     ARRANGE TIED OBSERVATIONS IN ORDER OF ESTIMATED VALUES

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 8/12/82
      INTEGER IR,IBK,IZ
      DIMENSION IR(1),IBK(1),D(1),W(1),IZ(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

      NB=1
      IF(L.NE.1)NB=IZ(L-1)
      NE=IZ(L)-2
      IF(NE.LT.NB) RETURN
      DO 30 JJ=NB,NE,2
      I=IBK(JJ)
      IF(I.GT.N) RETURN
      J=IBK(JJ+1)
C-----THE FOLLOWING FOUR STATEMENTS REPLACED BY THE NEXT FIVE 8/12/82
C     II=I+J-1
C     DO 10 K=I,II
C  10 W(K)=D(IR(K))
C     CALL SHEL9(W(I),IR(I),J)
      DO 10 I2=1,J
      K=I+I2-1
      W(I2)=D(IR(K))
   10 CONTINUE
      CALL SHEL9(W(1),IR(I),J)
   30 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE SES(IR,IBK,D,N,IZ,L)
C     ==============

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 08/04/82

C     PREPARES TIES FOR DISCRETE ORDINAL TRANSFORMATION
C     TAKE THE MEAN OF TIED OBSERVATIONS

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT WAS CHANGED FROM INTEGER*2 8/4/82
      INTEGER IR,IBK,IZ
      DIMENSION IR(1),IBK(1),D(1),IZ(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

      NB=1
      IF(L.NE.1)NB=IZ(L-1)
   21 NE=IZ(L)-2
      IF(NE.LT.NB) RETURN
      DO 30 II=NB,NE,2
      I=IBK(II)
      IF(I.GT.N) RETURN
      J=IBK(II+1)
      JJ=I+J-1
      S=0.0
      DO 10 K=I,JJ
   10 S=S+D(IR(K))
      S=S/J
      DO 12 K=I,JJ
   12 D(IR(K))=S
   30 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE TRS(DISP,DIST,IVEC,NELE)
C     ==============

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/23/82

C     COMPUTE KRUSKAL'S LEAST SQUARES MONOTONIC TRANSFORMATION

C     THIS ROUTINE DESTROYS THE DISTANCES IN DIST, THUS FORCING
C     DOUBLE STORAGE IN THE CALLING ROUTINE

C    WRITTEN BY ROSTYSLAW JAREMA LEWYCKYJ     MAY  1977

C----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 7/23/82
      INTEGER IVEC
      DIMENSION DIST(1),DISP(1),IVEC(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

C     PLACE DISTANCES IN ASCENDING ORDER IN DIST

      DO 141 I=1,NELE
  141 DIST(I)=DISP(IVEC(I))

C     PERFORM KRUSKAL'S LEAST SQUARES MONOTONIC TRANSFORMATION
C     PLACING RESULTING DISPARITIES IN DISP

      DO 70 I=2,NELE
      IM1=I-1

C     DETERMINE IF ORDER IS CORRECT

      IF(DIST(I).GE.DIST(IM1))GOTO 70
      SUM=DIST(I)
      FN=1.0

C     IF NOT DETERMINE BLOCK SIZE

      DO 10 J=1,IM1
      IMJ=I-J
      SUM=SUM+DIST(IMJ)
      FN=FN+1.0
      DISPT=SUM/FN
      IF(J.EQ.IM1)GOTO 20
      IF(DISPT.GE.DIST(IMJ-1))GOTO 20
   10 CONTINUE

C     SET BLOCK OF DISPARITIES EQUAL TO MEAN DISTANCE IN BLOCK

   20 DO 30 K=IMJ,I
   30 DIST(K)=DISPT
   70 CONTINUE

C     PLACE DISPARITIES IN DISP IN STANDARD ORDER

      DO 40 I=1,NELE
   40 DISP(I)=DIST(I)
      DO 143 I=1,NELE
  143 DIST(IVEC(I))=DISP(I)
      RETURN
      END


C***********************************************************************

      SUBROUTINE LENGTH(DIST,DISP,N)
C     =================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 11/10/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     NORMALIZE LENGTH OF DISPARITIES SO THAT THEY OPTIMIZE
C        NORMALIZED SSTRESS INSTED OF RAW STRESS.
C        ADJUSTMENT IS BASED ONLY ON DISTANCES WHICH CORRESPOND TO
C        ACTIVE DATA.  ADJUSTMENT IS ONLY MADE TO DISPARITIES WHICH
C        CORRESPOND TO ACTIVE DATA.  TYPE OF ADJUSTMENT DEPENDS ON
C        STRESS FORMULA BEING OPTIMIZED.

C-----------------------------------------------------------------------

      DIMENSION DIST(1),DISP(1)
      DOUBLE PRECISION SSQ,SPD,RATIO,DISTV,PROD,DISTM,DISPM
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL

      IF(NSIM.GT.3)GOTO 3

C     ADJUSTMENT FOR STRESS FORMULA 1

      SSQ=0.0
      SPD=0.0
      DO 1 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 1
      SSQ=SSQ+DIST(I)**2
      SPD=SPD+DIST(I)*DISP(I)
    1 CONTINUE
      IF(SPD.EQ.0.0)RETURN
      RATIO=SSQ/SPD
      DO 2 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 2
      DISP(I)=DISP(I)*RATIO
    2 CONTINUE
      RETURN

C     ADJUSTMENT FOR STRESS FORMULA 2

    3 DISTM=0.0
      DISPM=0.0
      NACT=0
      DO 4 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 4
      DISTM=DISTM+DIST(I)
      DISPM=DISPM+DISP(I)
      NACT=NACT+1
    4 CONTINUE
      DISTM=DISTM/NACT
      DISPM=DISPM/NACT
      PROD=0.0
      DISTV=0.0
      DO 5 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 5
      DISTV=DISTV+(DIST(I)-DISTM)**2
      PROD=PROD+(DIST(I)-DISTM)*(DISP(I)-DISPM)
    5 CONTINUE
      IF(PROD.EQ.0.0)RETURN
      RATIO=DISTV/PROD
      DO 6 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 6
      DISP(I)=DISP(I)*RATIO+DISTM*(1.0-RATIO)
    6 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE CATSES(WA,WC,IDY,N,CIV,NCAT)
C     =================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/23/82

C     ROUTINE TO COMPUTE DISCRETE NOMINAL TRANSFORMATION
C        (DISPARITIES ARE MEANS OF ALL DISTANCES IN CATEGORY)

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 7/23/82
      INTEGER IDY(1)
      DIMENSION WA(1),CIV(1),WC(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

      DO 10 I=1,NCAT
      CIV(I)=0.0
   10 WC(I)=0.0

C     AT THIS MOMENT WC IS BEING USED AS SCRATCH SPACE.
C     LATER IT WILL BE ASSIGNED NEW VALUES TO BE RETURNED TO THE
C     CALLING ROUTINE

      DO 11 I=1,N
      ID=IDY(I)
      CIV(ID)=CIV(ID)+WA(I)
   11 WC(ID)=WC(ID)+1.0
      DO 12 I=1,NCAT
      CIV(I)=CIV(I)/WC(I)
   12 CONTINUE

C     NOW WC IS ASSIGNED VALUES TO BE RETURNED TO THE CALLER

      DO 13 I=1,N
   13 WC(I)=CIV(IDY(I))
      RETURN
      END


C***********************************************************************

      SUBROUTINE POLYF(DISP,DIST,N,OBS,R,ALPH,NDEG,NCST,NAB)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/02/82

C     POLYNOMIAL FITTING

C-----------------------------------------------------------------------

      DIMENSION DISP(1),DIST(1),OBS(1)
      DIMENSION R(5,5),XY(5),ALPH(5)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSS2,NB,NS,NDTYP,NPS,NWC,MDEG,NT,NBS,NBNBNS

      NOR=NDEG
      LIN=0
      IF(NOR.EQ.1.AND.NCST.EQ.1) LIN=1
      IF(LIN.EQ.1) NOR=2
      NOR1=NOR+1
      DO 11 I=1,NOR1
      XY(I)=0.0
      DO 11 J=1,I
   11 R(I,J)=0.0
      R(1,1)=NAB
      DO 12 K=1,N
      IF(OBS(K).EQ.BIG) GOTO 12
      G=1.0
      XY(1)=XY(1)+DIST( K)
      L=0
      DO 13 I=2,NOR1
      I1=I-1
      DO 13 J=I1,I
      G=G*OBS(K)
      L=L+1
      R(I,J)=R(I,J)+G
      IF(L.GT.NOR) GOTO 13
      XY(L+1)=XY(L+1)+G*DIST( K)
   13 CONTINUE
   12 CONTINUE
      IF(NOR1.LT.3) GOTO 117
      DO 15 I=3,NOR1
      K=I-2
      DO 15 J=1,K
   15 R(I,J)=R(I-1,J+1)
  117 CONTINUE
      DO 16 I=1,NOR1
      DO 16 J=1,I
   16 R(J,I)=R(I,J)
      IF(NCST.EQ.1) GOTO 20
      DO 21 I=2,NOR1
      I1=I-1
      XY(I1)=XY(I)
      DO 21 J=2,NOR1
      J1=J-1
   21 R(I1,J1)=R(I,J)
      NOR1=NOR
   20 CONTINUE
      CALL MINV(R,NOR1,5)
      DO 17 I=1,NOR1
      ALPH(I)=0.0
      DO 17 J=1,NOR1
   17 ALPH(I)=ALPH(I)+R(I,J)*XY(J)
      IF(LIN.EQ.1) RETURN
      DO 25 K=1,N
      IF(OBS(K).EQ.BIG) GOTO 23
      DISP( K)=ALPH(1)
      IF(NCST.EQ.0) DISP( K)=ALPH(1)*OBS(K)
      IF(NOR1.LT.2) GOTO 18
      G=1.0
      IF(NCST.EQ.0) G=OBS(K)
      DO 19 I=2,NOR1
      G=G*OBS(K)
   19 DISP( K)=DISP( K)+G*ALPH(I)
      GOTO 18
   23 DISP(K)=BIG
   18 CONTINUE

C THE FOLLOWING LINE HAS BEEN ADDED JULY 1982.  BEFORE THIS CHANGE
C ALL NUMERICAL LEVELS OF MEASUREMENT WERE UNCONSTRAINED EXCEPT
C LINEAR INTERVAL, WHICH WAS CONSTRAINED TO YIELD POSITIVE ESTIMATES.
C AFTER THIS CHANGE ALL NUMERICAL TRANSFORMATIONS EXCEPT LINEAR INTERVAL
C USE PIECE-WISE LINEAR REGRESSION.  NEFATIVELY EXTIMATED DISPARITIES
C ARE REPLACED WITH ZERO ESTIMATES. CHANGE MADE BY YOUNG AND SARLE.
C IMPLEMENTED BY BROOKS.

      IF(DISP(K).LT.0.0)DISP(K)=0.0
   25 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE BLOC2(A,LL,IZ,N)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/23/82

C     FIND BLOCKS OF TIES

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 7/23/82
      INTEGER LL,IZ
      DIMENSION A(N),LL(N)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /IONUMS/NDP,NDQ,NDR,NDPP
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL

      II=1
    2 I=II
    3 II=II+1
      IF(II.GT.N) GOTO 7
      IF(A(I).EQ.A(II))GOTO 3
    7 IMI=II-I
      IF(IMI.EQ.1) GOTO 10
      IF(IZ.GE.NDMX-1)GOTO 4
      LL(IZ)=I
      LL(IZ+1)=IMI
      IZ=IZ+2
   10 IF(II.GE.N) RETURN
      GOTO 2
    4 WRITE(*,9900)
      WRITE(*,200)NDMX
 9900 FORMAT(/' ALSCAL FATAL ERROR:    COMPUTATIONS TERMINATED')
  200 FORMAT(8X,'MAXIMUM NUMBER OF TIE-BLOCKS EXCEEDS',I6)
      STOP
      END


C*********************************************************************

      SUBROUTINE INNER(CFL,W,X,WB,U11,U12,U22,UB1,UB2,XN,NB,NDIM,NS,
C     ================

     1 NDX,NBS,WS)
COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 03/22/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C  THIS ROUTINE COMPUTES STIMULUS COORDINATES

C----------------------------------------------------------------------

      DIMENSION CFL(NB,1),W(NS,1),WB(NBS,1),X(1),WS(NB,1)
      DIMENSION U11(NDX,1),U12(NDX,1),U22(NDX,1),UB1(1),UB2(1),XN(1)
      DOUBLE PRECISION T1,T2,T3,T4,TU11,TU12,TU22
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL
     +       /INICON/INITX,INITW,INITWS,INITXC

      NNS=NB*NS
      DO 15 M=1,NB
      IF(INITX.EQ.3.AND.M.GT.NCOL)GOTO 15
      IF(INITXC.EQ.3.AND.M.LE.NCOL)GOTO 15
      DO 13 I=1,NDIM
      T1=WS(M,I)
      DO 13 J=1,NDIM
      T2=T1*WS(M,J)
      TU11=0.0
      TU12=0.0
      TU22=0.0
      DO 14 L=1,NB
      IF(L.EQ.M)GOTO 14
      T3=0.0
      DO 24 K=1,NS
   24 T3=T3+W(K,I)*W(K,J)
      IF(J.GT.I)GOTO 12
      TU22=TU22+T3*T2
      TU11=TU11+T3*T2*CFL(L,I)*CFL(L,J)
   12 TU12=TU12+T3*T2*CFL(L,I)
   14 CONTINUE
      IF(J.GT.I)GOTO 13
      U11(I,J)=TU11
      U22(I,J)=TU22
   13 U12(I,J)=TU12
      JK=0
      DO 39 K=1,NS
      LJ=(M-1)*NB
      DO 39 J=1,NB
      JK=JK+1
      LJ=LJ+1
      IF(J.EQ.M)GOTO 153
      T1=WB(LJ,K)
      DO 40 I=1,NDIM
   40 T1=T1-W(K,I)*CFL(J,I)**2*WS(M,I)
      X(JK)=T1
      GOTO 39
  153 X(JK)=0.0
   39 CONTINUE
      DO 16 I=1,NDIM
      T1=0.0
      T2=0.0
      JK=0
      DO 18 K=1,NS
      T3=0.0
      T4=0.0
      DO 17 J=1,NB
      JK=JK+1
      T4=T4+X(JK)
   17 T3=T3+CFL(J,I)*X(JK)
      T1=T1+T3*W(K,I)
   18 T2=T2+T4*W(K,I)
      UB1(I)=T1*WS(M,I)
      UB2(I)=T2*WS(M,I)
   16 CONTINUE
      IF(NSIM.LT.2) GOTO 50
      IF(NWE.LT.2 ) GOTO 80
      DO 84 I=1,NDIM
      DO 84 J=1,NDIM
      TU11=0.0
      TU22=0.0
      TU12=0.0
      DO 83 L=1,NB
      IF(L.EQ.M)GOTO 83
      T1=0.0
      DO 82 K=1,NS
   82 T1=T1+W(K,I)*W(K,J)
      IF(J.GT.I)GOTO 81
      TU22=TU22+T1*(WS(L,I)*WS(L,J))
      TU11=TU11+T1*(WS(L,I)*WS(L,J))*CFL(L,I)*CFL(L,J)
   81 TU12=TU12+T1*(WS(L,I)*WS(L,J))*CFL(L,I)
   83 CONTINUE
      IF(J.GT.I)GOTO 84
      U11(I,J)=(U11(I,J)+TU11)*.5D0
      U22(I,J)=(U22(I,J)+TU22)*.5D0
   84 U12(I,J)=(U12(I,J)+TU12)*.5D0
   80 CONTINUE
      JK=NNS
      DO 59 K=1,NS
      LJ=M-NB
      DO 59 J=1,NB
      JK=JK+1
      LJ=LJ+NB
      IF(M.EQ.J)GOTO 253
  254 T1=WB(LJ,K)
      DO 60 I=1,NDIM
   60 T1=T1-W(K,I)*CFL(J,I)**2*WS(J,I)
      X(JK)=T1
      GOTO 59
  253 X(JK)=0.0
   59 CONTINUE
      DO 66 I=1,NDIM
      T1=0.0
      T2=0.0
      JK=NNS
      DO 67 K=1,NS
      T3=0.0
      T4=0.0
      DO 65 J=1,NB
      JK=JK+1
      T3=T3+X(JK)*WS(J,I)*CFL(J,I)
   65 T4=T4+X(JK)*WS(J,I)
      T1=T1+W(K,I)*T3
   67 T2=T2+W(K,I)*T4
      UB1(I)=(UB1(I)+T1)*.5D0
   66 UB2(I)=(UB2(I)+T2)*.5D0
   50 CONTINUE
      DO 68 I=1,NDIM
      DO 86 J=1,NDIM
      IF(J.GT.I)GOTO 86
      U11(I,J)=4.0*U11(I,J)
      U11(J,I)=U11(I,J)
      U22(J,I)=U22(I,J)
   86 U12(I,J)=-2.*U12(I,J)
   68 UB1(I)=-2.0*UB1(I)

      CALL COEF(U11,U12,U22,UB1,UB2,NDIM,CFL,M,XN,NB,NDX)

   15 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE COEF(U11,U12,U22,UB1,UB2,NDIM,CFL,M,XN,NB,NDX)
C     ===============

C COPYRIGHT DECEMBER 1977 BY FORREST W. YOUNG AND ROSTYLWW
C FINAL CHANGE 12/14/77
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     SUBROUTINE TO OBTAIN LEAST SSTRESS COORDINATES
C     BY SOLVING A SYSTEM OF CUBIC EQUATIONS FOR EACH COORDINATE
C     (TAKANE, YOUNG & DE LEEUW, PSYCHOMETRIKA, 1977)

C-----------------------------------------------------------------------

      DIMENSION U11(NDX,1),U12(NDX,1),U22(NDX,1),CFL(NB,1),
     *UB1(1),UB2(1),XN(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

      DO 20 LL=1,30
      DO 9 I=1,NDIM
    9 XN(I)=CFL(M,I)
      DO 10 L=1,NDIM
      A=2.0*U22(L,L)
      P=3.0*U12(L,L)/A
      Q=U11(L,L)-2.0*UB2(L)
      R=-UB1(L)
      DO 11 I=1,NDIM
      IF(I.EQ.L)GOTO 11
      Q=Q+2.0*CFL(M,I)*(U12(I,L)+U22(L,I)*CFL(M,I))
      R=R+CFL(M,I)*(U11(L,I)+U12(L,I)*CFL(M,I))
   11 CONTINUE
      Q=Q/A
      R=R/A
      CALL SCUBE(P,Q,R,CFL(M,L))
   10 CONTINUE
      DO 12 L=1,NDIM
      IF(ABS(CFL(M,L)-XN(L)).GT.0.0001) GOTO 20
   12 CONTINUE
      RETURN
   20 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE SCUBE(P,Q,R,CUBE)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 12/14/77
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     SOLUTION TO A CUBIC EQUATION

C-----------------------------------------------------------------------

      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DATA R3,R27/0.33333333,3.703704E-2/
      F(X)=(((X*0.25+P*R3)*X+Q*0.5)*X+R)*X

      A=(3.0*Q-P**2)*R3
      B=(P**2+P**2-9.0*Q)*P*R27+R
      D=B**2*0.25+A**3*R27
      IF(D.LT.0.0)GOTO 1
      E=SQRT(D)-.5*B
      CUBE=SIGN(ABS(E)**R3,E)-P*R3
      E= -(E+B)
      CUBE=SIGN(ABS(E)**R3,E)+CUBE
      RETURN
    1 T3=ACOS(-B/(2.0*SQRT(-A**3*R27)))*R3
      E=2.0*SQRT(-A*R3)
      C1=COS(T3)
      C2=COS(T3+4.1887902)
      C3=COS(T3+2.0943951)
      CUBE=AMIN1(C1,C2,C3)*E-P*R3
      CMAX=AMAX1(C1,C2,C3)*E-P*R3
      IF(F(CMAX).LT.F(CUBE))CUBE=CMAX
      RETURN
      END


C***********************************************************************

      SUBROUTINE MSTRS(N,DISP,DIST,STRSS1,RESID)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 11/08/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C ROUTINE TO CALCULATE THE STRESS OF TWO VECTORS
C AND TO ESTIMATE MISSING DATA

C     STRESS 1 OR 2 IS CALCULATED (NORMALIZED BY DISPARITIES, WHICH IS
C        THE SAME AS BY DISTANCES).  THE RESIDUAL VARIANCE (1-RSQ)
C        IS CALCULATED FOR NORMALIZATION OF DERIVED WEIGHTS.

C     ESTIMATES OF MISSING DATA ARE THEIR DISTANCES

C-----------------------------------------------------------------------

      DIMENSION DISP(1),DIST(1)
      DOUBLE PRECISION SV4,SV4SQ,SSQ,SXY,S4,S4SQ
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK1/NC,ND,BIG,NC2,NNC,NSC,
     + EPSI,NDIM,NDX,NDXS,NDXP,MAXIT,NADCT,NDCT,STRSO,
     + STRSS,STRSSB,NB,NS,NDTYP,NPS,NWC,NDEG,NT,NBS,NBNBNS
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL

      NACT=0
      SXY=0.0
      S4=0.0
      S4SQ=0.0
      SV4=0.0
      SV4SQ=0.0
      SSQ = 0
      DO 110 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 100
      S4=S4+DIST(I)
      SV4=SV4+DISP(I)
      SV4SQ=SV4SQ+DISP(I)*DISP(I)
      S4SQ=S4SQ+DIST(I)*DIST(I)
      SXY=SXY+DIST(I)*DISP(I)
      SSQ=SSQ+(DISP(I)-DIST(I))**2
      NACT=NACT+1
      GOTO 110
  100 DISP(I)=DIST(I)
  110 CONTINUE
      IF(NACT.LE.1)GOTO 120
      IF(NSIM.LT.4)GOTO 118
      SVAR=0.0
      DISPM=SV4/NACT
      DO 117 I=1,N
      IF(DISP(I).EQ.BIG)GOTO 117
      SVAR=SVAR+(DISP(I)-DISPM)**2
  117 CONTINUE
      STRSS1=SSQ/SVAR
      GOTO 119
  118 STRSS1=SSQ/SV4SQ
  119 RESID=1.0-(SXY-S4*SV4/NACT)**2/
     *((S4SQ-S4*S4/NACT)*(SV4SQ-SV4*SV4/NACT))
      RETURN
  120 STRSS1=0.0
      RESID=0.0
      RETURN
      END


C***********************************************************************

      SUBROUTINE NORMX(X,W,NB,NS,ND,NWE)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 02/12/80
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     NORMALIZE CONFIGURATION FOR PRINTING AT END OF JOB

C-----------------------------------------------------------------------

      DIMENSION X(NB,1),W(NS,1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION SUMSQ,SUM

      SUMSQ=0.0
      DO 4 J=1,ND

C     CENTER CONFIGURATION AT CENTROID

      SUM=0.0
      DO 1 I=1,NB
    1 SUM=SUM+X(I,J)
      SUM=SUM/NB
      DO 2 I=1,NB
      X(I,J)=X(I,J)-SUM
    2 SUMSQ=SUMSQ+X(I,J)**2
      IF (NWE.EQ.0) GOTO 4

C     FOR ALL OF THE WEIGHTED MODELS NORMALIZE SO THAT THE LENGTH OF
C     EACH STIMULUS DIMENSION IS UNITY

      SUMSQ=DSQRT(NB/SUMSQ)
      DO 3 I=1,NB
    3 X(I,J)=SUMSQ*X(I,J)
      SUMSQ=1.0/(SUMSQ**2)
      DO 6 I=1,NS
    6 W(I,J)=SUMSQ*W(I,J)
      SUMSQ=0.0
    4 CONTINUE
      IF(NWE.GT.0) RETURN

C     FOR THE UNWEIGHTED MODEL NORMALIZE LENGTH OF THE DIMENSIONS SO
C     THAT THEIR AVERAGE LENGTH IS UNITY.

      SUMSQ=DSQRT(NB*ND/SUMSQ)
      DO 5 J=1,ND
      DO 5 I=1,NB
    5 X(I,J)=SUMSQ*X(I,J)
      RETURN
      END


C***********************************************************************

      SUBROUTINE NORMW(W,N,ND,PHI,PHIROW,ICONFL)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 12/14/77
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     ROUTINE TO NORMALIZE WEIGHT MATRICES FOR OUTPUT

C-----------------------------------------------------------------------

      DIMENSION W(N,1),PHIROW(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION SUMSQ

C     NORMALIZE EITHER STIMULUS OR SUBJECT WEIGHTS SO THAT THEIR
C     LENGTH CORRESPONDS TO THE PROPORTION  OF VARIANCE IN THE
C     OPTIMALLY SCALED DATA ACCOUNTED FOR BY THE MODEL.

      SUMSQ=0.0
      DO 3 I=1,N
      DO 1 J=1,ND
    1 SUMSQ=SUMSQ+W(I,J)**2
      IF(ICONFL.EQ.1) GOTO 3

C     COME HERE WHEN CONDITIONALITY IS SUCH THAT EACH ROW OF WEIGHTS
C     IS NORMALIZED SEPARATELY.

      SUMSQ=DSQRT(PHIROW(I)/SUMSQ)
      DO 2 J=1,ND
    2 W(I,J)=SUMSQ*W(I,J)
      SUMSQ=0.0
    3 CONTINUE
      IF(ICONFL.EQ.0) RETURN

C     COME HERE WHEN CONDITIONALITY IS SUCH THAT ROWS ARE
C     NORMALIZED JOINTLY.

      SUMSQ=DSQRT(N*PHI/SUMSQ)
      DO 4 J=1,ND
      DO 4 I=1,N
    4 W(I,J)=SUMSQ*W(I,J)
      RETURN
      END


C***********************************************************************

      SUBROUTINE ARNGE(W,X,WS,NB,NS,NDIM,TR)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/23/82

C     SUBROUTINE TO REARRANGE WEIGHT AND STIMULUS MATRICES

C-----------------------------------------------------------------------

      DOUBLE PRECISION T
C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 7/23/82
      INTEGER NN
      DIMENSION W(NS,1),X(NB,1),TR(1),NN(6),WS(NB,1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
     +       /BLOCK2/NCST,NSIM,NWE,NDMX,NAB,NCOL

      NN(1)=1
      DO 11 J=1,NDIM
      T=0.0D0
      IF(NWE/2*2.EQ.NWE) GOTO 20
      DO 10 I=1,NS
   10 T=T+W(I,J)
      GOTO 11
   20 IF(NWE.LT.2) GOTO 30
      DO 40 I=1,NB
   40 T=T+WS(I,J)
      GOTO 11
   30 DO 22 I=1,NB
   22 T=T+X(I,J)**2
   11 TR(J)=(-T)

C     SORT TR INTO ASCENDING ORDER BY INSERTION SORT

      DO 55 J=2,NDIM
      T=TR(J)
      K=J-1
   50 IF(T.GE.TR(K))GOTO 51
      TR(K+1)=TR(K)
      NN(K+1)=NN(K)
      K=K-1
      IF(K.GE.1)GOTO 50
   51 TR(K+1)=T
      NN(K+1)=J
   55 CONTINUE

      DO 12 J=1,NDIM
      K=NN(J)
      IF(J.EQ.K) GOTO 12
      DO 13 I=1,NB
      WK=X(I,J)
      X(I,J)=X(I,K)
   13 X(I,K)=WK
      IF(NWE/2*2.EQ.NWE) GOTO 25
      DO 14 I=1,NS
      WK=W(I,J)
      W(I,J)=W(I,K)
   14 W(I,K)=WK
   25 IF(NWE.LT.2) GOTO 31
      DO 32 I=1,NB
      WK=WS(I,J)
      WS(I,J)=WS(I,K)
   32 WS(I,K)=WK
   31 CONTINUE
      JJ=J+1
      DO 15 I=JJ,NDIM
      IF(NN(I).EQ.J)GOTO 16
   15 CONTINUE
      GOTO 12
   16 NN(I)=K
   12 CONTINUE
      RETURN
      END



      SUBROUTINE EIGK(A,VALUE,N,NA)
C     ===============

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 02/21/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     EIGEN-ROUTINE BY KAISER (JK METHOD)

C    REF: "THE EIGENVALUES OF A REAL SYMMETRIC MATRIX"
C  BY H.F. KAISER .  IN  THE BRITISH COMPUTER JOURNAL
C  VOL 15 NO 3 PAGES 271-273.

C  THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A REAL
C  SYMMETRIC MATRIX. THE ABSOLUTE VALUES OF THE EIGENVALUES
C  ORDERED IN DESCENDING ORDER OF SIZE ARE RETURNED IN VALUE.
C  NORMALIZED EIGENVECTORS ARE RETURNED IN CORRESPONDING
C  COLUMNS OF A

C--------------------------------------------------------------------

      DIMENSION A(NA,1),VALUE(N)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION Q,HALFP

      IF(N-2)1,2,3
    3 Q=0.0
      DO 101 J=1,N
      DO 101 I=1,N
  101 Q=Q+A(I,J)*A(I,J)
      EPS=0.000001*Q/N
      NLESS1=(N-1)
      NN=(N-1)*N/2
      NCOUNT=NN
  116 DO 102 J=1,NLESS1
      JPLUS1=J+1
      DO 102 K=JPLUS1,N
      HALFP=0.0
      Q=0.0
      DO 103 I=1,N
      HALFP=HALFP+A(I,J)*A(I,K)
  103 Q=Q+(A(I,J)+A(I,K))*(A(I,J)-A(I,K))
      ABSP=DABS(HALFP+HALFP)
      IF(ABSP.LT.EPS.AND.Q.GE.0.0)GOTO 106
      ABSQ=DABS(Q)
      IF(ABSP.GT.ABSQ)GOTO 108
      TAN=ABSP/ABSQ
      COS=1.0/SQRT(1.0+TAN*TAN)
      SIN=TAN*COS
      GOTO 109
  108 CTN=ABSQ/ABSP
      SIN=1.0/SQRT(1.0+CTN*CTN)
      COS=CTN*SIN
  109 COS=SQRT((1.0+COS)/2.0)
      SIN=SIN/(COS+COS)
      IF(Q.GE.0.0)GOTO 111
      TEMP=COS
      COS=SIN
      SIN=TEMP
  111 IF(HALFP.LT.0.0)SIN=-SIN
      DO 114 I=1,N
      TEMP=A(I,J)
      A(I,J)=TEMP*COS+A(I,K)*SIN
  114 A(I,K)=-TEMP*SIN+A(I,K)*COS
      NCOUNT=NN
      GOTO 102
  106 NCOUNT=NCOUNT-1
      IF(NCOUNT.LE.0)GOTO 115
  102 CONTINUE
      GOTO 116
  115 DO 117 J=1,N
      VALUE(J)=0.0
      DO 118 I=1,N
  118 VALUE(J)=VALUE(J)+A(I,J)*A(I,J)
  117 VALUE(J)=SQRT(VALUE(J))
      DO 119 J=1,N
      DO 119 I=1,N
  119 A(I,J)=A(I,J)/VALUE(J)
      RETURN
    1 VALUE(1)=A(1,1)
      A(1,1)=1.0
      RETURN
    2 AD=A(1,1)+A(2,2)
      D=SQRT(AD*AD-4.0*(A(1,1)*A(2,2)-A(1,2)*A(2,1)))
      VALUE(1)=(AD+D)*0.5
      VALUE(2)= AD-VALUE(1)
      VL1=1.0/SQRT(A(1,2)**2+(A(1,1)-VALUE(1))**2)
      A(2,1)=(VALUE(1)-A(1,1))*VL1
      VL2=1.0/SQRT(A(1,2)**2+(A(1,1)-VALUE(2))**2)
      A(2,2)=(VALUE(2)-A(1,1))*VL2
      A(1,1)=A(1,2)*VL1
      A(1,2)=A(1,2)*VL2
      RETURN
      END


C***********************************************************************

      SUBROUTINE CJEIG(A,U,V,N,ND,B,W,ALAM,NFT,NB)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 12/14/77
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     EIGENVALUES AND VECTORS OF REAL SYMMETRIC MATRIX
C     BY CLINT AND JENNINGS' METHOD
C        WRITTEN BY YOSHIO TAKANE, 1974

C-----------------------------------------------------------------------

      DIMENSION A(N ,1),U(N ,1),V(N ,1),B(NB,1),W(N ,1),ALAM(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR
      DOUBLE PRECISION T1

      ND1=ND-1
      IF(NFT.NE.1) GOTO 11
      DO 36 J=1,ND
      DO 35 I=1,N
   35 U(I,J)=0.0
   36 U(J,J)=1.0

   11 DO 1  LL=1,30
      IF(LL.GT.30) RETURN
      DO 10 I=1,N
      DO 10 J=1,ND
      T1=0.0
      DO 9 K=1,N
    9 T1=T1+A(I,K)*U(K,J)
   10 V(I,J)=T1
      DO 13 I=1,ND
      DO 13 J=1,ND
      T1=0.0
      DO 12 K=1,N
   12 T1=T1+U(K,I)*V(K,J)
   13 B(I,J)=T1

      CALL EIGK(B,ALAM,ND,NB)

      DO 14 I=1,N
      DO 14 J=1,ND
      T1=0.0
      DO 41 K=1,ND
   41 T1=T1+V(I,K)*B(K,J)
   14 W(I,J)=T1
      DO 15 I=1,ND
      DO 15 J=1,ND
      T1=0.0
      DO 51 K=1,N
   51 T1=T1+W(K,I)*W(K,J)
   15 B(I,J)=T1

C     CHOLESKY FACTORIZATION

      DO 16 I=1,ND
      IF(B(I,I).LE.0.0) GOTO 18
      B(I,I)=SQRT(B(I,I))
      IF(I.GE.ND)GOTO 18
      J1=I+1
      DO 16 J=J1,ND
      B(I,J)=B(I,J)/B(I,I)
      DO 16 K=J1,J
   16 B(K,J)=B(K,J)-B(I,J)*B(I,K)
   18 CONTINUE

C     INVERSE OF CHOLESKY FACTOR

      DO 60 I=1,ND
      B(I,I)=1.0/B(I,I)
      IF(I.LE.1)GOTO 60
      J2=I-1
      DO 62 J=1,J2
      T1=0.0
      DO 63 K=J,J2
   63 T1=T1-B(K,I)*B(K,J)
      B(I,J)=T1*B(I,I)
   62 CONTINUE
   60 CONTINUE

      DO 30 J=1,ND
      DO 30 I=1,N
   30 V(I,J)=W(I,J)
      DO 19 I=1,N
      DO 19 J=1,ND
      T1=0.0
      DO 24 K=1,J
   24 T1=T1+V(I,K)*B(J,K)
   19 W(I,J)=T1

C     TEST OF CONVERGENCE

      DO 20 J=1,ND1
      DO 20 I=1,N
      IF(ABS(W(I,J)-U(I,J)).GT.1.0E-5)GOTO 21
   20 CONTINUE
      DO 22 J=1,ND
      DO 22 I=1,N
   22 U(I,J)=W(I,J)
      RETURN
   21 DO 23 J=1,ND
      DO 23 I=1,N
   23 U(I,J)=W(I,J)
    1 CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE MINV(A,N,NA)
C     ===============

C FINAL CHANGE 11/16/78
C EQUIVALENT TO ALSCAL82 VERSION (7/2/82)

C     OBTAIN THE INVERSE OF A SYMMETRIC MATRIX

C-----------------------------------------------------------------------

      DIMENSION A(1),L(30),M(30)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

C  IF N IS 1 OR 2 INVERT DIRECTLY

      IF(N-2)200,300,400
  200 A(1)=1.0/A(1)
      RETURN
  300 C=1.0/(A(1)*A(NA+2)-A(NA+1)*A(2))
      W=A(1)
      A(1)=A(NA+2)*C
      A(NA+2)=W*C
      A(NA+1)=-A(NA+1)*C
      A(2)=-A(2)*C
      RETURN

C  PACK THE MATRIX

  400 K=N
      I1=NA+1
      I2=N*NA
      DO 10 I3=I1,I2,NA
      I4=I3+N-1
      DO 10 I=I3,I4
      K=K+1
   10 A(K)=A(I)

C  SEARCH FOR LARGEST ELEMENT

      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
      IF(ABS(BIGA).LT.ABS(A(IJ)))GOTO 20
      BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE

C  INTERCHANGE ROWS

      J=L(K)
      IF(J.LE.K)GOTO 35
      KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI) =HOLD

C  INTERCHANGE COLUMNS

   35 I=M(K)
      IF(I.LE.K)GOTO 45
      JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI) =HOLD

C  DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
C  CONTAINED IN BIGA)

   45 IF(ABS(BIGA).GT.1.0E-20)GOTO 48
  900 PRINT 901
  901 FORMAT(/' ALSCAL FATAL ERROR:  INVERSE OF A SINGULAR MATRIX ATTEMPT
     *ED.')

C  CALL ERRTRA

      STOP
   48 DO 55 I=1,N
      IF(I.EQ.K)GOTO 55
      IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE

C  REDUCE MATRIX

      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I.EQ.K)GOTO 65
      IF(J.EQ.K)GOTO 65
      KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE

C  DIVIDE ROW BY PIVOT

      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J.NE.K)A(KJ)=A(KJ)/BIGA
   75 CONTINUE

C  PRODUCT OF PIVOTS
C  REPLACE PIVOT BY RECIPROCAL

      A(KK)=1.0/BIGA
   80 CONTINUE

C  FINAL ROW AND COLUMN INTERCHANGE

      K=N
  100 K=(K-1)
      IF(K.LE.0)GOTO 150
  105 I=L(K)
      IF(I.LE.K)GOTO 120
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI) =HOLD
  120 J=M(K)
      IF(J.LE.K)GOTO 100
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI) =HOLD
      GOTO 100

C  UNPACK THE MATRIX

  150 I3=N*(NA+1)+1
      I4=N*(N+1)+1
      DO 160 J=2,N
      I3=I3-NA
      I4=I4-N
      DO 160 I=1,N
  160 A(I3-I)=A(I4-I)
      RETURN
      END



C***********************************************************************

      SUBROUTINE SHEL9(A,C,NITEM)
C     ================

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/23/82

C     SUBROUTINE SHEL9   RICARDO DOBSON
C     A IS THE KEY VECTOR, C IS THE VECTOR TO SORT ON
C     NITEM IS THE NUMBER OF ITEMS IN THE TWO VECTORS
C     FIRST VECTOR ARGUMENT MUST BE REAL, SECOND VECTOR ARG
C     MUST BE INTEGER
C     SORT WILL BE IN ASCENDING ORDER

C-----------------------------------------------------------------------

C-----THE FOLLOWING STATEMENT CHANGED FROM INTEGER*2 7/23/82
      INTEGER C(1),KK
      REAL A(1)
      CHARACTER*72,COPR
      COMMON /CCOPR/COPR

      M=NITEM
20    M=M/2
      IF (M) 30,40,30
30    K=NITEM - M
      J=1
41    I=J
49    L=I+M
      IF(A(L)-A(I))50,60,60
50    B = A(I)
      A(I)= A(L)
      A(L)= B
      KK=C(I)
      C(I)=C(L)
      C(L)=KK
      I=I-M
      IF (I-1)60,49,49
60    J=J+1
      IF(J-K)41,41,20
40    CONTINUE
      RETURN
      END


C***********************************************************************

      SUBROUTINE LINT(DISP,DIST,N,OBS)
C     ===============

COPYRIGHT, 1977, FORREST W. YOUNG, YOSHIO TAKANE & ROSTYSLAW J LEWYCKYJ
C FINAL CHANGE 07/14/82

C     ROUTINE TO PERFORM CONSTRAINED LEAST SQUARES LINEAR REGRESSION.
C        IT FINDS THE BEST LINEAR TRANSFORMATION FOR QUANTITATIVE
C        (INTERVAL LEVEL) DATA UNDER THE CONSTRAINT THAT ALL
C        DISPARITIES (PREDICTIONS) BE NON-NEGATIVE.
C        THE NEWTON-RAPHSON METHOD IS USED.

C   ORIGINAL VERSION BY YOSHIO TAKANE AND FORREST YOUNG
C   REWRITTEN BY ROSTYSLAW LEWYCKYJ   MARCH 1977
C   MODIFIED ACCORDING TO SUGGESTIONS BY YOSHIO TAKANE JUNE 1977
C   REWRITTEN IN DOUBLE PRECISION BY YOSHIO TAKANE SEPTEMBER 1978
C   ERRORS CORRECTED BY YOSHIO TAKANE AND FORREST YOUNG, JANUARY 1981
C   ADAPTED FOR STAND ALONE VERSION OF ALSCAL, JULY 1982

C-----------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL DISP(1),DIST(1),OBS(1)
      REAL PRTMSG
      CHARACTER*72,COPR
C-----THE FOLLOWING STATEMENT IS REPLACED BY A DATA STATEMENT
C     COMMON /SLPMSG/ PRTMSG
C-----THE FOLLOWING STATEMENT WAS PUT IN TO TRY TO FIX UP A LINT PROBLEM
C-----   AND IS NOW TAKEN OUT BECAUSE THE FIX DID NOT WORK
C    +       /ITNUMR/ ITER
      COMMON /CCOPR/COPR
     +       /IONUMS/NDP,NDQ,NDR,NDPP
C-----THE FOLLOWING THREE STATEMENTS ADDED 7/14/82
      DOUBLE PRECISION CUT,STMIN
      INTEGER DEBUG,ICNSTR,NOULB
      COMMON /PRMBLK/CUT,STMIN,DEBUG,ICNSTR,NOULB
      DATA PRTMSG/0.0/
      DATA BIG/9.0E20/
C---THE FOLLOWING STATEMENT IS SAS ONLY
C     IOPT19=IOPT(19)
      FMIN=OBS(1)
      DO 10 K=1,N
   10 IF(FMIN.GT.OBS(K))FMIN=OBS(K)
      DH=0.0
      ODH=0.0
      O1=0.0
      O2=0.0
      O3=0.0
      O4=0.0
      DC=0.0
      D0=0.0
      D1=0.0
      D2=0.0

C COMPUTE INITIAL ESTIMATES BY LINEAR REGRESSION

      DO 22 K=1,N
      IF(OBS(K).EQ.BIG) GOTO 22
      DI=SQRT(DIST(K))
      DH=DH+DI
      OB1=OBS(K)
      ODH=ODH+OB1*DI
      DI=DIST(K)
      O1=O1+OB1
      DC=DC+1.0
      D0=D0+DI
      D1=D1+OB1*DI
      OB=OB1*OB1
      O2=O2+OB
      D2=D2+OB*DI
      OB=OB*OB1
      O3=O3+OB
      O4=O4+OB*OB1
   22 CONTINUE
      AO=(ODH*DC-DH*O1)/(O2*DC-O1*O1)
      BO=(DH-AO*O1)/DC
      FO=AO*AO*AO*AO*O4+6.0*AO*AO*BO*BO*O2+DC*BO*BO*BO*BO-2.0*AO*AO*D2-4
     1.0*AO*BO*D1-2.0*BO*BO*D0+4.0*AO*AO*AO*BO*O3+4.0*AO*BO*BO*BO*O1

C PERFORM UNCONSTRAINED NEWTON-RAPHSON ITERATIONS

    4 LL=0
    1 LL=LL+1
      IF(LL.GT.50)GOTO 2
      G1=((O4*AO+3.0*BO*O3)*AO+3.0*BO*BO*O2-D2)*AO+BO*(BO*BO*O1-D1)
      G2=((DC*BO+3.0*AO*O1)*BO+3.0*AO*AO*O2-D0)*BO+AO*(AO*AO*O3-D1)
      H11=D2-3.0*AO*AO*O4-6.0*AO*BO*O3-3.0*BO*BO*O2
      H22=D0-3.0*AO*AO*O2-6.0*AO*BO*O1-3.0*BO*BO*DC
      H12=D1-3.0*AO*AO*O3-6.0*AO*BO*O2-3.0*BO*BO*O1
      DET=H11*H22-H12*H12
      DET=1.0/DET
      H12=-H12*DET
      P=H11
      H11=H22*DET
      H22=P*DET
      S1=H11*G1+H12*G2
      S2=H12*G1+H22*G2
      STEP=1.0
      LLL=0
    3 LLL=LLL+1
      IF(LLL.GT.20) GOTO 2
      AN=AO+STEP*S1
      BN=BO+STEP*S2
      FN=AN*AN*AN*AN*O4+6.0*AN*AN*BN*BN*O2+DC*BN*BN*BN*BN-2.0*AN*AN*D2-4
     1.0*AN*BN*D1-2.0*BN*BN*D0+4.0*AN*AN*AN*BN*O3+4.0*AN*BN*BN*BN*O1
      IF(FN.LT.FO) GOTO 40
      IF(FN.EQ.FO) GOTO 2
      STEP=STEP*0.5
      GOTO 3
   40 GN=H11*G1*G1+2.0*H12*G1*G2+H22*G2*G2
      GN=DSQRT(-GN)
      FO=FN
      IF(GN.LE.0.0001) GOTO 2
      AO=AN
      BO=BN
      GOTO 1

C TEST FOR NEGATIVE SLOPE.

    2 IF(AN.LT.0.0) GOTO 50

C TEST FOR NEGATIVE ESTIMATES

  222 IF(AN*FMIN+BN.LT.0.0) GOTO 30

C COME HERE WHEN SLOPE AND ESTIMATES ARE POSITIVE

      DO 12 K=1,N
      IF(OBS(K).EQ.BIG)GOTO 13
      DISP(K)=(AN*OBS(K)+BN)*(AN*OBS(K)+BN)
      GOTO 12
   13 DISP(K)=BIG
   12 CONTINUE
C---THE FOLLOWING STATEMENT IS REPLACED BY THE NEXT
C     IF(IOPT19.EQ.1)WRITE(*,4738)AN,BN
      IF(DEBUG.GT.0)WRITE(*,4738)AN,BN
 4738 FORMAT(' LINT: POSITIVE EST+++SLOPE=',E15.5,'INTERCEPT=',E15.5)
      RETURN

C COME HERE WHEN THE SLOPE IS NEGATIVE.  SET THE INTERCEPT TO ZERO
C (RATIO LEVEL), PRINT WARNING, AND SOLVE FOR TRANSFORMATION.

   50 FMIN=0.0
      IF(PRTMSG.EQ.1)GOTO 30
      PRTMSG=1
      WRITE(*,9998)
      WRITE(*,9995)

C COME HERE WHEN NEGATIVE ESTIMATES DO EXIST AND IMPOSE CONSTRAINT

   30 T1=0.0
      T2=0.0
      DO 36 K=1,N
      IF(OBS(K).EQ.BIG)GOTO 36
      S=(OBS(K)-FMIN)**2
      T1=T1+S*DIST(K)
      T2=T2+S**2
   36 CONTINUE
      A=T1/T2
      DO 37 K=1,N
      IF(OBS(K).EQ.BIG)GOTO 38
      DISP(K)=A*(OBS(K)-FMIN)**2
      GOTO 37
   38 DISP(K)=BIG
   37 CONTINUE
C---THE FOLLOWING STATEMENT REPLACED BY THE NEXT 7/14/82
C     IF(IOPT19.EQ.1)WRITE(*,4748)A,FMIN
      IF(DEBUG.GT.0)WRITE(*,4748)A,FMIN
 4748 FORMAT(' LINT: NEGATIVE EST---SLOPE=',E15.5,'INTERCEPT=',E15.5)
      RETURN
 9995 FORMAT(T18,'THE MEASUREMENT LEVEL HAS BEEN CHANGED TO RATIO ON'/
     *,T18,'THIS ITERATION AND WILL BE CHANGED BACK TO INTERVAL ON'/
     *,T18,'THE NEXT.  THIS MAY RESULT IN DIVERGENCE.  IF SO,'/
     *,T18,'REANALYZE YOUR DATA WITH A DIFFERENT MEASUREMENT LEVEL.'/
     *,T18,'IN ALL CASES YOU SHOULD CHECK THE OPTIONS OF THE'/
     *,T18,'ANALYSIS, PARTICULARLY THE SIMILAR OPTION.')
 9998 FORMAT(/' ALSCAL WARNING: A LINEAR TRANSFORMATION OF YOUR DATA HAS
     + NEGATIVE SLOPE.')
      END



