C @(#)mdaif.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 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