Commit 9726449439e0e71bc27edbd07e5fd0c84f65efcb

Authored by daniau
1 parent b40f93a8d7

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

Showing 2 changed files with 17 additions and 17 deletions Inline Diff

No preview for this file type
%\documentclass[a4paper,10pt]{article} 1 1 %\documentclass[a4paper,10pt]{article}
\documentclass[a4paper,english]{article} 2 2 \documentclass[a4paper,english]{article}
3 3
\usepackage[utf8]{inputenc} 4 4 \usepackage[utf8]{inputenc}
\usepackage{a4wide} 5 5 \usepackage{a4wide}
\usepackage{eurosym} 6 6 \usepackage{eurosym}
\usepackage{url} 7 7 \usepackage{url}
\usepackage[colorlinks=true,hyperfigures=true]{hyperref} 8 8 \usepackage[colorlinks=true,hyperfigures=true]{hyperref}
%\usepackage{aeguill} 9 9 %\usepackage{aeguill}
10 10
\usepackage{graphicx} 11 11 \usepackage{graphicx}
\usepackage{babel} 12 12 \usepackage{babel}
\makeatother 13 13 \makeatother
14 14
15 15
%opening 16 16 %opening
\title{FVN Documentation} 17 17 \title{FVN Documentation}
\author{William Daniau} 18 18 \author{William Daniau}
19 19
20 20
\begin{document} 21 21 \begin{document}
22 22
\maketitle 23 23 \maketitle
24 24
%\begin{abstract} 25 25 %\begin{abstract}
26 26
%\end{abstract} 27 27 %\end{abstract}
\tableofcontents 28 28 \tableofcontents
29 29
\section{Whatis fvn,licence,disclaimer etc} 30 30 \section{Whatis fvn,licence,disclaimer etc}
\subsection{Whatis fvn} 31 31 \subsection{Whatis fvn}
fvn is a Fortran95 mathematical module. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, complex trigonometry etc. 32 32 fvn is a Fortran95 mathematical module. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, complex trigonometry etc.
33 33
Most of the work is done by interfacing Lapack \url{http://www.netlib.org/lapack} which means that Lapack and Blas \url{http://www.netlib.org/blas} must be available on your system for linking fvn. If you use an AMD microprocessor, the good idea is to use ACML ( AMD Core Math Library \url{http://developer.amd.com/acml.jsp} which contains an optimized Blas/Lapack. Fvn also contains a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack} for performing the numerical integration tasks. 34 34 Most of the work is done by interfacing Lapack \url{http://www.netlib.org/lapack} which means that Lapack and Blas \url{http://www.netlib.org/blas} must be available on your system for linking fvn. If you use an AMD microprocessor, the good idea is to use ACML ( AMD Core Math Library \url{http://developer.amd.com/acml.jsp} which contains an optimized Blas/Lapack. Fvn also contains a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack} for performing the numerical integration tasks.
35 35
This module has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in institute Femto-ST \url{http://www.femto-st.fr/}. 36 36 This module has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in institute Femto-ST \url{http://www.femto-st.fr/}.
37 37
\subsection{Licence} 38 38 \subsection{Licence}
The licence of fvn is free. You can do whatever you want with this code as far as you credit the authors. 39 39 The licence of fvn is free. You can do whatever you want with this code as far as you credit the authors.
40 40
\subsubsection*{Authors} 41 41 \subsubsection*{Authors}
As of the day this manuel is written there's only one author of fvn :\newline 42 42 As of the day this manuel is written there's only one author of fvn :\newline
William Daniau\newline 43 43 William Daniau\newline
william.daniau@femto-st.fr\newline 44 44 william.daniau@femto-st.fr\newline
45 45
\subsection{Disclaimer} 46 46 \subsection{Disclaimer}
The usual disclaimer applied : This software is provided AS IS in the hope it will be usefull. Use it at your own risks. The authors should not be taken responsible of anything that may result by the use of this software. 47 47 The usual disclaimer applied : This software is provided AS IS in the hope it will be usefull. Use it at your own risks. The authors should not be taken responsible of anything that may result by the use of this software.
48 48
\section{Naming scheme and convention} 49 49 \section{Naming scheme and convention}
The naming scheme of the routines is as follow : 50 50 The naming scheme of the routines is as follow :
\begin{verbatim} 51 51 \begin{verbatim}
fvn_x_name() 52 52 fvn_*_name()
\end{verbatim} 53 53 \end{verbatim}
where x can be s,d,c or z. 54 54 where * can be s,d,c or z.
\begin{itemize} 55 55 \begin{itemize}
\item s is for single precision real (real,real*4,real(4),real(kind=4)) 56 56 \item s is for single precision real (real,real*4,real(4),real(kind=4))
\item d for double precision real (double precision,real*8,real(8),real(kind=8)) 57 57 \item d for double precision real (double precision,real*8,real(8),real(kind=8))
\item c for single precision complex (complex,complex*8,complex(4),complex(kind=4)) 58 58 \item c for single precision complex (complex,complex*8,complex(4),complex(kind=4))
\item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8)) 59 59 \item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8))
\end{itemize} 60 60 \end{itemize}
In the following description of subroutines parameters, input parameters are followed by (in), output parameters by (out) and parameters which are used as input and modified by the subroutine are followed by (inout). 61 61 In the following description of subroutines parameters, input parameters are followed by (in), output parameters by (out) and parameters which are used as input and modified by the subroutine are followed by (inout).
62 62
\section{Linear algebra} 63 63 \section{Linear algebra}
The linear algebra routines of fvn are an interface to lapack, which make it easier to use. 64 64 The linear algebra routines of fvn are an interface to lapack, which make it easier to use.
\subsection{Matrix inversion} 65 65 \subsection{Matrix inversion}
\begin{verbatim} 66 66 \begin{verbatim}
call fvn_x_matinv(d,a,inva,status) 67 67 call fvn_*_matinv(d,a,inva,status) (*=s,d,c,z)
\end{verbatim} 68 68 \end{verbatim}
\begin{itemize} 69 69 \begin{itemize}
\item d (in) is an integer equal to the matrix rank 70 70 \item d (in) is an integer equal to the matrix rank
\item a (in) is a matrix of type x. It will remain untouched. 71 71 \item a (in) is a matrix of type x. It will remain untouched.
\item inva (out) is a matrix of type x which contain the inverse of a at the end of the routine 72 72 \item inva (out) is a matrix of type x which contain the inverse of a at the end of the routine
\item status (out) is an integer equal to zero if something went wrong 73 73 \item status (out) is an integer equal to zero if something went wrong
\end{itemize} 74 74 \end{itemize}
75 75
\subsubsection*{Example} 76 76 \subsubsection*{Example}
\begin{verbatim} 77 77 \begin{verbatim}
program inv 78 78 program inv
use fvn 79 79 use fvn
implicit none 80 80 implicit none
81 81
real(8),dimension(3,3) :: a,inva 82 82 real(8),dimension(3,3) :: a,inva
integer :: status 83 83 integer :: status
84 84
call random_number(a) 85 85 call random_number(a)
a=a*100 86 86 a=a*100
87 87
call fvn_d_matinv(3,a,inva,status) 88 88 call fvn_d_matinv(3,a,inva,status)
write (*,*) a 89 89 write (*,*) a
write (*,*) 90 90 write (*,*)
write (*,*) inva 91 91 write (*,*) inva
write (*,*) 92 92 write (*,*)
write (*,*) matmul(a,inva) 93 93 write (*,*) matmul(a,inva)
end program 94 94 end program
\end{verbatim} 95 95 \end{verbatim}
96 96
97 97
98 98
\subsection{Matrix determinants} 99 99 \subsection{Matrix determinants}
\begin{verbatim} 100 100 \begin{verbatim}
det=fvn_x_det(d,a,status) 101 101 det=fvn_*_det(d,a,status) (*=s,d,c,z)
\end{verbatim} 102 102 \end{verbatim}
\begin{itemize} 103 103 \begin{itemize}
\item d (in) is an integer equal to the matrix rank 104 104 \item d (in) is an integer equal to the matrix rank
\item a (in) is a matrix of type x. It will remain untouched. 105 105 \item a (in) is a matrix of type x. It will remain untouched.
\item status (out) is an integer equal to zero if something went wrong 106 106 \item status (out) is an integer equal to zero if something went wrong
\end{itemize} 107 107 \end{itemize}
108 108
\subsubsection*{Example} 109 109 \subsubsection*{Example}
\begin{verbatim} 110 110 \begin{verbatim}
program det 111 111 program det
use fvn 112 112 use fvn
implicit none 113 113 implicit none
114 114
real(8),dimension(3,3) :: a 115 115 real(8),dimension(3,3) :: a
real(8) :: deta 116 116 real(8) :: deta
integer :: status 117 117 integer :: status
118 118
call random_number(a) 119 119 call random_number(a)
a=a*100 120 120 a=a*100
121 121
deta=fvn_d_det(3,a,status) 122 122 deta=fvn_d_det(3,a,status)
write (*,*) a 123 123 write (*,*) a
write (*,*) 124 124 write (*,*)
write (*,*) "Det = ",deta 125 125 write (*,*) "Det = ",deta
end program 126 126 end program
127 127
\end{verbatim} 128 128 \end{verbatim}
129 129
130 130
131 131
\subsection{Matrix condition} 132 132 \subsection{Matrix condition}
\begin{verbatim} 133 133 \begin{verbatim}
call fvn_x_matcon(d,a,rcond,status) 134 134 call fvn_*_matcon(d,a,rcond,status) (*=s,d,c,z)
\end{verbatim} 135 135 \end{verbatim}
\begin{itemize} 136 136 \begin{itemize}
\item d (in) is an integer equal to the matrix rank 137 137 \item d (in) is an integer equal to the matrix rank
\item a (in) is a matrix of type x. It will remain untouched. 138 138 \item a (in) is a matrix of type x. It will remain untouched.
\item rcond (out) is a real of same kind as matrix a, it will contain the reciprocal condition number of the matrix 139 139 \item rcond (out) is a real of same kind as matrix a, it will contain the reciprocal condition number of the matrix
\item status (out) is an integer equal to zero if something went wrong 140 140 \item status (out) is an integer equal to zero if something went wrong
\end{itemize} 141 141 \end{itemize}
142 142
The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef} 143 143 The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}
\begin{equation} 144 144 \begin{equation}
R = \frac{1}{norm(A)*norm(invA)} 145 145 R = \frac{1}{norm(A)*norm(invA)}
\label{rconddef} 146 146 \label{rconddef}
\end{equation} 147 147 \end{equation}
148 148
The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm} 149 149 The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm}
\begin{equation} 150 150 \begin{equation}
L1 = max_j ( \sum_i{\mid A(i,j)\mid} ) 151 151 L1 = max_j ( \sum_i{\mid A(i,j)\mid} )
\label{l1norm} 152 152 \label{l1norm}
\end{equation} 153 153 \end{equation}
154 154
\subsubsection*{Example} 155 155 \subsubsection*{Example}
\begin{verbatim} 156 156 \begin{verbatim}
program cond 157 157 program cond
use fvn 158 158 use fvn
implicit none 159 159 implicit none
160 160
real(8),dimension(3,3) :: a 161 161 real(8),dimension(3,3) :: a
real(8) :: rcond 162 162 real(8) :: rcond
integer :: status 163 163 integer :: status
164 164
call random_number(a) 165 165 call random_number(a)
a=a*100 166 166 a=a*100
167 167
call fvn_d_matcon(3,a,rcond,status) 168 168 call fvn_d_matcon(3,a,rcond,status)
write (*,*) a 169 169 write (*,*) a
write (*,*) 170 170 write (*,*)
write (*,*) "Cond = ",rcond 171 171 write (*,*) "Cond = ",rcond
end program 172 172 end program
173 173
\end{verbatim} 174 174 \end{verbatim}
175 175
176 176
\subsection{Eigenvalues/Eigenvectors} 177 177 \subsection{Eigenvalues/Eigenvectors}
\begin{verbatim} 178 178 \begin{verbatim}
call fvn_x_matev(d,a,evala,eveca,status) 179 179 call fvn_*_matev(d,a,evala,eveca,status) (*=s,d,c,z)
\end{verbatim} 180 180 \end{verbatim}
\begin{itemize} 181 181 \begin{itemize}
\item d (in) is an integer equal to the matrix rank 182 182 \item d (in) is an integer equal to the matrix rank
\item a (in) is a matrix of type x. It will remain untouched. 183 183 \item a (in) is a matrix of type x. It will remain untouched.
\item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a 184 184 \item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a
\item eveca (out) is a complex matrix of same kind as a. Its columns are the eigenvectors of matrix a : eveca(:,j)=jth eigenvector associated with eigenvalue evala(j). 185 185 \item eveca (out) is a complex matrix of same kind as a. Its columns are the eigenvectors of matrix a : eveca(:,j)=jth eigenvector associated with eigenvalue evala(j).
\item status (out) is an integer equal to zero if something went wrong 186 186 \item status (out) is an integer equal to zero if something went wrong
\end{itemize} 187 187 \end{itemize}
188 188
\subsubsection*{Example} 189 189 \subsubsection*{Example}
\begin{verbatim} 190 190 \begin{verbatim}
program eigen 191 191 program eigen
use fvn 192 192 use fvn
implicit none 193 193 implicit none
194 194
real(8),dimension(3,3) :: a 195 195 real(8),dimension(3,3) :: a
complex(8),dimension(3) :: evala 196 196 complex(8),dimension(3) :: evala
complex(8),dimension(3,3) :: eveca 197 197 complex(8),dimension(3,3) :: eveca
integer :: status,i,j 198 198 integer :: status,i,j
199 199
call random_number(a) 200 200 call random_number(a)
a=a*100 201 201 a=a*100
202 202
call fvn_d_matev(3,a,evala,eveca,status) 203 203 call fvn_d_matev(3,a,evala,eveca,status)
write (*,*) a 204 204 write (*,*) a
write (*,*) 205 205 write (*,*)
do i=1,3 206 206 do i=1,3
write(*,*) "Eigenvalue ",i,evala(i) 207 207 write(*,*) "Eigenvalue ",i,evala(i)
write(*,*) "Associated Eigenvector :" 208 208 write(*,*) "Associated Eigenvector :"
do j=1,3 209 209 do j=1,3
write(*,*) eveca(j,i) 210 210 write(*,*) eveca(j,i)
end do 211 211 end do
write(*,*) 212 212 write(*,*)
end do 213 213 end do
214 214
end program 215 215 end program
216 216
\end{verbatim} 217 217 \end{verbatim}
218 218
219 219
\subsection{Sparse solving} 220 220 \subsection{Sparse solving}
By interfacing Tim Davis's SuiteSparse from university of Florida \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/} which is a reference for this kind of problems, fvn provides simple subroutines for solving linear sparse systems. 221 221 By interfacing Tim Davis's SuiteSparse from university of Florida \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/} which is a reference for this kind of problems, fvn provides simple subroutines for solving linear sparse systems.
222 222
The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form. 223 223 The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.
224 224
\begin{verbatim} 225 225 \begin{verbatim}
call fvn_*_sparse_solve(n,nz,T,Ti,Tj,B,x,status) where * is either zl, zi, dl or di 226 226 call fvn_*_sparse_solve(n,nz,T,Ti,Tj,B,x,status) (*=zl,zi,dl,di)
\end{verbatim} 227 227 \end{verbatim}
\begin{itemize} 228 228 \begin{itemize}
\item For this family of subroutine the two letters (zl,zi,dl,di) decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4) 229 229 \item For this family of subroutine the two letters (zl,zi,dl,di) decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)
\item n (in) is an integer equal to the matrix rank 230 230 \item n (in) is an integer equal to the matrix rank
\item nz (in) is an integer equal to the number of non-zero elements 231 231 \item nz (in) is an integer equal to the number of non-zero elements
\item T(nz) (in) is a complex/real array containing the non-zero elements 232 232 \item T(nz) (in) is a complex/real array containing the non-zero elements
\item Ti(nz),Tj(nz) (in) are the indexes of the corresponding element of T in the original matrix. 233 233 \item Ti(nz),Tj(nz) (in) are the indexes of the corresponding element of T in the original matrix.
\item B(n) (in) is a complex/real array containing the second member of the equation. 234 234 \item B(n) (in) is a complex/real array containing the second member of the equation.
\item x(n) (out) is a complex/real array containing the solution 235 235 \item x(n) (out) is a complex/real array containing the solution
\item status (out) is an integer which contain non-zero is something went wrong 236 236 \item status (out) is an integer which contain non-zero is something went wrong
\end{itemize} 237 237 \end{itemize}
238 238
\subsubsection*{Example} 239 239 \subsubsection*{Example}
\begin{verbatim} 240 240 \begin{verbatim}
program test_sparse 241 241 program test_sparse
242 242
use fvn 243 243 use fvn
implicit none 244 244 implicit none
245 245
integer(8), parameter :: nz=12 246 246 integer(8), parameter :: nz=12
integer(8), parameter :: n=5 247 247 integer(8), parameter :: n=5
complex(8),dimension(nz) :: A 248 248 complex(8),dimension(nz) :: A
integer(8),dimension(nz) :: Ti,Tj 249 249 integer(8),dimension(nz) :: Ti,Tj
complex(8),dimension(n) :: B,x 250 250 complex(8),dimension(n) :: B,x
integer(8) :: status 251 251 integer(8) :: status
252 252
A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),& 253 253 A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&
(1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /) 254 254 (1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /)
B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/) 255 255 B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /) 256 256 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /) 257 257 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
258 258
call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status) 259 259 call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(*,*) x 260 260 write(*,*) x
261 261
end program 262 262 end program
263 263
264 264
program test_sparse 265 265 program test_sparse
266 266
use fvn 267 267 use fvn
implicit none 268 268 implicit none
269 269
integer(4), parameter :: nz=12 270 270 integer(4), parameter :: nz=12
integer(4), parameter :: n=5 271 271 integer(4), parameter :: n=5
real(8),dimension(nz) :: A 272 272 real(8),dimension(nz) :: A
integer(4),dimension(nz) :: Ti,Tj 273 273 integer(4),dimension(nz) :: Ti,Tj
real(8),dimension(n) :: B,x 274 274 real(8),dimension(n) :: B,x
integer(4) :: status 275 275 integer(4) :: status
276 276
A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /) 277 277 A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
B = (/ 8., 45., -3., 3., 19./) 278 278 B = (/ 8., 45., -3., 3., 19./)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /) 279 279 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /) 280 280 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
281 281
call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status) 282 282 call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(*,*) x 283 283 write(*,*) x
284 284
end program 285 285 end program
286 286
287 287
288 288
\end{verbatim} 289 289 \end{verbatim}
290 290
\subsection{Identity matrix} 291 291 \subsection{Identity matrix}
\begin{verbatim} 292 292 \begin{verbatim}
I=fvn_x_ident(n) 293 293 I=fvn_*_ident(n) (*=s,d,c,z)
\end{verbatim} 294 294 \end{verbatim}
\begin{itemize} 295 295 \begin{itemize}
\item n (in) is an integer equal to the matrix rank 296 296 \item n (in) is an integer equal to the matrix rank
\end{itemize} 297 297 \end{itemize}
This function return the identity matrix of rank n, in the specified type (x = s,d,c,z). 298 298 This function return the identity matrix of rank n, in the specified type (x = s,d,c,z).
299 299
300 300
301 301
\section{Interpolation} 302 302 \section{Interpolation}
303 303
\subsection{Quadratic Interpolation} 304 304 \subsection{Quadratic Interpolation}
fvn provide function for interpolating values of a tabulated function of 1, 2 or 3 variables, for both single and double precision. 305 305 fvn provide function for interpolating values of a tabulated function of 1, 2 or 3 variables, for both single and double precision.
\subsubsection{One variable function} 306 306 \subsubsection{One variable function}
\begin{verbatim} 307 307 \begin{verbatim}
value=fvn_x_quad_interpol(x,n,xdata,ydata) 308 308 value=fvn_*_quad_interpol(x,n,xdata,ydata) (*=s,d)
\end{verbatim} 309 309 \end{verbatim}
\begin{itemize} 310 310 \begin{itemize}
\item x is the real where we want to evaluate the function 311 311 \item x is the real where we want to evaluate the function
\item n is the number of tabulated values 312 312 \item n is the number of tabulated values
\item xdata(n) contains the tabulated coordinates 313 313 \item xdata(n) contains the tabulated coordinates
\item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i)) 314 314 \item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))
\end{itemize} 315 315 \end{itemize}
xdata must be strictly increasingly ordered. 316 316 xdata must be strictly increasingly ordered.
x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation 317 317 x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation
\paragraph*{Example} 318 318 \paragraph*{Example}
\begin{verbatim} 319 319 \begin{verbatim}
program inter1d 320 320 program inter1d
321 321
use fvn 322 322 use fvn
implicit none 323 323 implicit none
324 324
integer(kind=4),parameter :: ndata=33 325 325 integer(kind=4),parameter :: ndata=33
integer(kind=4) :: i,nout 326 326 integer(kind=4) :: i,nout
real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata) 327 327 real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
real(kind=8) ::tv 328 328 real(kind=8) ::tv
329 329
intrinsic sin 330 330 intrinsic sin
331 331
f(x)=sin(x) 332 332 f(x)=sin(x)
333 333
xdata(1)=0. 334 334 xdata(1)=0.
fdata(1)=f(xdata(1)) 335 335 fdata(1)=f(xdata(1))
h=1./32. 336 336 h=1./32.
do i=2,ndata 337 337 do i=2,ndata
xdata(i)=xdata(i-1)+h 338 338 xdata(i)=xdata(i-1)+h
fdata(i)=f(xdata(i)) 339 339 fdata(i)=f(xdata(i))
end do 340 340 end do
call random_seed() 341 341 call random_seed()
call random_number(x) 342 342 call random_number(x)
343 343
q=fvn_d_quad_interpol(x,ndata,xdata,fdata) 344 344 q=fvn_d_quad_interpol(x,ndata,xdata,fdata)
345 345
tv=f(x) 346 346 tv=f(x)
write(*,*) "x ",x 347 347 write(*,*) "x ",x
write(*,*) "Calculated (real) value :",tv 348 348 write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q 349 349 write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv) 350 350 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
351 351
end program 352 352 end program
353 353
\end{verbatim} 354 354 \end{verbatim}
355 355
356 356
\subsubsection{Two variables function} 357 357 \subsubsection{Two variables function}
\begin{verbatim} 358 358 \begin{verbatim}
value=fvn_x_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) 359 359 value=fvn_*_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) (*=s,d)
\end{verbatim} 360 360 \end{verbatim}
\begin{itemize} 361 361 \begin{itemize}
\item x,y are the real coordinates where we want to evaluate the function 362 362 \item x,y are the real coordinates where we want to evaluate the function
\item nx is the number of tabulated values along x axis 363 363 \item nx is the number of tabulated values along x axis
\item xdata(nx) contains the tabulated x 364 364 \item xdata(nx) contains the tabulated x
\item ny is the number of tabulated values along y axis 365 365 \item ny is the number of tabulated values along y axis
\item ydata(ny) contains the tabulated y 366 366 \item ydata(ny) contains the tabulated y
\item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j)) 367 367 \item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))
\end{itemize} 368 368 \end{itemize}
xdata and ydata must be strictly increasingly ordered. 369 369 xdata and ydata must be strictly increasingly ordered.
(x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation 370 370 (x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
371 371
\paragraph*{Example} 372 372 \paragraph*{Example}
373 373
\begin{verbatim} 374 374 \begin{verbatim}
program inter2d 375 375 program inter2d
use fvn 376 376 use fvn
implicit none 377 377 implicit none
378 378
integer(kind=4),parameter :: nx=21,ny=42 379 379 integer(kind=4),parameter :: nx=21,ny=42
integer(kind=4) :: i,j 380 380 integer(kind=4) :: i,j
real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny) 381 381 real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
real(kind=8) :: tv 382 382 real(kind=8) :: tv
383 383
intrinsic dble,sin 384 384 intrinsic dble,sin
385 385
f(x,y)=sin(x+2.*y) 386 386 f(x,y)=sin(x+2.*y)
do i=1,nx 387 387 do i=1,nx
xdata(i)=dble(i-1)/dble(nx-1) 388 388 xdata(i)=dble(i-1)/dble(nx-1)
end do 389 389 end do
do i=1,ny 390 390 do i=1,ny
ydata(i)=dble(i-1)/dble(ny-1) 391 391 ydata(i)=dble(i-1)/dble(ny-1)
end do 392 392 end do
do i=1,nx 393 393 do i=1,nx
do j=1,ny 394 394 do j=1,ny
fdata(i,j)=f(xdata(i),ydata(j)) 395 395 fdata(i,j)=f(xdata(i),ydata(j))
end do 396 396 end do
end do 397 397 end do
call random_seed() 398 398 call random_seed()
call random_number(x) 399 399 call random_number(x)
call random_number(y) 400 400 call random_number(y)
401 401
q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata) 402 402 q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
tv=f(x,y) 403 403 tv=f(x,y)
404 404
write(*,*) "x y",x,y 405 405 write(*,*) "x y",x,y
write(*,*) "Calculated (real) value :",tv 406 406 write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q 407 407 write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv) 408 408 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
409 409
end program 410 410 end program
411 411
\end{verbatim} 412 412 \end{verbatim}
413 413
414 414
415 415
\subsubsection{Three variables function} 416 416 \subsubsection{Three variables function}
\begin{verbatim} 417 417 \begin{verbatim}
value=fvn_x_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) 418 418 value=fvn_*_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) (*=s,d)
\end{verbatim} 419 419 \end{verbatim}
\begin{itemize} 420 420 \begin{itemize}
\item x,y,z are the real coordinates where we want to evaluate the function 421 421 \item x,y,z are the real coordinates where we want to evaluate the function
\item nx is the number of tabulated values along x axis 422 422 \item nx is the number of tabulated values along x axis
\item xdata(nx) contains the tabulated x 423 423 \item xdata(nx) contains the tabulated x
\item ny is the number of tabulated values along y axis 424 424 \item ny is the number of tabulated values along y axis
\item ydata(ny) contains the tabulated y 425 425 \item ydata(ny) contains the tabulated y
\item nz is the number of tabulated values along z axis 426 426 \item nz is the number of tabulated values along z axis
\item zdata(ny) contains the tabulated z 427 427 \item zdata(ny) contains the tabulated z
\item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k)) 428 428 \item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k))
\end{itemize} 429 429 \end{itemize}
xdata, ydata and zdata must be strictly increasingly ordered. 430 430 xdata, ydata and zdata must be strictly increasingly ordered.
(x,y,z) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation 431 431 (x,y,z) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
432 432
\paragraph*{Example} 433 433 \paragraph*{Example}
\begin{verbatim} 434 434 \begin{verbatim}
program inter3d 435 435 program inter3d
use fvn 436 436 use fvn
437 437
implicit none 438 438 implicit none
439 439
integer(kind=4),parameter :: nx=21,ny=42,nz=18 440 440 integer(kind=4),parameter :: nx=21,ny=42,nz=18
integer(kind=4) :: i,j,k 441 441 integer(kind=4) :: i,j,k
real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz) 442 442 real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
real(kind=8) :: tv 443 443 real(kind=8) :: tv
444 444
intrinsic dble,sin 445 445 intrinsic dble,sin
446 446
f(x,y,z)=sin(x+2.*y+3.*z) 447 447 f(x,y,z)=sin(x+2.*y+3.*z)
do i=1,nx 448 448 do i=1,nx
xdata(i)=2.*(dble(i-1)/dble(nx-1)) 449 449 xdata(i)=2.*(dble(i-1)/dble(nx-1))
end do 450 450 end do
do i=1,ny 451 451 do i=1,ny
ydata(i)=2.*(dble(i-1)/dble(ny-1)) 452 452 ydata(i)=2.*(dble(i-1)/dble(ny-1))
end do 453 453 end do
do i=1,nz 454 454 do i=1,nz
zdata(i)=2.*(dble(i-1)/dble(nz-1)) 455 455 zdata(i)=2.*(dble(i-1)/dble(nz-1))
end do 456 456 end do
do i=1,nx 457 457 do i=1,nx
do j=1,ny 458 458 do j=1,ny
do k=1,nz 459 459 do k=1,nz
fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k)) 460 460 fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
end do 461 461 end do
end do 462 462 end do
end do 463 463 end do
call random_seed() 464 464 call random_seed()
call random_number(x) 465 465 call random_number(x)
call random_number(y) 466 466 call random_number(y)
call random_number(z) 467 467 call random_number(z)
468 468
q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata) 469 469 q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
tv=f(x,y,z) 470 470 tv=f(x,y,z)
471 471
write(*,*) "x y z",x,y,z 472 472 write(*,*) "x y z",x,y,z
write(*,*) "Calculated (real) value :",tv 473 473 write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q 474 474 write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv) 475 475 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
476 476
end program 477 477 end program
478 478
\end{verbatim} 479 479 \end{verbatim}
480 480
\subsubsection{Utility procedure} 481 481 \subsubsection{Utility procedure}
fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array. 482 482 fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.
\begin{verbatim} 483 483 \begin{verbatim}
call fvn_x_find_interval(x,i,xdata,n) 484 484 call fvn_*_find_interval(x,i,xdata,n) (*=s,d)
\end{verbatim} 485 485 \end{verbatim}
\begin{itemize} 486 486 \begin{itemize}
\item x (in) the real value to locate 487 487 \item x (in) the real value to locate
\item i (out) the resulting indice 488 488 \item i (out) the resulting indice
\item xdata(n) (in) increasingly ordered array 489 489 \item xdata(n) (in) increasingly ordered array
\item n (in) size of the array 490 490 \item n (in) size of the array
\end{itemize} 491 491 \end{itemize}
The resulting integer i is as : $xdata(i) <= x < xdata(i+1)$. If $x < xdata(1)$ then $i=0$ is returned. If $x > xdata(n)$ then $i=n$ is returned. Finally if $x=xdata(n)$ then $i=n-1$ is returned. 492 492 The resulting integer i is as : $xdata(i) <= x < xdata(i+1)$. If $x < xdata(1)$ then $i=0$ is returned. If $x > xdata(n)$ then $i=n$ is returned. Finally if $x=xdata(n)$ then $i=n-1$ is returned.
493 493
494 494
495 495
\subsection{Akima spline} 496 496 \subsection{Akima spline}
fvn provides Akima spline interpolation and evaluation for both single and double precision real. 497 497 fvn provides Akima spline interpolation and evaluation for both single and double precision real.
\subsubsection{Interpolation} 498 498 \subsubsection{Interpolation}
\begin{verbatim} 499 499 \begin{verbatim}
call fvn_x_akima(n,x,y,br,co) 500 500 call fvn_*_akima(n,x,y,br,co) (*=s,d)
\end{verbatim} 501 501 \end{verbatim}
\begin{itemize} 502 502 \begin{itemize}
\item n (in) is an integer equal to the number of points 503 503 \item n (in) is an integer equal to the number of points
\item x(n) (in) ,y(n) (in) are the known couples of coordinates 504 504 \item x(n) (in) ,y(n) (in) are the known couples of coordinates
\item br (out) on output contains a copy of x 505 505 \item br (out) on output contains a copy of x
\item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval. 506 506 \item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval.
\end{itemize} 507 507 \end{itemize}
508 508
\subsubsection{Evaluation} 509 509 \subsubsection{Evaluation}
\begin{verbatim} 510 510 \begin{verbatim}
y=fvn_x_spline_eval(x,n,br,co) 511 511 y=fvn_*_spline_eval(x,n,br,co) (*=s,d)
\end{verbatim} 512 512 \end{verbatim}
\begin{itemize} 513 513 \begin{itemize}
\item x (in) is the point where we want to evaluate 514 514 \item x (in) is the point where we want to evaluate
\item n (in) is the number of known points and br(n) (in), co(4,n) (in) \\ 515 515 \item n (in) is the number of known points and br(n) (in), co(4,n) (in) \\
are the outputs of fvn\_x\_akima(n,x,y,br,co) 516 516 are the outputs of fvn\_x\_akima(n,x,y,br,co)
\end{itemize} 517 517 \end{itemize}
518 518
\subsubsection{Example} 519 519 \subsubsection{Example}
In the following example we will use Akima splines to interpolate a sinus function with 30 points between -10 and 10. We then use the evaluation function to calculate the coordinates of 1000 points between -11 and 11, and write a 3 columns file containing : x, calculated sin(x), interpolation evaluation of sin(x). 520 520 In the following example we will use Akima splines to interpolate a sinus function with 30 points between -10 and 10. We then use the evaluation function to calculate the coordinates of 1000 points between -11 and 11, and write a 3 columns file containing : x, calculated sin(x), interpolation evaluation of sin(x).
521 521
One can see that the interpolation is very efficient even with only 30 points. Of course as soon as we leave the -10 to 10 interval, the values are extrapolated and thus can lead to very inacurrate values. 522 522 One can see that the interpolation is very efficient even with only 30 points. Of course as soon as we leave the -10 to 10 interval, the values are extrapolated and thus can lead to very inacurrate values.
523 523
\begin{verbatim} 524 524 \begin{verbatim}
program akima 525 525 program akima
use fvn 526 526 use fvn
implicit none 527 527 implicit none
528 528
integer :: nbpoints,nppoints,i 529 529 integer :: nbpoints,nppoints,i
real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d 530 530 real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
real(8),dimension(:,:),allocatable :: coeff_fvn_d 531 531 real(8),dimension(:,:),allocatable :: coeff_fvn_d
real(8) :: xstep_d,xp_d,ty_d,fvn_y_d 532 532 real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
533 533
open(2,file='fvn_akima_double.dat') 534 534 open(2,file='fvn_akima_double.dat')
open(3,file='fvn_akima_breakpoints_double.dat') 535 535 open(3,file='fvn_akima_breakpoints_double.dat')
nbpoints=30 536 536 nbpoints=30
allocate(x_d(nbpoints)) 537 537 allocate(x_d(nbpoints))
allocate(y_d(nbpoints)) 538 538 allocate(y_d(nbpoints))
allocate(breakpoints_d(nbpoints)) 539 539 allocate(breakpoints_d(nbpoints))
allocate(coeff_fvn_d(4,nbpoints)) 540 540 allocate(coeff_fvn_d(4,nbpoints))
541 541
xstep_d=20./dfloat(nbpoints) 542 542 xstep_d=20./dfloat(nbpoints)
do i=1,nbpoints 543 543 do i=1,nbpoints
x_d(i)=-10.+dfloat(i)*xstep_d 544 544 x_d(i)=-10.+dfloat(i)*xstep_d
y_d(i)=dsin(x_d(i)) 545 545 y_d(i)=dsin(x_d(i))
write(3,44) (x_d(i),y_d(i)) 546 546 write(3,44) (x_d(i),y_d(i))
end do 547 547 end do
close(3) 548 548 close(3)
549 549
call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d) 550 550 call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
551 551
nppoints=1000 552 552 nppoints=1000
xstep_d=22./dfloat(nppoints) 553 553 xstep_d=22./dfloat(nppoints)
do i=1,nppoints 554 554 do i=1,nppoints
xp_d=-11.+dfloat(i)*xstep_d 555 555 xp_d=-11.+dfloat(i)*xstep_d
ty_d=dsin(xp_d) 556 556 ty_d=dsin(xp_d)
fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d) 557 557 fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
write(2,44) (xp_d,ty_d,fvn_y_d) 558 558 write(2,44) (xp_d,ty_d,fvn_y_d)
end do 559 559 end do
560 560
close(2) 561 561 close(2)
562 562
44 FORMAT(4(1X,1PE22.14)) 563 563 44 FORMAT(4(1X,1PE22.14))
564 564
end program 565 565 end program
566 566
\end{verbatim} 567 567 \end{verbatim}
Results are plotted on figure \ref{akima} 568 568 Results are plotted on figure \ref{akima}
569 569
\begin{figure} 570 570 \begin{figure}
\begin{center} 571 571 \begin{center}
\includegraphics[width=0.9\textwidth]{akima.pdf} 572 572 \includegraphics[width=0.9\textwidth]{akima.pdf}
% akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720 573 573 % akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720
\caption{Akima Spline Interpolation} 574 574 \caption{Akima Spline Interpolation}
\label{akima} 575 575 \label{akima}
\end{center} 576 576 \end{center}
577 577
\end{figure} 578 578 \end{figure}
579 579
580 580
581 581
\section{Least square polynomial} 582 582 \section{Least square polynomial}
fvn provide a function to find a least square polynomial of a given degree, for real in single or double precision. It is performed using Lapack subroutine sgelss (dgelss), which solve this problem using singular value decomposition. 583 583 fvn provide a function to find a least square polynomial of a given degree, for real in single or double precision. It is performed using Lapack subroutine sgelss (dgelss), which solve this problem using singular value decomposition.
584 584
\begin{verbatim} 585 585 \begin{verbatim}
call fvn_x_lspoly(np,x,y,deg,coeff,status) 586 586 call fvn_*_lspoly(np,x,y,deg,coeff,status) (*=s,d)
\end{verbatim} 587 587 \end{verbatim}
\begin{itemize} 588 588 \begin{itemize}
\item np (in) is an integer equal to the number of points 589 589 \item np (in) is an integer equal to the number of points
\item x(np) (in),y(np) (in) are the known coordinates 590 590 \item x(np) (in),y(np) (in) are the known coordinates
\item deg (in) is an integer equal to the degree of the desired polynomial, it must be lower than np. 591 591 \item deg (in) is an integer equal to the degree of the desired polynomial, it must be lower than np.
\item coeff(deg+1) (out) on output contains the polynomial coefficients 592 592 \item coeff(deg+1) (out) on output contains the polynomial coefficients
\item status (out) is an integer containing 0 if a problem occured. 593 593 \item status (out) is an integer containing 0 if a problem occured.
\end{itemize} 594 594 \end{itemize}
595 595
\subsection*{Example} 596 596 \subsection*{Example}
Here's a simple example : we've got 13 measurement points and we want to find the least square degree 3 polynomial for these points : 597 597 Here's a simple example : we've got 13 measurement points and we want to find the least square degree 3 polynomial for these points :
\begin{verbatim} 598 598 \begin{verbatim}
program lsp 599 599 program lsp
use fvn 600 600 use fvn
implicit none 601 601 implicit none
602 602
integer,parameter :: npoints=13,deg=3 603 603 integer,parameter :: npoints=13,deg=3
integer :: status,i 604 604 integer :: status,i
real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc 605 605 real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
real(kind=8) :: coeff(deg+1) 606 606 real(kind=8) :: coeff(deg+1)
607 607
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. /) 608 608 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 /) 609 609 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 /)
610 610
open(2,file='fvn_lsp_double_mesure.dat') 611 611 open(2,file='fvn_lsp_double_mesure.dat')
open(3,file='fvn_lsp_double_poly.dat') 612 612 open(3,file='fvn_lsp_double_poly.dat')
613 613
do i=1,npoints 614 614 do i=1,npoints
write(2,44) xm(i),ym(i) 615 615 write(2,44) xm(i),ym(i)
end do 616 616 end do
close(2) 617 617 close(2)
618 618
619 619
call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status) 620 620 call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)
621 621
xstep=(xm(npoints)-xm(1))/1000. 622 622 xstep=(xm(npoints)-xm(1))/1000.
do i=1,1000 623 623 do i=1,1000
xc=xm(1)+(i-1)*xstep 624 624 xc=xm(1)+(i-1)*xstep
yc=poly(xc,coeff) 625 625 yc=poly(xc,coeff)
write(3,44) xc,yc 626 626 write(3,44) xc,yc
end do 627 627 end do
close(3) 628 628 close(3)
629 629
44 FORMAT(4(1X,1PE22.14)) 630 630 44 FORMAT(4(1X,1PE22.14))
631 631
contains 632 632 contains
function poly(x,coeff) 633 633 function poly(x,coeff)
implicit none 634 634 implicit none
real(8) :: x 635 635 real(8) :: x
real(8) :: coeff(deg+1) 636 636 real(8) :: coeff(deg+1)
real(8) :: poly 637 637 real(8) :: poly
integer :: i 638 638 integer :: i
639 639
poly=0. 640 640 poly=0.
641 641
do i=1,deg+1 642 642 do i=1,deg+1
poly=poly+coeff(i)*x**(i-1) 643 643 poly=poly+coeff(i)*x**(i-1)
end do 644 644 end do
645 645
end function 646 646 end function
end program 647 647 end program
\end{verbatim} 648 648 \end{verbatim}
The results are plotted on figure \ref{lsp} . 649 649 The results are plotted on figure \ref{lsp} .
650 650
\begin{figure} 651 651 \begin{figure}
\begin{center} 652 652 \begin{center}
\includegraphics[width=0.9\textwidth]{lsp.pdf} 653 653 \includegraphics[width=0.9\textwidth]{lsp.pdf}
\caption{Least Square Polynomial} 654 654 \caption{Least Square Polynomial}
\label{lsp} 655 655 \label{lsp}
\end{center} 656 656 \end{center}
\end{figure} 657 657 \end{figure}
658 658
659 659
660 660
\section{Zero finding} 661 661 \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}. 662 662 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}.
663 663
\begin{verbatim} 664 664 \begin{verbatim}
call fvn_z_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 665 665 call fvn_z_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
\end{verbatim} 666 666 \end{verbatim}
\begin{itemize} 667 667 \begin{itemize}
\item f (in) is the complex function (kind=8) for which we search zeros 668 668 \item f (in) is the complex function (kind=8) for which we search zeros
\item eps (in) is a real(8) corresponding to the first stopping criterion : let fp(z)=f(z)/p where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) and z(1),...,z(k-1) are previously found roots. if ((cdabs(f(z)).le.eps) .and. (cdabs(fp(z)).le.eps)), then z is accepted as a root. 669 669 \item eps (in) is a real(8) corresponding to the first stopping criterion : let fp(z)=f(z)/p where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) and z(1),...,z(k-1) are previously found roots. if ((cdabs(f(z)).le.eps) .and. (cdabs(fp(z)).le.eps)), then z is accepted as a root.
\item eps1 (in) is a real(8) corresponding to the second stopping criterion : a root is accepted if two successive approximations to a given root agree within eps1. Note that if either or both of the stopping criteria are fulfilled, the root is accepted. 670 670 \item eps1 (in) is a real(8) corresponding to the second stopping criterion : a root is accepted if two successive approximations to a given root agree within eps1. Note that if either or both of the stopping criteria are fulfilled, the root is accepted.
\item kn (in) is an integer equal to the number of known roots, which must be stored in x(1),...,x(kn), prior to entry in the subroutine. 671 671 \item kn (in) is an integer equal to the number of known roots, which must be stored in x(1),...,x(kn), prior to entry in the subroutine.
\item nguess (in) is the number of initial guesses provided. These guesses must be stored in x(kn+1),...,x(kn+nguess). nguess must be set equal to zero if no guesses are provided. 672 672 \item nguess (in) is the number of initial guesses provided. These guesses must be stored in x(kn+1),...,x(kn+nguess). nguess must be set equal to zero if no guesses are provided.
\item n (in) is an integer equal to the number of new roots to be found. 673 673 \item n (in) is an integer equal to the number of new roots to be found.
\item x (inout) is a complex(8) vector of length kn+n. x(1),...,x(kn) on input must contain any known roots. x(kn+1),..., x(kn+n) on input may, on user option, contain initial guesses for the n new roots which are to be computed. If the user does not provide an initial guess, zero is used. On output, x(kn+1),...,x(kn+n) contain the approximate roots found by the subroutine. 674 674 \item x (inout) is a complex(8) vector of length kn+n. x(1),...,x(kn) on input must contain any known roots. x(kn+1),..., x(kn+n) on input may, on user option, contain initial guesses for the n new roots which are to be computed. If the user does not provide an initial guess, zero is used. On output, x(kn+1),...,x(kn+n) contain the approximate roots found by the subroutine.
\item itmax (in) is an integer equal to the maximum allowable number of iterations per root. 675 675 \item itmax (in) is an integer equal to the maximum allowable number of iterations per root.
\item infer (out) is an integer vector of size kn+n. On output infer(j) contains the number of iterations used in finding the j-th root when convergence was achieved. If convergence was not obtained in itmax iterations, infer(j) will be greater than itmax 676 676 \item infer (out) is an integer vector of size kn+n. On output infer(j) contains the number of iterations used in finding the j-th root when convergence was achieved. If convergence was not obtained in itmax iterations, infer(j) will be greater than itmax
\item ier (out) is an integer used as an error parameter. ier = 33 indicates failure to converge within itmax iterations for at least one of the (n) new roots. 677 677 \item ier (out) is an integer used as an error parameter. ier = 33 indicates failure to converge within itmax iterations for at least one of the (n) new roots.
\end{itemize} 678 678 \end{itemize}
This subroutine always returns the last approximation for root j in x(j). if the convergence criterion is satisfied, then infer(j) is less than or equal to itmax. if the convergence criterion is not satisified, then infer(j) is set to either itmax+1 or itmax+k, with k greater than 1. infer(j) = itmax+1 indicates that muller did not obtain convergence in the allowed number of iterations. in this case, the user may wish to set itmax to a larger value. infer(j) = itmax+k means that convergence was obtained (on iteration k) for the deflated function fp(z) = f(z)/((z-z(1)...(z-z(j-1))) but failed for f(z). in this case, better initial guesses might help or, it might be necessary to relax the convergence criterion. 679 679 This subroutine always returns the last approximation for root j in x(j). if the convergence criterion is satisfied, then infer(j) is less than or equal to itmax. if the convergence criterion is not satisified, then infer(j) is set to either itmax+1 or itmax+k, with k greater than 1. infer(j) = itmax+1 indicates that muller did not obtain convergence in the allowed number of iterations. in this case, the user may wish to set itmax to a larger value. infer(j) = itmax+k means that convergence was obtained (on iteration k) for the deflated function fp(z) = f(z)/((z-z(1)...(z-z(j-1))) but failed for f(z). in this case, better initial guesses might help or, it might be necessary to relax the convergence criterion.
680 680
\subsection*{Example} 681 681 \subsection*{Example}
Example to find the ten roots of $x^{10}-1$ 682 682 Example to find the ten roots of $x^{10}-1$
\begin{verbatim} 683 683 \begin{verbatim}
program muller 684 684 program muller
use fvn 685 685 use fvn
implicit none 686 686 implicit none
687 687
integer :: i,info 688 688 integer :: i,info
complex(8),dimension(10) :: roots 689 689 complex(8),dimension(10) :: roots
integer,dimension(10) :: infer 690 690 integer,dimension(10) :: infer
complex(8), external :: f 691 691 complex(8), external :: f
692 692
call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info) 693 693 call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
694 694
write(*,*) "Error code :",info 695 695 write(*,*) "Error code :",info
do i=1,10 696 696 do i=1,10
write(*,*) roots(i),infer(i) 697 697 write(*,*) roots(i),infer(i)
enddo 698 698 enddo
end program 699 699 end program
700 700
function f(x) 701 701 function f(x)
complex(8) :: x,f 702 702 complex(8) :: x,f
f=x**10-1 703 703 f=x**10-1
end function 704 704 end function
705 705
\end{verbatim} 706 706 \end{verbatim}
707 707
708 708
\section{Trigonometry} 709 709 \section{Trigonometry}
\subsection{Complex Sine Arc} 710 710 \subsection{Complex Sine Arc}
( only complex(kind=8) version ) 711 711 ( only complex(kind=8) version )
\begin{verbatim} 712 712 \begin{verbatim}
y=fvn_z_asin(z) 713 713 y=fvn_z_asin(z)
\end{verbatim} 714 714 \end{verbatim}
This function return the complex arc sine of z. It is adapted from he c gsl library \url{http://www.gnu.org/software/gsl/}. 715 715 This function return the complex arc sine of z. It is adapted from he c gsl library \url{http://www.gnu.org/software/gsl/}.
716 716
717 717
\subsection{Complex Cosine Arc} 718 718 \subsection{Complex Cosine Arc}
( only complex(kind=8) version ) 719 719 ( only complex(kind=8) version )
\begin{verbatim} 720 720 \begin{verbatim}
y=fvn_z_acos(z) 721 721 y=fvn_z_acos(z)
\end{verbatim} 722 722 \end{verbatim}
This function return the complex arc cosine of z. It is adapted from he c gsl library \url{http://www.gnu.org/software/gsl/}. 723 723 This function return the complex arc cosine of z. It is adapted from he c gsl library \url{http://www.gnu.org/software/gsl/}.
724 724
\subsection{Real Sine Hyperbolic Arc} 725 725 \subsection{Real Sine Hyperbolic Arc}
( only real(kind=8) version ) 726 726 ( only real(kind=8) version )
\begin{verbatim} 727 727 \begin{verbatim}
y=fvn_d_asinh(x) 728 728 y=fvn_d_asinh(x)
\end{verbatim} 729 729 \end{verbatim}
This function return the arc hyperbolic sine of x. 730 730 This function return the arc hyperbolic sine of x.
731 731
\subsection{Real Cosine Hyperbolic Arc} 732 732 \subsection{Real Cosine Hyperbolic Arc}
( only real(kind=8) version ) 733 733 ( only real(kind=8) version )
\begin{verbatim} 734 734 \begin{verbatim}
y=fvn_d_acosh(x) 735 735 y=fvn_d_acosh(x)
\end{verbatim} 736 736 \end{verbatim}
This function return the arc hyperbolic cosine of x. In the current implementation error handling is ugly... it will stop program execution if argument is lesser than one. 737 737 This function return the arc hyperbolic cosine of x. In the current implementation error handling is ugly... it will stop program execution if argument is lesser than one.
738 738
\section{Numerical integration} 739 739 \section{Numerical integration}
Using an integrated slightly modified version of quadpack \url{http://www.netlib.org/quadpack}, fvn provide adaptative numerical integration (Gauss Kronrod) of real functions of 1 and 2 variables. fvn also provide a function to calculate Gauss-Legendre abscissas and weight, and a simple non adaptative integration subroutine. All routines exists only in fvn for double precision real. 740 740 Using an integrated slightly modified version of quadpack \url{http://www.netlib.org/quadpack}, fvn provide adaptative numerical integration (Gauss Kronrod) of real functions of 1 and 2 variables. fvn also provide a function to calculate Gauss-Legendre abscissas and weight, and a simple non adaptative integration subroutine. All routines exists only in fvn for double precision real.
741 741
\subsection{Gauss Legendre Abscissas and Weigth} 742 742 \subsection{Gauss Legendre Abscissas and Weigth}