      SUBROUTINE DECLAR(SUBR,DUMMY,RSUM)
 
************************************************************************
*                                                                      *
* Function    : Checks for enough memory.                              *
* Version     : 1.0                                                    *
* Date        : August 1996                                            *
* Author      : Patrick Groenen                                        *
* Calls       : SUBR                                                   *
* Used by     : Adapted for Prefmap3                                   *
*                                                                      *
************************************************************************
 
C-----------------------------------------------------------------------
C Parameters.
C
      INTEGER          RSUM
      EXTERNAL         SUBR
      REAL             DUMMY
C-----------------------------------------------------------------------
C Local variables.
C
      INTEGER          RMAXME
      PARAMETER        (RMAXME = 50000)
      REAL             A(RMAXME)
C-----------------------------------------------------------------------
C Program start.
C

      DO 10 I=1,RMAXME
        A(I) = 0
   10 CONTINUE
      IF (RSUM .LE. RMAXME) THEN
         CALL SUBR(A,RSUM)
      ELSE
         PRINT *,'  '
         PRINT *,' Not enough memory available for this task.'
      END IF
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     P R E F M A P - 3         VERSION 1.0          FEBRUARY 1985     C
C                                                                      C
C     E X T E R N A L   U N F O L D I N G                              C
C                                                                      C
C                                                                      C
C                                                                      C
C  INPUT DESCRIPTION AND OPTIONS                                       C
C                                                                      C
C                                                                      C
C    -CARD 1: TITLE (20A4)                                             C
C                                                                      C
C    -CARD 2: NROW,NCOL,IOPT,(IDIF(K),K=1,4) (7I5)                     C
C                                                                      C
C             NROW   = NUMBER OF ROW POINTS;                           C
C             NCOL   = NUMBER OF COLUMN POINTS                         C
C             IOPT   = 0: ALL ROWS SAME OPTIONS (DEFAULT)              C
C                    = 1: ALL ROWS DIFFERENT OPTIONS                   C
C                    = 2: INDICATED ROWS SAME OPTIONS, AS SPECIFIED IN C
C             IDIF   = STARTING POINTS FOR DIFFERENT OPTIONS           C
C                                                                      C
C    -CARD 3: NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT         C
C                                            (9I5,F10.8)               C
C                                                                      C
C             NDIM   = NUMBER OF DIMENSIONS                            C
C             NANA   = NUMBER OF SEPARATE ANALYSES                     C
C             IPRE   = 0: NO PRELIMINARY TRANSFORMATION (DEFAULT)      C
C                    = 1: CONFIGURATION WILL BE WEIGHTED               C
C                    = 2: CONFIGURATION WILL BE ROTATED AND WEIGHTED   C
C             ISTA   = 0: STANDARDIZATION EXTERNAL DATA (DEFAULT)      C
C                    = 1: CENTERING ONLY                               C
C                    = 2: NORMALIZATION ONLY                           C
C                    = 3: NONE OF THE ABOVE                            C
C             MOD4   = 1: VECTOR MODEL WILL BE APPLIED                 C
C             MOD3   = 1: UNFOLDING MODEL WILL BE APPLIED              C
C             MOD2   = 1: WEIGHTED UNFOLDING MODEL WILL BE APPLIED     C
C             MOD1   = 1: GENERAL UNFOLDING MODEL WILL BE APPLIED      C
C             MAX    = MAXIMUM NR OF NON-METRIC ITERATIONS (DEFAULT=50)C
C             CRIT   = CONVERGENCE CRITERION (DEFAULT=0.00001)         C
C                                                                      C
C    -CARD 4: IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,ISTAT,IWRI1,IWRI2 (9I5) C
C                                                                      C
C             IPRI0  = 0: NO PRINT OF INPUT                            C
C                    = 1: PRINT CONFIGURATION                          C
C                    = 2: PRINT EXTERNAL DATA                          C
C                    = 3: PRINT BOTH                                   C
C             IPRI1  = 0: NO RESULTS ACROSS OPTIONS                    C
C                    = 1: PRINT COMPLETE RESULTS ACROSS OPTIONS        C
C                    = 2: PRINT FIT ACROSS OPTIONS                     C
C                    = 3: PRINT FIT, CRITERION AND PREDICTED VALUES    C
C             IPRI2  = 0: NO RESULTS ACROSS ANALYSES                   C
C                    = 1: PRINT RESULTS ACROSS ANALYSES                C
C             ISHE   = 0: NO SCATTER PLOTS                             C
C                    = N: FOR N ROWS PLOT OF CRITERION                 C
C                         AND DISTANCES VERSUS THE DATA                C
C             IPLO   = 0: NO PLOT OF FINAL CONFIGURATION               C
C                    = 1: PLOT FIRST DIMENSION                         C
C                    = 2: PLOT FIRST TWO DIMENSIONS                    C
C                    = K: PLOT (PAIRWISE) K DIMENSIONS                 C
C             IHIS   = 0: NO PRINT OF HISTORY                          C
C                    = 1: PRINT HISTORY OF ITERATION                   C
C             IST    = 0: NO STATISTICS                                C
C                    = 1: PRINT STATISTICS FOR THE METRIC MODELS       C
C             IWRI1  = 0: NO STORAGE                                   C
C                    = 1: STORE COORDINATES                            C
C                    = 2: STORE 1 + WEIGHTS                            C
C                    = 3: STORE 2 + ROTATION MATRICES                  C
C             IWRI2  = 0: NO STORAGE                                   C
C                    = 1: PREDICETD VALUES                             C
C                    = 2: STORE 1 + CRITERION                          C
C                                                                      C
C    -CARD 5: INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4) (9I5)      C
C                                                                      C
C             INCON  = INPUT UNIT CONFIGURATION                        C
C             INDATA = INPUT UNIT EXTERNAL DATA                        C
C             INOPT  = INPUT UNIT OPTIONS                              C
C             ISCR1  = UNIT NUMBER SCRATCH FILE 1                      C
C             ISCR2  = UNIT NUMBER SCRATCH FILE 2                      C
C             LUN(1) = OUTPUT UNIT NUMBER COORDINATES                  C
C             LUN(2) = OUTPUT UNIT NUMBER WEIGHTS                      C
C             LUN(3) = OUTPUT UNIT NUMBER ROTATION MATRICES            C
C             LUN(4) = OUTPUT UNIT NUMBER PREDICTED VALUES             C
C                                                                      C
C                                                                      C
C    -CARD 6: FORMAT OF TARGET CONFIGURATION (20A4)                    C
C                                                                      C
C    -FOLLOWING CARDS: TARGET DATA                                     C
C                                                                      C
C    -FOLLOWING CARD : FORMAT OF EXTERNAL DATA (20A4)                  C
C                                                                      C
C    -FOLLOWING CARDS: EXTERNAL DATA                                   C
C                                                                      C
C     CARDS 1 - 5 AND THE TWO FORMAT CARDS ARE READ FROM INPUT MEDIUM  C
C     INPARA (IN THIS VERSION FIXED TO BE 5).                          C
C                                                                      C
C  TO ESTIMATE THE TOTAL NUMBER OF WORDS NEEDED FOR THE ARRAY AREA,    C
C  USE THE FORMULA NWORDS=NCOL*(2*NPRT+NDIM+5)+NPRT+2*NPRM+NCR1*NCR1   C
C                     +NDIM*NANA*(2*NPL1+NPL2+NDIM+2)+NDIM+NW2+NW3+NW4 C
C                     +NPR1**2+NPR2**2+NPR3**2+NPR4**2+NANA            C
C  IN THE STATIC ALLOCATION VERSION, NWORDS IS FIXED AT 25100.         C
C                                                                      C
C     AUTHORS : JACQUELINE MEULMAN & WILLEM HEISER                     C
C               BELL LABORATORIES                                      C
C               600 MOUNTAIN AVENUE - MURRAY HILL - NEW JERSEY 07974   C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      EXTERNAL MAIN2
      COMMON/CDECLA/NDIM,MAX,CRIT,IPRI1,IPRI2,IPLO,IHIS,ISHE,IWRI1,IST,
     X        IWRI2,IPRE,ISTA,NV,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2,
     X        NCR1,NW2,NW3,NW4,ND,NANA,IPRI0,MOD1,MOD2,MOD3,MOD4
      COMMON/CFORM /FORM1(20)
      COMMON/CHEAD /TITLE(20)
      COMMON/CPARAN/NROW,NCOL,IOPT,IDIF(4)
      COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4)
C
C   READER AND PRINTER SPECIFICATION; CHANGE HERE AND IN SUBROUTINE
C   DECLAR IF OTHER CONVENTIONS ARE USED ON YOUR MACHINE
C
      INPARA=5
      IWRITE=6
C
C
C
C  READ AND WRITE PARAMETER CARDS
C
      READ(INPARA,500) TITLE
      READ(INPARA,501) NROW,NCOL,IOPT,(IDIF(K),K=1,4)
      READ(INPARA,501) NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT
      READ(INPARA,501) IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,IST,IWRI1,IWRI2
      READ(INPARA,501) INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4)
      WRITE(IWRITE,600) TITLE
      WRITE(IWRITE,601)
      WRITE(IWRITE,602) NROW,NCOL,IOPT,(IDIF(K),K=1,4),
     1                 NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT,
     2                 IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,IST,IWRI1,IWRI2,
     3                  INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4)
C
C  SET DEFAULT VALUES AND EXTRA PARAMETERS
C
      IF (IOPT.GT.2) IOPT=0
      IDIF(1)=1
      IF (IOPT.EQ.1.AND.IDIF(2).EQ.0) IDIF(2)=NROW
      IF (IOPT.EQ.1.AND.IDIF(3).EQ.0) IDIF(3)=NROW
      IF (IOPT.EQ.1.AND.IDIF(4).EQ.0) IDIF(4)=NROW
      IF (IOPT.EQ.2.AND.IDIF(2).LE.IDIF(1)) IDIF(2)=NROW+1
      IF (IOPT.EQ.2.AND.IDIF(3).LE.IDIF(2)) IDIF(3)=NROW+1
      IF (IOPT.EQ.2.AND.IDIF(4).LE.IDIF(3)) IDIF(4)=NROW+1
      IF (INDATA.LE.0) INDATA=INPARA
      IF (INCON.LE.0) INCON=INPARA
      IF (INOPT.LE.0) INOPT=INPARA
      IF (LUN(1).LE.0.AND.IWRI1.GT.0) LUN(1)=IWRITE
      IF (LUN(2).LE.0.AND.IWRI1.GT.1) LUN(2)=IWRITE
      IF (LUN(3).LE.0.AND.IWRI1.GT.2) LUN(3)=IWRITE
      IF (LUN(4).LE.0.AND.IWRI2.GT.0) LUN(4)=IWRITE
      IF (NDIM.LE.0) NDIM=2
      IF (NANA.LE.0.OR.NANA.GT.4) NANA=4
      IF (MAX.LE.0) MAX=50
      IF (CRIT.LT.1.0E-10) CRIT=1.0E-5
      IF (IPRE.GT.2) IPRE=0
      IF (MOD4.NE.1) MOD4=0
      IF (MOD3.NE.1) MOD3=0
      IF (MOD2.NE.1) MOD2=0
      IF (MOD1.NE.1) MOD1=0
C
C  GIVE SUMMARY OF OPTIONS
C
      WRITE(IWRITE,603) NROW,NCOL,IOPT
      IN=1
      IR=NROW+1
      IF (IDIF(2).NE.IR) IN=IN+1
      IF (IDIF(3).NE.IR) IN=IN+1
      IF (IDIF(4).NE.IR) IN=IN+1
      IF (IOPT.EQ.2) WRITE(IWRITE,1603) (IDIF(K),K=1,IN)
      WRITE(IWRITE,2603)
      IF (NROW.EQ.0.OR.NCOL.EQ.0) WRITE(IWRITE,700)
      IF (NROW.EQ.0.OR.NCOL.EQ.0) STOP
      WRITE(IWRITE,604) NDIM,NANA,IPRE,ISTA
      WRITE(IWRITE,1604)
      WRITE(IWRITE,605) MOD4,MOD3,MOD2,MOD1,MAX,CRIT
      WRITE(IWRITE,606) IPRI0,IPRI1,IPRI2
      WRITE(IWRITE,1606) ISHE,IPLO
      WRITE(IWRITE,607) IHIS,IST,IWRI1,IWRI2
      WRITE(IWRITE,608) INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4)
C
      IF (ISCR1.EQ.0.OR.ISCR2.EQ.0) WRITE(IWRITE,700)
      IF (ISCR1.EQ.0.OR.ISCR2.EQ.0) STOP
C
C     COMPUTE THE NUMBER OF PREDICTOR VARIABLES FOR THE VARIOUS
C     MODELS. IF A MODEL WILL NOT BE APPLIED, THE NUMBER OF COLUMNS
C     OF ITS PREDICTOR MATRIX WILL BE DIMINISHED TO 1, UNLESS A PRE-
C     LIMINARY TRANSFORMATION IS ASKED FOR (NPR1 AND NPR4)
C     WHEN THE NUMBER OF COLUMNS IS NOT LARGE ENOUGH TO ESTIMATE
C     THE PARAMETERS OF A MODEL, THE MODEL WILL NOT BE APPLIED.
C     WHEN THE NUMBER OF COLUMNS IS NOT LARGE ENOUGH THE PRELIMI-
C     NARY TRANSFORMATION WILL NOT BE PERFORMED.
C     WHEN THE NUMBER OF DIMENSIONS EQUALS 1, THE GENERAL MODEL AND
C     THE PRELIMINARY ROTATION WILL NOT BE APPLIED
C
      NPRE=IPRE
      KPT1=NDIM+2
      KPT2=(NDIM*(NDIM+1))/2+1
      IF (KPT1.GT.NCOL) NPRE=0
      IF (KPT2.GT.NCOL.AND.IPRE.EQ.2) NPRE=0
      IF (NDIM.EQ.1) MOD1=0
      IF (NDIM.EQ.1.AND.IPRE.EQ.2) NPRE=0
      IF (NPRE.NE.IPRE) WRITE (IWRITE,701)
      IPRE=NPRE
C
      NPR1=(NDIM*(NDIM+3))/2+1
      IF (NPR1.GT.NCOL) MOD1=0
      KPR1=MAX0(MOD1,IPRE)
      IF (KPR1.GT.0) KPR1=1
      NPR1=KPR1*(NPR1-1)+1
      NPR2=2*NDIM+1
      IF (NPR2.GT.NCOL) MOD2=0
      NPR2=MOD2*(NPR2-1)+1
      NPR3=NDIM+2
      IF (NPR3.GT.NCOL) MOD3=0
      NPR3=MOD3*(NPR3-1)+1
      NPR4=NDIM+1
      IF (NPR4.GT.NCOL) MOD4=0
      KPR4=MAX0(MOD4,IPRE)
      IF (KPR4.GT.0) KPR4=1
      NPR4=KPR4*(NPR4-1)+1
      NPRT=NPR1+NPR2+NPR3+NPR4
      NPRM=MAX0(NPR1,NPR2,NPR3,NPR4)
      N=MAX0(IPLO,IPRI2)
      IF (N.NE.0) N=1
      NW2=MAX0((NROW+NCOL)*N,NCOL+NCOL,NDIM*NDIM)
      NW3=NW2
      NW4=NW2
      ND =NDIM*NANA
C
C  WHEN PLOTTING OF IDEAL-POINTS AND/OR VECTORS IS NOT DESIRED,
C  THE DIMENSIONS OF TWEIGH,TCOOR AND TROTAT WILL BE DIMINISHED
C
      N=MAX0(IPLO,IPRI2)
      IF (N.NE.0) N=1
      NPL1=(NROW-1)*N+1
      NPL2=(NROW-1)*NDIM*MOD1*N+NDIM
C
C  WHEN PREPROCESSING IS NOT DESIRED,
C  THE DIMENSION OF ARRAY CROSS WILL BE DIMINISHED
C
      N=0
      IF (IPRE.NE.0) N=1
      NCR1=(NCOL-1)*N+1
C
C  ALLOCATE STORAGE AND START COMPUTATIONS
C
C
C
C     CON    =    NCOL*NDIM
C     PR1    =    NCOL*NPR1
C     PR2    =    NCOL*NPR2
C     PR3    =    NCOL*NPR3
C     PR4    =    NCOL*NPR4
C     SIVA1  =    NPR1
C     SIVA2  =    NPR2
C     SIVA3  =    NPR3
C     SIVA4  =    NPR4
C     LSV1   =    NCOL*NPR1
C     LSV2   =    NCOL*NPR2
C     LSV3   =    NCOL*NPR3
C     LSV4   =    NCOL*NPR4
C     RSV1   =    NPR1*NPR1
C     RSV2   =    NPR2*NPR2
C     RSV3   =    NPR3*NPR3
C     RSV4   =    NPR4*NPR4
C     TCOOR  =    NPL1*NDIM*NANA
C     TWEIGH =    NPL1*NDIM*NANA
C     TROTAT =    NPL2*NDIM*NANA
C     COEF   =    NPRM
C     COOR   =    NDIM*NANA
C     WEIGH  =    NDIM*NANA
C     ROTAT  =    NDIM*NANA*NDIM
C     CROSS  =    NCR1*NCR1
C     WSIGN  =    NDIM
C     DATA   =    NCOL
C     DISP   =    NCOL
C     DSTAR  =    NCOL
C     IORD   =    NCOL
C     LBK    =    NCOL
C     WORK1  =    NPRM
C     WORK2  =    NW2
C     WORK3  =    NW3
C     WORK4  =    NW4
C     OPTION =    NANA
C     DANA   =    NCOL
C
C
      NWORDS=NCOL*(2*NPRT+NDIM+6)+NPRT+2*NPRM+NCR1*NCR1
     X      +NDIM*NANA*(2*NPL1+NPL2+NDIM+2)+NDIM+NW2+NW3+NW4
     X      +NPR1**2+NPR2**2+NPR3**2+NPR4**2+NANA
      CALL DECLAR(MAIN2,A,NWORDS)
      STOP
