lxsd / fvn

Blame view

fvn_test/test_lm.f90 1.9 KB
edit raw normal view history
31f2c259c   wdaniau   Added test exampl... 
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