C++*********************************************************************
C
C  HCLS.F
C               USED ALLOCATE JAN 2001 ARDEAN LEITH
C
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
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 the *
C=* 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 GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C  HCLS(LUNF,LUNT,LUND)
C
C  CALL TREE:  HCLS -> HCLP -> DIST_P
C                              CHAVA
C                              DENDRO
C                              DENDRO2
C                              ARBRE
C                              DENLST 
C
C--*********************************************************************

         SUBROUTINE HCLS(LUNF,LUNT,LUND)

         INCLUDE 'CMBLOCK.INC' 
         INCLUDE 'CMLIMIT.INC' 

         REAL, ALLOCATABLE, DIMENSION(:) ::  Q
         CHARACTER(LEN=MAXNAM)           :: FILNAM
         CHARACTER(LEN=MAXNAM)           :: FILNAMT
         CHARACTER(LEN=1)                :: NULL

         NULL = CHAR(0)

         WRITE(NOUT,*) ' YOU MAY USE A _IMC, _PIX, or _SEQ FILE'
         WRITE(NOUT,*) ' '

         CALL FILERD(FILNAM,NLET,NULL,
     &              'CORAN/PCA FILE (e.g. CORAN_01_IMC)~',IRTFLG)
         IF (IRTFLG .NE. 0) RETURN

         IT = INDEX(FILNAM,'_SEQ')
         IF (IT .GT. 0) THEN
            ITYPE = 1
            WRITE(NOUT,*)' For _SEQ coordinates file ---'
            CALL OPAUXFILE(.FALSE.,FILNAM,DATEXC,-LUNF,0,
     &                 'O',' ',.TRUE.,IRTFLG)
            IF (IRTFLG .NE. 0) RETURN

            READ(LUNF)  NKLA,NFAC

         ELSEIF(INDEX(FILNAM,'_IMC') .GT. 0) THEN
            CALL OPAUXFILE(.FALSE.,FILNAM,DATEXC,LUNF,0,
     &                 'O',' ',.TRUE.,IRTFLG)
            IF (IRTFLG .NE. 0) RETURN
            ITYPE = 2
            WRITE(NOUT,*)' In image coordinates file ---'
            READ(LUNF,*)  NKLA,NFAC

         ELSEIF(INDEX(FILNAM,'_PIX') .GT. 0) THEN
            CALL OPAUXFILE(.FALSE.,FILNAM,DATEXC,LUNF,0,
     &                 'O',' ',.TRUE.,IRTFLG)
            IF (IRTFLG .NE. 0) RETURN
            ITYPE = 3
            WRITE(NOUT,*)' In pixel coordinates file ---'
            READ(LUNF,*)  NKLA,NFAC

         ELSE
            CALL ERRT(101,'INVALID INPUT FILE TYPE',NE)
            RETURN
         ENDIF

         WRITE(NOUT,*)' Number of objects: ', NKLA
         WRITE(NOUT,*)' Number of factors: ', NFAC

         MAXFAC = NFAC
         MAXFA  = NFAC
	 MINFAC = 1
         CALL  RDPRAI(INUMBR,NFAC,MAXFA, MINFAC,MAXFAC,
     &	          'FACTOR NUMBERS',NULL,IER)
     
C        DIMENSION  W(MAXFA),COO(NFAC),COB(NFAC),INUM(MAXFA)

         K_W   = 1 
	 
C        CALCULATE NUMBER OF FACTORS
         MDIM   = NKLA * (NKLA-1) / 2
         JFIN   = 2 * NKLA - 1

         MEMD   = IPALIGN64(K_W+NFAC)
         MEMPK  = IPALIGN64(MEMD+MDIM)
         MEMVAL = IPALIGN64(MEMPK+JFIN)
         MEMLA  = IPALIGN64(MEMVAL+JFIN)
         MEMLB  = IPALIGN64(MEMLA+NKLA)
         MEMNT  = IPALIGN64(MEMLB+NKLA)
         MEMNO  = IPALIGN64(MEMNT+JFIN)
         MEMTOT = IPALIGN64(MEMNO+JFIN)
         MEMIV  = MEMD
         MEMIW  = IPALIGN64(MEMIV+NKLA)
         MEMV   = IPALIGN64(MEMIW+NKLA)
         MEMW   = IPALIGN64(MEMV+NKLA)
         MEMTOT = MAX0(MEMTOT,IPALIGN64(MEMW+NKLA))
	 
         ALLOCATE(Q(MEMTOT),STAT=IRTFLG)
	 IF (IRTFLG .NE. 0) THEN
            CALL ERRT(46,'HCLS; Q',MEMTOT)
            RETURN
         ENDIF
 
         WRITE(NOUT,*) ' DYNAMIC MEMORY ALLOCATION: ',MEMTOT
	    
         DO I=0,MAXFA-1
            Q(K_W+I) = 1.0
         ENDDO

         DO I=0,MAXFA-1
            CALL  RDPRM(W1,NOT_USED,'FACTOR WEIGHT')
            IF (W1 .EQ. 0.0)  EXIT
            Q(K_W+I) = W1
         ENDDO

201      WRITE(NOUT,*) '  FACTOR WEIGHTS USED:'
         WRITE(NOUT,23)  (Q(K_W+I),I=0,MAXFA-1)
23       FORMAT(10(2X,G10.3))

         CALL RDPRI1S(MODE,NOT_USED,
     &                 'CLUSTERING CRITERION (1-5)',IRTFLG)
         IF (IRTFLG .NE. 0) GOTO 9999

         CALL  HCLP(NKLA,MDIM,JFIN,
     &         Q(K_W),INUMBR,MAXFA,NFAC,
     &         Q(MEMD),Q(MEMPK),Q(MEMVAL),
     &         Q(MEMLA),Q(MEMLB),Q(MEMNT),
     &         Q(MEMNO),Q(MEMIV),Q(MEMIW),Q(MEMV),
     &         Q(MEMW),LUNF,LUND,LUNT,MODE,ITYPE)
 
9999     IF (ALLOCATED(Q)) DEALLOCATE(Q) 
	    
         WRITE (NDAT,2600)
2600     FORMAT (/' ',10('-'),'END OF HIERARCHICAL CLUSTERING',10('-')/)

         END