C
C  FORMAT SPECIFICATIONS
C
 500  FORMAT(20A4)
 501  FORMAT(9I5,F10.8)
 600  FORMAT(1H1//,5X,17HP R E F M A P - 3/
     1 1H ,47X,                31H  *** EXTERNAL UNFOLDING ***   ,25X,
     1 18HJACQUELINE MEULMAN /
     2 5X,15H VERSION - 1.0 ,28X,31H  ***         &          ***   ,25X,
     3 13HWILLEM HEISER/
     3 5X,14H FEBRUARY 1985, 29X,31H  *** PROPERTY FITTING   ***   ,25X,
     3 18HJ. DOUGLAS CARROLL/
     4 1H ,103X,17HBELL LABORATORIES/
     5 1H ,103X,15HMURRAY HILL, NJ////
     6 1H ,4X,13H1 JOB TITLE: ,20A4//)
 601  FORMAT(1H ,4X,24HECHO OF PARAMETER CARDS://)
 602  FORMAT(1H ,4X,2H2 ,7I5/
     1       1H ,4X,2H3 ,9I5,F10.8/
     2       1H ,4X,2H4 ,9I5/
     3       1H ,4X,2H5 ,9I5///)
 603  FORMAT(1H ,4X,22H2 DATA SPECIFICATIONS://
     1       1H ,8X,44HNUMBER OF ROW POINTS (THESE ARE FITTED) IS  ,I5/
     2       1H ,8X,44HNUMBER OF COLUMN POINTS (THESE ARE GIVEN) IS,I5/
     3       1H ,8X,44HOPTION SET SELECTION:                       ,I5/
     4       1H ,8X,44H         0 = ALL ROWS SAME OPTION SET       /
     5       1H ,8X,44H         1 = ALL ROWS DIFFERENT OPTION SETS /
     6       1H ,8X,44H         2 = SPECIFIED ROWS SAME OPTION SET )
1603  FORMAT(1H ,8X,27HOPTION SETS START WITH ROWS,17X,4I5)
2603  FORMAT(1H ,//)
 604  FORMAT(1H ,4X,26H3 ANALYSIS SPECIFICATIONS://
     1       1H ,8X,44HTHE NUMBER OF DIMENSIONS IS                 ,I5/
     2       1H ,8X,44HMAXIMUM NUMBER OF ANALYSES IN ANY OPTION SET,I5/
     3       1H ,8X,44HPRELIMINARY TRANSFORMATION OF CONFIGURATION ,I5/
     4       1H ,8X,44H         0 = REMAINS UNCHANGED              /
     5       1H ,8X,44H         1 = WEIGHTED                       /
     6       1H ,8X,44H         2 = ROTATED AND WEIGHTED           /
     7       1H ,8X,44HSTANDARDIZE EXTERNAL DATA                   ,I5/
     8       1H ,8X,44H         0 = YES (=1+2)                     /
     9       1H ,8X,44H         1 = CENTER ONLY                    )
 1604 FORMAT(1H ,8X,44H         2 = NORMALIZE ONLY                 /
     2       1H ,8X,44H         3 = NONE OF THE ABOVE              )
 605  FORMAT(1H ,8X,44HMODELS:  0 = NOT APPLIED, 1 = APPLIED       /
     1       1H ,8X,44H             VECTOR MODEL                   ,I5/
     2       1H ,8X,44H             UNFOLDING MODEL                ,I5/
     3       1H ,8X,44H             WEIGHTED UNFOLDING MODEL       ,I5/
     4       1H ,8X,44H             GENERAL UNFOLDING MODEL        ,I5/
     5       1H ,8X,44HTHE NUMBER OF NON-METRIC ITERATIONS IS      ,I5/
     6       1H ,8X,39HTHE CONVERGENCE CRITERION IS           ,E10.2///)
 606  FORMAT(1H ,4X,44H4 PRINT/PLOT OPTIONS: 0 = NONE              //
     1       1H ,8X,44HPRINT INPUT:                                ,I5/
     2       1H ,8X,44H         1 = TARGET CONFIGURATION           /
     3       1H ,8X,44H         2 = EXTERNAL DATA                  /
     4       1H ,8X,44H         3 = 1 & 2                          /
     5       1H ,8X,44HPRINT RESULTS ACROSS OPTION SETS            ,I5/
     6       1H ,8X,44H         1 = COMPLETE RESULTS               /
     7       1H ,8X,44H         2 = FIT                            /
     8       1H ,8X,44H         3 = 2 & CRITERION-PREDICTED VALUES /
     9       1H ,8X,44HSELECTED RESULTS ACROSS ANALYSES            ,I5)
1606  FORMAT(1H ,8X,44HSCATTER&TRANSFORMATION PLOT FOR N ROWS,  N =,I5/
     2       1H ,8X,44HPAIRWISE PLOTS IDEAL POINTS/VECTORS,   DIM =,I5)
 607  FORMAT(1H ,8X,44HHISTORY OF NON-METRIC REGRESSION            ,I5/
     1       1H ,8X,44HSTATISTICS FOR METRIC ANALYSES              ,I5/
     2       1H ,8X,44HSTORE OUTPUT:                               ,I5/
     3       1H ,8X,44H         1 = COORDINATES                    /
     4       1H ,8X,44H         2 = LIKE 1 + WEIGHTS               /
     5       1H ,8X,44H         3 = LIKE 2 + ROTATION MATRICES     /
     6       1H ,8X,44HSTORE OUTPUT:                               ,I5/
     8       1H ,8X,44H         1 = PREDICTED VALUES               /
     9       1H ,8X,44H         2 = LIKE 1 + CRITERION             ///)
 608  FORMAT(1H ,4X,31H5 UNIT NUMBERS FOR INPUT/OUTPUT//
     1       1H ,8X,44HTHE CONFIGURATION WILL BE READ FROM UNIT    ,I5/
     2       1H ,8X,44HTHE EXTERNAL DATA WILL BE READ FROM UNIT    ,I5/
     3       1H ,8X,44HTHE OPTIONS WILL BE READ FROM UNIT          ,I5/
     4       1H ,8X,44HUNIT NUMBER FOR SCRATCH FILE 1              ,I5/
     5       1H ,8X,44HUNIT NUMBER FOR SCRATCH FILE 2              ,I5/
     6       1H ,8X,44HOUTPUT UNIT FOR THE COORDINATES             ,I5/
     7       1H ,8X,44HOUTPUT UNIT FOR THE WEIGHTS                 ,I5/
     8       1H ,8X,44HOUTPUT UNIT FOR THE ROTATION MATRICES       ,I5/
     9       1H ,8X,44HOUTPUT UNIT FOR PREDICTED&CRITERION VALUES  ,I5/)
  700 FORMAT(1H0///,9X,43HTHE PROGRAM STOPS BECAUSE OF A FATAL ERROR ,
     X                 43HIN THE PARAMETER CARD JUST PRINTED.        /
     X       1H ,8X,   41HCHECK NUMBER OF ROWS OR COLUMNS, OR UNIT ,
     X                 43HNUMBERS FOR 2 OBLIGATORY SCRATCH FILES.    )
  701 FORMAT(1H0///,5X,43H*** WARNING *** THE PRELIMINARY TRANSFORMAT,
     X                 43HION WILL NOT BE PERFORMED, BECAUSE THE     /
     X       1H0,20X,  42HNUMBER OF TARGET POINTS IS NOT SUFFICIENT.)
C
      END
      SUBROUTINE ESTIMA(DATA,DSTAR,PRM,PSI,NCOL,NPR,COEF,RR,KK)
C     ******************************************************************
C     *                                                                *
C     *  E S T I M A                                                   *
C     *                                                                *
C     *  PURPOSE: COMPUTES THE REGRESSION COEFFICIENTS, THE PREDICTED  *
C     *           VALUES AND THE MULTIPLE CORRELATION                  *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION DATA(NCOL),DSTAR(NCOL),PRM(NCOL,NPR),PSI(NCOL,NPR),
     X          COEF(NPR)
C
C FOR THE VECTOR MODEL THE DISSIMILARITIES WILL BE REVERSED
C
      DSYM=1.
      IF (KK.EQ.4) DSYM=-1.0
C
      DO 10 J=1,NPR
         SUM=0.0
         DO 20 I=1,NCOL
            SUM=SUM+PSI(I,J)*DATA(I)*DSYM
   20    CONTINUE
         COEF(J)=SUM
   10 CONTINUE
C
      RMEAN=0.0
      DO 30 I=1,NCOL
         SUM=0.0
         DO 40 K=1,NPR
            SUM=SUM+PRM(I,K)*COEF(K)*DSYM
   40    CONTINUE
         DSTAR(I)=SUM
         RMEAN=RMEAN+SUM
   30 CONTINUE
      RMEAN=RMEAN/FLOAT(NCOL)
      STANDP=0.0
      DO 60 I=1,NCOL
         A=DSTAR(I)-RMEAN
         STANDP=STANDP+A*A
   60 CONTINUE
      STANDP=SQRT(STANDP)
C
      RMEAN=0.0
      STANDD=0.0
      DO 70 I=1,NCOL
         RMEAN=RMEAN+DATA(I)
   70 CONTINUE
      RMEAN=RMEAN/FLOAT(NCOL)
      DO 80 I=1,NCOL
         A=DATA(I)-RMEAN
         STANDD=STANDD+A*A
   80 CONTINUE
      STANDD=SQRT(STANDD)
C
      RR=STANDP/STANDD
      RETURN
      END
      SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR)
C     ******************************************************************
C     *                                                                *
C     *  I M T Q L 2                                                   *
C     *                                                                *
C     *  PURPOSE: DETERMINES THE EIGENVALUES AND EIGENVECTORS OF A     *
C     *           SYMMETRIC TRIDIAGONAL MATRIX USING THE IMPLICIT      *
C     *           QL METHOD.                                           *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     *  ADAPTED FROM EISPACK                                          *
C     *                                                                *
C     ******************************************************************
      DIMENSION D(N),E(N),Z(NM,N)
      REAL MACHEP
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPCIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
C
      MACHEP = 2.0**(-20)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      E(N) = 0.0
C
      DO 240 L = 1, N
         J = 0
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
  105    DO 110 M = L, N
            IF (M .EQ. N) GO TO 120
            IF (ABS(E(M)) .LE. MACHEP* (ABS(D(M)) + ABS(D(M+1))))
     X         GO TO 120
  110 CONTINUE
C
  120 P = D(L)
      IF (M .EQ. L) GO TO 240
      IF (J .EQ.30) GO TO 1000
      J = J + 1
C     ********** FORM SHIFT **********
      G = (D(L+1) - P) / (2.0 * E(L))
      R = SQRT(G*G+1.0)
      G = D(M) - P + E(L) / (G + SIGN(R,G))
      S = 1.0
      C = 1.0
      P = 0.0
      NML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
      DO 200 II = 1, NML
      I = M - II
      F = S * E(I)
      B = C * E(I)
      IF (ABS(F) .LT. ABS(G)) GO TO 150
      C = G / F
      R = SQRT(C*C+1.0)
      E(I+1) = F * R
      S = 1.0 / R
      C = C * S
      GO TO 160
  150 S = F / G
      R = SQRT(S*S+1.0)
      E(I+1) = G* R
      C = 1.0 / R
      S = S * C
  160 G = D(I+1) - P
      R = (D(I) - G) * S + 2.0 * C * B
      P = S * R
      D(I+1) = G + P
      G = C * R - B
C     ********** FORM VECTOR **********
      DO 180 K = 1, N
         F = Z(K,I+1)
         Z(K,I+1) = S * Z(K,I) + C * F
         Z(K,I) = C * Z(K,I) - S * F
  180 CONTINUE
C
  200 CONTINUE
C
      D(L) = D(L) - P
      E(L) = G
      E(M) = 0.0
      GO TO 105
  240 CONTINUE
C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
      DO 300 II = 2, N
      I = II - 1
      K = I
      P = D(I)
C
         DO 260 J = II, N
      IF (D(J) .LE. P ) GO TO 260
            K = J
            P = D(J)
  260    CONTINUE
C
         IF (K .EQ. I) GO TO 300
         D(K) = D(I)
         D(I) = P
C
         DO 280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
  280    CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
 1000 IERR = L
 1001 RETURN
      END
      SUBROUTINE MAIN2(A,NWORDS)
C     ******************************************************************
C     *                                                                *
C     *  M A I N 2                                                     *
C     *                                                                *
C     *  PURPOSE: COMPUTES THE RELATIVE ADRESSES FOR ALL ARRAYS IN     *
C     *           THE REST OF THE PROGRAM PREFMAP-III                  *
C     *                                                                *
C     *  SUBROUTINES CALLED: MAIN3                                     *
C     *                                                                *
C     ******************************************************************
      DIMENSION  A(NWORDS)
      COMMON/CDECLA/NDIM,MAX,CRIT,IPRI1,IPRI2,IPLO,IHIS,ISHE,IWRI1,IST,
     X        IWRI2,IPRE,ISTA,NV,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2,
     X        NCR1,NW2,NW3,NW4,ND,NANA,IPRI0,MOD1,MOD2,MOD3,MOD4
      COMMON/CPARAN/NROW,NCOL,IOPT,IDIF(4)
C
      ICON   =    1
      IPR1   =    1 + NCOL*NDIM
      IPR2   =    IPR1 + NCOL*NPR1
      IPR3   =    IPR2 + NCOL*NPR2
      IPR4   =    IPR3 + NCOL*NPR3
      ISV1   =    IPR4 + NCOL*NPR4
      ISV2   =    ISV1 + NPR1
      ISV3   =    ISV2 + NPR2
      ISV4   =    ISV3 + NPR3
      ILSV1  =    ISV4 + NPR4
      ILSV2  =    ILSV1 + NCOL*NPR1
      ILSV3  =    ILSV2 + NCOL*NPR2
      ILSV4  =    ILSV3 + NCOL*NPR3
      IRSV1  =    ILSV4 + NCOL*NPR4
      IRSV2  =    IRSV1 + NPR1*NPR1
      IRSV3  =    IRSV2 + NPR2*NPR2
      IRSV4  =    IRSV3 + NPR3*NPR3
      ITCOOR =    IRSV4+ NPR4*NPR4
      ITWEIG =    ITCOOR+ NPL1*NDIM*NANA
      ITROTA =    ITWEIG+ NPL1*NDIM*NANA
      ICOEF  =    ITROTA+ NPL2*NDIM*NANA
      ICOOR  =    ICOEF + NPRM
      IWEIGH =    ICOOR + NDIM*NANA
      IROTAT =    IWEIGH+ NDIM*NANA
      ICROSS =    IROTAT+ NDIM*NANA*NDIM
      IWSIGN =    ICROSS+ NCR1*NCR1
      IDATA  =    IWSIGN+ NDIM
      IDISP  =    IDATA + NCOL
      IDSTAR =    IDISP + NCOL
      IIORD  =    IDSTAR+ NCOL
      ILBK   =    IIORD + NCOL
      IW1    =    ILBK  + NCOL
      IW2    =    IW1   + NPRM
      IW3    =    IW2   + NW2
      IW4    =    IW3   + NW3
      IOO    =    IW4   + NW4
      IDANA  =    IOO   + NANA
C
      CALL MAIN3(A(ICON),A(IPR1),A(IPR2),A(IPR3),A(IPR4),
     X     A(ISV1),A(ISV2),A(ISV3),A(ISV4),
     X     A(ILSV1),A(ILSV2),A(ILSV3),A(ILSV4),
     X     A(IRSV1),A(IRSV2),A(IRSV3),A(IRSV4),A(ITCOOR),
     X     A(ITWEIG),A(ITROTA),A(ICOEF),A(ICOOR),A(IWEIGH),
     X     A(IROTAT),A(ICROSS),A(IWSIGN),A(IDATA),A(IDISP),
     X     A(IDSTAR),A(IIORD),A(ILBK),A(IW1),A(IW2),A(IW3),A(IW4),
     X     A(IOO),A(IDANA),NDIM,MAX,CRIT,
     X     IPRI0,IPRI1,IPRI2,IPLO,IHIS,IST,ISHE,IWRI1,IWRI2,IPRE,ISTA,
     X     NROW,NCOL,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2,NCR1,
     X     NW2,NW3,NW4,ND,NANA,IOPT,IDIF,MOD1,MOD2,MOD3,MOD4)
C
      RETURN
      END
      SUBROUTINE MAIN3(CON,PR1,PR2,PR3,PR4,SIVA1,SIVA2,SIVA3,SIVA4,
     X                 LSV1,LSV2,LSV3,LSV4,RSV1,RSV2,RSV3,RSV4,
     X                 TCOOR,TWEIGH,TROTAT,COEF,COOR,WEIGH,ROTAT,
     X                 CROSS,WSIGN,DATA,DISP,DSTAR,IORD,LBK,WORK1,
     X                 WORK2,WORK3,WORK4,OPTION,DANA,
     X                 NDIM,MAX,CRIT,IPRI0,IPRI1,IPRI2,
     X                 IPLO,IHIS,IST,ISHE,IWRI1,IWRI2,IPRE,ISTA,
     X                 NROW,NCOL,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2,
     X                 NCR1,NW2,NW3,NW4,ND,NANA,IOPT,IDIF,MOD1,MOD2,
     X                 MOD3,MOD4)
C     ******************************************************************
C     *                                                                *
C     *  M A I N 3                                                     *
C     *                                                                *
C     *  PURPOSE: CONTROLS THE FLOW OF THE PREFMAP-III PROGRAM         *
C     *                                                                *
C     *  SUBROUTINES CALLED: PREREC,ZERO,PREPRO,PRED1(-4),PSINV,ZEROI, *
C     *                      ESTIMA,MONO,UNRAV,RESUL2,STATIS,OUTP3     *
C     *                                                                *
C     ******************************************************************
      COMMON/CFORM /FORM1(20)
      COMMON/CHEAD /TITLE(20)
      COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4)
C
      DIMENSION CON(NCOL,NDIM),PR1(NCOL,NPR1),PR2(NCOL,NPR2),
     X          PR3(NCOL,NPR3),PR4(NCOL,NPR4)
      DIMENSION SIVA1(NPR1),SIVA2(NPR2),SIVA3(NPR3),SIVA4(NPR4),
     X          LSV1(NCOL,NPR1),LSV2(NCOL,NPR2),LSV3(NCOL,NPR3),
     X          LSV4(NCOL,NPR4),RSV1(NPR1,NPR1),RSV2(NPR2,NPR2),
     X          RSV3(NPR3,NPR3),RSV4(NPR4,NPR4),TCOOR(NPL1,ND),
     X          TWEIGH(NPL1,ND),TROTAT(NPL2,ND),COEF(NPRM),DANA(NCOL),
     X          COOR(ND),WEIGH(ND),ROTAT(NDIM,ND),CROSS(NCR1,NCR1),
     X          WSIGN(NDIM),OPTION(NANA),OPTLAB(12),OPTLBS(12),IDIF(4),
     X          DATA(NCOL),DISP(NCOL),DSTAR(NCOL),IORD(NCOL),LBK(NCOL),
     X          WORK1(NPRM),WORK2(NW2),WORK3(NW3),WORK4(NW4),FITT(4),
     X          HULP(4),HULP2(4),FITM(12),KCOUNT(12),VAF(12),PVAF(12),
     X          SFITM(12)
      REAL LSV1,LSV2,LSV3,LSV4
      INTEGER OPTION,OPTLAB,OPTLBS,HULP,HULP2
      DATA OPTLAB/'GM  ','WM  ','UM  ','VM  ','GP  ','WP  ','UP  ',
     X            'VP  ','GS  ','WS  ','US  ','VS  '/
      DATA OPTLBS/'gm  ','wm  ','um  ','vm  ','gp  ','wp  ','up  ',
     X            'vp  ','gs  ','ws  ','us  ','vs  '/
C
C     READ FORMAT AND TARGET CONFIGURATION
C
      READ(INPARA,500) FORM1
      WRITE(IWRITE,635) FORM1
      READ(INCON,FORM1) ((CON(J,M),M=1,NDIM),J=1,NCOL)
C
C CENTER TARGET CONFIGURATION
C
      DO 61 M=1,NDIM
         RMEAN=0.0
         DO 62 J=1,NCOL
            RMEAN=RMEAN+CON(J,M)
   62    CONTINUE
         RMEAN=RMEAN/FLOAT(NCOL)
         DO 63 J=1,NCOL
            CON(J,M)=CON(J,M)-RMEAN
   63    CONTINUE
         WRITE (IWRITE,629) M,RMEAN
   61 CONTINUE
C
      IF (IPRI0.EQ.1.OR.IPRI0.EQ.3) WRITE (IWRITE,600)
      IF (IPRI0.EQ.1.OR.IPRI0.EQ.3)
     X   CALL PREREC(CON,NCOL,NDIM,NCOL,IWRITE,0)
C
C READ THE EXTERNAL DATA AND STORE ON SCRATCH FILE. IF PREPROCESSING
C OF THE TARGET CONFIGURATION IS DESIRED, ACCUMULATE THE MATRIX OF
C SUMS OF SQUARES AND CROSS PRODUCTS.
C
      REWIND ISCR1
      REWIND ISCR2
C
C     DUMMY WRITING TO SCRATCH FILES (NECESSARY FOR VAX-UNIX)
C
      WRITE(ISCR1) (K,K=1,NCOL)
      WRITE(ISCR2) (K,K=1,NANA)
      REWIND ISCR1
      REWIND ISCR2
C
      CALL ZERO(CROSS,NCR1*NCR1)
         CALL ZERO(TWEIGH,NPL1*ND)
         CALL ZERO(TCOOR,NPL1*ND)
         CALL ZERO(TROTAT,NPL2*ND)
      READ (INPARA,500) FORM1
      WRITE(IWRITE,636) FORM1
      IF (IPRI0.EQ.0.OR.IPRI0.EQ.1) KN=MIN0(NROW,10)
      IF (IPRI0.EQ.2.OR.IPRI0.EQ.3) KN=NROW
      KNN=MIN0(NCOL,10)
      WRITE(IWRITE,604) KN,(K,K=1,KNN)
      NN=NROW
      DO 5 I=1,NROW
         READ (INDATA,FORM1) (DATA(J),J=1,NCOL)
         IF (IPRI0.EQ.2.OR.IPRI0.EQ.3.OR.I.LE.10)
     X      WRITE (IWRITE,605) I,(DATA(J),J=1,NCOL)
C
C NORMALIZE DATA AND/OR COMPUTE SSQ
C
C THE PRELIMINARY TRANSFORMATION AND THE METRIC ANALYSES ARE PERFORMED
C ACCORDING TO THE DESIRED NORMALIZATION OF THE EXTERNAL DATA.
C THE NONMETRIC ANALYSES ARE PERFORMED ON THE BASIS OF THE EXTERNAL
C DATA, NORMALIZED TO HAVE ZERO MEAN AND VARIANCE EQUAL TO ONE.
C THEREFORE THE RAW DATA WILL BE WRITTEN TO THE SCRATCH FILE (ARRAY
C DATA) WHILE ARRAY DANA WILL FIRST CONTAIN THE DATA, NORMALIZED
C ACCORDING TO OPTION ISTA. FOR NONMETRIC ANALYSES THE RAW DATA WILL
C BE NORMALIZED AS DESCRIBED ABOVE, NO MATTER WHAT VALUE THE PARAMETER
C ISTA TAKES. ARRAY DISP WILL THEN BE USED TO OBTAIN THE OPTIMALLY
C TRANSFORMED NORMALIZED DATA. IN THIS WAY THE NORMALIZATION OF THE ORI-
C GINAL DATA IS CONTROLLED (E.G. PRESERVING THE ORIGINAL VARIANCE), WHILE
C FOR NONMETRIC ANALYSES WE OBTAIN AN APPROPRIATE SCALING TO MAKE TRANS-
C FORMATION PLOTS.
C
         RMEAN=0.0
         STAND=0.0
         VAR=0.0
         DO 19 J=1,NCOL
            RMEAN=RMEAN+DATA(J)
   19    CONTINUE
         RMEAN=RMEAN/FLOAT(NCOL)
         DO 29 J=1,NCOL
            A=DATA(J)-RMEAN
            IF (ISTA.EQ.0.OR.ISTA.EQ.1) DANA(J)=A
            IF (ISTA.EQ.2.OR.ISTA.EQ.3) DANA(J)=DATA(J)
            STAND=STAND+DANA(J)*DANA(J)
            VAR=VAR+A*A
   29    CONTINUE
         IF (VAR.EQ.0.0) WRITE(IWRITE,900) I
         IF (VAR.EQ.0.0) GO TO 6
         IF (ISTA.EQ.1) GO TO 18
         IF (ISTA.EQ.3) GO TO 18
         STAND=SQRT(STAND/FLOAT(NCOL))
         DO 39 J=1,NCOL
            DANA(J)=DANA(J)/STAND
   39    CONTINUE
   18    CONTINUE
C
         IF (IPRE.LE.0) GO TO 6
            DO 7 J=1,NCOL
               DO 8 K=1,NCOL
                  CROSS(J,K)=CROSS(J,K)+DANA(J)*DANA(K)
    8          CONTINUE
    7       CONTINUE
    6    CONTINUE
         WRITE (ISCR1) (DATA(J),J=1,NCOL)
    5 CONTINUE
C     NROW=NN
      REWIND ISCR1
C
      DO 9 I=1,NDIM
         WSIGN(I)=1.0
    9 CONTINUE
C
      IF (IPRE.GT.0) CALL PREPRO(CON,CROSS,PR1,LSV1,PR4,LSV4,RSV1,RSV4,
     X                           WORK1,WSIGN,IPRE,IWRITE,NCOL,NDIM,NPR1,
     X                           NPR4,NPRM)
      IF (IPRE.GT.0) WRITE (IWRITE,601)
      IF (IPRE.GT.0) CALL PREREC(CON,NCOL,NDIM,NCOL,
     X                                          IWRITE,0)
C
C COMPUTE MAXIMAL DISTANCE FROM ORIGIN
C
      CMAX=0.0
      DO 52 I=1,NCOL
         A=0.0
         DO 53 J=1,NDIM
            A=A+CON(I,J)*CON(I,J)
   53    CONTINUE
         IF (A.GT.CMAX) CMAX=A
   52 CONTINUE
      CMAX=SQRT(CMAX)
C
C     CHECK MODELS THAT WILL BE USED
C     AND ACCORDINGLY CONSTRUCT PREDICTOR MATRICES
C
      IF (MOD1.EQ.1) CALL PRED1(CON,PR1,LSV1,NCOL,NDIM,NPR1)
      IF (MOD2.EQ.1) CALL PRED2(CON,PR2,LSV2,NCOL,NDIM,NPR2)
      IF (MOD3.EQ.1) CALL PRED3(CON,PR3,LSV3,WSIGN,NCOL,NDIM,NPR3)
      IF (MOD4.EQ.1) CALL PRED4(CON,PR4,LSV4,NCOL,NDIM,NPR4)
C
C
C     COMPUTE THE PSEUDO INVERSES OF THE PREDICTOR MATRICES
C     ON BASIS OF THEIR SINGULAR VALUE DECOMPOSITIONS
C
      IF (MOD1.EQ.1) CALL PSINV(LSV1,RSV1,SIVA1,NPR1,NCOL,WORK1)
      IF (MOD2.EQ.1) CALL PSINV(LSV2,RSV2,SIVA2,NPR2,NCOL,WORK1)
      IF (MOD3.EQ.1) CALL PSINV(LSV3,RSV3,SIVA3,NPR3,NCOL,WORK1)
      IF (MOD4.EQ.1) CALL PSINV(LSV4,RSV4,SIVA4,NPR4,NCOL,WORK1)
C
      CALL ZEROI(HULP2,4)
      CALL ZEROI(KCOUNT,12)
      CALL ZERO(FITM,12)
      CALL ZERO(SFITM,12)
      CALL ZERO(VAF,12)
      CALL ZERO(PVAF,12)
C
      IO=1
      II=1
      WRITE (IWRITE,632)
      WRITE (IWRITE,631) (J,J=1,NANA)
      DO 1 I=1,NROW
         IF (IOPT.LE.0.AND.I.GT.1) GO TO 17
         IF (IOPT.EQ.2.AND.I.NE.IDIF(IO)) GO TO 17
         IF (IOPT.GT.1.AND.IO.LT.4) IO=IO+1
         READ (INOPT,501) (HULP(J),J=1,NANA)
         WRITE (IWRITE,633) II,(HULP(J),J=1,NANA)
         II=II+1
         DO 2 K=1,NANA
            DO 3 J=1,12
               IF (HULP(K).EQ.OPTLAB(J)) GO TO 4
               IF (HULP(K).EQ.OPTLBS(J)) GO TO 4
    3       CONTINUE
            HULP(K)=0
            GO TO 2
    4       HULP(K)=J
            IF (J.GT.8) HULP(K)=HULP(K)+12
            IF (J.GT.4.AND.J.LT.9) HULP(K)=HULP(K)+6
            IF (MOD(HULP(K),10).EQ.1.AND.MOD1.EQ.0)
     X        WRITE(IWRITE,630) K
            IF (MOD(HULP(K),10).EQ.1.AND.MOD1.EQ.0) HULP(K)=0
            IF (MOD(HULP(K),10).EQ.2.AND.MOD2.EQ.0)
     X        WRITE(IWRITE,630) K
            IF (MOD(HULP(K),10).EQ.2.AND.MOD2.EQ.0) HULP(K)=0
            IF (MOD(HULP(K),10).EQ.3.AND.MOD3.EQ.0)
     X        WRITE(IWRITE,630) K
            IF (MOD(HULP(K),10).EQ.3.AND.MOD3.EQ.0) HULP(K)=0
            IF (MOD(HULP(K),10).EQ.4.AND.MOD4.EQ.0)
     X        WRITE(IWRITE,630) K
            IF (MOD(HULP(K),10).EQ.4.AND.MOD4.EQ.0) HULP(K)=0
    2   CONTINUE
   17   CONTINUE
        WRITE(ISCR2) (HULP(K),K=1,NANA)
    1 CONTINUE
      REWIND ISCR2
C
C MAIN LOOP ACROSS ROWS
C
C  FOR EACH ROW SEPARATELY :
C  REREAD THE EXTERNAL DATA AND THE MODEL DESIRED
C
      DO 10 I=1,NROW
         READ(ISCR1) (DATA(J),J=1,NCOL)
         READ(ISCR2) (OPTION(J),J=1,NANA)
C
         CALL ZERO(COOR,NANA*NDIM)
         CALL ZERO(WEIGH,NANA*NDIM)
         CALL ZERO(ROTAT,NANA*NDIM*NDIM)
         CALL ZERO(FITT,NANA)
C
C
C INNER LOOP ACROSS OPTIONS
C
         DO 20 J=1,NANA
C
C NORMALIZE EXTERNAL DATA
C
         RMEAN=0.0
         STAND=0.0
         VAR=0.0
         DO 119 M=1,NCOL
            RMEAN=RMEAN+DATA(M)
  119    CONTINUE
         RMEAN=RMEAN/FLOAT(NCOL)
         DO 129 M=1,NCOL
            A=DATA(M)-RMEAN
            IF (ISTA.EQ.0.OR.ISTA.EQ.1) DANA(M)=A
            IF (ISTA.EQ.2.OR.ISTA.EQ.3) DANA(M)=DATA(M)
            STAND=STAND+DANA(M)*DANA(M)
            VAR=VAR+A*A
  129    CONTINUE
         IF (VAR.EQ.0.0) GO TO 10
         STAND=SQRT(STAND/FLOAT(NCOL))
         VAR=VAR/FLOAT(NCOL)
         IF (ISTA.EQ.1) GO TO 118
         IF (ISTA.EQ.3) GO TO 118
         RMEAN=0.0
         VAR=0.0
         DO 139 M=1,NCOL
            DANA(M)=DANA(M)/STAND
            RMEAN=RMEAN+DANA(M)
  139    CONTINUE
         RMEAN=RMEAN/FLOAT(NCOL)
         DO 140 M=1,NCOL
            A=DANA(M)-RMEAN
            VAR=VAR+A*A
  140    CONTINUE
         VAR=VAR/FLOAT(NCOL)
  118    CONTINUE
C
C  COMPUTE REGRESSION COEFFICIENTS, PREDICTED VALUES AND MULTIPLE
C  CORRELATION FOR THE METRIC CASE
C
         KK=MOD(OPTION(J),10)
         IF (KK.EQ.0) GO TO 20
         IF (KK.EQ.1.OR.KK.EQ.2) HULP2(J)=HULP2(J)+1
         IF (IPRI1.NE.1) GO TO 38
         IF (KK.EQ.1) WRITE (IWRITE,611) I,J
         IF (KK.EQ.2) WRITE (IWRITE,612) I,J
         IF (KK.EQ.3) WRITE (IWRITE,613) I,J
         IF (KK.EQ.4) WRITE (IWRITE,614) I,J
         IF (OPTION(J).LT.11.AND.OPTION(J).GT.0) WRITE (IWRITE,608)
         IF (OPTION(J).GT.10) WRITE (IWRITE,609)
         IF (OPTION(J).GT.10.AND.OPTION(J).LT.20) WRITE (IWRITE,620)
         IF (OPTION(J).GT.20) WRITE (IWRITE,621)
C
   38 IF(KK.EQ.1) CALL ESTIMA(DANA,DSTAR,PR1,LSV1,NCOL,NPR1,COEF,FIT,KK)
      IF(KK.EQ.2) CALL ESTIMA(DANA,DSTAR,PR2,LSV2,NCOL,NPR2,COEF,FIT,KK)
      IF(KK.EQ.3) CALL ESTIMA(DANA,DSTAR,PR3,LSV3,NCOL,NPR3,COEF,FIT,KK)
      IF(KK.EQ.4) CALL ESTIMA(DANA,DSTAR,PR4,LSV4,NCOL,NPR4,COEF,FIT,KK)
C
         WF=(FIT**2)*VAR
         IK=OPTION(J)
         IF (IK.GT.10.AND.IK.LT.15) IK=IK-6
         IF (IK.GT.20)              IK=IK-12
         VAF(IK)=VAF(IK)+WF
         PVAF(IK)=PVAF(IK)+VAR
         KCOUNT(IK)=KCOUNT(IK)+1
         IF (IPRI1.EQ.1) WRITE (IWRITE,615) FIT,VAR,WF
         IF (IPRI1.GT.1) WRITE (IWRITE,617) I,J,FIT,VAR,WF
         FITT(J)=FIT
C
C COPY THE DATA, APPROPRIATELY NORMALIZED IN ARRAY DANA,
C INTO ARRAY DISP (FOR TRANSFORMATION PLOTS)
C
         DO 37 M=1,NCOL
            DISP(M)=DANA(M)
   37    CONTINUE
C
C  NON-METRIC REGRESSION; STANDARDIZE ORIGINAL DATA
C
         IF (OPTION(J).LE.4) GO TO 36
         RMEAN=0.0
         VAR=0.0
         DO 219 M=1,NCOL
            RMEAN=RMEAN+DATA(M)
  219    CONTINUE
         RMEAN=RMEAN/FLOAT(NCOL)
         DO 229 M=1,NCOL
            DANA(M)=DATA(M)-RMEAN
            VAR=VAR+DANA(M)*DANA(M)
  229    CONTINUE
         STAND=SQRT(VAR/FLOAT(NCOL))
         SSQ=FLOAT(NCOL)
         DO 239 M=1,NCOL
            DANA(M)=DANA(M)/STAND
            DISP(M)=DANA(M)
  239    CONTINUE
         IF (OPTION(J).LE.14) IAPP=1
         IF (OPTION(J).GE.21) IAPP=2
         IF (OPTION(J).EQ.11.OR.OPTION(J).EQ.21)
     X      CALL MONO (DANA,DISP,PR1,LSV1,IORD,LBK,NCOL,NPR1,MAX,CRIT,
     X                 IAPP,IWRITE,IHIS,WORK2,FIT,SSQ)
         IF (OPTION(J).EQ.12.OR.OPTION(J).EQ.22)
     X      CALL MONO (DANA,DISP,PR2,LSV2,IORD,LBK,NCOL,NPR2,MAX,CRIT,
     X                 IAPP,IWRITE,IHIS,WORK2,FIT,SSQ)
         IF (OPTION(J).EQ.13.OR.OPTION(J).EQ.23)
     X      CALL MONO (DANA,DISP,PR3,LSV3,IORD,LBK,NCOL,NPR3,MAX,CRIT,
     X                 IAPP,IWRITE,IHIS,WORK2,FIT,SSQ)
         IF (OPTION(J).EQ.14.OR.OPTION(J).EQ.24)
     X      CALL MONO (DANA,DISP,PR4,LSV4,IORD,LBK,NCOL,NPR4,MAX,CRIT,
     X                 IAPP,IWRITE,IHIS,WORK2,FIT,SSQ)
C
         IF (OPTION(J).EQ.11.OR.OPTION(J).EQ.21)
     X      CALL ESTIMA(DISP,DSTAR,PR1,LSV1,NCOL,NPR1,COEF,FIT,KK)
         IF (OPTION(J).EQ.12.OR.OPTION(J).EQ.22)
     X      CALL ESTIMA(DISP,DSTAR,PR2,LSV2,NCOL,NPR2,COEF,FIT,KK)
         IF (OPTION(J).EQ.13.OR.OPTION(J).EQ.23)
     X      CALL ESTIMA(DISP,DSTAR,PR3,LSV3,NCOL,NPR3,COEF,FIT,KK)
         IF (OPTION(J).EQ.14.OR.OPTION(J).EQ.24)
     X      CALL ESTIMA(DISP,DSTAR,PR4,LSV4,NCOL,NPR4,COEF,FIT,KK)
C
         IF (IPRI1.EQ.1) WRITE (IWRITE,616) FIT
         IF (IPRI1.GT.1) WRITE (IWRITE,618) FIT
C
   36    CONTINUE
C
         IK=OPTION(J)
         IF (IK.GT.10.AND.IK.LT.15) IK=IK-6
         IF (IK.GT.20)              IK=IK-12
         FITM(IK)=FITM(IK)+FIT
         SFITM(IK)=SFITM(IK)+FIT**2
C
C  UNRAVEL THE REGRESSION COEFFICIENTS INTO THE PARAMETERS
C  OF THE MODEL DESIRED
C
         JJ=J
         CALL UNRAV(COEF,NPRM,COOR,WEIGH,ROTAT,OPTION(J),NDIM,ND,JJ,I,
     X   IWRITE,WSIGN,IPRE,WORK1,WORK3,WORK2,WORK4,CMAX,IPRI1,SLOP,CEPT,
     X   FIT)
         CALL RESUL2 (DANA,DISP,DSTAR,WORK2,WORK3,I,J,OPTION(J),
     X   WORK4,NCOL,NCOL*2,IWRITE,IWRI2,LUN(4),ISHE,IPRI1,SLOP,CEPT)
   20   CONTINUE
C
C STORE COORDINATES,WEIGHTS,ROTATION MATRICES
C
      IF (IWRI1.EQ.0) GO TO 21
      LU=LUN(1)
      DO 23 J=1,NCOL
         IF (I.EQ.1) WRITE(LU,679) J,(CON(J,K),K=1,NDIM)
  23  CONTINUE
      KE=0
      DO 24 J=1,NANA
         KB=KE+1
         KE=KB+NDIM-1
         WRITE(LU,680) I,J,(COOR(K),K=KB,KE)
  24  CONTINUE
      IF (IWRI1.EQ.1) GO TO 21
      LU=LUN(2)
      KE=0
      DO 25 J=1,NANA
         KB=KE+1
         KE=KB+NDIM-1
         KK=MOD(OPTION(J),10)
         IF (KK.LT.4.AND.KK.GT.0) WRITE(LU,680) I,J,(WEIGH(K),K=KB,KE)
  25  CONTINUE
      IF (IWRI1.EQ.2) GO TO 21
      LU=LUN(3)
      KE=0
      DO 26 J=1,NANA
         KB=KE+1
         KE=KB+NDIM-1
         KK=MOD(OPTION(J),10)
         DO 27 M=1,NDIM
            IF (KK.EQ.1) WRITE(LU,685) I,J,M,(ROTAT(M,K),K=KB,KE)
  27     CONTINUE
  26  CONTINUE
C
C COMPUTE STATISTICS FOR METRIC ANALYSES
C
   21 IF(IST.EQ.1) CALL STATIS(OPTION,NROW,NCOL,NDIM,FITT,NANA,I,IWRITE)
C
C
        IF (NPL1.EQ.1.AND.NROW.GT.1) GO TO 10
        DO 30 N=1,ND
           TCOOR(I,N)=COOR(N)
           TWEIGH(I,N)=WEIGH(N)
   30   CONTINUE
        IF (NPL2.EQ.NDIM.AND.NROW.GT.1) GO TO 10
        DO 35 M=1,NDIM
           MM=M+(I-1)*NDIM
           DO 40 N=1,ND
              TROTAT(MM,N)=ROTAT(M,N)
   40      CONTINUE
   35   CONTINUE
C
C
   10 CONTINUE
C
      WRITE(IWRITE,634)
      DO 54 K=1,12
         IF (KCOUNT(K).EQ.0) GO TO 54
         SFITM(K)=SQRT(SFITM(K)/FLOAT(KCOUNT(K)))
         FITM(K)=FITM(K)/FLOAT(KCOUNT(K))
         WRITE(IWRITE,6343) OPTLAB(K),KCOUNT(K),FITM(K),SFITM(K)
         IF (K.GT.4) GO TO 54
         WRITE(IWRITE,6342) PVAF(K),VAF(K)
         PVAF(K)=VAF(K)/PVAF(K)
         WRITE(IWRITE,6341) PVAF(K)
   54 CONTINUE
      NTOT=NROW+NCOL
      N=MAX0(IPLO,IPRI2)
      IF (N.GT.0) CALL OUTP3(TCOOR,CON,WORK2,WORK3,WORK4,IPLO,NROW,
     X   NCOL,NDIM,NTOT,ND,NPL2,TWEIGH,TROTAT,NANA,HULP2,IPRI2,
     X   OPTION,CMAX)
C
  500 FORMAT(20A4)
  501 FORMAT(4A4)
  629 FORMAT(1H ,4X,42HTHE TARGET CONFIGURATION WILL BE CENTERED:,
     X     1X,23HMEAN ORIGINAL DIMENSION,I5,F11.3)
  600 FORMAT(1H0//,5X,31HTARGET CONFIGURATION (CENTERED)/,
     X       1H ,  4X,31H-------------------------------)
  601 FORMAT(1H0//,5X,41HTARGET CONFIGURATION AFTER TRANSFORMATION/
     X       1H ,  4X,41H-----------------------------------------)
  604 FORMAT(1H0//,5X,I5,26H ROWS OF THE EXTERNAL DATA/
     X         1H ,4X,31H-------------------------------/1H ,9X,10I9)
  605 FORMAT(1H ,5X,I5,2X,10F9.3/(1H ,12X,10F9.3))
  608 FORMAT(1H+,51X,18H METRIC REGRESSION/1H ,4X,
     X  52H====================================================,
     X  13H=============)
  609 FORMAT(1H+,51X,22H NON-METRIC REGRESSION)
  611 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2,
     X              25H GENERAL UNFOLDING MODEL )
  612 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2,
     X              25H WEIGHTED UNFOLDING MODEL)
  613 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2,
     X              25H UNFOLDING MODEL         )
  614 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2,
     X              25H VECTOR MODEL            )
  615 FORMAT(1H0/,16X,14HMETRIC FIT    ,F11.3,10H  VARIANCE,F11.3,
     X                                        10H    V.A.F.,F11.3)
  617 FORMAT(1H0,3X,3HROW,I5,7H OPTION,I2,14H METRIC FIT   ,F6.3,
     X                    10H  VARIANCE,F11.3,10H    V.A.F.,F11.3)
  618 FORMAT(1H0,16X,15H NON-METRIC FIT,3X,F6.3)
  616 FORMAT(1H0/,16X,14HNON-METRIC FIT,F11.3)
  620 FORMAT(1H+,74X,27H   PRIMARY APPROACH TO TIES/1H ,4X,
     X  52H====================================================,
     X  45H=============================================)
  621 FORMAT(1H+,74X,27H SECONDARY APPROACH TO TIES/1H ,4X,
     X  52H====================================================,
     X  45H=============================================)
  630 FORMAT(1H0/,5X,22H*** WARNING *** OPTION,I2,12H WILL NOT BE,
     X    51H APPLIED. EITHER THE MODEL HAS NOT BEEN REFERRED TO/
     X    1H0,20X,9HON CARD 3,
     X    49H OR THE NUMBER OF TARGET POINTS IS NOT SUFFICIENT,
     X    19H TO FIT THIS MODEL.)
  631 FORMAT(1H0//,19X,8HANALYSIS/1H ,4X,11HOPTION SET:,4I4)
  632 FORMAT(1H0/,5X,35HMODEL:      V  = VECTOR            /
     X       1H ,4X, 35H            U  = UNFOLDING         /
     X       1H ,4X, 35H            W  = WEIGHTED UNFOLDING/
     X       1H ,4X, 35H            G  = GENERAL UNFOLDING /
     X       1H ,4X, 35HREGRESSION:  M = METRIC            /
     X       1H ,4X, 29H             P = NON-METRIC, ,
     X       24HPRIMARY APPROACH TO TIES/
     X       1H ,4X, 29H             S = NON-METRIC, ,
     X       26HSECONDARY APPROACH TO TIES)
  633 FORMAT(1H ,10X,I5,2X,4A4)
  634 FORMAT(1H0//,5X,41HSUMMARY OF RESULTS: (PROPORTION OF TOTAL ,
     X 53HVARIANCE ACCOUNTED FOR ONLY GIVEN FOR METRIC OPTIONS)/)
 6343 FORMAT(1H0,4X,9HMODEL '  ,1A4,5H' N =,I5,15H  AVERAGE FIT =,F6.3/
     X                           18X,25H  ROOT MEAN SQUARED FIT =,F6.3)
 6342 FORMAT(1H ,4X,16HTOTAL VARIANCE =,F9.3,3X,14HTOTAL VARIANCE,
     X        16H ACCOUNTED FOR =,F9.3)
 6341 FORMAT(1H ,4X,44HPROPORTION OF TOTAL VARIANCE ACCOUNTED FOR =,
     X        F6.3)
  635 FORMAT(1H0//,5X,30HTHE CONFIGURATION WILL BE READ,
     X                13H WITH FORMAT ,20A4/(1H ,4X,20A4))
  636 FORMAT(1H0//,5X,30HTHE EXTERNAL DATA WILL BE READ,
     X                13H WITH FORMAT ,20A4/(1H ,4X,20A4))
  679 FORMAT(I10,5F12.5/(10X,5F12.5))
  680 FORMAT(2I5,5F12.5/(10X,5F12.5))
  685 FORMAT(3I5,4F12.5/(15X,4F12.5))
  900 FORMAT(1H0//,5X,3HROW,I5,35H WILL BE OMITTED FROM THE ANALYSIS ,
     X      26HSINCE IT HAS ZERO VARIANCE//)
C
C
      RETURN
      END
      SUBROUTINE MONO(DATA,DISP,PRM,LSV,IORD,LBK,NCOL,NPR,IMAX,CRIT,
     X                IAPP,IWRITE,IHIS,HELP,FIT,DSTAND)
C     ******************************************************************
C     *                                                                *
C     *  M O N O                                                       *
C     *                                                                *
C     *  PURPOSE: CONTROLS THE MONOTONE REGRESSION                     *
C     *                                                                *
C     *  SUBROUTINES CALLED: RANK1,PROJEC,MONOR                        *
C     *                                                                *
C     ******************************************************************
      DIMENSION DATA(NCOL),DISP(NCOL),IORD(NCOL),LBK(NCOL),HELP(NCOL)
      DIMENSION PRM(NCOL,NPR),LSV(NCOL,NPR)
      REAL LSV
C
      OLDFIT=FIT
      NITER=0
      IF (IHIS.EQ.1) WRITE (IWRITE,605)
      CALL RANK1(DATA,IORD,LBK,HELP,IAPP,NCOL)
      DO 10 NITER=1,IMAX
         CALL PROJEC(PRM,LSV,DISP,HELP,NCOL,NPR,STAND1)
         CALL MONOR(DISP,HELP,IORD,LBK,IAPP,NCOL,STAND2,DSTAND)
         FIT=STAND2/STAND1
         DIFF=FIT-OLDFIT
         IF (IHIS.EQ.1) WRITE(IWRITE,610) NITER,FIT,DIFF
         IF (ABS(DIFF).LT.CRIT) GO TO 60
         OLDFIT=FIT
   10 CONTINUE
   60 RETURN
C
  605 FORMAT (1H0/,20X,32HHISTORY OF NON-METRIC REGRESSION/,
     X       1H ,19X,32H--------------------------------/,
     X       1H ,50X,15HDIFFERENCE WITH/,
     X       1H ,19X,9HITERATION,7X,3HFIT,12X,19HPRECEDING ITERATION)
  610 FORMAT (1H0,25X,I4,1X,F12.5,3X,F12.5)
      END
      SUBROUTINE MONOR(DISP,HELP,IORD,LBK,IAPP,NCOL,FIT,DSTAND)
C     ******************************************************************
C     *                                                                *
C     *  M O N O R                                                     *
C     *                                                                *
C     *  PURPOSE: MONITOR FOR MONOTONE REGRESSION                      *
C     *                                                                *
C     *  SUBROUTINES CALLED: PRAR,SCAR,MRMNH                           *
C     *                                                                *
C     ******************************************************************
      DIMENSION DISP(NCOL),HELP(NCOL),IORD(NCOL),LBK(NCOL)
      REAL MEAN
      DO 10 I=1,NCOL
         HELP(I)=DISP(IORD(I))
   10 CONTINUE
      IF (IAPP.EQ.1) CALL PRAR(HELP,IORD,LBK,NCOL)
      IF (IAPP.EQ.2) CALL SCAR(HELP,LBK,NCOL)
      CALL MRMNH(HELP,NCOL)
      STAND=0.0
      MEAN=0.0
      FIT=0.0
      DO 15 I=1,NCOL
         STAND=STAND+HELP(I)*HELP(I)
         MEAN=MEAN+HELP(I)
   15 CONTINUE
      STAND=SQRT(DSTAND/STAND)
      MEAN=MEAN/FLOAT(NCOL)
      DO 20 I=1,NCOL
         DISP(IORD(I))=HELP(I)*STAND
         HELP(I)=HELP(I)-MEAN
         FIT=FIT+HELP(I)*HELP(I)
   20 CONTINUE
      FIT=SQRT(FIT)
      RETURN
      END
      SUBROUTINE MRMNH (DISP,N)
C     ******************************************************************
C     *                                                                *
C     *  M R M N H                                                     *
C     *                                                                *
C     *  PURPOSE: MONOTONE REGRESSION.                                 *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE.                                     *
C     *                                                                *
C     *  AUTHOR: ERNST VAN WANING.     RELEASED: JULY 1976             *
C     *                                                                *
C     ******************************************************************
C
C
      DIMENSION DISP(N)
C
      LOVBKH = 0
      LUPH   = 0
      DO 1 IUP = 2, N
         IF (DISP(IUP) .GE. DISP(IUP-1)) GOTO 1
            SDS = DISP(IUP)
            SW = 1.
            IDOWN = IUP
    5          IDOWN = IDOWN - 1
               IF (LUPH .EQ. IDOWN) GOTO 2
                  SDS = SDS + DISP(IDOWN)
                  SW = SW + 1.
               GOTO 3
    2             SDS = SDS + (LOVBKH * DISP(LUPH))
                  SW = SW + LOVBKH
                  IDOWN = IDOWN - LOVBKH + 1
    3          TRIALV = SDS / SW
               IF (IDOWN .EQ. 1) GOTO 4
            IF (DISP(IDOWN - 1) .GT. TRIALV) GOTO 5
    4       DO 6 J = IDOWN, IUP
    6          DISP(J) = TRIALV
            LOVBKH = IUP - IDOWN + 1
            LUPH = IUP
    1    CONTINUE
      RETURN
      END
      SUBROUTINE OUTP3(TCOOR,CON,H1,H2,IND,IPLO,NROW,NCOL,NDIM,NTOT,ND,
     X           NPL2,TWEIGH,TROTAT,NANA,IWE,IPRI2,OPTION,CMAX)
C     ******************************************************************
C     *                                                                *
C     *  O U T P 3                                                     *
C     *                                                                *
C     *  PURPOSE: PROVIDES FOR THE PRINTING, AND PLOTTING OF RESULTS   *
C     *                                                                *
C     *  SUBROUTINES CALLED: PLOT                                      *
C     *                                                                *
C     ******************************************************************
      COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4)
      DIMENSION TCOOR(NROW,ND),CON(NCOL,NDIM),H1(NTOT),H2(NTOT),IWE(4),
     X          IND(NTOT),TWEIGH(NROW,ND),TROTAT(NPL2,ND),
     X          OPTION(NANA),OPTLAB(13)
      INTEGER OPTION,OPTLAB
      DATA OPTLAB/'GM  ','WM  ','UM  ','VM  ','GP  ','WP  ','UP  ',
     X            'VP  ','GS  ','WS  ','US  ','VS  ','N.A.'/
      REWIND ISCR2
      IF (IPRI2.EQ.0) GO TO 21
      KE=0
      DO 1 J=1,NANA
         KB=KE+1
         KE=KB+NDIM-1
         DO 5 I=1,NROW
            READ (ISCR2) (OPTION(JN),JN=1,NANA)
            IF (I.EQ.1) WRITE(IWRITE,600) J
            IO=OPTION(J)
            IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
            IF (IO.GT.20)              IO=IO-12
            IF (IO.EQ.0)               IO=13
            WRITE (IWRITE,603) I,OPTLAB(IO),(TCOOR(I,K),K=KB,KE)
    5    CONTINUE
         REWIND ISCR2
         IC=0
         DO 10 I=1,NROW
            READ (ISCR2) (OPTION(JN),JN=1,NANA)
            IO=OPTION(J)
            IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
            IF (IO.GT.20)              IO=IO-12
            IF (IO.EQ.0)               IO=13
            IK=MOD(OPTION(J),10)
            IF (IK.LT.4.AND.IK.GT.0) IC=IC+1
            IF (IC.EQ.1) WRITE(IWRITE,601) J
            IF (IC.EQ.1) IC=IC+1
            IF (IK.LT.4.AND.IK.GT.0) WRITE(IWRITE,603)
     X                         I,OPTLAB(IO),(TWEIGH(I,K),K=KB,KE)
   10    CONTINUE
         REWIND ISCR2
         IC=0
         MM=1
         DO 15 I=1,NROW
            READ (ISCR2) (OPTION(JN),JN=1,NANA)
            IO=OPTION(J)
            IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
            IF (IO.GT.20)              IO=IO-12
            IF (IO.EQ.0)               IO=13
            IK=MOD(OPTION(J),10)
            IF (IK.LT.2.AND.IK.GT.0) IC=IC+1
            IF (IC.EQ.1) WRITE(IWRITE,602) J
            IF (IC.EQ.1) IC=IC+1
            DO 19 M=1,NDIM
               IF (IK.LT.2.AND.IK.GT.0.AND.M.EQ.1) WRITE(IWRITE,604)
     X                         I,OPTLAB(IO),M,(TROTAT(MM,K),K=KB,KE)
               IF (IK.LT.2.AND.IK.GT.0.AND.M.GT.1) WRITE(IWRITE,614)
     X                                      M,(TROTAT(MM,K),K=KB,KE)
               MM=MM+1
   19       CONTINUE
   15    CONTINUE
         REWIND ISCR2
    1 CONTINUE
   21 IF (IPLO.EQ.0) RETURN
      DO 25 M=1,NANA
         NLB=(M-1)*NDIM
         NDPL=MIN0(NDIM,IPLO)
         IF (NDPL.EQ.1) GO TO 3
         DO 20 I=2,NDPL
            IMIN=I-1
            DO 30 J=1,IMIN
               JB=J+NLB
               IB=I+NLB
               NA=0
 35            L=0
               NC=0
               NC2=0
               DO 40 K=1,NROW
                  READ (ISCR2) (OPTION(JN),JN=1,NANA)
                  IO=OPTION(M)
                  IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
                  IF (IO.GT.20)              IO=IO-12
                  IF (IO.EQ.0)               IO=13
                  NOUT=0
                  NT=1-NA
                  IF (TWEIGH(K,JB).LE.0.0.OR.TWEIGH(K,IB).LE.0.0)
     X               NT=NA
                  NC=NC+NT
                  NC2=NC2+NT
                  L=L+1
                  H1(L)=TCOOR(K,JB)*FLOAT(NT)
                  H2(L)=TCOOR(K,IB)*FLOAT(NT)
                  TDIS=SQRT(H1(L)**2+H2(L)**2)
                  IF (TDIS.GT.2*CMAX) NOUT=1
                  IF (NOUT.EQ.1) H1(L)=0.0
                  IF (NOUT.EQ.1) H2(L)=0.0
                  IF (NOUT.EQ.1) NC2=NC2-1
                  IF (NOUT.EQ.1) WRITE (IWRITE,605) M,K,
     X               OPTLAB(IO),TCOOR(K,JB),TCOOR(K,IB)
 40            CONTINUE
               REWIND ISCR2
               DO 50 K=1,NCOL
                  L=L+1
                  H1(L)=CON(K,J)
                  H2(L)=CON(K,I)
 50            CONTINUE
               IF (NC.EQ.0) GO TO 51
               IF (NC2.EQ.0) GO TO 51
               WRITE (IWRITE,110) M,J,I
               CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1)
               IF (NC.EQ.NROW) GO TO 30
 51            NA=1
               L=0
               NC=0
               NC2=0
               DO 52 K=1,NROW
                  READ (ISCR2) (OPTION(JN),JN=1,NANA)
                  IO=OPTION(M)
                  IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
                  IF (IO.GT.20)              IO=IO-12
                  IF (IO.EQ.0)               IO=13
                  NOUT=0
                  NT=1-NA
                  IF (TWEIGH(K,JB).LT.0.0.OR.TWEIGH(K,IB).LT.0.0)
     X               NT=NA
                  NC=NC+NT
                  NC2=NC2+NT
                  L=L+1
                  H1(L)=TCOOR(K,JB)*FLOAT(NT)
                  H2(L)=TCOOR(K,IB)*FLOAT(NT)
                  TDIS=SQRT(H1(L)**2+H2(L)**2)
                  IF (TDIS.GT.2*CMAX) NOUT=1
                  IF (NOUT.EQ.1) H1(L)=0.0
                  IF (NOUT.EQ.1) H2(L)=0.0
                  IF (NOUT.EQ.1) NC2=NC2-1
                  IF (NOUT.EQ.1) WRITE (IWRITE,605) M,K,
     X               OPTLAB(IO),TCOOR(K,JB),TCOOR(K,IB)
 52            CONTINUE
               REWIND ISCR2
               IF (NC.EQ.0) GO TO 30
               IF (NC2.EQ.0) GO TO 30
               WRITE (IWRITE,111) M,J,I
               CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1)
 30         CONTINUE
 20      CONTINUE
         IF (IWE(M).LE.1) GO TO 25
         DO 61 I=2,NDPL
            IMIN=I-1
            DO 62 J=1,IMIN
               JB=J+NLB
               IB=I+NLB
               WRITE (IWRITE,210) M,J,I
               CALL PRPLOT(TWEIGH(1,JB),TWEIGH(1,IB),IND,NROW,NROW,
     X                   1,IWRITE,1)
   62       CONTINUE
   61    CONTINUE
         GO TO 25
