lxsd / fvn | GitLab

Browse Code »

Commit 2919a9e2dae4c0b916072c9e7f155e93b02fa905

Authored by daniau 2008-10-01 17:08:24 +0200

1 parent 410c0dfaf8

Exists in master and in 3 other branches

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

Showing 15 changed files with 1377 additions and 1290 deletions Inline Diff

doc/fvn.pdf
doc/fvn.tex
fvn_test/test_akima.f90
fvn_test/test_det.f90
fvn_test/test_integ.f90
fvn_test/test_inter1d.f90
fvn_test/test_inter2d.f90
fvn_test/test_inter3d.f90
fvn_test/test_lsp.f90
fvn_test/test_matcon.f90
fvn_test/test_matev.f90
fvn_test/test_matinv.f90
fvn_test/test_muller.f90
fvn_test/test_operators.f90
fvn_test/test_sparse.f90

doc/fvn.pdf

Diff comments View file @ 2919a9e

No preview for this file type

doc/fvn.tex

Diff comments View file @ 2919a9e

%\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,license}	30	30	\section{Whatis fvn,license}
\subsection{Whatis fvn}	31	31	\subsection{Whatis fvn}
fvn is a Fortran95 mathematical library/module. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, special functions etc.	32	32	fvn is a Fortran95 mathematical library with several modules. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, special functions etc.
	33	33
Most of the work for linear algebra 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.	34	34	Most of the work for linear algebra 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.
	35	35
fvn include some integrated libraries : integration tasks uses a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack}, the fnlib library \url{http://www.netlib.org/fn} is used for special functions and sparse system resolution uses SuiteSparse \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/}.	36	36	fvn include some integrated libraries : integration tasks uses a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack}, the fnlib library \url{http://www.netlib.org/fn} is used for special functions and sparse system resolution uses SuiteSparse \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/}.
	37	37
This library has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in the Time and Frequency Department of institute Femto-ST \url{http://www.femto-st.fr/}.	38	38	This library has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in the Time and Frequency Department of institute Femto-ST \url{http://www.femto-st.fr/}.
	39	39
\subsection{License}	40	40	\subsection{License}
Your use or distribution of fvn or any modified version of	41	41	Your use or distribution of fvn or any modified version of
fvn implies that you agree to this License.	42	42	fvn implies that you agree to this License.
	43	43
This library is free software; you can redistribute it and/or	44	44	This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public	45	45	modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either	46	46	License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.	47	47	version 3 of the License, or (at your option) any later version.
	48	48
This library is distributed in the hope that it will be useful,	49	49	This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of	50	50	but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU	51	51	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.	52	52	Lesser General Public License for more details.
	53	53
You should have received a copy of the GNU Lesser General Public	54	54	You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software	55	55	License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301	56	56	Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
USA	57	57	USA
	58	58
Permission is hereby granted to use or copy this program under the	59	59	Permission is hereby granted to use or copy this program under the
terms of the GNU LGPL, provided that the Copyright, this License,	60	60	terms of the GNU LGPL, provided that the Copyright, this License,
and the Availability of the original version is retained on all copies.	61	61	and the Availability of the original version is retained on all copies.
User documentation of any code that uses this code or any modified	62	62	User documentation of any code that uses this code or any modified
version of this code must cite the Copyright, this License, the	63	63	version of this code must cite the Copyright, this License, the
Availability note, and "Used by permission." Permission to modify	64	64	Availability note, and "Used by permission." Permission to modify
the code and to distribute modified code is granted, provided the	65	65	the code and to distribute modified code is granted, provided the
Copyright, this License, and the Availability note are retained,	66	66	Copyright, this License, and the Availability note are retained,
and a notice that the code was modified is included.	67	67	and a notice that the code was modified is included.
	68	68
\subsubsection*{Authors}	69	69	\subsubsection*{Authors}
As of the day this manuel is written there's only one author of fvn :\newline	70	70	As of the day this manuel is written there's only one author of fvn :\newline
William Daniau\newline	71	71	William Daniau\newline
william.daniau@femto-st.fr\newline	72	72	william.daniau@femto-st.fr\newline
	73	73
\section{Naming scheme and convention}	74	74	\section{Naming scheme and convention}
The naming scheme of the routines is as follow :	75	75	The naming scheme of the routines is as follow :
\begin{verbatim}	76	76	\begin{verbatim}
fvn_*_name()	77	77	fvn_*_name()
\end{verbatim}	78	78	\end{verbatim}
where * can be s,d,c or z.	79	79	where * can be s,d,c or z.
\begin{itemize}	80	80	\begin{itemize}
\item s is for single precision real (real,real*4,real(4),real(kind=4))	81	81	\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))	82	82	\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))	83	83	\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))	84	84	\item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8))
\end{itemize}	85	85	\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).	86	86	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).
	87	87
For each routine, there is a generic interface (simply remove \verb'_*' in the name), so using the specific routine is not mandatory.	88	88	For each routine, there is a generic interface (simply remove \verb'_*' in the name), so using the specific routine is not mandatory.
	89	89
		90
		91	There's a general module called ``\verb'fvn''' that include all fvn submodules. So whatever part of the library is used in the program a ``\verb'use fvn''' will be sufficient instead of specifying the specific module.
		92
\section{Linear algebra}	90	93	\section{Linear algebra}
The linear algebra routines of fvn are an interface to lapack, which make it easier to use.	91	94	The linear algebra routines of fvn are an interface to lapack, which make it easier to use.
\subsection{Matrix inversion}	92	95	\subsection{Matrix inversion}
\begin{verbatim}	93	96	\begin{verbatim}
		97	Module : use fvn_linear
call fvn_matinv(d,a,inva,status)	94	98	call fvn_matinv(d,a,inva,status)
\end{verbatim}	95	99	\end{verbatim}
\begin{itemize}	96	100	\begin{itemize}
\item d (in) is an integer equal to the matrix rank	97	101	\item d (in) is an integer equal to the matrix rank
\item a (in) is a real or complex matrix. It will remain untouched.	98	102	\item a (in) is a real or complex matrix. It will remain untouched.
\item inva (out) is a real or complex matrix which contain the inverse of a at the end of the routine	99	103	\item inva (out) is a real or complex matrix which contain the inverse of a at the end of the routine
\item status (out) is an optional integer equal to zero if something went wrong	100	104	\item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}	101	105	\end{itemize}
	102	106
\subsubsection*{Example}	103	107	\subsubsection*{Example}
\begin{verbatim}	104	108	\begin{verbatim}
program inv	105	109	program inv
use fvn	106	110	use fvn_linear
implicit none	107	111	implicit none
	108	112
real(8),dimension(3,3) :: a,inva	109	113	real(8),dimension(3,3) :: a,inva
	110	114
call random_number(a)	111	115	call random_number(a)
a=a*100	112	116	a=a*100
	113	117
call fvn_matinv(3,a,inva)	114	118	call fvn_matinv(3,a,inva)
write (,) a	115	119	write (,) a
write (,)	116	120	write (,)
write (,) inva	117	121	write (,) inva
write (,)	118	122	write (,)
write (,) matmul(a,inva)	119	123	write (,) matmul(a,inva)
end program	120	124	end program
\end{verbatim}	121	125	\end{verbatim}
	122	126
	123	127
	124	128
\subsection{Matrix determinants}	125	129	\subsection{Matrix determinants}
\begin{verbatim}	126	130	\begin{verbatim}
		131	Module : use fvn_linear
det=fvn_det(d,a,status)	127	132	det=fvn_det(d,a,status)
\end{verbatim}	128	133	\end{verbatim}
\begin{itemize}	129	134	\begin{itemize}
\item d (in) is an integer equal to the matrix rank	130	135	\item d (in) is an integer equal to the matrix rank
\item a (in) is a real or complex matrix. It will remain untouched.	131	136	\item a (in) is a real or complex matrix. It will remain untouched.
\item status (out) is an optional integer equal to zero if something went wrong	132	137	\item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}	133	138	\end{itemize}
	134	139
\subsubsection*{Example}	135	140	\subsubsection*{Example}
\begin{verbatim}	136	141	\begin{verbatim}
program det	137	142	program det
use fvn	138	143	use fvn_linear
implicit none	139	144	implicit none
	140	145
real(8),dimension(3,3) :: a	141	146	real(8),dimension(3,3) :: a
real(8) :: deta	142	147	real(8) :: deta
integer :: status	143	148	integer :: status
	144	149
call random_number(a)	145	150	call random_number(a)
a=a*100	146	151	a=a*100
	147	152
deta=fvn_det(3,a,status)	148	153	deta=fvn_det(3,a,status)
write (,) a	149	154	write (,) a
write (,)	150	155	write (,)
write (,) "Det = ",deta	151	156	write (,) "Det = ",deta
end program	152	157	end program
	153	158
\end{verbatim}	154	159	\end{verbatim}
	155	160
	156	161
	157	162
\subsection{Matrix condition}	158	163	\subsection{Matrix condition}
\begin{verbatim}	159	164	\begin{verbatim}
		165	Module : use fvn_linear
call fvn_matcon(d,a,rcond,status)	160	166	call fvn_matcon(d,a,rcond,status)
\end{verbatim}	161	167	\end{verbatim}
\begin{itemize}	162	168	\begin{itemize}
\item d (in) is an integer equal to the matrix rank	163	169	\item d (in) is an integer equal to the matrix rank
\item a (in) is a real or complex matrix. It will remain untouched.	164	170	\item a (in) is a real or complex matrix. 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	165	171	\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 optional integer equal to zero if something went wrong	166	172	\item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}	167	173	\end{itemize}
	168	174
The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}	169	175	The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}
\begin{equation}	170	176	\begin{equation}
R = \frac{1}{norm(A)*norm(invA)}	171	177	R = \frac{1}{norm(A)*norm(invA)}
\label{rconddef}	172	178	\label{rconddef}
\end{equation}	173	179	\end{equation}
	174	180
The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm}	175	181	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}	176	182	\begin{equation}
L1 = max_j ( \sum_i{\mid A(i,j)\mid} )	177	183	L1 = max_j ( \sum_i{\mid A(i,j)\mid} )
\label{l1norm}	178	184	\label{l1norm}
\end{equation}	179	185	\end{equation}
	180	186
\subsubsection*{Example}	181	187	\subsubsection*{Example}
\begin{verbatim}	182	188	\begin{verbatim}
program cond	183	189	program cond
use fvn	184	190	use fvn_linear
implicit none	185	191	implicit none
	186	192
real(8),dimension(3,3) :: a	187	193	real(8),dimension(3,3) :: a
real(8) :: rcond	188	194	real(8) :: rcond
integer :: status	189	195	integer :: status
	190	196
call random_number(a)	191	197	call random_number(a)
a=a*100	192	198	a=a*100
	193	199
call fvn_d_matcon(3,a,rcond,status)	194	200	call fvn_d_matcon(3,a,rcond,status)
write (,) a	195	201	write (,) a
write (,)	196	202	write (,)
write (,) "Cond = ",rcond	197	203	write (,) "Cond = ",rcond
end program	198	204	end program
	199	205
\end{verbatim}	200	206	\end{verbatim}
	201	207
	202	208
\subsection{Eigenvalues/Eigenvectors}	203	209	\subsection{Eigenvalues/Eigenvectors}
\begin{verbatim}	204	210	\begin{verbatim}
		211	Module : use fvn_linear
call fvn_matev(d,a,evala,eveca,status)	205	212	call fvn_matev(d,a,evala,eveca,status)
\end{verbatim}	206	213	\end{verbatim}
\begin{itemize}	207	214	\begin{itemize}
\item d (in) is an integer equal to the matrix rank	208	215	\item d (in) is an integer equal to the matrix rank
\item a (in) is a real or complex matrix. It will remain untouched.	209	216	\item a (in) is a real or complex matrix. It will remain untouched.
\item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a	210	217	\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).	211	218	\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 optional integer equal to zero if something went wrong	212	219	\item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}	213	220	\end{itemize}
	214	221
		222
\subsubsection*{Example}	215	223	\subsubsection*{Example}
\begin{verbatim}	216	224	\begin{verbatim}
program eigen	217	225	program eigen
use fvn	218	226	use fvn_linear
implicit none	219	227	implicit none
	220	228
real(8),dimension(3,3) :: a	221	229	real(8),dimension(3,3) :: a
complex(8),dimension(3) :: evala	222	230	complex(8),dimension(3) :: evala
complex(8),dimension(3,3) :: eveca	223	231	complex(8),dimension(3,3) :: eveca
integer :: status,i,j	224	232	integer :: status,i,j
	225	233
call random_number(a)	226	234	call random_number(a)
a=a*100	227	235	a=a*100
	228	236
call fvn_matev(3,a,evala,eveca,status)	229	237	call fvn_matev(3,a,evala,eveca,status)
write (,) a	230	238	write (,) a
write (,)	231	239	write (,)
do i=1,3	232	240	do i=1,3
write(,) "Eigenvalue ",i,evala(i)	233	241	write(,) "Eigenvalue ",i,evala(i)
write(,) "Associated Eigenvector :"	234	242	write(,) "Associated Eigenvector :"
do j=1,3	235	243	do j=1,3
write(,) eveca(j,i)	236	244	write(,) eveca(j,i)
end do	237	245	end do
write(,)	238	246	write(,)
end do	239	247	end do
	240	248
end program	241	249	end program
	242	250
\end{verbatim}	243	251	\end{verbatim}
	244	252
	245	253
\subsection{Sparse solving}	246	254	\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.	247	255	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.
	248	256
The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.	249	257	The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.
	250	258
\begin{verbatim}	251	259	\begin{verbatim}
		260	Module : fvn_sparse
call fvn_sparse_solve(n,nz,T,Ti,Tj,B,x,status)	252	261	call fvn_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
\end{verbatim}	253	262	\end{verbatim}
\begin{itemize}	254	263	\begin{itemize}
\item For this family of subroutine the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)	255	264	\item For this family of subroutine the two letters (zl,zi,dl,di) of the specific interface name 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	256	265	\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	257	266	\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	258	267	\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.	259	268	\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.	260	269	\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	261	270	\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	262	271	\item status (out) is an integer which contain non-zero is something went wrong
\end{itemize}	263	272	\end{itemize}
	264	273
		274
\subsubsection*{Example}	265	275	\subsubsection*{Example}
\begin{verbatim}	266	276	\begin{verbatim}
program test_sparse	267	277	program test_sparse
	268	278
use fvn	269	279	use fvn_sparse
implicit none	270	280	implicit none
	271	281
integer(8), parameter :: nz=12	272	282	integer(8), parameter :: nz=12
integer(8), parameter :: n=5	273	283	integer(8), parameter :: n=5
complex(8),dimension(nz) :: A	274	284	complex(8),dimension(nz) :: A
integer(8),dimension(nz) :: Ti,Tj	275	285	integer(8),dimension(nz) :: Ti,Tj
complex(8),dimension(n) :: B,x	276	286	complex(8),dimension(n) :: B,x
integer(8) :: status	277	287	integer(8) :: status
	278	288
A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&	279	289	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.) /)	280	290	(1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /)
B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)	281	291	B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)	282	292	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 /)	283	293	Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	284	294
!specific routine that will be used here	285	295	!specific routine that will be used here
!call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	286	296	!call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	287	297	call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(,) x	288	298	write(,) x
	289	299
end program	290	300	end program
	291	301
	292	302
program test_sparse	293	303	program test_sparse
	294	304
use fvn	295	305	use fvn_sparse
implicit none	296	306	implicit none
	297	307
integer(4), parameter :: nz=12	298	308	integer(4), parameter :: nz=12
integer(4), parameter :: n=5	299	309	integer(4), parameter :: n=5
real(8),dimension(nz) :: A	300	310	real(8),dimension(nz) :: A
integer(4),dimension(nz) :: Ti,Tj	301	311	integer(4),dimension(nz) :: Ti,Tj
real(8),dimension(n) :: B,x	302	312	real(8),dimension(n) :: B,x
integer(4) :: status	303	313	integer(4) :: status
	304	314
A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)	305	315	A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
B = (/ 8., 45., -3., 3., 19./)	306	316	B = (/ 8., 45., -3., 3., 19./)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)	307	317	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 /)	308	318	Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	309	319
!specific routine that will be used here	310	320	!specific routine that will be used here
!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	311	321	!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	312	322	call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(,) x	313	323	write(,) x
	314	324
end program	315	325	end program
	316	326
	317	327
	318	328
\end{verbatim}	319	329	\end{verbatim}
	320	330
\subsection{Identity matrix}	321	331	\subsection{Identity matrix}
\begin{verbatim}	322	332	\begin{verbatim}
		333	Module : use fvn_linear
I=fvn__ident(n) (=s,d,c,z)	323	334	I=fvn__ident(n) (=s,d,c,z)
\end{verbatim}	324	335	\end{verbatim}
\begin{itemize}	325	336	\begin{itemize}
\item n (in) is an integer equal to the matrix rank	326	337	\item n (in) is an integer equal to the matrix rank
\end{itemize}	327	338	\end{itemize}
		339
This function return the identity matrix of rank n, in the specified type. No generic interface for this one.	328	340	This function return the identity matrix of rank n, in the specified type. No generic interface for this one.
	329	341
