module leastsqr
	interface
		subroutine leastsqrfit(x, y, sigma, kmax, func, sol, cov)
			implicit none
			real, dimension(:), intent(in) :: x, y, sigma
			integer, intent(in) :: kmax
			real, external :: func
			real, dimension(kmax), intent(out) :: sol
			real, dimension(kmax,kmax), intent(out) :: cov
		end subroutine
	end interface
end module

subroutine leastsqrfit(x, y, sigma, kmax, func, sol, cov)
	use matrix, only : qrdec, lineqnsystem, inverse
	implicit none
	real, dimension(:), intent(in) :: x, y, sigma
	integer, intent(in) :: kmax
	real, external :: func
	real, dimension(kmax), intent(out) :: sol
	real, dimension(kmax,kmax), intent(out) :: cov
	real, dimension(size(x),kmax) :: A, Q
	real, dimension(kmax,kmax) :: R, AE, QE, RE
	real, dimension(size(x)) :: b
	integer :: i, k

	! Lav matrix A og vector b:	
	b = y/sigma
	do k = 1,kmax
		do i = 1,size(x)
			A(i,k) = func(k,x(i))/sigma(i)
		enddo
	enddo
	
	! Solve linear equation-system using QR-decomposition:
	call qrdec(A, Q, R)
	sol = lineqnsystem(Q, R, b)
	
	! Calculate covariance matrix:
	AE = matmul(transpose(A), A)
	call qrdec(AE, QE, RE)
	cov = inverse(QE, RE)
end subroutine