C                                                                      C
C  ONE-DIMENSIONAL PLOT                                                C
C                                                                      C
 3       L=0
         NA=0
         NC=0
         NC2=0
         JB=NLB+1
         DO 63 K=1,NROW
            READ (ISCR2) (OPTION(JN),JN=1,NANA)
            IO=OPTION(M)
            IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
            IF (IO.GT.20)              IO=IO-12
            IF (IO.EQ.0)               IO=13
            NOUT=0
            NT=1-NA
            IF (TWEIGH(K,JB).LE.0.0) NT=NA
            L=L+1
            NC=NC+NT
            NC2=NC2+NT
            H1(L)=TCOOR(K,JB)*FLOAT(NT)
            H2(L)=0.0
            TDIS=ABS(H1(L))
            IF (TDIS.GT.2*CMAX) NOUT=1
            IF (NOUT.EQ.1) H1(L)=0.0
            IF (NOUT.EQ.1) NC2=NC2-1
            IF (NOUT.EQ.1) WRITE (IWRITE,606) M,K,
     X                     OPTLAB(IO),TCOOR(K,JB)
 63      CONTINUE
         REWIND ISCR2
         DO 64 K=1,NCOL
            L=L+1
            H1(L)=CON(K,1)
            H2(L)=0.0
 64      CONTINUE
         IF (NC.EQ.0) GO TO 26
         IF (NC2.EQ.0) GO TO 26
         WRITE (IWRITE,120) M
         CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1)
         IF (NC.EQ.NROW) GO TO 69
 26      NA=1
         L=0
         NC=0
         NC2=0
         DO 65 K=1,NROW
            READ (ISCR2) (OPTION(JN),JN=1,NANA)
            IO=OPTION(M)
            IF (IO.GT.10.AND.IO.LT.15) IO=IO-6
            IF (IO.GT.20)              IO=IO-12
            IF (IO.EQ.0)               IO=13
            NOUT=0
            NT=1-NA
            IF (TWEIGH(K,JB).LT.0.0) NT=NA
            L=L+1
            NC=NC+NT
            NC2=NC2+NT
            H1(L)=TCOOR(K,JB)*FLOAT(NT)
            H2(L)=0.0
            TDIS=ABS(H1(L))
            IF (TDIS.GT.2*CMAX) NOUT=1
            IF (NOUT.EQ.1) H1(L)=0.0
            IF (NOUT.EQ.1) NC2=NC2-1
            IF (NOUT.EQ.1) WRITE (IWRITE,606) M,K,
     X                     OPTLAB(IO),TCOOR(K,JB)
 65      CONTINUE
         REWIND ISCR2
         IF (NC.EQ.0) GO TO 69
         IF (NC2.EQ.0) GO TO 69
         WRITE(IWRITE,121) M
         CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1)
 69      IF (IWE(NANA).LE.1) GO TO 25
         DO 70 K=1,NROW
            H2(K)=0.0
 70      CONTINUE
         WRITE(IWRITE,220) M
         CALL PRPLOT(TWEIGH(1,JB),H2,IND,NROW,NROW,1,IWRITE,1)
 25   CONTINUE
