parameter (npmax=10, nxmax=10, ndmax= 1000) ! nxmax - max dimension of x vector ! npmax - max no of parameters ! ndmax - max data points that will be used character*80 dummy dimension x(ndmax, nxmax), f(ndmax, npmax), ft(npmax, ndmax) dimension a(npmax, npmax), y(ndmax), rhs(npmax), par(npmax) dimension xav(nxmax), fit_y(ndmax), index(npmax,3) open(20,file='fit.inp') READ(20,*) np, nx, nd DO i=1, nd READ(20,*) (x(i,j),j=1,nx), y(i) END DO READ(20,*) ishift ! ishift = 0 : Do not shift origin ! ishift = 1 : shift origin to x_av, y_av CLOSE(20) OPEN(30,file='fit.out') WRITE(30,120) nx, np, nd 120 format('No. of independant x , nx =', I3/, & 'No. of parametres in fit, np =', I3/, & 'No. of data sets , nd =', I3/) WRITE(30,115) DO j=1, nd WRITE(30,161) y(j), (x(j,i), i=1, nx) END DO 115 format(/' y_value | X-value -->') !________________________________________________________ ! Usually one would read in x and y from some file. ! one would also expect 'setf' routine that sets up ! F matrix from x. ! It is good practice to deal with as small values for ! x & y as possible. Hence IF(ishift.EQ.1) THEN CALL shift(x, y, ndmax, nd, nx, xav, yav) WRITE(30,116) yav, (xav(i), i=1, nx) WRITE(30,117) DO j=1, nd WRITE(30,161) y(j), (x(j,i), i=1, nx) END DO 117 format(/' y_shift | X_shift -->') 116 format(/'Y_av =', f10.3, 3x, 'X_av =', 5f9.3/) ELSE IF(ishift.EQ.0)THEN yav=0. DO i=1, nx xav(i) = 0.0 END DO ELSE WRITE(*,*) 'ishift cannot take a value =', ishift STOP END IF CALL setf(x, ndmax, nd, nx, f, np) CALL f_ft (f, ndmax, npmax, nd, np, ft) CALL mataxb (ft, f, a, npmax, ndmax, np, nd, np) CALL mataxb (ft, y, rhs, npmax, ndmax, np, nd, 1) !________________________________________________ WRITE(30,*) WRITE(30,*)'X-Matrix : {{X(i,j), j=1, np}, i=1, nd}' ! ft is actually F-Transpose do i=1, nd WRITE(30,160)(f(i,j), j=1,np) END DO WRITE(30,*) WRITE(30,*)'X^T * X' do j=1, np WRITE(30,160)(a(i,j), i=1,np) END DO WRITE(30,*) WRITE(30,*) 'X^T * y (i); i=1, np' WRITE(30,140) (rhs(i), i=1, np) 140 format(8f12.4) !________________________________________________ CALL matinv(a,npmax,index,np) WRITE(30,*) WRITE(30,*)'(X^T * X)^-1 ' do j=1, np WRITE(30,180)(a(i,j), i=1,np) 180 format(10e14.4) END DO CALL mataxb (a, rhs, par, npmax, npmax, np, np, 1) !________________________________________________ WRITE(30,*) WRITE(30,150) (par(i), i=1, np) 150 FORMAT(' Fitted parametrs'/ 10f8.4) open(40,file='repeat',status='old') do i=1, 1000 read(40,*,end=50) dummy end do 50 write(40,51) (par(i),i=1,np) 51 format(10f16.6) !________________________________________________ CALL mataxb (f, par, fit_y, ndmax, npmax, nd, np, 1) !________________________________________________ IF(ishift.EQ.1) THEN WRITE(30,155) ELSE WRITE(30,156) END IF 155 FORMAT(/'Actual and Fitted Values - Shifted'// & ' y_sh_actual y_sh_fitted | X_sh_values -->'/) 156 FORMAT(/'Actual and fitted values'// & /' y_actual y_fitted | X-values -->'/) avgsqr = 0. avgabs = 0. difmin= abs(y(1) - fit_y(1)) difmax= abs(y(1) - fit_y(1)) if(ishift.eq.0) yav=0 ymin = y(1)+yav ymax = y(1)+yav if(ishift.eq.1) WRITE(30,158) DO i=1, nd if(ishift.eq.0)WRITE(30,160) y(i), fit_y(i), (x(i,j), j=1, nx) if(ishift.eq.1)& WRITE(30,160)y(i)+yav,fit_y(i)+yav,(x(i,j)+xav(j),j=1,nx) diff= ( y(i) - fit_y(i) ) avgsqr = avgsqr+diff**2 avgabs = avgabs+abs(diff) if(abs(diff).lt.difmin)difmin=abs(diff) if(abs(diff).gt.difmax)difmax=abs(diff) if(y(i)+yav.lt.ymin) ymin=y(i)+yav if(y(i)+yav.gt.ymax) ymax=y(i)+yav END DO avgsqr= avgsqr / nd avgabs= avgabs / nd ravgsq=sqrt(avgsqr) WRITE(30, 170) avgsqr, ravgsq, avgabs, difmin,difmax,ymin,ymax 158 FORMAT(/'Fit was made by shifting the origin. Fit details '/& &'details after shifting origin back'/) 170 FORMAT(/'Average sqare Error in the fit =', E14.4/& &'Root average square error in the fit =', E14.4/& &'Average absolute error in the fit =', E14.4/& &'Minimum absolute error in the fit =', E14.4/& &'Maximum absoulte error in the fit =', E14.4/& &'Minimum y value =', E14.4/& &'Maximum y value =', E14.4/) 160 FORMAT(2f12.4, 8f10.3) 161 FORMAT( f12.4, 8f10.3) CLOSE(30) stop end !_______ SUBROUTINE shift(x, y, ndmax, nd, nx, xav, yav) DIMENSION x(ndmax, nx), y(ndmax), xav(nx) yav = 0. do j=1, nx xav(j) = 0. do i=1, nd xav(j) = xav(j) + x(i,j) IF(j.EQ.1) yav=yav+y(i) END DO xav(j) = xav(j)/nd END DO yav=yav/nd do j=1, nx do i=1, nd x(i,j) = x(i,j) - xav(j) IF(j.EQ.1) y(i) = y(i) - yav END DO END DO RETURN END !_________ SUBROUTINE f_ft (f, ndmax, npmax, nd, np, ft) DIMENSION ft(npmax, ndmax), f(ndmax, npmax) DO j=1, np DO i=1, nd ft(j,i) = f(i,j) END DO END DO RETURN END !________ subroutine fta(ft, nd, npmax, np, a, y, rhs) dimension ft(npmax, nd), a(npmax, npmax), y(nd), rhs(np) do i=1, np rhs(i) = 0. do j=1, np a(i,j) = 0. do k=1, nd a(i,j) = a(i,j) + ft(i,k) * ft (j,k) if(j.eq.1) rhs(i) = rhs(i) + ft(i,k)* y(k) END DO END DO END DO return end SUBROUTINE setf(x, ndmax, nd, nx, f, np) DIMENSION x(ndmax, nx), f(ndmax, np) DO i=1, nd f(i, 1) =1. f(i, 2) =x(i,1) !f(i, 3) =x(i,2) !f(i, 4) =x(i,2)*x(i,1) !f(i, 5) =x(i,1)**2 !f(i, 6) =x(i,2)**2 !f(i, 4) =x(i,3) !f(i, 5) =x(i,1)*x(i,2) !f(i, 6) =x(i,2)*x(i,3) !f(i, 7) =x(i,3)*x(i,1) !f(i, 8) =x(i,1)**2 !f(i, 9) =x(i,2)**2 !f(i, 10)=x(i,3)**2 END DO RETURN END SUBROUTINE matinv(vn,nmax,index,n) DIMENSION vn(nmax,nmax),index(nmax,3) EQUIVALENCE(irow,jrow),(icolum,jcolum),(amax,t,swap) DO j=1,n index(j,3)=0 END DO DO 130 i=1,n amax=0. DO 60 j=1,n IF (index(j,3)-1) 20,60,20 20 do 50 k=1,n IF (index(k,3)-1) 30,50,170 30 IF (amax-ABS(vn(j,k))) 40,50,50 40 irow=j icolum=k amax=ABS(vn(j,k)) 50 CONTINUE 60 CONTINUE index(icolum,3)=index(icolum,3)+1 index(i,1)=irow index(i,2)=icolum IF (irow-icolum) 70,90,70 70 determ=-determ DO 80 l=1,n swap=vn(irow,l) vn(irow,l)=vn(icolum,l) 80 vn(icolum,l)=swap 90 pivot=vn(icolum,icolum) vn(icolum,icolum)=1. DO 100 l=1,n 100 vn(icolum,l)=vn(icolum,l)/pivot DO 130 l1=1,n IF (l1-icolum) 110,130,110 110 t=vn(l1,icolum) vn(l1,icolum)=0. DO 120 l=1,n 120 vn(l1,l)=vn(l1,l)-vn(icolum,l)*t 130 CONTINUE DO 160 i=1,n l=n+1-i IF (index(l,1)-index(l,2)) 140,160,140 140 jrow=index(l,1) jcolum=index(l,2) DO 150 k=1,n swap=vn(k,jrow) vn(k,jrow)=vn(k,jcolum) vn(k,jcolum)=swap 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 RETURN END SUBROUTINE mataxb (a, b, c, iamax, ibmax, ia, ja, jb) DIMENSION a(iamax,iamax), b(ibmax,ibmax), c(iamax,ibmax) DO i=1, ia DO j=1, jb c(i,j) = 0. END DO END DO DO i=1,ia DO j=1,jb DO k=1, ja c(i,j) = c(i,j) + a(i,k) * b(k, j) END DO END DO END DO RETURN END