C                                            ***** STAT PACK *****
C     SUBROUTINE FOR PERFORMING DISCRIMINATE ANALYSIS
C     CALLING SEQUENCE: CALL DISCR(NV,NC,MV,MC,DATA,AV,JV,NAMES)
C     WHERE: NV - NUMBER OF VARIABLES
C            NC - NUMBER OF OBSERVATIONS
C            MV - MAXIMUM NUMBER OF VARIABLES POSSIBLE
C            MC - MAXIMUM NUMBER OF OBSERVATIONS POSSIBLE
C            DATA - MATRIX CONTAINING DATA
C            AV - EXTRA VECTOR AT LEAST NC LONG
C            JV - EXTRA VECTOR AT LEAST NC LONG
C            NAMES - VECTOR CONTAINING NAMES OF VARIABLES
C
C     DISCRIMINATE ANALYSIS WAS IN ORIGINAL REQUEST FOR ROUTINES
C     WHICH PROMPTED THE CREATION OF STAT PACK, SUBMITTED
C     BY AL LEADER (MANAGEMENT).  MORE RECENTLY SAM ANEMA
C     (COMPUTER CENTER FACULTY,STAFF COORDINATOR) EXPRESSED
C     AN OPINION THAT STAT PACK SHOULD CONTAIN A DISCRIMANNTE
C     ANALYSIS.  IN DISCUSSIONS WITH ULDIS SMMIDCHENS THE NEED
C     WAS ALSO INDICATED.  THIS ROUTINE IS ADOPTED FROM A PROGRAM
C     ORIGINALLY WRITTEN BY SAM ANEMA.
C
      SUBROUTINE DISCR(NV,NC,MV,MC,DATA,AV,JV,NAMES)
      DIMENSION XBAR(25,25),D(25,25),C(25,25),IP(25,25)
      DIMENSION OPT(8),IVAR(25),LOPT(8),RANGE(3,25),INPUT(50)
      DIMENSION IVALUE(15),SPACE(3),IV(25),XC(25),CON(25),XVAL(25)
      DIMENSION JV(1),AV(1),NAMES(1),DATA(MC,MV),NGRP(25)
      EQUIVALENCE (XBAR,IP),(INPUT,IV)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON /EXTRA/ HEDR(70),NSZ
      DATA LOPT/'DISCR','RANGE','MEANS','DISPM','INVDM','EVALF',
     1'XPROD','CLASS'/
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT('0ENTER OPTIONS SEPARATED BY COMMAS'/)
      READ(ICC,3,END=1)(IVAR(I),I=1,10)
3     FORMAT(10(A5,1X))
      IF(IVAR(1).EQ.'!') RETURN
      DO 13 I=1,8
13    OPT(I)=0
      DO 12 I=1,10
      IF(IVAR(I).EQ.' ') GO TO 20
      IF(IVAR(I).EQ.'HELP') GO TO 4
      IF(IVAR(I).NE.'HELP,') GO TO 6
4     WRITE(IDLG,5)
5     FORMAT('0THE DISCRIMINATE ANALYSIS MAY BE USED FOR UP TO'/
     1' 25 VARIABLES WITH UP TO 25 GROUPS.  PLACEMENT OF OBSERVATIONS'/
     2' INTO GROUPS IS BASED ON THE VALUE OF A BREAKDOWN VARIABLE.'/
     3' RANGES FOR GROUPING OBSERVATIONS ARE ASSUMED TO BE ENTERED BY'/
     4' THE USER, HOWEVER, OPTIONS ARE AVAILABLE TO FORM RANGES'/
     5' AUTOMATICALLY.  IF THE BREAKDOWN VARIABLE HAS 25 OR FEWER'/
     6' DISCRETE VALUES, A SEPARATE RANGE WILL BE CREATED FOR EACH'/
     7' INDIVIDUAL VALUE, OTHERWISE THE RANGE OF VALUES FOR THE '/
     8' BREAKDOWN VARIABLE WILL BE BROKEN INTO 25 SEPARATE AND EQUAL'/
     9' RANGES.'/
     1'0OPTIONS AVAILABLE ARE:'/
     2' "DISCR" - AUTOMATICALLY FORM GROUPINGS'/
     3' "RANGE" - LIST RANGES USED IN BREAKDOWNS (AVAILABLE ONLY IF'/
     4'           "DISCR" IS USED.)'/
     5' "AUTO"  - SAME AS "DISCR" EXCEPT SPECIFIED WHEN RANGES ARE'/
     6'           REQUESTED.'/
     7' "MEANS" - PRINT A MEANS TABLE FOR VARIABLE GROUP COMBINATIONS'/
     8' "XPROD" - PRINT THE CROSS PRODUCTS OF VARIANCE FROM MEAN FOR'/
     9'           EACH GROUP')
      WRITE(IDLG,16)
