program minq
!	programma generale per il metodo dei minimi quadrati
!   NB: la subroutine CALFUN che descrive la funzione del
!	    modello va cambiata per ogni problema  
implicit none
integer, parameter		:: n=5
integer, parameter		:: m=3
real, dimension(n)		:: x, y
real, dimension(m,m)	:: a, c, d
real, dimension(m)		:: b, u, diag, h, fi
integer					:: i

data x/1.,2.,3.,4.,5./y/1.,4.,9.,16.,25./
!     data x/1.,2.,3.,4.,5./y/1.1,4.1,9.1,16.1,25.1/
!     data x/1.,2.,3.,4.,5./y/1.,4.1,8.9,16.1,25./
!
write(*,*) 'dati in ingresso'
do i=1,n
        write(*,*) x(i),y(i)
enddo
!
call matdef(a,b,fi,m,x,y,n)
write(*,*) 'A e B definiscono le equazioni normali'
call matwrt(a,b,m)
!
call chodec(a,diag,m)
call chosol(a,b,u,diag,m)
write(*,*) 'vettore soluzione C'
call vecwrt(u,m)
!
stop
end program

subroutine matdef(a,b,fi,m,x,y,n)
! 
!      calcolo della matrice normale e del vettore noto
!      utilizzando le misure e le funzioni del modello
!
implicit none
integer					:: n, m, i, j, k
real, dimension(n)		:: x, y
real, dimension(m,m)	:: a
real, dimension(m)		:: b, fi
! acc a zero
do j=1,m
   b(j) = 0.
   do k=1,m
      a(j,k) = 0.
   enddo
enddo
!
do i=1,n
  call calfun(x(i),fi,m)
  do j=1,m
    b(j) = b(j) + y(i)*fi(j)
    do k=1,m
      a(j,k) = a(j,k) + fi(j)*fi(k)
    enddo
  enddo
enddo
!
return
end subroutine


subroutine calfun(xi,fi,m)
!
!      valutazione delle funzioni del modello nel nodo xi
!      Come esempio si mostra un modello polinomiale f(1)=1,
!      f(2)=x, ..., f(M)=x^(M-1).
!      Adattare al proprio modello.
!
implicit none
integer				:: m, i
real				:: xi
real, dimension(m)	:: fi       ! real :: fi(m)
       
do i=1,m
  fi(i) = xi**(i-1)
enddo
return
end subroutine

subroutine chodec(a,diag,n)
implicit none
integer					:: n, i, j, k
real					:: sum
real, dimension(n,n)	:: a
real, dimension(n)		:: diag
!
do i=1,n
!
  do j=i,n
!
    sum = a(i,j)
    do k = i-1,1,-1
      sum = sum - a(i,k) * a(j,k)
    enddo
    if (i .eq. j) then
      if (sum .le. 0.) pause 'matrice non simetrica e positiva'
      diag(i) = sqrt(sum)
    else
      a(j,i) = sum / diag(i)
    endif
!
  enddo
!
enddo
!
return
end subroutine

subroutine chosol(a,b,x,diag,n)
implicit none
integer					:: n, i, k
real					:: sum
real, dimension(n,n)	:: a
real, dimension(n)		:: b, diag, x
!
do i=1,n
  sum = b(i)
  do k=i-1,1,-1
    sum = sum - a(i,k) * x(k)
  enddo 
  x(i) = sum / diag(i)
enddo
!
do i=n,1,-1
  sum = x(i)
  do k = i+1,n
    sum = sum - a(k,i) * x(k)
  enddo
  x(i) = sum / diag(i)
enddo
!
return
end subroutine

subroutine matwrt(a,b,n)
implicit none
integer					:: n, i, j
real, dimension(n,n)	:: a
real, dimension(n)		:: b

do i=1,n
  write(*,'(1X,5f12.4)') (a(i,j), j=1,n)	! esempio di do implicito
enddo 
do i=1,n
  write(*,'(1X,f12.4)') b(i)
enddo
return
end subroutine

subroutine vecwrt(b,n)
implicit none
integer				:: n, i
real, dimension(n)	:: b
do i=1,n
  write(*,'(1X,f12.4)') b(i)
enddo
return
end subroutine
