C @(#)cmdsif.for	19.1 (ES0-DMD) 02/25/03 13:27:32
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 CMDSIF
C
C.AUTHOR
C
C  F. Murtagh, ST-ECF, Garching.             Version 1.0  17 June 1986
C  F. Murtagh                                Version 2.0  Oct. 1988
C                                            (New std. ifs.)
C
C.PURPOSE
C
C  Pass parameters to, and execute, CMDS
C  (Classical - Torgerson's, Gower's - Multidimensional Scaling,
C  Factor Analysis of a Distances Matrix, Principal Coordinates Analysis).
C
C.INPUT/OUTPUT
C
C  P1 - P3 contain parameters; these are: input table name, output
C  table name, number of principal coordinates wanted (optional, - default
C  = 3).
C
C.ALGORITHM
C
C . uses Midas Table interface routines;
C . input table is assumed to contain entries in single precision;
C . no select or null values 
C . all input table is sent to the CMDS routine;
C . description of the CMDS routine is in the corresponding program;
C . the output produced consists of a table, and descriptors associated
C   with this (eigenvalues): view with
C   READ/DESC outtable.TBL *
C . storage is limited only by the Midas Table system, and by the
C   overall system limitations;
C
C.MODIFICATIONS
C   NEW STANDARD INTERFACES, F. MURTAGH, AUG. 88
C   New Table File System,   M. Peron , SEP91
C
C-----------------------------------------------------------------------
        PROGRAM CMDSIF
C
        CHARACTER*60     NAMEIN, NAMEOUT
        CHARACTER*6     FORM 
        CHARACTER*16    UNIT, LABEL, LABEL2
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INTEGER         NACTV,ISTAT,TID1,NROW,NCOL,NSORTC
        INTEGER         NAC,NAR,IPTR,I,LEN,NSEL,NTOT,INULL
        INTEGER         IACTV,NCOLOUT,TID,IPTROUT
        INTEGER         IADD0,IADD1,IADD2,IADD3,KUNIT
        INCLUDE         'MID_INCLUDE:TABLES.INC/NOLIST'
        COMMON          /VMR/MADRID(1)
        INCLUDE         'MID_INCLUDE:TABLED.INC/NOLIST'
        DATA            FORM/'G14.6'/, UNIT/' '/
C
C ...   Assuming this is our first action ...
C
        CALL STSPRO('CMDSIF')
        DTYPE = D_R4_FORMAT
C
C ...   Get table name as a keyword passed in the command line.
C
        CALL STKRDC('P1',1,1,60,NACTV,NAMEIN,KUN,KNUL,ISTAT)
C
C ...   Read input table.
C
        CALL TBTOPN(NAMEIN,0,TID1,ISTAT)
        CALL TBIGET(TID1,NCOL,NROW,NSORTC,NAC,NAR,ISTAT)
        CALL TBCMAP(TID1,0,IPTR,ISTAT)
C
C ...   Some error checking on input.
C
        IF (NROW.LT.1.OR.NCOL.LT.1) THEN
           CALL STTPUT(' Nos. of rows/columns are less than 1.',ISTAT)
           CALL STTPUT(' What sort of a table is this ??',ISTAT)
           STOP
        ENDIF
C
        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)
              STOP
           ENDIF
        ENDDO
C
        CALL CHSEL(MADRID(IPTR),NROW,NSEL)
        IF (NSEL.NE.NROW) THEN
           CALL STTPUT(' Not all rows are SELECTed. ',ISTAT)
           CALL STTPUT(' In current implementation, MUST select all.',
     X                                                     ISTAT)
           STOP
        ENDIF
C
        CALL TBCMAP(TID1,1,IPTR,ISTAT)
        NTOT = NROW*NCOL
        CALL CHNULL(MADRID(IPTR),NTOT,INULL)
        IF (INULL.NE.0) THEN
           CALL STTPUT
     X     (' Null entries in the table are not allowed.',ISTAT)
           CALL STTPUT
     X     (' Use SELECT, and then construct another table.',
     X                                                    ISTAT)
           STOP
        ENDIF
