lxsd / fvn

Download as
Browse Code ยป

Commit 9158e74d6be4ac1ce78ef091480d45ff0793becc

Authored by daniau
1 parent 15adbf52ba
Exists in master and in 3 other branches geevx, geevx_2020, multi

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

Showing 3 changed files with 867 additions and 0 deletions Inline Diff

No preview for this file type
File was created 1 %\documentclass[a4paper,10pt]{article}
2 \documentclass[a4paper,english]{article}
3
4 \usepackage[utf8]{inputenc}
5 \usepackage{a4wide}
6 \usepackage{eurosym}
7 \usepackage{url}
8 %\usepackage{aeguill}
9
10 \usepackage{graphicx}
11 \usepackage{babel}
12 \makeatother
13
14
15 %opening
16 \title{FVN Documentation}
17 \author{William Daniau}
18
19
20 \begin{document}
21
22 \maketitle
23
24 %\begin{abstract}
25
26 %\end{abstract}
27 \tableofcontents
28
29 \section{Whatis fvn,licence,disclaimer etc}
30 \subsection{Whatis fvn}
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
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
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
37 \subsection{Licence}
38 The licence of fvn is free. You can do whatever you want with this code as far as you credit the authors.
39
40 \subsubsection*{Authors}
41 As of the day this manuel is written there's only one author of fvn :\newline
42 William Daniau\newline
43 william.daniau@femto-st.fr\newline
44
45 \subsection{Disclaimer}
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
48 \section{Naming scheme and convention}
49 The naming scheme of the routines is as follow :
50 \begin{verbatim}
51 fvn_x_name()
52 \end{verbatim}
53 where x can be s,d,c or z.
54 \begin{itemize}
55 \item s is for single precision real (real,real*4,real(4),real(kind=4))
56 \item d for double precision real (double precision,real*8,real(8),real(kind=8))
57 \item c for single precision complex (complex,complex*8,complex(4),complex(kind=4))
58 \item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8))
59 \end{itemize}
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
62 \section{Linear algebra}
63 The linear algebra routines of fvn are an interface to lapack, which make it easier to use.
64 \subsection{Matrix inversion}
65 \begin{verbatim}
66 call fvn_x_matinv(d,a,inva,status)
67 \end{verbatim}
68 \begin{itemize}
69 \item d (in) is an integer equal to the matrix rank
70 \item a (in) is a matrix of type x. It will remain untouched.
71 \item inva (out) is a matrix of type x which contain the inverse of a at the end of the routine
72 \item status (out) is an integer equal to zero if something went wrong
73 \end{itemize}
74
75 \subsubsection*{Example}
76 \begin{verbatim}
77 program inv
78 use fvn
79 implicit none
80
81 real(8),dimension(3,3) :: a,inva
82 integer :: status
83
84 call random_number(a)
85 a=a*100
86
87 call fvn_d_matinv(3,a,inva,status)
88 write (*,*) a
89 write (*,*)
90 write (*,*) inva
91 write (*,*)
92 write (*,*) matmul(a,inva)
93 end program
94 \end{verbatim}
95
96
97
98 \subsection{Matrix determinants}
99 \begin{verbatim}
100 det=fvn_x_det(d,a,status)
101 \end{verbatim}
102 \begin{itemize}
103 \item d (in) is an integer equal to the matrix rank
104 \item a (in) is a matrix of type x. It will remain untouched.
105 \item status (out) is an integer equal to zero if something went wrong
106 \end{itemize}
107
108 \subsubsection*{Example}
109 \begin{verbatim}
110 program det
111 use fvn
112 implicit none
113
114 real(8),dimension(3,3) :: a
115 real(8) :: deta
116 integer :: status
117
118 call random_number(a)
119 a=a*100
120
121 deta=fvn_d_det(3,a,status)
122 write (*,*) a
123 write (*,*)
124 write (*,*) "Det = ",deta
125 end program
126
127 \end{verbatim}
128
129
130
131 \subsection{Matrix condition}
132 \begin{verbatim}
133 call fvn_x_matcon(d,a,rcond,status)
134 \end{verbatim}
135 \begin{itemize}
136 \item d (in) is an integer equal to the matrix rank
137 \item a (in) is a matrix of type x. It will remain untouched.
138 \item rcond (out) is a real of same kind as matrix a, it will contain the reciprocal condition number of the matrix
139 \item status (out) is an integer equal to zero if something went wrong
140 \end{itemize}
141
142 The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}
143 \begin{equation}
144 R = \frac{1}{norm(A)*norm(invA)}
145 \label{rconddef}
146 \end{equation}
147
148 The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm}
149 \begin{equation}
150 L1 = max_j ( \sum_i{\mid A(i,j)\mid} )
151 \label{l1norm}
152 \end{equation}
153
154 \subsubsection*{Example}
155 \begin{verbatim}
156 program cond
157 use fvn
158 implicit none
159
160 real(8),dimension(3,3) :: a
161 real(8) :: rcond
162 integer :: status
163
164 call random_number(a)
165 a=a*100
166
167 call fvn_d_matcon(3,a,rcond,status)
168 write (*,*) a
169 write (*,*)
170 write (*,*) "Cond = ",rcond
171 end program
172
173 \end{verbatim}
174
175
176 \subsection{Eigenvalues/Eigenvectors}
177 \begin{verbatim}
178 call fvn_x_matev(d,a,evala,eveca,status)
179 \end{verbatim}
180 \begin{itemize}
181 \item d (in) is an integer equal to the matrix rank
182 \item a (in) is a matrix of type x. It will remain untouched.
183 \item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a
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).
185 \item status (out) is an integer equal to zero if something went wrong
186 \end{itemize}
187
188 \subsubsection*{Example}
189 \begin{verbatim}
190 program eigen
191 use fvn
192 implicit none
193
194 real(8),dimension(3,3) :: a
195 complex(8),dimension(3) :: evala
196 complex(8),dimension(3,3) :: eveca
197 integer :: status,i,j
198
199 call random_number(a)
200 a=a*100
201
202 call fvn_d_matev(3,a,evala,eveca,status)
203 write (*,*) a
204 write (*,*)
205 do i=1,3
206 write(*,*) "Eigenvalue ",i,evala(i)
207 write(*,*) "Associated Eigenvector :"
208 do j=1,3
209 write(*,*) eveca(j,i)
210 end do
211 write(*,*)
212 end do
213
214 end program
215
216 \end{verbatim}
217
218
219 \subsection{Sparse solving}
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
222 The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.
223
224 \begin{verbatim}
225 call fvn_*_sparse_solve(n,nz,T,Ti,Tj,B,x,status) where * is either zl, zi, dl or di
226 \end{verbatim}
227 \begin{itemize}
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)
229 \item n (in) is an integer equal to the matrix rank
230 \item nz (in) is an integer equal to the number of non-zero elements
231 \item T(nz) (in) is a complex/real array containing the non-zero elements
232 \item Ti(nz),Tj(nz) (in) are the indexes of the corresponding element of T in the original matrix.
233 \item B(n) (in) is a complex/real array containing the second member of the equation.
234 \item x(n) (out) is a complex/real array containing the solution
235 \item status (out) is an integer which contain non-zero is something went wrong
236 \end{itemize}
237
238 \subsubsection*{Example}
239 \begin{verbatim}
240 program test_sparse
241
242 use fvn
243 implicit none
244
245 integer(8), parameter :: nz=12
246 integer(8), parameter :: n=5
247 complex(8),dimension(nz) :: A
248 integer(8),dimension(nz) :: Ti,Tj
249 complex(8),dimension(n) :: B,x
250 integer(8) :: status
251
252 A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&
253 (1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /)
254 B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
255 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
256 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
257
258 call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
259 write(*,*) x
260
261 end program
262
263
264 program test_sparse
265
266 use fvn
267 implicit none
268
269 integer(4), parameter :: nz=12
270 integer(4), parameter :: n=5
271 real(8),dimension(nz) :: A
272 integer(4),dimension(nz) :: Ti,Tj
273 real(8),dimension(n) :: B,x
274 integer(4) :: status
275
276 A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
277 B = (/ 8., 45., -3., 3., 19./)
278 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
279 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
280
281 call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
282 write(*,*) x
283
284 end program
285
286
287
288 \end{verbatim}
289
290
291
292 \section{Interpolation}
293
294 \subsection{Quadratic Interpolation}
295 fvn provide function for interpolating values of a tabulated function of 1, 2 or 3 variables, for both single and double precision.
296 \subsubsection{One variable function}
297 \begin{verbatim}
298 value=fvn_x_quad_interpol(x,n,xdata,ydata)
299 \end{verbatim}
300 \begin{itemize}
301 \item x is the real where we want to evaluate the function
302 \item n is the number of tabulated values
303 \item xdata(n) contains the tabulated coordinates
304 \item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))
305 \end{itemize}
306 xdata must be strictly increasingly ordered.
307 x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation
308 \paragraph*{Example}
309 \begin{verbatim}
310 program inter1d
311
312 use fvn
313 implicit none
314
315 integer(kind=4),parameter :: ndata=33
316 integer(kind=4) :: i,nout
317 real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
318 real(kind=8) ::tv
319
320 intrinsic sin
321
322 f(x)=sin(x)
323
324 xdata(1)=0.
325 fdata(1)=f(xdata(1))
326 h=1./32.
327 do i=2,ndata
328 xdata(i)=xdata(i-1)+h
329 fdata(i)=f(xdata(i))
330 end do
331 call random_seed()
332 call random_number(x)
333
334 q=fvn_d_quad_interpol(x,ndata,xdata,fdata)
335
336 tv=f(x)
337 write(*,*) "x ",x
338 write(*,*) "Calculated (real) value :",tv
339 write(*,*) "fvn interpolation :",q
340 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
341
342 end program
343
344 \end{verbatim}
345
346
347 \subsubsection{Two variables function}
348 \begin{verbatim}
349 value=fvn_x_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
350 \end{verbatim}
351 \begin{itemize}
352 \item x,y are the real coordinates where we want to evaluate the function
353 \item nx is the number of tabulated values along x axis
354 \item xdata(nx) contains the tabulated x
355 \item ny is the number of tabulated values along y axis
356 \item ydata(ny) contains the tabulated y
357 \item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))
358 \end{itemize}
359 xdata and ydata must be strictly increasingly ordered.
360 (x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
361
362 \paragraph*{Example}
363
364 \begin{verbatim}
365 program inter2d
366 use fvn
367 implicit none
368
369 integer(kind=4),parameter :: nx=21,ny=42
370 integer(kind=4) :: i,j
371 real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
372 real(kind=8) :: tv
373
374 intrinsic dble,sin
375
376 f(x,y)=sin(x+2.*y)
377 do i=1,nx
378 xdata(i)=dble(i-1)/dble(nx-1)
379 end do
380 do i=1,ny
381 ydata(i)=dble(i-1)/dble(ny-1)
382 end do
383 do i=1,nx
384 do j=1,ny
385 fdata(i,j)=f(xdata(i),ydata(j))
386 end do
387 end do
388 call random_seed()
389 call random_number(x)
390 call random_number(y)
391
392 q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
393 tv=f(x,y)
394
395 write(*,*) "x y",x,y
396 write(*,*) "Calculated (real) value :",tv
397 write(*,*) "fvn interpolation :",q
398 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
399
400 end program
401
402 \end{verbatim}
403
404
405
406 \subsubsection{Three variables function}
407 \begin{verbatim}
408 value=fvn_x_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
409 \end{verbatim}
410 \begin{itemize}
411 \item x,y,z are the real coordinates where we want to evaluate the function
412 \item nx is the number of tabulated values along x axis
413 \item xdata(nx) contains the tabulated x
414 \item ny is the number of tabulated values along y axis
415 \item ydata(ny) contains the tabulated y
416 \item nz is the number of tabulated values along z axis
417 \item zdata(ny) contains the tabulated z
418 \item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k))
419 \end{itemize}
420 xdata, ydata and zdata must be strictly increasingly ordered.
421 (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
423 \paragraph*{Example}
424 \begin{verbatim}
425 program inter3d
426 use fvn
427
428 implicit none
429
430 integer(kind=4),parameter :: nx=21,ny=42,nz=18
431 integer(kind=4) :: i,j,k
432 real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
433 real(kind=8) :: tv
434
435 intrinsic dble,sin
436
437 f(x,y,z)=sin(x+2.*y+3.*z)
438 do i=1,nx
439 xdata(i)=2.*(dble(i-1)/dble(nx-1))
440 end do
441 do i=1,ny
442 ydata(i)=2.*(dble(i-1)/dble(ny-1))
443 end do
444 do i=1,nz
445 zdata(i)=2.*(dble(i-1)/dble(nz-1))
446 end do
447 do i=1,nx
448 do j=1,ny
449 do k=1,nz
450 fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
451 end do
452 end do
453 end do
454 call random_seed()
455 call random_number(x)
456 call random_number(y)
457 call random_number(z)
458
459 q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
460 tv=f(x,y,z)
461
462 write(*,*) "x y z",x,y,z
463 write(*,*) "Calculated (real) value :",tv
464 write(*,*) "fvn interpolation :",q
465 write(*,*) "Relative fvn error :",abs((q-tv)/tv)
466
467 end program
468
469 \end{verbatim}
470
471 \subsubsection{Utility procedure}
472 fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.
473 \begin{verbatim}
474 call fvn_x_find_interval(x,i,xdata,n)
475 \end{verbatim}
476 \begin{itemize}
477 \item x (in) the real value to locate
478 \item i (out) the resulting indice
479 \item xdata(n) (in) increasingly ordered array
480 \item n (in) size of the array
481 \end{itemize}
482 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
484
485
486 \subsection{Akima spline}
487 fvn provides Akima spline interpolation and evaluation for both single and double precision real.
488 \subsubsection{Interpolation}
489 \begin{verbatim}
490 call fvn_x_akima(n,x,y,br,co)
491 \end{verbatim}
492 \begin{itemize}
493 \item n (in) is an integer equal to the number of points
494 \item x(n) (in) ,y(n) (in) are the known couples of coordinates
495 \item br (out) on output contains a copy of x
496 \item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval.
497 \end{itemize}
498
499 \subsubsection{Evaluation}
500 \begin{verbatim}
501 y=fvn_x_spline_eval(x,n,br,co)
502 \end{verbatim}
503 \begin{itemize}
504 \item x (in) is the point where we want to evaluate
505 \item n (in) is the number of known points and br(n) (in), co(4,n) (in) \\
506 are the outputs of fvn\_x\_akima(n,x,y,br,co)
507 \end{itemize}
508
509 \subsubsection{Example}
510 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
512 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
514 \begin{verbatim}
515 program akima
516 use fvn
517 implicit none
518
519 integer :: nbpoints,nppoints,i
520 real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
521 real(8),dimension(:,:),allocatable :: coeff_fvn_d
522 real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
523
524 open(2,file='fvn_akima_double.dat')
525 open(3,file='fvn_akima_breakpoints_double.dat')
526 nbpoints=30
527 allocate(x_d(nbpoints))
528 allocate(y_d(nbpoints))
529 allocate(breakpoints_d(nbpoints))
530 allocate(coeff_fvn_d(4,nbpoints))
531
532 xstep_d=20./dfloat(nbpoints)
533 do i=1,nbpoints
534 x_d(i)=-10.+dfloat(i)*xstep_d
535 y_d(i)=dsin(x_d(i))
536 write(3,44) (x_d(i),y_d(i))
537 end do
538 close(3)
539
540 call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
541
542 nppoints=1000
543 xstep_d=22./dfloat(nppoints)
544 do i=1,nppoints
545 xp_d=-11.+dfloat(i)*xstep_d
546 ty_d=dsin(xp_d)
547 fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
548 write(2,44) (xp_d,ty_d,fvn_y_d)
549 end do
550
551 close(2)
552
553 44 FORMAT(4(1X,1PE22.14))
554
555 end program
556
557 \end{verbatim}
558 Results are plotted on figure \ref{akima}
559
560 \begin{figure}
561 \begin{center}
562 \includegraphics[width=0.9\textwidth]{akima.pdf}
563 % akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720
564 \caption{Akima Spline Interpolation}
565 \label{akima}
566 \end{center}
567
568 \end{figure}
569
570
571
572 \section{Least square polynomial}
573 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
575 \begin{verbatim}
576 call fvn_x_lspoly(np,x,y,deg,coeff,status)
577 \end{verbatim}
578 \begin{itemize}
579 \item np (in) is an integer equal to the number of points
580 \item x(np) (in),y(np) (in) are the known coordinates
581 \item deg (in) is an integer equal to the degree of the desired polynomial, it must be lower than np.
582 \item coeff(deg+1) (out) on output contains the polynomial coefficients
583 \item status (out) is an integer containing 0 if a problem occured.
584 \end{itemize}
585
586 \subsection*{Example}
587 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 :
588 \begin{verbatim}
589 program lsp
590 use fvn
591 implicit none
592
593 integer,parameter :: npoints=13,deg=3
594 integer :: status,i
595 real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
596 real(kind=8) :: coeff(deg+1)
597
598 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. /)
599 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
601 open(2,file='fvn_lsp_double_mesure.dat')
602 open(3,file='fvn_lsp_double_poly.dat')
603
604 do i=1,npoints
605 write(2,44) xm(i),ym(i)
606 end do
607 close(2)
608
609
610 call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)
611
612 xstep=(xm(npoints)-xm(1))/1000.
613 do i=1,1000
614 xc=xm(1)+(i-1)*xstep
615 yc=poly(xc,coeff)
616 write(3,44) xc,yc
617 end do
618 close(3)
619
620 44 FORMAT(4(1X,1PE22.14))
621
622 contains
623 function poly(x,coeff)
624 implicit none
625 real(8) :: x
626 real(8) :: coeff(deg+1)
627 real(8) :: poly
628 integer :: i
629
630 poly=0.
631
632 do i=1,deg+1
633 poly=poly+coeff(i)*x**(i-1)
634 end do
635
636 end function
637 end program
638 \end{verbatim}
639 The results are plotted on figure \ref{lsp} .
640
641 \begin{figure}
642 \begin{center}
643 \includegraphics[width=0.9\textwidth]{lsp.pdf}
644 \caption{Least Square Polynomial}
645 \label{lsp}
646 \end{center}
647 \end{figure}
648
649
650
651 \section{Zero finding}
652 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
654 \begin{verbatim}
655 call fvn_z_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
656 \end{verbatim}
657 \begin{itemize}
658 \item f (in) is the complex function (kind=8) for which we search zeros
659 \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.
660 \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.
661 \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.
662 \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.
663 \item n (in) is an integer equal to the number of new roots to be found.
664 \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.
665 \item itmax (in) is an integer equal to the maximum allowable number of iterations per root.
666 \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
667 \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.
668 \end{itemize}
669 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
671 \subsection*{Example}
672 Example to find the ten roots of $x^{10}-1$
673 \begin{verbatim}
674 program muller
675 use fvn
676 implicit none
677
678 integer :: i,info
679 complex(8),dimension(10) :: roots
680 integer,dimension(10) :: infer
681 complex(8), external :: f
682
683 call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
684
685 write(*,*) "Error code :",info
686 do i=1,10
No preview for this file type