16    FORMAT(' "DISPM" - PRINT THE DISPERSION MATRIX'/
     1' "INVDM" - PRINT THE INVERSE OF THE DISPERSION MATRIX'/
     2' "EVALF" - EVALUATE THE DISCRIMINATE FUNCTION FOR EACH OBSER.'/
     4' "CLASS" - PRINT THE CLASSIFICATION MATIRX'/
     5' "OUTPT" - "MEANS", "DISPM", "XPROD", "INVDM", "EVALF", AND'/
     6'           "CLASS"'/
     7' IF NO OPTIONS ARE DESIRED TYPE A CARRIAGE RETURN')
      GO TO 1
6     IF(IVAR(I).EQ.'AUTO') GO TO 7
      IF(IVAR(I).NE.'AUTO,') GO TO 9
7     WRITE(IDLG,8)
8     FORMAT(' "AUTO" IS ONLY SPECIFIED WHEN RANGES ARE REQUESTED')
      GO TO 1
9     IOPT=0
      DO 10 J=1,8
      IF(IVAR(I).NE.LOPT(J)) GO TO 10
      OPT(J)=1
      GO TO 12
10    CONTINUE
      IF(IVAR(I).NE.'OUTPT') GO TO 15
      DO 14 J=3,8
14    OPT(J)=1
      GO TO 12
15    WRITE(IDLG,11) IVAR(I)
11    FORMAT(' OPTION "',A5,'" DOES NOT EXIST')
      GO TO 1
12    CONTINUE
C
C     OPTIONS TAKEN CRE OF NOW FIND VARIABLES TO BE ANALYSED
C
20    IF(ICC.NE.2) WRITE(IDLG,21)
21    FORMAT(' WHICH VARIABLES ARE TO BE ANALYSED?'/)
      CALL ALPHA(IVAR,25,NN,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 20
      IF(IHELP.EQ.1) GO TO 20
      IF(NN.LT.NV) GO TO 23
      WRITE(IDLG,22)
22    FORMAT(' MORE VARIABLES SPECIFIED THAN POSSIBLE FOR THIS DATA')
      GO TO 20
23    DO 24 I=1,NN-1
      IF(IVAR(I).LT.0) GO TO 24
      DO 24 J=I,NN
      IF(IVAR(I).NE.IVAR(N)) GO TO 24
      WRITE(IDLG,25)NAMES(IVAR(I))
25    FORMAT(' VARIABLE "',A5,'" SPECIFIED TWICE')
      GO TO 20
24    CONTINUE
30    IF(ICC.NE.2) WRITE(IDLG,31)
31    FORMAT(' ENTER THE BREAKDOWN VARIABLE? ',$)
      CALL ALPHA(IBREAK,1,I,IRET,IHELP,IERR,NAMES,NV)
      IF(IERR.EQ.1) GO TO 30
      IF(IRET.EQ.1) RETURN
      IF(IHELP.EQ.1) GO TO 30
      IF(IBREAK.GT.0) GO TO 33
      WRITE(IDLG,32)
32    FORMAT(' "*" AND "ALL" MAY NOT BE USED FOR A BREAKDOWN')
      GO TO 30
33    DO 34 I=1,NN
      IF(IBREAK.NE.IVAR(I)) GO TO 34
      WRITE(IDLG,35)
35    FORMAT(' THE BREAKDOWN VARIABLE WAS ALSO SPECIFIED AS ONE OF'/
     1' THE VARIALBES TO BE ANALYSED')
      GO TO 30
34    CONTINUE
C
C     NOW WORK ON BREAKDOWN VARIABLE
C
      DO 40 J=1,NC
      JV(J)=J
40    AV(J)=DATA(J,IBREAK)
      CALL SORT(AV,JV,NC)
      IF(OPT(1).EQ.1) GO TO 100
C
C     RANGES ENTERED BY USER
C
70    IF(ICC.NE.2) WRITE(IDLG,71)
71    FORMAT(' ENTER RANGES FOR BREAKDOWNS'/)
      I=1
72    IF(ICC.NE.2) WRITE(IDLG,73)
73    FORMAT('+? ',$)
      READ(ICC,74,END=60) INPUT
74    FORMAT(50A1)
      IF(INPUT(1).EQ.'!') RETURN
      IF(INPUT(1).EQ.' ') GO TO 60
      IF((INPUT(1).EQ.'A').AND.(INPUT(2).EQ.'U').AND.(INPUT(3).EQ.'T')
     1.AND.(INPUT(4).EQ.'O')) GO TO 100
      IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'P')
     1.OR.(INPUT(4).NE.'P')) GO TO 65
      WRITE(IDLG,75)
