Anonymous user
Polynomial regression: Difference between revisions
quick and dirty f90 code
m (Moved to Math cat) |
(quick and dirty f90 code) |
||
Line 11:
This task is intended as a subtask for [[Measure relative performance of sorting algorithms implementations]].
=={{header|Fortran}}==
{{libheader|LAPACK}}
<pre>module fitting
contains
function polyfit(vx, vy, d)
implicit none
integer, intent(in) :: d
real(8), dimension(d+1) :: polyfit
real(8), dimension(:), intent(in) :: vx, vy
real(8), dimension(:,:), allocatable :: X
real(8), dimension(:,:), allocatable :: XT
real(8), dimension(:,:), allocatable :: XTX
integer :: i, j
integer :: n, lda, lwork
integer :: info
integer, dimension(:), allocatable :: ipiv
real(8), dimension(:), allocatable :: work
n = d+1
lda = n
lwork = n
allocate(ipiv(n))
allocate(work(lwork))
allocate(XT(n, size(vx)))
allocate(X(size(vx), n))
allocate(XTX(n, n))
! prepare the matrix
do i = 0, d
do j = 1, size(vx)
X(j, i+1) = vx(j)**i
end do
end do
XT = transpose(X)
XTX = matmul(XT, X)
! calls to LAPACK subs DGETRF and DGETRI
call DGETRF(n, n, XTX, lda, ipiv, info)
if ( info /= 0 ) then
print *, "problem"
return
end if
call DGETRI(n, XTX, lda, ipiv, work, lwork, info)
if ( info /= 0 ) then
print *, "problem"
return
end if
polyfit = matmul( matmul(XTX, XT), vy)
deallocate(ipiv)
deallocate(work)
deallocate(X)
deallocate(XT)
deallocate(XTX)
end function
end module</pre>
===Example===
<pre>program PolynomalFitting
use fitting
implicit none
! let us test it
integer :: i
real(8), dimension(11) :: x = (/ (i,i=0,10) /)
real(8), dimension(11) :: y = (/ 1, 6, 17, 34, &
57, 86, 121, 162, &
209, 262, 321 /)
real(8), dimension(3) :: a
a = polyfit(x, y, 2)
print *,a
end program</pre>
Output (lower powers first, so this seems the opposite of the Python output):
<pre>
0.999999999999813 2.00000000000002 2.99999999999998
</pre>
=={{header|Python}}==
| |||