Added documentation for fvn_lm

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

Added documentation for fvn_lm
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@80 b657c933-2333-4658-acf2-d3c7c2708721
wdaniau
1 parent 31f2c259c6
Showing 2 changed files with 169 additions and 4 deletions Side-by-side Diff
doc/fvn.pdf
doc/fvn.tex
@@ -295,7 +295,7 @@
  \item \texttt{n} (in) is an integer equal to the matrix rank
  \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
  \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
- \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
+ \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
  \item \texttt{B(n) or Bx(n),Bz(n)} (in) is a real array or two real arrays for real/imaginary in case of a complex system) containing the second member of the equation.
  \item \texttt{x(n)} (out) is a complex/real array containing the solution
  \item \texttt{status} (out) is an integer which contain non-zero is something went wrong
@@ -387,7 +387,7 @@
  \item \texttt{n} (in) is an integer equal to the matrix rank
  \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
  \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
- \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
+ \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
  \item \texttt{det} (out), a real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
 \item \texttt{status} (out) is an integer which contain non-zero is something went wrong
 \end{itemize}
@@ -913,7 +913,7 @@
 \end{verbatim}
 The results are plotted on figure \ref{lsp} .
  
-\begin{figure}
+\begin{figure}[ht]
  \begin{center}
  \includegraphics[width=0.9\textwidth]{lsp.pdf}
  \caption{Least Square Polynomial}
  
  
  
@@ -921,10 +921,175 @@
  \end{center}
 \end{figure}
  
