C @(#)partif.for 17.1.1.1 (ES0-DMD) 01/25/02 17:16:49 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