C @(#)carif.for	17.1.1.1 (ES0-DMD) 01/25/02 17:16:49
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 CARIF
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
C.PURPOSE
C
C  Pass parameters to, and execute, CA (Correspondence Analysis).
C  The rows, only, are studied (i.e. in the parameter space).
C
C.INPUT/OUTPUT
C
C  P1 - P4 contain parameters; these are: input table name, output
C  table name, number of factors wanted (optional, - default
C  = 3), and whether or not a second output table is required, to store
C  the factor coordinates.
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 CA routine;
C . description of the CA 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 STANDARD INTERFACES, F. MURTAGH, AUG. 1988
C  New Table File System,   M. Peron,   SEP91 
C-----------------------------------------------------------------------
        PROGRAM CARIF
C
        CHARACTER*60     NAMEIN, NAMEOUT, TABEVEC
        CHARACTER*6     FORM
        CHARACTER*16    UNIT, LABEL, LABEL2
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INTEGER         ISTAT,NACTV,NCOLIN,NROW,TID1
        INTEGER         NSORTC,NAC,NAR,I,IPTRIN
        INTEGER         LEN,NSEL,IWPTRIN,NTOT,INULL
        INTEGER         IACTV,TID,NCOLOUT,NNCOL,NNROW
        INTEGER         NSC1,NAROUT,IWPTROUT,IPTROUT
        INTEGER         IEVEC,TID2,NNAC,NAROUT2,IWPTROUT2
        INTEGER         IPTROUT2,IADD0,IADD1,IADD2,IADD3
        INTEGER         IADD4,IADD5,IADD6,IERR,KUNIT 
C
        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('CARIF')
        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)
        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
C
        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. of
C       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 output (projections) table.
C
        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 = 'FAC00'
        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)
        IF (TABEVEC.NE.'NO') THEN
           IEVEC = 1
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,NSORTC,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)
        CALL GETSTOR(NROW,IADD5)
        CALL GETSTOR(NCOLIN,IADD6)
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),
     X     MADRID(IADD5),MADRID(IADD6),IERR)
        CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NROW,NCOLOUT)

        ELSE
        CALL APPL2(NROW,NCOLIN,NCOLOUT,MADRID(IWPTRIN),MADRID(IWPTROUT),
     X     MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),MADRID(IADD3),
     X     MADRID(IADD4),
     X     MADRID(IWPTROUT2),MADRID(IADD5),MADRID(IADD6),IERR)
        CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NROW,NCOLOUT)
        CALL MAPBG(MADRID(IWPTROUT2),MADRID(IPTROUT2),
     1             NAROUT2,NCOLIN,NCOLOUT)

        ENDIF
        IF (IERR.EQ.1) THEN
           CALL STTPUT(' IERR = 1 on return from CA.',ISTAT)
           CALL STTPUT(' No convergence in eigenroutine.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.EQ.2) THEN
           CALL STTPUT(' IERR = 2 on return from CA.',ISTAT)
           CALL STTPUT(' All values in a row or column = 0.',ISTAT)
           CALL STTPUT
     X          (' Please remove this row or column and redo.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from CA.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Write out descriptors with output table.
C
        CALL DSCROUT(NAMEOUT,NCOLIN,MADRID(IADD3),TID,KUNIT)
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)
        CALL RELSTOR(NROW,IADD5)
        CALL RELSTOR(NCOLIN,IADD6)
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 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                          RR,CC,IERR)
C
C ...   This is now the application routine - note the setting of
C       IPRINT which in the future may be passed if desired
C       on the command line.
C
        INTEGER NR,NC1,NC2,IERR,I,IPRINT,J
        REAL*4  AIN(NR,NC1), AOUT(NR,NC2), ACREA(NR,NC1),A1(NC1,NC1),
     X          A2(NC1,NC1),W1(NC1),W2(NC1),RR(NR),CC(NC1)
C
        IPRINT = 1
C
        DO I = 1, NR
           DO J = 1, NC1
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        CALL CA(NR,NC1,ACREA,IPRINT,A1,W1,W2,A2,RR,CC,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, NR
           DO J = 1, NC2
              AOUT(I,J) = ACREA(I,J)
           ENDDO
        ENDDO
C
        RETURN
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL2(NR,NC1,NC2,AIN,AOUT,ACREA,A1,A2,W1,W2,
     X             EVEC,RR,CC,IERR)
C
C ...   This is now the application routine - note the setting of
C       IPRINT which in the future may be passed if desired
C       on the command line.
C
        INTEGER   NR,NC1,NC2,IERR,IPRINT,I,J
        REAL*4  AIN(NR,NC1), AOUT(NR,NC2), ACREA(NR,NC1),A1(NC1,NC1),
     X          A2(NC1,NC1),W1(NC1),W2(NC1),EVEC(NC1,NC2),RR(NR),CC(NC1)
C
        IPRINT = 1
C
        DO I = 1, NR
           DO J = 1, NC1
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        CALL CA(NR,NC1,ACREA,IPRINT,A1,W1,W2,A2,RR,CC,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, NR
           DO J = 1, NC2
              AOUT(I,J) = ACREA(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.
        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 GETSTOR(NVALS,IPTR)
        INTEGER NVALS,IPTR,ISTAT
C ...   Allocate storage space for NVALS real values.
        CALL TDMGET(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE RELSTOR(NVALS,IPTR)
        INTEGER NVALS,IPTR,ISTAT
C ...   Release storage space (already allocated) for NVALS real values.
        CALL TDMFRE(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C---------------------------------------------------------------------
        SUBROUTINE DSCROUT(TAB,M,VALS,TID1,KUNIT)
C ...   Output descriptors, with table.
        INTEGER M,TID1,KUNIT,I,INDX,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
C       CALL PUTDSCR_ST(TABNAME,'EIGENVALUES','R',VALS,1,M,ISTAT)
        CALL STDWRR(TID1,'EIGENVALUES',VALS,1,M,KUNIT,ISTAT)
C
        RETURN
        END