\subsection{Operators}	330	342	\subsection{Operators}
fvn defines some linear operators similar to those defined in IMSL\textregistered (\url{http://www.vni.com/products/imsl/}), that can be used for matrix operations.	331	343	fvn defines some linear operators similar to those defined in IMSL\textregistered (\url{http://www.vni.com/products/imsl/}), that can be used for matrix operations.
		344	\begin{verbatim}
		345	Module : use fvn_linear
		346	\end{verbatim}
	332	347
\subsubsection{Unary operators}	333	348	\subsubsection{Unary operators}
\paragraph{.i.}	334	349	\paragraph{.i.}
This operator gives the inverse matrix of the argument which must be a square matrix. The status of the operation can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).	335	350	This operator gives the inverse matrix of the argument which must be a square matrix. The status of the operation can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[	336	351	\[
\textrm{b=.i.a}~\Longleftrightarrow~b=a^{-1}	337	352	\textrm{b=.i.a}~\Longleftrightarrow~b=a^{-1}
\]	338	353	\]
	339	354
\paragraph{.t.}	340	355	\paragraph{.t.}
This operator gives the transpose matrix of the argument.	341	356	This operator gives the transpose matrix of the argument.
\[	342	357	\[
\textrm{b=.t.a}~\Longleftrightarrow~b= {}^t \! a	343	358	\textrm{b=.t.a}~\Longleftrightarrow~b= {}^t \! a
\]	344	359	\]
	345	360
	346	361
\paragraph{.h.}	347	362	\paragraph{.h.}
This operator gives the conjugate transpose matrix of the argument (also called Hermitian transpose or adjoint matrix).	348	363	This operator gives the conjugate transpose matrix of the argument (also called Hermitian transpose or adjoint matrix).
\[	349	364	\[
\textrm{b=.h.a}~\Longleftrightarrow~b=a^*=\overline{({}^t\!a)}={}^t \overline{a}	350	365	\textrm{b=.h.a}~\Longleftrightarrow~b=a^*=\overline{({}^t\!a)}={}^t \overline{a}
\]	351	366	\]
	352	367
	353	368
\subsubsection{Binary operators}	354	369	\subsubsection{Binary operators}
\paragraph{.x.}	355	370	\paragraph{.x.}
This operator gives the matrix product of the two operands.	356	371	This operator gives the matrix product of the two operands.
\[	357	372	\[
\textrm{c=a.x.b}~\Longleftrightarrow~c=ab	358	373	\textrm{c=a.x.b}~\Longleftrightarrow~c=ab
\]	359	374	\]
	360	375
	361	376
\paragraph{.ix.}	362	377	\paragraph{.ix.}
This operator gives the matrix product of the inverse of the first operand and the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).	363	378	This operator gives the matrix product of the inverse of the first operand and the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[	364	379	\[
\textrm{c=a.ix.b}~\Longleftrightarrow~c=a^{-1}b	365	380	\textrm{c=a.ix.b}~\Longleftrightarrow~c=a^{-1}b
\]	366	381	\]
	367	382
\paragraph{.xi.}	368	383	\paragraph{.xi.}
This operator gives the matrix product of the first operand and the inverse of the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).	369	384	This operator gives the matrix product of the first operand and the inverse of the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[	370	385	\[
\textrm{c=a.xi.b}~\Longleftrightarrow~c=ab^{-1}	371	386	\textrm{c=a.xi.b}~\Longleftrightarrow~c=ab^{-1}
\]	372	387	\]
	373	388
\paragraph{.tx.}	374	389	\paragraph{.tx.}
This operator gives the matrix product of the transpose matrix of the first operand and the second one.	375	390	This operator gives the matrix product of the transpose matrix of the first operand and the second one.
\[	376	391	\[
\textrm{c=a.tx.b}~\Longleftrightarrow~c={}^t\!ab	377	392	\textrm{c=a.tx.b}~\Longleftrightarrow~c={}^t\!ab
\]	378	393	\]
	379	394
	380	395
\paragraph{.xt.}	381	396	\paragraph{.xt.}
This operator gives the matrix product of the first operand and the transpose matrix of the second one.	382	397	This operator gives the matrix product of the first operand and the transpose matrix of the second one.
\[	383	398	\[
\textrm{c=a.xt.b}~\Longleftrightarrow~c=a{}^t\!b	384	399	\textrm{c=a.xt.b}~\Longleftrightarrow~c=a{}^t\!b
\]	385	400	\]
	386	401
\paragraph{.hx.}	387	402	\paragraph{.hx.}
This operator gives the matrix product of the conjugate transpose matrix of the first operand and the second one.	388	403	This operator gives the matrix product of the conjugate transpose matrix of the first operand and the second one.
\[	389	404	\[
\textrm{c=a.hx.b}~\Longleftrightarrow~c=a^*b	390	405	\textrm{c=a.hx.b}~\Longleftrightarrow~c=a^*b
\]	391	406	\]
	392	407
\paragraph{.xh.}	393	408	\paragraph{.xh.}
This operator gives the matrix product of thee first operand and the conjugate transpose matrix of the second one.	394	409	This operator gives the matrix product of thee first operand and the conjugate transpose matrix of the second one.
\[	395	410	\[
\textrm{c=a.xh.b}~\Longleftrightarrow~c=ab^*	396	411	\textrm{c=a.xh.b}~\Longleftrightarrow~c=ab^*
\]	397	412	\]
	398	413
	399	414
\section{Interpolation}	400	415	\section{Interpolation}
	401	416
\subsection{Quadratic Interpolation}	402	417	\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.	403	418	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}	404	419	\subsubsection{One variable function}
\begin{verbatim}	405	420	\begin{verbatim}
		421	Module : use fvn_interpol
value=fvn_quad_interpol(x,n,xdata,ydata)	406	422	value=fvn_quad_interpol(x,n,xdata,ydata)
\end{verbatim}	407	423	\end{verbatim}
\begin{itemize}	408	424	\begin{itemize}
\item x is the real where we want to evaluate the function	409	425	\item x is the real where we want to evaluate the function
\item n is the number of tabulated values	410	426	\item n is the number of tabulated values
\item xdata(n) contains the tabulated coordinates	411	427	\item xdata(n) contains the tabulated coordinates
\item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))	412	428	\item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))
\end{itemize}	413	429	\end{itemize}
xdata must be strictly increasingly ordered.	414	430	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	415	431	x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation
\paragraph*{Example}	416	432	\paragraph*{Example}
\begin{verbatim}	417	433	\begin{verbatim}
program inter1d	418	434	program inter1d
	419	435
use fvn	420	436	use fvn_interpol
implicit none	421	437	implicit none
	422	438
integer(kind=4),parameter :: ndata=33	423	439	integer(kind=4),parameter :: ndata=33
integer(kind=4) :: i,nout	424	440	integer(kind=4) :: i,nout
real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)	425	441	real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
real(kind=8) ::tv	426	442	real(kind=8) ::tv
	427	443
intrinsic sin	428	444	intrinsic sin
	429	445
f(x)=sin(x)	430	446	f(x)=sin(x)
	431	447
xdata(1)=0.	432	448	xdata(1)=0.
fdata(1)=f(xdata(1))	433	449	fdata(1)=f(xdata(1))
h=1./32.	434	450	h=1./32.
do i=2,ndata	435	451	do i=2,ndata
xdata(i)=xdata(i-1)+h	436	452	xdata(i)=xdata(i-1)+h
fdata(i)=f(xdata(i))	437	453	fdata(i)=f(xdata(i))
end do	438	454	end do
call random_seed()	439	455	call random_seed()
call random_number(x)	440	456	call random_number(x)
	441	457
q=fvn_d_quad_interpol(x,ndata,xdata,fdata)	442	458	q=fvn_d_quad_interpol(x,ndata,xdata,fdata)
	443	459
tv=f(x)	444	460	tv=f(x)
write(,) "x ",x	445	461	write(,) "x ",x
write(,) "Calculated (real) value :",tv	446	462	write(,) "Calculated (real) value :",tv
write(,) "fvn interpolation :",q	447	463	write(,) "fvn interpolation :",q
write(,) "Relative fvn error :",abs((q-tv)/tv)	448	464	write(,) "Relative fvn error :",abs((q-tv)/tv)
	449	465
end program	450	466	end program
	451	467
\end{verbatim}	452	468	\end{verbatim}
	453	469
	454	470
\subsubsection{Two variables function}	455	471	\subsubsection{Two variables function}
\begin{verbatim}	456	472	\begin{verbatim}
		473	Module : use fvn_interpol
value=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)	457	474	value=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
\end{verbatim}	458	475	\end{verbatim}
\begin{itemize}	459	476	\begin{itemize}
\item x,y are the real coordinates where we want to evaluate the function	460	477	\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	461	478	\item nx is the number of tabulated values along x axis
\item xdata(nx) contains the tabulated x	462	479	\item xdata(nx) contains the tabulated x
\item ny is the number of tabulated values along y axis	463	480	\item ny is the number of tabulated values along y axis
\item ydata(ny) contains the tabulated y	464	481	\item ydata(ny) contains the tabulated y
\item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))	465	482	\item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))
\end{itemize}	466	483	\end{itemize}
xdata and ydata must be strictly increasingly ordered.	467	484	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	468	485	(x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
	469	486
\paragraph*{Example}	470	487	\paragraph*{Example}
	471	488
\begin{verbatim}	472	489	\begin{verbatim}
program inter2d	473	490	program inter2d
use fvn	474	491	use fvn_interpol
implicit none	475	492	implicit none
	476	493
integer(kind=4),parameter :: nx=21,ny=42	477	494	integer(kind=4),parameter :: nx=21,ny=42
integer(kind=4) :: i,j	478	495	integer(kind=4) :: i,j
real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)	479	496	real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
real(kind=8) :: tv	480	497	real(kind=8) :: tv
	481	498
intrinsic dble,sin	482	499	intrinsic dble,sin
	483	500
f(x,y)=sin(x+2.*y)	484	501	f(x,y)=sin(x+2.*y)
do i=1,nx	485	502	do i=1,nx
xdata(i)=dble(i-1)/dble(nx-1)	486	503	xdata(i)=dble(i-1)/dble(nx-1)
end do	487	504	end do
do i=1,ny	488	505	do i=1,ny
ydata(i)=dble(i-1)/dble(ny-1)	489	506	ydata(i)=dble(i-1)/dble(ny-1)
end do	490	507	end do
do i=1,nx	491	508	do i=1,nx
do j=1,ny	492	509	do j=1,ny
fdata(i,j)=f(xdata(i),ydata(j))	493	510	fdata(i,j)=f(xdata(i),ydata(j))
end do	494	511	end do
end do	495	512	end do
call random_seed()	496	513	call random_seed()
call random_number(x)	497	514	call random_number(x)
call random_number(y)	498	515	call random_number(y)
	499	516
q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)	500	517	q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
tv=f(x,y)	501	518	tv=f(x,y)
	502	519
write(,) "x y",x,y	503	520	write(,) "x y",x,y
write(,) "Calculated (real) value :",tv	504	521	write(,) "Calculated (real) value :",tv
write(,) "fvn interpolation :",q	505	522	write(,) "fvn interpolation :",q
write(,) "Relative fvn error :",abs((q-tv)/tv)	506	523	write(,) "Relative fvn error :",abs((q-tv)/tv)
	507	524
end program	508	525	end program
	509	526
\end{verbatim}	510	527	\end{verbatim}
	511	528
	512	529
	513	530
\subsubsection{Three variables function}	514	531	\subsubsection{Three variables function}
\begin{verbatim}	515	532	\begin{verbatim}
		533	Module : use fvn_interpol
value=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)	516	534	value=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
\end{verbatim}	517	535	\end{verbatim}
\begin{itemize}	518	536	\begin{itemize}
\item x,y,z are the real coordinates where we want to evaluate the function	519	537	\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	520	538	\item nx is the number of tabulated values along x axis
\item xdata(nx) contains the tabulated x	521	539	\item xdata(nx) contains the tabulated x
\item ny is the number of tabulated values along y axis	522	540	\item ny is the number of tabulated values along y axis
\item ydata(ny) contains the tabulated y	523	541	\item ydata(ny) contains the tabulated y
\item nz is the number of tabulated values along z axis	524	542	\item nz is the number of tabulated values along z axis
\item zdata(ny) contains the tabulated z	525	543	\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))	526	544	\item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k))
\end{itemize}	527	545	\end{itemize}
xdata, ydata and zdata must be strictly increasingly ordered.	528	546	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	529	547	(x,y,z) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
	530	548
\paragraph*{Example}	531	549	\paragraph*{Example}
\begin{verbatim}	532	550	\begin{verbatim}
program inter3d	533	551	program inter3d
use fvn	534	552	use fvn_interpol
	535	553
implicit none	536	554	implicit none
	537	555
integer(kind=4),parameter :: nx=21,ny=42,nz=18	538	556	integer(kind=4),parameter :: nx=21,ny=42,nz=18
integer(kind=4) :: i,j,k	539	557	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)	540	558	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	541	559	real(kind=8) :: tv
	542	560
intrinsic dble,sin	543	561	intrinsic dble,sin
	544	562
f(x,y,z)=sin(x+2.y+3.z)	545	563	f(x,y,z)=sin(x+2.y+3.z)
do i=1,nx	546	564	do i=1,nx
xdata(i)=2.*(dble(i-1)/dble(nx-1))	547	565	xdata(i)=2.*(dble(i-1)/dble(nx-1))
end do	548	566	end do
do i=1,ny	549	567	do i=1,ny
ydata(i)=2.*(dble(i-1)/dble(ny-1))	550	568	ydata(i)=2.*(dble(i-1)/dble(ny-1))
end do	551	569	end do
do i=1,nz	552	570	do i=1,nz
zdata(i)=2.*(dble(i-1)/dble(nz-1))	553	571	zdata(i)=2.*(dble(i-1)/dble(nz-1))
end do	554	572	end do
do i=1,nx	555	573	do i=1,nx
do j=1,ny	556	574	do j=1,ny
do k=1,nz	557	575	do k=1,nz
fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))	558	576	fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
end do	559	577	end do
end do	560	578	end do
end do	561	579	end do
call random_seed()	562	580	call random_seed()
call random_number(x)	563	581	call random_number(x)
call random_number(y)	564	582	call random_number(y)
call random_number(z)	565	583	call random_number(z)
	566	584
q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)	567	585	q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
tv=f(x,y,z)	568	586	tv=f(x,y,z)
	569	587
write(,) "x y z",x,y,z	570	588	write(,) "x y z",x,y,z
write(,) "Calculated (real) value :",tv	571	589	write(,) "Calculated (real) value :",tv
write(,) "fvn interpolation :",q	572	590	write(,) "fvn interpolation :",q
write(,) "Relative fvn error :",abs((q-tv)/tv)	573	591	write(,) "Relative fvn error :",abs((q-tv)/tv)
	574	592
end program	575	593	end program
	576	594
\end{verbatim}	577	595	\end{verbatim}
	578	596
\subsubsection{Utility procedure}	579	597	\subsubsection{Utility procedure}
fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.	580	598	fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.
\begin{verbatim}	581	599	\begin{verbatim}
		600	Module : use fvn_interpol
call fvn_find_interval(x,i,xdata,n)	582	601	call fvn_find_interval(x,i,xdata,n)
\end{verbatim}	583	602	\end{verbatim}
\begin{itemize}	584	603	\begin{itemize}
\item x (in) the real value to locate	585	604	\item x (in) the real value to locate
\item i (out) the resulting indice	586	605	\item i (out) the resulting indice
\item xdata(n) (in) increasingly ordered array	587	606	\item xdata(n) (in) increasingly ordered array
\item n (in) size of the array	588	607	\item n (in) size of the array
\end{itemize}	589	608	\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.	590	609	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.
	591	610
	592	611
	593	612
\subsection{Akima spline}	594	613	\subsection{Akima spline}
fvn provides Akima spline interpolation and evaluation for both single and double precision real.	595	614	fvn provides Akima spline interpolation and evaluation for both single and double precision real.
\subsubsection{Interpolation}	596	615	\subsubsection{Interpolation}
\begin{verbatim}	597	616	\begin{verbatim}
		617	Module : use fvn_interpol
call fvn_akima(n,x,y,br,co)	598	618	call fvn_akima(n,x,y,br,co)
\end{verbatim}	599	619	\end{verbatim}
\begin{itemize}	600	620	\begin{itemize}
\item n (in) is an integer equal to the number of points	601	621	\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	602	622	\item x(n) (in) ,y(n) (in) are the known couples of coordinates
\item br (out) on output contains a copy of x	603	623	\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.	604	624	\item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval.
\end{itemize}	605	625	\end{itemize}
	606	626
\subsubsection{Evaluation}	607	627	\subsubsection{Evaluation}
\begin{verbatim}	608	628	\begin{verbatim}
		629	Module : use fvn_interpol
y=fvn_spline_eval(x,n,br,co)	609	630	y=fvn_spline_eval(x,n,br,co)
\end{verbatim}	610	631	\end{verbatim}
\begin{itemize}	611	632	\begin{itemize}
\item x (in) is the point where we want to evaluate	612	633	\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) \\	613	634	\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)	614	635	are the outputs of fvn\_x\_akima(n,x,y,br,co)
\end{itemize}	615	636	\end{itemize}
	616	637
\subsubsection{Example}	617	638	\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).	618	639	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).
	619	640
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.	620	641	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.
	621	642
\begin{verbatim}	622	643	\begin{verbatim}
program akima	623	644	program akima
use fvn	624	645	use fvn_interpol
implicit none	625	646	implicit none
	626	647
integer :: nbpoints,nppoints,i	627	648	integer :: nbpoints,nppoints,i
real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d	628	649	real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
real(8),dimension(:,:),allocatable :: coeff_fvn_d	629	650	real(8),dimension(:,:),allocatable :: coeff_fvn_d
real(8) :: xstep_d,xp_d,ty_d,fvn_y_d	630	651	real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
	631	652
open(2,file='fvn_akima_double.dat')	632	653	open(2,file='fvn_akima_double.dat')
open(3,file='fvn_akima_breakpoints_double.dat')	633	654	open(3,file='fvn_akima_breakpoints_double.dat')
nbpoints=30	634	655	nbpoints=30
allocate(x_d(nbpoints))	635	656	allocate(x_d(nbpoints))
allocate(y_d(nbpoints))	636	657	allocate(y_d(nbpoints))
allocate(breakpoints_d(nbpoints))	637	658	allocate(breakpoints_d(nbpoints))
allocate(coeff_fvn_d(4,nbpoints))	638	659	allocate(coeff_fvn_d(4,nbpoints))
	639	660
