C @(#)pcarif.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 PCARIF 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, PCA (Principal Components Analysis). C The rows, only, are studied (i.e. in the parameter space). C C.INPUT/OUTPUT C C P1 - P5 contain parameters; these are: input table name, output C table name, analysis option (3 = correlations, 2 = covariances, 1 = C 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 selects 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 Modified for new interface calls, F. Murtagh, August 1988. C Modified for new Table File System M.Peron Sep 91 C----------------------------------------------------------------------- PROGRAM PCARIF C CHARACTER*60 NAMEIN, NAMEOUT, TABEVEC CHARACTER*6 FORM CHARACTER*16 UNIT, LABEL, LABEL2 INTEGER MADRID, KUN(1), KNUL, DTYPE INTEGER NACTV,ISTAT,TID1,NCOLIN,NROW,NSORTC,NAC INTEGER NAR,IPTRIN,I,LEN,NSEL,IWPTRIN,NTOT INTEGER INULL,IACTV,TID,IOPTION,IWPTROUT2 INTEGER NNCOL,NNROW,NSC1,NAROUT,NCOLOUT,IPTROUT INTEGER IWPTROUT,TID2,IEVEC,NNAC,NAROUT2 INTEGER IPTROUT2,IADD0,IADD1,IADD2,IADD3 INTEGER IADD4,IERR C C ... Must do the following ... C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:TABLED.INC' C C C ... Assuming this is our first action ... C CALL STSPRO('PCARIF') 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,F_IO_MODE,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 CALL CHSEL(TID1,NROW,NSEL) IF (NSEL.NE.NROW) THEN CALL STTPUT(' Not all rows are SELECTed. ',ISTAT) CALL STTPUT X (' 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 PCA analysis option 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,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,NCOLIN,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 = '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*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(1:3).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,I,ISTAT) ENDDO 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,NROW,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,NROW,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,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 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 REAL*4 AIN(NR,NC1), AOUT(NR,NC2), ACREA(NR,NC1),A1(NC1,NC1), X A2(NC1,NC1),W1(NC1),W2(NC1) INTEGER METHOD,NR,NC1,NC2,IERR,IPRINT,I,J 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, 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,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 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) INTEGER METHOD,NR,NC1,NC2,IERR,IPRINT,I,J 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, 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. REAL*4 X(NT), TRNULL INTEGER TINULL INTEGER NULL,I,NT 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. REAL*4 VALS(M),V INTEGER I,N,M,TID2,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 = MIN0(N,M) CALL STDWRR(TID2,'EIGENVALUES',VALS,1,MM,KUNIT,ISTAT) C RETURN END