C @(#)tdregp.for 17.1.1.1 (ES0-DMD) 01/25/02 17:47:16 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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 16:02 - 19 NOV 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.IDENTIFICATION C C program TDREGP C C.PURPOSE C C Execute the command C REGRESSION/POL table dep_var ind_var1[,ind_var2] deg1[,deg2] [outkey] C OR C REGRESSION/POL table dep_var ind_var1[,ind_var2,...] ? [outkey] C C.KEYWORDS C C Tables, Multiple linear regression, POLYNOMIAL FIT C C.ALGORITHM C C uses table interface routines C fitting algorithm using Householder transformation C C----------------------------------------------------------- C SUBROUTINE TDREGP IMPLICIT NONE C LOGICAL NULL(9), SEL, IPOL, NEXT C INTEGER PARVAL INTEGER IVAR(9),IDEG(2),IPAR(13) INTEGER NROW,NCOL, NAC, NAR, I, ISTAT, STATUS INTEGER IAV, TID, NSC1, NN, NVAR, I1, I2, K, L, N, LL, N1 INTEGER IROW, KMIN, NPAR, NM1, K1 INTEGER INDEX, KUN, KNUL, TYPE, MAXPAR, INDKEY INTEGER NOELM,BYTELM,DEGMAX PARAMETER (MAXPAR=100) C REAL ERROR,RPAR1(4),RPAR(13) C DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX DOUBLE PRECISION VAL(9) DOUBLE PRECISION G(MAXPAR,MAXPAR),X(MAXPAR) C CHARACTER*16 MSG CHARACTER*80 TABLE,LINE,SELE CHARACTER*90 SDUM CHARACTER*17 COLUMN(9) CHARACTER*8 FORM,TYPKEY CHARACTER*20 CPAR, KEYNAM, INPNAM C INCLUDE 'MID_INCLUDE:TABLES.INC' INCLUDE 'MID_INCLUDE:TABLED.INC' DATA MSG/'ERR:TREGRMUL'/ DATA PARVAL/5/ DATA CPAR/' '/ DATA IPAR/13*0/ DATA RPAR/13*0./ C C ... GET INPUT PARAMETERS C DO 10 I = 1,9 NULL(I) = .FALSE. VAL(I) = 0 10 CONTINUE CALL TDPGET(PARVAL,NPAR,ISTAT) IF (ISTAT.NE.0) GO TO 70 IF (NPAR.LT.4) THEN ISTAT = ERRPAR GO TO 70 END IF TABLE = TPARBF(1) KEYNAM = TPARBF(5) INDKEY = INDEX(KEYNAM,' ') - 1 DO 15 I = 1, MAXPAR X(I) = 0.D0 15 CONTINUE INPNAM = KEYNAM(1:INDKEY)//'D' CALL STKFND(INPNAM,TYPKEY,NOELM,BYTELM,ISTAT) CALL STKWRD(INPNAM,X,1,NOELM,KUN,ISTAT) CALL STKRDR('INPUTR',1,2,IAV,RPAR1,KUN,KNUL,ISTAT) IDEG(1) = RPAR1(1) IDEG(2) = RPAR1(2) DEGMAX = 0 IF (IDEG(2).GT.0) DEGMAX = IDEG(1) IF (DEGMAX.GT.8) THEN CALL STETER(9,'REGRES/POLY: degree along X is limited to 8') ENDIF IF (IDEG(2).GT.8) THEN CALL STETER(9,'REGRES/POLY: degree along Y is limited to 8') ENDIF C C ... READ INPUT TABLE C CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT) CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT) CALL TDRSEL(TID,SELE,ISTAT) C C ... DECODE PARAMETERS C COLUMN(1) = TPARBF(2) LINE = ' ' LINE = TPARBF(3) NN = INDEX(LINE,' ') - 1 NVAR = 1 I1 = 1 SDUM = LINE//',' 20 CONTINUE I2 = INDEX(SDUM,',') SDUM(I2:I2) = ';' NVAR = NVAR + 1 IF (NVAR.GT.9) THEN ISTAT = ERRPAR GO TO 70 END IF COLUMN(NVAR) = LINE(I1:I2-1) I1 = I2 + 1 IF (I1.LT.NN) GO TO 20 K = IDEG(1) L = IDEG(2) IF (K.EQ.0 .AND. L.EQ.0) THEN IPOL = .FALSE. N = NVAR NM1 = N - 1 ELSE IPOL = .TRUE. N = (K+1)* (L+1) END IF IF (N.GT.NOELM) THEN ISTAT = ERRPAR GO TO 70 END IF C C ... SEARCH FOR COLUMNS AND CHECK FORMAT C DO 30 I = 1,NVAR CALL TBCSER(TID,COLUMN(I),IVAR(I),ISTAT) IF (ISTAT.NE.0) GO TO 70 IF (IVAR(I).EQ.-1) THEN ISTAT = ERRPAR GO TO 70 END IF CALL TBFGET(TID,IVAR(I),FORM,LL,TYPE,ISTAT) IF (TYPE.EQ.D_C_FORMAT) THEN ISTAT = ERRFMT GO TO 70 END IF 30 CONTINUE CALL STTPUT(TABLE,ISTAT) C C ... SOLVE LEAST SQUARE PROBLEM C XMIN = 1.D30 YMIN = 1.D30 XMAX = -XMIN YMAX = -YMIN LL = 0 N1 = 0 DO 50 IROW = 1,NROW CALL TBSGET(TID,IROW,SEL,ISTAT) IF (SEL) THEN CALL TBRRDD(TID,IROW,NVAR,IVAR,VAL,NULL,ISTAT) NEXT = ( .NOT. NULL(1)) .AND. ( .NOT. NULL(2)) .AND. + ( .NOT. NULL(3)) .AND. ( .NOT. NULL(4)) .AND. + ( .NOT. NULL(5)) .AND. ( .NOT. NULL(6)) .AND. + ( .NOT. NULL(7)) .AND. ( .NOT. NULL(8)) .AND. + ( .NOT. NULL(9)) IF (NEXT) THEN N1 = N1 + 1 IF (IPOL) THEN XMIN = MIN(XMIN,VAL(2)) XMAX = MAX(XMAX,VAL(2)) YMIN = MIN(YMIN,VAL(3)) YMAX = MAX(YMAX,VAL(3)) CALL TDSET2(LL+1,VAL(2),VAL(3),VAL(1),K,L,G,X,N, + MAXPAR) ELSE CALL TDSET1(LL+1,VAL(2),VAL(1),NM1,G,X,N,MAXPAR) END IF IF (LL.NE.0) THEN KMIN = MIN(LL,N+1) DO 40 K1 = 1,KMIN CALL TDHHTR(K1,LL+1,G,X,N,MAXPAR) 40 CONTINUE END IF LL = MIN(LL+1,N+1) END IF END IF 50 CONTINUE CALL TDSOLV(G,X,N,MAXPAR) ERROR = ABS(G(N+1,N+1)) C C ... COPY RESULT OF REGRESSION C CPAR(9:16) = TABLE IPAR(1) = N1 IPAR(2) = NVAR - 1 IPAR(3) = IVAR(1) IF (IPOL) THEN CPAR(17:20) = 'MULT' IPAR(4) = IVAR(2) IPAR(5) = IVAR(3) IPAR(6) = IDEG(1) IPAR(7) = IDEG(2) ELSE CPAR(17:20) = 'LINE' DO 60 I = 1,NVAR - 1 IPAR(I+3) = IVAR(I+1) 60 CONTINUE END IF RPAR(1) = XMIN RPAR(2) = XMAX RPAR(3) = YMIN RPAR(4) = YMAX RPAR(5) = SQRT((ERROR*ERROR)/IPAR(1)) INPNAM = KEYNAM(1:INDKEY)//'C' CALL STKWRC(INPNAM,1,CPAR,1,20,KUN,ISTAT) INPNAM = KEYNAM(1:INDKEY)//'I' CALL STKWRI(INPNAM,IPAR,1,13,KUN,ISTAT) INPNAM = KEYNAM(1:INDKEY)//'R' CALL STKWRR(INPNAM,RPAR,1,5,KUN,ISTAT) INPNAM = KEYNAM(1:INDKEY)//'D' CALL STKWRD(INPNAM,X,1,N,KUN,ISTAT) C C ... PRINT RESULT C IF (SELE(1:1).NE.'-') THEN LINE = ' SELECT '//SELE CALL STTPUT(LINE,ISTAT) END IF RPAR(5) = ERROR CALL TDRDIS(CPAR,IPAR,RPAR,X,ISTAT) IF (ISTAT.NE.0) GO TO 70 C C ... END C CALL DSCUPT(TID,TID,' ',ISTAT) CALL TBTCLO(TID,ISTAT) 70 IF (ISTAT.NE.0) THEN WRITE (MSG(13:16),9000) ISTAT CALL TDERRR(ISTAT,MSG,STATUS) END IF RETURN 9000 FORMAT (I4) END