xstep_d=20./dfloat(nbpoints)	640	661	xstep_d=20./dfloat(nbpoints)
do i=1,nbpoints	641	662	do i=1,nbpoints
x_d(i)=-10.+dfloat(i)*xstep_d	642	663	x_d(i)=-10.+dfloat(i)*xstep_d
y_d(i)=dsin(x_d(i))	643	664	y_d(i)=dsin(x_d(i))
write(3,44) (x_d(i),y_d(i))	644	665	write(3,44) (x_d(i),y_d(i))
end do	645	666	end do
close(3)	646	667	close(3)
	647	668
call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)	648	669	call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
	649	670
nppoints=1000	650	671	nppoints=1000
xstep_d=22./dfloat(nppoints)	651	672	xstep_d=22./dfloat(nppoints)
do i=1,nppoints	652	673	do i=1,nppoints
xp_d=-11.+dfloat(i)*xstep_d	653	674	xp_d=-11.+dfloat(i)*xstep_d
ty_d=dsin(xp_d)	654	675	ty_d=dsin(xp_d)
fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)	655	676	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)	656	677	write(2,44) (xp_d,ty_d,fvn_y_d)
end do	657	678	end do
	658	679
close(2)	659	680	close(2)
	660	681
44 FORMAT(4(1X,1PE22.14))	661	682	44 FORMAT(4(1X,1PE22.14))
	662	683
end program	663	684	end program
	664	685
\end{verbatim}	665	686	\end{verbatim}
Results are plotted on figure \ref{akima}	666	687	Results are plotted on figure \ref{akima}
	667	688
\begin{figure}	668	689	\begin{figure}
\begin{center}	669	690	\begin{center}
\includegraphics[width=0.9\textwidth]{akima.pdf}	670	691	\includegraphics[width=0.9\textwidth]{akima.pdf}
% akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720	671	692	% akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720
\caption{Akima Spline Interpolation}	672	693	\caption{Akima Spline Interpolation}
\label{akima}	673	694	\label{akima}
\end{center}	674	695	\end{center}
	675	696
\end{figure}	676	697	\end{figure}
	677	698
	678	699
	679	700
\section{Least square polynomial}	680	701	\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 sgels (dgels), which solve this problem.	681	702	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 sgels (dgels), which solve this problem.
	682	703
\begin{verbatim}	683	704	\begin{verbatim}
		705	Module : use fvn_linear
call fvn_lspoly(np,x,y,deg,coeff,status)	684	706	call fvn_lspoly(np,x,y,deg,coeff,status)
\end{verbatim}	685	707	\end{verbatim}
\begin{itemize}	686	708	\begin{itemize}
\item np (in) is an integer equal to the number of points	687	709	\item np (in) is an integer equal to the number of points
\item x(np) (in),y(np) (in) are the known coordinates	688	710	\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.	689	711	\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	690	712	\item coeff(deg+1) (out) on output contains the polynomial coefficients
\item status (out) is an integer containing 0 if a problem occured.	691	713	\item status (out) is an integer containing 0 if a problem occured.
\end{itemize}	692	714	\end{itemize}
	693	715
\subsection*{Example}	694	716	\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 :	695	717	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}	696	718	\begin{verbatim}
program lsp	697	719	program lsp
use fvn	698	720	use fvn_linear
implicit none	699	721	implicit none
	700	722
integer,parameter :: npoints=13,deg=3	701	723	integer,parameter :: npoints=13,deg=3
integer :: status,i	702	724	integer :: status,i
real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc	703	725	real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
real(kind=8) :: coeff(deg+1)	704	726	real(kind=8) :: coeff(deg+1)
	705	727
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. /)	706	728	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 /)	707	729	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 /)
	708	730
open(2,file='fvn_lsp_double_mesure.dat')	709	731	open(2,file='fvn_lsp_double_mesure.dat')
open(3,file='fvn_lsp_double_poly.dat')	710	732	open(3,file='fvn_lsp_double_poly.dat')
	711	733
do i=1,npoints	712	734	do i=1,npoints
write(2,44) xm(i),ym(i)	713	735	write(2,44) xm(i),ym(i)
end do	714	736	end do
close(2)	715	737	close(2)
	716	738
	717	739
call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)	718	740	call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)
	719	741
xstep=(xm(npoints)-xm(1))/1000.	720	742	xstep=(xm(npoints)-xm(1))/1000.
do i=1,1000	721	743	do i=1,1000
xc=xm(1)+(i-1)*xstep	722	744	xc=xm(1)+(i-1)*xstep
yc=poly(xc,coeff)	723	745	yc=poly(xc,coeff)
write(3,44) xc,yc	724	746	write(3,44) xc,yc
end do	725	747	end do
close(3)	726	748	close(3)
	727	749
44 FORMAT(4(1X,1PE22.14))	728	750	44 FORMAT(4(1X,1PE22.14))
	729	751
contains	730	752	contains
function poly(x,coeff)	731	753	function poly(x,coeff)
implicit none	732	754	implicit none
real(8) :: x	733	755	real(8) :: x
real(8) :: coeff(deg+1)	734	756	real(8) :: coeff(deg+1)
real(8) :: poly	735	757	real(8) :: poly
integer :: i	736	758	integer :: i
	737	759
poly=0.	738	760	poly=0.
	739	761
do i=1,deg+1	740	762	do i=1,deg+1
poly=poly+coeff(i)x*(i-1)	741	763	poly=poly+coeff(i)x*(i-1)
end do	742	764	end do
	743	765
end function	744	766	end function
end program	745	767	end program
\end{verbatim}	746	768	\end{verbatim}
The results are plotted on figure \ref{lsp} .	747	769	The results are plotted on figure \ref{lsp} .
	748	770
\begin{figure}	749	771	\begin{figure}
\begin{center}	750	772	\begin{center}
\includegraphics[width=0.9\textwidth]{lsp.pdf}	751	773	\includegraphics[width=0.9\textwidth]{lsp.pdf}
\caption{Least Square Polynomial}	752	774	\caption{Least Square Polynomial}
\label{lsp}	753	775	\label{lsp}
\end{center}	754	776	\end{center}
\end{figure}	755	777	\end{figure}
	756	778
	757	779
	758	780
\section{Zero finding}	759	781	\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}.	760	782	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}.
	761	783
\begin{verbatim}	762	784	\begin{verbatim}
		785	Module : use fvn_misc
call fvn_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)	763	786	call fvn_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
\end{verbatim}	764	787	\end{verbatim}
\begin{itemize}	765	788	\begin{itemize}
\item f (in) is the complex function (kind=8) for which we search zeros	766	789	\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.	767	790	\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.	768	791	\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.	769	792	\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.	770	793	\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.	771	794	\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.	772	795	\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.	773	796	\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	774	797	\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.	775	798	\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}	776	799	\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.	777	800	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.
	778	801
\subsection*{Example}	779	802	\subsection*{Example}
Example to find the ten roots of $x^{10}-1$	780	803	Example to find the ten roots of $x^{10}-1$
\begin{verbatim}	781	804	\begin{verbatim}
program muller	782	805	program muller
use fvn	783	806	use fvn_misc
implicit none	784	807	implicit none
	785	808
integer :: i,info	786	809	integer :: i,info
complex(8),dimension(10) :: roots	787	810	complex(8),dimension(10) :: roots
integer,dimension(10) :: infer	788	811	integer,dimension(10) :: infer
complex(8), external :: f	789	812	complex(8), external :: f
	790	813
call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)	791	814	call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
	792	815
write(,) "Error code :",info	793	816	write(,) "Error code :",info
do i=1,10	794	817	do i=1,10
write(,) roots(i),infer(i)	795	818	write(,) roots(i),infer(i)
enddo	796	819	enddo
end program	797	820	end program
	798	821
function f(x)	799	822	function f(x)
complex(8) :: x,f	800	823	complex(8) :: x,f
f=x**10-1	801	824	f=x**10-1
end function	802	825	end function
	803	826
\end{verbatim}	804	827	\end{verbatim}
	805	828
	806	829
	807	830
\section{Numerical integration}	808	831	\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.	809	832	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.
	810	833
\subsection{Gauss Legendre Abscissas and Weigth}	811	834	\subsection{Gauss Legendre Abscissas and Weigth}
This subroutine was inspired by Numerical Recipes routine gauleg.	812	835	This subroutine was inspired by Numerical Recipes routine gauleg.
\begin{verbatim}	813	836	\begin{verbatim}
		837	Module : use fvn_integ
call fvn_gauss_legendre(n,qx,qw)	814	838	call fvn_gauss_legendre(n,qx,qw)
\end{verbatim}	815	839	\end{verbatim}
\begin{itemize}	816	840	\begin{itemize}
\item n (in) is an integer equal to the number of Gauss Legendre points	817	841	\item n (in) is an integer equal to the number of Gauss Legendre points
\item qx (out) is a real(8) vector of length n containing the abscissas.	818	842	\item qx (out) is a real(8) vector of length n containing the abscissas.
\item qw (out) is a real(8) vector of length n containing the weigths.	819	843	\item qw (out) is a real(8) vector of length n containing the weigths.
\end{itemize}	820	844	\end{itemize}
This subroutine computes n Gauss-Legendre abscissas and weigths	821	845	This subroutine computes n Gauss-Legendre abscissas and weigths
	822	846
\subsection{Gauss Legendre Numerical Integration}	823	847	\subsection{Gauss Legendre Numerical Integration}
\begin{verbatim}	824	848	\begin{verbatim}
		849	Module : use fvn_integ
call fvn_gl_integ(f,a,b,n,res)	825	850	call fvn_gl_integ(f,a,b,n,res)
\end{verbatim}	826	851	\end{verbatim}
\begin{itemize}	827	852	\begin{itemize}
\item f (in) is a real(8) function to integrate	828	853	\item f (in) is a real(8) function to integrate
\item a (in) and b (in) are real(8) respectively lower and higher bound of integration	829	854	\item a (in) and b (in) are real(8) respectively lower and higher bound of integration
\item n (in) is an integer equal to the number of Gauss Legendre points to use	830	855	\item n (in) is an integer equal to the number of Gauss Legendre points to use
\item res (out) is a real(8) containing the result	831	856	\item res (out) is a real(8) containing the result
\end{itemize}	832	857	\end{itemize}
This function is a simple Gauss Legendre integration subroutine, which evaluate the integral of function f as in equation \ref{intsple} using n Gauss-Legendre pairs.	833	858	This function is a simple Gauss Legendre integration subroutine, which evaluate the integral of function f as in equation \ref{intsple} using n Gauss-Legendre pairs.
	834	859
\subsection{Gauss Kronrod Adaptative Integration}	835	860	\subsection{Gauss Kronrod Adaptative Integration}
This kind of numerical integration is an iterative procedure which try to achieve a given precision.	836	861	This kind of numerical integration is an iterative procedure which try to achieve a given precision.
\subsubsection{Numerical integration of a one variable function}	837	862	\subsubsection{Numerical integration of a one variable function}
\begin{verbatim}	838	863	\begin{verbatim}
		864	Module : use fvn_integ
call fvn_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)	839	865	call fvn_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
\end{verbatim}	840	866	\end{verbatim}
This routine evaluate the integral of function f as in equation \ref{intsple}	841	867	This routine evaluate the integral of function f as in equation \ref{intsple}
\begin{itemize}	842	868	\begin{itemize}
\item f (in) is an external real(8) function of one variable	843	869	\item f (in) is an external real(8) function of one variable
\item a (in) and b (in) are real(8) respectively lower an higher bound of integration	844	870	\item a (in) and b (in) are real(8) respectively lower an higher bound of integration
\item epsabs (in) and epsrel (in) are real(8) respectively desired absolute and relative error	845	871	\item epsabs (in) and epsrel (in) are real(8) respectively desired absolute and relative error
\item key (in) is an integer between 1 and 6 correspondind to the Gauss-Kronrod rule to use :	846	872	\item key (in) is an integer between 1 and 6 correspondind to the Gauss-Kronrod rule to use :
\begin{itemize}	847	873	\begin{itemize}
\item 1 : 7 - 15 points	848	874	\item 1 : 7 - 15 points
\item 2 : 10 - 21 points	849	875	\item 2 : 10 - 21 points
\item 3 : 15 - 31 points	850	876	\item 3 : 15 - 31 points
\item 4 : 20 - 41 points	851	877	\item 4 : 20 - 41 points
\item 5 : 25 - 51 points	852	878	\item 5 : 25 - 51 points
\item 6 : 30 - 61 points	853	879	\item 6 : 30 - 61 points
\end{itemize}	854	880	\end{itemize}
\item res (out) is a real(8) containing the estimation of the integration.	855	881	\item res (out) is a real(8) containing the estimation of the integration.
\item abserr (out) is a real(8) equal to the estimated absolute error	856	882	\item abserr (out) is a real(8) equal to the estimated absolute error
\item ier (out) is an integer used as an error flag	857	883	\item ier (out) is an integer used as an error flag
\begin{itemize}	858	884	\begin{itemize}
\item 0 : no error	859	885	\item 0 : no error
\item 1 : maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yield no improvement it is advised to analyze the integrand in order to determine the integration difficulaties. If the position of a local difficulty can be determined (i.e.singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved.	860	886	\item 1 : maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yield no improvement it is advised to analyze the integrand in order to determine the integration difficulaties. If the position of a local difficulty can be determined (i.e.singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved.
\item 2 : the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved.	861	887	\item 2 : the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved.
\item 3 : extremely bad integrand behaviour occurs at some points of the integration interval.	862	888	\item 3 : extremely bad integrand behaviour occurs at some points of the integration interval.
\item 6 : the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit4. result, abserr, neval, last are set to zero. Except when lenw is invalid, iwork(1), work(limit2+1) and work(limit3+1) are set to zero, work(1) is set to a and work(limit+1) to b.	863	889	\item 6 : the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit4. result, abserr, neval, last are set to zero. Except when lenw is invalid, iwork(1), work(limit2+1) and work(limit3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
\end{itemize}	864	890	\end{itemize}
\item limit (in) is an optional integer equal to maximum number of subintervals in the partition of the given integration interval (a,b). If the parameter is not present a default value of 500 will be used.	865	891	\item limit (in) is an optional integer equal to maximum number of subintervals in the partition of the given integration interval (a,b). If the parameter is not present a default value of 500 will be used.
\end{itemize}	866	892	\end{itemize}
	867	893
\begin{equation}	868	894	\begin{equation}
\int_a^b f(x)~dx	869	895	\int_a^b f(x)~dx
\label{intsple}	870	896	\label{intsple}
\end{equation}	871	897	\end{equation}
	872	898
	873	899
	874	900
	875	901
\subsubsection{Numerical integration of a two variable function}	876	902	\subsubsection{Numerical integration of a two variable function}
\begin{verbatim}	877	903	\begin{verbatim}
		904	Module : use fvn_integ
call fvn_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)	878	905	call fvn_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
\end{verbatim}	879	906	\end{verbatim}
This function evaluate the integral of a function f(x,y) as defined in equation \ref{intdble}. The parameters of same name as in the previous paragraph have exactly the same function and behaviour thus only what differs is decribed here	880	907	This function evaluate the integral of a function f(x,y) as defined in equation \ref{intdble}. The parameters of same name as in the previous paragraph have exactly the same function and behaviour thus only what differs is decribed here
\begin{itemize}	881	908	\begin{itemize}
\item a (in) and b (in) are real(8) corresponding respectively to lower and higher bound of integration for the x variable.	882	909	\item a (in) and b (in) are real(8) corresponding respectively to lower and higher bound of integration for the x variable.
\item g(x) (in) and h(x) (in) are external functions describing the lower and higher bound of integration for the y variable as a function of x.	883	910	\item g(x) (in) and h(x) (in) are external functions describing the lower and higher bound of integration for the y variable as a function of x.
\end{itemize}	884	911	\end{itemize}
	885	912
\begin{equation}	886	913	\begin{equation}
\int_a^b \int_{g(x)}^{h(x)} f(x,y)~dy~dx	887	914	\int_a^b \int_{g(x)}^{h(x)} f(x,y)~dy~dx
\label{intdble}	888	915	\label{intdble}
\end{equation}	889	916	\end{equation}
	890	917
\subsubsection*{Example}	891	918	\subsubsection*{Example}
\begin{verbatim}	892	919	\begin{verbatim}
program integ	893	920	program integ
use fvn	894	921	use fvn_integ
implicit none	895	922	implicit none
	896	923
real(8), external :: f1,f2,g,h	897	924	real(8), external :: f1,f2,g,h
real(8) :: a,b,epsabs,epsrel,abserr,res	898	925	real(8) :: a,b,epsabs,epsrel,abserr,res
integer :: key,ier	899	926	integer :: key,ier
	900	927
a=0.	901	928	a=0.
b=1.	902	929	b=1.
epsabs=1d-8	903	930	epsabs=1d-8
epsrel=1d-8	904	931	epsrel=1d-8
key=2	905	932	key=2
call fvn_d_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier,500)	906	933	call fvn_d_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier,500)
write(,) "Integration of x*x between 0 and 1 : "	907	934	write(,) "Integration of x*x between 0 and 1 : "
write(,) res	908	935	write(,) res
	909	936
call fvn_d_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,500)	910	937	call fvn_d_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,500)
write(,) "Integration of x*y between 0 and 1 on both x and y : "	911	938	write(,) "Integration of x*y between 0 and 1 on both x and y : "
write(,) res	912	939	write(,) res
	913	940
	914	941
end program	915	942	end program
	916	943
function f1(x)	917	944	function f1(x)
implicit none	918	945	implicit none
real(8) :: x,f1	919	946	real(8) :: x,f1
f1=x*x	920	947	f1=x*x
end function	921	948	end function
	922	949
function f2(x,y)	923	950	function f2(x,y)
implicit none	924	951	implicit none
real(8) :: x,y,f2	925	952	real(8) :: x,y,f2
f2=x*y	926	953	f2=x*y
end function	927	954	end function
	928	955