C
C ...   OUTPUT TABLE - PREPARE
C
C ...   First get the output table name and no. cols. via keywords.
C
        CALL STKRDC('P2',1,1,60,IACTV,NAMEOUT,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',1,1,IACTV,NCOLOUT,KUN,KNUL,ISTAT)
C
C ...   Now create the table.
C
        CALL TBTINI(NAMEOUT,F_TRANS,17,NCOLOUT,NROW,TID,ISTAT)
        LABEL = 'NEW00'
        LABEL2 = '10001'
        DO I = 1, NCOLOUT
           WRITE (LABEL2(1:5),100) I+10000
           LABEL(4:5) = LABEL2(4:5)
           CALL TBCINI(TID,DTYPE,1,FORM,UNIT,LABEL,NSORTC,ISTAT)
        ENDDO
  100   FORMAT(I5)
C
        CALL TBCMAP(TID,1,IPTROUT,ISTAT)
C
C ...   ALLOCATE STORAGE
C
        CALL GETSTOR(NROW*NROW,IADD0)
        CALL GETSTOR(NROW*NROW,IADD1)
        CALL GETSTOR(NROW,IADD2)
        CALL GETSTOR(NROW,IADD3)
C
C ...   DO THE WORK
C
        CALL APPL(NROW,NCOL,NCOLOUT,MADRID(IPTR),MADRID(IPTROUT),
     X  MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),MADRID(IADD3))
C
C ...   WRITE DESCRIPTORS
C
        CALL DSCROUT(NAMEOUT,NROW,MADRID(IADD2),TID,KUNIT)
C
C ...   FINISH UP
C
C
        CALL RELSTOR(NROW*NROW,IADD0)
        CALL RELSTOR(NROW*NROW,IADD1)
        CALL RELSTOR(NROW,IADD2)
        CALL RELSTOR(NROW,IADD3)
C
        CALL TBTCLO(TID1,ISTAT)
C
        CALL TBIPUT(TID,NCOLOUT,NROW,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
C
        CALL STSEPI
C
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL(NR,NC1,NC2,AIN,AOUT,ACREA,A,W1,W2)
        INTEGER NR,NC1,NC2,IERR,IPRINT,I,J,ISTAT
        REAL*4  AIN(NR,NC1), AOUT(NR,NC2), ACREA(NR,NR),A(NR,NR),
     X                  W1(NR),W2(NR)
        IERR = 0
        IPRINT = 1
        DO I = 1, NR
           DO J = 1, NC1
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        IF (NR.NE.NC1) THEN
           CALL STTPUT
     X          (' A symmetric distances matrix is expected.',ISTAT)
           CALL STTPUT(' The input matrix is not symmetric !',ISTAT)
           STOP
        ENDIF
C
        CALL CMDS(NR,ACREA,IPRINT,W1,W2,A,IERR)
C
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from CMDS.',ISTAT)
           STOP
        ENDIF
C
        DO I = 1, NR
           DO J = 1, NC2
              AOUT(I,J) = ACREA(I,J)
           ENDDO
        ENDDO
C
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHSEL(MASK,NROW,NSEL)
C ...   Count table rows which are SELECTed.
        INTEGER I,NSEL,NROW
        REAL*4  MASK(NROW),TBLSEL
        DOUBLE PRECISION TDTRUE,TDFALS
        CALL TBMCON(TBLSEL,TDTRUE,TDFALS)
        NSEL = 0
        DO I = 1, NROW
           IF (MASK(I).EQ.TBLSEL) NSEL = NSEL + 1
        ENDDO
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHNULL(X,NT,NULL)
C ...   Check if null values are present.
        INTEGER I,NULL,NT
        REAL*4  X(NT), TRNULL
        INTEGER TINULL
        DOUBLE PRECISION TDNULL
        CALL TBMNUL(TINULL,TRNULL,TDNULL)
        NULL = 0
        DO I = 1, NT
           IF (X(I).EQ.TRNULL) NULL = NULL + 1
           IF (NULL.GT.0) RETURN
        ENDDO
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE GETSTOR(NVALS,IPTR)
C ...   Allocate storage space for NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMGET(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE RELSTOR(NVALS,IPTR)
C ...   Release storage space of NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMFRE(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE DSCROUT(TAB,N,VALS,TID,KUNIT)
C ...   Output descriptors, with table.
        INTEGER I,N,TID,KUNIT,INDX,ISTAT
        REAL*4 VALS(N),V
        CHARACTER*8  TAB 
        CHARACTER*12 TABNAME
C
        DO I = 1,(N-1)/2
           V = VALS(I)
           VALS(I) = VALS(N-I)
           VALS(N-I) = V
        ENDDO
        INDX = INDEX(TAB//' ',' ')-1
        TABNAME = TAB(1:INDX)//'.TBL'
        CALL STDWRR(TID,'EIGENVALUES',VALS,1,N-1,KUNIT,ISTAT)
        RETURN
        END
