C @(#)pcacif.for	17.1.1.1 (ES0-DMD) 01/25/02 17:16:50
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 PCACIF
C
C.AUTHOR
C
C  F. Murtagh, ST-ECF, Garching.             Version 1.0  9 May 1986
C  F. Murtagh                                Version 2.0  Oct. 1988
C                                            (New std. ifs.)
C  M. Peron                                  Version 3.0  Sep 91
C.PURPOSE
C
C  Pass parameters to, and execute, PCA (Principal Components Analysis).
C  The columns, only, are studied (i.e. in the objects space).
C
C.INPUT/OUTPUT
C
C  P1 - P5 contain parameters; these are: input table name, output
C  table name, PCA option (3 = on correlations, 2 = on covariances, 1 =
C  on SSCPs), number of principal components wanted (optional, - default
C  = 3), and whether or not a second output table is required, to store
C  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 nulls;
C . all columns of the input table are sent to the PCA routine;
C . description of the PCA 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: NROW > 7 or NCOLIN (this isn't a restrictive assumption).
C
C.MODIFICATIONS
C
C   New MIDAS interface calls, F. Murtagh August 1988.
C
C-----------------------------------------------------------------------
        PROGRAM PCACIF
C
        CHARACTER*60     NAMEIN, NAMEOUT, TABEVEC
        CHARACTER*6     FORM
        CHARACTER*16    UNIT, LABEL, LABEL2
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INTEGER         NACTV,ISTAT,TID1,NCOLIN,NROW,NSORTC
        INTEGER         NAC,NAR,I,LEN,NSEL,IWPTRIN,IPTRIN
        INTEGER         NTOT,INULL,IACTV,IERR,TID,TID2
        INTEGER         IADD0,IADD1,IADD2,IADD3,IADD4
        INTEGER         IPTROUT,IWPTROUT,IWPTROUT2
        INTEGER         NNCOL,NNROW,NSC1,NAROUT,IEVEC
        INTEGER         NCOLOUT,IOPTION,IPTROUT2,NNAC
        INTEGER         NAROUT2 
        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('PCACIF')
        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        CALL TBCMAP(TID1,0,IPTRIN,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)
        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
C
C ...   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,IOPTION,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',2,1,IACTV,NCOLOUT,KUN,KNUL,ISTAT)
C
C ...   Now create the output (projections) table.
C
        CALL TBTINI(NAMEOUT,F_TRANS,17,NCOLOUT,NCOLIN,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 = '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,I,ISTAT)
        ENDDO
  100   FORMAT(I5)
C
        CALL TBIGET(TID,NNCOL,NNROW,NSC1,NAC,NAROUT,ISTAT)
        CALL TDMGET(4*NCOLOUT*NCOLIN,IWPTROUT,ISTAT)
        CALL TBCMAP(TID,1,IPTROUT,ISTAT)
        CALL MAPSM(MADRID(IPTROUT),MADRID(IWPTROUT),
     1              NAROUT,NCOLIN,NCOLOUT)
C
        IEVEC = 0
        CALL STKRDC('INPUTC',1,1,60,IACTV,TABEVEC,KUN,KNUL,ISTAT)
        IF (TABEVEC.NE.'NO') THEN
           IEVEC = 1
           IF (NROW.LT.NCOLIN) THEN
           CALL STTPUT(' Warning.',ISTAT)
           CALL STTPUT
     X     (' Number of rows is usually > number of columns;',ISTAT)
           CALL STTPUT
     X     (' Eigenvector coords. in NCOLIN-dimension space are filed.',
     X                                                          ISTAT)
           CALL STTPUT('  ',ISTAT)
           ENDIF
C
C ...      Create the output (eigenvectors) table.
C
           CALL TBTINI(TABEVEC,F_TRANS,17,NCOLOUT,NCOLIN,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,I,ISTAT)
           ENDDO
C
           CALL TBIGET(TID2,NNCOL,NNROW,NSC1,NNAC,NAROUT2,ISTAT)
           CALL TDMGET(4*NCOLOUT*NCOLIN,IWPTROUT2,ISTAT)
           CALL TBCMAP(TID2,1,IPTROUT2,ISTAT)
           CALL MAPSM(MADRID(IPTROUT2),MADRID(IWPTROUT2),
     1                NAROUT2,NCOLIN,NCOLOUT) 
        ENDIF
C

C
C ...   Allocate storage.
C
        CALL GETSTOR(NROW*NCOLIN,IADD0)
        CALL GETSTOR(NCOLIN*NCOLIN,IADD1)
        CALL GETSTOR(NCOLIN*NCOLIN,IADD2)
        CALL GETSTOR(NCOLIN,IADD3)
        CALL GETSTOR(NCOLIN,IADD4)
C
C ...   We now have all the storage we need; APPL will look after
C       a few further odds and ends - e.g. setting a few 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,NCOLIN,NCOLOUT,MADRID(IWPTRIN),MADRID(IWPTROUT),
     X     MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),MADRID(IADD3),
     X     MADRID(IADD4),IOPTION,IERR)
        CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NCOLIN,NCOLOUT)
        ELSE
        CALL APPL2(NROW,NCOLIN,NCOLOUT,MADRID(IWPTRIN),MADRID(IWPTROUT),
     X     MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),MADRID(IADD3),
     X     MADRID(IADD4),MADRID(IWPTROUT2),IOPTION,IERR)
        CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NCOLIN,NCOLOUT)
        CALL MAPBG(MADRID(IWPTROUT2),MADRID(IPTROUT2),
     1             NAROUT2,NCOLIN,NCOLOUT)
        ENDIF
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from PCA.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Write out descriptors with output table.
C
        CALL DSCROUT(NAMEOUT,NROW,NCOLIN,MADRID(IADD3),TID)