function g(x)	929	956	function g(x)
implicit none	930	957	implicit none
real(8) :: x,g	931	958	real(8) :: x,g
g=0.	932	959	g=0.
end function	933	960	end function
	934	961
function h(x)	935	962	function h(x)
implicit none	936	963	implicit none
real(8) :: x,h	937	964	real(8) :: x,h
h=1.	938	965	h=1.
end function	939	966	end function
\end{verbatim}	940	967	\end{verbatim}
	941	968
	942	969
\section{Special functions}	943	970	\section{Special functions}
Specials functions are available in fvn by using an implementation of fnlib \url{http://www.netlib.org/fn} with some additions. This can be used separatly from the rest of fvn by using the module \verb'fvn_fnlib' and linking the library \verb'libfvn_fnlib.a' . The module provides a generic interfaces to all the routines. Specific names of the routines are given in the description.	944	971	Specials functions are available in fvn by using an implementation of fnlib \url{http://www.netlib.org/fn} with some additions. This can be used separatly from the rest of fvn by using the module \verb'fvn_fnlib' and linking the library \verb'libfvn_fnlib.a' . The module provides a generic interfaces to all the routines. Specific names of the routines are given in the description.
	945	972
		973	\begin{verbatim}
		974	Module : use fvn_fnlib
		975	\end{verbatim}
		976
\paragraph{Important Note}	946	977	\paragraph{Important Note}
Due to the addition of fnlib to fvn, some functions that were in fvn and are redondant will be removed from fvn, so update your code now and replace them with the fnlib version. These are listed here after :	947	978	Due to the addition of fnlib to fvn, some functions that were in fvn and are redondant are now removed from fvn, so update your code now and replace them with the fnlib version. These are listed here after :
\begin{itemize}	948	979	\begin{itemize}
\item \verb'fvn_z_acos' replaced by \verb'acos'	949	980	\item \verb'fvn_z_acos' replaced by \verb'acos'
\item \verb'fvn_z_asin' replaced by \verb'asin'	950	981	\item \verb'fvn_z_asin' replaced by \verb'asin'
\item \verb'fvn_d_asinh' replaced by \verb'asinh'	951	982	\item \verb'fvn_d_asinh' replaced by \verb'asinh'
\item \verb'fvn_d_acosh' replaced by \verb'acosh'	952	983	\item \verb'fvn_d_acosh' replaced by \verb'acosh'
\item \verb'fvn_s_csevl' replaced by \verb'csevl'	953	984	\item \verb'fvn_s_csevl' replaced by \verb'csevl'
\item \verb'fvn_d_csevl' replaced by \verb'csevl'	954	985	\item \verb'fvn_d_csevl' replaced by \verb'csevl'
\item \verb'fvn_d_factorial' replaced by \verb'fac'	955	986	\item \verb'fvn_d_factorial' replaced by \verb'fac'
\item \verb'fvn_d_lngamma' replaced by \verb'alngam'	956	987	\item \verb'fvn_d_lngamma' replaced by \verb'alngam'
\end{itemize}	957	988	\end{itemize}
	958	989
	959	990
\subsection{Elementary functions}	960	991	\subsection{Elementary functions}
\subsubsection{carg}	961	992	\subsubsection{carg}
\begin{verbatim}	962	993	\begin{verbatim}
carg(z)	963	994	carg(z)
\end{verbatim}	964	995	\end{verbatim}
\begin{itemize}	965	996	\begin{itemize}
\item z (in) is a complex	966	997	\item z (in) is a complex
\end{itemize}	967	998	\end{itemize}
This function evaluates the argument of the complex z. That is $\theta$ for $z=\rho e^{i\theta}$.	968	999	This function evaluates the argument of the complex z. That is $\theta$ for $z=\rho e^{i\theta}$.
	969	1000
Specific interfaces : \verb'carg,zarg'	970	1001	Specific interfaces : \verb'carg,zarg'
	971	1002
	972	1003
\subsubsection{cbrt}	973	1004	\subsubsection{cbrt}
\begin{verbatim}	974	1005	\begin{verbatim}
cbrt(x)	975	1006	cbrt(x)
\end{verbatim}	976	1007	\end{verbatim}
\begin{itemize}	977	1008	\begin{itemize}
\item x is a real or complex	978	1009	\item x is a real or complex
\end{itemize}	979	1010	\end{itemize}
This function evaluates the cubic root of the argument x.	980	1011	This function evaluates the cubic root of the argument x.
	981	1012
Specific interfaces : \verb'cbrt,dcbrt,ccbrt,zcbrt'	982	1013	Specific interfaces : \verb'cbrt,dcbrt,ccbrt,zcbrt'
	983	1014
	984	1015
\subsubsection{exprl}	985	1016	\subsubsection{exprl}
\begin{verbatim}	986	1017	\begin{verbatim}
exprl(x)	987	1018	exprl(x)
\end{verbatim}	988	1019	\end{verbatim}
\begin{itemize}	989	1020	\begin{itemize}
\item x is a real or complex	990	1021	\item x is a real or complex
\end{itemize}	991	1022	\end{itemize}
This function evaluates ${e^x-1}\over x$.	992	1023	This function evaluates ${e^x-1}\over x$.
	993	1024
Specific interfaces : \verb'exprel,dexprl,cexprl,zexprl'	994	1025	Specific interfaces : \verb'exprel,dexprl,cexprl,zexprl'
	995	1026
\subsubsection{log10}	996	1027	\subsubsection{log10}
\begin{verbatim}	997	1028	\begin{verbatim}
log10(x)	998	1029	log10(x)
\end{verbatim}	999	1030	\end{verbatim}
\begin{itemize}	1000	1031	\begin{itemize}
\item x is a real or complex	1001	1032	\item x is a real or complex
\end{itemize}	1002	1033	\end{itemize}
This function is an extension of the intrinsic function log10 to complex arguments.	1003	1034	This function is an extension of the intrinsic function log10 to complex arguments.
	1004	1035
Specific interfaces : \verb'clog10,zlog10'	1005	1036	Specific interfaces : \verb'clog10,zlog10'
	1006	1037
	1007	1038
\subsubsection{alnrel}	1008	1039	\subsubsection{alnrel}
\begin{verbatim}	1009	1040	\begin{verbatim}
alnrel(x)	1010	1041	alnrel(x)
\end{verbatim}	1011	1042	\end{verbatim}
\begin{itemize}	1012	1043	\begin{itemize}
\item x is a real or complex	1013	1044	\item x is a real or complex
\end{itemize}	1014	1045	\end{itemize}
This function evaluates $ln(1+x)$.	1015	1046	This function evaluates $ln(1+x)$.
	1016	1047
Specific interfaces : \verb'alnrel,dlnrel,clnrel,zlnrel'	1017	1048	Specific interfaces : \verb'alnrel,dlnrel,clnrel,zlnrel'
	1018	1049
	1019	1050
\subsection{Trigonometry}	1020	1051	\subsection{Trigonometry}
\subsubsection{tan}	1021	1052	\subsubsection{tan}
\begin{verbatim}	1022	1053	\begin{verbatim}
tan(x)	1023	1054	tan(x)
\end{verbatim}	1024	1055	\end{verbatim}
\begin{itemize}	1025	1056	\begin{itemize}
\item x is a real or complex	1026	1057	\item x is a real or complex
\end{itemize}	1027	1058	\end{itemize}
This function evaluates the tangent of the argument. It is an extension of the intrinsic function tan to complex arguments.	1028	1059	This function evaluates the tangent of the argument. It is an extension of the intrinsic function tan to complex arguments.
	1029	1060
Specific interfaces : \verb'ctan,ztan'	1030	1061	Specific interfaces : \verb'ctan,ztan'
	1031	1062
	1032	1063
\subsubsection{cot}	1033	1064	\subsubsection{cot}
\begin{verbatim}	1034	1065	\begin{verbatim}
cot(x)	1035	1066	cot(x)
\end{verbatim}	1036	1067	\end{verbatim}
\begin{itemize}	1037	1068	\begin{itemize}
\item x is a real or complex	1038	1069	\item x is a real or complex
\end{itemize}	1039	1070	\end{itemize}
This function evaluate the cotangent of the argument.	1040	1071	This function evaluate the cotangent of the argument.
	1041	1072
Specific interfaces : \verb'cot,dcot,ccot,zcot'	1042	1073	Specific interfaces : \verb'cot,dcot,ccot,zcot'
	1043	1074
\subsubsection{sindg}	1044	1075	\subsubsection{sindg}
\begin{verbatim}	1045	1076	\begin{verbatim}
sindg(x)	1046	1077	sindg(x)
\end{verbatim}	1047	1078	\end{verbatim}
\begin{itemize}	1048	1079	\begin{itemize}
\item x is a real	1049	1080	\item x is a real
\end{itemize}	1050	1081	\end{itemize}
This function evaluate the sinus of the argument expressed in degrees.	1051	1082	This function evaluate the sinus of the argument expressed in degrees.
	1052	1083
Specific interfaces : \verb'sindg,dsindg'	1053	1084	Specific interfaces : \verb'sindg,dsindg'
	1054	1085
	1055	1086
\subsubsection{cosdg}	1056	1087	\subsubsection{cosdg}
\begin{verbatim}	1057	1088	\begin{verbatim}
cosdg(x)	1058	1089	cosdg(x)
\end{verbatim}	1059	1090	\end{verbatim}
\begin{itemize}	1060	1091	\begin{itemize}
\item x is a real	1061	1092	\item x is a real
\end{itemize}	1062	1093	\end{itemize}
This function evaluate the cosinus of the argument expressed in degrees.	1063	1094	This function evaluate the cosinus of the argument expressed in degrees.
	1064	1095
Specific interfaces : \verb'cosdg,dcosdg'	1065	1096	Specific interfaces : \verb'cosdg,dcosdg'
	1066	1097
	1067	1098
\subsubsection{asin}	1068	1099	\subsubsection{asin}
\begin{verbatim}	1069	1100	\begin{verbatim}
asin(x)	1070	1101	asin(x)
\end{verbatim}	1071	1102	\end{verbatim}
\begin{itemize}	1072	1103	\begin{itemize}
\item x is a real or complex	1073	1104	\item x is a real or complex
\end{itemize}	1074	1105	\end{itemize}
This function evaluates the arc sine of the argument. It is an extension of the intrinsic function asin to complex arguments.	1075	1106	This function evaluates the arc sine of the argument. It is an extension of the intrinsic function asin to complex arguments.
	1076	1107
Specific interfaces : \verb'casin,zasin'	1077	1108	Specific interfaces : \verb'casin,zasin'
	1078	1109
\subsubsection{acos}	1079	1110	\subsubsection{acos}
\begin{verbatim}	1080	1111	\begin{verbatim}
acos(x)	1081	1112	acos(x)
\end{verbatim}	1082	1113	\end{verbatim}
\begin{itemize}	1083	1114	\begin{itemize}
\item x is a real or complex	1084	1115	\item x is a real or complex
\end{itemize}	1085	1116	\end{itemize}
This function evaluates the arc cosine of the argument. It is an extension of the intrinsic function acos to complex arguments.	1086	1117	This function evaluates the arc cosine of the argument. It is an extension of the intrinsic function acos to complex arguments.
	1087	1118
Specific interfaces : \verb'cacos,zacos'	1088	1119	Specific interfaces : \verb'cacos,zacos'
	1089	1120
	1090	1121
\subsubsection{atan}	1091	1122	\subsubsection{atan}
\begin{verbatim}	1092	1123	\begin{verbatim}
atan(x)	1093	1124	atan(x)
\end{verbatim}	1094	1125	\end{verbatim}
\begin{itemize}	1095	1126	\begin{itemize}
\item x is a real or complex	1096	1127	\item x is a real or complex
\end{itemize}	1097	1128	\end{itemize}
This function evaluates the arc tangent of the argument. It is an extension of the intrinsic function atan to complex arguments.	1098	1129	This function evaluates the arc tangent of the argument. It is an extension of the intrinsic function atan to complex arguments.
	1099	1130
Specific interfaces : \verb'catan,zatan'	1100	1131	Specific interfaces : \verb'catan,zatan'
	1101	1132
\subsubsection{atan2}	1102	1133	\subsubsection{atan2}
\begin{verbatim}	1103	1134	\begin{verbatim}
atan2(x,y)	1104	1135	atan2(x,y)
\end{verbatim}	1105	1136	\end{verbatim}
\begin{itemize}	1106	1137	\begin{itemize}
\item x,y are real or complex	1107	1138	\item x,y are real or complex
\end{itemize}	1108	1139	\end{itemize}
This function evaluates the arc tangent of $x \over y$. It is an extension of the intrinsic function atan2 to complex arguments.	1109	1140	This function evaluates the arc tangent of $x \over y$. It is an extension of the intrinsic function atan2 to complex arguments.
	1110	1141
Specific interfaces : \verb'catan2,zatan2'	1111	1142	Specific interfaces : \verb'catan2,zatan2'
	1112	1143
\subsubsection{sinh}	1113	1144	\subsubsection{sinh}
\begin{verbatim}	1114	1145	\begin{verbatim}
sinh(x)	1115	1146	sinh(x)
\end{verbatim}	1116	1147	\end{verbatim}
\begin{itemize}	1117	1148	\begin{itemize}
\item x is a real or complex	1118	1149	\item x is a real or complex
\end{itemize}	1119	1150	\end{itemize}
This function evaluates the hyperbolic sine of the argument. It is an extension of the intrinsic function sinh to complex arguments.	1120	1151	This function evaluates the hyperbolic sine of the argument. It is an extension of the intrinsic function sinh to complex arguments.
	1121	1152
Specific interfaces : \verb'csinh,zsinh'	1122	1153	Specific interfaces : \verb'csinh,zsinh'
	1123	1154
	1124	1155
\subsubsection{cosh}	1125	1156	\subsubsection{cosh}
\begin{verbatim}	1126	1157	\begin{verbatim}
cosh(x)	1127	1158	cosh(x)
\end{verbatim}	1128	1159	\end{verbatim}
\begin{itemize}	1129	1160	\begin{itemize}
\item x is a real or complex	1130	1161	\item x is a real or complex
\end{itemize}	1131	1162	\end{itemize}
This function evaluates the hyperbolic cosine of the argument. It is an extension of the intrinsic function cosh to complex arguments.	1132	1163	This function evaluates the hyperbolic cosine of the argument. It is an extension of the intrinsic function cosh to complex arguments.
	1133	1164
Specific interfaces : \verb'ccosh,zcosh'	1134	1165	Specific interfaces : \verb'ccosh,zcosh'
	1135	1166
\subsubsection{tanh}	1136	1167	\subsubsection{tanh}
\begin{verbatim}	1137	1168	\begin{verbatim}
tanh(x)	1138	1169	tanh(x)
\end{verbatim}	1139	1170	\end{verbatim}
This function evaluates the hyperbolic tangent of the argument. It is an extension of the intrinsic function tanh to complex arguments.	1140	1171	This function evaluates the hyperbolic tangent of the argument. It is an extension of the intrinsic function tanh to complex arguments.
	1141	1172
Specific interfaces : \verb'ctanh,ztanh'	1142	1173	Specific interfaces : \verb'ctanh,ztanh'
	1143	1174
\subsubsection{asinh}	1144	1175	\subsubsection{asinh}
\begin{verbatim}	1145	1176	\begin{verbatim}
asinh(x)	1146	1177	asinh(x)
\end{verbatim}	1147	1178	\end{verbatim}
\begin{itemize}	1148	1179	\begin{itemize}
\item x is a real or complex	1149	1180	\item x is a real or complex
\end{itemize}	1150	1181	\end{itemize}
This function evaluates the arc hyperbolic sine of the argument.	1151	1182	This function evaluates the arc hyperbolic sine of the argument.
	1152	1183
Specific interfaces : \verb'asinh,dasinh,casinh,zasinh'	1153	1184	Specific interfaces : \verb'asinh,dasinh,casinh,zasinh'
	1154	1185
\subsubsection{acosh}	1155	1186	\subsubsection{acosh}
\begin{verbatim}	1156	1187	\begin{verbatim}
acosh(x)	1157	1188	acosh(x)
\end{verbatim}	1158	1189	\end{verbatim}
\begin{itemize}	1159	1190	\begin{itemize}
\item x is a real or complex	1160	1191	\item x is a real or complex
\end{itemize}	1161	1192	\end{itemize}
This function evaluates the arc hyperbolic cosine of the argument.	1162	1193	This function evaluates the arc hyperbolic cosine of the argument.
	1163	1194
Specific interfaces : \verb'acosh,dacosh,cacosh,zacosh'	1164	1195	Specific interfaces : \verb'acosh,dacosh,cacosh,zacosh'
	1165	1196
\subsubsection{atanh}	1166	1197	\subsubsection{atanh}
\begin{verbatim}	1167	1198	\begin{verbatim}
atanh(x)	1168	1199	atanh(x)
\end{verbatim}	1169	1200	\end{verbatim}
\begin{itemize}	1170	1201	\begin{itemize}
\item x is a real or complex	1171	1202	\item x is a real or complex
\end{itemize}	1172	1203	\end{itemize}
This function evaluates the arc hyperbolic tangent of the argument.	1173	1204	This function evaluates the arc hyperbolic tangent of the argument.
	1174	1205
Specific interfaces : \verb'atanh,datanh,catanh,zatanh'	1175	1206	Specific interfaces : \verb'atanh,datanh,catanh,zatanh'
	1176	1207
\subsection{Exponential Integral and related}	1177	1208	\subsection{Exponential Integral and related}
\subsubsection{ei}	1178	1209	\subsubsection{ei}
\begin{verbatim}	1179	1210	\begin{verbatim}
ei(x)	1180	1211	ei(x)
\end{verbatim}	1181	1212	\end{verbatim}
\begin{itemize}	1182	1213	\begin{itemize}
\item x is a real	1183	1214	\item x is a real
\end{itemize}	1184	1215	\end{itemize}
This function evaluates the exponential integral for argument greater then 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{ei} for $x \neq 0$.	1185	1216	This function evaluates the exponential integral for argument greater then 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{ei} for $x \neq 0$.
\begin{equation}	1186	1217	\begin{equation}
\label{ei}	1187	1218	\label{ei}
ei(x)= - \int _{-x} ^\infty {e^{-t}\over t}dt	1188	1219	ei(x)= - \int _{-x} ^\infty {e^{-t}\over t}dt
\end{equation}	1189	1220	\end{equation}
	1190	1221
Specific interfaces : \verb'ei,dei'	1191	1222	Specific interfaces : \verb'ei,dei'
	1192	1223
	1193	1224
\subsubsection{e1}	1194	1225	\subsubsection{e1}
\begin{verbatim}	1195	1226	\begin{verbatim}
e1(x)	1196	1227	e1(x)
\end{verbatim}	1197	1228	\end{verbatim}
\begin{itemize}	1198	1229	\begin{itemize}
\item x is a real	1199	1230	\item x is a real
\end{itemize}	1200	1231	\end{itemize}
This function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.	1201	1232	This function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
\begin{equation}	1202	1233	\begin{equation}
\label{e1}	1203	1234	\label{e1}
e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt	1204	1235	e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt
\end{equation}	1205	1236	\end{equation}
	1206	1237
Specific interfaces : \verb'e1,de1'	1207	1238	Specific interfaces : \verb'e1,de1'
	1208	1239
\subsubsection{ali}	1209	1240	\subsubsection{ali}
\begin{verbatim}	1210	1241	\begin{verbatim}
ali(x)	1211	1242	ali(x)
\end{verbatim}	1212	1243	\end{verbatim}
\begin{itemize}	1213	1244	\begin{itemize}
\item x is a real	1214	1245	\item x is a real
\end{itemize}	1215	1246	\end{itemize}
This function evaluates the logarithm integral. it is define by equation \ref{ali} for $x > 0$ and $x \neq 1$.	1216	1247	This function evaluates the logarithm integral. it is define by equation \ref{ali} for $x > 0$ and $x \neq 1$.
\begin{equation}	1217	1248	\begin{equation}
\label{ali}	1218	1249	\label{ali}
ali(x)= - \int _0 ^x {dt \over ln(x)}	1219	1250	ali(x)= - \int _0 ^x {dt \over ln(x)}
\end{equation}	1220	1251	\end{equation}
	1221	1252
Specific interfaces : \verb'ali,dli'	1222	1253	Specific interfaces : \verb'ali,dli'
	1223	1254
\subsubsection{si}	1224	1255	\subsubsection{si}
\begin{verbatim}	1225	1256	\begin{verbatim}
si(x)	1226	1257	si(x)
\end{verbatim}	1227	1258	\end{verbatim}
\begin{itemize}	1228	1259	\begin{itemize}
\item x is a real	1229	1260	\item x is a real
\end{itemize}	1230	1261	\end{itemize}
This function evaluates the sine integral defined by equation \ref{si}.	1231	1262	This function evaluates the sine integral defined by equation \ref{si}.
\begin{equation}	1232	1263	\begin{equation}
\label{si}	1233	1264	\label{si}
si(x)= \int _0 ^x {sin(t) \over t }dt	1234	1265	si(x)= \int _0 ^x {sin(t) \over t }dt
\end{equation}	1235	1266	\end{equation}
	1236	1267
Specific interfaces : \verb'si,dsi'	1237	1268	Specific interfaces : \verb'si,dsi'
	1238	1269
	1239	1270
\subsubsection{ci}	1240	1271	\subsubsection{ci}
\begin{verbatim}	1241	1272	\begin{verbatim}
ci(x)	1242	1273	ci(x)
\end{verbatim}	1243	1274	\end{verbatim}
\begin{itemize}	1244	1275	\begin{itemize}
\item x is a real	1245	1276	\item x is a real
\end{itemize}	1246	1277	\end{itemize}
This function evaluates the cosine integral defined by equation \ref{ci} where $\gamma \approx 0.57721566$ represent Euler's constant.	1247	1278	This function evaluates the cosine integral defined by equation \ref{ci} where $\gamma \approx 0.57721566$ represent Euler's constant.
\begin{equation}	1248	1279	\begin{equation}
\label{ci}	1249	1280	\label{ci}
ci(x)= \gamma + ln(x) + \int _0 ^x {{1-cos(t)} \over t} dt	1250	1281	ci(x)= \gamma + ln(x) + \int _0 ^x {{1-cos(t)} \over t} dt
\end{equation}	1251	1282	\end{equation}
	1252	1283
Specific interfaces : \verb'ci,dci'	1253	1284	Specific interfaces : \verb'ci,dci'
	1254	1285
\subsubsection{cin}	1255	1286	\subsubsection{cin}
\begin{verbatim}	1256	1287	\begin{verbatim}
cin(x)	1257	1288	cin(x)
\end{verbatim}	1258	1289	\end{verbatim}
\begin{itemize}	1259	1290	\begin{itemize}
\item x is a real	1260	1291	\item x is a real
\end{itemize}	1261	1292	\end{itemize}
This function evaluates the cosine integral alternate definition given by equation \ref{cin}.	1262	1293	This function evaluates the cosine integral alternate definition given by equation \ref{cin}.
\begin{equation}	1263	1294	\begin{equation}
\label{cin}	1264	1295	\label{cin}
cin(x)= \int _0 ^x {{1-cos(t)} \over t} dt	1265	1296	cin(x)= \int _0 ^x {{1-cos(t)} \over t} dt
\end{equation}	1266	1297	\end{equation}
	1267	1298
Specific interface : \verb'cin,dcin'	1268	1299	Specific interface : \verb'cin,dcin'
	1269	1300
\subsubsection{shi}	1270	1301	\subsubsection{shi}
\begin{equation}	1271	1302	\begin{equation}
shi(x)	1272	1303	shi(x)
\end{equation}	1273	1304	\end{equation}
\begin{itemize}	1274	1305	\begin{itemize}
\item x is a real	1275	1306	\item x is a real
\end{itemize}	1276	1307	\end{itemize}
This function evaluates the hyperbolic sine integral defined by equation \ref{shi}.	1277	1308	This function evaluates the hyperbolic sine integral defined by equation \ref{shi}.
\begin{equation}	1278	1309	\begin{equation}
\label{shi}	1279	1310	\label{shi}
shi(x) = \int _0 ^x {sinh(t) \over t}dt	1280	1311	shi(x) = \int _0 ^x {sinh(t) \over t}dt
\end{equation}	1281	1312	\end{equation}
	1282	1313
Specific interfaces : \verb'shi,dshi'	1283	1314	Specific interfaces : \verb'shi,dshi'
	1284	1315
\subsubsection{chi}	1285	1316	\subsubsection{chi}
\begin{verbatim}	1286	1317	\begin{verbatim}
chi(x)	1287	1318	chi(x)
\end{verbatim}	1288	1319	\end{verbatim}
\begin{itemize}	1289	1320	\begin{itemize}
\item x is a real	1290	1321	\item x is a real
\end{itemize}	1291	1322	\end{itemize}
This function evaluates the hyperbolic cosine integral defined by equation \ref{chi} where $\gamma \approx 0.57721566$ represent Euler's constant.	1292	1323	This function evaluates the hyperbolic cosine integral defined by equation \ref{chi} where $\gamma \approx 0.57721566$ represent Euler's constant.
\begin{equation}	1293	1324	\begin{equation}
\label{chi}	1294	1325	\label{chi}
chi(x)= \gamma + ln(x) + \int _0 ^x {{cosh(t) -1} \over t}dt	1295	1326	chi(x)= \gamma + ln(x) + \int _0 ^x {{cosh(t) -1} \over t}dt
\end{equation}	1296	1327	\end{equation}
	1297	1328
Specific interfaces : chi,dchi	1298	1329	Specific interfaces : chi,dchi
	1299	1330
	1300	1331
\subsubsection{cinh}	1301	1332	\subsubsection{cinh}
\begin{verbatim}	1302	1333	\begin{verbatim}
cinh(x)	1303	1334	cinh(x)
\end{verbatim}	1304	1335	\end{verbatim}
\begin{itemize}	1305	1336	\begin{itemize}
\item x is a real	1306	1337	\item x is a real
\end{itemize}	1307	1338	\end{itemize}
This function evaluates the hyperbolic cosine integral alternate definition given by equation \ref{cinh}.	1308	1339	This function evaluates the hyperbolic cosine integral alternate definition given by equation \ref{cinh}.
\begin{equation}	1309	1340	\begin{equation}
\label{cinh}	1310	1341	\label{cinh}
cinh(x) = \int _0 ^x {{cosh(t) -1} \over t}dt	1311	1342	cinh(x) = \int _0 ^x {{cosh(t) -1} \over t}dt
\end{equation}	1312	1343	\end{equation}
	1313	1344
Specific interfaces : cinh,dcinh	1314	1345	Specific interfaces : cinh,dcinh
	1315	1346
	1316	1347
\subsection{Gamma function and related}	1317	1348	\subsection{Gamma function and related}
\subsubsection{fac}	1318	1349	\subsubsection{fac}
\begin{verbatim}	1319	1350	\begin{verbatim}
fac(n)	1320	1351	fac(n)
dfac(n)	1321	1352	dfac(n)
\end{verbatim}	1322	1353	\end{verbatim}
\begin{itemize}	1323	1354	\begin{itemize}
\item n is an integer	1324	1355	\item n is an integer
\end{itemize}	1325	1356	\end{itemize}
This function return $n!$ as a real(4) or real(8) for dfac. There's no generic interface for this one.	1326	1357	This function return $n!$ as a real(4) or real(8) for dfac. There's no generic interface for this one.
	1327	1358
Specific interfaces : \verb'fac,dfac'	1328	1359	Specific interfaces : \verb'fac,dfac'
	1329	1360
\subsubsection{binom}	1330	1361	\subsubsection{binom}
\begin{verbatim}	1331	1362	\begin{verbatim}
binom(n,m)	1332	1363	binom(n,m)
dbinom(n,m)	1333	1364	dbinom(n,m)
\end{verbatim}	1334	1365	\end{verbatim}
\begin{itemize}	1335	1366	\begin{itemize}
\item n,m are integers	1336	1367	\item n,m are integers
\end{itemize}	1337	1368	\end{itemize}
This function return the binomial coefficient defined by equation \ref{binom} with $n \geq m \geq 0$. binom returns a real(4), dbinom a real(8). There's no generic interface for this one.	1338	1369	This function return the binomial coefficient defined by equation \ref{binom} with $n \geq m \geq 0$. binom returns a real(4), dbinom a real(8). There's no generic interface for this one.
\begin{equation}	1339	1370	\begin{equation}
\label{binom}	1340	1371	\label{binom}
binom(n,m) = C_n^m = {{n!} \over {m!(n-m)!}}	1341	1372	binom(n,m) = C_n^m = {{n!} \over {m!(n-m)!}}
\end{equation}	1342	1373	\end{equation}
	1343	1374
Specific interfaces : \verb'binom,dbinom'	1344	1375	Specific interfaces : \verb'binom,dbinom'
	1345	1376
	1346	1377
\subsubsection{gamma}	1347	1378	\subsubsection{gamma}
\begin{verbatim}	1348	1379	\begin{verbatim}
gamma(x)	1349	1380	gamma(x)
\end{verbatim}	1350	1381	\end{verbatim}
\begin{itemize}	1351	1382	\begin{itemize}
\item x is a real or complex	1352	1383	\item x is a real or complex
\end{itemize}	1353	1384	\end{itemize}
This function evaluates $ \Gamma (x) $ defined by equation \ref{gamma}.	1354	1385	This function evaluates $ \Gamma (x) $ defined by equation \ref{gamma}.
\begin{equation}	1355	1386	\begin{equation}
\label{gamma}	1356	1387	\label{gamma}
\Gamma (x) = \int _0 ^{\infty} t^{x-1}e^{-t}dt	1357	1388	\Gamma (x) = \int _0 ^{\infty} t^{x-1}e^{-t}dt
\end{equation}	1358	1389	\end{equation}
Note that $n!=\Gamma (n+1)$.	1359	1390	Note that $n!=\Gamma (n+1)$.
	1360	1391
Specific interfaces :\verb'gamma,dgamma,cgamma,zgamm'	1361	1392	Specific interfaces :\verb'gamma,dgamma,cgamma,zgamm'
	1362	1393
\subsubsection{gamr}	1363	1394	\subsubsection{gamr}
\begin{verbatim}	1364	1395	\begin{verbatim}
gamr(x)	1365	1396	gamr(x)
\end{verbatim}	1366	1397	\end{verbatim}
\begin{itemize}	1367	1398	\begin{itemize}
\item x is a real or complex	1368	1399	\item x is a real or complex
\end{itemize}	1369	1400	\end{itemize}
This function evaluates the reciprocal gamma function $gamr(x)= {1 \over \Gamma(x)}$	1370	1401	This function evaluates the reciprocal gamma function $gamr(x)= {1 \over \Gamma(x)}$
	1371	1402
	1372	1403
\subsubsection{alngam}	1373	1404	\subsubsection{alngam}
\begin{verbatim}	1374	1405	\begin{verbatim}
alngam(x)	1375	1406	alngam(x)
\end{verbatim}	1376	1407	\end{verbatim}
\begin{itemize}	1377	1408	\begin{itemize}
\item x is a real or complex	1378	1409	\item x is a real or complex
\end{itemize}	1379	1410	\end{itemize}
This function evaluates $ln(\|\Gamma(x)\|)$	1380	1411	This function evaluates $ln(\|\Gamma(x)\|)$
	1381	1412
Specific interfaces : \verb'alngam,dlngam,clngam,zlngam'	1382	1413	Specific interfaces : \verb'alngam,dlngam,clngam,zlngam'
	1383	1414
	1384	1415
\subsubsection{algams}	1385	1416	\subsubsection{algams}
\begin{verbatim}	1386	1417	\begin{verbatim}
call algams(x,algam,sgngam)	1387	1418	call algams(x,algam,sgngam)
\end{verbatim}	1388	1419	\end{verbatim}
\begin{itemize}	1389	1420	\begin{itemize}
\item x (in) is a real	1390	1421	\item x (in) is a real
\item algam (out) is a real	1391	1422	\item algam (out) is a real
\item sgngam (out) is a real	1392	1423	\item sgngam (out) is a real
\end{itemize}	1393	1424	\end{itemize}
This subroutine evaluates the logarithm of the absolute value of gamma and the sign of gamma.	1394	1425	This subroutine evaluates the logarithm of the absolute value of gamma and the sign of gamma.
$algam=ln(\|\Gamma(x)\|)$ and $sgngam=1.0$ or $-1.0$ according to the sign of $\Gamma(x)$.	1395	1426	$algam=ln(\|\Gamma(x)\|)$ and $sgngam=1.0$ or $-1.0$ according to the sign of $\Gamma(x)$.
	1396	1427
Specific interfaces : \verb'algams,dlgams'	1397	1428	Specific interfaces : \verb'algams,dlgams'
	1398	1429
\subsubsection{gami}	1399	1430	\subsubsection{gami}
\begin{verbatim}	1400	1431	\begin{verbatim}
gami(a,x)	1401	1432	gami(a,x)
\end{verbatim}	1402	1433	\end{verbatim}
\begin{itemize}	1403	1434	\begin{itemize}
\item x is a positive real	1404	1435	\item x is a positive real
\item a is a strictly positive real	1405	1436	\item a is a strictly positive real
\end{itemize}	1406	1437	\end{itemize}
This function evaluates the incomplete gamma function defined by equation \ref{gami}.	1407	1438	This function evaluates the incomplete gamma function defined by equation \ref{gami}.
\begin{equation}	1408	1439	\begin{equation}
\label{gami}	1409	1440	\label{gami}
gami(a,x)=\gamma(a,x)=\int _0 ^x t^{a-1} e^{-t}dt	1410	1441	gami(a,x)=\gamma(a,x)=\int _0 ^x t^{a-1} e^{-t}dt
\end{equation}	1411	1442	\end{equation}
	1412	1443
Specific interfaces : \verb'gami,dgami'	1413	1444	Specific interfaces : \verb'gami,dgami'
	1414	1445
\subsubsection{gamic}	1415	1446	\subsubsection{gamic}
\begin{verbatim}	1416	1447	\begin{verbatim}
gamic(a,x)	1417	1448	gamic(a,x)
\end{verbatim}	1418	1449	\end{verbatim}
\begin{itemize}	1419	1450	\begin{itemize}
\item x is a positive real	1420	1451	\item x is a positive real
\item a is a real	1421	1452	\item a is a real
\end{itemize}	1422	1453	\end{itemize}
This function evaluates the complementary incomplete gamma function defined by equation \ref{gamic}.	1423	1454	This function evaluates the complementary incomplete gamma function defined by equation \ref{gamic}.
\begin{equation}	1424	1455	\begin{equation}
\label{gamic}	1425	1456	\label{gamic}
gamic(a,x)=\Gamma(a,x)=\int _x ^\infty t^{a-1} e^{-t}dt	1426	1457	gamic(a,x)=\Gamma(a,x)=\int _x ^\infty t^{a-1} e^{-t}dt
\end{equation}	1427	1458	\end{equation}
	1428	1459
Specific interfaces : \verb'gamic,dgamic'	1429	1460	Specific interfaces : \verb'gamic,dgamic'
	1430	1461
\subsubsection{gamit}	1431	1462	\subsubsection{gamit}
\begin{verbatim}	1432	1463	\begin{verbatim}
gamit(a,x)	1433	1464	gamit(a,x)
\end{verbatim}	1434	1465	\end{verbatim}
\begin{itemize}	1435	1466	\begin{itemize}
\item x is a positive real	1436	1467	\item x is a positive real
\item a is a real	1437	1468	\item a is a real
\end{itemize}	1438	1469	\end{itemize}
This function evaluates the Tricomi's incomplete gamma function defined by equation \ref{gamit}.	1439	1470	This function evaluates the Tricomi's incomplete gamma function defined by equation \ref{gamit}.
\begin{equation}	1440	1471	\begin{equation}
\label{gamit}	1441	1472	\label{gamit}
gamit(a,x)=\gamma^* (a,x)= {{x^{-a}\gamma(a,x)}\over \Gamma(a)}	1442	1473	gamit(a,x)=\gamma^* (a,x)= {{x^{-a}\gamma(a,x)}\over \Gamma(a)}
\end{equation}	1443	1474	\end{equation}
	1444	1475
Specific interfaces : \verb'gamit,dgamit'	1445	1476	Specific interfaces : \verb'gamit,dgamit'
	1446	1477
	1447	1478
\subsubsection{psi}	1448	1479	\subsubsection{psi}
\begin{verbatim}	1449	1480	\begin{verbatim}
psi(x)	1450	1481	psi(x)
\end{verbatim}	1451	1482	\end{verbatim}
\begin{itemize}	1452	1483	\begin{itemize}
\item x is a real or complex	1453	1484	\item x is a real or complex
\end{itemize}	1454	1485	\end{itemize}
This function evaluates the psi function which is the logarithm derivative of the gamma function as defined in equation \ref{psi}.	1455	1486	This function evaluates the psi function which is the logarithm derivative of the gamma function as defined in equation \ref{psi}.
\begin{equation}	1456	1487	\begin{equation}
\label{psi}	1457	1488	\label{psi}
psi(x)= \psi(x) = {d\over dx} ln(\Gamma(x))	1458	1489	psi(x)= \psi(x) = {d\over dx} ln(\Gamma(x))
\end{equation}	1459	1490	\end{equation}
x must not be zero or a negative integer.	1460	1491	x must not be zero or a negative integer.
	1461	1492
Specific interfaces : \verb'psi,dpsi,cpsi,zpsi'	1462	1493	Specific interfaces : \verb'psi,dpsi,cpsi,zpsi'
	1463	1494
	1464	1495
\subsubsection{poch}	1465	1496	\subsubsection{poch}
\begin{verbatim}	1466	1497	\begin{verbatim}
poch(a,x)	1467	1498	poch(a,x)
\end{verbatim}	1468	1499	\end{verbatim}
\begin{itemize}	1469	1500	\begin{itemize}
\item x is a real	1470	1501	\item x is a real
\item a is a real	1471	1502	\item a is a real
\end{itemize}	1472	1503	\end{itemize}
This function evaluates a generalization of Pochhammer's symbol.	1473	1504	This function evaluates a generalization of Pochhammer's symbol.
	1474	1505
Pochhammer's symbol for n a positive integer is given by equation \ref{poch_int}	1475	1506	Pochhammer's symbol for n a positive integer is given by equation \ref{poch_int}
\begin{equation}	1476	1507	\begin{equation}
\label{poch_int}	1477	1508	\label{poch_int}
(a)_n = a(a-1)(a-2)...(a-n+1)	1478	1509	(a)_n = a(a-1)(a-2)...(a-n+1)
\end{equation}	1479	1510	\end{equation}
	1480	1511
The generalization of Pochhammer's symbol is given by equation \ref{poch}	1481	1512	The generalization of Pochhammer's symbol is given by equation \ref{poch}
\begin{equation}	1482	1513	\begin{equation}
\label{poch}	1483	1514	\label{poch}
poch(a,x)= (a)_x = {\Gamma(a+x) \over \Gamma(a)}	1484	1515	poch(a,x)= (a)_x = {\Gamma(a+x) \over \Gamma(a)}
\end{equation}	1485	1516	\end{equation}
	1486	1517
Specific interfaces : \verb'poch,dpoch'	1487	1518	Specific interfaces : \verb'poch,dpoch'
	1488	1519
	1489	1520
\subsubsection{poch1}	1490	1521	\subsubsection{poch1}
\begin{verbatim}	1491	1522	\begin{verbatim}
poch1(a,x)	1492	1523	poch1(a,x)
\end{verbatim}	1493	1524	\end{verbatim}
\begin{itemize}	1494	1525	\begin{itemize}
\item x is a real	1495	1526	\item x is a real
\item a is a real	1496	1527	\item a is a real
\end{itemize}	1497	1528	\end{itemize}
This function is defined by equation \ref{poch1}. It is usefull for certains situations, especially when x is small.	1498	1529	This function is defined by equation \ref{poch1}. It is usefull for certains situations, especially when x is small.
	1499	1530
\begin{equation}	1500	1531	\begin{equation}
\label{poch1}	1501	1532	\label{poch1}
poch1(a,x)={{(a)_x-1} \over x}	1502	1533	poch1(a,x)={{(a)_x-1} \over x}
\end{equation}	1503	1534	\end{equation}
	1504	1535
Specific interfaces : \verb'poch1,dpoch1'	1505	1536	Specific interfaces : \verb'poch1,dpoch1'
	1506	1537
\subsubsection{beta}	1507	1538	\subsubsection{beta}
\begin{verbatim}	1508	1539	\begin{verbatim}
beta(a,b)	1509	1540	beta(a,b)
\end{verbatim}	1510	1541	\end{verbatim}
\begin{itemize}	1511	1542	\begin{itemize}
\item a,b are real positive or complex	1512	1543	\item a,b are real positive or complex
\end{itemize}	1513	1544	\end{itemize}
This function evaluates $\beta$ function defined by equation \ref{beta}.	1514	1545	This function evaluates $\beta$ function defined by equation \ref{beta}.
\begin{equation}	1515	1546	\begin{equation}
\label{beta}	1516	1547	\label{beta}
beta(a,b)=\beta(a,b)={ {\Gamma(a) \Gamma(b)} \over \Gamma(a+b) }	1517	1548	beta(a,b)=\beta(a,b)={ {\Gamma(a) \Gamma(b)} \over \Gamma(a+b) }
\end{equation}	1518	1549	\end{equation}
	1519	1550
Specific interfaces : \verb'beta,dbeta,cbeta,zbeta'	1520	1551	Specific interfaces : \verb'beta,dbeta,cbeta,zbeta'
	1521	1552
	1522	1553
\subsubsection{albeta}	1523	1554	\subsubsection{albeta}
\begin{verbatim}	1524	1555	\begin{verbatim}
albeta(a,b)	1525	1556	albeta(a,b)
\end{verbatim}	1526	1557	\end{verbatim}
\begin{itemize}	1527	1558	\begin{itemize}
\item a,b are real positive or complex	1528	1559	\item a,b are real positive or complex
\end{itemize}	1529	1560	\end{itemize}
This function evaluates the natural logarithm of beta function : $ln(\beta(a,b))$	1530	1561	This function evaluates the natural logarithm of beta function : $ln(\beta(a,b))$
	1531	1562
Specific interfaces : \verb'albeta,dlbeta,clbeta,zlbeta'	1532	1563	Specific interfaces : \verb'albeta,dlbeta,clbeta,zlbeta'
	1533	1564
\subsubsection{betai}	1534	1565	\subsubsection{betai}
\begin{verbatim}	1535	1566	\begin{verbatim}
betai(x,pin,qin)	1536	1567	betai(x,pin,qin)
\end{verbatim}	1537	1568	\end{verbatim}
\begin{itemize}	1538	1569	\begin{itemize}
\item x is a real in [0,1]	1539	1570	\item x is a real in [0,1]
\item pin and qin are strictly positive real	1540	1571	\item pin and qin are strictly positive real
\end{itemize}	1541	1572	\end{itemize}
This function evaluates the incomplete beta function ratio, that is the probability that a random variable from a beta distribution having parameters pin and qin will be less than or equal to x. It is defined by equation \ref{betai}.	1542	1573	This function evaluates the incomplete beta function ratio, that is the probability that a random variable from a beta distribution having parameters pin and qin will be less than or equal to x. It is defined by equation \ref{betai}.
	1543	1574
\begin{equation}	1544	1575	\begin{equation}
\label{betai}	1545	1576	\label{betai}
betai(x,pin,qin)=I_x(pin,qin)={1 \over \beta(pin,qin)} \int _0 ^x t^{pin-1}(1-t)^{qin-1}dt	1546	1577	betai(x,pin,qin)=I_x(pin,qin)={1 \over \beta(pin,qin)} \int _0 ^x t^{pin-1}(1-t)^{qin-1}dt
\end{equation}	1547	1578	\end{equation}
	1548	1579
Specific interfaces : \verb'betai,dbetai'	1549	1580	Specific interfaces : \verb'betai,dbetai'
	1550	1581
\subsection{Error function and related}	1551	1582	\subsection{Error function and related}
\subsubsection{erf}	1552	1583	\subsubsection{erf}
\begin{verbatim}	1553	1584	\begin{verbatim}
erf(x)	1554	1585	erf(x)
\end{verbatim}	1555	1586	\end{verbatim}
\begin{itemize}	1556	1587	\begin{itemize}
\item x is a real	1557	1588	\item x is a real
\end{itemize}	1558	1589	\end{itemize}
This function evaluates the error function defined by equation \ref{erf}.	1559	1590	This function evaluates the error function defined by equation \ref{erf}.
\begin{equation}	1560	1591	\begin{equation}
\label{erf}	1561	1592	\label{erf}
erf(x)={2\over \sqrt{ \pi}} \int _0 ^x e^{-t^2}dt	1562	1593	erf(x)={2\over \sqrt{ \pi}} \int _0 ^x e^{-t^2}dt
\end{equation}	1563	1594	\end{equation}
	1564	1595
Specific interfaces : \verb'erf,derf'	1565	1596	Specific interfaces : \verb'erf,derf'
	1566	1597
\subsubsection{erfc}	1567	1598	\subsubsection{erfc}
\begin{verbatim}	1568	1599	\begin{verbatim}
erfc(x)	1569	1600	erfc(x)
\end{verbatim}	1570	1601	\end{verbatim}
\begin{itemize}	1571	1602	\begin{itemize}
\item x is a real	1572	1603	\item x is a real
\end{itemize}	1573	1604	\end{itemize}
This function evaluates the complimentary error function defined by equation \ref{erfc}.	1574	1605	This function evaluates the complimentary error function defined by equation \ref{erfc}.
\begin{equation}	1575	1606	\begin{equation}
\label{erfc}	1576	1607	\label{erfc}
erfc(x)={2\over \sqrt{ \pi}} \int _x ^\infty e^{-t^2}dt	1577	1608	erfc(x)={2\over \sqrt{ \pi}} \int _x ^\infty e^{-t^2}dt
\end{equation}	1578	1609	\end{equation}
	1579	1610
Specific interfaces : \verb'erfc,derfc'	1580	1611	Specific interfaces : \verb'erfc,derfc'
	1581	1612
	1582	1613
\subsubsection{daws}	1583	1614	\subsubsection{daws}
\begin{verbatim}	1584	1615	\begin{verbatim}
daws(x)	1585	1616	daws(x)
\end{verbatim}	1586	1617	\end{verbatim}
\begin{itemize}	1587	1618	\begin{itemize}
\item x is a real	1588	1619	\item x is a real
\end{itemize}	1589	1620	\end{itemize}
This function evaluates Dawson's function defined by equation \ref{daws}.	1590	1621	This function evaluates Dawson's function defined by equation \ref{daws}.
\begin{equation}	1591	1622	\begin{equation}
\label{daws}	1592	1623	\label{daws}
daws(x)=e^{-x^2} \int _0 ^x e^{t^2}dt	1593	1624	daws(x)=e^{-x^2} \int _0 ^x e^{t^2}dt
\end{equation}	1594	1625	\end{equation}
	1595	1626
Specific interfaces : \verb'daws,ddaws'	1596	1627	Specific interfaces : \verb'daws,ddaws'
	1597	1628
\subsection{Bessel functions and related}	1598	1629	\subsection{Bessel functions and related}
\subsubsection{bsj0}	1599	1630	\subsubsection{bsj0}
\begin{verbatim}	1600	1631	\begin{verbatim}
bsj0(x)	1601	1632	bsj0(x)
\end{verbatim}	1602	1633	\end{verbatim}
\begin{itemize}	1603	1634	\begin{itemize}
\item x is a real	1604	1635	\item x is a real
\end{itemize}	1605	1636	\end{itemize}
This function evaluates Bessel function of the first kind of order 0 defined by equation \ref{bsj0}.	1606	1637	This function evaluates Bessel function of the first kind of order 0 defined by equation \ref{bsj0}.
\begin{equation}	1607	1638	\begin{equation}
\label{bsj0}	1608	1639	\label{bsj0}
bsj0(x)=J_0(x)= {1 \over \pi} \int _0 ^\pi cos(x sin(\theta)) d\theta	1609	1640	bsj0(x)=J_0(x)= {1 \over \pi} \int _0 ^\pi cos(x sin(\theta)) d\theta
\end{equation}	1610	1641	\end{equation}
	1611	1642
Specific interfaces : \verb'besj0,dbesj0'	1612	1643	Specific interfaces : \verb'besj0,dbesj0'
	1613	1644
\subsubsection{bsj1}	1614	1645	\subsubsection{bsj1}
\begin{verbatim}	1615	1646	\begin{verbatim}
bsj1(x)	1616	1647	bsj1(x)
\end{verbatim}	1617	1648	\end{verbatim}
\begin{itemize}	1618	1649	\begin{itemize}
\item x is a real	1619	1650	\item x is a real
\end{itemize}	1620	1651	\end{itemize}
This function evaluates Bessel function of the first kind of order 1 defined by equation \ref{bsj1}.	1621	1652	This function evaluates Bessel function of the first kind of order 1 defined by equation \ref{bsj1}.
\begin{equation}	1622	1653	\begin{equation}
\label{bsj1}	1623	1654	\label{bsj1}
bsj1(x)=J_1(x)={1 \over \pi} \int _0 ^\pi cos(x sin(\theta)-\theta) d\theta	1624	1655	bsj1(x)=J_1(x)={1 \over \pi} \int _0 ^\pi cos(x sin(\theta)-\theta) d\theta
\end{equation}	1625	1656	\end{equation}
	1626	1657
Specific interfaces : \verb'besj1,dbesj1'	1627	1658	Specific interfaces : \verb'besj1,dbesj1'
	1628	1659
\subsubsection{bsjn}	1629	1660	\subsubsection{bsjn}
\begin{verbatim}	1630	1661	\begin{verbatim}
bsjn(n,x,factor,big)	1631	1662	bsjn(n,x,factor,big)
\end{verbatim}	1632	1663	\end{verbatim}
\begin{itemize}	1633	1664	\begin{itemize}
\item n is an integer	1634	1665	\item n is an integer
\item x is a real	1635	1666	\item x is a real
\item factor is an optional integer	1636	1667	\item factor is an optional integer
\item big is an optional real	1637	1668	\item big is an optional real
\end{itemize}	1638	1669	\end{itemize}
This function evaluates Bessel function of the first kind of order n (plotted in figure \ref{besseljnfamily}). These functions satisfy the recurrent relation \ref {bsjn}.	1639	1670	This function evaluates Bessel function of the first kind of order n (plotted in figure \ref{besseljnfamily}). These functions satisfy the recurrent relation \ref {bsjn}.
\begin{equation}	1640	1671	\begin{equation}
\label{bsjn}	1641	1672	\label{bsjn}
J_{n+1}(x)={{2n}\over x} J_n(x) - J_{n-1}(x)	1642	1673	J_{n+1}(x)={{2n}\over x} J_n(x) - J_{n-1}(x)
\end{equation}	1643	1674	\end{equation}
This relation is directly used in upward direction to compute $J_n(x)$ for $x>n$. However it is unstable for $x<n$, therefore a Miller's Algorithm is used. The principle of this method is to use the reccurent relation downward from an arbitrary higher than n order with an arbitrary seed and then normalize the solution with \ref{bsjnnorm}	1644	1675	This relation is directly used in upward direction to compute $J_n(x)$ for $x>n$. However it is unstable for $x<n$, therefore a Miller's Algorithm is used. The principle of this method is to use the reccurent relation downward from an arbitrary higher than n order with an arbitrary seed and then normalize the solution with \ref{bsjnnorm}
\begin{equation}	1645	1676	\begin{equation}
\label{bsjnnorm}	1646	1677	\label{bsjnnorm}
1=J_0+2J_2+2J_4+2J_6+...	1647	1678	1=J_0+2J_2+2J_4+2J_6+...
\end{equation}	1648	1679	\end{equation}
The optional parameters \verb'factor' and \verb'big' can be used to modify the behaviour of the algorithm. \verb'factor' is used in determining the arbitrary starting order ( an even integer near $n+\sqrt{factor~n}$), the default $factor$ value is 40. \verb'big' is a real determining the threshold for which anti-overflow counter measures has to be taken, default value is $1.10^{10}$	1649	1680	The optional parameters \verb'factor' and \verb'big' can be used to modify the behaviour of the algorithm. \verb'factor' is used in determining the arbitrary starting order ( an even integer near $n+\sqrt{factor~n}$), the default $factor$ value is 40. \verb'big' is a real determining the threshold for which anti-overflow counter measures has to be taken, default value is $1.10^{10}$
	1650	1681
By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsj0(x)$ or $bsj1(x)$ is actually performed.	1651	1682	By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsj0(x)$ or $bsj1(x)$ is actually performed.
	1652	1683
\begin{figure}	1653	1684	\begin{figure}
\begin{center}	1654	1685	\begin{center}
\includegraphics[width=0.9\textwidth]{besselj-mono.pdf}	1655	1686	\includegraphics[width=0.9\textwidth]{besselj-mono.pdf}
\caption{Bessel $J_n$ functions family}	1656	1687	\caption{Bessel $J_n$ functions family}
\label{besseljnfamily}	1657	1688	\label{besseljnfamily}
\end{center}	1658	1689	\end{center}
\end{figure}	1659	1690	\end{figure}
	1660	1691
	1661	1692
Specific interfaces : \verb'besjn,dbesjn'	1662	1693	Specific interfaces : \verb'besjn,dbesjn'
	1663	1694
\subsubsection{bsy0}	1664	1695	\subsubsection{bsy0}
\begin{verbatim}	1665	1696	\begin{verbatim}
bsy0(x)	1666	1697	bsy0(x)
\end{verbatim}	1667	1698	\end{verbatim}
\begin{itemize}	1668	1699	\begin{itemize}
\item x is a strictly positive real	1669	1700	\item x is a strictly positive real
\end{itemize}	1670	1701	\end{itemize}
This function evaluates the Bessel function of the second kind of order 0 defined by equation \ref{bsy0}	1671	1702	This function evaluates the Bessel function of the second kind of order 0 defined by equation \ref{bsy0}
\begin{equation}	1672	1703	\begin{equation}
\label{bsy0}	1673	1704	\label{bsy0}
bsy0(x)=Y_0(x)={1 \over \pi} \int _0 ^\pi sin(x sin(\theta))d\theta -{2 \over \pi} \int _0 ^\infty e^{-x sinh(t)}dt	1674	1705	bsy0(x)=Y_0(x)={1 \over \pi} \int _0 ^\pi sin(x sin(\theta))d\theta -{2 \over \pi} \int _0 ^\infty e^{-x sinh(t)}dt
\end{equation}	1675	1706	\end{equation}
	1676	1707
Specific interfaces : \verb'besy0,dbesy0'	1677	1708	Specific interfaces : \verb'besy0,dbesy0'
	1678	1709
\subsubsection{bsy1}	1679	1710	\subsubsection{bsy1}
\begin{verbatim}	1680	1711	\begin{verbatim}
bsy1(x)	1681	1712	bsy1(x)
\end{verbatim}	1682	1713	\end{verbatim}
\begin{itemize}	1683	1714	\begin{itemize}
\item x is a strictly positive real	1684	1715	\item x is a strictly positive real
\end{itemize}	1685	1716	\end{itemize}
This function evaluates the Bessel function of the second kind of order 1 defined by equation \ref{bsy1}.	1686	1717	This function evaluates the Bessel function of the second kind of order 1 defined by equation \ref{bsy1}.
\begin{equation}	1687	1718	\begin{equation}
\label{bsy1}	1688	1719	\label{bsy1}
bsy1(x)=Y_1(x)=-{1 \over \pi} \int _0 ^\pi sin (\theta - x sin(\theta)) d\theta	1689	1720	bsy1(x)=Y_1(x)=-{1 \over \pi} \int _0 ^\pi sin (\theta - x sin(\theta)) d\theta
- {1 \over \pi} \int _0 ^\infty (e^t -e^{-t})e^{-x sinh(t)}dt	1690	1721	- {1 \over \pi} \int _0 ^\infty (e^t -e^{-t})e^{-x sinh(t)}dt
\end{equation}	1691	1722	\end{equation}
	1692	1723
Specific interfaces : \verb'besy1,dbesy1'	1693	1724	Specific interfaces : \verb'besy1,dbesy1'
	1694	1725
	1695	1726
\subsubsection{bsyn}	1696	1727	\subsubsection{bsyn}
\begin{verbatim}	1697	1728	\begin{verbatim}
bsyn(n,x)	1698	1729	bsyn(n,x)
\end{verbatim}	1699	1730	\end{verbatim}
\begin{itemize}	1700	1731	\begin{itemize}
\item n is an integer	1701	1732	\item n is an integer
\item x is a strictly positive real	1702	1733	\item x is a strictly positive real
\end{itemize}	1703	1734	\end{itemize}
This function evaluates the Bessel function of the second kind of order n (plotted in figure \ref{besselynfamily}). These functions satisfy the recurrent relation \ref {bsyn}.	1704	1735	This function evaluates the Bessel function of the second kind of order n (plotted in figure \ref{besselynfamily}). These functions satisfy the recurrent relation \ref {bsyn}.
\begin{equation}	1705	1736	\begin{equation}
\label{bsyn}	1706	1737	\label{bsyn}
Y_{n+1}(x)={{2n}\over x} Y_n(x) - Y_{n-1}(x)	1707	1738	Y_{n+1}(x)={{2n}\over x} Y_n(x) - Y_{n-1}(x)
\end{equation}	1708	1739	\end{equation}
This recurrent relation is directly used in the upward direction to compute $Y_n(x)$.	1709	1740	This recurrent relation is directly used in the upward direction to compute $Y_n(x)$.
	1710	1741
By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsy0(x)$ or $bsy1(x)$ is actually performed.	1711	1742	By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsy0(x)$ or $bsy1(x)$ is actually performed.
	1712	1743
\begin{figure}	1713	1744	\begin{figure}
\begin{center}	1714	1745	\begin{center}
\includegraphics[width=0.9\textwidth]{bessely-mono.pdf}	1715	1746	\includegraphics[width=0.9\textwidth]{bessely-mono.pdf}
\caption{Bessel $Y_n$ functions family}	1716	1747	\caption{Bessel $Y_n$ functions family}
\label{besselynfamily}	1717	1748	\label{besselynfamily}
\end{center}	1718	1749	\end{center}
\end{figure}	1719	1750	\end{figure}
	1720	1751
	1721	1752
Specific interfaces : \verb'besyn,dbesyn'	1722	1753	Specific interfaces : \verb'besyn,dbesyn'
	1723	1754
\subsubsection{bsi0}	1724	1755	\subsubsection{bsi0}
\begin{verbatim}	1725	1756	\begin{verbatim}
bsi0(x)	1726	1757	bsi0(x)
\end{verbatim}	1727	1758	\end{verbatim}
\begin{itemize}	1728	1759	\begin{itemize}
\item x is a real	1729	1760	\item x is a real
\end{itemize}	1730	1761	\end{itemize}
This function evaluates the Bessel function of the third kind of order 0 defined by equation \ref{bsi0}.	1731	1762	This function evaluates the Bessel function of the third kind of order 0 defined by equation \ref{bsi0}.
\begin{equation}	1732	1763	\begin{equation}

fvn_test/test_akima.f90

Diff comments View file @ 2919a9e

program akima	1	1	program akima
use fvn	2	2	use fvn_interpol
implicit none	3	3	implicit none
integer :: nbpoints,nppoints,i	4	4	integer :: nbpoints,nppoints,i
real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d	5	5	real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
real(8),dimension(:,:),allocatable :: coeff_fvn_d	6	6	real(8),dimension(:,:),allocatable :: coeff_fvn_d
real(8) :: xstep_d,xp_d,ty_d,fvn_y_d	7	7	real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
open(2,file='fvn_akima_double.dat')	8	8	open(2,file='fvn_akima_double.dat')
open(3,file='fvn_akima_breakpoints_double.dat')	9	9	open(3,file='fvn_akima_breakpoints_double.dat')
nbpoints=30	10	10	nbpoints=30
allocate(x_d(nbpoints))	11	11	allocate(x_d(nbpoints))
allocate(y_d(nbpoints))	12	12	allocate(y_d(nbpoints))
allocate(breakpoints_d(nbpoints))	13	13	allocate(breakpoints_d(nbpoints))
allocate(coeff_fvn_d(4,nbpoints))	14	14	allocate(coeff_fvn_d(4,nbpoints))
xstep_d=20./dfloat(nbpoints)	15	15	xstep_d=20./dfloat(nbpoints)
do i=1,nbpoints	16	16	do i=1,nbpoints
x_d(i)=-10.+dfloat(i)*xstep_d	17	17	x_d(i)=-10.+dfloat(i)*xstep_d
y_d(i)=dsin(x_d(i))	18	18	y_d(i)=dsin(x_d(i))
write(3,44) x_d(i),y_d(i)	19	19	write(3,44) x_d(i),y_d(i)
end do	20	20	end do
close(3)	21	21	close(3)
call fvn_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)	22	22	call fvn_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
nppoints=1000	23	23	nppoints=1000
xstep_d=22./dfloat(nppoints)	24	24	xstep_d=22./dfloat(nppoints)
do i=1,nppoints	25	25	do i=1,nppoints
xp_d=-11.+dfloat(i)*xstep_d	26	26	xp_d=-11.+dfloat(i)*xstep_d
ty_d=dsin(xp_d)	27	27	ty_d=dsin(xp_d)
fvn_y_d=fvn_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)	28	28	fvn_y_d=fvn_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
write(2,44) xp_d,ty_d,fvn_y_d	29	29	write(2,44) xp_d,ty_d,fvn_y_d
end do	30	30	end do
close(2)	31	31	close(2)
deallocate(coeff_fvn_d,breakpoints_d,y_d,x_d)	32	32	deallocate(coeff_fvn_d,breakpoints_d,y_d,x_d)
write(6,*) "All done, plot results with gnuplot using command :"	33	33	write(6,*) "All done, plot results with gnuplot using command :"
write(6,*) "pl 'fvn_akima_double.dat' u 1:2 w l,'fvn_akima_breakpoints_double.dat' w p"	34	34	write(6,*) "pl 'fvn_akima_double.dat' u 1:2 w l,'fvn_akima_breakpoints_double.dat' w p"
	35	35