C
         RETURN
C
  600 FORMAT(1H0///,12X,8HANALYSIS,I2,12H COORDINATES/1H ,11X,
     X       22H----------------------)
  601 FORMAT(1H0//,12X,8HANALYSIS,I2,8H WEIGHTS/1H ,11X,
     X       18H------------------)
  602 FORMAT(1H0//,12X,8HANALYSIS,I2,18H ROTATION MATRICES/1H ,11X,
     X       28H----------------------------)
  603 FORMAT(1H ,25X,I6,5X,A4,10F9.3/(1H ,40X,10F9.3))
  604 FORMAT(1H ,25X,I6,5X,A4,I2,10F7.3/(1H ,42X,10F7.3))
  614 FORMAT(1H ,40X,I2,10F7.3/(1H ,42X,10F7.3))
  605 FORMAT(1H0//,6X,8HANALYSIS,I2,20H *** NOT PLOTTED ***,I6,5X,A4,
     X       2F9.3,30H *** SEVERELY OUT OF RANGE ***)
  606 FORMAT(1H0//,6X,8HANALYSIS,I2,20H *** NOT PLOTTED ***,I6,5X,A4,
     X        F9.3,30H *** SEVERELY OUT OF RANGE ***)
 110  FORMAT(1H1,3X,38HIDEAL POINTS AND/OR VECTORS (INTEGERS),
     X              26H IN TARGET CONFIGURATION. ,
     X       8HANALYSIS,I2,4H DIM,I2,16H (X-AXIS) VS DIM,I2,9H (Y-AXIS))
 111  FORMAT(1H1,3X,42HANTI-IDEAL AND/OR SADDLE POINTS (INTEGERS),
     X              26H IN TARGET CONFIGURATION. ,
     X       8HANALYSIS,I2,4H DIM,I2,16H (X-AXIS) VS DIM,I2,9H (Y-AXIS))
 120  FORMAT(1H1,3X,43HONE-DIMENSIONAL PLOT OF IDEAL POINTS AND/OR,
     X       32H VECTORS (INTEGERS) FOR ANALYSIS,I2)
 121  FORMAT(1H1,3X,48HONE-DIMENSIONAL PLOT OF ANTI-IDEAL AND/OR SADDLE,
     X       20H POINTS FOR ANALYSIS,I2)
 210  FORMAT(1H1,3X,28HPLOT OF WEIGHTS FOR ANALYSIS,
     1       I2,10H DIMENSION,I2,26H (X-AXIS) VERSUS DIMENSION,I2,
     2        9H (Y-AXIS))
 220  FORMAT(1H1,3X,31HONE-DIMENSIONAL PLOT OF WEIGHTS,
     X       13H FOR ANALYSIS,I2)
      END
      SUBROUTINE PRAR(DIST,IORD,LBK,N)
