C @(#)tcluster.for	19.1 (ES0-DMD) 02/25/03 13:27:33
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C
C.IDENTIFICATION
C
C  program TCLUSTER.FOR
C  F.Murtagh,J.D.Ponz           ESO/ST-ECF - Garching
C                               VERSION 1.0     22 JAN 85
C  F. Murtagh                   Version 2.0 New std. ifs. Oct. 1988
C
C  M. Peron                     Version 3.0 Sep 91
C.PURPOSE
C
C  Execute the command
C  CLUSTER/TABLE  intable outable [method]
C
C.KEYWORDS
C
C  Tables, Cluster Analysis
C
C.INPUT/OUTPUT
C
C  P1 - P3  contain input parameters
C
C.ALGORITHM
C
C . uses table interface routines
C . input table is assumed to contain entries in single precision,
C . No selection or null values
C . description of clustering algorithm in the corresponding routine
C . the output of the command consists of
C   .  9 columns with the class
C      number for each entry in  2 cluster classification, 3 cluster
C      and so on. Column labels :nCL where n=2,...,9
C   .  auxiliar descriptors
C      IORDER(9)
C      HEIGHT(9)
C      CRITVAL(9)
C      CRIT(9)
C
C.MODIFICATIONS
C
C F. Murtagh, Aug. 88 - new interface calls implemented.
C
C-----------------------------------------------------------
C
        IMPLICIT INTEGER*4      (A-Z)
        CHARACTER       INTAB*60, OUTAB*60, METHOD*4
        CHARACTER*80    LINE
        CHARACTER*16    LABEL,UNIT
        CHARACTER*6     FORM
        INTEGER*4       IORDER(9), HEIGHT(9), HVALS(9)
        REAL*4          CRITVAL(9)
C       REAL*4          AUX(9)
        LOGICAL         INULL
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INCLUDE         'MID_INCLUDE:TABLES.INC/NOLIST'
        COMMON          /VMR/MADRID(1)
        INCLUDE         'MID_INCLUDE:TABLED.INC/NOLIST'
        DATA            LABEL/'2CL'/,UNIT/' '/,FORM/'I4'/
C
C ... GET INTO MIDAS
C
        CALL STSPRO('TCLUSTER')
        DTYPE = D_R4_FORMAT
C
C ... GET INPUT PARAMETERS
C
        CALL STKRDC('P1',1,1,60,IAV,INTAB,KUN,KNUL,ISTAT)
        CALL STKRDC('P2',1,1,60,IAV,OUTAB,KUN,KNUL,ISTAT)
        CALL STKRDC('P3',1,1,4,IAV,METHOD,KUN,KNUL,ISTAT)
C        CALL UPCAS(METHOD,METHOD)
C
C ... READ INPUT TABLE
C
        CALL TBTOPN(INTAB, F_I_MODE, TID1, ISTAT)
        CALL TBIGET(TID1, NCOL,NROW,NSC1,NAC,NAR,ISTAT)

        CALL CHSEL(TID1,NROW,NSEL)
        IF (NSEL.LE.0) THEN
          CALL STTPUT(' Empty table, check SELECT', ISTAT)
          GOTO 1000
        ENDIF
        IF (NSEL.NE.NROW) THEN
          WRITE(LINE,100) NSEL,NROW
          CALL STTPUT(LINE, ISTAT)
          STOP
        ENDIF
C
        LENDIS = (NROW*(NROW-1))/2
C
C ... CREATE OUTPUT TABLE
C
        NCOLOUT = NCOL + 8
        NAC = NCOLOUT
        CALL TBTINI(OUTAB, F_TRANS, 17, NAC, NROW, TID, ISTAT)
        DO I = 1, NCOL
          CALL TBFGET(TID1, I, FORM, LEN, DTYPE, ISTAT)
          IF (DTYPE.NE.D_R4_FORMAT) THEN
            CALL STTPUT(' Illegal format', ISTAT)
            CALL STTPUT(' Only R*4 column type allowed' , ISTAT)
            GOTO 1000
          ENDIF
          CALL TBCMAP(TID1,I,IADD1,ISTAT)
          CALL TBLGET(TID1, I, LABEL, ISTAT)
          CALL TBUGET(TID1, I, UNIT, ISTAT)
          CALL TBCINI(TID, DTYPE, 1 , FORM, UNIT, LABEL, I1, ISTAT)
          CALL TBCMAP(TID,I1,IADD2,ISTAT)
          CALL CPYCOL(MADRID(IADD1), MADRID(IADD2), NROW)
        ENDDO
C
C ... CREATE CLUSTER COLUMNS
C
        LABEL = 'CL2'
        UNIT  = ' '
        FORM  = 'I4'
        DO I = 2, 9
          WRITE(LABEL(3:3),150) I
          CALL TBCINI(TID, DTYPE, 1, FORM, UNIT, LABEL, I1, ISTAT)
        ENDDO
C
C ... GET ADDRESSES
C
C     DATA ARRAY ADDRESS
        CALL TDMGET(4*NCOL*NSEL,IWADD1,ISTAT)
        CALL TBIGET(TID,NNCOL,NNROW,NSC1,NNAC,NAR,ISTAT)
        CALL TBCMAP(TID, 1, IADD1, ISTAT)
        CALL MAPSM(MADRID(IADD1),MADRID(IWADD1),
     1              NAR,NSEL,NCOL)

C
C ... CHECK IF NULL VALUE
C
        NDAT = NCOL*NROW
        CALL CHECK(MADRID(IWADD1),NDAT,INULL)
        IF (INULL) THEN
          CALL STTPUT(' Null entries in the table are not allowed'
     .                     , ISTAT)
          CALL STTPUT(' Use the SELECT command ' , ISTAT)
          GOTO 1000
        ENDIF
C     CLASS ARRAY ADDRESS
        CALL TDMGET(32*NSEL,IWADD2,ISTAT)
        CALL TBCMAP(TID, NCOL+1, IADD2, ISTAT)
C        CALL MIPSM(MADRID(IADD2),MADRID(IWADD2),
C     1              NAR,NSEL,8)
C     ADDRESS OF 5 VECTORS USED AS WORK SPACE
        CALL TDMGET(4*NROW,IWORK1,ISTAT)
        CALL TDMGET(4*NROW,IWORK2,ISTAT)
        CALL TDMGET(4*NROW,IWORK3,ISTAT)
        CALL TDMGET(4*NROW,IWORK4,ISTAT)
        CALL TDMGET(4*NROW,IWORK5,ISTAT)
        CALL TDMGET(4*NROW,IWORK6,ISTAT)
        CALL TDMGET(4*NROW,IWORK7,ISTAT)
        IF (METHOD.NE.'MNVR') CALL TDMGET(4*LENDIS,IWORK8,ISTAT)
C
C ... ALGORITHM
C
        IF (METHOD.EQ.'MNVR')
     .  CALL CLUS(NROW, NCOL, MADRID(IWADD1), MADRID(IWORK1),
     .          MADRID(IWORK3), MADRID(IWORK2), MADRID(IWORK4),
     .          MADRID(IWORK2), MADRID(IWORK1), MADRID(IWORK5))
        IF (METHOD.EQ.'MVAR')
     .  CALL CLST(NROW, NCOL, LENDIS, 1, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'SLNK')
     .  CALL CLST(NROW, NCOL, LENDIS, 2, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'CLNK')
     .  CALL CLST(NROW, NCOL, LENDIS, 3, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'ALNK')
     .  CALL CLST(NROW, NCOL, LENDIS, 4, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'WLNK')
     .  CALL CLST(NROW, NCOL, LENDIS, 5, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'MEDN')
     .  CALL CLST(NROW, NCOL, LENDIS, 6, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
        IF (METHOD.EQ.'CNTR')
     .  CALL CLST(NROW, NCOL, LENDIS, 7, MADRID(IWADD1),
     .          MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWORK4), MADRID(IWORK5), MADRID(IWORK6),
     .          MADRID(IWORK7), MADRID(IWORK8))
C
        CALL SEECL(NROW, MADRID(IWORK1), MADRID(IWORK3), MADRID(IWORK2),
     .          MADRID(IWADD2), HVALS, IORDER, CRITVAL, HEIGHT)
C
        CALL MIPBG(MADRID(IWADD2),MADRID(IADD2),
     .             NAR,NSEL,8)
        CALL DENDR(IORDER, HEIGHT, CRITVAL)
C
        CALL TDMFRE(4*NROW,IWORK1,ISTAT)
        CALL TDMFRE(4*NROW,IWORK2,ISTAT)
        CALL TDMFRE(4*NROW,IWORK3,ISTAT)
        CALL TDMFRE(4*NROW,IWORK4,ISTAT)
        CALL TDMFRE(4*NROW,IWORK5,ISTAT)
        CALL TDMFRE(4*NROW,IWORK6,ISTAT)
        CALL TDMFRE(4*NROW,IWORK7,ISTAT)
        CALL TDMFRE(32*NSEL,IWADD2,ISTAT)
        IF (METHOD.NE.'MNVR') CALL TDMFRE(4*LENDIS,IWORK8,ISTAT)
C
C ... CONVERT VALUES IN IADD2 FROM INTG*4 INTO REAL*4
C
        NTOT = 8*NROW
C        CALL CONVERT(MADRID(IADD2),MADRID(IADD2),NTOT)
C
C ... SAVE PARAMETERS AS DESCRIPTORS
C
C       CALL COPY(MADRID(IWORK2),AUX,NROW)
C       INDX      = INDEX(OUTAB//' ',' ') -1
C       TABLEFULL = OUTAB(1:INDX)//'.TBL'
C       CALL STDWRI(TID,'IORDER',IORDER,1,9,KUNIT,ISTAT)
C       CALL STDWRI(TID,'HEIGHT',HEIGHT,1,9,KUNIT,ISTAT)
C       CALL STDWRR(TID,'CRITVAL',CRITVAL,1,9,KUNIT,ISTAT)
C       CALL STDWRR(TID,'CRIT',AUX,1,9,KUNIT,ISTAT)
C
C ... COPY ORIGINAL ARRAY BACK INTO THE OUTPUT TABLE
C
C        CALL TBCMAP(TID1,0,IADD0,ISTAT)
C        DO I = 1, NCOL
C          CALL TBCMAP(TID1,I,IADD1,ISTAT)
C          CALL TBCMAP(TID,I,IADD2,ISTAT)
C          CALL CPYCOL(MADRID(IADD1), MADRID(IADD2), NROW)
C        ENDDO
C
C ... END
C
        CALL TBIPUT(TID,NCOLOUT,NROW,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
        CALL TBTCLO(TID1,ISTAT)
 1000   CONTINUE
        CALL STSEPI
 100    FORMAT(' No SELECT allowed; create new table.',I2)
 150    FORMAT(I1)
 200    FORMAT(' Selected ',I6,' rows from ',I6,' in total')
        END




C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        SUBROUTINE COPY(INPUT,AUX,NSEL)
C
C       COPY N-1 TO N-9 ELEMENTS FRON INPUT INTO AUX
C
        INTEGER  I,NSEL
        REAL*4  INPUT(NSEL),AUX(9)
C
        DO I = 1, 9
          AUX(I) = INPUT(NSEL-I)
        ENDDO
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        SUBROUTINE CONVERT(INPUT,OUTPUT,NTOT)
C
C       CONVERT FROM INTEGER*4 TO REAL*4
C
        INTEGER I,NTOT
        REAL*4          OUTPUT(NTOT)
        INTEGER*4       INPUT(NTOT)
C
        DO I = 1, NTOT
          OUTPUT(I) = FLOAT(INPUT(I))
        ENDDO
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        SUBROUTINE CHECK(X,NDAT,INULL)
C
C       CHECK IF NULL VALUE
C
        INTEGER         I,NDAT
        REAL*4          X(NDAT), TRNULL
        INTEGER         TINULL
        DOUBLE PRECISION TDNULL
        LOGICAL         INULL
        CALL TBMNUL(TINULL,TRNULL,TDNULL)
C
        INULL = .FALSE.
        DO I = 1, NDAT
          IF(X(I).EQ.TRNULL) GOTO 10
        ENDDO
        RETURN
10      INULL = .TRUE.
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        SUBROUTINE CPYCOL(A,B,N)   ! COPY VECT. TO VECT. / REAL ASSUMED.
        INTEGER I,N
        REAL*4  A(N), B(N)
        DO 100 I = 1, N
           B(I) = A(I)
  100   CONTINUE
        RETURN
        END








