PB3D  [2.45]
Ideal linear high-n MHD stability in 3-D
num_utilities.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------!
3 !------------------------------------------------------------------------------!
4 module num_utilities
5 #include <PB3D_macros.h>
6 #include <wrappers.h>
7  use str_utilities
8  use messages
9  use num_vars, only: dp, iu, max_str_ln, pi
10 
11  implicit none
12  private
13  public calc_ext_var, calc_det, calc_int, add_arr_mult, c, &
14  &check_deriv, calc_inv, calc_mult, calc_aux_utilities, derivs, &
15  &con2dis, dis2con, round_with_tol, conv_mat, is_sym, con, &
16  &calc_coeff_fin_diff, fac, d, m, f, bubble_sort, gcd, lcm, &
17  &order_per_fun, shift_f, spline, solve_vand
18 #if ldebug
19  public debug_con2dis_reg, debug_calc_coeff_fin_diff, calc_d2_smooth
20 #endif
21 
22  ! global variables
23  integer, allocatable :: d(:,:,:)
24  integer, allocatable :: m(:,:)
25  integer, allocatable :: f(:,:)
26 #if ldebug
27  logical :: debug_con2dis_reg = .false.
28  logical :: debug_calc_coeff_fin_diff = .false.
29 #endif
30 
31  ! interfaces
32 
39  interface add_arr_mult
41  module procedure add_arr_mult_3_3
43  module procedure add_arr_mult_3_1
45  module procedure add_arr_mult_1_1
46  end interface
47 
63  interface calc_det
65  module procedure calc_det_0d
67  module procedure calc_det_3d
68  end interface
69 
81  interface calc_inv
83  module procedure calc_inv_0d
85  module procedure calc_inv_3d
86  end interface
87 
95  interface calc_mult
97  module procedure calc_mult_0d_real
99  module procedure calc_mult_3d_real
101  module procedure calc_mult_3d_complex
102  end interface
103 
123  interface conv_mat
125  module procedure conv_mat_3d
127  module procedure conv_mat_3d_complex
129  module procedure conv_mat_0d
131  module procedure conv_mat_0d_complex
132  end interface
133 
160  interface calc_int
162  module procedure calc_int_eqd
164  module procedure calc_int_reg
165  end interface
166 
169  interface round_with_tol
171  module procedure round_with_tol_ind
173  module procedure round_with_tol_arr
174  end interface
175 
188  interface con2dis
190  module procedure con2dis_eqd
192  module procedure con2dis_reg
193  end interface
194 
206  interface dis2con
208  module procedure dis2con_eqd
210  module procedure dis2con_reg
211  end interface
212 
221  interface con
223  module procedure con_3d
225  module procedure con_2d
227  module procedure con_1d
229  module procedure con_0d
230  end interface
231 
237  interface bubble_sort
239  module procedure bubble_sort_int
241  module procedure bubble_sort_real
242  end interface
243 
248  interface order_per_fun
250  module procedure order_per_fun_1
252  module procedure order_per_fun_2
253  end interface
254 
276  interface spline
278  module procedure spline_real
280  module procedure spline_complex
281  end interface
282 
283 contains
285  integer function calc_int_reg(var,x,var_int) result(ierr)
286  character(*), parameter :: rout_name = 'calc_int_reg'
287 
288  ! input / output
289  real(dp), intent(inout) :: var_int(:)
290  real(dp), intent(in) :: var(:)
291  real(dp), intent(in) :: x(:)
292 
293  ! local variables
294  integer :: n_max
295  integer :: kd
296  character(len=max_str_ln) :: err_msg ! error message
297 
298  ! initialize ierr
299  ierr = 0
300 
301  ! set n_max
302  n_max = size(var)
303 
304  ! tests
305  if (size(x).ne.n_max) then
306  err_msg = 'The arrays of points and values are not of the same &
307  &length'
308  ierr = 1
309  chckerr(err_msg)
310  else if (n_max.lt.2) then
311  err_msg = 'Asking to integrate with '//trim(i2str(n_max))&
312  &//' points. Need at least 2'
313  ierr = 1
314  chckerr(err_msg)
315  end if
316 
317  ! calculate integral for all points
318  var_int = 0.0_dp
319 
320  do kd = 2, n_max
321  var_int(kd) = var_int(kd-1) + &
322  &(var(kd)+var(kd-1))/2 * (x(kd)-x(kd-1))
323  end do
324  end function calc_int_reg
326  integer function calc_int_eqd(var,step_size,var_int) result(ierr)
327  character(*), parameter :: rout_name = 'calc_int_eqd'
328 
329  ! input / output
330  real(dp), intent(inout) :: var_int(:)
331  real(dp), intent(in) :: var(:)
332  real(dp), intent(in) :: step_size
333 
334  ! local variables
335  integer :: n_max
336  integer :: kd
337  character(len=max_str_ln) :: err_msg ! error message
338 
339  ! initialize ierr
340  ierr = 0
341 
342  ! set n_max
343  n_max = size(var)
344 
345  ! tests
346  if (n_max.lt.2) then
347  err_msg = 'Asking to integrate with '//trim(i2str(n_max))&
348  &//' points. Need at least 2'
349  ierr = 1
350  chckerr(err_msg)
351  end if
352 
353  ! calculate integral for all points
354  var_int = 0.0_dp
355 
356  do kd = 2, n_max
357  var_int(kd) = var_int(kd-1) + (var(kd)+var(kd-1))/2 * step_size
358  end do
359  end function calc_int_eqd
360 
362  integer function add_arr_mult_3_3(arr_1,arr_2,arr_3,deriv) result(ierr)
363  character(*), parameter :: rout_name = 'add_arr_mult_3_3'
364 
365  ! input / output
366  real(dp), intent(in) :: arr_1(1:,1:,1:,0:,0:,0:)
367  real(dp), intent(in) :: arr_2(1:,1:,1:,0:,0:,0:)
368  real(dp), intent(out) :: arr_3(1:,1:,1:)
369  integer, intent(in) :: deriv(3)
370 
371  ! local variables
372  real(dp) :: bin_fac(3) ! binomial factor for norm., pol. and tor. sum
373  integer :: r, m, n ! current degree of norm., pol., tor. derivatives
374  integer :: kd ! normal counter
375  character(len=max_str_ln) :: err_msg ! error message
376 
377  ! initialize ierr
378  ierr = 0
379 
380  ! tests
381  if (size(arr_1,1).ne.size(arr_2,1) .or. size(arr_1,2).ne.size(arr_2,2) &
382  &.or. size(arr_1,3).ne.size(arr_2,3)) then
383  err_msg = 'Arrays 1 and 2 need to have the same size'
384  ierr = 1
385  chckerr(err_msg)
386  end if
387  if (size(arr_1,1).ne.size(arr_3,1) .or. size(arr_1,2).ne.size(arr_3,2) &
388  &.or. size(arr_1,3).ne.size(arr_3,3)) then
389  err_msg = 'Arrays 1 and 2 need to have the same size as the &
390  &resulting array 3'
391  ierr = 1
392  chckerr(err_msg)
393  end if
394  do kd = 1,3
395  if (size(arr_1,3+kd).lt.deriv(kd)+1 .or. &
396  &size(arr_2,3+kd).lt.deriv(kd)+1) then
397  err_msg = 'Arrays 1 and 2 do not provide deriv. orders for a &
398  &derivative of order ['//trim(i2str(deriv(1)))//','//&
399  &trim(i2str(deriv(2)))//','//trim(i2str(deriv(3)))//'&
400  &] in coordinate '//trim(i2str(kd))
401  ierr = 1
402  chckerr(err_msg)
403  end if
404  end do
405 
406  ! distribute normal and angular derivatives using binomial theorem
407  ! normal derivatives
408  do r = 0,deriv(1)
409  if (r.eq.0) then ! first term in sum
410  bin_fac(1) = 1.0_dp
411  else
412  bin_fac(1) = bin_fac(1)*(deriv(1)-(r-1))/r
413  end if
414  ! poloidal derivatives
415  do m = 0,deriv(2)
416  if (m.eq.0) then ! first term in sum
417  bin_fac(2) = 1.0_dp
418  else
419  bin_fac(2) = bin_fac(2)*(deriv(2)-(m-1))/m
420  end if
421  ! toroidal derivatives
422  do n = 0,deriv(3)
423  if (n.eq.0) then ! first term in sum
424  bin_fac(3) = 1.0_dp
425  else
426  bin_fac(3) = bin_fac(3)*(deriv(3)-(n-1))/n
427  end if
428 
429  ! current term in the tripple summation
430  arr_3 = arr_3 + &
431  &bin_fac(1)*bin_fac(2)*bin_fac(3) &
432  &* arr_1(:,:,:,r,m,n)&
433  &* arr_2(:,:,:,deriv(1)-r,deriv(2)-m,deriv(3)-n)
434  end do
435  end do
436  end do
437  end function add_arr_mult_3_3
440  integer function add_arr_mult_3_1(arr_1,arr_2,arr_3,deriv) result(ierr)
441  character(*), parameter :: rout_name = 'add_arr_mult_3_1'
442 
443  ! input / output
444  real(dp), intent(in) :: arr_1(1:,1:,1:,0:,0:,0:)
445  real(dp), intent(in) :: arr_2(1:,0:)
446  real(dp), intent(out) :: arr_3(1:,1:,1:)
447  integer, intent(in) :: deriv(3)
448 
449  ! local variables
450  real(dp) :: bin_fac ! binomial factor for norm. sum
451  integer :: r ! current degree of norm. derivative
452  integer :: kd ! normal counter
453  character(len=max_str_ln) :: err_msg ! error message
454 
455  ! initialize ierr
456  ierr = 0
457 
458  ! tests
459  if (size(arr_1,3).ne.size(arr_2,1)) then
460  err_msg = 'Arrays 1 and 2 need to have the same size in the flux &
461  &coord.'
462  ierr = 1
463  chckerr(err_msg)
464  end if
465  if (size(arr_1,1).ne.size(arr_3,1) .or. size(arr_1,2).ne.size(arr_3,2) &
466  &.or. size(arr_1,3).ne.size(arr_3,3)) then
467  err_msg = 'Array 1 needs to have the same size as the resulting &
468  &array 3'
469  ierr = 1
470  chckerr(err_msg)
471  end if
472 
473  ! distribute normal derivatives using binomial theorem, angular
474  ! derivatives only operate on the first array
475  ! normal derivatives
476  do r = 0,deriv(1)
477  if (r.eq.0) then ! first term in sum
478  bin_fac = 1.0_dp
479  else
480  bin_fac = bin_fac*(deriv(1)-(r-1))/r
481  end if
482 
483  ! current term in the tripple summation
484  do kd = 1, size(arr_1,3)
485  arr_3(:,:,kd) = arr_3(:,:,kd) + bin_fac * &
486  &arr_1(:,:,kd,r,deriv(2),deriv(3))* arr_2(kd,deriv(1)-r)
487  end do
488  end do
489  end function add_arr_mult_3_1
491  integer function add_arr_mult_1_1(arr_1,arr_2,arr_3,deriv) result(ierr)
492  character(*), parameter :: rout_name = 'add_arr_mult_1_1'
493 
494  ! input / output
495  real(dp), intent(in) :: arr_1(1:,0:)
496  real(dp), intent(in) :: arr_2(1:,0:)
497  real(dp), intent(out) :: arr_3(1:)
498  integer, intent(in) :: deriv(3)
499 
500  ! local variables
501  real(dp) :: bin_fac ! binomial factor for norm. sum
502  integer :: r ! current degree of norm. derivative
503  character(len=max_str_ln) :: err_msg ! error message
504 
505  ! initialize ierr
506  ierr = 0
507 
508  ! tests
509  if (size(arr_1,1).ne.size(arr_2,1)) then
510  err_msg = 'Arrays 1 and 2 need to have the same size in the &
511  &flux coord.'
512  ierr = 1
513  chckerr(err_msg)
514  end if
515  if (size(arr_1,1).ne.size(arr_3,1)) then
516  err_msg = 'Array 1 needs to have the same size as the resulting &
517  &array 3'
518  ierr = 1
519  chckerr(err_msg)
520  end if
521 
522  if (deriv(2).gt.0 .or. deriv(3).gt.0) then ! no addition to arr_3
523  return
524  else
525  ! distribute normal derivatives using binomial theorem, angular
526  ! derivatives only operate on the first array
527  ! normal derivatives
528  do r = 0,deriv(1)
529  if (r.eq.0) then ! first term in sum
530  bin_fac = 1.0_dp
531  else
532  bin_fac = bin_fac*(deriv(1)-(r-1))/r
533  end if
534 
535  ! current term in the tripple summation
536  arr_3 = arr_3 + bin_fac * arr_1(:,r)* arr_2(:,deriv(1)-r)
537  end do
538  end if
539  end function add_arr_mult_1_1
540 
542  integer recursive function calc_det_3d(deta,a,n) result (ierr)
543  character(*), parameter :: rout_name = 'calc_det_3D'
544 
545  ! input / output
546  real(dp), intent(inout) :: deta(:,:,:)
547  real(dp), intent(in) :: a(:,:,:,:)
548  integer, intent(in) :: n
549 
550  ! local variables
551  logical :: sym ! .true. if matrix symmetric
552  integer :: id, jd, kd ! counter
553  integer :: nn ! nr. of elements in matrix
554  real(dp) :: sgn ! holds either plus or minus one
555  integer, allocatable :: slct(:) ! 0 in between 1's selects which column to delete
556  integer, allocatable :: idx(:) ! counts from 1 to size(A)
557  character(len=max_str_ln) :: err_msg ! error message
558  real(dp), allocatable :: work(:,:,:) ! work array
559  integer, allocatable :: c_sub(:) ! indices of submatrix in storage convention of eq_vars.eq_2_type
560  integer, allocatable :: k_sub(:) ! first indices of submatrix in 2D
561 
562  ! initialize ierr
563  ierr = 0
564 
565  ! tests
566  if (size(deta,1).ne.size(a,1) .or. size(deta,2).ne.size(a,2) .or. &
567  &size(deta,3).ne.size(a,3)) then
568  err_msg = 'The output and input matrix have to have the same sizes'
569  ierr = 1
570  chckerr(err_msg)
571  end if
572 
573  ! set total number of elements in matrix
574  nn = size(a,4)
575 
576  ! set sym
577  ierr = is_sym(n,nn,sym)
578  chckerr('')
579 
580  ! intialize other quantities
581  sgn = 1.0_dp
582  allocate (idx(n))
583  idx = [(id,id=1,n)]
584  allocate (slct(n))
585  slct = 1
586  allocate (work(size(a,1),size(a,2),size(a,3)))
587 
588  deta = 0.0_dp
589 
590  if (n.eq.1) then ! shouldn't be used
591  deta = a(:,:,:,c([1,1],sym,n))
592  else if (n.eq.2) then
593  deta = a(:,:,:,c([1,1],sym,n))*a(:,:,:,c([2,2],sym,n))-&
594  &a(:,:,:,c([1,2],sym,n))*a(:,:,:,c([2,1],sym,n))
595  else
596  ! allocate coordinates of (n-1)x(n-1) submatrix
597  allocate(c_sub((n-1)**2))
598  allocate(k_sub(n**2)); k_sub = 0
599  ! iterate over indices in second dimension
600  do id = 1, n
601  slct(id) = 0
602  ! set up coordinates of submatrix (2:n,pack(idx,slct.gt.0))
603  ! (in general not symmetric, even if parent matrix is)
604  k_sub = pack(idx,slct.gt.0)
605  do kd = 1,n-1
606  do jd = 1,n-1
607  c_sub((kd-1)*(n-1)+jd) = c([jd+1,k_sub(kd)],sym,n)
608  end do
609  end do
610  ierr = calc_det_3d(work,a(:,:,:,c_sub),n-1)
611  chckerr('')
612  deta = deta + a(:,:,:,c([1,id],sym,n))*sgn*work
613  sgn = -sgn
614  slct(id) = 1
615  end do
616  ! deallocate c_sub and k_sub
617  deallocate(c_sub,k_sub)
618  end if
619  end function calc_det_3d
621  integer function calc_det_0d(det_0D,A) result(ierr)
622  character(*), parameter :: rout_name = 'calc_det_0D'
623 
624  ! input / output
625  real(dp), intent(inout) :: det_0d
626  real(dp), intent(in) :: a(:,:)
627 
628  ! local variables
629  integer :: n
630  real(dp) :: sgn
631  integer :: id, info
632  integer, allocatable :: ipiv(:)
633  character(len=max_str_ln) :: err_msg ! error message
634  real(dp), allocatable :: a_loc(:,:) ! copy of A
635 
636  ! initialize ierr
637  ierr = 0
638 
639  ! tests
640  if (size(a,1).ne.size(a,2)) then
641  err_msg = 'The matrix A has to be square'
642  ierr = 1
643  chckerr(err_msg)
644  end if
645 
646  ! initialize
647  n = size(a,1)
648  allocate(ipiv(n))
649  ipiv = 0
650 
651  ! set up local A
652  allocate(a_loc(n,n))
653  a_loc = a
654 
655  call dgetrf(n,n,a_loc,n,ipiv,info)
656 
657  det_0d = 1.0_dp
658  do id = 1, n
659  det_0d = det_0d*a_loc(id,id)
660  end do
661 
662  sgn = 1.0_dp
663  do id = 1, n
664  if(ipiv(id).ne.id) then
665  sgn = -sgn
666  end if
667  end do
668  det_0d = sgn*det_0d
669 
670  ! deallocate local variables
671  deallocate(a_loc)
672  end function calc_det_0d
673 
675  integer function calc_inv_3d(inv_3D,A,n) result(ierr)
676  use num_vars, only: tol_zero
677 
678  character(*), parameter :: rout_name = 'calc_inv_3D'
679 
680  ! input / output
681  real(dp), intent(inout) :: inv_3d(:,:,:,:)
682  real(dp), intent(in) :: a(:,:,:,:)
683  integer, intent(in) :: n
684 
685  ! local variables
686  logical :: sym ! .true. if matrix symmetric
687  real(dp), allocatable :: deta(:,:,:) ! determinant of a submatrix of A
688  integer :: id, jd, kd, ld ! counters
689  integer :: nn ! nr. of elements in matrix
690  integer, allocatable :: slct(:,:) ! 0 in between 1's selects which column to delete
691  integer, allocatable :: idx(:) ! counts from 1 to size(A)
692  character(len=max_str_ln) :: err_msg ! error message
693  integer, allocatable :: c_sub(:) ! indices of submatrix in storage convention of eq_vars.eq_2_type
694  integer, allocatable :: l_sub(:), k_sub(:) ! first and second indices of submatrix in 2D
695  integer :: kd_min ! min. of kd
696 
697  ! initialize ierr
698  ierr = 0
699 
700  ! tests
701  if (size(inv_3d,1).ne.size(a,1) .or. size(inv_3d,2).ne.size(a,2) .or. &
702  &size(inv_3d,3).ne.size(a,3) .or. size(inv_3d,4).ne.size(a,4)) then
703  err_msg = 'The output and input matrix have to have the same sizes'
704  ierr = 1
705  chckerr(err_msg)
706  end if
707 
708  ! set total number of elements in matrix
709  nn = size(a,4)
710 
711  ! set sym
712  ierr = is_sym(n,nn,sym)
713  chckerr('')
714 
715  ! initialize other quantitites
716  allocate (idx(n))
717  idx = [(id,id=1,n)]
718  allocate (slct(n,2))
719  slct = 1
720  allocate(deta(size(a,1),size(a,2),size(a,3)))
721  deta = 0.0_dp
722 
723  ! allocate coordinates of (n-1)x(n-1) submatrix
724  allocate(c_sub((n-1)**2))
725  allocate(l_sub(n**2)); l_sub = 0
726  allocate(k_sub(n**2)); k_sub = 0
727 
728  ! set up kd_min
729  kd_min = 1
730 
731  ! calculate cofactor matrix, replacing original elements
732  ! use is made of detA
733  do ld = 1,n
734  ! adjust kd_min if symmetric
735  if (sym) kd_min = ld
736  do kd = kd_min,n
737  ! set up coordinates of submatrix (strike out row i, column j)
738  slct(ld,1) = 0 ! take out column ld
739  slct(kd,2) = 0 ! take out row kd
740  l_sub = pack(idx,slct(:,1).gt.0)
741  k_sub = pack(idx,slct(:,2).gt.0)
742  do jd = 1,n-1
743  do id = 1,n-1
744  c_sub((jd-1)*(n-1)+id) = c([l_sub(id),k_sub(jd)],sym,n) ! taking the transpose!
745  end do
746  end do
747  ierr = calc_det_3d(deta,a(:,:,:,c_sub),n-1)
748  chckerr('')
749  inv_3d(:,:,:,c([kd,ld],sym,n)) = (-1.0_dp)**(kd+ld)*deta
750  slct(ld,1) = 1
751  slct(kd,2) = 1
752  end do
753  end do
754 
755  ! deallocate c_sub l_sub and k_sub
756  deallocate(c_sub,l_sub,k_sub)
757 
758  ! calculate determinant in detA
759  ierr = calc_det(deta,a,n)
760  chckerr('')
761 
762  ! limit determinant from being zero
763  deta = sign(max(abs(deta),tol_zero),deta)
764 
765  ! divide by determinant
766  do kd = 1,nn
767  inv_3d(:,:,:,kd) = inv_3d(:,:,:,kd) / deta
768  end do
769  end function calc_inv_3d
771  integer function calc_inv_0d(inv_0D,A) result(ierr)
772  character(*), parameter :: rout_name = 'calc_inv_0D'
773 
774  ! input / output
775  real(dp), intent(inout) :: inv_0d(:,:)
776  real(dp), intent(in) :: a(:,:)
777 
778  ! local variables
779  integer :: n
780  integer, allocatable :: ipiv(:) ! pivot variable, used by Lapack
781  real(dp), allocatable :: w(:) ! work variable
782  character(len=max_str_ln) :: err_msg ! error message
783 
784  ! initialize ierr
785  ierr = 0
786 
787  ! tests
788  if (size(a,1).ne.size(a,1)) then
789  err_msg = 'The input matrix A has to be square'
790  ierr = 1
791  chckerr(err_msg)
792  end if
793  if (size(inv_0d,1).ne.size(a,1) .or. size(inv_0d,2).ne.size(a,2)) then
794  err_msg = 'The output and input matrix have to have the same sizes'
795  ierr = 1
796  chckerr(err_msg)
797  end if
798 
799  ! initialize
800  n = size(a,1)
801  allocate(ipiv(n))
802  ipiv = 0
803  allocate(w(n))
804 
805  inv_0d = a
806 
807  call dgetrf(n,n,inv_0d,n,ipiv,ierr) ! LU factorization
808  err_msg = 'Lapack couldn''t find the LU factorization'
809  chckerr(err_msg)
810 
811  call dgetri(n,inv_0d,n,ipiv,w,n,ierr) ! inverse of LU
812  chckerr('Lapack couldn''t compute the inverse')
813  end function calc_inv_0d
814 
816  integer function calc_mult_3d_real(A,B,AB,n,transp) result(ierr)
817  character(*), parameter :: rout_name = 'calc_mult_3D_real'
818 
819  ! input / output
820  real(dp), intent(in) :: a(:,:,:,:)
821  real(dp), intent(in) :: b(:,:,:,:)
822  real(dp), intent(inout) :: ab(:,:,:,:)
823  integer, intent(in) :: n
824  logical, intent(in), optional :: transp(2)
825 
826  ! local variables
827  logical :: sym(3) ! .true. if matrices symmetric
828  integer :: nn(3) ! nr. of elements in matrices
829  character(len=max_str_ln) :: err_msg ! error message
830  integer :: id, jd, kd ! counters
831  integer :: id_min ! min. of id
832  integer :: c1, c2, c3 ! coordinates
833  integer :: ind(2,3) ! indices
834  logical :: transp_loc(2) ! local copy of transp
835 
836  ! initialize ierr
837  ierr = 0
838 
839  ! tests
840  if (size(a,1).ne.size(b,1) .or. size(a,2).ne.size(b,2) .or. &
841  &size(a,3).ne.size(b,3) .or. size(a,1).ne.size(ab,1) .or. &
842  &size(a,2).ne.size(ab,2) .or. size(a,3).ne.size(ab,3)) then
843  err_msg = 'The output and input matrices have to defined on the &
844  &same grid'
845  ierr = 1
846  chckerr(err_msg)
847  end if
848 
849  ! set local transp
850  transp_loc = .false.
851  if (present(transp)) transp_loc = transp
852 
853  ! set total number of elements in matrices
854  nn(1) = size(a,4)
855  nn(2) = size(b,4)
856  nn(3) = size(ab,4)
857 
858  ! set sym
859  do id = 1,3
860  ierr = is_sym(n,nn(id),sym(id))
861  chckerr('')
862  end do
863 
864  ! initialize AB
865  ab = 0._dp
866 
867  ! set up id_min
868  id_min = 1
869 
870  ! loop over all 2D rows and columns of AB
871  do jd = 1,n
872  ! adjust id_min if AB symmetric
873  if (sym(3)) id_min = jd
874  do id = id_min,n
875  do kd = 1,n
876  ind(:,3) = [id,jd]
877  c3 = c(ind(:,3),sym(3),n)
878  if (transp_loc(1)) then
879  ind(:,1) = [kd,id]
880  else
881  ind(:,1) = [id,kd]
882  end if
883  c1 = c(ind(:,1),sym(1),n)
884  if (transp_loc(2)) then
885  ind(:,2) = [jd,kd]
886  else
887  ind(:,2) = [kd,jd]
888  end if
889  c2 = c(ind(:,2),sym(2),n)
890  ab(:,:,:,c3) = ab(:,:,:,c3) + &
891  &a(:,:,:,c1)*b(:,:,:,c2)
892  end do
893  end do
894  end do
895  end function calc_mult_3d_real
897  integer function calc_mult_3d_complex(A,B,AB,n,transp) result(ierr)
898  character(*), parameter :: rout_name = 'calc_mult_3D_complex'
899 
900  ! input / output
901  complex(dp), intent(in) :: a(:,:,:,:)
902  complex(dp), intent(in) :: b(:,:,:,:)
903  complex(dp), intent(inout) :: ab(:,:,:,:)
904  integer, intent(in) :: n
905  logical, intent(in), optional :: transp(2)
906 
907  ! local variables
908  logical :: sym(3) ! .true. if matrices symmetric
909  integer :: nn(3) ! nr. of elements in matrices
910  character(len=max_str_ln) :: err_msg ! error message
911  integer :: id, jd, kd ! counters
912  integer :: id_min ! min. of id
913  integer :: c1, c2, c3 ! coordinates
914  logical :: transp_loc(2) ! local copy of transp
915  integer :: ind(2,3) ! indices
916  integer :: d(3) ! dimension of matrices A and B
917 
918  ! initialize ierr
919  ierr = 0
920 
921  ! tests
922  if (size(a,1).ne.size(b,1) .or. size(a,2).ne.size(b,2) .or. &
923  &size(a,3).ne.size(b,3) .or. size(a,1).ne.size(ab,1) .or. &
924  &size(a,2).ne.size(ab,2) .or. size(a,3).ne.size(ab,3)) then
925  err_msg = 'The output and input matrices have to defined on the &
926  &same grid'
927  ierr = 1
928  chckerr(err_msg)
929  end if
930 
931  ! set local transp
932  transp_loc = .false.
933  if (present(transp)) transp_loc = transp
934 
935  ! set d
936  d = [size(a,1),size(a,2),size(a,3)]
937 
938  ! set total number of elements in matrices
939  nn(1) = size(a,4)
940  nn(2) = size(b,4)
941  nn(3) = size(ab,4)
942 
943  ! set sym
944  do id = 1,3
945  ierr = is_sym(n,nn(id),sym(id))
946  chckerr('')
947  end do
948 
949  ! initialize AB
950  ab = 0._dp
951 
952  ! set up id_min
953  id_min = 1
954 
955  ! loop over all 2D rows and columns of AB
956  do jd = 1,n
957  ! adjust id_min if AB symmetric
958  if (sym(3)) id_min = jd
959  do id = id_min,n
960  do kd = 1,n
961  ind(:,3) = [id,jd]
962  c3 = c(ind(:,3),sym(3),n)
963  if (transp_loc(1)) then
964  ind(:,1) = [kd,id]
965  else
966  ind(:,1) = [id,kd]
967  end if
968  c1 = c(ind(:,1),sym(1),n)
969  if (transp_loc(2)) then
970  ind(:,2) = [jd,kd]
971  else
972  ind(:,2) = [kd,jd]
973  end if
974  c2 = c(ind(:,2),sym(2),n)
975  ab(:,:,:,c3) = ab(:,:,:,c3) + &
976  &con(a(:,:,:,c1),ind(:,1),sym(1),d)*&
977  &con(b(:,:,:,c2),ind(:,2),sym(2),d)
978  end do
979  end do
980  end do
981  end function calc_mult_3d_complex
983  integer function calc_mult_0d_real(A,B,AB,n,transp) result(ierr)
984  character(*), parameter :: rout_name = 'calc_mult_0D_real'
985 
986  ! input / output
987  real(dp), intent(in) :: a(:)
988  real(dp), intent(in) :: b(:)
989  real(dp), intent(inout) :: ab(:)
990  integer, intent(in) :: n
991  logical, intent(in), optional :: transp(2)
992 
993  ! local variables
994  real(dp), allocatable :: loc_ab(:,:,:,:) ! local copy of C
995 
996  ! initialize ierr
997  ierr = 0
998 
999  ! allocate local AB
1000  allocate(loc_ab(1,1,1,size(ab)))
1001 
1002  ! calculate using 2D version
1003  ierr = calc_mult_3d_real(reshape(a,[1,1,1,size(a)]),&
1004  &reshape(b,[1,1,1,size(b)]),loc_ab,n,transp)
1005 
1006  ! update AB with local value
1007  ab = loc_ab(1,1,1,:)
1008 
1009  ! deallocate local variables
1010  deallocate(loc_ab)
1011  end function calc_mult_0d_real
1012 
1014  integer function conv_mat_3d(A,B,n,transp) result(ierr)
1015  character(*), parameter :: rout_name = 'conv_mat_3D'
1016 
1017  ! input / output
1018  real(dp), intent(in) :: a(:,:,:,:)
1019  real(dp), intent(inout) :: b(:,:,:,:)
1020  integer, intent(in) :: n
1021  logical, intent(in), optional :: transp
1022 
1023  ! local variables
1024  logical :: sym(2) ! .true. if matrices symmetric
1025  logical :: transp_loc ! local transp
1026  integer :: nn(2) ! nr. of elements in matrices
1027  integer :: id, jd ! counters
1028  integer :: id_min ! minimum value of id
1029  integer :: ind(2,2) ! indices
1030  real(dp), allocatable :: a_loc(:,:,:,:) ! local copy of A
1031 
1032  ! initialize ierr
1033  ierr = 0
1034 
1035  ! set total number of elements in matrix
1036  nn(1) = size(a,4)
1037  nn(2) = size(b,4)
1038 
1039  ! set local transp
1040  transp_loc = .false.
1041  if (present(transp)) transp_loc = transp
1042 
1043  ! set local A
1044  allocate(a_loc(size(a,1),size(a,2),size(a,3),size(a,4)))
1045  a_loc = a
1046 
1047  ! set sym
1048  do id = 1,2
1049  ierr = is_sym(n,nn(id),sym(id))
1050  chckerr('')
1051  end do
1052 
1053  ! copy A into B
1054  if (sym(1).and.sym(2) .or. &
1055  &.not.sym(1).and..not.sym(2).and..not.transp_loc) then ! both symmetric or both unsymmetric and no transpose needed
1056  b = a ! symmetric matrix equal to transpose or no tranpose needed
1057  else ! at least one is not symmetric
1058  ! loop over row
1059  do jd = 1,n
1060  ! set id_min
1061  if (sym(2)) then
1062  id_min = jd
1063  else
1064  id_min = 1
1065  end if
1066  ! loop over column
1067  do id = id_min,n
1068  ! set ind
1069  if (transp_loc) then
1070  ind(:,1) = [jd,id]
1071  else
1072  ind(:,1) = [id,jd]
1073  end if
1074  ind(:,2) = [id,jd]
1075  ! copy elements of A into B
1076  b(:,:,:,c(ind(:,2),sym(2),n)) = &
1077  &a_loc(:,:,:,c(ind(:,1),sym(1),n))
1078  end do
1079  end do
1080  end if
1081  end function conv_mat_3d
1083  integer function conv_mat_0d(A,B,n,transp) result(ierr)
1084  character(*), parameter :: rout_name = 'conv_mat_0D'
1085 
1086  ! input / output
1087  real(dp), intent(in) :: a(:)
1088  real(dp), intent(inout) :: b(:)
1089  integer, intent(in) :: n
1090  logical, intent(in), optional :: transp
1091 
1092  ! local variables
1093  real(dp), allocatable :: b_loc(:,:,:,:) ! local B
1094 
1095  ! initialize ierr
1096  ierr = 0
1097 
1098  ! set up local B
1099  allocate(b_loc(1,1,1,size(b)))
1100 
1101  ! call 3D version
1102  ierr = conv_mat_3d(reshape(a,[1,1,1,size(a)]),b_loc,n,transp)
1103  chckerr('')
1104 
1105  ! copy to B
1106  b = b_loc(1,1,1,:)
1107  end function conv_mat_0d
1109  integer function conv_mat_3d_complex(A,B,n,transp) result(ierr)
1110  character(*), parameter :: rout_name = 'conv_mat_3D_complex'
1111 
1112  ! input / output
1113  complex(dp), intent(in) :: a(:,:,:,:)
1114  complex(dp), intent(inout) :: b(:,:,:,:)
1115  integer, intent(in) :: n
1116  logical, intent(in), optional :: transp
1117 
1118  ! local variables
1119  logical :: sym(2) ! .true. if matrices symmetric
1120  logical :: transp_loc ! local transp
1121  integer :: nn(2) ! nr. of elements in matrices
1122  integer :: id, jd ! counters
1123  integer :: id_min ! minimum value of id
1124  integer :: ind(2,2) ! indices
1125  complex(dp), allocatable :: a_loc(:,:,:,:) ! local copy of A
1126 
1127  ! initialize ierr
1128  ierr = 0
1129 
1130  ! set total number of elements in matrix
1131  nn(1) = size(a,4)
1132  nn(2) = size(b,4)
1133 
1134  ! set local transp
1135  transp_loc = .false.
1136  if (present(transp)) transp_loc = transp
1137 
1138  ! set local A
1139  allocate(a_loc(size(a,1),size(a,2),size(a,3),size(a,4)))
1140  a_loc = a
1141 
1142  ! set sym
1143  do id = 1,2
1144  ierr = is_sym(n,nn(id),sym(id))
1145  chckerr('')
1146  end do
1147 
1148  ! copy A into B
1149  if (sym(1).and.sym(2) .or. &
1150  &.not.sym(1).and..not.sym(2).and..not.transp_loc) then ! both symmetric or both unsymmetric and no transpose needed
1151  b = a ! symmetric matrix equal to transpose or no tranpose needed
1152  else ! at least one is not symmetric
1153  ! loop over row
1154  do jd = 1,n
1155  ! set id_min
1156  if (sym(2)) then
1157  id_min = jd
1158  else
1159  id_min = 1
1160  end if
1161  ! loop over column
1162  do id = id_min,n
1163  ! set ind
1164  if (transp_loc) then
1165  ind(:,1) = [jd,id]
1166  else
1167  ind(:,1) = [id,jd]
1168  end if
1169  ind(:,2) = [id,jd]
1170  ! copy elements of A into B
1171  b(:,:,:,c(ind(:,2),sym(2),n)) = &
1172  &a_loc(:,:,:,c(ind(:,1),sym(1),n))
1173  end do
1174  end do
1175  end if
1176  end function conv_mat_3d_complex
1178  integer function conv_mat_0d_complex(A,B,n,transp) result(ierr)
1179  character(*), parameter :: rout_name = 'conv_mat_0D_complex'
1180 
1181  ! input / output
1182  complex(dp), intent(in) :: a(:)
1183  complex(dp), intent(inout) :: b(:)
1184  integer, intent(in) :: n
1185  logical, intent(in), optional :: transp
1186 
1187  ! local variables
1188  complex(dp), allocatable :: b_loc(:,:,:,:) ! local B
1189 
1190  ! initialize ierr
1191  ierr = 0
1192 
1193  ! set up local B
1194  allocate(b_loc(1,1,1,size(b)))
1195 
1196  ! call 3D version
1197  ierr = conv_mat_3d_complex(reshape(a,[1,1,1,size(a)]),b_loc,n,transp)
1198  chckerr('')
1199 
1200  ! copy to B
1201  b = b_loc(1,1,1,:)
1202  end function conv_mat_0d_complex
1203 
1205  integer function con2dis_eqd(pt_c,pt_d,lim_c,lim_d) result(ierr)
1206  !character(*), parameter :: rout_name = 'con2dis_eqd'
1207 
1208  ! input / output
1209  real(dp), intent(in) :: pt_c
1210  real(dp), intent(inout) :: pt_d
1211  real(dp), intent(in) :: lim_c(2)
1212  integer, intent(in) :: lim_d(2)
1213 
1214  ! initialize ierr
1215  ierr = 0
1216 
1217  ! Calculate, using the formula:
1218  ! pt_c - min_c pt_d - min_d
1219  ! ------------- = ------------- ,
1220  ! max_c - min_c max_d - min_d
1221  pt_d = lim_d(1) + (lim_d(2)-lim_d(1)) * (pt_c-lim_c(1)) / &
1222  &(lim_c(2)-lim_c(1))
1223  end function con2dis_eqd
1225  integer function con2dis_reg(pt_c,pt_d,var_c) result(ierr)
1226  character(*), parameter :: rout_name = 'con2dis_reg'
1227 
1228  ! input / output
1229  real(dp), intent(in) :: pt_c
1230  real(dp), intent(inout) :: pt_d
1231  real(dp), intent(in) :: var_c(:)
1232 
1233  ! local variables
1234  real(dp), allocatable :: var_c_loc(:) ! local copy of var
1235  real(dp), allocatable :: var_c_inv(:) ! inverted var_c
1236  integer :: size_c ! size of var_c
1237  integer :: ind_lo, ind_hi ! lower and upper index comprising pt_c
1238  integer :: id ! counter
1239  real(dp) :: pt_c_loc ! local pt_c
1240  character(len=max_str_ln) :: err_msg ! error message
1241  real(dp), parameter :: tol = 1.e-8_dp ! tolerance for comparisons
1242 
1243  ! initialize ierr
1244  ierr = 0
1245 
1246  ! set up size_c
1247  size_c = size(var_c)
1248 
1249 #if ldebug
1250  if (debug_con2dis_reg) &
1251  &call writo('finding the discrete index of '//trim(r2str(pt_c)))
1252 #endif
1253 
1254  ! set up local var_c and pt_c_loc
1255  allocate(var_c_loc(size_c))
1256  if (var_c(1).lt.var_c(size_c)) then
1257  var_c_loc = var_c
1258  pt_c_loc = pt_c
1259  else
1260 #if ldebug
1261  if (debug_con2dis_reg) &
1262  &call writo('setting var_c_loc = - var_c and pt_c_loc = -pt_c')
1263 #endif
1264  var_c_loc = - var_c ! local var should be rising
1265  pt_c_loc = - pt_c ! invert pt_c as well
1266  end if
1267 
1268  ! set up inverse var_c
1269  allocate(var_c_inv(size_c))
1270  var_c_inv = var_c_loc(size_c:1:-1)
1271 
1272  ! find the lower and upper index comprising local pt_c
1273  ind_lo = 0
1274  ind_hi = size_c+1
1275  do id = 1,size_c
1276  if (pt_c_loc+tol*abs(pt_c_loc) .ge. var_c_loc(id)) then
1277  ind_lo = id
1278 #if ldebug
1279  if (debug_con2dis_reg) call writo('for iteration '//&
1280  &trim(i2str(id))//'/'//trim(i2str(size_c))//', '//&
1281  &trim(r2strt(pt_c_loc))//' + '//trim(r2strt(tol))//&
1282  &'*abs('//trim(r2strt(pt_c_loc))//') >= '//&
1283  &trim(r2strt(var_c_loc(id)))//' => ['//&
1284  &trim(i2str(ind_lo))//','//trim(i2str(ind_hi))//']')
1285 #endif
1286  end if
1287  if (pt_c_loc-tol*abs(pt_c_loc) .le. var_c_inv(id)) then
1288  ind_hi = size_c+1-id
1289 #if ldebug
1290  if (debug_con2dis_reg) call writo('for iteration '//&
1291  &trim(i2str(id))//'/'//trim(i2str(size_c))//', '//&
1292  &trim(r2strt(pt_c_loc))//' - '//trim(r2strt(tol))//&
1293  &'*abs('//trim(r2strt(pt_c_loc))//') < '//&
1294  &trim(r2strt(var_c_inv(id)))//' => ['//&
1295  &trim(i2str(ind_lo))//','//trim(i2str(ind_hi))//']')
1296 #endif
1297  end if
1298  end do
1299 #if ldebug
1300  if (debug_con2dis_reg) then
1301  call writo('final ind_lo = '//trim(i2str(ind_lo))//', ind_hi = '//&
1302  &trim(i2str(ind_hi)))
1303  !!call print_ex_2D('var_c_loc, var_c_inv','',&
1304  !!&reshape([var_c_loc,var_c_inv],[size_c,2]))
1305  end if
1306 #endif
1307 
1308  ! tests
1309  if (ind_lo.eq.0 .or. ind_hi.eq.size_c+1) then ! not within range
1310  pt_d = -1._dp
1311  ierr = 1
1312  err_msg = 'pt_c '//trim(r2str(pt_c))//' not within range '//&
1313  &trim(r2str(minval(var_c)))//'..'//trim(r2str(maxval(var_c)))
1314  chckerr(err_msg)
1315  end if
1316 
1317  ! set output
1318  if (ind_lo.lt.ind_hi) then ! valid output that does not correspond to a point on grid
1319  pt_d = ind_lo + (pt_c_loc-var_c_loc(ind_lo))/&
1320  &(var_c_loc(ind_hi)-var_c_loc(ind_lo))
1321  else if (ind_lo.eq.ind_hi) then ! valid output that does correspond to a point on grid
1322  pt_d = ind_lo
1323  else ! invalid output
1324  pt_d = -1._dp
1325  ierr = 1
1326  err_msg = 'ind_lo cannot be higher than ind_hi'
1327  chckerr(err_msg)
1328  end if
1329 
1330  ! deallocate
1331  deallocate(var_c_loc,var_c_inv)
1332  end function con2dis_reg
1333 
1335  integer function dis2con_eqd(pt_d,pt_c,lim_d,lim_c) result(ierr)
1336  !character(*), parameter :: rout_name = 'dis2con_eqd'
1337 
1338  ! input / output
1339  integer, intent(in) :: pt_d
1340  real(dp), intent(inout) :: pt_c
1341  integer, intent(in) :: lim_d(2)
1342  real(dp), intent(in) :: lim_c(2)
1343 
1344  ! initialize ierr
1345  ierr = 0
1346 
1347  ! Calculate, using the formula:
1348  ! pt_c - min_c pt_d - min_d
1349  ! ------------- = ------------- ,
1350  ! max_c - min_c max_d - min_d
1351  pt_c = lim_c(1) + (lim_c(2)-lim_c(1)) * (pt_d-lim_d(1)) / &
1352  &(lim_d(2)-lim_d(1))
1353  end function dis2con_eqd
1355  integer function dis2con_reg(pt_d,pt_c,var_c) result(ierr)
1356  character(*), parameter :: rout_name = 'dis2con_reg'
1357 
1358  ! input / output
1359  integer, intent(in) :: pt_d
1360  real(dp), intent(inout) :: pt_c
1361  real(dp), intent(in) :: var_c(:)
1362 
1363  ! local variables
1364  character(len=max_str_ln) :: err_msg ! error message
1365 
1366  ! initialize ierr
1367  ierr = 0
1368 
1369  ! Check whether the discrete value lies inside the range
1370  if (pt_d.lt.1 .or. pt_d.gt.size(var_c)) then
1371  pt_c = -1._dp
1372  ierr = 1
1373  err_msg = 'pt_c not within range'
1374  chckerr(err_msg)
1375  end if
1376 
1377  ! Return the continuous variable
1378  pt_c = var_c(pt_d)
1379  end function dis2con_reg
1380 
1382  integer function round_with_tol_arr(vals,lim_lo,lim_hi,tol) result(ierr)
1383  character(*), parameter :: rout_name = 'round_with_tol_arr'
1384 
1385  ! input / output
1386  real(dp), intent(inout) :: vals(:)
1387  real(dp), intent(in) :: lim_lo
1388  real(dp), intent(in) :: lim_hi
1389  real(dp), intent(in), optional :: tol
1390 
1391  ! local variables
1392  real(dp) :: loc_tol ! local value of tol
1393  character(len=max_str_ln) :: err_msg ! error message
1394 
1395  ! initialize ierr
1396  ierr = 0
1397 
1398  ! set loc_tol
1399  if (present(tol)) then
1400  loc_tol = tol
1401  else
1402  loc_tol = 1e-5_dp
1403  end if
1404 
1405  ! set possible error message
1406  err_msg = 'value has to lie within 0..1'
1407 
1408  if (minval(vals-lim_lo+loc_tol).lt.0) then
1409  ierr = 1
1410  chckerr(err_msg)
1411  else if (maxval(vals-lim_hi-loc_tol).gt.0) then
1412  ierr = 1
1413  chckerr(err_msg)
1414  else
1415  where (vals.lt.lim_lo) vals = lim_lo
1416  where (vals.gt.lim_hi) vals = lim_hi
1417  end if
1418  end function round_with_tol_arr
1420  integer function round_with_tol_ind(val,lim_lo,lim_hi,tol) result(ierr)
1421  character(*), parameter :: rout_name = 'round_with_tol_ind'
1422 
1423  ! input / output
1424  real(dp), intent(inout) :: val
1425  real(dp), intent(in) :: lim_lo
1426  real(dp), intent(in) :: lim_hi
1427  real(dp), intent(in), optional :: tol
1428 
1429  ! local variables
1430  real(dp) :: vals(1) ! array copy of val
1431 
1432  ! initialize ierr
1433  ierr = 0
1434 
1435  vals(1) = val
1436 
1437  ierr = round_with_tol_arr(vals,lim_lo,lim_hi,tol)
1438  chckerr('')
1439 
1440  val = vals(1)
1441  end function round_with_tol_ind
1442 
1444  function con_3d(A,c,sym,d) result(B)
1445  complex(dp), intent(in) :: a(:,:,:)
1446  integer, intent(in) :: c(2)
1447  logical, intent(in) :: sym
1448  integer, intent(in) :: d(3)
1449  complex(dp) :: b(d(1),d(2),d(3))
1450 
1451  ! determine whether to take complex conjugate or not
1452  if (c(2).gt.c(1) .and. sym) then ! upper triangular matrix
1453  b = conjg(a)
1454  else
1455  b = a
1456  end if
1457  end function con_3d
1459  function con_2d(A,c,sym,d) result(B)
1460  complex(dp), intent(in) :: a(:,:)
1461  integer, intent(in) :: c(2)
1462  logical, intent(in) :: sym
1463  integer, intent(in) :: d(2)
1464  complex(dp) :: b(d(1),d(2))
1465 
1466  ! determine whether to take complex conjugate or not
1467  if (c(2).gt.c(1) .and. sym) then ! upper triangular matrix
1468  b = conjg(a)
1469  else
1470  b = a
1471  end if
1472  end function con_2d
1474  function con_1d(A,c,sym,d) result(B)
1475  complex(dp), intent(in) :: a(:)
1476  integer, intent(in) :: c(2)
1477  logical, intent(in) :: sym
1478  integer, intent(in) :: d(1)
1479  complex(dp) :: b(d(1))
1480 
1481  ! determine whether to take complex conjugate or not
1482  if (c(2).gt.c(1) .and. sym) then ! upper triangular matrix
1483  b = conjg(a)
1484  else
1485  b = a
1486  end if
1487  end function con_1d
1489  function con_0d(A,c,sym) result(B)
1490  complex(dp), intent(in) :: a
1491  integer, intent(in) :: c(2)
1492  logical, intent(in) :: sym
1493  complex(dp) :: b
1494 
1495  ! determine whether to take complex conjugate or not
1496  if (c(2).gt.c(1) .and. sym) then ! upper triangular matrix
1497  b = conjg(a)
1498  else
1499  b = a
1500  end if
1501  end function con_0d
1502 
1504  subroutine bubble_sort_real(a,piv)
1505  ! input / output
1506  real(dp), intent(inout) :: a(:)
1507  integer, intent(inout), optional :: piv(:)
1508 
1509  ! local variables
1510  real(dp) :: temp ! temporary value
1511  integer :: id, jd ! counter
1512  logical :: swapped ! value swapped with next
1513  logical :: save_piv ! save pivots
1514 
1515  save_piv = .false.
1516  if (present(piv)) then
1517  if (size(piv).eq.size(a)) then
1518  save_piv = .true.
1519  piv = [(jd,jd=1,size(a))]
1520  end if
1521  end if
1522 
1523  do jd = size(a)-1,1,-1
1524  swapped = .false.
1525  do id = 1, jd
1526  if (a(id) .gt. a(id+1)) then
1527  temp = a(id)
1528  a(id) = a(id+1)
1529  a(id+1) = temp
1530  piv(id:id+1) = piv(id+1:id:-1)
1531  swapped = .true.
1532  end if
1533  end do
1534  if (.not. swapped) exit
1535  end do
1536  end subroutine bubble_sort_real
1538  subroutine bubble_sort_int(a,piv)
1539  ! input / output
1540  integer, intent(inout) :: a(:)
1541  integer, intent(inout), optional :: piv(:)
1542 
1543  ! local variables
1544  integer :: temp ! temporary value
1545  integer :: id, jd ! counter
1546  logical :: swapped ! value swapped with next
1547  logical :: save_piv ! save pivots
1548 
1549  save_piv = .false.
1550  if (present(piv)) then
1551  if (size(piv).eq.size(a)) then
1552  save_piv = .true.
1553  piv = [(jd,jd=1,size(a))]
1554  end if
1555  end if
1556 
1557  do jd = size(a)-1,1,-1
1558  swapped = .false.
1559  do id = 1, jd
1560  if (a(id) .gt. a(id+1)) then
1561  temp = a(id)
1562  a(id) = a(id+1)
1563  a(id+1) = temp
1564  if (present(piv)) piv(id:id+1) = piv(id+1:id:-1)
1565  swapped = .true.
1566  end if
1567  end do
1568  if (.not. swapped) exit
1569  end do
1570  end subroutine bubble_sort_int
1571 
1573  integer function order_per_fun_1(x,y,x_out,y_out,overlap,tol) result(ierr)
1574  use num_vars, only: tol_zero
1575 
1576  character(*), parameter :: rout_name = 'order_per_fun_1'
1577 
1578  ! input / output
1579  real(dp), intent(in) :: x(:)
1580  real(dp), intent(in) :: y(:)
1581  real(dp), intent(inout), allocatable :: x_out(:)
1582  real(dp), intent(inout), allocatable :: y_out(:)
1583  integer, intent(in) :: overlap
1584  real(dp), intent(in), optional :: tol
1585 
1586  ! local variables
1587  real(dp) :: tol_loc ! local tolerance
1588  real(dp), allocatable :: x_fund(:) ! x on fundamental interval 0..2pi
1589  integer :: ml_x ! location of maximum of x
1590  integer :: lim(2) ! limits in x
1591  integer :: lim_loc(2) ! limits in loc_x
1592  character(len=max_str_ln) :: err_msg ! error message
1593 
1594  ! initialize ierr
1595  ierr = 0
1596 
1597  ! set local tol
1598  tol_loc = tol_zero
1599  if (present(tol)) tol_loc = tol
1600 
1601  ! tests
1602  if (abs(cos(x(size(x)))-cos(x(1))).lt.tol_loc .and. &
1603  &abs(sin(x(size(x)))-sin(x(1))).lt.tol_loc) then
1604  ierr = 1
1605  err_msg = 'First and last point cannot be identical. Remove one.'
1606  chckerr(err_msg)
1607  end if
1608 
1609  ! initialize
1610  allocate(x_fund(size(x)))
1611  x_fund = mod(x,2*pi)
1612  where(x_fund.lt.0._dp) x_fund = x_fund+2*pi
1613  ml_x = maxloc(x_fund,1)
1614 
1615  allocate(x_out(size(x)+2*overlap))
1616  allocate(y_out(size(x)+2*overlap))
1617 
1618  ! left overlap
1619  lim = [max(ml_x-overlap+1,1),ml_x] ! all possible points to the left of ml_x, u
1620  lim_loc = [overlap-lim(2)+lim(1),overlap]
1621  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))-2*pi
1622  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1623  lim_loc = [1,lim_loc(1)-1]
1624  lim = [size(x)-(lim_loc(2)-lim_loc(1)),size(x)]
1625  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))-2*pi
1626  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1627 
1628  ! bulk
1629  lim = [ml_x+1,size(x)] ! points after max, until last
1630  lim_loc = [overlap+1,overlap+lim(2)-lim(1)+1]
1631  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))
1632  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1633  lim = [1,ml_x] ! points from start, until max
1634  lim_loc = lim_loc(2) + [1,lim(2)-lim(1)+1]
1635  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))
1636  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1637 
1638  ! right overlap
1639  lim = [ml_x+1,min(ml_x+overlap,size(x))] ! all possible poinst to the right of ml_X
1640  lim_loc = lim_loc(2) + [1,lim(2)-lim(1)+1]
1641  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))+2*pi
1642  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1643  lim_loc = [lim_loc(2)+1,size(x)+2*overlap] ! remaining point loop around
1644  lim = [1,lim_loc(2)-lim_loc(1)+1]
1645  x_out(lim_loc(1):lim_loc(2)) = x_fund(lim(1):lim(2))+2*pi
1646  y_out(lim_loc(1):lim_loc(2)) = y(lim(1):lim(2))
1647  end function order_per_fun_1
1649  integer function order_per_fun_2(xy,xy_out,overlap,tol) result(ierr)
1650  character(*), parameter :: rout_name = 'order_per_fun_2'
1651 
1652  ! input / output
1653  real(dp), intent(in) :: xy(:,:)
1654  real(dp), intent(inout), allocatable :: xy_out(:,:)
1655  integer, intent(in) :: overlap
1656  real(dp), intent(in), optional :: tol
1657 
1658  ! local variables
1659  real(dp), allocatable :: x_out(:) ! local x_out
1660  real(dp), allocatable :: y_out(:) ! local y_out
1661 
1662  ! initialize ierr
1663  ierr = 0
1664 
1665  ierr = order_per_fun_1(xy(:,1),xy(:,2),x_out,y_out,overlap,tol)
1666  chckerr('')
1667 
1668  allocate(xy_out(size(x_out),2))
1669  xy_out(:,1) = x_out
1670  xy_out(:,2) = y_out
1671  end function order_per_fun_2
1672 
1674  integer function spline_real(x,y,xnew,ynew,ord,deriv,bcs,bcs_val,extrap) &
1675  &result(ierr)
1676 
1677  use ezspline_obj
1678  use ezspline
1679 
1680  character(*), parameter :: rout_name = 'spline_real'
1681 
1682  ! input / output
1683  real(dp), intent(in), target :: x(:)
1684  real(dp), intent(in) :: y(:)
1685  real(dp), intent(in), target :: xnew(:)
1686  real(dp), intent(out) :: ynew(:)
1687  integer, intent(in), optional :: ord
1688  integer, intent(in), optional :: deriv
1689  integer, intent(in), optional :: bcs(2)
1690  real(dp), intent(in), optional :: bcs_val(2)
1691  logical, intent(in), optional :: extrap
1692 
1693  ! local variables
1694  type(ezspline1_r8) :: f_spl ! spline object
1695  integer :: ord_loc ! local order
1696  integer :: deriv_loc ! local deriv
1697  integer :: bcs_loc(2) ! local bcs
1698  integer :: kd ! counter
1699  integer :: n ! size of x, y
1700  integer :: nnew ! size of output x, y
1701  integer :: nnew_interp ! size of interpolated x, y, rest is extrapolated
1702  integer :: il(2) ! limits of interp., what falls outside is extrapolated
1703  real(dp), pointer :: x_loc(:) ! x or -x
1704  real(dp), pointer :: xnew_loc(:) ! xnew or -xnew
1705  real(dp) :: bcs_val_loc(2) ! local bcs_val
1706  real(dp) :: lim_vals(0:3) ! limit values at boundaries for extrapolation
1707  real(dp), allocatable :: xnew_ez(:) ! xnew in EZ spline doubles
1708  real(dp), allocatable :: ynew_ez(:) ! ynew in EZ spline doubles
1709  character(len=max_str_ln) :: err_msg ! error message
1710  logical :: flip_x ! whether x axis is flipped
1711  logical :: extrap_loc ! local extrap
1712 
1713  ! initialize ierr
1714  ierr = 0
1715 
1716  ! set local variables
1717  ord_loc = 3
1718  if (present(ord)) ord_loc = ord
1719  deriv_loc = 0
1720  if (present(deriv)) deriv_loc = deriv
1721  bcs_loc = [0,0]
1722  if (present(bcs)) bcs_loc = bcs
1723  if (present(bcs_val)) bcs_val_loc = bcs_val
1724  extrap_loc = .false.
1725  if (present(extrap)) extrap_loc = extrap
1726 
1727  ! set other variables
1728  n = size(x)
1729  nnew = size(xnew)
1730 
1731  ! set local x and xnew, possibly flipped if negative
1732  flip_x = (x(2) .lt. x(1))
1733  if (flip_x) then
1734  allocate(x_loc(n))
1735  allocate(xnew_loc(nnew))
1736  x_loc = -x
1737  xnew_loc = -xnew
1738  else
1739  x_loc => x
1740  xnew_loc => xnew
1741  end if
1742 
1743 #if ldebug
1744  ! check monotony
1745  do kd = 2,n
1746  if (x_loc(kd).le.x_loc(kd-1)) then
1747  ierr = 1
1748  err_msg = '|x| is not monotonously increasing'
1749  chckerr(err_msg)
1750  end if
1751  end do
1752 
1753  ! check array sizes
1754  if (size(y).ne.n) then
1755  ierr = 1
1756  err_msg = 'x and y need to have the same size'
1757  chckerr(err_msg)
1758  end if
1759  if (size(ynew).ne.nnew) then
1760  ierr = 1
1761  err_msg = 'xnew and ynew need to have the same size'
1762  chckerr(err_msg)
1763  end if
1764 
1765  ! check order
1766  if (ord_loc.lt.1 .or. ord_loc.gt.3) then
1767  ierr = 1
1768  err_msg = 'Only order 1 (linear), 2 (akima hermite) or 3 (cubic) &
1769  &possible'
1770  chckerr(err_msg)
1771  end if
1772 
1773  ! check derviative
1774  select case (deriv_loc)
1775  case (:-1)
1776  ierr = 1
1777  err_msg = 'only nonnegative degrees of derivative are possible'
1778  chckerr(err_msg)
1779  case (0:1)
1780  ! do nothing
1781  case (2:3)
1782  if (ord_loc.eq.1 .or. ord_loc.eq.2) then
1783  ierr = 1
1784  err_msg = 'Derivative of degree '//trim(i2str(deriv_loc))//&
1785  &' not possible for order '//trim(i2str(ord_loc))
1786  chckerr(err_msg)
1787  end if
1788  case (4:)
1789  ierr = 1
1790  err_msg = 'maximum degree of derivative is 2 for order 4'
1791  chckerr(err_msg)
1792  end select
1793 
1794  ! check boundary conditions
1795  select case (ord_loc)
1796  case (1)
1797  if (present(bcs) .or. present(bcs_val)) then
1798  call writo('for order 1, no boundary conditions can be &
1799  &prescribed',alert=.true.)
1800  end if
1801  case (2:3)
1802  ! check possibility
1803  do kd = 1,2
1804  if (bcs_loc(kd).lt.-1 .or. bcs_loc(kd).gt.ord_loc-1) then
1805  ierr = 1
1806  err_msg = 'For order '//trim(i2str(ord_loc))//&
1807  &' bcs has to be -1..'//trim(i2str(ord_loc-1))
1808  chckerr(err_msg)
1809  end if
1810  end do
1811 
1812  ! check whether derivatives are provided
1813  if ((any(bcs_loc.eq.1) .or. any(bcs_loc.eq.2)) &
1814  &.and. .not.present(bcs_val)) then
1815  ierr = 1
1816  err_msg = 'When prescribing first or second derviatives, &
1817  &need to provide bcs_val'
1818  chckerr(err_msg)
1819  end if
1820  end select
1821 #endif
1822 
1823  ! set up interpolation limits
1824  do kd = 1,nnew
1825  if (xnew_loc(kd).ge.x_loc(1)) exit
1826  end do
1827  il(1) = kd
1828  do kd = nnew,1,-1
1829  if (xnew_loc(kd).le.x_loc(n)) exit
1830  end do
1831  il(2) = kd
1832  nnew_interp = il(2)-il(1)+1
1833 
1834  ! check for extrapolation
1835  if ((il(1).ne.1 .or. il(2).ne.nnew) .and. &
1836  &.not.extrap_loc) then
1837  ierr = 1
1838  call writo('xnew = ['//trim(r2str(minval(xnew)))//'..'//&
1839  &trim(r2str(maxval(xnew)))//']')
1840  call writo('but x = ['//trim(r2str(minval(x)))//'..'//&
1841  &trim(r2str(maxval(x)))//']')
1842  err_msg = 'Extrapolation needed, but not allowed'
1843  chckerr(err_msg)
1844  end if
1845 
1846  ! initialize
1847  if (ord_loc.eq.1) then
1848  call ezlinear_init(f_spl,n,ierr)
1849  call ezspline_error(ierr)
1850  chckerr('')
1851  else
1852  call ezspline_init(f_spl,n,bcs_loc,ierr)
1853  call ezspline_error(ierr)
1854  chckerr('')
1855  if (ord_loc.eq.2) f_spl%isHermite = 1
1856  end if
1857 
1858  ! set grid
1859  f_spl%x1 = x_loc
1860 
1861  ! set boundary condition
1862  if (present(bcs_val)) then
1863  f_spl%bcval1min = bcs_val_loc(1)
1864  f_spl%bcval1max = bcs_val_loc(2)
1865  if (flip_x) then ! if x-axis flipped, also first derivative flipped
1866  if (bcs_loc(1) .eq. 1) f_spl%bcval1min = -f_spl%bcval1min
1867  if (bcs_loc(2) .eq. 1) f_spl%bcval1max = -f_spl%bcval1max
1868  end if
1869  end if
1870 
1871  ! set up
1872  call ezspline_setup(f_spl,y,ierr,exact_dim=.true.) ! match exact dimensions
1873  call ezspline_error(ierr)
1874  chckerr('')
1875  ! interpolated part
1876  if (il(1).le.il(2)) then
1877  allocate(xnew_ez(nnew_interp))
1878  allocate(ynew_ez(nnew_interp))
1879  xnew_ez = xnew_loc(il(1):il(2))
1880 
1881  if (deriv_loc.eq.0) then
1882  ! interpolate
1883  call ezspline_interp(f_spl,nnew_interp,xnew_ez,ynew_ez,ierr)
1884  call ezspline_error(ierr)
1885  chckerr('')
1886  else
1887  call ezspline_derivative(f_spl,deriv_loc,nnew_interp,xnew_ez,&
1888  &ynew_ez,ierr)
1889  call ezspline_error(ierr)
1890  chckerr('')
1891  end if
1892 
1893  ynew(il(1):il(2)) = ynew_ez
1894  deallocate(xnew_ez,ynew_ez)
1895  end if
1896 
1897  ! extrapolated part
1898  if (extrap_loc) then
1899  ! extrapolate left
1900  if (il(1).gt.1) then
1901  ! set up limit values at first point
1902  ierr = setup_lim_vals(x_loc(1),lim_vals)
1903  chckerr('')
1904 
1905  ! calculate extrapolation
1906  call calc_extrap(xnew_loc(1:il(1)-1),x_loc(1),lim_vals,&
1907  &ynew(1:il(1)-1))
1908  end if
1909 
1910  ! extrapolate right
1911  if (il(2).lt.nnew) then
1912  ! set up limit values at last point
1913  ierr = setup_lim_vals(x_loc(n),lim_vals)
1914  chckerr('')
1915 
1916  ! calculate extrapolation
1917  call calc_extrap(xnew_loc(il(2)+1:nnew),x_loc(n),lim_vals,&
1918  &ynew(il(2)+1:nnew))
1919  end if
1920  end if
1921 
1922  ! free
1923  call ezspline_free(f_spl,ierr)
1924  chckerr('')
1925  call ezspline_error(ierr)
1926  chckerr('')
1927  if (flip_x) then
1928  deallocate(x_loc)
1929  deallocate(xnew_loc)
1930  end if
1931  nullify(x_loc)
1932  nullify(xnew_loc)
1933  contains
1937  integer function setup_lim_vals(xb,lim_vals) result(ierr)
1938  character(*), parameter :: rout_name = 'setup_lim_vals'
1939 
1940  ! input / output
1941  real(dp), intent(in) :: xb ! x of boundary
1942  real(dp), intent(out) :: lim_vals(0:3) ! limit values
1943 
1944  ! local variables
1945  integer :: kd ! counter
1946  integer :: max_deriv ! maximum degree of derivative
1947 
1948  ! initialize ierr
1949  ierr = 0
1950 
1951  ! initialize
1952  lim_vals = 0._dp
1953  select case(ord_loc)
1954  case (1:2)
1955  max_deriv = 1
1956  case (3)
1957  max_deriv = 3
1958  end select
1959 
1960  ! interpolation
1961  call ezspline_interp(f_spl,xb,lim_vals(0),ierr)
1962  call ezspline_error(ierr)
1963  chckerr('')
1964 
1965  ! derivatives
1966  do kd = 1,max_deriv
1967  call ezspline_derivative(f_spl,kd,xb,lim_vals(kd),ierr)
1968  call ezspline_error(ierr)
1969  chckerr('')
1970  end do
1971  end function setup_lim_vals
1972 
1974  subroutine calc_extrap(xnew,xb,lim_vals,ynew)
1975  ! input / output
1976  real(dp), intent(in) :: xnew(:) ! xnew where to extrapolate
1977  real(dp), intent(in) :: xb ! x of boundary
1978  real(dp), intent(in) :: lim_vals(0:2) ! limit values
1979  real(dp), intent(out) :: ynew(:) ! y at xnew
1980 
1981  ! local variables
1982  real(dp), allocatable :: xdel(:) ! delta x for extrapolation
1983 
1984  allocate(xdel(size(xnew)))
1985  xdel = xnew-xb
1986 
1987  select case (deriv_loc)
1988  case (0)
1989  ynew = lim_vals(0) + &
1990  &lim_vals(1)*xdel + &
1991  &lim_vals(2)*xdel**2*0.5_dp
1992  case (1)
1993  ynew = lim_vals(1) + &
1994  &lim_vals(2)*xdel
1995  case (2)
1996  ynew = lim_vals(2)
1997  end select
1998  end subroutine calc_extrap
1999  end function spline_real
2001  integer function spline_complex(x,y,xnew,ynew,ord,deriv,bcs,bcs_val,&
2002  &extrap) result(ierr)
2003 
2004  use ezspline_obj
2005  use ezspline
2006 
2007  character(*), parameter :: rout_name = 'spline_complex'
2008 
2009  ! input / output
2010  real(dp), intent(in) :: x(:)
2011  complex(dp), intent(in) :: y(:)
2012  real(dp), intent(in) :: xnew(:)
2013  complex(dp), intent(out) :: ynew(:)
2014  integer, intent(in), optional :: ord
2015  integer, intent(in), optional :: deriv
2016  integer, intent(in), optional :: bcs(2)
2017  complex(dp), intent(in), optional :: bcs_val(2)
2018  logical, intent(in), optional :: extrap
2019 
2020  ! local variables
2021  real(dp), allocatable :: ynew_loc(:,:) ! local ynew
2022 
2023  ! set up local variables
2024  allocate(ynew_loc(size(ynew),2))
2025 
2026  ! call real version for real part
2027  if (present(bcs_val)) then
2028  ierr = spline_real(x,rp(y),xnew,ynew_loc(:,1),ord,deriv,bcs,&
2029  &bcs_val=rp(bcs_val),extrap=extrap)
2030  chckerr('')
2031  else
2032  ierr = spline_real(x,rp(y),xnew,ynew_loc(:,1),ord,deriv,bcs,&
2033  &extrap=extrap)
2034  chckerr('')
2035  end if
2036 
2037  ! call real version for complex part
2038  if (present(bcs_val)) then
2039  ierr = spline_real(x,ip(y),xnew,ynew_loc(:,2),ord,deriv,bcs,&
2040  &bcs_val=ip(bcs_val),extrap=extrap)
2041  chckerr('')
2042  else
2043  ierr = spline_real(x,ip(y),xnew,ynew_loc(:,2),ord,deriv,bcs,&
2044  &extrap=extrap)
2045  chckerr('')
2046  end if
2047 
2048  ! save
2049  ynew = ynew_loc(:,1) + iu*ynew_loc(:,2)
2050  end function spline_complex
2051 
2062  subroutine calc_aux_utilities(n_mod)
2063  use num_vars, only: max_deriv
2064 
2065  ! input / output
2066  integer, intent(in), optional :: n_mod
2067 
2068  ! local variables
2069  integer :: id, jd, kd ! counters
2070 
2071  ! common for all program stles
2072 
2073  ! derivatives d from 0 to max_deriv+1
2074  allocate(d(0:max_deriv+1,0:max_deriv+1,0:max_deriv+1))
2075  do kd = 0,max_deriv+1
2076  do jd = 0,max_deriv+1
2077  do id = 0,max_deriv+1
2078  if (id+jd+kd.le.max_deriv+1) then ! valid derivatives
2079  d(id,jd,kd) = calc_derivs_1d_id([id,jd,kd],3)
2080  else ! derivatives too high
2081  d(id,jd,kd) = 0
2082  end if
2083  end do
2084  end do
2085  end do
2086 
2087  ! metrics m from 1 to 3
2088  allocate(m(1:3,1:3))
2089  do jd = 1,3
2090  do id = 1,3
2091  m(id,jd) = calc_derivs_1d_id([id,jd],2)
2092  end do
2093  end do
2094 
2095  ! Fourier modes from 1 to n_mod
2096  if (present(n_mod)) then
2097  allocate(f(1:n_mod,1:n_mod))
2098  do jd = 1,n_mod
2099  do id = 1,n_mod
2100  f(id,jd) = calc_derivs_1d_id([id,jd],2)
2101  end do
2102  end do
2103  end if
2104  end subroutine calc_aux_utilities
2105 
2145  integer function calc_derivs_1d_id(deriv,dims) result(res)
2146  ! input / output
2147  integer, intent(in) :: deriv(:)
2148  integer, intent(in) :: dims
2149 
2150  ! local variables
2151  integer :: id, jd ! counters
2152  integer :: id_max ! index of last nonzero element in deriv
2153  integer :: tot_deriv ! total derivative
2154  integer, allocatable :: deriv_loc(:) ! local extended deriv
2155 
2156  ! set up local deriv
2157  allocate(deriv_loc(0:size(deriv)))
2158  deriv_loc(0) = 0
2159  deriv_loc(1:size(deriv)) = deriv
2160 
2161  ! set up total derivative
2162  tot_deriv = sum(deriv)
2163 
2164  ! initialize res
2165  res = 1
2166 
2167  ! calculate displacement if total deriv not equal to zero (assuming
2168  ! deriv > 0)
2169  if (tot_deriv.gt.0) then
2170  ! get index of last nonzero element in deriv
2171  id_max = 0
2172  do id = 1,size(deriv)
2173  if (deriv(id).gt.0) id_max = id
2174  end do
2175 
2176  ! iterate over dimensions
2177  do jd = 0,id_max-1
2178  ! iterate over derivatives of dimension jd (where deriv(jd) = 0)
2179  do id = 0,tot_deriv-sum(deriv_loc(0:jd))-1
2180  res = res + fac(dims-jd+id-1)/(fac(id)*fac(dims-jd-1))
2181  end do
2182  end do
2183  end if
2184  end function calc_derivs_1d_id
2185 
2203  recursive subroutine calc_derivs_of_order(res,order,dims)
2204  ! input / output
2205  integer, intent(inout), allocatable :: res(:,:)
2206  integer, intent(in) :: order
2207  integer, intent(in), optional :: dims
2208 
2209  ! local variables
2210  integer :: id ! counter
2211  integer, allocatable :: derivs_loc(:,:) ! derivs of local subblock
2212  integer :: id_tot ! counter for total dimension of local derivs.
2213  integer :: dims_loc ! local version of dims
2214 
2215  ! set up local dims
2216  dims_loc = 3
2217  if (present(dims)) dims_loc = dims
2218 
2219  ! allocate resulting derivs
2220  if (allocated(res)) deallocate(res)
2221  allocate(res(&
2222  &dims_loc,fac(order+dims_loc-1)/(fac(order)*fac(dims_loc-1))))
2223 
2224  ! iterate over the first index if order greater than 0
2225  if (order.gt.0) then
2226  id_tot = 0
2227  do id = order,0,-1
2228  ! calculate the derivs of local subblock if more than 1
2229  ! dimension
2230  if (dims_loc.gt.1) then
2231  call calc_derivs_of_order(derivs_loc,order-id,dims_loc-1) ! calculate iteratively subblock
2232  res(1,id_tot+1:id_tot+size(derivs_loc,2)) = id ! set first dimension
2233  res(2:dims_loc,id_tot+1:id_tot+size(derivs_loc,2)) &
2234  &= derivs_loc ! set next dimensions
2235  id_tot = id_tot+size(derivs_loc,2)
2236  else
2237  res = order
2238  end if
2239  end do
2240  else
2241  res = 0
2242  end if
2243  end subroutine calc_derivs_of_order
2245  function derivs(order,dims) result(res)
2246  ! input / output
2247  integer, intent(in) :: order
2248  integer, intent(in), optional :: dims
2249  integer, allocatable :: res(:,:)
2250 
2251  ! call the subroutine
2252  call calc_derivs_of_order(res,order,dims)
2253  end function derivs
2254 
2259  integer function check_deriv(deriv,max_deriv,sr_name) result(ierr)
2260  character(*), parameter :: rout_name = 'check_deriv'
2261 
2262  ! input / output
2263  integer, intent(in) :: deriv(3)
2264  integer, intent(in) :: max_deriv
2265  character(len=*), intent(in) :: sr_name
2266 
2267  ! local variables
2268  integer :: id
2269  character(len=max_str_ln) :: err_msg ! error message
2270 
2271  ! initialize ierr
2272  ierr = 0
2273 
2274  ! test the derivatives
2275  do id = 1, 3
2276  if (deriv(id).gt.max_deriv) then
2277  err_msg = 'Asking '//trim(sr_name)//' for a derivative&
2278  & in the '//trim(i2str(id))//'th dimension of order '&
2279  &//trim(i2str(deriv(id)))//', while the maximum order is '&
2280  &//trim(i2str(max_deriv))
2281  ierr = 1
2282  chckerr(err_msg)
2283  end if
2284  end do
2285  end function
2286 
2323  integer function calc_ext_var(ext_var,var,var_points,ext_point,deriv_in) &
2324  &result(ierr)
2325  character(*), parameter :: rout_name = 'ext_var'
2326 
2327  ! input / output
2328  real(dp), intent(inout) :: ext_var
2329  real(dp), intent(in) :: var(:)
2330  real(dp), intent(in) :: var_points(:)
2331  real(dp), intent(in) :: ext_point
2332  integer, intent(in), optional :: deriv_in
2333 
2334  ! local variables
2335  integer :: id, jd
2336  integer :: istat
2337  integer :: pol_deg
2338  integer :: deriv
2339  real(dp), allocatable :: a(:,:), lu(:,:)
2340  real(dp), allocatable :: a_i(:)
2341  real(dp) :: fact ! multiplicative factor, see below
2342  character(len=max_str_ln) :: err_msg ! error message
2343 
2344  ! initialize ierr
2345  ierr = 0
2346 
2347  ! determine the degree of the polynomial ext_var
2348  pol_deg = size(var)
2349 
2350  ! tests
2351  if (size(var).ne.size(var_points)) then
2352  err_msg = 'The size of the arrays containing the variable and &
2353  &the points do not match'
2354  ierr = 1
2355  chckerr(err_msg)
2356  end if
2357  if (present(deriv_in)) then
2358  deriv = deriv_in
2359  if (deriv.ge.pol_deg) then ! order of derivative too hight
2360  err_msg = 'Asking ext_var for '//trim(i2str(deriv_in))&
2361  &//'th derivative, while using a polynomial of degree '&
2362  &//trim(i2str(pol_deg))
2363  ierr = 1
2364  chckerr(err_msg)
2365  else if (deriv.lt.0) then
2366  err_msg = 'Asking ext_var for a derivative of order '&
2367  &//trim(i2str(deriv_in))
2368  ierr = 1
2369  chckerr(err_msg)
2370  end if
2371  else
2372  deriv = 0
2373  end if
2374 
2375  ! fill matrix
2376  allocate(a(pol_deg,pol_deg))
2377  allocate(lu(pol_deg,pol_deg))
2378  a(:,1) = 1.0_dp
2379  col: do jd = 2, pol_deg
2380  row: do id = 1, pol_deg
2381  a(id,jd) = var_points(id)**(jd-1)
2382  end do row
2383  end do col
2384 
2385  ! solve the system for a_i
2386  allocate(a_i(pol_deg))
2387  a_i = var
2388  lu = a
2389  call dgesv( pol_deg, 1, lu, pol_deg, [(id,id=1,pol_deg)], &
2390  &a_i, pol_deg, istat)
2391 
2392 
2393  if (istat.ne.0) then
2394  err_msg = 'The solution could not be found. Are the extrapolating &
2395  &points independent?'
2396  ierr = 1
2397  if (ierr.ne.0) then
2398  call writo('The matrix equation to be solved was Ax=b, A = ')
2399  call print_ar_2(a)
2400  call writo('and b = ')
2401  call print_ar_1(var)
2402  end if
2403  chckerr(err_msg)
2404  end if
2405 
2406  ext_var = 0.0_dp
2407  ! for deriv of order i:
2408  ! d^i/dx^i f = sum_j = i ^ pol_deg ( j*(j-1)*...*(j-i) a_j x^(j-i) )
2409  polynom: do id = 1+deriv, pol_deg
2410  ! construct the factor
2411  fact = 1.0_dp
2412  d_dx: do jd = 1,deriv
2413  fact = fact * float((id-1)-(jd-1))
2414  end do d_dx
2415  ! update ext_var with current polynomial term a_i(id)
2416  ext_var = ext_var + fact*a_i(id)*ext_point**(id-1-deriv)
2417  end do polynom
2418  end function calc_ext_var
2419 
2448  integer function calc_coeff_fin_diff(deriv,nr,ind,coeff) result(ierr)
2449  character(*), parameter :: rout_name = 'calc_coeff_fin_diff'
2450 
2451  ! input / output
2452  integer, intent(in) :: deriv
2453  integer, intent(in) :: nr
2454  integer, intent(in) :: ind
2455  real(dp), intent(inout), allocatable :: coeff(:)
2456 
2457  ! local variables
2458  character(len=max_str_ln) :: err_msg ! error message
2459  real(dp), allocatable :: mat_loc(:) ! elements of local Vandermonde matrix
2460  real(dp), allocatable :: rhs_loc(:) ! local right-hand side
2461  integer :: id ! counter
2462 
2463  ! initialize ierr
2464  ierr = 0
2465 
2466  ! test whether number of points nr and position ind positive and
2467  ! consistent
2468  if (deriv.le.0) then
2469  ierr = 1
2470  err_msg = 'The degree of the derivative has to be positive'
2471  chckerr(err_msg)
2472  end if
2473  if (nr.le.0) then
2474  ierr = 1
2475  err_msg = 'The number of points in the finite differences has to &
2476  &be positive'
2477  chckerr(err_msg)
2478  end if
2479  if (ind.lt.1 .or. ind.gt.nr) then
2480  ierr = 1
2481  err_msg = 'The position of the derivative is not in range.'
2482  chckerr(err_msg)
2483  end if
2484 
2485  ! test whether number of points enough for derivative
2486  if (deriv.gt.nr-1) then
2487  ierr = 1
2488  err_msg = 'For derivatives of degree '//trim(i2str(deriv))//&
2489  &' minimally '//trim(i2str(deriv+1))//' points are necessary'
2490  chckerr(err_msg)
2491  end if
2492 
2493  ! set up output
2494  if (allocated(coeff)) deallocate(coeff)
2495  allocate(coeff(nr))
2496  allocate(mat_loc(nr))
2497  allocate(rhs_loc(nr))
2498 
2499  ! set up matrix and rhs
2500  do id = 1,nr
2501  mat_loc(id) = (id-ind)*1._dp
2502  end do
2503  rhs_loc = 0._dp
2504  rhs_loc(deriv+1) = 1._dp ! looking for this derivative
2505  call solve_vand(nr,mat_loc,rhs_loc,coeff,transp=.true.)
2506  coeff = coeff*fac(deriv)
2507  end function calc_coeff_fin_diff
2508 
2555  integer function c(ij,sym,n,lim_n)
2556  ! input / output
2557  integer, intent(in) :: ij(2)
2558  logical, intent(in) :: sym
2559  integer, intent(in), optional :: n
2560  integer, intent(in), optional :: lim_n(2,2)
2561 
2562  ! local variables
2563  integer :: n_loc ! local n
2564  integer :: js ! j* = min(0,j,min(1)-min(2)+1)
2565  integer :: ij_loc(2) ! local ij
2566  integer :: kd ! dummy integer
2567 
2568  ! set local n
2569  n_loc = 3
2570  if (present(n)) n_loc = n
2571 
2572  ! set c depending on symmetry
2573  if (sym) then ! symmetric
2574  ! set local ij, refering to the total indices
2575  ij_loc = ij
2576  if (present(lim_n)) ij_loc = ij + lim_n(1,:) - 1 ! displace with minimum indices
2577 
2578  ! swap indices if upper diagonal
2579  if (ij_loc(2).gt.ij_loc(1)) then
2580  kd = ij_loc(2)
2581  ij_loc(2) = ij_loc(1)
2582  ij_loc(1) = kd
2583  end if
2584 
2585  ! get c for whole matrix
2586  c = nint((ij_loc(2)-1)*n_loc + ij_loc(1) - &
2587  &(ij_loc(2)-1)*ij_loc(2)*0.5)
2588 
2589  if (present(lim_n)) then
2590  ! subtract if submatrix
2591  js = max(0,min(ij_loc(2)-lim_n(1,2)+1,lim_n(1,1)-lim_n(1,2)+1))
2592  c = c - nint((lim_n(1,2)-1)*(n_loc+1-lim_n(1,2)*0.5) &
2593  &+js*(lim_n(1,1)-lim_n(1,2)+0.5-js*0.5) &
2594  &+(n_loc-lim_n(2,1))*(ij_loc(2)-lim_n(1,2)))
2595  end if
2596  else ! asymmetric
2597  if (present(lim_n)) then
2598  ! set c with local size
2599  c = (ij(2)-1)*(lim_n(2,1)-lim_n(1,1)+1) + ij(1)
2600  else
2601  ! set c with total size
2602  c = (ij(2)-1)*n_loc + ij(1)
2603  end if
2604  end if
2605  end function c
2606 
2608  recursive function fac(n) result(fact)
2609  ! input / output
2610  integer, intent(in) :: n
2611  integer :: fact
2612 
2613  ! calculate recursively
2614  if (n.eq.0) then
2615  fact = 1
2616  else
2617  fact = n * fac(n-1)
2618  end if
2619  end function fac
2620 
2625  integer function is_sym(n,nn,sym) result(ierr)
2626  character(*), parameter :: rout_name = 'is_sym'
2627 
2628  ! input / output
2629  integer, intent(in) :: n
2630  integer, intent(in) :: nn
2631  logical, intent(inout) :: sym
2632 
2633  ! local variables
2634  character(len=max_str_ln) :: err_msg ! error message
2635 
2636  ! initialize ierr
2637  ierr = 0
2638 
2639  ! set sym
2640  if (nn.eq.n**2) then
2641  sym = .false.
2642  else if (nn.eq.n*(n+1)/2) then
2643  sym = .true.
2644  else
2645  ierr = 1
2646  err_msg = 'The matrix corresponds neither to a symmetric nor to a &
2647  &normal matrix'
2648  chckerr('')
2649  end if
2650  end function is_sym
2651 
2656  recursive function lcm(u, v) result(res)
2657  integer, intent(in) :: u
2658  integer, intent(in) :: v
2659  integer :: res
2660 
2661  res = u*v / gcd(u,v)
2662  end function lcm
2663 
2668  recursive function gcd(u, v) result(res)
2669  integer, intent(in) :: u
2670  integer, intent(in) :: v
2671  integer :: res
2672 
2673  if (mod(u, v) /= 0) then
2674  res = gcd(v, mod(u, v))
2675  else
2676  res = v
2677  end if
2678  end function gcd
2679 
2695  subroutine shift_f(Al,Bl,Cl,A,B,C)
2696  ! input / output
2697  integer, intent(in) :: al(2)
2698  integer, intent(in) :: bl(2)
2699  integer, intent(in) :: cl(2)
2700  real(dp), intent(in) :: a(al(1):al(2),2)
2701  real(dp), intent(in) :: b(bl(1):bl(2),2)
2702  real(dp), intent(inout) :: c(cl(1):cl(2),2)
2703 
2704  ! local variables
2705  integer :: i_a, i_b, i_c ! indices in A, B and C
2706 
2707  ! initialize C
2708  c = 0._dp
2709 
2710  ! loop over A
2711  do i_a = al(1),al(2)
2712  do i_b = bl(1),bl(2)
2713  ! contribution to i_A + i_B
2714  i_c = i_a + i_b
2715 
2716  if (i_c.ge.cl(1) .and. i_c.le.cl(2)) then
2717  c(i_c,1) = c(i_c,1) + 0.5* a(i_a,1)*b(i_b,1)
2718  c(i_c,2) = c(i_c,2) + 0.5* a(i_a,1)*b(i_b,2)
2719  c(i_c,2) = c(i_c,2) + 0.5* a(i_a,2)*b(i_b,1)
2720  c(i_c,1) = c(i_c,1) - 0.5* a(i_a,2)*b(i_b,2)
2721  end if
2722 
2723  ! contribution to i_A + i_B
2724  i_c = i_a - i_b
2725 
2726  if (i_c.ge.cl(1) .and. i_c.le.cl(2)) then
2727  c(i_c,1) = c(i_c,1) + 0.5* a(i_a,1)*b(i_b,1)
2728  c(i_c,2) = c(i_c,2) - 0.5* a(i_a,1)*b(i_b,2)
2729  c(i_c,2) = c(i_c,2) + 0.5* a(i_a,2)*b(i_b,1)
2730  c(i_c,1) = c(i_c,1) + 0.5* a(i_a,2)*b(i_b,2)
2731  end if
2732  end do
2733  end do
2734  end subroutine shift_f
2735 
2750  subroutine solve_vand(n,a,b,x,transp)
2751  ! input / output
2752  integer, intent(in) :: n
2753  real(dp), intent(in) :: a(n)
2754  real(dp), intent(in) :: b(n)
2755  real(dp), intent(inout) :: x(n)
2756  logical, intent(in), optional :: transp
2757 
2758  ! local variables
2759  integer :: jd, kd ! counters
2760  logical :: transp_loc ! local transp
2761 
2762  ! initialize
2763  transp_loc = .false.
2764  if (present(transp)) transp_loc = transp
2765  x(1:n) = b(1:n)
2766 
2767  if (transp_loc) then
2768  do kd = 1, n - 1
2769  do jd = n, kd + 1, -1
2770  x(jd) = x(jd) - a(kd) * x(jd-1)
2771  end do
2772  end do
2773 
2774  do kd = n - 1, 1, -1
2775  do jd = kd + 1, n
2776  x(jd) = x(jd) / ( a(jd) - a(jd-kd) )
2777  end do
2778  do jd = kd, n - 1
2779  x(jd) = x(jd) - x(jd+1)
2780  end do
2781  end do
2782  else
2783  do kd = 1, n - 1
2784  do jd = n, kd + 1, -1
2785  x(jd) = ( x(jd) - x(jd-1) ) / ( a(jd) - a(jd-kd) )
2786  end do
2787  end do
2788 
2789  do kd = n - 1, 1, -1
2790  do jd = kd, n - 1
2791  x(jd) = x(jd) - a(kd) * x(jd+1)
2792  end do
2793  end do
2794  end if
2795  end subroutine solve_vand
2796 
2797 #if ldebug
2798 
2828  integer function calc_d2_smooth(fil_N,x,y,D2y,style) result(ierr)
2829  character(*), parameter :: rout_name = 'calc_D2_smooth'
2830 
2831  ! input / output
2832  integer, intent(in) :: fil_n
2833  real(dp), intent(in) :: x(:)
2834  real(dp), intent(in) :: y(:)
2835  real(dp), intent(inout) :: d2y(:)
2836  integer, intent(in) :: style
2837 
2838  ! local variables
2839  integer :: n ! size of arrays
2840  integer :: id, kd ! counters
2841  integer :: sym_m ! center of symmetry
2842  real(dp), allocatable :: s(:) ! coeffcients s
2843  real(dp), allocatable :: ab(:) ! alpha (style 1) or beta (style 2)
2844  character(len=max_str_ln) :: err_msg ! error message
2845 
2846  ! initialize ierr
2847  ierr = 0
2848 
2849  ! tests
2850  if (fil_n.lt.5) then
2851  ierr = 1
2852  err_msg = 'filter length must be at least 5'
2853  chckerr(err_msg)
2854  end if
2855  if (mod(fil_n,2).ne.1) then
2856  ierr = 1
2857  err_msg = 'filter length must be odd'
2858  chckerr(err_msg)
2859  end if
2860 
2861  ! set size of arrays
2862  sym_m = (fil_n-1)/2
2863  n = size(x)
2864  allocate(s(0:sym_m+1))
2865  allocate(ab(sym_m))
2866 
2867  ! calculate coefficients s
2868  s(sym_m+1) = 0._dp
2869  s(sym_m) = 1._dp
2870  do kd = sym_m-1,0,-1
2871  s(kd) = ((2*fil_n-10)*s(kd+1) - (fil_n+2*kd+3)*s(kd+2)) / &
2872  &(fil_n-2*kd-1)
2873  end do
2874 
2875  ! calculate for all n-fil_N interior (style 1) or left (style 2) points
2876  ! The second derivative at other points are set to zero.
2877  select case (style)
2878  case (1) ! central
2879  d2y(1:sym_m) = 0._dp
2880  d2y(n-sym_m+1:n) = 0._dp
2881  do id = sym_m+1,n-sym_m
2882  ! calculate alpha
2883  do kd = 1,sym_m
2884  ab(kd) = 4*kd**2*s(kd)*(x(id+kd)-x(id-kd))**(-2)
2885  end do
2886  ! set D2y
2887  d2y(id) = sum(ab*(y(id+1:id+sym_m))) + &
2888  &sum(ab*y(id-1:id-sym_m:-1)) - &
2889  &2*y(id)*sum(ab)
2890  end do
2891  case (2) ! left
2892  d2y(1:2*sym_m) = 0._dp
2893  do id = 2*sym_m+1,n
2894  ! calculate beta
2895  do kd = 1,sym_m
2896  ab(kd) = 4*kd**2*s(kd)*(x(id+kd-sym_m)-x(id-kd-sym_m))**(-2)
2897  end do
2898  ! set D2y
2899  d2y(id) = sum(ab*(y(id+1-sym_m:id))) + &
2900  &sum(ab*y(id-1-sym_m:id-2*sym_m:-1)) - &
2901  &2*y(id-sym_m)*sum(ab)
2902  end do
2903  case default
2904  err_msg = 'No style associated with '//trim(i2str(style))
2905  ierr = 1
2906  chckerr(err_msg)
2907  end select
2908 
2909  ! rescale
2910  d2y = d2y * 2._dp**(3._dp-fil_n)
2911  end function calc_d2_smooth
2912 #endif
2913 end module num_utilities
num_utilities::c
integer function, public c(ij, sym, n, lim_n)
Convert 2-D coordinates (i,j) to the storage convention used in matrices.
Definition: num_utilities.f90:2556
num_utilities::gcd
recursive integer function, public gcd(u, v)
Returns least denominator using the GCD.
Definition: num_utilities.f90:2669
num_utilities::calc_inv_0d
integer function calc_inv_0d(inv_0D, A)
private constant version
Definition: num_utilities.f90:772
num_utilities::calc_int
Integrates a function using the trapezoidal rule.
Definition: num_utilities.f90:160
num_utilities::m
integer, dimension(:,:), allocatable, public m
1-D array indices of metric indices
Definition: num_utilities.f90:24
num_utilities::debug_con2dis_reg
logical, public debug_con2dis_reg
plot debug information for con2dis_reg()
Definition: num_utilities.f90:27
num_utilities::round_with_tol
Rounds an arry of values to limits, with a tolerance that can optionally be modified.
Definition: num_utilities.f90:169
num_vars::dp
integer, parameter, public dp
double precision
Definition: num_vars.f90:46
num_utilities::con
Either takes the complex conjugate of a square matrix element A, defined on a 3-D grid,...
Definition: num_utilities.f90:221
num_utilities::calc_ext_var
integer function, public calc_ext_var(ext_var, var, var_points, ext_point, deriv_in)
Extrapolates a function.
Definition: num_utilities.f90:2325
num_vars
Numerical variables used by most other modules.
Definition: num_vars.f90:4
num_utilities::f
integer, dimension(:,:), allocatable, public f
1-D array indices of Fourier mode combination indices
Definition: num_utilities.f90:25
num_vars::max_str_ln
integer, parameter, public max_str_ln
maximum length of strings
Definition: num_vars.f90:50
num_utilities::con2dis
Convert between points from a continuous grid to a discrete grid.
Definition: num_utilities.f90:188
num_utilities::calc_aux_utilities
subroutine, public calc_aux_utilities(n_mod)
Initialize utilities for fast future reference, depending on program style.
Definition: num_utilities.f90:2063
num_utilities::calc_d2_smooth
integer function, public calc_d2_smooth(fil_N, x, y, D2y, style)
Calculate second derivative with smoothing formula by Holoborodko, .
Definition: num_utilities.f90:2829
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_utilities::order_per_fun
Order a periodic function to include and an overlap.
Definition: num_utilities.f90:248
str_utilities
Operations on strings.
Definition: str_utilities.f90:4
num_utilities::solve_vand
subroutine, public solve_vand(n, a, b, x, transp)
Solve a Vandermonde system .
Definition: num_utilities.f90:2751
num_utilities::calc_mult
Calculate matrix multiplication of two square matrices .
Definition: num_utilities.f90:95
num_utilities::calc_coeff_fin_diff
integer function, public calc_coeff_fin_diff(deriv, nr, ind, coeff)
Calculate the coefficients for finite differences.
Definition: num_utilities.f90:2449
num_utilities::spline
Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the m...
Definition: num_utilities.f90:276
num_utilities::calc_inv
Calculate inverse of square matrix A.
Definition: num_utilities.f90:81
num_utilities::shift_f
subroutine, public shift_f(Al, Bl, Cl, A, B, C)
Calculate multiplication through shifting of fourier modes A and B into C.
Definition: num_utilities.f90:2696
num_utilities::derivs
integer function, dimension(:,:), allocatable, public derivs(order, dims)
Returns derivatives of certain order.
Definition: num_utilities.f90:2246
num_utilities::lcm
recursive integer function, public lcm(u, v)
Returns common multiple using the Euclid's algorithm.
Definition: num_utilities.f90:2657
num_utilities::conv_mat
Converts a (symmetric) matrix A with the storage convention described in eq_vars.eq_2_type.
Definition: num_utilities.f90:123
num_utilities::d
integer, dimension(:,:,:), allocatable, public d
1-D array indices of derivatives
Definition: num_utilities.f90:23
num_utilities::check_deriv
integer function, public check_deriv(deriv, max_deriv, sr_name)
checks whether the derivatives requested for a certain subroutine are valid
Definition: num_utilities.f90:2260
num_vars::tol_zero
real(dp), public tol_zero
tolerance for zeros
Definition: num_vars.f90:133
num_utilities::dis2con
Convert between points from a discrete grid to a continuous grid.
Definition: num_utilities.f90:206
num_vars::max_deriv
integer, parameter, public max_deriv
highest derivatives for metric factors in Flux coords.
Definition: num_vars.f90:52
num_utilities
Numerical utilities.
Definition: num_utilities.f90:4
num_utilities::is_sym
integer function, public is_sym(n, nn, sym)
Determines whether a matrix making use of the storage convention in eq_vars.eq_2_type is symmetric or...
Definition: num_utilities.f90:2626
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_utilities::calc_derivs_1d_id
integer function calc_derivs_1d_id(deriv, dims)
Returns the 1D indices for derivatives of a certain order in certain dimensions.
Definition: num_utilities.f90:2146
num_utilities::fac
recursive integer function, public fac(n)
Calculate factorial.
Definition: num_utilities.f90:2609
num_utilities::bubble_sort
Sorting with the bubble sort routine.
Definition: num_utilities.f90:237
num_utilities::add_arr_mult
Add to an array (3) the product of arrays (1) and (2).
Definition: num_utilities.f90:39
setup_lim_vals
integer function setup_lim_vals(xb, lim_vals)
/private set up limit values
Definition: num_utilities.f90:1938
num_utilities::debug_calc_coeff_fin_diff
logical, public debug_calc_coeff_fin_diff
plot debug information for calc_coeff_fin_diff()
Definition: num_utilities.f90:28
num_utilities::calc_det
Calculate determinant of a matrix.
Definition: num_utilities.f90:63
this documentation was created with Doxygen 1.8.17 by Toon Weyens