44 FORMAT(4(1X,1PE22.14))	36	36	44 FORMAT(4(1X,1PE22.14))
end program	37	37	end program

fvn_test/test_det.f90

Diff comments View file @ 2919a9e

program test_det	1	1	program test_det
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
real(8),dimension(3,3) :: a	4	4	real(8),dimension(3,3) :: a
real(8) :: deta	5	5	real(8) :: deta
integer :: status,i	6	6	integer :: status,i
	7	7
call init_random_seed()	8	8	call init_random_seed()
call random_number(a)	9	9	call random_number(a)
a=a*100	10	10	a=a*100
deta=fvn_det(3,a,status)	11	11	deta=fvn_det(3,a,status)
do i=1,3	12	12	do i=1,3
write (*,'(3(f15.5))') a(i,:)	13	13	write (*,'(3(f15.5))') a(i,:)
end do	14	14	end do
write (,)	15	15	write (,)
write (,) "Det = ",deta,"Status = ",status	16	16	write (,) "Det = ",deta,"Status = ",status
end program	17	17	end program

fvn_test/test_integ.f90

Diff comments View file @ 2919a9e

program integ	1	1	program integ
use fvn	2	2	use fvn_integ
implicit none	3	3	implicit none
real(8), external :: f1,f2,g,h	4	4	real(8), external :: f1,f2,g,h
real(8) :: a,b,epsabs,epsrel,abserr,res	5	5	real(8) :: a,b,epsabs,epsrel,abserr,res
integer :: key,ier	6	6	integer :: key,ier
a=0.	7	7	a=0.
b=1.	8	8	b=1.
epsabs=1d-8	9	9	epsabs=1d-8
epsrel=1d-8	10	10	epsrel=1d-8
key=2	11	11	key=2
call fvn_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier)	12	12	call fvn_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier)
write(,) "Integration of x*x between 0 and 1 : "	13	13	write(,) "Integration of x*x between 0 and 1 : "
write(,) res	14	14	write(,) res
call fvn_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier)	15	15	call fvn_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier)
write(,) "Integration of x*y between 0 and 1 on both x and y : "	16	16	write(,) "Integration of x*y between 0 and 1 on both x and y : "
write(,) res	17	17	write(,) res
	18	18