75    FORMAT('0ENTER RANGES 1 PER LINE IN RESPONSE TO A QUESTION MARK.'/
     1' EACH RANGE IS COMPRISED OF A MINIMUM, A MAXIMUM, AND AN'/
     2' OPTIONAL RANGE NAME (UP TO 5 CHARACTERS).  EACH ENTRY IN'/
     3' THE RANGE IS SEPARATED FROM THE NEXT ENTRY BY A COMMA.'/
     4' FOR EXAMPLE:  1,7,LOW; THE MINIMUM IS 1, THE MAXIMUM FOR THE'/
     5' RANGE IS 7, AND THE NAME OF THAT GROUPING IS LOW.  WHEN THE'/
     6' LAST RANGE HAS BEEN ENTERED TYPE A CONTROL Z (^Z) OR A BLANK'/
     7' LINE.'//)
      GO TO 72
C
C     MINIMUM OF RANGE
C
65    DO 76 J=1,15
76    IVALUE(J)=' '
      K=1
      J=1
77    IF(INPUT(J).EQ.' ') GO TO 80
      IF(INPUT(J).EQ.',') GO TO 82
      IF(INPUT(J).EQ.'E') GO TO 79
      IF(INPUT(J).EQ.'-') GO TO 79
      IF(INPUT(J).EQ.'.') GO TO 79
      IF((INPUT(J).LE.'9').AND.(INPUT(J).GE.'0')) GO TO 79
      WRITE(IDLG,78) INPUT(J)
78    FORMAT(' ILLEGAL CHARACTER "',A1,'" IN MINIMUM OF RANGE'/)
      GO TO 72
79    IF(K.GT.15) GO TO 80
      IVALUE(K)=INPUT(J)
      K=K+1
      J=J+1
      GO TO 77
80    WRITE(IDLG,81)
81    FORMAT(' NO COMMA BETWEEN MINIMUM AND MAXIMUM OR EXTRA SPACES'/)
      GO TO 72
82    ENCODE(15,74,SPACE) IVALUE
      DECODE(15,83,SPACE) RANGE(1,I)
83    FORMAT(G)
      DO 84 K=1,15
84    IVALUE(K)=' '
      K=1
      J=J+1
C
C     MAXIMUM OF RANGE
C
85    IF(INPUT(J).EQ.' ') GO TO 90
      IF(INPUT(J).EQ.',') GO TO 90
      IF(INPUT(J).EQ.'E') GO TO 87
      IF(INPUT(J).EQ.'-') GO TO 87
      IF(INPUT(J).EQ.'.') GO TO 87
      IF((INPUT(J).LE.'9').AND.(INPUT(J).GE.'0')) GO TO 87
      WRITE(IDLG,86) INPUT(J)
86    FORMAT(' ILLEGAL CHARACTER "',A1,'" IN MAXIMUM OF RANGE'/)
      GO TO 72
87    IF(K.GT.15) GO TO 88
      IVALUE(K)=INPUT(J)
      K=K+1
      J=J+1
      GO TO 85
88    WRITE(IDLG,89)
89    FORMAT(' NO COMMA BETWEEN MAXIMUM AND NAME OR TOO MANY CHAR.'/)
      GO TO 72
90    ENCODE(15,74,SPACE) IVALUE
      DECODE(15,83,SPACE) RANGE(2,I)
      IF(INPUT(J).NE.',') GO TO 91
      IF(INPUT(J+1).EQ.' ') GO TO 91
      ENCODE(5,74,RANGE(3,I))(INPUT(K),K=J+1,J+5)
      GO TO 93
91    ENCODE(5,92,RANGE(3,I)) I
92    FORMAT('GRP',I2)
93    I=I+1
      IF(I.LE.25) GO TO 72
      WRITE(IDLG,94)
94    FORMAT(' MAXIMUM OF 25 RANGES')
      GO TO 60
60    NBRK=I-1
      GO TO 108
C
C     DISCR OR AUTO USED
C
100   I=1
      J=2
      RANGE(1,I)=AV(1)
      RANGE(2,I)=AV(1)
      ENCODE(5,92,RANGE(3,I)) I
101   IF(RANGE(1,I).EQ.AV(J)) GO TO 102
      I=I+1
      IF(I.GT.25) GO TO 103
      RANGE(1,I)=AV(J)
      RANGE(2,I)=AV(J)
      ENCODE(5,92,RANGE(3,I)) I
102   J=J+1
      IF(J.LE.NC) GO TO 101
      NBRK=I
      GO TO 105
103   XNC=(AV(NC)-AV(1))/25.
      RANGE(1,I)=AV(1)
      DO 104 I=1,24
      RANGE(2,I)=RANGE(1,I)+XNC
      RANGE(1,I+1)=RANGE(2,I)
104   CONTINUE
      RANGE(2,25)=AV(NC)
      NBRK=25
C
C     OK RANGES ESTABLISHED SEE IF USER WANTS LIST
C
105   IF(OPT(2).NE.1) GO TO 108
      WRITE(IDLG,95)
95    FORMAT(' RANGES USED:')
      DO 106 I=1,NBRK
106   WRITE(IDLG,107) RANGE(1,I),RANGE(2,I),RANGE(3,I)
107   FORMAT(1X,G,',',G,',',A5)
C
C     NOW ARRANGE GROUPS
C
108   DO 109 I=1,NBRK
109   NGRP(I)=0
      DO 110 I=1,NC
      DO 111 J=1,NBRK
      IF(AV(I).LT.RANGE(1,J)) GO TO 111
      IF(AV(I).GT.RANGE(2,J)) GO TO 111
      AV(I)=J
      NGRP(J)=NGRP(J)+1
      GO TO 110
111   CONTINUE
115   AV(I)=99
110   CONTINUE
      CALL SORT(AV,JV,NC)
      I=1
112   IF(NGRP(I).NE.0) GO TO 117
      WRITE(IDLG,113) RANGE(3,I)
113   FORMAT(' RANGE "',A5,'" HAD NO OBSERVATIONS - DISCARDED')
      IF(I.EQ.NBRK) GO TO 116
      DO 114 J=I+1,NBRK
      NGRP(J-1)=NGRP(J)
      DO 114 K=1,3
114   RANGE(K,J-1)=RANGE(K,J)
116   NBRK=NBRK-1
      IF(I.LE.NBRK) GO TO 112
      GO TO 118
117   I=I+1
      IF(I.LE.NBRK) GO TO 112
118   IF(NBRK.GT.1) GO TO 120
      WRITE(IDLG,119)
119   FORMAT(' ONLY 1 GROUP - NO DISCRIMINATE ANALYSIS PERFORMED')
      RETURN
C
C     OK NOW GET VARIABLES TO BE USED.
C
120   K=1
      DO 121 I=1,NN
      IV(I)=IVAR(I)
      IF(IVAR(I).GT.0) GO TO 121
      IV(I)=K
      K=K+1
121   CONTINUE
      GO TO 135
130   J=NN
131   IF(IVAR(J).GT.0) GO TO 132
      IV(J)=IV(J)+1
      IF(IV(J).LE.NV) GO TO 133
132   J=J-1
      IF(J.GE.1) GO TO 131
      RETURN
133   K=IV(J)
      IF(J.EQ.NN) GO TO 135
      DO 134 I=J+1,NN
      IF(IVAR(I).GT.0) GO TO 134
      K=K+1
      IF(K.GT.NV) GO TO 132
      IV(I)=K
134   CONTINUE
135   DO 136 I=1,NN-1
      DO 136 K=I+1,NN
      IF(IV(I).EQ.IV(K)) GO TO 130
136   CONTINUE
C
C     BEGIN ANALYSIS
C
      IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(J),J=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,137) NAMES(IBREAK),(NAMES(IV(I)),I=1,NN)
