PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
num_ops.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module num_ops
5 #include <PB3D_macros.h>
6  use str_utilities
7  use messages
8  use num_vars, only: dp, iu, max_str_ln, pi
9  use num_utilities
10  use output_ops
11 
12  implicit none
13  private
15 #if ldebug
16  public debug_calc_zero
17 #endif
18 
19  ! global variables
20 #if ldebug
21 
22  logical :: debug_calc_zero = .false.
23 #endif
24 
25  ! interfaces
26 
35  interface calc_zero_hh
37  module procedure calc_zero_hh_0d
39  module procedure calc_zero_hh_3d
40  end interface
41 
42 contains
48  function calc_zero_hh_0d(zero,fun,ord,guess,max_nr_backtracks,output) &
49  &result(err_msg)
51 
52  ! input / output
53  real(dp), intent(inout) :: zero
54  interface
55  function fun(x,ord) ! the function
56  use num_vars, only: dp
57  real(dp), intent(in) :: x
58  integer, intent(in) :: ord
59  real(dp) :: fun
60  end function fun
61  end interface
62  integer, intent(in) :: ord
63  real(dp), intent(in) :: guess
64  integer, intent(in), optional :: max_nr_backtracks
65  logical, intent(in), optional :: output
66  character(len=max_str_ln) :: err_msg
67 
68  ! local variables
69  integer :: id, jd, kd ! counters
70  integer :: max_nr_backtracks_loc ! local max_nr_backtracks
71  logical :: output_loc ! local output
72  logical :: relaxed_enough ! whether step was relaxed enough to get a smaller new value
73  real(dp) :: zero_new ! new zero calculated
74  real(dp) :: fun_new ! function value at new zero
75  real(dp) :: corr ! correction
76  real(dp), allocatable :: fun_vals(:) ! function values
77 #if ldebug
78  real(dp), allocatable :: corrs(:) ! corrections for all steps
79  real(dp), allocatable :: values(:) ! values for all steps
80 #endif
81 
82  ! initialize error message
83  err_msg = ''
84 
85  ! set up local output
86  output_loc = .false.
87  if (present(output)) output_loc = output
88 
89  ! set up zero
90  zero = guess
91 
92  ! set up local max_nr_backtracks
93  max_nr_backtracks_loc = max_nr_backtracks_hh
94  if (present(max_nr_backtracks)) &
95  &max_nr_backtracks_loc = max_nr_backtracks
96 
97  ! initialize function values
98  allocate(fun_vals(0:ord))
99 
100 #if ldebug
101  if (debug_calc_zero) then
102  ! set up corrs
103  allocate(corrs(max_it_zero))
104  allocate(values(max_it_zero))
105  end if
106 #endif
107 
108  ! loop over iterations
109  hh: do jd = 1,max_it_zero
110  ! calculate function values
111  do kd = 0,ord
112  fun_vals(kd) = fun(zero,kd)
113  end do
114 
115  ! correction to theta_HH
116  select case (ord)
117  case (1) ! Newton-Rhapson
118  corr = -fun_vals(0)/fun_vals(1)
119  case (2) ! Halley
120  corr = -fun_vals(0)*fun_vals(1)/&
121  &(fun_vals(1)**2 - 0.5_dp*&
122  &fun_vals(0)*fun_vals(2))
123  case (3) ! 3rd order
124  corr = -(6*fun_vals(0)*fun_vals(1)**2-&
125  &3*fun_vals(0)**2*fun_vals(2))/&
126  &(6*fun_vals(1)**3 - 6*fun_vals(0)*&
127  &fun_vals(1)*fun_vals(2)+&
128  &fun_vals(0)**2*fun_vals(3))
129  end select
130 #if ldebug
131  if (debug_calc_zero) then
132  corrs(jd) = corr
133  values(jd) = zero
134  end if
135 #endif
136  relaxed_enough = .false.
137  do id = 1,max_nr_backtracks_loc
138  zero_new = zero + corr ! propose a new zero
139 
140  fun_new = fun(zero_new,0) ! calculate function value there
141 
142  if (abs(fun_new)-abs(fun_vals(0)).le.0._dp) then ! error has decreased
143  relaxed_enough = .true.
144  zero = zero_new
145  exit
146  end if
147  corr = corr*0.5_dp
148  end do
149 
150  ! check for relaxation
151  if (.not.relaxed_enough) then
152  err_msg = trim(i2str(max_nr_backtracks_hh))//' backtracks were &
153  &not sufficient to get a better guess for the zero'
154  zero = 0.0_dp
155  else
156  ! output
157  if (output_loc) call writo(trim(i2str(jd))//' / '&
158  &//trim(i2str(max_it_zero))//&
159  &' - relative error: '//&
160  &trim(r2str(abs(corr)))//' (tol: '//&
161  &trim(r2str(tol_zero))//')')
162 
163  ! check for convergence
164  if (abs(corr).lt.tol_zero) then
165 #if ldebug
166  if (debug_calc_zero) &
167  &call plot_evolution(corrs(1:jd),values(1:jd))
168 #endif
169  return
170  else if (jd .eq. max_it_zero) then
171 #if ldebug
172  if (debug_calc_zero) then
173  call plot_evolution(corrs,values)
174  deallocate(corrs)
175  deallocate(values)
176  end if
177 #endif
178  err_msg = 'Not converged after '//trim(i2str(jd))//&
179  &' iterations, with residual '//&
180  &trim(r2strt(corr))//' and final value '//&
181  &trim(r2strt(zero))
182  zero = 0.0_dp
183  end if
184  end if
185  end do hh
186 #if ldebug
187  contains
188  ! plots corrections
190  subroutine plot_evolution(corrs,values)
191  ! input / output
192  real(dp), intent(in) :: corrs(:) ! corrections
193  real(dp), intent(in) :: values(:) ! values
194 
195  ! local variables
196  real(dp) :: min_x, max_x ! min. and max. x for which to plot the function
197  real(dp) :: extra_x = 1._dp ! how much bigger to take the plot interval than just min and max
198  real(dp) :: delta_x ! interval between zero and guess
199  integer :: n_x = 200 ! how many x values to plot
200  real(dp), allocatable :: x_out(:,:) ! output x of plot
201  real(dp), allocatable :: y_out(:,:) ! output y of plot
202 
203  ! set up min and max x
204  min_x = min(zero,guess)
205  max_x = max(zero,guess)
206  delta_x = max_x-min_x
207  min_x = min_x - delta_x*extra_x
208  max_x = max_x + delta_x*extra_x
209 
210  ! set up output of plot
211  allocate(x_out(n_x,2))
212  allocate(y_out(n_x,2))
213  x_out(:,1) = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
214  x_out(:,2) = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
215  y_out(:,1) = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),0),&
216  &jd=1,n_x)]
217  y_out(:,2) = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),1),&
218  &jd=1,n_x)]
219 
220  ! plot function
221  call writo('Initial guess: '//trim(r2str(guess)))
222  call writo('Resulting zero: '//trim(r2str(zero)))
223  call print_ex_2d(['function ','derivative'],'',y_out,x=x_out)
224 
225  ! user output
226  call writo('Last correction '//trim(r2strt(corrs(size(corrs))))&
227  &//' with tolerance '//trim(r2strt(tol_zero)))
228 
229  ! plot corrections and values
230  call print_ex_2d('corrs','',corrs)
231  call print_ex_2d('values','',values)
232  end subroutine plot_evolution
233 #endif
234  end function calc_zero_hh_0d
241  function calc_zero_hh_3d(dims,zero,fun,ord,guess,max_nr_backtracks,output) &
242  &result(err_msg)
244 
245  ! input / output
246  integer, intent(in) :: dims(3)
247  real(dp), intent(inout) :: zero(dims(1),dims(2),dims(3))
248  interface
249  function fun(dims,x,ord) ! the function
250  use num_vars, only: dp
251  integer, intent(in) :: dims(3)
252  real(dp), intent(in) :: x(dims(1),dims(2),dims(3))
253  integer, intent(in) :: ord
254  real(dp) :: fun(dims(1),dims(2),dims(3))
255  end function fun
256  end interface
257  integer, intent(in) :: ord
258  real(dp), intent(in) :: guess(dims(1),dims(2),dims(3))
259  integer, intent(in), optional :: max_nr_backtracks
260  logical, intent(in), optional :: output
261  character(len=max_str_ln) :: err_msg
262 
263  ! local variables
264  integer :: id, jd, kd ! counters
265  integer :: max_nr_backtracks_loc ! local max_nr_backtracks
266  integer :: mc_ind(3) ! index of maximum correction
267  logical :: output_loc ! local output
268  logical :: relaxed_enough ! relaxed enough
269  real(dp) :: zero_new(dims(1),dims(2),dims(3)) ! new zeros calculated
270  real(dp) :: fun_new(dims(1),dims(2),dims(3)) ! function values at new zero
271  real(dp) :: corr(dims(1),dims(2),dims(3)) ! correction
272  real(dp), allocatable :: fun_vals(:,:,:,:) ! function values
273 #if ldebug
274  real(dp), allocatable :: corrs(:,:,:,:) ! corrections for all steps
275  real(dp), allocatable :: values(:,:,:,:) ! values for all steps
276 #endif
277 
278  ! initialize error message
279  err_msg = ''
280 
281  ! test
282  if (ord.lt.1 .or. ord.gt.3) then
283  zero = 0.0_dp
284  err_msg = 'only orders 1 (Newton-Rhapson), 2 (Halley) and 3 &
285  &implemented, not '//trim(i2str(ord))
286  return
287  end if
288 
289  ! set up local output
290  output_loc = .false.
291  if (present(output)) output_loc = output
292 
293  ! set up zero
294  zero = guess
295 
296  ! set up local max_nr_backtracks
297  max_nr_backtracks_loc = max_nr_backtracks_hh
298  if (present(max_nr_backtracks)) &
299  &max_nr_backtracks_loc = max_nr_backtracks
300 
301  ! initialize function values
302  allocate(fun_vals(dims(1),dims(2),dims(3),0:ord))
303 
304 #if ldebug
305  if (debug_calc_zero) then
306  ! set up corrs
307  allocate(corrs(dims(1),dims(2),dims(3),max_it_zero))
308  allocate(values(dims(1),dims(2),dims(3),max_it_zero))
309  end if
310 #endif
311 
312  ! loop over iterations
313  hh: do jd = 1,max_it_zero
314  ! calculate function values
315  do kd = 0,ord
316  fun_vals(:,:,:,kd) = fun(dims,zero,kd)
317  end do
318 
319  ! correction to theta_HH
320  select case (ord)
321  case (1) ! Newton-Rhapson
322  corr = -fun_vals(:,:,:,0)/fun_vals(:,:,:,1)
323  case (2) ! Halley
324  corr = -fun_vals(:,:,:,0)*fun_vals(:,:,:,1)/&
325  &(fun_vals(:,:,:,1)**2 - 0.5_dp*&
326  &fun_vals(:,:,:,0)*fun_vals(:,:,:,2))
327  case (3) ! 3rd order
328  corr = -(6*fun_vals(:,:,:,0)*fun_vals(:,:,:,1)**2-&
329  &3*fun_vals(:,:,:,0)**2*fun_vals(:,:,:,2))/&
330  &(6*fun_vals(:,:,:,1)**3 - 6*fun_vals(:,:,:,0)*&
331  &fun_vals(:,:,:,1)*fun_vals(:,:,:,2)+&
332  &fun_vals(:,:,:,0)**2*fun_vals(:,:,:,3))
333  end select
334 #if ldebug
335  if (debug_calc_zero) then
336  corrs(:,:,:,jd) = corr
337  values(:,:,:,jd) = zero
338  end if
339 #endif
340  relaxed_enough = .false.
341  do id = 1,max_nr_backtracks_loc
342  zero_new = zero ! otherwise some values might stay unitialized
343  where (abs(corr).gt.tol_zero) zero_new = zero + corr ! propose a new zero
344 
345  fun_new = fun(dims,zero_new,0) ! calculate function value there
346 
347 #if ldebug
348  if (debug_calc_zero) &
349  &call writo('backtrack '//trim(i2str(id))//' of max. '//&
350  &trim(i2str(max_nr_backtracks_loc))//' - criterion: '//&
351  &trim(r2str(maxval(abs(fun_new)-abs(fun_vals(:,:,:,0)))))//&
352  &' <= 0?')
353 #endif
354  if (maxval(abs(fun_new)-abs(fun_vals(:,:,:,0))).le.0._dp) then ! error has decreased
355  relaxed_enough = .true.
356  zero = zero_new
357  exit
358  end if
359  corr = corr*0.5_dp
360  end do
361 
362  ! check for relaxation
363  if (.not.relaxed_enough) then
364  err_msg = trim(i2str(max_nr_backtracks_hh))//' backtracks were &
365  &not sufficient to get a better guess for the zero'
366  zero = 0.0_dp
367  exit
368  else
369  ! output
370  if (output_loc) call writo(trim(i2str(jd))//' / '&
371  &//trim(i2str(max_it_zero))//&
372  &' - maximum relative error: '//&
373  &trim(r2str(maxval(abs(corr))))//' (tol: '//&
374  &trim(r2str(tol_zero))//')')
375 
376  ! check for convergence
377  if (maxval(abs(corr)).lt.tol_zero) then
378 #if ldebug
379  if (debug_calc_zero) call &
380  &plot_evolution(corrs(:,:,:,1:jd),values(:,:,:,1:jd))
381 #endif
382  return
383  else if (jd .eq. max_it_zero) then
384 #if ldebug
385  if (debug_calc_zero) then
386  call plot_evolution(corrs,values)
387  deallocate(corrs)
388  deallocate(values)
389  end if
390 #endif
391  mc_ind = maxloc(abs(corr))
392  err_msg = 'Not converged after '//trim(i2str(jd))//&
393  &' iterations, with maximum residual '//&
394  &trim(r2strt(corr(mc_ind(1),mc_ind(2),mc_ind(3))))&
395  &//' and final value '//trim(r2strt(&
396  &zero(mc_ind(1),mc_ind(2),mc_ind(3))))
397  zero = 0.0_dp
398  end if
399  end if
400  end do hh
401 #if ldebug
402  contains
403  ! plots corrections
405  subroutine plot_evolution(corrs,values)
406  ! input / output
407  real(dp), intent(in) :: corrs(:,:,:,:) ! corrections
408  real(dp), intent(in) :: values(:,:,:,:) ! values
409 
410  ! local variables
411  integer :: n_corrs ! number of corrections
412  character(len=max_str_ln) :: var_name(1) ! names of variable
413  character(len=max_str_ln) :: file_name ! name of file
414 
415  ! initialize
416  n_corrs = size(corrs,4)
417  mc_ind = maxloc(abs(corrs(:,:,:,n_corrs)))
418 
419  ! user output
420  call writo('Last maximum correction '//&
421  &trim(r2strt(corrs(mc_ind(1),mc_ind(2),mc_ind(3),n_corrs)))&
422  &//' and final value '//trim(r2strt(&
423  &values(mc_ind(1),mc_ind(2),mc_ind(3),n_corrs)))//&
424  &' with tolerance '//trim(r2strt(tol_zero)))
425 
426  ! plot corrs
427  file_name = 'corrs'
428  var_name = 'var'
429 
430  call plot_hdf5(var_name,file_name,corrs,col=2,&
431  &descr='corrections')
432 
433  ! plot values
434  file_name = 'values'
435  var_name = 'var'
436 
437  call plot_hdf5(var_name,file_name,values,col=2,&
438  &descr='values')
439  end subroutine plot_evolution
440 #endif
441  end function calc_zero_hh_3d
442 
461  function calc_zero_zhang(zero,fun,x_int_in) result(err_msg)
463 
464  ! input / output
465  real(dp), intent(inout) :: zero
466  interface
467  function fun(x) ! the function
468  use num_vars, only: dp
469  real(dp), intent(in) :: x
470  real(dp) :: fun
471  end function fun
472  end interface
473  real(dp), intent(in) :: x_int_in(2)
474  character(len=max_str_ln) :: err_msg
475 
476  ! local variables
477  logical :: converged ! whether converged
478  integer :: jd ! counter
479  integer :: id_loc ! local id (either 1 or 2)
480  real(dp) :: x_int(2) ! points of current interval [a,b]
481  real(dp) :: x_mid ! x at midpoint [c]
482  real(dp) :: x_rule ! x found with rule [s]
483  real(dp) :: x_temp ! temporary x
484  real(dp) :: fun_int(2) ! function at x_int
485  real(dp) :: fun_mid ! function at x_mid
486  real(dp) :: fun_rule ! function at x_rule
487 #if ldebug
488  real(dp), allocatable :: x_ints(:,:) ! bounds for all steps
489 #endif
490 
491  ! initialize error message and zero
492  err_msg = ''
493  zero = 0.0_dp
494 
495  ! initialize from input
496  x_int(1) = minval(x_int_in)
497  x_int(2) = maxval(x_int_in)
498  fun_int = [fun(x_int(1)),fun(x_int(2))]
499 
500  ! test whether already zero
501  if (abs(fun_int(1)).lt.tol_zero) then
502  zero = x_int(1)
503  return
504  end if
505  if (abs(fun_int(2)).lt.tol_zero) then
506  zero = x_int(2)
507  return
508  end if
509 
510  ! test whether different sign
511  if (product(fun_int).gt.0) then
512  err_msg = 'Function values on starting interval have same sign'
513  return
514  end if
515 
516  ! intialize converged
517  converged = .false.
518 
519 #if ldebug
520  if (debug_calc_zero) then
521  ! set up x_ints
522  allocate(x_ints(max_it_zero,2))
523  end if
524 #endif
525 
526  ! iterations
527  do jd = 1,max_it_zero
528  ! calculate midpoint and function value
529  x_mid = 0.5_dp*sum(x_int)
530  fun_mid = fun(x_mid)
531 
532  if (abs(fun_int(1)-fun_mid).gt.tol_zero .and. &
533  &abs(fun_int(2)-fun_mid).gt.tol_zero) then
534  ! three different function values: use inverse quadratic
535  ! interpolation
536  x_rule = fun_mid*fun_int(2)*x_int(1)/&
537  &((fun_int(1)-fun_mid)*(fun_int(1)-fun_int(2))) + &
538  &fun_int(1)*fun_int(2)*x_mid/&
539  &((fun_mid-fun_int(1))*(fun_mid-fun_int(2))) + &
540  &fun_int(1)*fun_mid*x_int(2)/&
541  &((fun_int(2)-fun_mid)*(fun_int(2)-fun_mid))
542 
543  if (x_int(1).lt.x_rule .and. x_rule.lt.x_int(2)) then
544  ! found x within interval
545  else
546  ! use bisection instead
547  if (fun_int(1)*fun_mid.lt.0) then ! [a,c] contains sign change
548  id_loc = 1
549  else if (fun_int(2)*fun_mid.lt.0) then ! [b,c] contains sign change
550  id_loc = 2
551  else
552  err_msg = 'Something went wrong...'
553  converged = .true.
554  end if
555  if (.not.converged) x_rule = 0.5_dp*(x_int(id_loc)+x_mid)
556  end if
557  else
558  ! two function values overlap so use secant rule
559  if (fun_int(1)*fun_mid.lt.0) then ! [a,c] contains sign change
560  id_loc = 1
561  else if (fun_int(2)*fun_mid.lt.0) then ! [b,c] contains sign change
562  id_loc = 2
563  else
564  err_msg = 'Something went wrong...'
565  converged = .true.
566  end if
567  if (.not.converged) x_rule = x_int(id_loc) - &
568  &fun_int(id_loc)*(x_int(id_loc)-x_mid)/&
569  &(fun_int(id_loc)-fun_mid)
570  end if
571 
572  ! calculate fun(s)
573  fun_rule = fun(x_rule)
574 
575  ! make sure x_mid < x_rule
576  if (x_mid.gt.x_rule) then
577  x_temp = x_rule
578  x_rule = x_mid
579  x_mid = x_temp
580  end if
581 
582  ! check whether [c,s] contains root
583  if (fun_mid*fun_rule.lt.0) then
584  x_int(1) = x_mid
585  x_int(2) = x_rule
586  else
587  if (fun_int(1)*fun_mid.lt.0) then
588  x_int(2) = x_mid
589  else
590  x_int(1) = x_rule
591  end if
592  end if
593 
594 #if ldebug
595  if (debug_calc_zero) then
596  x_ints(jd,:) = x_int
597  end if
598 #endif
599 
600  ! check for convergence
601  fun_int = [fun(x_int(1)),fun(x_int(2))]
602  if (abs(fun_int(1)).lt.tol_zero) then
603  zero = x_int(1)
604  converged = .true.
605  end if
606  if (abs(fun_int(2)).lt.tol_zero) then
607  zero = x_int(2)
608  converged = .true.
609  end if
610  if ((x_int(2)-x_int(1)).lt.tol_zero) then
611  zero = 0.5_dp*sum(x_int)
612  converged = .true.
613  end if
614 
615  if (converged) then
616 #if ldebug
617  if (debug_calc_zero) call plot_evolution(x_ints(1:jd,:))
618 #endif
619  return
620  end if
621  end do
622 #if ldebug
623  contains
624  ! plots corrections
626  subroutine plot_evolution(x_ints)
627  ! input / output
628  real(dp), intent(in) :: x_ints(:,:) ! x_intervals
629 
630  ! local variables
631  real(dp) :: min_x, max_x ! min. and max. x for which to plot the function
632  integer :: n_x = 200 ! how many x values to plot
633  integer :: n_ints ! nr. of intervals
634  real(dp), allocatable :: x_out(:) ! output x of plot
635  real(dp), allocatable :: y_out(:) ! output y of plot
636 
637  ! set up nr. of intervals
638  n_ints = size(x_ints,1)
639 
640  ! set up min and max x
641  min_x = x_int_in(1)
642  max_x = x_int_in(2)
643 
644  ! set up output of plot
645  allocate(x_out(n_x))
646  allocate(y_out(n_x))
647  x_out = [(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1),jd=1,n_x)]
648  y_out = [(fun(min_x+(max_x-min_x)*(jd-1._dp)/(n_x-1)),jd=1,n_x)]
649 
650  ! plot function
651  call writo('Initial interval: ['//trim(r2str(x_int_in(1)))//&
652  &'..'//trim(r2str(x_int_in(2)))//']')
653  call writo('Resulting zero: '//trim(r2str(zero)))
654  call print_ex_2d('function','',y_out,x=x_out)
655 
656  ! user output
657  call writo('Final interval ['//trim(r2str(x_ints(n_ints,1)))//&
658  &'..'//trim(r2str(x_ints(n_ints,2)))//'] with tolerance '//&
659  &trim(r2strt(tol_zero)))
660 
661  ! plot intervals
662  call print_ex_2d(['x_intervals'],'',x_ints)
663  end subroutine plot_evolution
664 #endif
665  end function calc_zero_zhang
666 end module num_ops
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
num_vars::max_it_zero
integer, public max_it_zero
maximum number of iterations to find zeros
Definition: num_vars.f90:131
num_ops
Numerical operations.
Definition: num_ops.f90:4
str_utilities::i2str
elemental character(len=max_str_ln) function, public i2str(k)
Convert an integer to string.
Definition: str_utilities.f90:18
num_vars::iu
complex(dp), parameter, public iu
complex unit
Definition: num_vars.f90:85
num_ops::debug_calc_zero
logical, public debug_calc_zero
plot debug information for calc_zero
Definition: num_ops.f90:22
output_ops::print_ex_2d
Print 2-D output on a file.
Definition: output_ops.f90:47
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
str_utilities::r2strt
elemental character(len=max_str_ln) function, public r2strt(k)
Convert a real (double) to string.
Definition: str_utilities.f90:54
output_ops::plot_hdf5
Prints variables vars with names var_names in an HDF5 file with name c file_name and accompanying XDM...
Definition: output_ops.f90:137
num_vars::max_nr_backtracks_hh
integer, public max_nr_backtracks_hh
maximum number of backtracks for Householder, relax. factors
Definition: num_vars.f90:132
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
num_ops::calc_zero_hh
Finds the zero of a function using Householder iteration.
Definition: num_ops.f90:35
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
str_utilities::r2str
elemental character(len=max_str_ln) function, public r2str(k)
Convert a real (double) to string.
Definition: str_utilities.f90:42
messages::writo
subroutine, public writo(input_str, persistent, error, warning, alert)
Write output to file identified by output_i.
Definition: messages.f90:275
messages
Numerical utilities related to giving output.
Definition: messages.f90:4
num_vars::pi
real(dp), parameter, public pi
Definition: num_vars.f90:83
num_ops::calc_zero_zhang
character(len=max_str_ln) function, public calc_zero_zhang(zero, fun, x_int_in)
Finds the zero of a function using Zhang's method, which is simpler than Brent's method.
Definition: num_ops.f90:462
output_ops
Operations concerning giving output, on the screen as well as in output files.
Definition: output_ops.f90:5