end program	19	19	end program
function f1(x)	20	20	function f1(x)
implicit none	21	21	implicit none
real(8) :: x,f1	22	22	real(8) :: x,f1
f1=x*x	23	23	f1=x*x
end function	24	24	end function
function f2(x,y)	25	25	function f2(x,y)
implicit none	26	26	implicit none
real(8) :: x,y,f2	27	27	real(8) :: x,y,f2
f2=x*y	28	28	f2=x*y
end function	29	29	end function
function g(x)	30	30	function g(x)
implicit none	31	31	implicit none
real(8) :: x,g	32	32	real(8) :: x,g
g=0.	33	33	g=0.
end function	34	34	end function
function h(x)	35	35	function h(x)
implicit none	36	36	implicit none
real(8) :: x,h	37	37	real(8) :: x,h
h=1.	38	38	h=1.
end function	39	39	end function

fvn_test/test_inter1d.f90

Diff comments View file @ 2919a9e

program inter1d	1	1	program inter1d
use fvn	2	2	use fvn_interpol
implicit none	3	3	implicit none
integer(kind=4),parameter :: ndata=33	4	4	integer(kind=4),parameter :: ndata=33
integer(kind=4) :: i,nout	5	5	integer(kind=4) :: i,nout
real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)	6	6	real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
real(kind=8) ::tv	7	7	real(kind=8) ::tv
intrinsic sin	8	8	intrinsic sin
f(x)=sin(x)	9	9	f(x)=sin(x)
xdata(1)=0.	10	10	xdata(1)=0.
fdata(1)=f(xdata(1))	11	11	fdata(1)=f(xdata(1))
h=1./32.	12	12	h=1./32.
do i=2,ndata	13	13	do i=2,ndata
xdata(i)=xdata(i-1)+h	14	14	xdata(i)=xdata(i-1)+h
fdata(i)=f(xdata(i))	15	15	fdata(i)=f(xdata(i))
end do	16	16	end do
call init_random_seed()	17	17	call init_random_seed()
call random_number(x)	18	18	call random_number(x)
q=fvn_quad_interpol(x,ndata,xdata,fdata)	19	19	q=fvn_quad_interpol(x,ndata,xdata,fdata)
tv=f(x)	20	20	tv=f(x)
write(*,'("x y z ",1(f8.5))') x	21	21	write(*,'("x y z ",1(f8.5))') x
write(*,'("Calculated (real) value :",f8.5)') tv	22	22	write(*,'("Calculated (real) value :",f8.5)') tv
write(*,'("fvn interpolation : ",f8.5)') q	23	23	write(*,'("fvn interpolation : ",f8.5)') q
write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)	24	24	write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
end program	25	25	end program

