C 2 C PROGRAM 'LSAPROX' LEAST SQUARES APPROXIMATION 3 C C NUMBR ==> NUMBER OF VALUES TO BE PREDICTD C NEGATIVE SIGN ==> BACK PREDICT C NO SIGN ==> FORWARD PREDICT C 4 C 63 IMPLICIT REAL*8 (A-H,O-Z) 64 REAL * 8NORM(4,4) DIMENSION X(100),FX(100),NAME(20),C(50),W(100),FPHI(50) 66 DIMENSION AM(4,25) DIMENSION B1(2,25) COMMON SING 67 EXTERNAL WEIGHT,PHI 68 SING=1.0D-25 69 NMAX=4 MMAX=100 72 XSTOP = -9999.0D0 C 73 C READ AND PRINT JOB TITLE 74 C 75 READ 1000,NAME 76 PRINT 1010,NAME 77 C 78 C READ NUMBER OF BASE FUNCTIONS AND WEIGHT PARAMETER 79 C 80 READ 1020,N,IW,INUM 81 C 82 C READ X,F(X),W(X) 83 C 84 M = 1 10 READ 1030,X(M),FX(M),W(M) IF (X(M).LT.-8999.0D0) GOTO 20 M = M + 1 GO TO 10 20 M = M - 1 C 91 C MAKE CHECKS 92 C 93 IF(N.LT.1)GO TO 150 94 IF(N.GT.NMAX)GO TO 160 95 IF(M.GT.MMAX)GO TO 170 96 IF(N.GT.M)GO TO 180 97 C 98 C IF WEIGHT NOT ASSIGNED COMPUTE IT 99 C 100 IF(IW.EQ.0)GO TO 35 101 DO 30 I=1,M 102 30 W(I)=WEIGHT(X(I)) 103 C 104 C FORM THE NORMAL EQUATIONS 105 C 106 35 DO 50 I=1,N 107 DO 50 J=I,N 109 NORM(I,J)=0.0D0 110 DO 40 K=1,M 111 WW=W(K) 112 XX=X(K) 113 PHII=PHI(I,XX) 114 PHIJ=PHI(J,XX) 115 NORM(I,J)=NORM(I,J)+WW*PHII*PHIJ 116 40 CONTINUE 50 NORM(J,I)=NORM(I,J) 118 C 119 C INVERT THE NORMALS - ERROR EXIT FOR DETERMINANT LESS THAN SING 120 C 121 CALL CHOLD(NORM,NMAX,N,DETA,&190) 122 C 123 C C M = INV * PHI(I) C DO 54 I=1,N DO 54 L=1,M 54 AM(I,L) = 0.0D0 C DO 55 I=1,N DO 55 J=1,N DO 66 L=1,M XX = X(L) PHII = PHI(J,XX) 66 AM(I,L) = AM(I,L) + NORM(I,J) * PHII 55 CONTINUE C C *F* = M * PHI(UNKOWN) * F(KNOWN) C C PRINT INPUT DATA PRINT 1040 DO 110 I=1,M 110 PRINT 1050,X(I),FX(I),W(I) PRINT 700 700 FORMAT ('1') C K = M - 1 XBACK = X(1) - 1 XFORWD = X(M) + 1 DO 79 I=1,M SUM1 = 0.0D0 SUM2 = 0.0D0 DO 77 J=1,N SUM1 = AM(J,I) * PHI(J,XBACK) + SUM1 SUM2 = AM(J,I) * PHI(J,XFORWD) + SUM2 77 CONTINUE B1(1,I) = SUM1 B1(2,I) = SUM2 79 CONTINUE PRINT,B1 C J = 1 IF (INUM.GT.0) J = 2 KOUNT = 1 94 CONTINUE PREDIC = 0.0D0 DO 88 I=1,M 88 PREDIC = B1(J,I) * FX(I) + PREDIC IF (J.EQ.2) GO TO 323 DO 91 I=1,K 91 FPHI(I + 1) = FX(I) DO 93 I=1,K 93 FX(I + 1) = FPHI(I + 1) FX(1) = PREDIC IPRINT = -KOUNT PRINT 1035,IPRINT,PREDIC KOUNT = KOUNT + 1 IF (KOUNT.GT.IABS(INUM)) STOP GOTO 94 323 CONTINUE DO 95 I=1,K 95 FX(I) = FX(I + 1) FX(K + 1) = PREDIC IPRINT = K + KOUNT PRINT 1035,IPRINT,PREDIC KOUNT = KOUNT + 1 IF (KOUNT.GT.INUM) STOP GOTO 94 150 PRINT 1130 173 STOP 174 160 PRINT 1140,N,NMAX 175 STOP 176 170 PRINT 1150,M,MMAX 177 STOP 178 180 PRINT 1160,M,N 179 STOP 180 190 PRINT 1170 181 1000 FORMAT(20A4) 183 1010 FORMAT('1'////20X,20A4////) 184 1020 FORMAT(3I4) 1030 FORMAT(3F10.4) 1035 FORMAT(3(/),20X,'PREDICTED VALUE AT ',1X,I3 ,1X,' = ',F10.5) 1040 FORMAT(45X,'INPUT CHECK'//22X,'X',21X,'F(X)',23X,'W(X)'//) 187 1050 FORMAT('0', 5X,F20.5, 5X,F20.5, 5X,F20.5) 188 1130 FORMAT('0',5X,'ERROR - DEGREE OF FUNCTION < 1 CHECK DATA CARD 2') 197 1140 FORMAT('0',5X,'ERROR - N=',I4,5X,'NMAX=',I4,5X,'CHECK DATA CARD 2 198 1OR INCREASE NMAX AND DIMENSION') 199 1150 FORMAT('0',5X,'ERROR - M=',I4,5X,'MMAX=',I4,5X,'CHECK XSTOP OR INC 200 1REASE MMAX AND DIMENSION') 201 1160 FORMAT('0',5X,'ERROR - M=',I4,5X,'N=',I4,5X,'CHECK INPUT OR INCREA 202 1SE NO. OF DATA POINTS OR DECREASE NO. OF BASE FUNCTIONS') 203 1170 FORMAT('0',5X,'ERROR - DETERMINANT OF GRAMM''S MATRIX .LT.SING, DE 204 1CREASE SING OR CHECK BASE FUNCTIONS OR WEIGHT FUNCTION') 205 END 206 FUNCTION WEIGHT(X) REAL*8 WEIGHT,X WEIGHT = 1.0D0 RETURN END FUNCTION PHI ( I,HR ) IMPLICIT REAL*8(A-H,O-Z) PI = 3.141592653589 C OMT = (HR *PI)/12.0D0 GOTO (10,20,30,40),I 10 PHI = 1.0D0 RETURN 20 PHI = HR RETURN 30 PHI = DCOS(OMT) RETURN 40 PHI = DSIN(OMT) RETURN END SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) 220 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 221 DIMENSION A(IRDA,NA) 222 COMMON SING 223 SQRT(SUM) = DSQRT(SUM) 224 ABS(DETA) = DABS(DETA) 225 DETA = A(1,1) 226 A(1,1) = SQRT(A(1,1)) 227 IF(NA .EQ. 1) GO TO 6 228 DO 1 I = 2,NA 229 1 A(I,1) = A(I,1) / A(1,1) 230 DO 5 J = 2,NA 231 SUM = 0. 232 J1 = J - 1 233 DO 2 K = 1,J1 234 2 SUM = SUM + A(J,K) ** 2 235 DETA=DETA*(A(J,J)-SUM) 236 A(J,J) = SQRT(A(J,J) - SUM) 237 IF(J .EQ. NA) GO TO 5 238 J2 = J + 1 239 DO 4 I = J2,NA 240 SUM = 0. 241 DO 3 K = 1,J1 242 3 SUM = SUM + A(I,K) * A(J,K) 243 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 244 5 CONTINUE 245 6 IF(ABS(DETA).LT.SING)GO TO 16 246 DO 7 I = 1,NA 247 7 A(I,I) = 1. / A(I,I) 248 IF(NA .EQ. 1) GO TO 10 249 N1 = NA - 1 250 DO 9 J = 1,N1 251 J2 = J + 1 252 DO 9 I = J2,NA 253 SUM = 0. 254 I1 = I - 1 255 DO 8 K = J,I1 256 8 SUM = SUM + A(I,K) * A(K,J) 257 9 A(I,J) = - A(I,I) * SUM 258 10 DO 15 J = 1,NA 259 IF(J .EQ. 1) GO TO 12 260 J1 = J - 1 261 DO 11 I = 1,J1 262 11 A(I,J) = A(J,I) 263 12 DO 14 I = J,NA 264 SUM = 0. 265 DO 13 K = I,NA 266 13 SUM = SUM + A(K,I) * A(K,J) 267 14 A(I,J) = SUM 268 15 CONTINUE 269 RETURN 270 16 RETURN 1 271 END 272