C @(#)mdaif.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 MDAIF
C
C.AUTHOR
C
C  F. Murtagh, ST-ECF, Garching.             Version 1.0  16 June 1986
C  F. Murtagh                                Version 2.0  Oct. 1988
C                                            (New std. ifs.)
C
C.PURPOSE
C
C  Pass parameters to, and execute, MDA (Multiple Discriminant Analysis)
C
C.INPUT/OUTPUT
C
C  P1 - P3 contain parameters; these are: input table name, output
C  table name, and whether or not a second output table is required, to
C  store the eigenvectors.
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 columns of the input table are sent to the MDA routine;
C . description of the MDA routine is in the corresponding program;
C . the output produced consists of one (or two) tables, and descriptors
C   to the first of these (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 . assumed is: NG < M < NROW in the MDA routine.
C . note that NG eigenvalues are output, although the last is in theory 0;
C   this will indicate the limited (but usually adequate) accuracy achieved.
C
C.MODIFICATIONS
C
C  STANDARD INTERFACES, F. MURTAGH, AUG. 1988
C  New Table File System M.Peron SEP91
C-----------------------------------------------------------------------
        PROGRAM MDAIF
C
        CHARACTER*6     FORM
        CHARACTER*60     NAMEIN, NAMEOUT, TABEVEC
        CHARACTER*16    UNIT, LABEL, LABEL2
        CHARACTER*80    LINE
        INTEGER         MADRID, KUN, KNUL, DTYPE,TID,TID1,TID2
        INTEGER         NACTV,ISTAT,NCOLIN,NROW,I,NAC,NAR
        INTEGER         LEN,NSEL,IPTRIN,IWPTRIN
        INTEGER         NTOT,INULL,IPTRGP,IMISS,NG
        INTEGER         IACTV,NCOLOUT,NNCOL,NNROW
        INTEGER         NSC1,NAROUT,IWPTROUT,IPTROUT
        INTEGER         IEVEC,M,NNAC,NAROUT2,IWPTROUT2
        INTEGER         IADD0,IADD3,IADD4,IADD5,IADD6
        INTEGER         IADD7,IADD8,IADD9,IADD2
        INTEGER         IADD10,IADD11,IERR,IPRINT,IADD1,KUNIT 
        INTEGER         NSORTC,IPTROUT2
        INCLUDE         'MID_INCLUDE:TABLES.INC/NOLIST'
        COMMON          /VMR/MADRID(1)
        INCLUDE         'MID_INCLUDE:TABLED.INC/NOLIST'
C
C
C ...   Assuming this is our first action ...
C
        CALL STSPRO('MDAIF')
        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,16,TID1,ISTAT)
        CALL TBIGET(TID1,NCOLIN,NROW,NSORTC,NAC,NAR,ISTAT)
C
C ...   Some error checking on input.
C
        IF (NROW.LT.1.OR.NCOLIN.LT.1) THEN
           CALL STTPUT(' Nos. of rows/columns are less than 1.',ISTAT)
           CALL STTPUT(' What sort of a table is this ??',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        DO I = 1, NCOLIN
           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)
              CALL STTPUT(' Aborting.',ISTAT)
              GOTO 9000
           ENDIF
        ENDDO
        CALL CHSEL(TID1,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)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        CALL TDMGET(4*NROW*NCOLIN,IWPTRIN,ISTAT)
        NSEL = NROW
        CALL TBCMAP(TID1,1,IPTRIN,ISTAT)
        CALL MAPSM(MADRID(IPTRIN),MADRID(IWPTRIN),
     1              NAR,NROW,NCOLIN)
        NTOT = NROW*NCOLIN
        CALL CHNULL(MADRID(IWPTRIN),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)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        CALL TBCMAP(TID1,NCOLIN,IPTRGP,ISTAT)
        CALL TBL_CHGP(MADRID(IPTRGP),NROW,IMISS,NG)
C
        WRITE (LINE,90) NG
        CALL STTPUT(LINE,ISTAT)
  90    FORMAT(' Number of groups:',I6)
C
C       IF (NG.GE.NROW.OR.NG.GE.NCOLIN-1) THEN
        IF (NG.GT.NROW.OR.NG.GT.NCOLIN) THEN
           CALL STTPUT(' Unacceptable no. of groups present.',ISTAT)
           CALL STTPUT(' (NG .GT. NROW or NCOL)',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        IF (IMISS.GT.0) THEN
           CALL STTPUT(' Unacceptable group memberships:',ISTAT)
           CALL STTPUT(' Check last col. of input table',ISTAT)
           CALL STTPUT(' for missing sequence numbers.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Get the output table name.
C
        CALL STKRDC('P2',1,1,60,IACTV,NAMEOUT,KUN,KNUL,ISTAT)
C
C ...   Now create the output (projections) table.
C
C        NCOLOUT = NG-1
        NCOLOUT = NG
        CALL TBTINI(NAMEOUT,F_TRANS,17,NCOLOUT,NROW,TID,ISTAT)
        FORM = 'G14.6'
        UNIT = ' '
C       (The messing around with LABEL2 which follows is so that
C       leading zeros in a column label appear as such, and not as
C       blanks.)
        LABEL = 'DIS00'
        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 TBIGET(TID,NNCOL,NNROW,NSC1,NAC,NAROUT,ISTAT)
        CALL TDMGET(4*NCOLOUT*NROW,IWPTROUT,ISTAT)
        CALL TBCMAP(TID,1,IPTROUT,ISTAT)
        CALL MAPSM(MADRID(IPTROUT),MADRID(IWPTROUT),
     1              NAROUT,NROW,NCOLOUT)

C
        IEVEC = 0
        CALL STKRDC('INPUTC',1,1,60,IACTV,TABEVEC,KUN,KNUL,ISTAT)

        M = NCOLIN - 1
        IF (TABEVEC.NE.'NO') THEN
           IEVEC = 1
C
C ...      Create the output (eigenvectors) table.
C
           M = NCOLIN - 1
           CALL TBTINI(TABEVEC,F_TRANS,17,NCOLOUT,M,TID2,ISTAT)
C          (The messing around with LABEL2 which follows is so that
C          leading zeros in a column label appear as such, and not as
C          blanks.)
           LABEL = 'EVC00'
           LABEL2 = '10001'
           DO I = 1, NCOLOUT
              WRITE (LABEL2(1:5),100) I+10000
              LABEL(4:5) = LABEL2(4:5)
              CALL TBCINI(TID2,DTYPE,1,FORM,UNIT,LABEL,NSORTC,ISTAT)
           ENDDO
C
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)
           CALL TDMGET(4*NCOLOUT*M,IWPTROUT2,ISTAT)
           CALL TBCMAP(TID2,1,IPTROUT2,ISTAT)
           CALL MAPSM(MADRID(IPTROUT2),MADRID(IWPTROUT2),
     1                NAROUT2,M,NCOLOUT)
        ENDIF
C
C ...   Allocate storage.
C
        CALL GETSTOR(NROW*M+NG+M,IADD0)
C        CALL GETSTOR(NROW*M,IADD0)
C        CALL GETSTOR(NG,IADD1)
C        CALL GETSTOR(M,IADD2)
        IADD1 = IADD0+NROW*M
        IADD2 = IADD1+NG
         CALL GETSTOR(NG*M,IADD3)
        CALL GETSTOR(M*M,IADD4)
        CALL GETSTOR(NG*NG,IADD5)
        CALL GETSTOR(NG*NG,IADD6)
        CALL GETSTOR(M*NG,IADD7)
        CALL GETSTOR(M,IADD8)
        CALL GETSTOR(M,IADD9)
        CALL GETSTOR(NROW,IADD10)
        CALL GETSTOR(NG,IADD11)
C
C ...   We now have all the storage we need; APPL will look after
C       a few further odds and ends - e.g. setting hard-wired
C       options, rearranging output, etc.  Then it will call the
C       application program.
C
        IERR = 0
        IF (IEVEC.EQ.0) THEN
           CALL APPL1(NROW,M,NG,MADRID(IWPTRIN),MADRID(IPTRGP),
     X     MADRID(IADD0),IPRINT,MADRID(IWPTROUT),
     X     MADRID(IADD1),MADRID(IADD2),MADRID(IADD3),MADRID(IADD4),
     X     MADRID(IADD5),
     X     MADRID(IADD6),MADRID(IADD7),MADRID(IADD8),MADRID(IADD9),
     X     MADRID(IADD10),
     X     MADRID(IADD11),IERR)
           CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NROW,NCOLOUT)

        ELSE
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)
           CALL APPL2(NROW,M,NG,MADRID(IWPTRIN),MADRID(IPTRGP),
     X     MADRID(IADD0),IPRINT,MADRID(IWPTROUT),
     X     MADRID(IADD1),MADRID(IADD2),MADRID(IADD3),MADRID(IADD4),
     X     MADRID(IADD5),
     X     MADRID(IADD6),MADRID(IADD7),MADRID(IADD8),MADRID(IADD9),
     X     MADRID(IADD10),
     X     MADRID(IADD11),MADRID(IWPTROUT2),IERR)
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)
           CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NROW,NCOLOUT)
           CALL MAPBG(MADRID(IWPTROUT2),MADRID(IPTROUT2),
     1             NAROUT2,M,NCOLOUT)
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)


        ENDIF
        IF (IERR.EQ.1) THEN
           CALL STTPUT(' IERR = 1 on return from MDA.',ISTAT)
           CALL STTPUT
     X  (' No convergence in eigenroutine TQL2, called from MDA.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.EQ.2) THEN
           CALL STTPUT(' IERR = 2 on return from MDA.',ISTAT)
           CALL STTPUT(' Matrix for inversion was singular.',ISTAT)
           CALL STTPUT(' Check that Nrows > Ncols. > Ngps., or',ISTAT)
           CALL STTPUT
     X          (' that no cols. of same values are present.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from MDA.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Write out descriptors with output table.
C
        CALL DSCROUT(NAMEOUT,NCOLOUT+1,MADRID(IADD8),TID,KUNIT)
C
C ...   Finish up: write output file to disk, and release memory space.
C
        CALL RELSTOR(NROW*M+NG+M,IADD0)
        CALL RELSTOR(NG*M,IADD3)
        CALL RELSTOR(M*M,IADD4)
        CALL RELSTOR(NG*NG,IADD5)
        CALL RELSTOR(NG*NG,IADD6)
        CALL RELSTOR(M*NG,IADD7)
        CALL RELSTOR(M,IADD8)
        CALL RELSTOR(M,IADD9)
        CALL RELSTOR(NROW,IADD10)
        CALL RELSTOR(NG,IADD11)
C
        CALL TBTCLO(TID1,ISTAT)
C
        CALL TBIPUT(TID,NCOLOUT,NROW,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
        IF (IEVEC.NE.0) THEN
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)
           CALL TBIPUT(TID2,NCOLOUT,M,ISTAT)
           CALL TBSINI(TID2,ISTAT)
           CALL TBTCLO(TID2,ISTAT)
        ENDIF
C
 9000   CONTINUE
        CALL STSEPI
C
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL1(N,M,NG,AIN,GPR,ACREA,IPRINT,AOUT,GPS,MEAN,GPM,
     X                      TOT,BETW,BETW2,CPROJ,W1,W2,GP,NOG,IERR)
C
C ...   This is now the application routine - note the setting of
C       IPRINT, which in future may be passed if desired
C       on the command line.
C
        INTEGER M,N,NG,IPRINT,IERR,I,J
        REAL*4  AIN(N,M), AOUT(N,NG), ACREA(N,M),GPR(N),GPS(NG),MEAN(M),
     X          TOT(M,M), GPM(NG,M), BETW(NG,NG), BETW2(NG,NG),
     X          CPROJ(M,NG), W1(M), W2(M)
        INTEGER*4       GP(N), NOG(NG)
C
        IPRINT = 3
C
        DO I = 1, N
           GP(I) = GPR(I)
           DO J = 1, M
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
        DO I = 1, NG
           NOG(I) = GPS(I)
        ENDDO
C
        CALL MDA(N,M,NG,ACREA,GP,IPRINT,NOG,MEAN,GPM,TOT,BETW,
     X           BETW2,CPROJ,W1,W2,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, N
           DO J = 1, NG
              AOUT(I,J) = ACREA(I,J)
           ENDDO
        ENDDO
C
        RETURN
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL2(N,M,NG,AIN,GPR,ACREA,IPRINT,AOUT,GPS,MEAN,GPM,
     X            TOT,BETW,BETW2,CPROJ,W1,W2,GP,NOG,AOUT2,IERR)
C
C ...   This is now the application routine - note the setting of
C       IPRINT, which in future may be passed if desired
C       on the command line.
C
        INTEGER N,M,NG,IPRINT,I,J,IERR
        REAL*4  AIN(N,M), AOUT(N,NG), ACREA(N,M),GPR(N),GPS(NG),MEAN(M),
     X          TOT(M,M), GPM(NG,M), BETW(NG,NG), BETW2(NG,NG),
     X          CPROJ(M,NG), W1(M), W2(M), AOUT2(M,NG)
        INTEGER*4       GP(N), NOG(NG)
C
        IPRINT = 3
C
        DO I = 1, N
           GP(I) = GPR(I)
           DO J = 1, M
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
        DO I = 1, NG
           NOG(I) = GPS(I)
        ENDDO
C
        CALL MDA(N,M,NG,ACREA,GP,IPRINT,NOG,MEAN,GPM,TOT,BETW,
     X           BETW2,CPROJ,W1,W2,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, N
           DO J = 1, NG
              AOUT(I,J) = ACREA(I,J)
           ENDDO
        ENDDO
C
        DO I = 1, M
           DO J = 1, NG
              AOUT2(I,J) = CPROJ(I,J)
           ENDDO
        ENDDO
C
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHNULL(X,NT,NULL)
C ...   Check if null values are present.
        INTEGER I,NT,NULL
        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 TBL_CHGP(X,NT,IMISS,NG)
C ...   Determine no. of groups and check for consistent sequencing.
        INTEGER IMISS,NG,NT,I2,I
        REAL*4  X(NT),G
        IMISS = 0
        NG = 0
        G = NG
        DO I = 1, NT
           IF (X(I).GT.G) THEN
              G = X(I)
              NG = G
           ENDIF
        ENDDO
C
        DO I = 1, NG
           G = I
           DO I2 = 1, NT
              IF (X(I2).EQ.G) GOTO 200
           ENDDO
           IMISS = I
           IF (IMISS.GT.0) RETURN
  200      CONTINUE
        ENDDO
C
        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 (already allocated) for NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMFRE(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE DSCROUT(TAB,M,VALS,TID,KUNIT)
C ...   Output descriptors, with table.
        REAL*4 VALS(M),V
        INTEGER TID,M,I,KUNIT,ISTAT,INDX,MM
        CHARACTER TAB*8, TABNAME*12
C
C ...   The following values are rearranged in decreasing rather
C       than increasing order.
        DO I = 1, M/2
           V = VALS(I)
           VALS(I) = VALS(M-I+1)
           VALS(M-I+1) = V
        ENDDO
C
        INDX = INDEX(TAB//' ',' ')-1
        TABNAME = TAB(1:INDX)//'.TBL'
C
        MM = M-1
        CALL STDWRR(TID,'EIGENVALUES',VALS,1,MM,KUNIT,ISTAT)
C
        RETURN
        END