C     ******************************************************************
C     *                                                                *
C     *  P R A R                                                       *
C     *                                                                *
C     *  PURPOSE:  TO SORT THE DISTANCES WITHIN TIEBLOCKS AND REORDER  *
C     *            THEIR INDICES IN IORD ACCORDINGLY                   *
C     *                                                                *
C     *  SUBROUTINES CALLED: SHEL9                                     *
C     *                                                                *
C     ******************************************************************
      DIMENSION DIST(N),IORD(N),LBK(N)
C
      IDLBK=0
      JB=1
 1       IDLBK=IDLBK+1
         NLBK=LBK(IDLBK)
         IF (NLBK.GT.1) CALL SHEL9(DIST(JB),IORD(JB),NLBK)
         JB=JB+NLBK
      IF (JB.LT.N) GO TO 1
      RETURN
      END
      SUBROUTINE PRED1(CON,PRM,LSV1,NCOL,NDIM,NPR1)
C     ******************************************************************
C     *                                                                *
C     *  P R E D 1                                                     *
C     *                                                                *
C     *  PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE              *
C     *           GENERAL UNFOLDING MODEL                              *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR1),LSV1(NCOL,NPR1)
      REAL LSV1
C
      DO 10 J=1,NCOL
         PRM(J,1)=1.0
         LSV1(J,1)=PRM(J,1)
         DO 20 I=1,NDIM
            L=I+1
            LL=L+NDIM
            PRM(J,L)=CON(J,I)
            LSV1(J,L)=PRM(J,L)
            PRM(J,LL)=CON(J,I)**2
            LSV1(J,LL)=PRM(J,LL)
   20    CONTINUE
         LB=2
         LE=NDIM-1
         L=2*NDIM+2
         DO 30 M=1,LE
            DO 40 I=LB,NDIM
               PRM(J,L)=CON(J,M)*CON(J,I)
               LSV1(J,L)=PRM(J,L)
               L=L+1
   40       CONTINUE
            LB=LB+1
   30    CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE PRED2(CON,PRM,LSV2,NCOL,NDIM,NPR2)
C     ******************************************************************
C     *                                                                *
C     *  P R E D 2                                                     *
C     *                                                                *
C     *  PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE              *
C     *           WEIGHTED UNFOLDING MODEL                             *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR2),LSV2(NCOL,NPR2)
      REAL LSV2
C
      DO 10 J=1,NCOL
         PRM(J,1)=1.0
         LSV2(J,1)=1.0
         DO 20 I=1,NDIM
            L=I+1
            LL=L+NDIM
            PRM(J,L)=CON(J,I)
            LSV2(J,L)=PRM(J,L)
            PRM(J,LL)=CON(J,I)**2
            LSV2(J,LL)=PRM(J,LL)
   20    CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE PRED3(CON,PRM,LSV3,WSIGN,NCOL,NDIM,NPR3)
C     ******************************************************************
C     *                                                                *
C     *  P R E D 3                                                     *
C     *                                                                *
C     *  PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE              *
C     *           UNFOLDING MODEL                                      *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR3),LSV3(NCOL,NPR3)
      DIMENSION WSIGN(NDIM)
      REAL LSV3
C
      DO 10 J=1,NCOL
         PRM(J,1)=1.0
         LSV3(J,1)=1.0
         SSQ=0.
         DO 20 I=1,NDIM
            L=I+1
            PRM(J,L)=CON(J,I)
            LSV3(J,L)=PRM(J,L)
            SSQ=SSQ+WSIGN(I)*CON(J,I)**2
   20    CONTINUE
         L=NDIM+2
         PRM(J,L)=SSQ
         LSV3(J,L)=SSQ
   10 CONTINUE
      RETURN
      END
      SUBROUTINE PRED4(CON,PRM,LSV4,NCOL,NDIM,NPR4)
C     ******************************************************************
C     *                                                                *
C     *  P R E D 4                                                     *
C     *                                                                *
C     *  PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE              *
C     *           VECTOR MODEL                                         *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR4),LSV4(NCOL,NPR4)
      REAL LSV4
C
      DO 10 J=1,NCOL
         PRM(J,1)=1.0
         LSV4(J,1)=1.0
         DO 20 I=1,NDIM
            L=I+1
            PRM(J,L)=CON(J,I)
            LSV4(J,L)=PRM(J,L)
   20    CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE PREPRO(CON,CROSS,U,UH,V,PSEU,WORK1,WORK4,W,WNEG,
     X                 IPRE,IWRITE,NCOL,NDIM,NPR1,NPR4,NPRM)
C     ******************************************************************
C     *                                                                *
C     *  P R E P R O                                                   *
C     *                                                                *
C     *  PURPOSE: PREPROCESSING THE TARGET CONFIGURATION               *
C     *           IF IPRE EQUALS 1, ONLY WEIGHTING OF AXES             *
C     *           IF IPRE EQUALS 2, CANONICAL REFERENCE FRAME          *
C     *                                                                *
C     *  SUBROUTINES CALLED: PRED4,PSINV,TRED2,IMTQL2,PREREC           *
C     *                                                                *
C     ******************************************************************
      DIMENSION CON(NCOL,NDIM),CROSS(NCOL,NCOL),U(NCOL,NPR1),
     X          UH(NCOL,NPR1),V(NCOL,NPR4),PSEU(NCOL,NPR4),
     X          WORK1(NPR1,NPR1),WORK4(NPR4,NPR4),W(NPRM),WNEG(NDIM)
C
C  CHECK OPTION AND SET UP U,V AND PSEU
C
      IF (IPRE.EQ.1) NRAV=NDIM
      IF (IPRE.EQ.2) NRAV=NPR1-(NDIM+1)
      CALL PRED4(CON,V,PSEU,NCOL,NDIM,NPR4)
      CALL PSINV(PSEU,WORK4,W,NPR4,NCOL,WORK1)
      DO 10 J=1,NCOL
         DO 20 L=1,NDIM
            U(J,L)=CON(J,L)*CON(J,L)
 20      CONTINUE
         IF (IPRE.EQ.1) GO TO 10
            ME=NDIM-1
            KB=2
            L=NDIM+1
            DO 30 M=1,ME
               DO 40 K=KB,NDIM
                  U(J,L)=2*CON(J,M)*CON(J,K)
                  L=L+1
 40            CONTINUE
               KB=KB+1
 30         CONTINUE
 10   CONTINUE
C
C  FORM UH
C
      DO 50 I=1,NCOL
         DO 60 L=1,NRAV
            UH(I,L)=0.0
            DO 70 J=1,NCOL
               DO 80 K=1,NPR4
                  UH(I,L)=UH(I,L)+PSEU(I,K)*V(J,K)*U(J,L)
 80            CONTINUE
 70         CONTINUE
            UH(I,L)=U(I,L)-UH(I,L)
 60      CONTINUE
 50   CONTINUE
C
C
C  USE U AS WORKING ARRAY FOR INVERSE SQUARE ROOT OF UHU
C
      DO 90 L=1,NRAV
         DO 100 M=1,NRAV
            WORK1(L,M)=0.0
            DO 110 I=1,NCOL
               WORK1(L,M)=WORK1(L,M)+UH(I,L)*UH(I,M)
 110        CONTINUE
 100     CONTINUE
 90   CONTINUE
      CALL TRED2(NPR1,NRAV,W,WORK4,WORK1)
      CALL IMTQL2(NPR1,NRAV,W,WORK4,WORK1,IERR)
      IF (IERR.GT.0) GO TO 900
      SMALL=0.0000001
      DO 120 L=1,NRAV
         IF (W(L).LE.SMALL) W(L)=0.0
         IF (W(L).GT.SMALL) W(L)=1.0/SQRT(W(L))
 120  CONTINUE
      DO 130 L=1,NRAV
         DO 140 M=1,NRAV
            U(L,M)=0.0
            DO 150 N=1,NRAV
               U(L,M)=U(L,M)+WORK1(L,N)*W(N)*WORK1(M,N)
 150        CONTINUE
 140     CONTINUE
 130  CONTINUE
C
C
C  FORM MATRIX FOR GENERALIZED EIGENVECTOR PROBLEM
C
      DO 160 L=1,NRAV
         DO 170 M=1,NRAV
            WORK1(L,M)=0.0
            DO 180 K1=1,NRAV
               DO 190 I=1,NCOL
                  DO 200 J=1,NCOL
                     DO 210 K2=1,NRAV
                        P=U(L,K1)*UH(I,K1)*CROSS(I,J)*UH(J,K2)*U(K2,M)
                        WORK1(L,M)=WORK1(L,M)+P
 210                 CONTINUE
 200              CONTINUE
 190           CONTINUE
 180        CONTINUE
 170     CONTINUE
 160  CONTINUE
C
C
C  FIND SOLUTION FOR W
C
      CALL TRED2(NPR1,NRAV,W,WORK4,WORK1)
      CALL IMTQL2(NPR1,NRAV,W,WORK4,WORK1,IERR)
      IF (IERR.GT.0) GO TO 900
      SUM=0.0
      DO 220 L=1,NRAV
         W(L)=0.0
         DO 230 M=1,NRAV
            W(L)=W(L)+U(L,M)*WORK1(M,1)
 230     CONTINUE
         IF (L.LE.NDIM) SUM=SUM+W(L)
 220  CONTINUE
      DO 235 L=1,NRAV
         W(L)=W(L)*SIGN(1.0,SUM)
  235 CONTINUE
C
      IF (IPRE.EQ.2) GO TO 260
C
C  FINISH HERE IF WEIGHTING ONLY
C
      STAND=0.0
      DO 236 L=1,NRAV
         IF (L.LE.NDIM) STAND=STAND+ABS(W(L))
 236  CONTINUE
      STAND=FLOAT(NDIM)/STAND
      DO 240 L=1,NDIM
         WNEG(L)=SIGN(1.,W(L))
         W(L)=SQRT(WNEG(L)*W(L)*STAND)
         DO 250 I=1,NCOL
            CON(I,L)=CON(I,L)*W(L)
 250     CONTINUE
 240  CONTINUE
      WRITE (IWRITE,601)
      WRITE (IWRITE,602) (W(L),L=1,NDIM)
      RETURN