fvn_test/test_inter2d.f90

Diff comments View file @ 2919a9e

program inter2d	1	1	program inter2d
use fvn	2	2	use fvn_interpol
implicit none	3	3	implicit none
integer(kind=4),parameter :: nx=21,ny=42	4	4	integer(kind=4),parameter :: nx=21,ny=42
integer(kind=4) :: i,j	5	5	integer(kind=4) :: i,j
real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)	6	6	real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
real(kind=8) :: tv	7	7	real(kind=8) :: tv
intrinsic dble,sin	8	8	intrinsic dble,sin
f(x,y)=sin(x+2.*y)	9	9	f(x,y)=sin(x+2.*y)
do i=1,nx	10	10	do i=1,nx
xdata(i)=dble(i-1)/dble(nx-1)	11	11	xdata(i)=dble(i-1)/dble(nx-1)
end do	12	12	end do
do i=1,ny	13	13	do i=1,ny
ydata(i)=dble(i-1)/dble(ny-1)	14	14	ydata(i)=dble(i-1)/dble(ny-1)
end do	15	15	end do
do i=1,nx	16	16	do i=1,nx
do j=1,ny	17	17	do j=1,ny
fdata(i,j)=f(xdata(i),ydata(j))	18	18	fdata(i,j)=f(xdata(i),ydata(j))
end do	19	19	end do
end do	20	20	end do
call init_random_seed()	21	21	call init_random_seed()
call random_number(x)	22	22	call random_number(x)
call random_number(y)	23	23	call random_number(y)
q=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)	24	24	q=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
tv=f(x,y)	25	25	tv=f(x,y)
write(*,'("x y z ",2(f8.5))') x,y	26	26	write(*,'("x y z ",2(f8.5))') x,y
write(*,'("Calculated (real) value :",f8.5)') tv	27	27	write(*,'("Calculated (real) value :",f8.5)') tv
write(*,'("fvn interpolation : ",f8.5)') q	28	28	write(*,'("fvn interpolation : ",f8.5)') q
write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)	29	29	write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
end program	30	30	end program

fvn_test/test_inter3d.f90

Diff comments View file @ 2919a9e

program test_inter3d	1	1	program test_inter3d
use fvn	2	2	use fvn_interpol
implicit none	3	3	implicit none
integer(kind=4),parameter :: nx=21,ny=42,nz=18	4	4	integer(kind=4),parameter :: nx=21,ny=42,nz=18
integer(kind=4) :: i,j,k	5	5	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)	6	6	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	7	7	real(kind=8) :: tv
intrinsic dble,sin	8	8	intrinsic dble,sin
f(x,y,z)=sin(x+2.y+3.z)	9	9	f(x,y,z)=sin(x+2.y+3.z)
do i=1,nx	10	10	do i=1,nx
xdata(i)=2.*(dble(i-1)/dble(nx-1))	11	11	xdata(i)=2.*(dble(i-1)/dble(nx-1))
end do	12	12	end do
do i=1,ny	13	13	do i=1,ny
ydata(i)=2.*(dble(i-1)/dble(ny-1))	14	14	ydata(i)=2.*(dble(i-1)/dble(ny-1))
end do	15	15	end do
do i=1,nz	16	16	do i=1,nz
zdata(i)=2.*(dble(i-1)/dble(nz-1))	17	17	zdata(i)=2.*(dble(i-1)/dble(nz-1))
end do	18	18	end do
do i=1,nx	19	19	do i=1,nx
do j=1,ny	20	20	do j=1,ny
do k=1,nz	21	21	do k=1,nz
fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))	22	22	fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
end do	23	23	end do
end do	24	24	end do
end do	25	25	end do
call init_random_seed()	26	26	call init_random_seed()
call random_number(x)	27	27	call random_number(x)
call random_number(y)	28	28	call random_number(y)
call random_number(z)	29	29	call random_number(z)
q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)	30	30	q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
tv=f(x,y,z)	31	31	tv=f(x,y,z)
write(*,'("x y z ",3(f8.5))') x,y,z	32	32	write(*,'("x y z ",3(f8.5))') x,y,z
write(*,'("Calculated (real) value :",f8.5)') tv	33	33	write(*,'("Calculated (real) value :",f8.5)') tv
write(*,'("fvn interpolation : ",f8.5)') q	34	34	write(*,'("fvn interpolation : ",f8.5)') q
write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)	35	35	write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
end program	36	36	end program
	37	37

