C 2 C PROGRAM 'LSAPROX' LEAST SQUARES APPROXIMATION 3 C 4 C INPUT: CARD 1: JOB TITLE FORMAT(20A4) 5 C CARD 2: N=NUMBER OF BASE FUNCTIONS 6 C IW=WEIGHTING CODE ==> IW=0 READ WEIGHT ON CARD 3 7 C IW=1 COMPUTE WEIGHT IN FUNCTION 8 C WEIGHT 9 C FORMAT(2I4) 10 C CARD TYPE 3: X - VALUE OF BASE POINTS 11 C FUNCTIONAL VALUES OF X 12 C WEIGHT OF X (IF IW=1 LEAVE BLANK) 13 C FORMAT(3F20.0) 14 C NOTE: CARD TYPE 3 IS REPEATED ONCE FOR EACH DATA POINT 15 C CARD TYPE 4: MARKS END OF CARD TYPE 3 -- TYPE -99999.0 IN 16 C COLUMNS 1 - 8 17 C CARD TYPE 5: VALUES OF 'X' AT WHICH APPROXIMATED VALUE IS 18 C DESIRED 19 C NOTE: IF ONLY COEFFICIENTS ARE DESIRED CARD TYPE 5 IS NOT 20 C REQUIRED 21 C 22 C NOTES: MAXIMUM DEGREE OF APPROXIMATION FUNCTION IS 50 - TO INCREASE 23 C CHANGE 'NMAX' AND DIMENSION STATEMENT 24 C MAXIMUM NUMBER OF DATA POINTS IS 100 - TO INCREASE CHANGE 25 C 'MMAX' AND DIMENSION STATEMENT 26 C 27 C THE FOLLOWING CHECKS ARE PROVIDED 28 C (1) N.LT.1 29 C (2) N.GT.NMAX 30 C (3) M.GT.MMAX 31 C (4) N.GT.M 32 C (5) DETERMINANT OF GRAMM'S MATRIX .LT.SING IE. BASE FUNCTIONS ARE 33 C NOT LINEARLY INDEPENDENT 34 C 35 C THE USER MUST SUPPLY TWO DOUBLE PRECISION FUNCTIONS 36 C (1) FUNCTION WEIGHT -- DESCRIBES THE WEIGHT 37 C (2) FUNCTION PHI -- DESCRIBES THE BASE FUNCTIONS 38 C 39 C THE FOLLOWING MINIMUM REQUIREMENTS ARE SUPPLIED 40 C 41 C FUNCTION WEIGHT(X) 42 C REAL*8 WEIGHT,X 43 C WEIGHT=1.0D0 44 C RETURN 45 C END 46 C 47 C FUNCTION PHI(I,X) 48 C REAL*8 PHI,X 49 C PHI=1.0D0 50 C RETURN 51 C END 52 C 53 C IF THESE ARE NOT SUFFICIENT ADD 54 C 55 C (1) WEIGHT=F(X) 56 C (2) IF(I.EQ.1)PHI=F(X) 57 C IF(I.EQ.2)PHI=F(X) 58 C ETC 59 C 60 C C. CHAMBERLAIN MARCH 1976 61 C 62 C 63 IMPLICIT REAL*8 (A-H,O-Z) 64 REAL*8 NORM(50,50) 65 DIMENSION X(100),FX(100),NAME(20),C(50),W(100),FPHI(50) 66 DIMENSION EX(100) COMMON SING 67 EXTERNAL WEIGHT,PHI 68 SING=1.0D-25 69 NMAX=50 71 MMAX=100 72 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 81 C 82 C READ X,F(X),W(X) 83 C 84 M = 1 10 READ 1030,X(M),FX(M),W(M) X(M) = X(M) - 1.0D0 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 FPHI(I)=0.0D0 108 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 IF(I.EQ.J)FPHI(I)=FPHI(I)+WW*PHII*FX(K) 117 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 COMPUTE THE COEFFICIENTS C(I) 124 C 125 DO 60 I=1,N 126 C(I)=0.0D0 127 DO 60 J=1,N 128 60 C(I)=C(I)+NORM(J,I)*FPHI(J) 129 C 130 C COMPUTE LEAST SQUARES DISTANCE,VARIANCE AND STANDARD DEVIATION 131 C 132 RHO=0.0D0 133 DO 80 I=1,M 134 WW=W(I) 135 XX=X(I) 136 PN=0.0D0 137 DO 70 J=1,N 138 70 PN=PN+C(J)*PHI(J,XX) 139 80 RHO=RHO+WW*(FX(I)-PN)**2 140 MN=M-N 141 IF(MN.EQ.0)GO TO 90 142 VAR=RHO/DFLOAT(MN) 143 STD=DSQRT(VAR) 144 TVAR = NORM(2,2) * VAR TSTD = DSQRT (TVAR) TSTD = TSTD * X(N) * 0.3048D0 TRD = C(2) * X(N) * 0.3048D0 INTER=0 145 GO TO 100 146 90 VAR=0.0D0 147 STD=0.0D0 148 INTER=1 149 C 150 C PRINT ALL RESULTS 151 C 152 100 CONTINUE 153 PRINT 1040 154 DO 110 I=1,M 155 110 PRINT 1050,X(I),FX(I),W(I) 156 PRINT 1060 157 DO 120 I=1,N 158 120 PRINT 1070,I,C(I) 159 IF(INTER.EQ.0)PRINT 1080,RHO,VAR,STD 160 IF(INTER.EQ.1)PRINT 1090 161 C C READ VALUES AT WHICH PN(X) IS TO BE EVALUATED AND EVALUATE C PRINT 1100 125 READ(5,1110,END=140)VAL VALUE = 0.0D0 DO 130 I=1,N 130 VALUE = VALUE + C(I) * PHI(I,VAL) PRINT 1120,VAL,VALUE GO TO 125 140 STOP 172 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 STOP 182 1000 FORMAT(20A4) 183 1010 FORMAT('1'////20X,20A4////) 184 1020 FORMAT(2I4) 185 1030 FORMAT (3F10.3) 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 1060 FORMAT('-',40X,'COMPUTED COEFFICIENTS') 189 1070 FORMAT('0',30X,'C(',I2,')',10X,F20.10) 190 1080 FORMAT('-',10X,'RHO=',F20.10,10X,'VAR=',F20.10,10X,'STD=',F20.10) 191 1090 FORMAT('-',10X,'M=N : RHO=VAR=STD= 0.0D0') 192 1100 FORMAT('1',20X,'APPROXIMATIONS USING ABOVE COEFFICIENTS AND BASE F 193 1UNCTIONS IN FUNCTION PHI'///39X,'X',32X,'PN(X)') 194 1110 FORMAT(F20.0) 195 1120 FORMAT('0',25X,F20.10,15X,F20.10) 196 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 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