C
C  FINISH HERE FOR FULL TRANSFORMATION
C
 260  CONTINUE
      DO 270 L=1,NDIM
         WORK4(L,L)=W(L)
 270  CONTINUE
      ME=NDIM-1
      KB=2
      L=NDIM+1
      DO 280 M=1,ME
         DO 290 K=KB,NDIM
            WORK4(M,K)=W(L)
            WORK4(K,M)=W(L)
            L=L+1
 290     CONTINUE
         KB=KB+1
 280  CONTINUE
      CALL TRED2(NPR4,NDIM,W,WORK1,WORK4)
      CALL IMTQL2(NPR4,NDIM,W,WORK1,WORK4,IERR)
      IF (IERR.GT.0) GO TO 900
C
C  TRANSMIT SIGNS OF AXES THROUGH WNEG AND PERFORM REGULAR
C  TRANSFORMATION HERE.
C
      STAND=0.0
      DO 295 L=1,NDIM
         STAND=STAND+ABS(W(L))
 295  CONTINUE
      STAND=FLOAT(NDIM)/STAND
      DO 300 L=1,NDIM
         WNEG(L)=SIGN(1.,W(L))
         W(L)=SQRT(WNEG(L)*W(L)*STAND)
 300  CONTINUE
      WRITE (IWRITE,601)
      WRITE (IWRITE,602) (W(L),L=1,NDIM)
      WRITE (IWRITE,603)
      CALL PREREC(WORK4,NDIM,NDIM,NPR4,IWRITE,0)
      DO 310 I=1,NCOL
         DO 320 L=1,NDIM
            WORK1(1,L)=0.0
            DO 330 M=1,NDIM
               WORK1(1,L)=WORK1(1,L)+CON(I,M)*WORK4(M,L)
 330        CONTINUE
 320     CONTINUE
         DO 340 L=1,NDIM
            CON(I,L)=WORK1(1,L)*W(L)
 340     CONTINUE
 310  CONTINUE
      RETURN