137   FORMAT('0',15X,'***** DISCRIMINATE ANALYSIS *****'/
     1'0VARIABLE ',A5,' IS USED TO DETERMINE GROUPINGS.  THE '/
     2' VARIABLES BEING ANALYSED ARE:',7(1X,A5)/11(1X,A5))
      LINES=9
      LINPO=(NN+7)/8+1
      IF(NN.EQ.8) LINPO=3
      LINPOB=(NBRK+7)/8+1
      IF(NBRK.EQ.8) LINPOB=3
      DO 140 J=1,NN
      XC(J)=0
      DO 141 K=1,NN
141   D(J,K)=0.
      DO 142 K=1,NBRK
142   XBAR(J,K)=0.
140   CONTINUE
      N=0
      NTOT=0
      DO 150 I=1,NBRK
      NTOT=NTOT+NGRP(I)
      DO 151 J=1,NN
      DO 151 K=1,NN
151   C(J,K)=0.
      XN=NGRP(I)
      DO 152 J=N+1,N+NGRP(I)
      JJ=JV(J)
      DO 153 K=1,NN
      X=DATA(JJ,IV(K))
      XBAR(K,I)=XBAR(K,I)+X
      XC(K)=XC(K)+X
      DO 154 L=K,NN
154   C(K,L)=C(K,L)+X*DATA(JJ,IV(L))
153   CONTINUE
152   CONTINUE
      DO 155 K=1,NN