+\section{General least square fitting}
+fvn provides a routine performing a general least square fitting using Levenberg-Marquardt algorithm. It uses routine \verb'lmdif' from \verb'MINPACK' (\url{http://www.netlib.org/minpack}). The used version has been converted to fortran90 by Alan Miller \verb'amiller @ bigpond.net.au'.
+\newline
+The purpose of \verb'fvn_lm' is to minimize the sum of the squares of m nonlinear
+functions in n variables by a modification of the Levenberg-Marquardt
+algorithm.  This is done by using the more general least-squares
+solver lmdif.  The user must provide a subroutine which calculates the
+functions.  The jacobian is then calculated by a forward-difference
+approximation.
  
+\begin{verbatim}
+Module : use fvn_misc
+call fvn_lm(fcn,m,n,a,info,tol)
+\end{verbatim}
+\begin{itemize}
+\item fcn is the user-supplied subroutine which calculates
+        the functions.  fcn must follow the following interface that must
+        be declared in the calling subroutine :
+\begin{verbatim}
+         interface
+             subroutine fcn(m,n,a,fvec,iflag)
+                 use fvn_common
+                 integer(ip_kind), intent(in) :: m
+                 integer(ip_kind), intent(in) :: n
+                 real(dp_kind), dimension(:), intent(in) :: a
+                 real(dp_kind), dimension(:), intent(inout) :: fvec
+                 integer(ip_kind), intent(inout) :: iflag
+             end subroutine
+         end interface
+ 
+\end{verbatim}
  
+          This is the function which calculate the differences for which which square sum 
+          will be minimized outputing this difference in vector fvec.
+           Parameters of fcn are as follows :
+\begin{itemize}
+ \item m : positive integer input variable set to the number of functions
+               (number of measurement points)
+ \item n : positive integer input variable set to the number of variables
+               (number of parameters in the function to fit)
+ \item a : vector of length n containing parameters for which fcn should
+               perform the calculation
+ \item fvec : vector of length m containing the resulting evaluation
+ \item iflag : integer normally not used, can be used to exit the
+               the algorithm by setting it to a negative value
+\end{itemize}
+
+\item m : positive integer input variable set to the number of functions
+               (number of measurement points)
+\item n : positive integer input variable set to the number of variables
+               (number of parameters in the function to fit)
+\item a(n) : vector of length n, on input must contains an initial guess (or unity vector)
+           and on output the solution vector
+\item info : is an output positive integer
+\begin{itemize}
+ \item info = 0  improper input parameters.
+ \item info = 1  algorithm estimates that the relative error
+                 in the sum of squares is at most tol.
+ \item info = 2  algorithm estimates that the relative error
+                between x and the solution is at most tol.
+ \item info = 3  conditions for info = 1 and info = 2 both hold.
+ \item info = 4  fvec is orthogonal to the columns of the
+                jacobian to machine precision.
+ \item info = 5  number of calls to fcn has reached or exceeded 200*(n+1).
+ \item info = 6  tol is too small. no further reduction in
+                the sum of squares is possible.
+ \item info = 7  tol is too small.  No further improvement in
+                the approximate solution x is possible.
+\end{itemize}
+         
+\item tol : is an optional positive value. Termination occurs when the
+      algorithm estimates either that the relative error in the sum of
+      squares is at most tol or that the relative error between x and the
+      solution is at most tol. If not provided default value is :
+               sqrt(epsilon(0.0d0))
+\end{itemize}
+
+\subsection*{Example}
+Here's a simple example solving the same problem as for the least square polynomial example but using \verb'fvn_lm' instead :
+\begin{verbatim}
+module excursion
+real(8), dimension(:), pointer :: xm => NULL()
+real(8), dimension(:), pointer :: ym => NULL()
+end module
+
+
+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
+\end{verbatim}
+Note the need to use a supplementary module \verb'excursion', to allow the subroutine \verb'zef' to have access to the measurement points.
+
 \section{Zero finding}
-fvn provide a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.
+fvn provides a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.
  
 \begin{verbatim}
  Module : use fvn_misc
...	...	@@ -295,7 +295,7 @@
295	295	\item \texttt{n} (in) is an integer equal to the matrix rank
296	296	\item \texttt{nz} (in) is an integer equal to the number of non-zero elements
297	297	\item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
298		- \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
	298	+ \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
299	299	\item \texttt{B(n) or Bx(n),Bz(n)} (in) is a real array or two real arrays for real/imaginary in case of a complex system) containing the second member of the equation.
300	300	\item \texttt{x(n)} (out) is a complex/real array containing the solution
301	301	\item \texttt{status} (out) is an integer which contain non-zero is something went wrong
...	...	@@ -387,7 +387,7 @@
387	387	\item \texttt{n} (in) is an integer equal to the matrix rank
388	388	\item \texttt{nz} (in) is an integer equal to the number of non-zero elements
389	389	\item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
390		- \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
	390	+ \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
391	391	\item \texttt{det} (out), a real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
392	392	\item \texttt{status} (out) is an integer which contain non-zero is something went wrong
393	393	\end{itemize}
...	...	@@ -913,7 +913,7 @@
913	913	\end{verbatim}
914	914	The results are plotted on figure \ref{lsp} .
915	915
916		-\begin{figure}
	916	+\begin{figure}[ht]
917	917	\begin{center}
918	918	\includegraphics[width=0.9\textwidth]{lsp.pdf}
919	919	\caption{Least Square Polynomial}
920	920
921	921
922	922
...	...	@@ -921,10 +921,175 @@
921	921	\end{center}
922	922	\end{figure}
923	923
	924	+\section{General least square fitting}
	925	+fvn provides a routine performing a general least square fitting using Levenberg-Marquardt algorithm. It uses routine \verb'lmdif' from \verb'MINPACK' (\url{http://www.netlib.org/minpack}). The used version has been converted to fortran90 by Alan Miller \verb'amiller @ bigpond.net.au'.
	926	+\newline
	927	+The purpose of \verb'fvn_lm' is to minimize the sum of the squares of m nonlinear
	928	+functions in n variables by a modification of the Levenberg-Marquardt
	929	+algorithm. This is done by using the more general least-squares
	930	+solver lmdif. The user must provide a subroutine which calculates the
	931	+functions. The jacobian is then calculated by a forward-difference
	932	+approximation.
924	933
	934	+\begin{verbatim}
	935	+Module : use fvn_misc
	936	+call fvn_lm(fcn,m,n,a,info,tol)
	937	+\end{verbatim}
	938	+\begin{itemize}
	939	+\item fcn is the user-supplied subroutine which calculates
	940	+ the functions. fcn must follow the following interface that must
	941	+ be declared in the calling subroutine :
	942	+\begin{verbatim}
	943	+ interface
	944	+ subroutine fcn(m,n,a,fvec,iflag)
	945	+ use fvn_common
	946	+ integer(ip_kind), intent(in) :: m
	947	+ integer(ip_kind), intent(in) :: n
	948	+ real(dp_kind), dimension(:), intent(in) :: a
	949	+ real(dp_kind), dimension(:), intent(inout) :: fvec
	950	+ integer(ip_kind), intent(inout) :: iflag
	951	+ end subroutine
	952	+ end interface
	953	+
	954	+\end{verbatim}
925	955
	956	+ This is the function which calculate the differences for which which square sum
	957	+ will be minimized outputing this difference in vector fvec.
	958	+ Parameters of fcn are as follows :
	959	+\begin{itemize}
	960	+ \item m : positive integer input variable set to the number of functions
	961	+ (number of measurement points)
	962	+ \item n : positive integer input variable set to the number of variables
	963	+ (number of parameters in the function to fit)
	964	+ \item a : vector of length n containing parameters for which fcn should
	965	+ perform the calculation
	966	+ \item fvec : vector of length m containing the resulting evaluation
	967	+ \item iflag : integer normally not used, can be used to exit the
	968	+ the algorithm by setting it to a negative value
	969	+\end{itemize}
	970	+
	971	+\item m : positive integer input variable set to the number of functions
	972	+ (number of measurement points)
	973	+\item n : positive integer input variable set to the number of variables
	974	+ (number of parameters in the function to fit)
	975	+\item a(n) : vector of length n, on input must contains an initial guess (or unity vector)
	976	+ and on output the solution vector
	977	+\item info : is an output positive integer
	978	+\begin{itemize}
	979	+ \item info = 0 improper input parameters.
	980	+ \item info = 1 algorithm estimates that the relative error
	981	+ in the sum of squares is at most tol.
	982	+ \item info = 2 algorithm estimates that the relative error
	983	+ between x and the solution is at most tol.
	984	+ \item info = 3 conditions for info = 1 and info = 2 both hold.
	985	+ \item info = 4 fvec is orthogonal to the columns of the
	986	+ jacobian to machine precision.
	987	+ \item info = 5 number of calls to fcn has reached or exceeded 200*(n+1).
	988	+ \item info = 6 tol is too small. no further reduction in
	989	+ the sum of squares is possible.
	990	+ \item info = 7 tol is too small. No further improvement in
	991	+ the approximate solution x is possible.
	992	+\end{itemize}
	993	+
	994	+\item tol : is an optional positive value. Termination occurs when the
	995	+ algorithm estimates either that the relative error in the sum of
	996	+ squares is at most tol or that the relative error between x and the
	997	+ solution is at most tol. If not provided default value is :
	998	+ sqrt(epsilon(0.0d0))
	999	+\end{itemize}
	1000	+
	1001	+\subsection*{Example}
	1002	+Here's a simple example solving the same problem as for the least square polynomial example but using \verb'fvn_lm' instead :
	1003	+\begin{verbatim}
	1004	+module excursion
	1005	+real(8), dimension(:), pointer :: xm => NULL()
	1006	+real(8), dimension(:), pointer :: ym => NULL()
	1007	+end module
	1008	+
	1009	+
	1010	+program gls
	1011	+use fvn_misc
	1012	+use excursion
	1013	+implicit none
	1014	+
	1015	+ interface
	1016	+ subroutine zef(m, n, x, fvec, iflag)
	1017	+ integer(4), intent(in) :: m,n
	1018	+ real(8), dimension(:), intent(in) :: x
	1019	+ real(8), dimension(:), intent(inout) :: fvec
	1020	+ integer(4), intent(inout) :: iflag
	1021	+ end subroutine
	1022	+ end interface
	1023	+
	1024	+
	1025	+ integer,parameter :: npoints=13,deg=3
	1026	+ integer :: status,i
	1027	+ real(kind=dp_kind) :: xstep,xc,yc
	1028	+ real(kind=dp_kind) :: coeff(deg+1)
	1029	+ real(kind=dp_kind) :: fvec(npoints)
	1030	+ integer(4) :: iwa(deg+1)
	1031	+ integer(4) :: info
	1032	+ real(8) :: tol
	1033	+
	1034	+ allocate(xm(npoints),ym(npoints))
	1035	+
	1036	+ 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. /)
	1037	+ 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 /)
	1038	+ open(2,file='fvn_lsp_double_mesure.dat')
	1039	+ open(3,file='fvn_lm.dat')
	1040	+ do i=1,npoints
	1041	+ write(2,44) xm(i),ym(i)
	1042	+ end do
	1043	+ close(2)
	1044	+
	1045	+ coeff=1.
	1046	+ call fvn_lm(zef,npoints,deg+1,coeff,info)
	1047	+ write(,) "info : ",info
	1048	+
	1049	+ xstep=(xm(npoints)-xm(1))/1000.
	1050	+ do i=1,1000
	1051	+ xc=xm(1)+(i-1)*xstep
	1052	+ yc=poly(xc,coeff)
	1053	+ write(3,44) xc,yc
	1054	+ end do
	1055	+ close(3)
	1056	+write(,) "All done, plot results with gnuplot using command :"
	1057	+write(,) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lm.dat' u 1:2 w l"
	1058	+44 FORMAT(4(1X,1PE22.14))
	1059	+contains
	1060	+function poly(x,coeff)
	1061	+ use fvn_common
	1062	+ implicit none
	1063	+ real(kind=dp_kind) :: x
	1064	+ real(kind=dp_kind) :: coeff(deg+1)
	1065	+ real(kind=dp_kind) :: poly
	1066	+ integer :: i
	1067	+ poly=0.
	1068	+ do i=1,deg+1
	1069	+ poly=poly+coeff(i)x*(i-1)
	1070	+ end do
	1071	+end function
	1072	+end program
	1073	+
	1074	+subroutine zef(m, n, x, fvec, iflag)
	1075	+ use excursion
	1076	+ implicit none
	1077	+ integer(4), intent(in) :: m,n
	1078	+ real(8), dimension(:), intent(in) :: x
	1079	+ real(8), dimension(:), intent(inout) :: fvec
	1080	+ integer(4), intent(inout) :: iflag
	1081	+
	1082	+ integer(4) :: i
	1083	+
	1084	+ do i=1,m
	1085	+ fvec(i)=ym(i)-(x(1)+x(2)xm(i)+x(3)xm(i)*2+x(4)xm(i)**3)
	1086	+ enddo
	1087	+end subroutine
	1088	+\end{verbatim}
	1089	+Note the need to use a supplementary module \verb'excursion', to allow the subroutine \verb'zef' to have access to the measurement points.
	1090	+
926	1091	\section{Zero finding}
927		-fvn provide a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.
	1092	+fvn provides a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.
928	1093
929	1094	\begin{verbatim}
930	1095	Module : use fvn_misc