fvn_test/test_lsp.f90

Diff comments View file @ 2919a9e

program lsp	1	1	program lsp
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
integer,parameter :: npoints=13,deg=3	4	4	integer,parameter :: npoints=13,deg=3
integer :: status,i	5	5	integer :: status,i
real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc	6	6	real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
real(kind=8) :: coeff(deg+1)	7	7	real(kind=8) :: coeff(deg+1)
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. /)	8	8	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 /)	9	9	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')	10	10	open(2,file='fvn_lsp_double_mesure.dat')
open(3,file='fvn_lsp_double_poly.dat')	11	11	open(3,file='fvn_lsp_double_poly.dat')
do i=1,npoints	12	12	do i=1,npoints
write(2,44) xm(i),ym(i)	13	13	write(2,44) xm(i),ym(i)
end do	14	14	end do
close(2)	15	15	close(2)
call fvn_lspoly(npoints,xm,ym,deg,coeff,status)	16	16	call fvn_lspoly(npoints,xm,ym,deg,coeff,status)
xstep=(xm(npoints)-xm(1))/1000.	17	17	xstep=(xm(npoints)-xm(1))/1000.
do i=1,1000	18	18	do i=1,1000
xc=xm(1)+(i-1)*xstep	19	19	xc=xm(1)+(i-1)*xstep
yc=poly(xc,coeff)	20	20	yc=poly(xc,coeff)
write(3,44) xc,yc	21	21	write(3,44) xc,yc
end do	22	22	end do
close(3)	23	23	close(3)
write(,) "All done, plot results with gnuplot using command :"	24	24	write(,) "All done, plot results with gnuplot using command :"
write(,) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lsp_double_poly.dat' u 1:2 w l"	25	25	write(,) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lsp_double_poly.dat' u 1:2 w l"
44 FORMAT(4(1X,1PE22.14))	26	26	44 FORMAT(4(1X,1PE22.14))
contains	27	27	contains
function poly(x,coeff)	28	28	function poly(x,coeff)
implicit none	29	29	implicit none
real(8) :: x	30	30	real(8) :: x
real(8) :: coeff(deg+1)	31	31	real(8) :: coeff(deg+1)
real(8) :: poly	32	32	real(8) :: poly
integer :: i	33	33	integer :: i
poly=0.	34	34	poly=0.
do i=1,deg+1	35	35	do i=1,deg+1
poly=poly+coeff(i)x*(i-1)	36	36	poly=poly+coeff(i)x*(i-1)
end do	37	37	end do
end function	38	38	end function
end program	39	39	end program
	40	40

fvn_test/test_matcon.f90

Diff comments View file @ 2919a9e

program test_matcon	1	1	program test_matcon
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
real(8),dimension(3,3) :: a	4	4	real(8),dimension(3,3) :: a
real(8) :: rcond	5	5	real(8) :: rcond
integer :: status,i	6	6	integer :: status,i
call init_random_seed()	7	7	call init_random_seed()
call random_number(a)	8	8	call random_number(a)
a=a*100	9	9	a=a*100
call fvn_matcon(3,a,rcond,status)	10	10	call fvn_matcon(3,a,rcond,status)
write(,) "Reasonnably conditionned matrix"	11	11	write(,) "Reasonnably conditionned matrix"
do i=1,3	12	12	do i=1,3
write (*,'(3(e12.5))') a(i,:)	13	13	write (*,'(3(e12.5))') a(i,:)
end do	14	14	end do
write (,)	15	15	write (,)
write (,) "Cond = ",rcond	16	16	write (,) "Cond = ",rcond
write (,)	17	17	write (,)
write (,)	18	18	write (,)
a(1,1)=a(1,1)*1d9	19	19	a(1,1)=a(1,1)*1d9
write(,) "Badly conditionned matrix"	20	20	write(,) "Badly conditionned matrix"
do i=1,3	21	21	do i=1,3
write (*,'(3(e12.5))') a(i,:)	22	22	write (*,'(3(e12.5))') a(i,:)
end do	23	23	end do
call fvn_matcon(3,a,rcond,status)	24	24	call fvn_matcon(3,a,rcond,status)
write (,)	25	25	write (,)
write (,) "Cond = ",rcond	26	26	write (,) "Cond = ",rcond
	27	27
end program	28	28	end program

fvn_test/test_matev.f90

Diff comments View file @ 2919a9e

program test_matev	1	1	program test_matev
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
real(8),dimension(3,3) :: a	4	4	real(8),dimension(3,3) :: a
complex(8),dimension(3) :: evala	5	5	complex(8),dimension(3) :: evala
complex(8),dimension(3,3) :: eveca	6	6	complex(8),dimension(3,3) :: eveca
integer :: status,i,j	7	7	integer :: status,i,j
call init_random_seed()	8	8	call init_random_seed()
call random_number(a)	9	9	call random_number(a)
a=a*100	10	10	a=a*100
call fvn_matev(3,a,evala,eveca,status)	11	11	call fvn_matev(3,a,evala,eveca,status)
	12	12
write(,) "The matrix :"	13	13	write(,) "The matrix :"
write (*,'(3(e12.5))') a	14	14	write (*,'(3(e12.5))') a
write (,)	15	15	write (,)
do i=1,3	16	16	do i=1,3
write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)	17	17	write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
write(,) "Associated Eigenvector :"	18	18	write(,) "Associated Eigenvector :"
do j=1,3	19	19	do j=1,3
write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)	20	20	write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
end do	21	21	end do
write(,)	22	22	write(,)
end do	23	23	end do
end program	24	24	end program

fvn_test/test_matinv.f90

Diff comments View file @ 2919a9e

program test_matinv	1	1	program test_matinv
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
integer(4), parameter :: n=3	4	4	integer(4), parameter :: n=3
integer(4) :: status,i	5	5	integer(4) :: status,i
complex(8),dimension(n,n) :: m,im,prod	6	6	complex(8),dimension(n,n) :: m,im,prod
real(8),dimension(n,n) :: rtmp,itmp	7	7	real(8),dimension(n,n) :: rtmp,itmp
character(len=80) :: fmreal,fmcmplx	8	8	character(len=80) :: fmreal,fmcmplx
	9	9
fmcmplx='(3("(",f8.5,",",f8.5,") "))'	10	10	fmcmplx='(3("(",f8.5,",",f8.5,") "))'
	11	11
! initialize pseudo random generator	12	12	! initialize pseudo random generator
call init_random_seed()	13	13	call init_random_seed()
! fill real and imaginary part	14	14	! fill real and imaginary part
call random_number(rtmp)	15	15	call random_number(rtmp)
call random_number(itmp)	16	16	call random_number(itmp)
! create the complex matrix (fvn_i is defined in the fvn module)	17	17	! create the complex matrix (fvn_i is defined in the fvn module)
m=rtmp+fvn_i*itmp	18	18	m=rtmp+fvn_i*itmp
write(,) "Matrix M"	19	19	write(,) "Matrix M"
do i=1,n	20	20	do i=1,n
write(*,fmcmplx) m(i,:)	21	21	write(*,fmcmplx) m(i,:)
end do	22	22	end do
	23	23
! Invertion	24	24	! Invertion
call fvn_matinv(n,m,im)	25	25	call fvn_matinv(n,m,im)
write(,) "Inverse of M"	26	26	write(,) "Inverse of M"
do i=1,n	27	27	do i=1,n
write(*,fmcmplx) im(i,:)	28	28	write(*,fmcmplx) im(i,:)
end do	29	29	end do
	30	30
! Result should be identity matrix	31	31	! Result should be identity matrix
write(,) "Product of M and inverse of M :"	32	32	write(,) "Product of M and inverse of M :"
prod=matmul(m,im)	33	33	prod=matmul(m,im)
do i=1,n	34	34	do i=1,n

fvn_test/test_muller.f90

Diff comments View file @ 2919a9e

program muller	1	1	program muller
use fvn	2	2	use fvn_misc
implicit none	3	3	implicit none
integer :: i,info	4	4	integer :: i,info
complex(8),dimension(10) :: roots	5	5	complex(8),dimension(10) :: roots
integer,dimension(10) :: infer	6	6	integer,dimension(10) :: infer
complex(8), external :: f	7	7	complex(8), external :: f
	8	8
real(8) :: eps1	9	9	real(8) :: eps1
eps1=1.d-10	10	10	eps1=1.d-10
call fvn_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)	11	11	call fvn_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
write(,) "Error code :",info	12	12	write(,) "Error code :",info
do i=1,10	13	13	do i=1,10
write(,) roots(i),infer(i)	14	14	write(,) roots(i),infer(i)
enddo	15	15	enddo
end program	16	16	end program
	17	17
function f(x)	18	18	function f(x)
complex(8) :: x,f	19	19	complex(8) :: x,f
f=x**10-1	20	20	f=x**10-1

fvn_test/test_operators.f90

Diff comments View file @ 2919a9e

program test_matinv	1	1	program test_matinv
use fvn	2	2	use fvn_linear
implicit none	3	3	implicit none
	4	4
integer, parameter :: n=3	5	5	integer, parameter :: n=3
complex(8),dimension(n,n) :: m1,m2,res	6	6	complex(8),dimension(n,n) :: m1,m2,res
real(8),dimension(n,n) :: rtmp,itmp	7	7	real(8),dimension(n,n) :: rtmp,itmp
character(len=80) :: fmcmplx	8	8	character(len=80) :: fmcmplx
integer :: i	9	9	integer :: i
	10	10
fmcmplx='(3("(",f8.5,",",f8.5,") "))'	11	11	fmcmplx='(3("(",f8.5,",",f8.5,") "))'
! initialize pseudo random generator	12	12	! initialize pseudo random generator
!call init_random_seed()	13	13	!call init_random_seed()
! fill real and imaginary part	14	14	! fill real and imaginary part
call random_number(rtmp)	15	15	call random_number(rtmp)
call random_number(itmp)	16	16	call random_number(itmp)
! create the complex matrix (fvn_i is defined in the fvn module)	17	17	! create the complex matrix (fvn_i is defined in the fvn module)
m1=rtmp+fvn_i*itmp	18	18	m1=rtmp+fvn_i*itmp
write(,) "Matrix M1"	19	19	write(,) "Matrix M1"
do i=1,n	20	20	do i=1,n
write(*,fmcmplx) m1(i,:)	21	21	write(*,fmcmplx) m1(i,:)
end do	22	22	end do
	23	23
call random_number(rtmp)	24	24	call random_number(rtmp)
call random_number(itmp)	25	25	call random_number(itmp)
m2=rtmp+fvn_i*itmp	26	26	m2=rtmp+fvn_i*itmp
write(,)	27	27	write(,)
write(,) "Matrix M2"	28	28	write(,) "Matrix M2"
do i=1,n	29	29	do i=1,n
write(*,fmcmplx) m2(i,:)	30	30	write(*,fmcmplx) m2(i,:)
end do	31	31	end do
	32	32
	33	33
write(,)	34	34	write(,)
write(,) "M1.x.M2"	35	35	write(,) "M1.x.M2"
res=m1.x.m2	36	36	res=m1.x.m2
call write_res()	37	37	call write_res()
write(,)	38	38	write(,)
write(,) "M1.x.M2 standard"	39	39	write(,) "M1.x.M2 standard"
res=matmul(m1,m2)	40	40	res=matmul(m1,m2)
call write_res()	41	41	call write_res()
	42	42
write(,)	43	43	write(,)
write(,) ".i.M1"	44	44	write(,) ".i.M1"
res=.i.m1	45	45	res=.i.m1
call write_res()	46	46	call write_res()
write(,)	47	47	write(,)
write(,) ".i.M1 standard"	48	48	write(,) ".i.M1 standard"
call fvn_matinv(3,m1,res)	49	49	call fvn_matinv(3,m1,res)
call write_res()	50	50	call write_res()
	51	51
write(,)	52	52	write(,)
write(,) "M1.ix.M2"	53	53	write(,) "M1.ix.M2"
res=m1.ix.m2	54	54	res=m1.ix.m2
call write_res()	55	55	call write_res()
write(,)	56	56	write(,)
write(,) "M1.ix.M2 standard"	57	57	write(,) "M1.ix.M2 standard"
call fvn_matinv(3,m1,res)	58	58	call fvn_matinv(3,m1,res)
res=matmul(res,m2)	59	59	res=matmul(res,m2)
call write_res()	60	60	call write_res()
	61	61
write(,)	62	62	write(,)
write(,) "M1.xi.M2"	63	63	write(,) "M1.xi.M2"
res=m1.xi.m2	64	64	res=m1.xi.m2
call write_res()	65	65	call write_res()
write(,)	66	66	write(,)
write(,) "M1.xi.M2 standard"	67	67	write(,) "M1.xi.M2 standard"
res=m1.xi.m2	68	68	res=m1.xi.m2
call fvn_matinv(3,m2,res)	69	69	call fvn_matinv(3,m2,res)
res=matmul(m1,res)	70	70	res=matmul(m1,res)
call write_res()	71	71	call write_res()
	72	72
write(,)	73	73	write(,)
write(,) ".t.M1"	74	74	write(,) ".t.M1"
res=.t.m1	75	75	res=.t.m1
call write_res()	76	76	call write_res()
write(,)	77	77	write(,)
write(,) ".t.M1 standard"	78	78	write(,) ".t.M1 standard"
res=transpose(m1)	79	79	res=transpose(m1)
call write_res()	80	80	call write_res()
	81	81
write(,)	82	82	write(,)
write(,) "M1.tx.M2"	83	83	write(,) "M1.tx.M2"
res=m1.tx.m2	84	84	res=m1.tx.m2
call write_res()	85	85	call write_res()
write(,)	86	86	write(,)
write(,) "M1.tx.M2 standard"	87	87	write(,) "M1.tx.M2 standard"
res=matmul(transpose(m1),m2)	88	88	res=matmul(transpose(m1),m2)
call write_res()	89	89	call write_res()
	90	90
write(,)	91	91	write(,)
write(,) "M1.xt.M2"	92	92	write(,) "M1.xt.M2"
res=m1.xt.m2	93	93	res=m1.xt.m2
call write_res()	94	94	call write_res()
write(,)	95	95	write(,)
write(,) "M1.xt.M2 standard"	96	96	write(,) "M1.xt.M2 standard"
res=matmul(m1,transpose(m2))	97	97	res=matmul(m1,transpose(m2))
call write_res()	98	98	call write_res()
	99	99
write(,)	100	100	write(,)
write(,) ".h.M1"	101	101	write(,) ".h.M1"
res=.h.m1	102	102	res=.h.m1
call write_res()	103	103	call write_res()
write(,)	104	104	write(,)
write(,) ".h.M1 standard"	105	105	write(,) ".h.M1 standard"
res=transpose(conjg(m1))	106	106	res=transpose(conjg(m1))
call write_res()	107	107	call write_res()
	108	108
write(,)	109	109	write(,)
write(,) "M1.hx.M2"	110	110	write(,) "M1.hx.M2"
res=m1.hx.m2	111	111	res=m1.hx.m2
call write_res()	112	112	call write_res()
write(,)	113	113	write(,)
write(,) "M1.hx.M2 standard"	114	114	write(,) "M1.hx.M2 standard"
res=matmul(transpose(conjg(m1)),m2)	115	115	res=matmul(transpose(conjg(m1)),m2)
call write_res()	116	116	call write_res()
	117	117
write(,)	118	118	write(,)
write(,) "M1.xh.M2"	119	119	write(,) "M1.xh.M2"

fvn_test/test_sparse.f90

Diff comments View file @ 2919a9e

program test_sparse	1	1	program test_sparse
use fvn	2	2	use fvn_sparse
implicit none	3	3	implicit none
integer(4), parameter :: nz=12	4	4	integer(4), parameter :: nz=12
integer(4), parameter :: n=5	5	5	integer(4), parameter :: n=5
real(8),dimension(nz) :: A	6	6	real(8),dimension(nz) :: A
real(8),dimension(n,n) :: As	7	7	real(8),dimension(n,n) :: As
integer(4),dimension(nz) :: Ti,Tj	8	8	integer(4),dimension(nz) :: Ti,Tj
real(8),dimension(n) :: B,x	9	9	real(8),dimension(n) :: B,x
integer(4) :: status,i	10	10	integer(4) :: status,i
! Description of the matrix in triplet form	11	11	! Description of the matrix in triplet form
A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)	12	12	A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
B = (/ 8., 45., -3., 3., 19./)	13	13	B = (/ 8., 45., -3., 3., 19./)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)	14	14	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 /)	15	15	Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	16	16
! Reconstruction of the matrix in standard form	17	17	! Reconstruction of the matrix in standard form
! just needed for printing the matrix here	18	18	! just needed for printing the matrix here
As=0.	19	19	As=0.
do i=1,nz	20	20	do i=1,nz
As(Ti(i),Tj(i))=A(i)	21	21	As(Ti(i),Tj(i))=A(i)
end do	22	22	end do
write(,) "Matrix in standard representation :"	23	23	write(,) "Matrix in standard representation :"
do i=1,5	24	24	do i=1,5
write(*,'(5f8.4)') As(i,:)	25	25	write(*,'(5f8.4)') As(i,:)
end do	26	26	end do
write(,)	27	27	write(,)
write(*,'("Right hand side :",5f8.4)') B	28	28	write(*,'("Right hand side :",5f8.4)') B
	29	29
!specific routine that will be used here	30	30	!specific routine that will be used here
!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	31	31	!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)	32	32	call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(*,'("Solution :",5f8.4)') x	33	33	write(*,'("Solution :",5f8.4)') x
write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)	34	34	write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)