155   XBAR(K,I)=XBAR(K,I)/NGRP(I)
      DO 156 K=1,NN
      DO 156 L=K,NN
      C(K,L)=C(K,L)-XN*XBAR(K,I)*XBAR(L,I)
156   C(L,K)=C(K,L)
      IF(OPT(7).EQ.0) GO TO 165
      IF(IOUT.NE.21) GO TO 400
      LINES=LINES+4+LINPO
      IF(LINES.LE.(LINPP-LINPO)) GO TO 400
      CALL PRNTHD
      LINES=6+LINPO
400   WRITE(IOUT,157) RANGE(3,I)
157   FORMAT(//'0MATRIX OF CROSS PRODUCTS OF DEVIATIONS FROM THE MEANS'
     1,' FOR GROUP: ',A5)
      IF(IOUT.NE.21) GO TO 401
      WRITE(IOUT,160) (NAMES(IV(K)),K=1,NN)
160   FORMAT('0',7X,8(2X,A5,8X)/(8X,8(2X,A5,8X)))
      GO TO 402
401   WRITE(IOUT,161) (NAMES(IV(K)),K=1,NN)
161   FORMAT('0',7X,5(2X,A5,6X)/(8X,5(2X,A5,6X)))
402   DO 158 K=1,NN
      IF(IOUT.NE.21) GO TO 404
      LINES=LINES+LINPO
      IF(LINES.LE.LINPP) GO TO 403
      CALL PRNTHD
      WRITE(IOUT,160)(NAMES(IV(KLP)),KLP=1,NN)
      LINES=2+2*LINPO
403   WRITE(IOUT,159) NAMES(IV(K)),(C(K,L),L=1,NN)
      GO TO 158
159   FORMAT('0',A5,2X,8G15.7/(8X,8G15.7))
404   WRITE(IOUT,162) NAMES(IV(K)),(C(K,L),L=1,NN)
162   FORMAT('0',A5,1X,5G13.7/(7X,5G13.7))
158   CONTINUE
165   DO 166 K=1,NN
      DO 166 L=1,NN
166   D(K,L)=D(K,L)+C(K,L)
      N=N+NGRP(I)
150   CONTINUE
      XN=N-NBRK
      IF(OPT(3).NE.1) GO TO 172
      IF(IOUT.NE.21) GO TO 405
      LINES=LINES+4+LINPOB
      IF(LINES.LE.(LINPP-LINPOB)) GO TO 405
      CALL PRNTHD
      LINES=6+LINPOB
