diff --git a/fvn_test/Makefile b/fvn_test/Makefile index 97ea79d..381b4f7 100644 --- a/fvn_test/Makefile +++ b/fvn_test/Makefile @@ -6,7 +6,8 @@ test_det$(exext) test_matcon$(exext) test_matev$(exext) test_inter1d$(exext) \ test_inter2d$(exext) test_inter3d$(exext) test_akima$(exext) test_lsp$(exext) test_muller$(exext) \ test_integ$(exext) test_bsyn$(exext) test_bsjn$(exext) test_bskn$(exext) test_bsin$(exext) test_operators$(exext) test_ze1$(exext) \ test_besri$(exext) test_besrj$(exext) test_bestime$(exext) \ -test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext) +test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext) \ +test_lm$(exext) prog:$(programs) @@ -21,3 +22,8 @@ clean: $(F95) $(F95FLAGS) -c $< -o $@ $(programs): init_random_seed.o + +test_lm.o: test_lm.f90 excursion.o + +test_lm$(exext): test_lm.o excursion.o + $(LINK) $(LINKFLAGS) test_lm.o excursion.o $(LINKFVN) -o $@ diff --git a/fvn_test/excursion.f90 b/fvn_test/excursion.f90 new file mode 100644 index 0000000..db38fa6 --- /dev/null +++ b/fvn_test/excursion.f90 @@ -0,0 +1,6 @@ +module excursion + +real(8), dimension(:), pointer :: xm => NULL() +real(8), dimension(:), pointer :: ym => NULL() + +end module diff --git a/fvn_test/test_lm.f90 b/fvn_test/test_lm.f90 new file mode 100644 index 0000000..74e54f9 --- /dev/null +++ b/fvn_test/test_lm.f90 @@ -0,0 +1,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 \ No newline at end of file