C
C  TROUBLE REPORT
C
 900  WRITE(IWRITE,901)
 901  FORMAT(1H1,33HT R O U B L E                    ,
     X      /1H0,34HONE OF THE EIGENVALUE ROUTINES    ,
     X      /1H0,34HIN TRANSFORMATION STEP BREAKS DOWN,
     X      /1H0,34HCONFIGURATION LEFT UNCHANGED      )
 601  FORMAT(1H0//,5X,32HNORMALIZED WEIGHTS PER DIMENSION,
     X      /1H ,  4X,32H--------------------------------)
 602  FORMAT(1H ,12X,12F9.3)
 603  FORMAT(1H0///5X,15HROTATION MATRIX/1H ,4X,15H---------------)
      RETURN
      END
      SUBROUTINE PREREC(X,M,N,MM,IWRITE,NR)
C     ******************************************************************
C     *                                                                *
C     *  P R E F R E C                                                 *
C     *                                                                *
C     *  PURPOSE: PRINTS A M*N ARRAY WITH ROW- AND COLUMN IDENTIFI-    *
C     *           CATION, EACH LINE CONTAINING AT MOST 10 VALUES (F9.3)*
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION X(MM,N)
C                                                                      C
C  CHECK HOW MANY TIMES THE NUMBER OF COLUMNS EXCEEDS 10               C
C                                                                      C
      NTRU=N/10
      NREM=N-NTRU*10
      IF (NREM.GT.0) GO TO 1
         NREM=10
         NTRU=NTRU-1
 1    NREP=NTRU+1
C                                                                      C
C  PRINT NREP TIMES A BLOCK OF ELEMENTS                                C
C                                                                      C
      NK=10
      KB=1
      DO 2 I=1,NREP
         IF (I.EQ.NREP) NK=NREM
         KE=KB+NK-1
         WRITE(IWRITE,60) (K,K=KB,KE)
         WRITE(IWRITE,61)
         DO 3 L=1,M
            IF (NR.EQ.0) WRITE(IWRITE,62) L,(X(L,K),K=KB,KE)
            IF (NR.NE.0) WRITE(IWRITE,62) NR,(X(L,K),K=KB,KE)
 3       CONTINUE
         KB=KB+NK
 2    CONTINUE
      RETURN
 60   FORMAT(1H0,9X,10I9)
 61   FORMAT(1H )
 62   FORMAT(1H ,I10,2X,10F9.3)
      END
      SUBROUTINE PROJEC(PRM,LSV,DISP,WORK,NCOL,NPR,STAND)
C     ******************************************************************
C     *                                                                *
C     *  P R O J E C T                                                 *
C     *                                                                *
C     *  PURPOSE: COMPUTES THE PROJECTION OF THE PREDICTAND ON THE     *
C     *           PREDICTOR                                            *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION DISP(NCOL),PRM(NCOL,NPR),LSV(NCOL,NPR),WORK(NCOL)
      REAL LSV,MEAN
C
      DO 5 I=1,NCOL
         WORK(I)=0.0
    5 CONTINUE
      DO 10 K=1,NCOL
         DO 20 I=1,NCOL
            SUM=0.0
            DO 30 J=1,NPR
               SUM=SUM+PRM(K,J)*LSV(I,J)*DISP(K)
   30       CONTINUE
            WORK(I)=WORK(I)+SUM
   20    CONTINUE
   10 CONTINUE
      MEAN=0.0
      STAND=0.0
      DO 40 I=1,NCOL
         DISP(I)=WORK(I)
         MEAN=MEAN+DISP(I)
   40 CONTINUE
      MEAN=MEAN/FLOAT(NCOL)
      DO 50 I=1,NCOL
         A=WORK(I)-MEAN
         STAND=STAND+A*A
   50 CONTINUE
      STAND=SQRT(STAND)
      RETURN
      END
      SUBROUTINE PRPLOT(X,Y,IND,NITEM,NRC,IDENT,IWRITE,NOR)
C     ******************************************************************
C     *                                                                *
C     *  P R P L O T                                                   *
C     *                                                                *
C     *  PURPOSE: PRINTPLOT OF Y VERSUS X                              *
C     *           IF IDENT.GT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY*
C     *           INTEGER NUMBERS 1-9,0,1-9... , THE REST OF THE POINTS*
C     *           BY CHARACTERS A-I,J,A-I...                           *
C     *           IF IDENT.EQ.0, ALL POINTS ARE PRINTED AS STARS       *
C     *           IF IDENT.LT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY*
C     *           STARS, THE OTHERS BY A 'D'                           *
C     *           IF NOR.EQ.1, THE ORIGIN WILL BE PLOTTED              *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     *  AUTHOR: WILLEM HEISER       RELEASED: MARCH 1978              *
C     *                              ADAPTED : NOVEMBER 1982           *
C     *                              ADAPTED : DECEMBER 1984           *
C     *                                                                *
C     ******************************************************************
C
C TO ADJUST THE PLOT CHANGE THE CONSTANTS UNDER WHICH
C ## IS PLACED ACCORDING TO THE COMMENTS
C
      DIMENSION LINE(71),XPR(11),X(NITEM),Y(NITEM),LABEL(20),IND(NITEM)
C                            ##=NN+1 (NN=INT(NC/7))
C                    ##=NUMBER OF COLUMNS PER LINE=NC
      INTEGER BLANK,STAR
      DATA BLANK,STAR,MORE/1H ,1H*,1HM/,NL/56/,NC/71/,NN/10/
C                                                 ##=NC##=NN
C                                          ##=NUMBER OF LINES PER PAGE
      DATA LABEL/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
     X           1HJ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI/
C
 1    FORMAT(16X,3H.--,10(7H+------),5H+---.)
C                      ##=NN
 2    FORMAT(8X,F7.3,1X,1HI,2X,71A1,3X,1HI)
C                              ##=NC
C                           ##=NC+5
 4    FORMAT(15X,11F7.3)
C                ##=NN+1
C
         DO 10 I=1,NITEM
 10         IND(I)=I
C
C  THE INDICES ARE ORDERED SUCH THAT THEY WOULD ORDER Y IN DESCENDING
C  ORDER
C
      M=NITEM
 20      M=M/2
         IF (M.LT.1) GO TO 40
            K=NITEM-M
            J=1
 41            I=J
 49               L=I+M
                  IF (Y(IND(I)).GE.Y(IND(L))) GO TO 60
                     II=IND(I)
                     IND(I)=IND(L)
                     IND(L)=II
                  I=I-M
               IF (I.GE.1) GO TO 49
 60         J=J+1
      IF (J-K) 41,41,20
 40   CONTINUE
C
C  THE PLOT IS ADJUSTED TO THE RANGE OF THE 'LONGEST' AXIS
C
      XMAX=X(1)
      XMIN=X(1)
      DO 70 I=2,NITEM
         IF (X(I).GT.XMAX) XMAX=X(I)
         IF (X(I).LT.XMIN) XMIN=X(I)
 70      CONTINUE
      IF (XMIN.GT. 0.0.AND.NOR.EQ.1) XMIN=0.0
      IF (XMAX.LT. 0.0.AND.NOR.EQ.1) XMAX=0.0
      SPANX=XMAX-XMIN
      YMAX=Y(IND(1))
      IF (YMAX.LT. 0.0.AND.NOR.EQ.1) YMAX=0.0
      YMIN=Y(IND(NITEM))
      IF (YMIN.GT. 0.0.AND.NOR.EQ.1) YMIN=0.0
      SPANY=YMAX-YMIN
      IF (SPANY.LE.SPANX) GO TO 75
         XMIN=XMIN-(SPANY-SPANX)/2.0
         SPANX=SPANY
 75   STEPX=SPANX/(7.0*NN)
      STEPY=(STEPX*NC)/(NL+0.)
      HSTEP=STEPY/2.0
      TOP=(YMIN+YMAX+SPANX)/2.0
C
C  THE POINTS ARE PLOTTED LINE AFTER LINE
C
      WRITE(IWRITE,1)
      L=1
      I=1
      YPR=TOP+STEPY
      IF (NOR.EQ.1) IP=(0.0-XMIN)/STEPX+1.0
 80      YPR=YPR-STEPY
 90            DO 95 J=1,NC
 95               LINE(J)=BLANK
         IF ((YPR-HSTEP.LE.0.0.AND.YPR+HSTEP.GT.0.0).AND.NOR.EQ.1)
     X      LINE(IP)=STAR
         IF (I.GT.NITEM) GO TO 110
            IF ((YPR-Y(IND(I))).GT.HSTEP) GO TO 110
 100           JP=(X(IND(I))-XMIN)/STEPX+1.0
               IF (LINE(JP).EQ.BLANK) GO TO 101
                  LINE(JP)=MORE
                  GO TO 103
 101           IF (IDENT.EQ.0) GO TO 102
               IF (IDENT.GT.0) GO TO 104
                  IF (IND(I).LE.NRC) LINE(JP)=STAR
                  IF (IND(I).GT.NRC) LINE(JP)=LABEL(15)
                  GO TO 103
 104              IF (IND(I).LE.NRC) LINE(JP)=LABEL(MOD(IND(I),10)+1)
                  IF (IND(I).GT.NRC)
     X               LINE(JP)=LABEL(MOD((IND(I)-NRC),10)+11)
                  GO TO 103
 102           LINE(JP)=STAR
 103           I=I+1
               IF (I.GT.NITEM) GO TO 105
                  IF ((YPR-Y(IND(I))).LT.HSTEP) GO TO 100
 105           WRITE(IWRITE,2) YPR,LINE
               GO TO 115
 110     WRITE(IWRITE,2) YPR,LINE
 115     L=L+1
      IF (L.LE.NL) GO TO 80
         WRITE(IWRITE,1)
C
C  VALUES OF X-AXIS ARE PRINTED
C
      XPR(1)=XMIN
      DO 130 J=1,10
 130     XPR(J+1)=XPR(J)+STEPX*7.0
      WRITE(IWRITE,4) XPR
      RETURN
      END
      SUBROUTINE PSINV(LSV,RSV,SIVA,NPR,NCOL,WORK)
C     ******************************************************************
C     *                                                                *
C     *  P S I N V                                                     *
C     *                                                                *
C     *  PURPOSE: COMPUTES THE PSEUDO INVERSE                          *
C     *                                                                *
C     *                -1                                              *
C     *           (A'A)  A'    ON THE BASIS OF THE                     *
C     *                                                                *
C     *           SINGULAR VALUE DECOMPOSITION A = KSL' BY TAKING      *
C     *                                                                *
C    *            ********************                                 *
C     *                 -1        -1                                   *
C     *            (A'A)  A' =  LS  K'                                 *
C     *            ************                                        *
C     *                                                                *
C     *  SUBROUTINES CALLED: SIVAD                                     *
C     *                                                                *
C     ******************************************************************
      DIMENSION LSV(NCOL,NPR),RSV(NPR,NPR),SIVA(NPR),WORK(NPR)
      REAL LSV
C
      IF (NPR.GT.1) CALL SIVAD
     X (NCOL,NPR,SIVA,LSV,NCOL*NPR,RSV,NPR*NPR,WORK,ICV)
C
      DO 5 I=1,NPR
         SMALL=.000001
         IF(SIVA(I).GE.SMALL) SIVA(I)=1.0/SIVA(I)
         IF(SIVA(I).LT.SMALL) SIVA(I)=0.
    5 CONTINUE
C
C
      DO 10 K=1,NCOL
         DO 20 I=1,NPR
            SUM=0.0
            DO 30 J=1,NPR
               SUM=SUM+RSV(I,J)*SIVA(J)*LSV(K,J)
   30       CONTINUE
            WORK(I)=SUM
   20    CONTINUE
         DO 40 I=1,NPR
            LSV(K,I)=WORK(I)
   40    CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE RANK1(DATA,IORD,LBK,HELP,IAPP,NCOL)
C     ******************************************************************
C     *                                                                *
C     * R A N K 1                                                      *
C     *                                                                *
C     *  PURPOSE: CONSTRUCT IORD AND LBK FROM DATA.                    *
C     *                                                                *
C     * SUBROUTINES CALLED: SHEL9, TIEBL                               *
C     *                                                                *
C     ******************************************************************
      DIMENSION DATA(NCOL),IORD(NCOL),LBK(NCOL),HELP(NCOL)
      DO 10 I=1,NCOL
         HELP(I)=DATA(I)
         IORD(I)=I
         LBK(I)=0
   10 CONTINUE
      CALL SHEL9(HELP,IORD,NCOL)
      CALL TIEBL(HELP,LBK,NCOL,IDLBK)
      IF (IDLBK.EQ.NCOL) IAPP=0
      RETURN
      END
      SUBROUTINE RESUL1(COOR,WEIGH,ROTAT,NDIM,ND,NOP,IWRITE,OPT,NR,IP)
C     ******************************************************************
C     *                                                                *
C     *  R E S U L 1                                                   *
C     *                                                                *
C     *  PURPOSE: PRINTS RESULTS PER SUBJECT                           *
C     *                                                                *
C     *  SUBROUTINES CALLED: PREREC                                   *
C     *                                                                *
C     ******************************************************************
      DIMENSION COOR(ND),WEIGH(ND),ROTAT(NDIM,ND)
      INTEGER OPT
C
      NOPT=NOP-1
C
      WRITE(IWRITE,600)
      ILB=NOPT*NDIM+1
      IF (OPT.NE.3) WRITE(IWRITE,606)
      IF (OPT.NE.3) GO TO 5
      IF (IP.EQ.0.AND.WEIGH(ILB).GE.0.0) WRITE(IWRITE,603)
      IF (IP.EQ.0.AND.WEIGH(ILB).LT.0.0) WRITE(IWRITE,604)
      IF (IP.EQ.0) GO TO 5
      WRITE (IWRITE,605)
      K=ILB
      KE=ILB+NDIM-1
      WRITE(IWRITE,607) (WEIGH(K),K=ILB,KE)
      WRITE(IWRITE,608)
    5 CALL PREREC(COOR(ILB),1,NDIM,1,IWRITE,NR)
      IF (OPT.LT.3) WRITE(IWRITE,601)
      IF (OPT.LT.3) CALL PREREC(WEIGH(ILB),1,NDIM,1,IWRITE,NR)
      IF (OPT.LT.2) WRITE(IWRITE,602)
      IF (OPT.LT.2) CALL PREREC(ROTAT(1,ILB),NDIM,NDIM,NDIM,IWRITE,0)
C
  600 FORMAT(1H0/,16X,11HCOORDINATES)
  606 FORMAT(1H ,15X,11H-----------)
  601 FORMAT(1H0/,16X,18HNORMALIZED WEIGHTS/
     X       1H , 15X,18H------------------)
  602 FORMAT(1H0/,16X,8HROTATION/1H ,15X,8H--------)
  603 FORMAT(1H+,27X,11HIDEAL POINT/1H ,15X,23H-----------------------)
  604 FORMAT(1H+,27X,16HANTI-IDEAL POINT/1H ,15X,
     X       28H----------------------------)
C
  605 FORMAT(1H+,27X,13HWITH WEIGHTS:)
  607 FORMAT(1H+,40X,20F4.0)
  608 FORMAT(1H ,15X,25H-------------------------)
      RETURN
      END
      SUBROUTINE RESUL2(DATA,DISP,DSTAR,H1,H2,NI,NJ,OPTION,IND,NCOL,
     X                  NC,IWRITE,IW,LU,ISHE,IPRI1,SLOPE,CEPT)
C     ******************************************************************
C     *                                                                *
C     *  R E S U L 2                                                   *
C     *                                                                *
C     *  PURPOSE: PLOTS DATA VERSUS DISP AND DSTAR (OPTIONAL)          *
C     *           STORES DSTAR AND DISP (OPTIONAL)                     *
C     *                                                                *
C     *  SUBROUTINES CALLED: PLOT                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION DATA(NCOL),DISP(NCOL),DSTAR(NCOL),H1(NC),H2(NC),IND(NC)
      INTEGER OPTION,OPT
C
      OPT=MOD(OPTION,10)
      IPRI=0
      IF (IPRI1.EQ.1.OR.IPRI1.EQ.3) IPRI=1
      IF (IPRI.EQ.1) WRITE (IWRITE,601)
      IF (IPRI.EQ.1) CALL PREREC(DISP,1,NCOL,1,IWRITE,NI)
      IF (IPRI.EQ.1) WRITE (IWRITE,603)
      IF (IPRI.EQ.1) CALL PREREC(DSTAR,1,NCOL,1,IWRITE,NI)
      IF (IPRI.EQ.1) WRITE (IWRITE,604) SLOPE,CEPT
      IF (ISHE.LT.NI) GO TO 30
      L=0
      DO 10 I=1,NCOL
         L=L+1
         H1(L)=DATA(I)
         H2(L)=DISP(I)
   10 CONTINUE
      DO 20 I=1,NCOL
         L=L+1
         H1(L)=DATA(I)
         H2(L)=DSTAR(I)
   20 CONTINUE
      WRITE(IWRITE,600) NI
      IF (OPT.EQ.1) WRITE(IWRITE,611)
      IF (OPT.EQ.2) WRITE(IWRITE,612)
      IF (OPT.EQ.3) WRITE(IWRITE,613)
      IF (OPT.EQ.4) WRITE(IWRITE,614)
      CALL PRPLOT(H1,H2,IND,NC,NCOL,-1,IWRITE,0)
   30 IF (IW.EQ.0) RETURN
      DO 31 J=1,NCOL
         IF (IW.EQ.1) WRITE (LU,620) NI,NJ,J,DSTAR(J)
         IF (IW.EQ.2) WRITE (LU,620) NI,NJ,J,DSTAR(J),DATA(J),DISP(J)
   31 CONTINUE
C
  600 FORMAT(1H1,3X,44HSCATTER PLOT. DATA (X-AXIS) VS CRITERION (*),
     X      38H AND PREDICTED VALUES (D) (M=*+D), ROW,I5)
  611 FORMAT(1H+,90X,25H GENERAL UNFOLDING MODEL )
  612 FORMAT(1H+,90X,25H WEIGHTED UNFOLDING MODEL)
  613 FORMAT(1H+,90X,25H UNFOLDING MODEL         )
  614 FORMAT(1H+,90X,25H VECTOR MODEL            )
  601 FORMAT(1H0/,16X,16HCRITERION VALUES,
     X   59H (EXTERNAL DATA, POSSIBLY TRANSFORMED BY MONOTONE FUNCTION)/
     X   1H ,15X,16H----------------)
  603 FORMAT(1H0/,16X,16HPREDICTED VALUES/1H ,15X,16H----------------)
  604 FORMAT(1H0/,16X,5HSLOPE,F11.5,3X,9HINTERCEPT,F11.5)
  620 FORMAT(3I5,3F13.6)
      RETURN
      END
      SUBROUTINE SCAR(DIST,LBK,N)
C     ******************************************************************
C     *                                                                *
C     *  S C A R                                                       *
C     *                                                                *
C     *  PURPOSE:  POOL DISTANCES WITHIN TIEBLOCKS                     *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION DIST(N),LBK(N)
C
      IDLBK=0
      JB=0
 3       IDLBK=IDLBK+1
         NLBK=LBK(IDLBK)
         MAX=JB+NLBK
         MIN=JB+1
         DMEAN=0.0
         DO 1 K=MIN,MAX
            DMEAN=DMEAN+DIST(K)
 1          CONTINUE
         DMEAN=DMEAN/FLOAT(NLBK)
         DO 2 K=MIN,MAX
            DIST(K)=DMEAN
 2          CONTINUE
         JB=MAX
      IF (JB.LT.N) GO TO 3
      RETURN
      END
      SUBROUTINE SHEL9(DIS,IND,N)
C     ******************************************************************
C
      DIMENSION DIS(N) ,IND(N)
C
      M=N
    4 M=M/2
      IF (M .LT. 1) RETURN
      K=N-M
      J=1
    3 I=J
    2    L=I+M
         IF ( DIS(L) .GE. DIS(I)) GOTO 1
            B = DIS(I)
            DIS(I) = DIS(L)
            DIS(L) = B
            KK = IND(I)
            IND(I) = IND(L)
            IND(L) = KK
            I=I-M
         IF (I .GE. 1) GOTO 2
    1 J=J+1
      IF (J-K) 3,3,4
      END
      SUBROUTINE SIVAD(M,N,Q,U,IU,V,IV,E,ICV)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C****                                                                  C
C**** NAME:         SIVAD(M,N,Q,U,IU,V,IV,E,ICV)                       C
C****                                                                  C
C****                                                                  C
C**** PURPOSE    :  COMPUTES AND ORDERS THE SINGULAR VALUES OF A REAL  C
C****               MATRIX A AND GIVES A COMPLETE ORTHOGONAL DECOMPO-  C
C****               SITION:                                            C
C****                           A = U * Q * V(T)                       C
C****                                                                  C
C****                                                                  C
C**** PARAMETERS :  M         NUMBER OF ROWS OF THE MATRIX A           C
C****                                                                  C
C****               N         NUMBER OF COLUMNS OF THE MATRIX A        C
C****                                                                  C
C****               IU        LENGTH OF THE INPUT- AND OUTPUTVECTOR    C
C****                         U,SO IU EQUALS M*N                       C
C****                                                                  C
C****               IV        LENGTH OF THE OUTPUTVECTOR V,SO IV       C
C****                         EQUALS N*N                               C
C****                                                                  C
C****               ICV       OUTPUT PARAMETER,CONTAINS THE MAXIMUM    C
C****                         NUMBER OF ITERATIONS USED,IF ICV = 0     C
C****                         THERE WAS NO CONVERGENCE FOR ONE OF THE  C
C****                         SINGULAR VALUES AND THE SUBROUTINE       C
C****                         HAS STOPPED                              C
C****                                                                  C
C****                                                                  C
C**** DIMENSIONS:   U(IU)     CONTAINS ON INPUT THE M*N MATRIX TO BE   C
C****                         COMPOSED,ON OUTPUT IT CONTAINS THE M*N   C
C****                         ORTHONORMALIZED MATRIX U                 C
C****                                                                  C
C****               V(IV)     REPRESENTS THE ORTHOGONAL MATRIX V       C
C****                                                                  C
C****               Q(N)      A VECTOR HOLDING THE NONNEGATIVE         C
C****                         SINGULAR VALUES OF A,IN DECREASING       C
C****                         SEQUENCE                                 C
C****                                                                  C
C****               E(N)      A VECTOR USED IN THE SUBROUTINE,IT       C
C****                         SHOULD BE DECLARED IN THE MAIN PROGRAM   C
C****                                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C****
C**** * * * * * *   START OF PROGRAM   * * * * * * * * * * * *
C****
C****
      DIMENSION U(IU),V(IV),Q(N),E(N)
      DATA ZERO,ONE,TWO/0.0,1.0,2.0/
      EPS=2.0**(-24)
      TOL=10.0**(-38)
C****
C****
C**** * * * * * *   HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM
C****
C****
      G=ZERO
      X=ZERO
      DO 130 I=1,N
      E(I)=G
      S=ZERO
      L=I+1
      DO 20 J=I,M
      II=J+(I-1)*M
      S=S+U(II)*U(II)
   20 CONTINUE
      IF(S.GE.TOL) GOTO 30
      G=ZERO
      GOTO 60
   30 CONTINUE
      JJ=I+(I-1)*M
      F=U(JJ)
      G=-SQRT(S)
      IF(F.LT.ZERO) G=-G
      H=F*G-S
      U(JJ)=F-G
      IF(L.GT.N) GOTO 60
      DO 50 J=L,N
      S=ZERO
      DO 40 K=I,M
      S=S+U(K+(I-1)*M)*U(K+(J-1)*M)
   40 CONTINUE
      F=S/H
      DO 50 K=I,M
      KJ=K+(J-1)*M
      U(KJ)=U(KJ)+F*U(K+(I-1)*M)
   50 CONTINUE
   60 CONTINUE
      Q(I)=G
      S=ZERO
      IF(L.GT.N) GOTO 75
      DO 70 J=L,N
      IJ=I+(J-1)*M
      S=S+U(IJ)*U(IJ)
   70 CONTINUE
   75 CONTINUE
      IF(S.GE.TOL) GOTO 80
      G=ZERO
      GOTO 120
   80 CONTINUE
      I1=(M+1)*I
      F=U(I1)
      G=-SQRT(S)
      IF(F.LT.ZERO) G=-G
      H=F*G-S
      U(I1)=F-G
      DO 90 J=L,N
      E(J)=U(I+(J-1)*M)/H
   90 CONTINUE
      DO 110 J=L,M
      S=ZERO
      DO 100 K=L,N
      S=S+U(J+(K-1)*M)*U(I+(K-1)*M)
  100 CONTINUE
      DO 110 K=L,N
      JK=J+(K-1)*M
      U(JK)=U(JK)+S*E(K)
  110 CONTINUE
  120 CONTINUE
      Y=ABS(Q(I))+ABS(E(I))
      IF(Y.GT.X) X=Y
  130 CONTINUE
C****
C****
C**** * * * * * *   ACCUMULATION OF RIGHT HAND TRANSFORMATIONS
C****
C****
      DO 210 II=1,N
      I=(N+1)-II
      IF(G.NE.ZERO) GOTO 150
      GOTO 190
  150 CONTINUE
      IPLUS=(M+1)*I
      H=U(IPLUS)*G
      DO 160 J=L,N
      V(J+(I-1)*N)=U(I+(J-1)*M)/H
  160 CONTINUE
      DO 180 J=L,N
      S=ZERO
      DO 170 K=L,N
      S=S+U(I+(K-1)*M)*V(K+(J-1)*N)
  170 CONTINUE
      DO 180 K=L,N
      KJ=K+(J-1)*N
      V(KJ)=V(KJ)+S*V(K+(I-1)*N)
  180 CONTINUE
  190 CONTINUE
      IF(L.GT.N) GOTO 205
      DO 200 J=L,N
      V(J+(I-1)*N)=ZERO
      V(I+(J-1)*N)=ZERO
  200 CONTINUE
  205 CONTINUE
      V(I+(I-1)*N)=ONE
      G=E(I)
      L=I
  210 CONTINUE
C****
C****
C**** * * * * * *   ACCUMULATION OF LEFT HAND TRANSFORMATIONS
C****
C****
      DO 320 II=1,N
      I=(N+1)-II
      JJ=I+(I-1)*M
      L=I+1
      G=Q(I)
      IF(L.GT.N) GOTO 245
      DO 240 J=L,N
      U(I+(J-1)*M)=ZERO
  240 CONTINUE
  245 CONTINUE
      IF(G.NE.ZERO) GOTO 250
      GOTO 290
  250 CONTINUE
      H=U(JJ)*G
      IF(L.GT.N) GOTO 275
      DO 270 J=L,N
      S=ZERO
      DO 260 K=L,M
      S=S+U(K+(I-1)*M)*U(K+(J-1)*M)
  260 CONTINUE
      F=S/H
      DO 270 K=I,M
      KJ=K+(J-1)*M
      U(KJ)=U(KJ)+F*U(K+(I-1)*M)
  270 CONTINUE
  275 CONTINUE
      DO 280 J=I,M
      JI=J+(I-1)*M
      U(JI)=U(JI)/G
  280 CONTINUE
      GOTO 310
  290 CONTINUE
      DO 300 J=I,M
      U(J+(I-1)*M)=ZERO
  300 CONTINUE
  310 CONTINUE
      U(JJ)=U(JJ)+ONE
  320 CONTINUE
C****
C****
C**** * * * * * *   DIAGONALIZATION OF THE BIDIAGONAL FORM   * *
C****
C****
      EPS=EPS*X
      ICV=0
      ICONV=1
      DO 470 KK=1,N
      K=(N+1)-KK
      IF(ICV.LT.ICONV) ICV=ICONV
      ICONV=0
  340 CONTINUE
      ICONV=ICONV+1
      IF(ICONV.GT.50) GOTO 550
      DO 350 LL=1,K
      L=(K+1)-LL
      IF(ABS(E(L)).LE.EPS) GOTO 390
      LMIN=L-1
      IF(ABS(Q(LMIN)).LE.EPS) GOTO 360
  350 CONTINUE
  360 CONTINUE
C****
C****
C**** * * * * * *   CANCELLATION OF E(L) IF L>1   * * * * * * *
C****
C****
      C=ZERO
      S=ONE
      L1=L-1
      DO 380 I=L,K
      F=S*E(I)
      E(I)=C*E(I)
      IF(ABS(F).LE.EPS) GOTO 390
      G=Q(I)
      IF(ABS(F).LT.ABS(G)) GOTO 362
      IF (F.EQ.ZERO) GOTO 365
      H=G/F
      H=ABS(F)*SQRT(ONE+H*H)
      GOTO 368
  362 CONTINUE
      H=F/G
      H=ABS(G)*SQRT(ONE+H*H)
      GOTO 368
  365 CONTINUE
      H=ZERO
      GOTO 369
  368 C=G/H
      S=-F/H
  369 Q(I)=H
      DO 370 J=1,M
      JL=J+(L1-1)*M
      JI=J+(I-1)*M
      Y=U(JL)
      Z=U(JI)
      U(JL)=Y*C+Z*S
      U(JI)=-Y*S+Z*C
  370 CONTINUE
  380 CONTINUE
  390 CONTINUE
C****
C****
C**** * * * * * *   TEST CONVERGENCE   * * * * * * * * * * * *
C****
C****
      Z=Q(K)
      IF(L.EQ.K) GOTO 450
C****
C****
C**** * * * * * *   SHIFT FROM BOTTOM 2 * 2 MINOR   * * * * * *
C****
C****
      X=Q(L)
      KMIN=K-1
      Y=Q(KMIN)
      G=E(KMIN)
      H=E(K)
      F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)
      G=SQRT(F*F+ONE)
      R=F+G
      IF(F.LT.ZERO) R=F-G
      F=((X-Z)*(X+Z)+H*(Y/R-H))/X
C****
C****
C**** * * * * * *   NEXT QR TRANSFORMATION   * * * * * * * * *
C****
C****
      C=ONE
      S=ONE
      L2=L+1
      DO 430 I=L2,K
      G=E(I)
      Y=Q(I)
      H=S*G
      G=C*G
      IMIN=I-1
      IF(ABS(F).LT.ABS(H)) GOTO 392
      IF(F.EQ.ZERO) GOTO 395
      Z=H/F
      Z=ABS(F)*SQRT(ONE+Z*Z)
      GOTO 398
  392 CONTINUE
      Z=F/H
      Z=ABS(H)*SQRT(ONE+Z*Z)
      GOTO 398
  395 CONTINUE
      Z=ZERO
      GOTO 399
  398 C=F/Z
      S=H/Z
  399 E(IMIN)=Z
      F=X*C+G*S
      G=-X*S+G*C
      H=Y*S
      Y=Y*C
      DO 400 J=1,N
      JM=J+(IMIN-1)*N
      JI=J+(I-1)*N
      X=V(JM)
      Z=V(JI)
      V(JM)=X*C+Z*S
      V(JI)=-X*S+Z*C
  400 CONTINUE
      IF(ABS(F).LT.ABS(H)) GOTO 412
      IF(F.EQ.ZERO) GOTO 415
      Z=H/F
      Z=ABS(F)*SQRT(ONE+Z*Z)
      GOTO 418
  412 CONTINUE
      Z=F/H
      Z=ABS(H)*SQRT(ONE+Z*Z)
      GOTO 418
  415 CONTINUE
      Z=ZERO
      GOTO 419
  418 C=F/Z
      S=H/Z
  419 Q(IMIN)=Z
      F=C*G+S*Y
      X=-S*G+C*Y
      DO 420 J=1,M
      JM=J+(IMIN-1)*M
      JI=J+(I-1)*M
      Y=U(JM)
      Z=U(JI)
      U(JM)=Y*C+Z*S
      U(JI)=-Y*S+Z*C
  420 CONTINUE
  430 CONTINUE
      E(L)=ZERO
      E(K)=F
      Q(K)=X
      GOTO 340
C****
C****
C**** * * * * * *   CONVERGENCE   * * * * * * * * * * * * * * *
C****
C****
  450 CONTINUE
      IF(Z.GE.ZERO) GOTO 470
C****
C****
C**** * * * * * *   Q(K) IS MADE NON-NEGATIVE   * * * * * * * *
C****
C****
      Q(K)=-Z
      DO 460 J=1,N
      JK=J+(K-1)*N
      V(JK)=-V(JK)
  460 CONTINUE
  470 CONTINUE
C****
C****
C**** * * * * * *   ORDERING OF SINGULAR VALUES   * * * * * * *
C****
C****
      NMIN=N-1
      DO 540 I=1,NMIN
      IPLUS=I+1
  480 CONTINUE
      DO 530 J=IPLUS,N
      IF(Q(I).GE.Q(J)) GOTO 530
      R=Q(I)
      Q(I)=Q(J)
      Q(J)=R
      DO 490 K=1,N
      KI=K+(I-1)*N
      KJ=K+(J-1)*N
      R=V(KI)
      V(KI)=V(KJ)
      V(KJ)=R
  490 CONTINUE
      DO 510 K=1,M
      KI=K+(I-1)*M
      KJ=K+(J-1)*M
      R=U(KI)
      U(KI)=U(KJ)
      U(KJ)=R
  510 CONTINUE
      JPLUS=J+1
      IF(JPLUS.GT.N) GOTO 540
      GOTO 480
  530 CONTINUE
  540 CONTINUE
      GOTO 570
  550 CONTINUE
      ICV=0
  570 CONTINUE
      RETURN
      END
      SUBROUTINE STATIS(OPTION,NROW,NCOL,NDIM,FITT,NANA,K,IWRITE)
C     ******************************************************************
C     *                                                                *
C     * S T A T I S                                                    *
C     *                                                                *
C     *  PURPOSE: COMPUTE STATISTICS                                   *
C     *                                                                *
C     * SUBROUTINES CALLED: SHEL9                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION OPTION(NANA),IORD(4),HELP(4),FITT(NANA)
      INTEGER OPTION
      WRITE(IWRITE,640) K
      DO 5 J=1,NANA
C
C WARNING
C
         IF (OPTION(J).GT.4)  WRITE (IWRITE,660)
         IF (OPTION(J).GT.4) GO TO 6
    5 CONTINUE
    6 DO 10 J=1,NANA
         FITT(J)=FITT(J)**2
         HELP(J)=FLOAT(MOD(OPTION(J),10))
         IORD(J)=J
   10 CONTINUE
      CALL SHEL9(HELP,IORD,NANA)
      DO 15 J=1,NANA
         R=HELP(J)
         IF (R.EQ.0.) GO TO 15
         IF (FITT(IORD(J)).GT. .9995) WRITE(IWRITE,650) IORD(J)
         IF (FITT(IORD(J)).GT. .9995) GO TO 15
         F1=FITT(IORD(J))/(1.-FITT(IORD(J)))
         ID=NCOL-1
         IF (J.EQ.NANA) RR=5
         IF (J.EQ.NANA) GO TO 16
         RR=HELP(J+1)
         IF (RR.EQ.R) GO TO 15
         F2=(FITT(IORD(J))-FITT(IORD(J+1)))/(1.-FITT(IORD(J)))
   16    IF (R.EQ.4.) GO TO 40
         IF (R.EQ.3.) GO TO 30
         IF (R.EQ.2.) GO TO 20
         IDF1=IFIX((FLOAT(NDIM)/2)*(NDIM+3))
         IF (RR.EQ.2.) IDF12=.5*NDIM*(NDIM-1)
         IF (RR.EQ.3.) IDF12=.5*NDIM*(NDIM+1)-1
         IF (RR.EQ.4.) IDF12=.5*NDIM*(NDIM+1)
         GO TO 50
   20    IDF1=2*NDIM
         IF (RR.EQ.3.) IDF12=NDIM-1
         IF (RR.EQ.4.) IDF12=NDIM
         GO TO 50
   30    IDF1=NDIM+1
         IF (RR.EQ.4.) IDF12=1
         GO TO 50
   40    IDF1=NDIM
         IDF2=ID-IDF1
         GO TO 60
   50    IDF2=ID-IDF1
         IF (RR.EQ.0.) GO TO 60
         IF (RR.EQ.5.) GO TO 60
         FF2=(FLOAT(IDF2)/FLOAT(IDF12))*F2
   60    FF1=(FLOAT(IDF2)/FLOAT(IDF1))*F1
         FIT=FITT(IORD(J))
         IF (R.EQ.1. .AND.RR.NE.1.) WRITE (IWRITE,600) FIT,FF1,IDF1,IDF2
         IF (R.EQ.1. .AND.RR.EQ.2.) WRITE (IWRITE,601) FF2,IDF12,IDF2
         IF (R.EQ.1. .AND.RR.EQ.3.) WRITE (IWRITE,602) FF2,IDF12,IDF2
         IF (R.EQ.1. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2
         IF (R.EQ.2. .AND.RR.NE.2.) WRITE (IWRITE,604) FIT,FF1,IDF1,IDF2
         IF (R.EQ.2. .AND.RR.EQ.3.) WRITE (IWRITE,602) FF2,IDF12,IDF2
         IF (R.EQ.2. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2
         IF (R.EQ.3. .AND.RR.NE.3.) WRITE (IWRITE,605) FIT,FF1,IDF1,IDF2
         IF (R.EQ.3. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2
         IF (R.EQ.4. .AND.RR.NE.4.) WRITE (IWRITE,606) FIT,FF1,IDF1,IDF2
   15 CONTINUE
      WRITE(IWRITE,607)
C
  600 FORMAT(1H0,12X,34HGENERAL UNFOLDING MODEL,  P.VAF. =,F6.3,5H, F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  601 FORMAT(1H ,15X,42HVERSUS WEIGHTED UNFOLDING MODEL        F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  602 FORMAT(1H ,15X,42HVERSUS UNFOLDING MODEL                 F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  603 FORMAT(1H ,15X,42HVERSUS VECTOR MODEL                    F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  604 FORMAT(1H0,12X,34HWEIGHTED UNFOLDING MODEL, P.VAF. =,F6.3,5H, F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  605 FORMAT(1H0,12X,34HUNFOLDING MODEL,          P.VAF. =,F6.3,5H, F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  606 FORMAT(1H0,12X,34HVECTOR MODEL,             P.VAF. =,F6.3,5H, F =,
     X        F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM)
  607 FORMAT(1H0//)
  640 FORMAT(1H0//,13X,33HSTATISTICS ACROSS OPTIONS FOR ROW,I5/
     X       1H ,12X,38H--------------------------------------)
  650 FORMAT(1H0/,13X,22HPERFECT FIT FOR OPTION,I2,16H. NO STATISTICS.)
  660 FORMAT(1H ,12X,38HTHE STATISTICS GIVEN BELOW DEPEND ONLY,
     X               46H ON THE METRIC FIT, NOT ON THE NON-METRIC FIT.)
      RETURN
      END
      SUBROUTINE TEST(P1,N1,N2)
C     ******************************************************************
C     *                                                                *
C     *  T E S T                                                       *
C     *                                                                *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION P1(N1,N2)
C
      DO 10 I=1,N1
         WRITE(6,600) (P1(I,J),J=1,N2)
   10 CONTINUE
      WRITE(6,601) N1,N2
  600 FORMAT(10F11.4)
  601 FORMAT(10X,2I10)
      RETURN
      END
      SUBROUTINE TIEBL (DISM, LBK, N, IDLBK )
C     ******************************************************************
C     *                                                                *
C     *  T I E B L                                                     *
C     *                                                                *
C     *  PURPOSE:- FIND  A VECTOR LBK SO THAT EACH ELEMENT OF IT       *
C     *            REPRESENTS THE LENGTH OF A TIEBLOCK FOUND IN DISM.  *
C     *          - IDLBK GIVES THE NUMBER OF TIEBLOCKS                 *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE.                                     *
C     *                                                                *
C     *  SEE ACCOMPANYING PAPER SECTIONS: 2.2.                         *
C     *                                                                *
C     *  AUTHOR: ERNST VAN WANING.     RELEASED: JULY 1976             *
C     *                                                                *
C     ******************************************************************
C
C
      DIMENSION DISM(N),LBK(N)
C
      IDLBK = 0
      JBASE = 1
    3    K = 1
    2    IDDISM = JBASE + K
         IF (IDDISM .GT. N) GOTO 1
         IF (DISM(JBASE) .NE. DISM(IDDISM)) GOTO 1
            K = K + 1
            GOTO 2
    1    IDLBK = IDLBK + 1
         LBK(IDLBK) = K
         JBASE = IDDISM
      IF (JBASE .LE. N) GOTO 3
      RETURN
      END
      SUBROUTINE TRED2(NM,N,D,E,Z)
C     ******************************************************************
C     *                                                                *
C     *  T R E D 2                                                     *
C     *                                                                *
C     *  PURPOSE: REDUCES A REAL SYMMETRIC MATRIX TO A SYMMETRIC TRI-  *
C     *           DIAGONAL MATRIX USING ORTHOGONAL SIMILARITY TRANS-   *
C     *           FORMATIONS; THE ACCUMULATED ORTHOGONAL TRANSFORMA-   *
C     *           TIONS ARE RETAINED IN Z.                             *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     *  ADAPTED FROM EISPACK                                          *
C     *                                                                *
C     ******************************************************************
      DIMENSION Z(NM,N),D(N),E(N)
C
      IF (N .EQ. 1) GO TO 320
C     ********** FOR I=N STEP -1 UNTIL 2 DO -- **********
      DO 300 II =2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0
         SCALE = 0.0
         IF (L .LT.2) GO TO 130
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
         DO 120 K =1, L
  120    SCALE = SCALE + ABS(Z(I,K))
C
         IF (SCALE .NE. 0.0) GO TO 140
  130    E(I) = Z(I,L)
         GO TO 290
C
  140    DO 150 K = 1, L
            Z(I,K) = Z(I,K) / SCALE
            H = H + Z(I,K) * Z(I,K)
  150    CONTINUE
C
         F = Z(I,L)
         G = -SIGN(SQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         Z(I,L) = F - G
         F = 0.0
C
         DO 240 J = 1, L
            Z(J,I) = Z(I,J) / (SCALE * H)
            G = 0.0
C     ********** FORM ELEMENT OF A*U **********
            DO 180 K = 1, J
  180       G = G + Z(J,K) * Z(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + Z(K,J) * Z(I,K)
C     ********** FORM ELEMENT OF P **********
  220       E(J) = G / H
            F = F + E(J) * Z(I,J)
  240    CONTINUE
C
         HH = F / (H + H)
C     ********** FORM REDUCED A **********
         DO 260 J =1, L
            F = Z(I,J)
            G = E(J) - HH * F
            E(J) = G
C
            DO 260 K = 1, J
               Z(J,K) = Z(J,K) - F * E(K) -G * Z(I,K)
  260    CONTINUE
C
         DO 280 K = 1, L
  280    Z(I,K) = SCALE * Z(I,K)
C
  290    D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0
      E(1) = 0.0
C     ********** ACCUMULATION OF TRANSFORMATION MATRICES **********
      DO 500 I = 1, N
         L = I - 1
         IF (D(I) .EQ. 0.0) GO TO 380
C
C
         DO 360 J = 1, L
            G = 0.0
            DO 340 K =1, L
  340       G = G + Z(I,K) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * Z(K,I)
  360    CONTINUE
C
  380    D(I) = Z(I,I)
         Z(I,I) = 1.0
         IF (L .LT. 1) GO TO 500
C
         DO 400 J =1, L
            Z(I,J) = 0.0
            Z(J,I) = 0.0
  400    CONTINUE
C
  500 CONTINUE
C
      RETURN
      END
      SUBROUTINE UNRAV(COEF,NPRM,COOR,WEIGH,ROTAT,OPTION,NDIM,ND,JJ,NR,
     X IWRITE,WSIGN,IPRE,WORK1,WORK3,WORK2,WORK4,CMAX,IPRI1,SLOPE,CEPT,
     X FIT)
C     ******************************************************************
C     *                                                                *
C     *  U N R A V E L                                                 *
C     *                                                                *
C     *  PURPOSE: UNRAVELS THE REGRESSION COEFFICIENTS AND COMPUTES    *
C     *           THE COORDINATES FOR THE IDEAL POINTS AND/OR VECTORS  *
C     *                                                                *
C     *  SUBROUTINES CALLED: UNRA1,UNRA2,UNRA3,UNRA4,RESUL1            *
C     *                                                                *
C     ******************************************************************
      DIMENSION COEF(NPRM),COOR(ND),WEIGH(ND),ROTAT(NDIM,ND),
     X          WSIGN(NDIM),WORK1(NDIM),WORK4(NDIM),WORK2(NDIM,NDIM),
     X          WORK3(NDIM,NDIM)
      INTEGER OPTION,OPT
C
      OPT=MOD(OPTION,10)
      IF (OPT.EQ.1) CALL UNRA1(COEF,NPRM,COOR,WEIGH,ROTAT,NDIM,ND,JJ,
     X                         WORK1,WORK3,WORK2,WORK4,SLOPE,CEPT)
      IF (OPT.EQ.2) CALL UNRA2(COEF,NPRM,COOR,WEIGH,NDIM,ND,JJ,
     X                         SLOPE,CEPT)
      IF (OPT.EQ.3) CALL UNRA3(COEF,NPRM,COOR,WEIGH,WSIGN,NDIM,ND,JJ,
     X                         SLOPE,CEPT)
      IF (OPT.EQ.4) CALL UNRA4(COEF,NPRM,COOR,WEIGH,NDIM,ND,JJ,CMAX,
     X                         SLOPE,CEPT,FIT)
      IF (OPT.NE.0.AND.IPRI1.EQ.1)
     X    CALL RESUL1(COOR,WEIGH,ROTAT,NDIM,ND,JJ,IWRITE,OPT,NR,IPRE)
C
      RETURN
      END
      SUBROUTINE UNRA1(COEF,NPR,COOR,WEIGH,ROTAT,NDIM,ND,NOP,
     X                 WORK1,LSV,RSV,SIVA,STAND,CEPT)
C     ******************************************************************
C     *                                                                *
C     *  U N R A 1                                                     *
C     *                                                                *
C     *  PURPOSE: UNRAVELS FOR THE GENERAL UNFOLDING MODEL             *
C     *                                                                *
C     *  SUBROUTINES CALLED: TRED2,IMTQL2                              *
C     *                                                                *
C     ******************************************************************
      DIMENSION COEF(NPR),COOR(ND),LSV(NDIM,NDIM),RSV(NDIM,NDIM),
     X          SIVA(NDIM),WORK1(NDIM),WEIGH(ND),ROTAT(NDIM,ND)
      REAL LSV
C
C
      NOPT=NOP-1
C
C     NOPT WILL INDICATE ON WHAT POSITIONS ARRAY COOR HAS TO BE FILLED
C
C
      N=NDIM+2
      K=2*NDIM+2
      L=2
      M=NDIM-1
      DO 1 I=1,M
         DO 3 J=L,NDIM
            RSV(I,J)=.5*COEF(K)
            RSV(J,I)=RSV(I,J)
            K=K+1
    3    CONTINUE
         L=L+1
    1 CONTINUE
      DO 5 I=1,NDIM
         RSV(I,I)=COEF(N)
         N=N+1
    5 CONTINUE
C
C
C   COMPUTE THE EIGENVALUE DECOMPOSITION OF THE TRANSFORMATION MATRIX
C   RSV.
C   ON THE BASIS THEREOF COMPUTE (IN LSV) THE INVERSE (FOR THE IDEAL
C   POINTS); THE EIGENVALUES (IN SIVA) WILL GIVE THE WEIGHTS;
C   THE EIGENVECTORS (IN RSV) WILL GIVE THE ROTATION.
C
      CALL TRED2(NDIM,NDIM,SIVA,WORK1,RSV)
      CALL IMTQL2(NDIM,NDIM,SIVA,WORK1,RSV,IERR)
C
C COLLECT WEIGHTS FROM ARRAY SIVA AND ROTATION FROM ARRAY RSV
C STANDARDIZE WEIGHTS AND COMPUTE SLOPE
C
      STAND=0.0
      DO 64 I=1,NDIM
         STAND=STAND+SIVA(I)*SIVA(I)
   64 CONTINUE
      STAND=SQRT(STAND/NDIM)
      DO 65 I=1,NDIM
         IL=NOPT*NDIM+I
         WEIGH(IL)=SIVA(I)/STAND
   65 CONTINUE
      DO 70 I=1,NDIM
         DO 80 J=1,NDIM
            IL=NOPT*NDIM+J
            ROTAT(I,IL)=RSV(I,J)
   80    CONTINUE
   70 CONTINUE
C
      DO 6 I=1,NDIM
         SMALL=.0000001
         IF (ABS(SIVA(I)).GT.SMALL) SIVA(I)=1.0/SIVA(I)
         IF (ABS(SIVA(I)).LE.SMALL) SIVA(I)=0.0
    6 CONTINUE
C
      DO 10 K=1,NDIM
         DO 20 I=1,NDIM
            SUM=0.0
            DO 30 J=1,NDIM
               SUM=SUM+RSV(I,J)*RSV(K,J)*SIVA(J)
   30       CONTINUE
            WORK1(I)=SUM
   20    CONTINUE
         DO 40 I=1,NDIM
            LSV(K,I)=WORK1(I)
   40    CONTINUE
   10 CONTINUE
C
C
C  COMPUTE THE IDEAL POINTS AND INTERCEPT
C
      CEPT=0.0
      DO 50 I=1,NDIM
         IL=NOPT*NDIM+I
         SUM=0.0
         DO 60 J=1,NDIM
            L=J+1
            SUM=SUM+COEF(L)*LSV(J,I)
   60    CONTINUE
         L=I+1
         COOR(IL)=-.5*SUM
         ASIGN=SIGN(1.,COOR(IL))
         IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*ASIGN
         CEPT=CEPT+COOR(IL)*COEF(L)
   50 CONTINUE
      CEPT=COEF(1)+.5*CEPT
C
      RETURN
      END
      SUBROUTINE UNRA2(COEF,NPR,COOR,WEIGH,NDIM,ND,NOP,STAND,CEPT)
C     ******************************************************************
C     *                                                                *
C     *  U N R A 2                                                     *
C     *                                                                *
C     *  PURPOSE: UNRAVELS FOR THE WEIGHTED UNFOLDING MODEL            *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION COEF(NPR),COOR(ND),WEIGH(ND)
C
      NOPT=NOP-1
C
C     FOR COMMENT SEE UNRA1
C
C COMPUTE SLOPE AND STANDARDIZE WEIGHTS
C
      STAND=0.0
      DO 5 I=1,NDIM
         K=I+NDIM+1
         STAND=STAND+COEF(K)*COEF(K)
    5 CONTINUE
      STAND=SQRT(STAND/NDIM)
      CEPT=0.0
      DO 10 I=1,NDIM
         IL=NOPT*NDIM+I
         L=I+1
         K=L+NDIM
         WEIGH(IL)=COEF(K)/STAND
         COOR(IL)=-.5*COEF(L)/COEF(K)
         ASIGN=SIGN(1.,COOR(IL))
         IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*ASIGN
         CEPT=CEPT+COOR(IL)*COEF(L)
   10 CONTINUE
      CEPT=COEF(1)+.5*CEPT
      RETURN
      END
      SUBROUTINE UNRA3(COEF,NPR,COOR,WEIGH,WSIGN,NDIM,ND,NOP,STAND,CEPT)
C     ******************************************************************
C     *                                                                *
C     *  U N R A 3                                                     *
C     *                                                                *
C     *  PURPOSE: UNRAVELS FOR THE UNFOLDING MODEL                     *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION COEF(NPR),COOR(ND),WEIGH(ND),WSIGN(NDIM)
C
C
      NOPT=NOP-1
C
C     SEE UNRA1 FOR COMMENT
C
      K=NDIM+2
      STAND=ABS(COEF(K))
      CEPT=0.0
      ASIGN=1.
      DO 10 I=1,NDIM
         L=I+1
         IL=NOPT*NDIM+I
         COOR(IL)=-.5*COEF(L)/COEF(K)*WSIGN(I)
         CEPT=CEPT+COOR(IL)*COEF(L)
         IF (COEF(K).LT.0.) ASIGN=-1.
         WEIGH(IL)=WSIGN(I)*ASIGN
         BSIGN=SIGN(1.,COOR(IL))
         IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*BSIGN
   10 CONTINUE
      CEPT=COEF(1)+.5*CEPT
      RETURN
      END
      SUBROUTINE UNRA4(COEF,NPR,COOR,WEIGH,NDIM,ND,NOP,CMAX,STAND,CEPT,
     X                 FIT)
C     ******************************************************************
C     *                                                                *
C     *  U N R A 4                                                     *
C     *                                                                *
C     *  PURPOSE: UNRAVELS FOR THE VECTOR MODEL                        *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION COEF(NPR),COOR(ND),WEIGH(ND)
C
      NOPT=NOP-1
C
C
C     SEE UNRA1 FOR COMMENT
C
      SUM=0.0
      K=NDIM+1
      DO 10 J=2,K
         SUM=SUM+COEF(J)*COEF(J)
   10 CONTINUE
      STAND=SQRT(SUM)
      DO 20 I=1,NDIM
         IL=NOPT*NDIM+I
         L=I+1
         COOR(IL)=FIT*CMAX*COEF(L)/STAND
         WEIGH(IL)=1.
   20 CONTINUE
      CEPT=COEF(1)
      STAND=STAND/(FIT*CMAX)
      RETURN
      END
      SUBROUTINE ZERO(A,L)
C     ******************************************************************
C     *                                                                *
C     *  Z E R O                                                       *
C     *                                                                *
C     *  PURPOSE: FILLS ARRAY A WITH ZERO'S                            *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION A(L)
C
      DO 10 J=1,L
         A(J)=0.0
   10 CONTINUE
      RETURN
      END
      SUBROUTINE ZEROI(N,L)
C     ******************************************************************
C     *                                                                *
C     *  Z E R O I                                                     *
C     *                                                                *
C     *  PURPOSE: FILLS ARRAY N WITH ZERO'S                            *
C     *                                                                *
C     *  SUBROUTINES CALLED: NONE                                      *
C     *                                                                *
C     ******************************************************************
      DIMENSION N(L)
C
      DO 10 J=1,L
         N(J)=0
   10 CONTINUE
      RETURN
      END

