C++*********************************************************************
C
C  HCLS.F
C               USED ALLOCATE JAN 2001 ARDEAN LEITH
C
C **********************************************************************
C=*                                                                    *
C=* This file is part of:   SPIDER - Modular Image Processing System.  *
C=* SPIDER System Authors:  Joachim Frank & ArDean Leith               *
C=* Copyright 1985-2010  Health Research Inc.,                         *
C=* Riverview Center, 150 Broadway, Suite 560, Menands, NY 12204.      *
C=* Email: spider@wadsworth.org                                        *
C=*                                                                    *
C=* SPIDER 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=* SPIDER 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=* You should have received a copy of the GNU General Public License  *
C=* along with this program. If not, see <http://www.gnu.org/licenses> *
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     :: Q(:)
         CHARACTER(LEN=MAXNAM) :: FILNAM
         CHARACTER(LEN=MAXNAM) :: FILNAMT
         CHARACTER(LEN=1)      :: NULL

         INTEGER *8            :: MDIM8,HUG8
         INTEGER               :: MDIM 

         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
         MDIM8  = NKLA 
         MDIM8  = MDIM8 * (NKLA-1) / 2   ! FOR BIG NKLA
         HUG8   = HUGE(IT) / 2

         IF (MDIM8 > HUG8) THEN
             WRITE(NOUT,*) '*** ALLOCATION NEEDED:',MDIM8
             WRITE(6,*) ' *** ALLOCATION NEEDED:',MDIM8
             CALL ERRT(101, 'EXCESSIVE ALLOCATION FOR DISTANCES',NE)
             RETURN
         ENDIF

         MDIM   = MDIM8

         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 = MAX(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