lxsd / fvn

Download as
Browse Code ยป

Commit 31f2c259c63cf0ca76232d3a16f53b27036658d1

Authored by wdaniau
1 parent 20ce530d4f
Exists in master and in 2 other branches geevx_2020, multi

Added test example for fvn_lm 

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@79 b657c933-2333-4658-acf2-d3c7c2708721

Showing 3 changed files with 91 additions and 1 deletions Side-by-side Diff

... ... @@ -6,7 +6,8 @@
6 6 test_inter2d$(exext) test_inter3d$(exext) test_akima$(exext) test_lsp$(exext) test_muller$(exext) \
7 7 test_integ$(exext) test_bsyn$(exext) test_bsjn$(exext) test_bskn$(exext) test_bsin$(exext) test_operators$(exext) test_ze1$(exext) \
8 8 test_besri$(exext) test_besrj$(exext) test_bestime$(exext) \
9   -test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext)
  9 +test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext) \
  10 +test_lm$(exext)
10 11  
11 12 prog:$(programs)
12 13  
... ... @@ -21,4 +22,9 @@
21 22 $(F95) $(F95FLAGS) -c $< -o $@
22 23  
23 24 $(programs): init_random_seed.o
  25 +
  26 +test_lm.o: test_lm.f90 excursion.o
  27 +
  28 +test_lm$(exext): test_lm.o excursion.o
  29 + $(LINK) $(LINKFLAGS) test_lm.o excursion.o $(LINKFVN) -o $@
fvn_test/excursion.f90
Diff comments   View file @ 31f2c25
  1 +module excursion
  2 +
  3 +real(8), dimension(:), pointer :: xm => NULL()
  4 +real(8), dimension(:), pointer :: ym => NULL()
  5 +
  6 +end module
fvn_test/test_lm.f90
Diff comments   View file @ 31f2c25
  1 +program gls
  2 +use fvn_misc
  3 +use excursion
  4 +implicit none
  5 +
  6 + interface
  7 + subroutine zef(m, n, x, fvec, iflag)
  8 + integer(4), intent(in) :: m,n
  9 + real(8), dimension(:), intent(in) :: x
  10 + real(8), dimension(:), intent(inout) :: fvec
  11 + integer(4), intent(inout) :: iflag
  12 + end subroutine
  13 + end interface
  14 +
  15 +
  16 + integer,parameter :: npoints=13,deg=3
  17 + integer :: status,i
  18 + real(kind=dp_kind) :: xstep,xc,yc
  19 + real(kind=dp_kind) :: coeff(deg+1)
  20 + real(kind=dp_kind) :: fvec(npoints)
  21 + integer(4) :: iwa(deg+1)
  22 + integer(4) :: info
  23 + real(8) :: tol
  24 +
  25 + allocate(xm(npoints),ym(npoints))
  26 +
  27 + 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. /)
  28 + 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 /)
  29 + open(2,file='fvn_lsp_double_mesure.dat')
  30 + open(3,file='fvn_lm.dat')
  31 + do i=1,npoints
  32 + write(2,44) xm(i),ym(i)
  33 + end do
  34 + close(2)
  35 +
  36 + coeff=1.
  37 + call fvn_lm(zef,npoints,deg+1,coeff,info)
  38 + write(*,*) "info : ",info
  39 +
  40 + xstep=(xm(npoints)-xm(1))/1000.
  41 + do i=1,1000
  42 + xc=xm(1)+(i-1)*xstep
  43 + yc=poly(xc,coeff)
  44 + write(3,44) xc,yc
  45 + end do
  46 + close(3)
  47 +write(*,*) "All done, plot results with gnuplot using command :"
  48 +write(*,*) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lm.dat' u 1:2 w l"
  49 +44 FORMAT(4(1X,1PE22.14))
  50 +contains
  51 +function poly(x,coeff)
  52 + use fvn_common
  53 + implicit none
  54 + real(kind=dp_kind) :: x
  55 + real(kind=dp_kind) :: coeff(deg+1)
  56 + real(kind=dp_kind) :: poly
  57 + integer :: i
  58 + poly=0.
  59 + do i=1,deg+1
  60 + poly=poly+coeff(i)*x**(i-1)
  61 + end do
  62 +end function
  63 +end program
  64 +
  65 +subroutine zef(m, n, x, fvec, iflag)
  66 + use excursion
  67 + implicit none
  68 + integer(4), intent(in) :: m,n
  69 + real(8), dimension(:), intent(in) :: x
  70 + real(8), dimension(:), intent(inout) :: fvec
  71 + integer(4), intent(inout) :: iflag
  72 +
  73 + integer(4) :: i
  74 +
  75 + do i=1,m
  76 + fvec(i)=ym(i)-(x(1)+x(2)*xm(i)+x(3)*xm(i)**2+x(4)*xm(i)**3)
  77 + enddo
  78 +end subroutine