C @(#)partif.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
C.IDENTIFICATION
C
C  program PARTIF.FOR
C  F.Murtagh,J.D.Ponz           ESO/ST-ECF - Garching
C                               VERSION 1.0     22 JAN 85
C  F.Murtagh.                   VERSION 2.0     16 JUN 86
C  F.Murtagh                    Version 3.0     3  AUG 88 (new std. i/f's.)
C  M.Peron                      Version 4.0     9 SEP 91
C.PURPOSE
C
C  Execute the command
C  PARTITION/TABLE  intable outable nclasses method ngrps seed
C
C.KEYWORDS
C
C  Tables, Cluster Analysis, Partitioning, Minimal Distance Algorithm,
C  Exchange Algorithm.
C
C.INPUT/OUTPUT
C
C  P1 - P6  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 select or null values
C . description of clustering algorithm in the corresponding routine
C . the output of the command consists of
C   one column with the class number for each entry.
C
C-----------------------------------------------------------
        PROGRAM         PARTIF
C
        IMPLICIT INTEGER*4      (A-Z)
        INTEGER*4       NCLASS,METHOD,ITER,NNCOL,NROW
        REAL*4          CTOT
        CHARACTER*6     FORM
        CHARACTER*60     INTAB, OUTAB
        CHARACTER*16    LABEL,UNIT
        CHARACTER*80    LINE
        LOGICAL         INULL
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INCLUDE         'MID_INCLUDE:TABLES.INC/NOLIST'
        COMMON          /VMR/MADRID(1)
        INCLUDE         'MID_INCLUDE:TABLED.INC/NOLIST'
        DATA            FORM/'I4'/,LABEL/'CL'/,UNIT/' '/ 
C
C ... GET INTO MIDAS
C
        CALL STSPRO('PARTIF')
        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 STKRDI('INPUTI',1,1,IAV,NCLASS,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',2,1,IAV,METHOD,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',3,1,IAV,NGP0,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',4,1,IAV,ISEED,KUN,KNUL,ISTAT)
C
C ... READ INPUT TABLE
C
        CALL TBTOPN(INTAB, 16, TID1, ISTAT)
        CALL TBIGET(TID1,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
C        CALL TBCMAP(TID1,0,IADD0,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
        NSEL = NROW
C
C ... CREATE OUTPUT TABLE
C
        NCOLOUT = NCOL + 1
        CALL TBTINI(OUTAB, F_TRANS, 17, NCOLOUT, NSEL, 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 COLUMN
C
        LABEL = 'CL'
        UNIT  = ' '
        FORM  = 'I4'
        CALL TBCINI(TID, DTYPE, 1, FORM, UNIT, LABEL, I1, ISTAT)
C
C ... GET ADDRESSES
C
C     DATA ARRAY ADDRESS
        CALL TDMGET(4*NCOL*NSEL,IWADD1,ISTAT)
        CALL TBIGET(TID,NNCOL,NNROW,NSC1,NAC,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*NSEL
        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 TBCMAP(TID, NCOL+1, IADD2, ISTAT)
C
        LENSTOR = (NCLASS*NCOL + NCLASS + NSEL)
        CALL TDMGET(4*LENSTOR,IWORK1,ISTAT)
        CALL TDMGET(4*NCLASS*NCOL,IWORK2,ISTAT)
        CALL TDMGET(4*NCLASS,IWORK3,ISTAT)
C
C ... ALGORITHM
C
        IWORK2 = IWORK1 +4*NCLASS*NCOL
        IWORK3 = IWORK2 +4*NCLASS
        CALL RANDP(NSEL,NCLASS,MADRID(IADD2),ISEED)
        IF (METHOD.EQ.2) THEN
           CALL READIMA(MADRID(IWADD1),NSEL,NCOL)
           CALL EXCH(MADRID(IWADD1),NSEL,NCOL,MADRID(IADD2),
     X               NGP0,MADRID(IWORK3),MADRID(IWORK1),NCLASS,
     X               MADRID(IWORK2),
     X               CTOT,ITER,IERR)
        ELSE
           CALL MINDST(MADRID(IWADD1),NSEL,NCOL,MADRID(IADD2),
     X               NGP0,MADRID(IWORK3),MADRID(IWORK1),NCLASS,
     X               MADRID(IWORK2),
     X               CTOT,ITER,IERR)
        ENDIF
C
        IF (IERR.EQ.2) THEN
           CALL STTPUT
     X(' A group has < minimum allowed number of members.',ISTAT)
           CALL STTPUT
     X(' Reduce number of groups and try again.',ISTAT)
           STOP
        ENDIF
        IF (IERR.EQ.1) THEN
           CALL STTPUT(' Invalid group assignment has been detected',
     X                                                          ISTAT)
           CALL STTPUT(' (i.e. < 1 or > no. of groups).',ISTAT)
           CALL STTPUT(' Is the no. of gps. correctly specified?',
     X                                                          ISTAT)
           STOP
        ENDIF
        WRITE(LINE,140) CTOT, ITER
  140   FORMAT(' Sum of variances:',F16.4,'; No. of iterations:',I8)
        CALL STTPUT(LINE,ISTAT)

        CALL TDMFRE(4*LENSTOR,IWORK1,ISTAT)
        CALL TDMFRE(4*NCLASS*NCOL,IWORK2,ISTAT)
        CALL TDMFRE(4*NCLASS,IWORK3,ISTAT)
C
C ... CONVERT VALUES IN IADD2 FROM INTG*4 INTO REAL*4
C
        NTOT = NSEL
        CALL CONVERT(MADRID(IADD2),MADRID(IADD2),NTOT)
C
C ... COPY ORIGINAL ARRAY BACK INTO THE OUTPUT TABLE
C
        CALL TBCMAP(TID1,0,IADD0,ISTAT)
        DO I = 1, NCOL
          CALL TBCMAP(TID1,I,IADD1,ISTAT)
          CALL TBCMAP(TID,I,IADD2,ISTAT)
          CALL CPYCOL(MADRID(IADD1), MADRID(IADD2), NROW)
        ENDDO
C
C ... END
C
        CALL TBIPUT(TID,NCOLOUT,NSEL,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
        CALL TBTCLO(TID1,ISTAT)
1000    CONTINUE
        CALL STSEPI
100     FORMAT(' Must SELECT all rows',2I1)
200     FORMAT(' Selected ',I6,' rows from ',I6,' in total')
        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 COLUMN VECTOR
        INTEGER I,N
        REAL*4          A(N),B(N)
        DO I = 1,N
           B(I)=A(I)
        ENDDO
        RETURN
        END
        SUBROUTINE READIMA(A,NSEL,NCOL)
        INTEGER NSEL,NCOL,I,J
        REAL*4 A(NSEL,NCOL)
        DO 1000 J = 1,NCOL
          DO 2000 I = 1,NSEL
             WRITE(6,*)A(I,J)
2000      CONTINUE 
1000    CONTINUE
        RETURN
        END
