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