405   WRITE(IOUT,171)
171   FORMAT(//'0','MEANS')
      IF(IOUT.NE.21) GO TO 407
406   WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
      GO TO 408
407   WRITE(IOUT,161)(RANGE(3,K),K=1,NBRK)
408   DO 170 I=1,NN
      IF(IOUT.NE.21) GO TO 410
      LINES=LINES+LINPOB
      IF(LINES.LE.LINPP) GO TO 409
      CALL PRNTHD
      WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
      LINES=2+2*LINPOB
409   WRITE(IOUT,159) NAMES(IV(I)),(XBAR(I,J),J=1,NBRK)
      GO TO 170
410   WRITE(IOUT,162) NAMES(IV(I)),(XBAR(I,J),J=1,NBRK)
170   CONTINUE
172   DO 175 K=1,NN
      DO 175 L=1,NN
175   C(K,L)=D(K,L)/XN
      IF(OPT(4).NE.1) GO TO 176
      IF(IOUT.NE.21) GO TO 411
      LINES=LINES+4+LINPO
      IF(LINES.LE.(LINPP-LINPO)) GO TO 411
      CALL PRNTHD
      LINES=6+LINPO
411   WRITE(IOUT,173)
173   FORMAT(//'0',5X,'DISPERSION MATRIX')
      IF(IOUT.NE.21) GO TO 413
412   WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
      GO TO 414
413   WRITE(IOUT,161)(NAMES(IV(K)),K=1,NN)
414   DO 174 I=1,NN
      IF(IOUT.NE.21) GO TO 416
      LINES=LINES+LINPO
      IF(LINES.LE.LINPP) GO TO 415
      CALL PRNTHD
      WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
      LINES=2+2*LINPO
415   WRITE(IOUT,159) NAMES(IV(I)),(C(I,J),J=1,NN)
      GO TO 174
416   WRITE(IOUT,162) NAMES(IV(I)),(C(I,J),J=1,NN)
174   CONTINUE
C
C
C     INVERSE OF DISPERSION MATRIX
C
176   DO 202 I=1,NN
      DO 201 J=1,NN
201   D(I,J)=0.
202   D(I,I)=1.0
      DO 210 I=1,NN
      IF(((C(I,I)+100.)-100.).NE.0.0) GO TO 220
      IF(I.EQ.N) GO TO 300
      DO 211 J=I+1,NN
      IF(((C(J,I)+100.)-100.).NE.0.0) GO TO 212
211   CONTINUE
      GO TO 300
212   DO 213 K=1,NN
      D(I,K)=D(I,K)+D(J,K)
213   C(I,K)=C(I,K)+C(J,K)
220   G=C(I,I)
      DO 221 J=1,NN
      D(I,J)=D(I,J)/G
221   C(I,J)=C(I,J)/G
      DO 230 L=1,NN
      IF(L.EQ.I) GO TO 230
      G=C(L,I)
      DO 231 J=1,NN
      D(L,J)=D(L,J)-G*D(I,J)
231   C(L,J)=C(L,J)-G*C(I,J)
230   CONTINUE
210   CONTINUE
C
C
      IF(OPT(5).NE.1) GO TO 242
      IF(IOUT.NE.21) GO TO 417
      LINES=LINES+4+LINPO
      IF(LINES.LE.(LINPP-LINPO)) GO TO 417
      CALL PRNTHD
      LINES=6+LINPO
417   WRITE(IOUT,240)
240   FORMAT(//'0','INVERSE OF DISPERSION MATRIX')
      IF(IOUT.NE.21) GO TO 419
418   WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
      GO TO 420
419   WRITE(IOUT,161)(NAMES(IV(K)),K=1,NN)
420   DO 241 I=1,NN
      IF(IOUT.NE.21) GO TO 422
      LINES=LINES+LINPO
      IF(LINES.LE.LINPP) GO TO 421
      CALL PRNTHD
      WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
      LINES=2+2*LINPO
421   WRITE(IOUT,159) NAMES(IV(I)),(D(I,J),J=1,NN)
      GO TO 241
422   WRITE(IOUT,162) NAMES(IV(I)),(D(I,J),J=1,NN)
241   CONTINUE
242   DO 243 I=1,NN
243   XC(I)=XC(I)/NTOT
      DSQ=0.
      DO 244 I=1,NN
      DO 244 J=1,NN
      DS=0
      DO 245 K=1,NBRK
245   DS=DS+NGRP(K)*(XBAR(I,K)-XC(I))*(XBAR(J,K)-XC(J))
244   DSQ=DSQ+DS*D(I,J)
      IDF=NN*(NBRK-1)
      PROB=FISHER(IDF,1000,DSQ)
      IF(IOUT.NE.21) GO TO 423
      LINES=LINES+8
      IF(LINES.LE.LINPP) GO TO 423
      CALL PRNTHD
      LINES=10
423   WRITE(IOUT,246) DSQ,DSQ,IDF,NBRK,NN,PROB
246   FORMAT(//'0GENERALIZED MAHALANOBIS D- SQUARE =',G/
     1'0USING THE VALUE ',G,' AS A CHI-SQUARE WITH ',I4/
     1' DEGREES OF FREEDOM TO TEST THE HYPOTHESIS THAT THE MEAN'/
     1' VALUES ARE THE SAME IN ALL ',I3,' GROUPS FOR THESE ',I3/
     4' VARIABLES, THE PROBABILITY IS ',F5.2)
      DO 250 I=1,NBRK
      DSS=0
      DO 251 J=1,NN
      DS=0
      DO 252 K=1,NN
      DS=DS+D(J,K)*XBAR(K,I)
252   DSS=DSS+D(J,K)*XBAR(J,I)*XBAR(K,I)
251   C(J,I)=DS
250   CON(I)=DSS*-.5
      IF(IOUT.NE.21) GO TO 420
      LINES=LINES+4+LINPOB
      IF(LINES.LE.(LINPP-LINPOB)) GO TO 424
      CALL PRNTHD
      LINES=6+LINPOB
424   WRITE(IOUT,260)
260   FORMAT(//'0FUNCTION COEFFICIENTS')
      IF(IOUT.NE.21) GO TO 425
      WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
      GO TO 426
425   WRITE(IOUT,161)(RANGE(3,K),K=1,NBRK)
426   DO 261 I=1,NN
      IF(IOUT.NE.21) GO TO 428
      LINES=LINES+LINPOB
      IF(LINES.LE.LINPP) GO TO 427
      CALL PRNTHD
      WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
      LINES=2+2*LINPOB
427   WRITE(IOUT,159) NAMES(IV(I)),(C(I,J),J=1,NBRK)
      GO TO 261
428   WRITE(IOUT,162) NAMES(IV(I)),(C(I,J),J=1,NBRK)
261   CONTINUE
      LINES=LINES+1+LINPOB
      IF(LINES.LE.LINPP) GO TO 429
      CALL PRNTHD
      WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
      LINES=2+2*LINPOB
429   WRITE(IOUT,262)
262   FORMAT(/)
      CONST='CONST'
      IF(IOUT.EQ.21) WRITE(IOUT,159) CONST,(CON(I),I=1,NBRK)
      IF(IOUT.NE.21) WRITE(IOUT,162) CONST,(CON(I),I=1,NBRK)
      IF((OPT(6).NE.1).AND.(OPT(8).NE.1)) GO TO 130
      IF(OPT(6).NE.1) GO TO 430
      LINES=LINES+9
      LINPOB=(NBRK+9)/10
      IF(NBRK.EQ.10) LINPOB=2
      IF(LINES.LE.(LINPP-2-LINPOB)) GO TO 431
      CALL PRNTHD
      LINES=11
431   WRITE(IOUT,272)
272   FORMAT(//'0','EVALUATION OF CLASSIFICATION FUNCTIONS',
     1' FOR EACH CASE')
      WRITE(IOUT,285)
285   FORMAT('0',9X,'FUNCTION'/12X,'WITH'/10X,'LARGEST',4X,'LARGEST'/
     32X,'OBS.',6X,'PROB.',6X,'PROB.')
430   DO 270 I=1,NBRK
      DO 270 J=1,NBRK
270   IP(I,J)=0
      N=0
      DO 271 I=1,NBRK
      IF(IOUT.NE.21) GO TO 432
      LINES=LINES+2
      IF(LINES.LE.(LINPP-LINPOB)) GO TO 432
      CALL PRNTHD
      WRITE(IOUT,285)
      LINES=9
432   WRITE(IOUT,273) RANGE(3,I)
273   FORMAT('0GROUP: ',A5)
      DO 274 J=1+N,NGRP(I)+N
      M=JV(J)
      XX=-1.E30
      DO 275 K=1,NBRK
      E=CON(K)
      DO 276 L=1,NN
276   E=E+C(L,K)*DATA(M,IV(L))
      IF(E.LE.XX) GO TO 275
      XX=E
      IX=K
275   XVAL(K)=E
      DS=0
      DO 277 JJ=1,NBRK
      XVAL(JJ)=EXP(XVAL(JJ)-XX)
277   DS=DS+XVAL(JJ)
      XX=1./DS
      DO 278 JJ=1,NBRK
278   XVAL(JJ)=XVAL(JJ)/DS
      IF(OPT(6).NE.1) GO TO 274
      IF(IOUT.NE.21) GO TO 434
      LINES=LINES+LINPOB
      IF(LINES.LE.LINPP) GO TO 433
      CALL PRNTHD
      WRITE(IOUT,285)
      LINES=7+LINPOB
433   IF(IOUT.EQ.21) WRITE(IOUT,279) JV(J),IX,XX,(XVAL(JJ),JJ=1,NBRK)
279   FORMAT(1X,I5,7X,I2,4X,F8.5,4X,10(2X,F8.5)/(31X,10(2X,F8.5)))
      GO TO 274
434   WRITE(IOUT,280) JV(J),IX,XX,(XVAL(JJ),JJ=1,NBRK)
280   FORMAT(1X,I5,7X,I2,4X,F8.5,4X,5(1X,F7.4)/(31X,5(F7.4)))
274   IP(I,IX)=IP(I,IX)+1
      N=N+NGRP(I)
271   CONTINUE
      IF(OPT(8).NE.1) GO TO 130
      LINPOB=(NBRK+19)/20
      IF(NBRK.EQ.20) LINPOB=2
      IF(IOUT.NE.21) GO TO 435
      LINES=LINES+7+LINPOB
      IF(LINES.LE.(LINPP-LINPOB-1)) GO TO 435
      CALL PRNTHD
      LINES=9+LINPOB
435   WRITE(IOUT,284)
284   FORMAT(//'0CLASSIFICATION MATRIX'/8X,'GROUP PLACED IN AS A',
     1' RESULT OF EVALUATING THE FUNCTION')
      IF(IOUT.EQ.21) WRITE(IOUT,281)(RANGE(3,I),I=1,NBRK)
281   FORMAT('0GROUP  ',20(1X,A5)/(8X,20(1X,A5)))
      IF(IOUT.NE.21) WRITE(IOUT,282) (RANGE(3,I),I=1,NBRK)
282   FORMAT('0GROUP  ',10(1X,A5)/(8X,10(1X,A5)))
      WRITE(IOUT,283)
283   FORMAT(1X)
      DO 290 I=1,NBRK
      IF(IOUT.NE.21) GO TO 437
      LINES=LINES+LINPOB
      IF(LINES.LE.LINPP) GO TO 436
      CALL PRNTHD
      WRITE(IOUT,281)(RANGE(3,J),J=1,NBRK)
      WRITE(IOUT,283)
      LINES=4+LINPOB*2
436   WRITE(IOUT,291) RANGE(3,I),(IP(I,J),J=1,NBRK)
291   FORMAT(1X,A5,2X,20(1X,I5)/(8X,20(1X,I5)))
      GO TO 290
437   WRITE(IOUT,292) RANGE(3,I),(IP(I,J),J=1,NBRK)
292   FORMAT(1X,A5,2X,10(1X,I5)/(8X,10(1X,I5)))
290   CONTINUE
      GO TO 130
300   WRITE(IDLG,301)
      IF(IOUT.EQ.21) WRITE(IOUT,301)
301   FORMAT(' INVERSE OF DISPERSION MATRIX DOES NOT EXIST')
      GO TO 130
      END
      SUBROUTINE SORT(AV,SP,NC)
      DIMENSION SP(1),AV(1),IL(20),IU(20)
      M=1
      II=1
      J=NC
71    IF(II.GE.J) GO TO 78
72    K=II
      IJ=(J+II)/2
      T=AV(IJ)
      S=SP(IJ)
      IF(AV(II).GT.T) GO TO 83
      IF(AV(II).LT.T) GO TO 73
      IF(SP(II).LE.S) GO TO 73
83    AV(IJ)=AV(II)
      AV(II)=T
      T=AV(IJ)
      SP(IJ)=SP(II)
      SP(II)=S
      S=SP(IJ)
73    LL=J
      IF(AV(J).LT.T) GO TO 84
      IF(AV(J).GT.T) GO TO 75
      IF(SP(J).GE.S) GO TO 75
84    AV(IJ)=AV(J)
      AV(J)=T
      T=AV(IJ)
      SP(IJ)=SP(J)
      SP(J)=S
      S=SP(IJ)
      IF(AV(II).GT.T) GO TO 85
      IF(AV(II).LT.T) GO TO 75
      IF(SP(II).LE.S) GO TO 75
85    AV(IJ)=AV(II)
      AV(II)=T
      T=AV(IJ)
      SP(IJ)=SP(II)
      SP(II)=S
      S=SP(IJ)
      GO TO 75
74    AV(LL)=AV(K)
      AV(K)=TT
      SP(LL)=SP(K)
      SP(K)=SS
75    LL=LL-1
      IF(AV(LL).GT.T) GO TO 75
      IF((AV(LL).EQ.T).AND.(SP(LL).GT.S)) GO TO 75
      TT=AV(LL)
      SS=SP(LL)
76    K=K+1
      IF(AV(K).LT.T) GO TO 76
      IF((AV(K).EQ.T).AND.(SP(K).LT.S)) GO TO 76
      IF(K.LE.LL) GO TO74
      IF((LL-II).LE.(J-K)) GO TO 77
      IL(M)=II
      IU(M)=LL
      II=K
      M=M+1
      GO TO 79
77    IL(M)=K
      IU(M)=J
      J=LL
      M=M+1
      GO TO 79
78    M=M-1
      IF(M.EQ.0) RETURN
      II=IL(M)
      J=IU(M)
79    IF((J-II).GE.11) GO TO 72
      IF(II.EQ.1) GO TO 71
      II=II-1
80    II=II+1
      IF(II.EQ.J) GO TO 78
      S=SP(II+1)
      T=AV(II+1)
      IF(AV(II).LT.T) GO TO 80
      IF((AV(II).EQ.T).AND.(SP(II).LE.S)) GO TO 80
      K=II
81    AV(K+1)=AV(K)
      SP(K+1)=SP(K)
      K=K-1
      IF(T.LT.AV(K)) GO TO 81
      IF((T.EQ.AV(K)).AND.(S.LT.SP(K))) GO TO 81
      AV(K+1)=T
      SP(K+1)=S
      GO TO 80
      END