C
C ...   Finish up: write output file to disk, and release memory space.
C
        CALL RELSTOR(NROW*NCOLIN,IADD0)
        CALL RELSTOR(NCOLIN*NCOLIN,IADD1)
        CALL RELSTOR(NCOLIN*NCOLIN,IADD2)
        CALL RELSTOR(NCOLIN,IADD3)
        CALL RELSTOR(NCOLIN,IADD4)
C
        CALL TBTCLO(TID1,ISTAT)
C
        CALL TBIPUT(TID,NCOLOUT,NCOLIN,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
        IF (IEVEC.NE.0) THEN
           CALL TBIPUT(TID2,NCOLOUT,NCOLIN,ISTAT)
           CALL TBSINI(TID2,ISTAT)
           CALL TBTCLO(TID2,ISTAT)
        ENDIF
C
 9000   CONTINUE
        CALL STSEPI
C
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL1(NR,NC1,NC2,AIN,AOUT,ACREA,A1,A2,W1,W2,
     X                  METHOD,IERR)
C
C ...   This is now the application routine - note the setting of
C       parameters, which in a future epoch may be passed if desired
C       on the command line.
C
        INTEGER NR,NC1,NC2,METHOD,IERR,IPRINT,I,J
        REAL*4  AIN(NR,NC1), AOUT(NC1,NC2), ACREA(NR,NC1),A1(NC1,NC1),
     X          A2(NC1,NC1),W1(NC1),W2(NC1)
C
        IPRINT = 3
C
        DO I = 1, NR
           DO J = 1, NC1
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        CALL PCA(NR,NC1,ACREA,METHOD,IPRINT,A1,W1,W2,A2,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, NC1
           DO J = 1, NC2
              AOUT(I,J) = A1(I,J)
           ENDDO
        ENDDO
C
        RETURN
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL2(NR,NC1,NC2,AIN,AOUT,ACREA,A1,A2,W1,W2,
     X             EVEC,METHOD,IERR)
C
C ...   This is now the application routine - note the setting of
C       parameters, which in a future epoch may be passed if desired
C       on the command line.
C
        INTEGER NR,NC1,NC2,I,J,IERR,IPRINT,METHOD
        REAL*4  AIN(NR,NC1), AOUT(NC1,NC2), ACREA(NR,NC1),A1(NC1,NC1),
     X          A2(NC1,NC1),W1(NC1),W2(NC1),EVEC(NC1,NC2)
C
        IPRINT = 1
C
        DO I = 1, NR
           DO J = 1, NC1
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        CALL PCA(NR,NC1,ACREA,METHOD,IPRINT,A1,W1,W2,A2,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, NC1
           DO J = 1, NC2
              AOUT(I,J) = A1(I,J)
           ENDDO
        ENDDO
C
        DO I = 1, NC1
           DO J = 1, NC2
              EVEC(I,J) = A2(I,NC1-J+1)
           ENDDO
        ENDDO
C
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHNULL(X,NT,NULL)
C ...   Check if null values are present.
        REAL*4  X(NT), TRNULL
        INTEGER TINULL,NULL,NT,I
        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 (already allocated) for NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMFRE(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE DSCROUT(TAB,N,M,VALS,TID2)
C ...   Output descriptors, with table.
        INTEGER N,M,TID2,I,INDX,MM
        INTEGER KUNIT,ISTAT
        REAL*4 VALS(M),V
        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 = MIN0(N,M)
        CALL STDWRR(TID2,'EIGENVALUES',VALS,1,MM,KUNIT,ISTAT)
C
        RETURN
        END