moindres_carres Subroutine

   subroutine moindres_carres(nb_var, nb_pts, hij, vec_xy, beta, f, df, typ, eps, relax, nb_var_der, info)
   !================================================================================================
   !< @note Function that returns the parameters of a function that approximates a data set. The parameters
   !<        determination is achieved by non linear least squares approximation.
   !<
   !<  @endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in   )                             :: nb_var      !! *number of parameters to be determined*
   integer(kind=I4), intent(in   )                             :: nb_pts      !! *number of points for function evaluation*
   integer(kind=I4), intent(in   )                             :: nb_var_der  !! *number of derivatives*
   real   (kind=R8), intent(in   ), dimension(1:nb_pts)        :: hij         !! *vector of evaluation points*
   real   (kind=R8), intent(in   ), dimension(1:nb_pts, 1:2)   :: vec_xy      !! *x and y coordinates of evaluation points*
   real   (kind=R8), intent(inout), dimension(1:nb_var)        :: beta        !! *parameters vector*
   real   (kind=R8), intent(in   )                             :: eps         !! *stop criterion*
   real   (kind=R8), intent(in   )                             :: relax       !! *relaxation parameter*
   character(len=*), intent(in   )                             :: typ         !! *kind of function used*
   integer(kind=I4), intent(out  )                             :: info        !! *information from Cholesky resolution*

   interface
      real(kind=R8) function f(xi, yi, var, nb_var, typ)
         use data_arch, only : I4, R8
         implicit none
         integer(kind=I4), intent(in   )                      :: nb_var
         real   (kind=R8), intent(in   )                      :: xi
         real   (kind=R8), intent(in   )                      :: yi
         real   (kind=R8), intent(inout), dimension(1:nb_var) :: var
         character(len=*), intent(in)                         :: typ
      endfunction f
      real(kind=R8) function df(xi, yi, var, nb_var, ivar, typ)
         use data_arch, only : I4, R8
         implicit none
         integer(kind=I4), intent(in   )                      :: nb_var
         integer(kind=I4), intent(in   )                      :: ivar
         real   (kind=R8), intent(in   )                      :: xi
         real   (kind=R8), intent(in   )                      :: yi
         real   (kind=R8), intent(inout), dimension(1:nb_var) :: var
         character(len=*), intent(in)                         :: typ
      endfunction df
   endinterface

      integer(kind=I4) :: i, ii, j, compteur
      integer(kind=I4) :: n, pseudo_newton
      real   (kind=R8) :: xi, yi, m_abs_scd

      real   (kind=R8),              dimension(1:nb_var_der)               :: scd
      real   (kind=R8),              dimension(1:nb_var_der, 1:nb_var_der) :: mat, sav_mat
      real   (kind=R8), allocatable, dimension(:,:)                        :: Jf, JfT
      real   (kind=R8), allocatable, dimension(:)                          :: vec_f, r, pt

      n = nint( sqrt(real(nb_pts, kind=R8)) )

      allocate(  Jf(1:nb_pts, 1:nb_var_der) )
      allocate( JfT(1:nb_var_der, 1:nb_pts) )
      allocate( vec_f(1:nb_pts) )
      allocate(     r(1:nb_pts) )
      allocate( pt(1:nb_var_der))

      compteur      = 0
      pseudo_newton = 10

w:    do
         do ii = 1, nb_pts
            xi = vec_xy(ii, 1)
            yi = vec_xy(ii, 2)
            vec_f(ii) = f(xi = xi,              & !
                          yi = yi,              & !
                         var = beta(1:nb_var),  & !
                      nb_var = nb_var,          & !
                         typ = typ)               !
         enddo

         if (mod(compteur, pseudo_newton)==0) then
            ! http://en.wikipedia.org/wiki/Least_squares
            do j = 1, nb_var_der
               do ii = 1, nb_pts
                  xi = vec_xy(ii, 1)
                  yi = vec_xy(ii, 2)
                  Jf(ii, j) = df(xi = xi,             & !
                                 yi = yi,             & !
                                var = beta(1:nb_var), & !
                             nb_var = nb_var,         & !
                               ivar = j,              & !
                                typ = typ)              !
               enddo
            enddo
            JfT(1:nb_var_der, 1:nb_pts) = transpose( Jf(1:nb_pts, 1:nb_var_der) )
            sav_mat(1:nb_var_der, 1:nb_var_der) = matmul( JfT(1:nb_var_der, 1:nb_pts), Jf(1:nb_pts, 1:nb_var_der) )
         endif
         r(1:nb_pts) = hij(1:nb_pts) -vec_f(1:nb_pts)
         scd(1:nb_var_der) = matmul( JfT(1:nb_var_der, 1:nb_pts), r(1:nb_pts) )

         mat(1:nb_var_der, 1:nb_var_der) = sav_mat(1:nb_var_der, 1:nb_var_der)
         forall(i=1:nb_var_der) mat(i,i) = 1.5*mat(i,i) ! type levenberg-marquardt

         call choldc(a = mat(1:nb_var_der, 1:nb_var_der),   & !
                     n = nb_var_der,                        & !
                    np = nb_var_der,                        & !
                     p = pt(1:nb_var_der),                  & !
                  info = info)                                !

         call cholsl(a = mat(1:nb_var_der, 1:nb_var_der),   & !
                     n = nb_var_der,                        & !
                    np = nb_var_der,                        & !
                     p = pt(1:nb_var_der),                  & !
                     b = scd(1:nb_var_der),                 & !
                     x = scd(1:nb_var_der),                 & !
                  info = info)                                !

         beta(1:nb_var_der) = beta(1:nb_var_der) +relax*scd(1:nb_var_der)
         m_abs_scd = 100*maxval( abs(scd(1:nb_var_der)/abs(beta(1:nb_var_der))) )

         if ((m_abs_scd) < eps) exit w
         compteur = compteur +1
         if (mod(compteur,1000) == 0) then
            if (verbose) write(*,*) compteur, m_abs_scd
         endif

      enddo w

      deallocate( Jf, JfT, vec_f, r, pt )

   return
   endsubroutine moindres_carres