-
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@79 b657c933-2333-4658-acf2-d3c7c2708721
test_lm.f90
1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
program gls
use fvn_misc
use excursion
implicit none
interface
subroutine zef(m, n, x, fvec, iflag)
integer(4), intent(in) :: m,n
real(8), dimension(:), intent(in) :: x
real(8), dimension(:), intent(inout) :: fvec
integer(4), intent(inout) :: iflag
end subroutine
end interface
integer,parameter :: npoints=13,deg=3
integer :: status,i
real(kind=dp_kind) :: xstep,xc,yc
real(kind=dp_kind) :: coeff(deg+1)
real(kind=dp_kind) :: fvec(npoints)
integer(4) :: iwa(deg+1)
integer(4) :: info
real(8) :: tol
allocate(xm(npoints),ym(npoints))
xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5 /)
open(2,file='fvn_lsp_double_mesure.dat')
open(3,file='fvn_lm.dat')
do i=1,npoints
write(2,44) xm(i),ym(i)
end do
close(2)
coeff=1.
call fvn_lm(zef,npoints,deg+1,coeff,info)
write(*,*) "info : ",info
xstep=(xm(npoints)-xm(1))/1000.
do i=1,1000
xc=xm(1)+(i-1)*xstep
yc=poly(xc,coeff)
write(3,44) xc,yc
end do
close(3)
write(*,*) "All done, plot results with gnuplot using command :"
write(*,*) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lm.dat' u 1:2 w l"
44 FORMAT(4(1X,1PE22.14))
contains
function poly(x,coeff)
use fvn_common
implicit none
real(kind=dp_kind) :: x
real(kind=dp_kind) :: coeff(deg+1)
real(kind=dp_kind) :: poly
integer :: i
poly=0.
do i=1,deg+1
poly=poly+coeff(i)*x**(i-1)
end do
end function
end program
subroutine zef(m, n, x, fvec, iflag)
use excursion
implicit none
integer(4), intent(in) :: m,n
real(8), dimension(:), intent(in) :: x
real(8), dimension(:), intent(inout) :: fvec
integer(4), intent(inout) :: iflag
integer(4) :: i
do i=1,m
fvec(i)=ym(i)-(x(1)+x(2)*xm(i)+x(3)*xm(i)**2+x(4)*xm(i)**3)
enddo
end subroutine