Octopus
math.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22module math_oct_m
23 use debug_oct_m
24 use global_oct_m
25 use, intrinsic :: iso_fortran_env
30
31 implicit none
32
33 private
34 public :: &
35 is_close, &
37 ylmr_cmplx, &
38 ylmr_real, &
39 weights, &
40 factorial, &
41 hermite, &
46 ddelta, &
47 member, &
51 even, &
52 odd, &
56 pad, &
57 pad_pow2, &
58 log2, &
60 phi1, &
61 phi2, &
63 is_prime, &
75
76 interface is_close
77 module procedure dis_close_scalar, zis_close_scalar
78 end interface
79
80 interface diagonal_matrix
82 end interface diagonal_matrix
83
87
91
94 interface interpolate
97 end interface interpolate
98
99 interface symmetrize_matrix
100 module procedure dsymmetrize_matrix, zsymmetrize_matrix
101 end interface symmetrize_matrix
102
103 interface log2
104 module procedure dlog2, ilog2, llog2
105 end interface log2
106
107 interface pad
108 module procedure pad4, pad8, pad48, pad88
109 end interface pad
110
111contains
112
131 elemental logical function dis_close_scalar(x, y, rtol, atol)
132 real(real64), intent(in) :: x, y
133 real(real64), optional, intent(in) :: rtol
134 real(real64), optional, intent(in) :: atol
135 real(real64) :: atol_, rtol_
136
137 atol_ = optional_default(atol, 1.e-8_real64)
138 rtol_ = optional_default(rtol, 1.e-5_real64)
139
140 dis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
141 end function dis_close_scalar
142
144 elemental logical function zis_close_scalar(x, y, rtol, atol)
145 complex(real64), intent(in) :: x, y
146 real(real64), optional, intent(in) :: rtol
147 real(real64), optional, intent(in) :: atol
148 real(real64) :: atol_, rtol_
149
150 atol_ = optional_default(atol, 1.e-8_real64)
151 rtol_ = optional_default(rtol, 1.e-5_real64)
152
153 zis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
154 end function zis_close_scalar
155
156
157 ! ---------------------------------------------------------
160 pure function idiagonal_matrix(dim, diag) result(matrix)
161 integer, intent(in) :: dim
162 integer, intent(in) :: diag
163 integer :: matrix(dim, dim)
164
165 integer :: ii
166
167 matrix = 0
168 do ii = 1, dim
169 matrix(ii, ii) = diag
170 end do
171
172 end function idiagonal_matrix
174 ! ---------------------------------------------------------
175 recursive function hermite(n, x) result (h)
176 integer, intent(in) :: n
177 real(real64), intent(in) :: x
178
179 real(real64) :: h
180
181 ! no push_sub, called too frequently
182
183 if (n <= 0) then
184 h = m_one
185 elseif (n == 1) then
186 h = m_two*x
187 else
188 h = m_two*x*hermite(n-1,x) - m_two*(n-1)*hermite(n-2,x)
189 end if
190
191 end function hermite
193
194 ! ---------------------------------------------------------
195 recursive function factorial (n) RESULT (fac)
196 integer, intent(in) :: n
197
198 integer :: fac
199
200 ! no push_sub for recursive function
201
202 if (n <= 1) then
203 fac = 1
204 else
205 fac = n * factorial(n-1)
206 end if
207 end function factorial
208
209 ! ---------------------------------------------------------
211 subroutine ylmr_cmplx(xx, li, mi, ylm)
212 real(real64), intent(in) :: xx(3)
213 integer, intent(in) :: li, mi
214 complex(real64), intent(out) :: ylm
215
216 integer :: i
217 real(real64) :: dx, dy, dz, r, plm, cosm, sinm, cosmm1, sinmm1, cosphi, sinphi
218 real(real64), parameter :: tiny = 1.e-15_real64
219
220 ! no push_sub: this routine is called too many times
221
222 ! if l=0, no calculations are required
223 if (li == 0) then
224 ylm = cmplx(sqrt(m_one/(m_four*m_pi)), m_zero, real64)
225 return
226 end if
227
228 r = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
229
230 ! if r=0, direction is undefined => make ylm=0 except for l=0
231 if (r <= tiny) then
232 ylm = m_z0
233 return
234 end if
235 dx = xx(1)/r
236 dy = xx(2)/r
237 dz = xx(3)/r
238
239 ! get the associated Legendre polynomial (including the normalization factor)
240 ! We need to avoid getting outside [-1:1] due to rounding errors
241 if(abs(dz) > m_one) dz = sign(m_one, dz)
242 plm = loct_legendre_sphplm(li, abs(mi), dz)
243
244 ! compute sin(|m|*phi) and cos(|m|*phi)
245 r = sqrt(dx*dx + dy*dy)
246 if(r > tiny) then
247 cosphi = dx/r
248 sinphi = dy/r
249 else ! In this case the values are ill defined so we choose \phi=0
250 cosphi = m_zero
251 sinphi = m_zero
252 end if
254 cosm = m_one
255 sinm = m_zero
256 do i = 1, abs(mi)
257 cosmm1 = cosm
258 sinmm1 = sinm
259 cosm = cosmm1*cosphi - sinmm1*sinphi
260 sinm = cosmm1*sinphi + sinmm1*cosphi
261 end do
262
263 !And now ylm
264 ylm = plm*cmplx(cosm, sinm, real64)
265
266 if (mi < 0) then
267 ylm = conjg(ylm)
268 do i = 1, abs(mi)
269 ylm = -ylm
270 end do
271 end if
272
273 end subroutine ylmr_cmplx
274
275
276 ! ---------------------------------------------------------
282 subroutine ylmr_real(xx, li, mi, ylm)
283 real(real64), intent(in) :: xx(3)
284 integer, intent(in) :: li, mi
285 real(real64), intent(out) :: ylm
286
287 integer, parameter :: lmaxd = 20
288 real(real64), parameter :: tiny = 1.e-15_real64
289 integer :: i, ilm0, l, m, mabs
290 integer, save :: lmax = -1
291
292 real(real64) :: cmi, cosm, cosmm1, cosphi, dphase, dplg, fac, &
293 fourpi, plgndr, phase, pll, pmm, pmmp1, sinm, &
294 sinmm1, sinphi, rsize, rx, ry, rz, xysize
295 real(real64), save :: c(0:(lmaxd+1)*(lmaxd+1))
296
297 ! no push_sub: called too frequently
298
299 ! evaluate normalization constants once and for all
300 if (li > lmax) then
301 fourpi = m_four*m_pi
302 do l = 0, li
303 ilm0 = l*l + l
304 do m = 0, l
305 fac = (2*l+1)/fourpi
306 do i = l - m + 1, l + m
307 fac = fac/i
308 end do
309 c(ilm0 + m) = sqrt(fac)
310 ! next line because ylm are real combinations of m and -m
311 if (m /= 0) c(ilm0 + m) = c(ilm0 + m)*sqrt(m_two)
312 c(ilm0 - m) = c(ilm0 + m)
313 end do
314 end do
315 lmax = li
316 end if
317
318 ! if l=0, no calculations are required
319 if (li == 0) then
320 ylm = c(0)
321 return
322 end if
323
324 ! if r=0, direction is undefined => make ylm=0 except for l=0
325 rsize = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
326 if (rsize < tiny) then
327 ylm = m_zero
328 return
329 end if
330
331 rx = xx(1)/rsize
332 ry = xx(2)/rsize
333 rz = xx(3)/rsize
334
335 ! explicit formulas for l=1 and l=2
336 if (li == 1) then
337 select case (mi)
338 case (-1)
339 ylm = (-c(1))*ry
340 case (0)
341 ylm = c(2)*rz
342 case (1)
343 ylm = (-c(3))*rx
344 end select
345 return
346 end if
347
348 if (li == 2) then
349 select case (mi)
350 case (-2)
351 ylm = c(4)*6.0_real64*rx*ry
352 case (-1)
353 ylm = (-c(5))*m_three*ry*rz
354 case (0)
355 ylm = c(6)*m_half*(m_three*rz*rz - m_one)
356 case (1)
357 ylm = (-c(7))*m_three*rx*rz
358 case (2)
359 ylm = c(8)*m_three*(rx*rx - ry*ry)
360 end select
361 return
362 end if
363
364 ! general algorithm based on routine plgndr of numerical recipes
365 mabs = abs(mi)
366 xysize = sqrt(max(rx*rx + ry*ry, tiny))
367 cosphi = rx/xysize
368 sinphi = ry/xysize
369 cosm = m_one
370 sinm = m_zero
371 do m = 1, mabs
372 cosmm1 = cosm
373 sinmm1 = sinm
374 cosm = cosmm1*cosphi - sinmm1*sinphi
375 sinm = cosmm1*sinphi + sinmm1*cosphi
376 end do
377
378 if (mi < 0) then
379 phase = sinm
380 dphase = mabs*cosm
381 else
382 phase = cosm
383 dphase = (-mabs)*sinm
384 end if
385
386 pmm = m_one
387 fac = m_one
388
389 if (mabs > m_zero) then
390 do i = 1, mabs
391 pmm = (-pmm)*fac*xysize
392 fac = fac + m_two
393 end do
394 end if
395
396 if (li == mabs) then
397 plgndr = pmm
398 dplg = (-li)*rz*pmm/(xysize**2)
399 else
400 pmmp1 = rz*(2*mabs + 1)*pmm
401 if (li == mabs + 1) then
402 plgndr = pmmp1
403 dplg = -((li*rz*pmmp1 - (mabs + li)*pmm)/(xysize**2))
404 else
405 do l = mabs + 2, li
406 pll = (rz*(2*l - 1)*pmmp1 - (l + mabs - 1)*pmm)/(l - mabs)
407 pmm = pmmp1
408 pmmp1 = pll
409 end do
410 plgndr = pll
411 dplg = -((li*rz*pll - (l + mabs - 1)*pmm)/(xysize**2))
412 end if
413 end if
414
415 ilm0 = li*li + li
416 cmi = c(ilm0 + mi)
417 ylm = cmi*plgndr*phase
418
419 end subroutine ylmr_real
420
421
422 ! ---------------------------------------------------------
432 subroutine weights(N, M, cc, side)
433 integer, intent(in) :: n, m
434 real(real64), intent(out) :: cc(0:,0:,0:)
435 integer, optional, intent(in) :: side
436
437 integer :: i, j, k, mn, side_
438 real(real64) :: c1, c2, c3, c4, c5, xi
439 real(real64), allocatable :: x(:)
440
441 push_sub(weights)
442
443 safe_allocate(x(0:m))
444
445 if (present(side)) then
446 side_ = side
447 else
448 side_ = 0
449 end if
450
451 select case (side_)
452 case (-1)
453 ! grid-points for left-side finite-difference formulas on an equi-spaced grid
454 mn = m
455 x(:) = (/(-i,i=0,mn)/)
456 case (+1)
457 ! grid-points for right-side finite-difference formulas on an equi-spaced grid
458 mn = m
459 x(:) = (/(i,i=0,mn)/)
460 case default
461 ! grid-points for centered finite-difference formulas on an equi-spaced grid
462 mn = m/2
463 x(:) = (/0,(-i,i,i=1,mn)/)
464 end select
465
466
467
468 xi = m_zero ! point at which the approx. are to be accurate
469
470 cc = m_zero
471 cc(0,0,0) = m_one
472
473 c1 = m_one
474 c4 = x(0) - xi
475
476 do j = 1, m
477 mn = min(j,n)
478 c2 = m_one
479 c5 = c4
480 c4 = x(j) - xi
481
482 do k = 0, j - 1
483 c3 = x(j) - x(k)
484 c2 = c2*c3
485
486 if (j <= n) cc(k, j - 1, j) = m_zero
487 cc(k, j, 0) = c4*cc(k, j - 1, 0)/c3
488
489 do i = 1, mn
490 cc(k, j, i) = (c4*cc(k, j - 1, i) - i*cc(k, j - 1, i - 1))/c3
491 end do
492
493 cc(j, j, 0) = -c1*c5*cc(j - 1, j - 1, 0) / c2
494 end do
495
496 do i = 1, mn
497 cc(j, j, i) = c1*(i*cc(j - 1, j - 1, i - 1) - c5*cc(j - 1, j - 1, i))/c2
498 end do
499
500 c1 = c2
501 end do
502
503 safe_deallocate_a(x)
504 pop_sub(weights)
505 end subroutine weights
506
507 ! ---------------------------------------------------------
508 real(real64) pure function ddelta(i, j)
509 integer, intent(in) :: i
510 integer, intent(in) :: j
511
512 ! no push_sub in pure function
513
514 if (i == j) then
515 ddelta = m_one
516 else
517 ddelta = m_zero
518 end if
519
520 end function ddelta
521
522
523 ! ---------------------------------------------------------
526 subroutine make_idx_set(n, out, length, in)
527 integer, intent(in) :: n
528 integer, allocatable, intent(out) :: out(:)
529 integer, intent(out) :: length
530 integer, optional, intent(in) :: in(:)
531
532 integer :: i
533
534 push_sub(make_idx_set)
535
536 if (present(in)) then
537 length = ubound(in, 1)
538 safe_allocate(out(1:length))
539 out(:) = in(:)
540 else
541 length = n
542 safe_allocate(out(1:length))
543 do i = 1, length
544 out(i) = i
545 end do
546 end if
547
548 pop_sub(make_idx_set)
549 end subroutine make_idx_set
550
551
552 ! ---------------------------------------------------------
555 logical function member(n, a)
556 integer, intent(in) :: n
557 integer, intent(in) :: a(:)
558
559 integer :: i
560
561 ! no push_sub, called too frequently
562
563 member = .false.
564
565 do i = 1, ubound(a, 1)
566 if (a(i) == n) then
567 member = .true.
568 exit
569 end if
570 end do
571
572 end function member
573
574
575 ! ---------------------------------------------------------
576 subroutine interpolation_coefficients(nn, xa, xx, cc)
577 integer, intent(in) :: nn
578 real(real64), intent(in) :: xa(:)
579 real(real64), intent(in) :: xx
580 real(real64), intent(out) :: cc(:)
581
582 integer :: ii, kk
583
584 ! no push_sub, called too frequently
585
586 do ii = 1, nn
587 cc(ii) = m_one
588 do kk = 1, nn
589 if (kk == ii) cycle
590 cc(ii) = cc(ii)*(xx - xa(kk))/(xa(ii) - xa(kk))
591 end do
592 end do
593
594 end subroutine interpolation_coefficients
595
596
597 ! ---------------------------------------------------------
599 logical pure function even(n)
600 integer, intent(in) :: n
602 even = (mod(n, 2) == 0)
603
604 end function even
605
606
607 ! ---------------------------------------------------------
609 logical pure function odd(n)
610 integer, intent(in) :: n
611
612 odd = .not. even(n)
613
614 end function odd
615
616 ! ---------------------------------------------------------
620 real(real64), intent(in) :: x(:)
621 real(real64), intent(out) :: u(:)
622
623 integer :: n, k, j
624 real(real64) :: sumx2
625 real(real64), allocatable :: xx(:)
626
628
629 n = size(x)
630 assert(n>1)
631 assert(size(u) == n-1)
632
633 ! These lines make the code less machine-dependent.
634 safe_allocate(xx(1:n))
635 do j = 1, n
636 if (abs(x(j))<1.0e-8_real64) then
637 xx(j) = m_zero
638 else
639 xx(j) = x(j)
640 end if
641 end do
642
643 u = m_zero
644 do k = 1, n-1
645 sumx2 = m_zero
646 do j = k+1, n
647 sumx2 = sumx2 + xx(j)**2
648 end do
649 if (abs(sumx2) <= m_epsilon .and. abs(xx(k)) <= m_epsilon) exit
650 if (k < n - 1) then
651 u(k) = atan2(sqrt(sumx2), xx(k))
652 else
653 u(n-1) = atan2(xx(n), xx(n-1))
654 end if
655 end do
656
657
659 end subroutine cartesian2hyperspherical
660 ! ---------------------------------------------------------
661
662
663 ! ---------------------------------------------------------
666 subroutine hyperspherical2cartesian(u, x)
667 real(real64), intent(in) :: u(:)
668 real(real64), intent(out) :: x(:)
670 integer :: n, j, k
671
673
674 n = size(x)
675 assert(n>1)
676 assert(size(u) == n-1)
677
678 if (n == 2) then
679 x(1) = cos(u(1))
680 x(2) = sin(u(1))
681 elseif (n == 3) then
682 x(1) = cos(u(1))
683 x(2) = sin(u(1))*cos(u(2))
684 x(3) = sin(u(1))*sin(u(2))
685 else
686 x(1) = cos(u(1))
687 x(2) = sin(u(1))*cos(u(2))
688 x(3) = sin(u(1))*sin(u(2))*cos(u(3))
689 do j = 4, n - 1
690 x(j) = m_one
691 do k = 1, j - 1
692 x(j) = x(j) * sin(u(k))
693 end do
694 x(j) = x(j) * cos(u(j))
695 end do
696 x(n) = m_one
697 do k = 1, n - 2
698 x(n) = x(n) * sin(u(k))
699 end do
700 x(n) = x(n) * sin(u(n-1))
701 end if
704 end subroutine hyperspherical2cartesian
705 ! ---------------------------------------------------------
706
707
708 ! ---------------------------------------------------------
712 subroutine hypersphere_grad_matrix(grad_matrix, r, x)
713 real(real64), intent(out) :: grad_matrix(:,:)
714 real(real64), intent(in) :: r
715 real(real64), intent(in) :: x(:)
716
717 integer :: n, l, m
718
720
721 n = size(x)+1 ! the dimension of the matrix is (n-1)x(n)
722
723 ! --- l=1 ---
724 grad_matrix = m_one
725 grad_matrix(1,1) = -r*sin(x(1))
726 do m = 2, n-1
727 grad_matrix(m,1) = m_zero
728 end do
729
730 ! --- l=2..(n-1) ---
731 do l=2, n-1
732 do m=1, n-1
733 if (m == l) then
734 grad_matrix(m,l) = -r*grad_matrix(m,l)*product(sin(x(1:l)))
735 elseif (m < l) then
736 grad_matrix(m,l) = grad_matrix(m,l)*r*cos(x(m))*cos(x(l))*product(sin(x(1:m-1)))*product(sin(x(m+1:l-1)))
737 else
738 grad_matrix(m,l) = m_zero
739 end if
740 end do
741 end do
742
743 ! --- l=n ---
744 do m=1, n-1
745 grad_matrix(m,n) = r*cos(x(m))*grad_matrix(m,n)*product(sin(x(1:m-1)))*product(sin(x(m+1:n-1)))
746 end do
747
749 end subroutine hypersphere_grad_matrix
750
751 ! ---------------------------------------------------------
752 integer(int64) pure function pad88(size, blk)
753 integer(int64), intent(in) :: size
754 integer(int64), intent(in) :: blk
755
756 integer(int64) :: mm
757
758 mm = mod(size, blk)
759 if (mm == 0) then
760 pad88 = size
761 else
762 pad88 = size + blk - mm
763 end if
764 end function pad88
765
766 ! ---------------------------------------------------------
767 integer(int64) pure function pad48(size, blk)
768 integer, intent(in) :: size
769 integer(int64), intent(in) :: blk
770
771 pad48 = pad88(int(size, int64), blk)
772 end function pad48
773
774 ! ---------------------------------------------------------
775 integer(int64) pure function pad8(size, blk)
776 integer(int64), intent(in) :: size
777 integer, intent(in) :: blk
778
779 pad8 = pad88(size, int(blk, int64))
780 end function pad8
781
782 ! ---------------------------------------------------------
783 integer pure function pad4(size, blk)
784 integer, intent(in) :: size
785 integer, intent(in) :: blk
786
787 pad4 = int(pad88(int(size, int64), int(blk, int64)), int32)
788 end function pad4
789
790 ! ---------------------------------------------------------
791
796 integer pure function pad_pow2(size)
797 integer, intent(in) :: size
798
799 integer :: mm, mm2
800
801 mm = size
802 pad_pow2 = 1
803
804 ! loop below never terminates otherwise! just pick 1 as value.
805 if (size == 0) return
806
807 ! first we divide by two all the times we can, so we catch exactly
808 ! the case when size is already a power of two
809 do
810 mm2 = mm/2
811 if (mm - 2*mm2 /= 0) exit
813 mm = mm2
814 end do
815
816 ! the rest is handled by log2
817 if (mm /= 1) then
818 pad_pow2 = pad_pow2*2**(ceiling(log(real(mm, real64) )/log(2.0_8)))
819 end if
820
821 end function pad_pow2
822
823 ! -------------------------------------------------------
824
825 real(real64) pure function dlog2(xx)
826 real(real64), intent(in) :: xx
827
828 dlog2 = log(xx)/log(m_two)
829 end function dlog2
830
831 ! -------------------------------------------------------
832
833 integer pure function ilog2(xx)
834 integer, intent(in) :: xx
835
836 ilog2 = nint(log2(real(xx, real64) ))
837 end function ilog2
838
839 ! -------------------------------------------------------
840
841 integer(int64) pure function llog2(xx)
842 integer(int64), intent(in) :: xx
843
844 llog2 = nint(log2(real(xx, real64) ), kind=int64)
845 end function llog2
846
847 ! -------------------------------------------------------
849 complex(real64) pure function exponential(z)
850 complex(real64), intent(in) :: z
851
852 exponential = exp(z)
853 end function exponential
854
855 ! -------------------------------------------------------
861 complex(real64) pure function phi1(z)
862 complex(real64), intent(in) :: z
863
864 real(real64), parameter :: cut = 0.002_real64
865
866 if (abs(z) < cut) then
867 phi1 = m_one + m_half*z*(m_one + z/m_three*(m_one + m_fourth*z*(m_one + z/m_five)))
868 else
869 phi1 = (exp(z) - m_z1)/z
870 end if
871 end function phi1
872
873 ! -------------------------------------------------------
880 complex(real64) pure function phi2(z)
881 complex(real64), intent(in) :: z
882
883 real(real64), parameter :: cut = 0.002_real64
884
885 if (abs(z) < cut) then
886 phi2 = m_half + z/6.0_real64*(m_one + m_fourth*z*(m_one + z/m_five*(m_one + z/6.0_real64)))
887 else
888 phi2 = (exp(z) - z - m_z1)/z**2
889 end if
890 end function phi2
891
892 ! -------------------------------------------------------
894 real(real64) pure function square_root(x)
895 real(real64), intent(in) :: x
896
897 square_root = sqrt(x)
898 end function square_root
899
900 ! -------------------------------------------------------
901
902 logical function is_prime(n)
903 integer, intent(in) :: n
904
905 integer :: i, root
906
907 push_sub(is_prime)
908
909 if (n < 1) then
910 message(1) = "Internal error: is_prime does not take negative numbers."
911 call messages_fatal(1)
912 end if
913 if (n == 1) then
914 is_prime = .false.
915 pop_sub(is_prime)
916 return
917 end if
919 root = nint(sqrt(real(n, real64) ))
920 do i = 2, root
921 if (mod(n,i) == 0) then
922 is_prime = .false.
923 pop_sub(is_prime)
924 return
925 end if
926 end do
927
928 is_prime = .true.
929 pop_sub(is_prime)
930 end function is_prime
931
932 ! ---------------------------------------------------------
937 subroutine generate_rotation_matrix(R, ff, tt)
938 real(real64), intent(out) :: r(:,:)
939 real(real64), intent(in) :: ff(:)
940 real(real64), intent(in) :: tt(:)
941
942 integer :: dim, i, j
943 real(real64) :: th, uv, uu, vv, ft
944 real(real64), allocatable :: axis(:), u(:), v(:), f(:), t(:), p(:)
945
947
948 dim = size(ff,1)
949
950 assert((dim < 3) .or. (dim > 2))
951 assert(size(tt,1) == dim)
952 assert((size(r,1) == dim) .and. (size(r,2) == dim))
953
954 safe_allocate(u(1:dim))
955 safe_allocate(v(1:dim))
956 safe_allocate(f(1:dim))
957 safe_allocate(t(1:dim))
958
959
960 !normalize
961 f = ff / norm2(ff)
962 t = tt / norm2(tt)
963
964 ft = dot_product(f,t)
965
966 if (abs(ft) < m_one) then
967 select case (dim)
968 case (2)
969 th = acos(ft)
970 r(1,1) = cos(th)
971 r(1,2) = -sin(th)
972
973 r(2,1) = sin(th)
974 r(2,2) = cos(th)
975
976 case (3)
977 if (.false.) then
978 !Old implementation
979 safe_allocate(axis(1:dim))
980 th = acos(ft)
981
982 u = f / dot_product(f,f)
983 v = t /dot_product(t,t)
984
985 axis = dcross_product(u,v)
986 axis = axis / norm2(axis)
988 r(1,1) = cos(th) + axis(1)**2 * (1 - cos(th))
989 r(1,2) = axis(1)*axis(2)*(1-cos(th)) + axis(3)*sin(th)
990 r(1,3) = axis(1)*axis(3)*(1-cos(th)) - axis(2)*sin(th)
991
992 r(2,1) = axis(2)*axis(1)*(1-cos(th)) - axis(3)*sin(th)
993 r(2,2) = cos(th) + axis(2)**2 * (1 - cos(th))
994 r(2,3) = axis(2)*axis(3)*(1-cos(th)) + axis(1)*sin(th)
996 r(3,1) = axis(3)*axis(1)*(1-cos(th)) + axis(2)*sin(th)
997 r(3,2) = axis(3)*axis(2)*(1-cos(th)) - axis(1)*sin(th)
998 r(3,3) = cos(th) + axis(3)**2 * (1 - cos(th))
999
1000 safe_deallocate_a(axis)
1001 end if
1002
1003 if (.true.) then
1004 !Naive implementation
1005 th = acos(ft)
1006 u = dcross_product(f,t)
1007 u = u / norm2(u)
1008
1009 r(1,1) = u(1)**2 + (1-u(1)**2)*cos(th)
1010 r(1,2) = u(1)*u(2)*(1-cos(th)) - u(3)*sin(th)
1011 r(1,3) = u(1)*u(3) + u(2)*sin(th)
1012
1013 r(2,1) = u(1)*u(2)*(1-cos(th)) + u(3)*sin(th)
1014 r(2,2) = u(2)**2 + (1-u(2)**2)*cos(th)
1015 r(2,3) = u(2)*u(3)*(1-cos(th)) - u(1)*sin(th)
1016
1017 r(3,1) = u(1)*u(3)*(1-cos(th)) - u(2)*sin(th)
1018 r(3,2) = u(2)*u(3)*(1-cos(th)) + u(1)*sin(th)
1019 r(3,3) = u(3)**2 + (1-u(3)**2)*cos(th)
1020 end if
1021
1022 if (.false.) then
1023 !Fast
1024 safe_allocate(p(1:dim))
1025
1026 if (abs(f(1)) <= abs(f(2)) .and. abs(f(1)) < abs(f(3))) then
1027 p = (/m_one, m_zero, m_zero/)
1028 else if (abs(f(2)) < abs(f(1)) .and. abs(f(2)) <= abs(f(3))) then
1029 p = (/m_zero, m_one, m_zero/)
1030 else if (abs(f(3)) <= abs(f(1)) .and. abs(f(3)) < abs(f(2))) then
1031 p = (/m_zero, m_zero, m_one/)
1032 end if
1033
1034 u = p - f
1035 v = p - t
1036
1037 uu = dot_product(u,u)
1038 vv = dot_product(v,v)
1039 uv = dot_product(u,v)
1040
1041 do i = 1,3
1042 do j = 1,3
1043
1044 r(i,j) = ddelta(i,j) - m_two * u(i)*u(j)/uu - m_two * v(i)*v(j)/vv &
1045 + 4_real64*uv * v(i)*u(j) /(uu*vv)
1046
1047 end do
1048 end do
1049
1050 safe_deallocate_a(p)
1051 end if
1052
1053
1054 end select
1055
1056 else
1057
1058 r = m_zero
1059 do i = 1,dim
1060 r(i,i) = m_one
1061 end do
1062
1063 end if
1064
1065
1066 safe_deallocate_a(u)
1067 safe_deallocate_a(v)
1068
1070 end subroutine generate_rotation_matrix
1071
1072 ! ---------------------------------------------------------
1081 subroutine numder_ridders(x, h, res, err, f)
1082 implicit none
1083 real(real64), intent(in) :: x, h
1084 real(real64), intent(out) :: res, err
1085 interface
1086 subroutine f(x, fx)
1087 use, intrinsic :: iso_fortran_env
1088 implicit none
1089 real(real64), intent(in) :: x
1090 real(real64), intent(inout) :: fx
1091 end subroutine f
1092 end interface
1093
1094 real(real64), parameter :: con = 1.4_real64, big = 1.0e30_real64, safe = m_two
1095 integer, parameter :: ntab = 20
1096 integer :: i, j
1097 real(real64) :: errt, fac, hh, fx1, fx2
1098 real(real64), allocatable :: a(:, :)
1099
1100 push_sub(numder_ridders)
1101
1102 if (abs(h) <= m_epsilon) then
1103 message(1) = "h must be nonzero in numder_ridders"
1104 call messages_fatal(1)
1105 end if
1106
1107 safe_allocate(a(1:ntab, 1:ntab))
1108
1109 hh = h
1110 call f(x+hh, fx1)
1111 call f(x-hh, fx2)
1112 a(1,1) = (fx1-fx2) / (m_two*hh)
1113 err = big
1114 do i = 2, ntab
1115 hh = hh / con
1116 call f(x+hh, fx1)
1117 call f(x-hh, fx2)
1118 a(1,i) = (fx1-fx2) / (m_two*hh)
1119 fac = con**2
1120 do j = 2, i
1121 a(j,i) = (a(j-1,i)*fac-a(j-1,i-1)) / (fac-m_one)
1122 fac = con**2*fac
1123 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
1124 if (errt .le. err) then
1125 err = errt
1126 res = a(j,i)
1127 end if
1128 end do
1129 if (abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) return
1130 end do
1131
1132 safe_deallocate_a(a)
1133 pop_sub(numder_ridders)
1134 end subroutine numder_ridders
1135
1136 ! ---------------------------------------------------------
1137 pure function dzcross_product(a, b) result(c)
1138 real(real64), intent(in) :: a(:)
1139 complex(real64), intent(in) :: b(:)
1140
1141 complex(real64) :: c(1:3)
1142
1143 c(1) = a(2)*b(3) - a(3)*b(2)
1144 c(2) = a(3)*b(1) - a(1)*b(3)
1145 c(3) = a(1)*b(2) - a(2)*b(1)
1146
1147 end function dzcross_product
1148
1149 ! ---------------------------------------------------------
1150 pure function zdcross_product(a, b) result(c)
1151 complex(real64), intent(in) :: a(:)
1152 real(real64), intent(in) :: b(:)
1153
1154 complex(real64) :: c(1:3)
1155
1156 c(1) = a(2)*b(3) - a(3)*b(2)
1157 c(2) = a(3)*b(1) - a(1)*b(3)
1158 c(3) = a(1)*b(2) - a(2)*b(1)
1159
1160 end function zdcross_product
1161
1162
1163!*****************************************************************************80
1164!
1165!! LM_POLYNOMIAL evaluates Laguerre polynomials Lm(n,m,x).
1166!
1167! First terms:
1168!
1169! M = 0
1170!
1171! Lm(0,0,X) = 1
1172! Lm(1,0,X) = -X + 1
1173! Lm(2,0,X) = X^2 - 4 X + 2
1174! Lm(3,0,X) = -X^3 + 9 X^2 - 18 X + 6
1175! Lm(4,0,X) = X^4 - 16 X^3 + 72 X^2 - 96 X + 24
1176! Lm(5,0,X) = -X^5 + 25 X^4 - 200 X^3 + 600 X^2 - 600 x + 120
1177! Lm(6,0,X) = X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720
1178!
1179! M = 1
1180!
1181! Lm(0,1,X) = 0
1182! Lm(1,1,X) = -1,
1183! Lm(2,1,X) = 2 X - 4,
1184! Lm(3,1,X) = -3 X^2 + 18 X - 18,
1185! Lm(4,1,X) = 4 X^3 - 48 X^2 + 144 X - 96
1186!
1187! M = 2
1188!
1189! Lm(0,2,X) = 0
1190! Lm(1,2,X) = 0,
1191! Lm(2,2,X) = 2,
1192! Lm(3,2,X) = -6 X + 18,
1193! Lm(4,2,X) = 12 X^2 - 96 X + 144
1194!
1195! M = 3
1196!
1197! Lm(0,3,X) = 0
1198! Lm(1,3,X) = 0,
1199! Lm(2,3,X) = 0,
1200! Lm(3,3,X) = -6,
1201! Lm(4,3,X) = 24 X - 96
1202!
1203! M = 4
1204!
1205! Lm(0,4,X) = 0
1206! Lm(1,4,X) = 0
1207! Lm(2,4,X) = 0
1208! Lm(3,4,X) = 0
1209! Lm(4,4,X) = 24
1210!
1211! Recursion:
1212!
1213! Lm(0,M,X) = 1
1214! Lm(1,M,X) = (M+1-X)
1215!
1216! if 2 <= N:
1217!
1218! Lm(N,M,X) = ( (M+2*N-1-X) * Lm(N-1,M,X)
1219! + (1-M-N) * Lm(N-2,M,X) ) / N
1220!
1221! Special values:
1222!
1223! For M = 0, the associated Laguerre polynomials Lm(N,M,X) are equal
1224! to the Laguerre polynomials L(N,X).
1225!
1226! Licensing:
1227!
1228! This code is distributed under the GNU LGPL license.
1229!
1230! Modified:
1231!
1232! 08 February 2003
1233!
1234! Author:
1235!
1236! John Burkardt
1237!
1238! Reference:
1239!
1240! Milton Abramowitz, Irene Stegun,
1241! Handbook of Mathematical Functions,
1242! National Bureau of Standards, 1964,
1243! ISBN: 0-486-61272-4,
1244! LC: QA47.A34.
1245!
1246! Parameters:
1247!
1248! Input, integer ( kind = 4 ) MM, the number of evaluation points.
1249!
1250! Input, integer ( kind = 4 ) N, the highest order polynomial to compute.
1251! Note that polynomials 0 through N will be computed.
1252!
1253! Input, integer ( kind = 4 ) M, the parameter. M must be nonnegative.
1254!
1255! Input, real ( kind = rk ) X(MM), the evaluation points.
1256!
1257! Output, real ( kind = rk ) CX(MM,0:N), the associated Laguerre polynomials
1258! of degrees 0 through N evaluated at the evaluation points.
1259!
1260! Taken from
1261! https://people.sc.fsu.edu/~jburkardt/f_src/laguerre_polynomial/laguerre_polynomial.html
1262! and adapted to Octopus by N. Tancogne-Dejean
1263 subroutine generalized_laguerre_polynomial ( np, nn, mm, xx, cx )
1264 integer, intent(in) :: np
1265 integer, intent(in) :: nn, mm
1266 real(real64), intent(in) :: xx(np)
1267 real(real64), intent(out) :: cx(np)
1268
1269 integer ii
1270 real(real64), allocatable :: cx_tmp(:,:)
1271
1272 assert(mm >= 0)
1273
1274 safe_allocate(cx_tmp(1:np, 0:nn))
1275
1276 if (nn < 0) then
1277 cx = m_zero
1278 return
1279 end if
1280
1281 if (nn == 0) then
1282 cx(1:np) = m_one
1283 return
1284 end if
1285
1286 cx_tmp(1:np,0) = m_one
1287 cx_tmp(1:np,1) = real(mm + 1, real64) - xx(1:np)
1288
1289 do ii = 2, nn
1290 cx_tmp(1:np, ii) = &
1291 (( real(mm + 2*ii - 1, real64) - xx(1:np)) * cx_tmp(1:np,ii-1) &
1292 + real(-mm -ii + 1, real64) * cx_tmp(1:np,ii-2)) / real(ii, real64)
1293 end do
1294
1295 cx(1:np) = cx_tmp(1:np, nn)
1296
1297 safe_deallocate_a(cx_tmp)
1298
1299 end subroutine generalized_laguerre_polynomial
1300
1301 subroutine dupper_triangular_to_hermitian(nn, aa)
1302 integer, intent(in) :: nn
1303 real(real64), intent(inout) :: aa(:, :)
1304
1305 integer :: ii, jj
1306
1307 do jj = 1, nn
1308 do ii = 1, jj-1
1309 aa(jj, ii) = aa(ii, jj)
1310 end do
1311 end do
1312 end subroutine dupper_triangular_to_hermitian
1313
1314 subroutine zupper_triangular_to_hermitian(nn, aa)
1315 integer, intent(in) :: nn
1316 complex(real64), intent(inout) :: aa(:, :)
1317
1318 integer :: ii, jj
1319
1320 do jj = 1, nn
1321 do ii = 1, jj-1
1322 aa(jj, ii) = conjg(aa(ii, jj))
1323 end do
1324 aa(jj, jj) = real(aa(jj, jj), real64)
1325 end do
1326 end subroutine zupper_triangular_to_hermitian
1327
1328 subroutine dlower_triangular_to_hermitian(nn, aa)
1329 integer, intent(in) :: nn
1330 real(real64), intent(inout) :: aa(:, :)
1331
1332 integer :: ii, jj
1333
1334 do jj = 1, nn
1335 do ii = 1, jj-1
1336 aa(ii, jj) = aa(jj, ii)
1337 end do
1338 end do
1339 end subroutine dlower_triangular_to_hermitian
1340
1341 subroutine zlower_triangular_to_hermitian(nn, aa)
1342 integer, intent(in) :: nn
1343 complex(real64), intent(inout) :: aa(:, :)
1344
1345 integer :: ii, jj
1346
1347 do jj = 1, nn
1348 do ii = 1, jj-1
1349 aa(ii, jj) = conjg(aa(jj, ii))
1350 end do
1351 aa(jj, jj) = real(aa(jj, jj), real64)
1352 end do
1353 end subroutine zlower_triangular_to_hermitian
1354
1356 subroutine dsymmetrize_matrix(nn, aa)
1357 integer, intent(in) :: nn
1358 real(real64), intent(inout) :: aa(:, :)
1359
1360 integer :: ii, jj
1361
1362 do ii = 1, nn
1363 do jj = ii+1, nn
1364 aa(ii, jj) = m_half * (aa(ii, jj) + aa(jj, ii))
1365 aa(jj, ii) = aa(ii, jj)
1366 end do
1367 end do
1368 end subroutine dsymmetrize_matrix
1369
1371 subroutine zsymmetrize_matrix(n, a)
1372 integer, intent(in) :: n
1373 complex(real64), intent(inout) :: a(:, :)
1374
1375 integer :: i,j
1376
1377 do i = 1, n
1378 do j = i+1, n
1379 a(i, j) = m_half * (a(i, j) + conjg(a(j, i)))
1380 a(j, i) = conjg(a(i, j))
1381 end do
1382 end do
1383 end subroutine zsymmetrize_matrix
1384
1385
1386 subroutine dzero_small_elements_matrix(nn, aa, tol)
1387 integer, intent(in) :: nn
1388 real(real64), intent(inout) :: aa(:, :)
1389 real(real64) :: tol
1390
1391 where(abs(aa)<tol)
1392 aa = m_zero
1393 end where
1395
1397 logical function is_symmetric(a, tol)
1398 real(real64), intent(in) :: a(:, :)
1399 real(real64), optional, intent(in) :: tol
1400
1401 real(real64) :: tolerance
1402
1403 push_sub(is_symmetric)
1404
1405 tolerance = optional_default(tol, 1.0e-8_real64)
1406 ! a must be square
1407 assert(size(a, 1) == size(a, 2))
1408
1409 is_symmetric = all(abs(a - transpose(a)) < tolerance)
1410
1411 pop_sub(is_symmetric)
1412
1413 end function is_symmetric
1414
1416 logical pure function is_diagonal(dim, matrix, tol)
1417 integer, intent(in) :: dim
1418 real(real64), intent(in) :: matrix(:, :)
1419 real(real64), intent(in) :: tol
1420
1421 integer :: ii, jj
1422 real(real64) :: max_diag
1423
1424 max_diag = m_zero
1425 do ii = 1, dim
1426 max_diag = max(abs(matrix(ii, ii)), max_diag)
1427 end do
1428
1429 is_diagonal = .true.
1430 do ii = 1, dim
1431 do jj = 1, dim
1432 if (ii == jj) cycle
1433 if (abs(matrix(ii, jj)) > tol*max_diag) is_diagonal = .false.
1434 end do
1435 end do
1436
1437 end function is_diagonal
1438
1440 subroutine convert_to_base(num, base, converted)
1441 integer, intent(in) :: num
1442 integer, intent(in) :: base
1443 integer, intent(out) :: converted(:)
1444
1445 integer :: i, num_, length
1446
1447 ! no PUSH_SUB, called too often
1448
1449 converted = 0
1450 if (num > 0) then
1451 ! length of converted number
1452 length = floor(log(real(num, real64)) / log(real(base, real64))) + 1
1453 assert(size(converted) >= length)
1454
1455 num_ = num
1456 ! do conversion by repeated division
1457 do i = 1, size(converted)
1458 converted(i) = mod(num_, base)
1459 num_ = num_ / base
1460 end do
1461 end if
1462
1463 end subroutine convert_to_base
1466 subroutine convert_from_base(converted, base, num)
1467 integer, intent(in) :: converted(:)
1468 integer, intent(in) :: base
1469 integer, intent(out) :: num
1470
1471 integer :: i
1472
1473 ! no PUSH_SUB, called too often
1474
1475 num = m_zero
1476 do i = 1, size(converted)
1477 num = num + converted(i) * base**(i-1)
1478 end do
1480 end subroutine convert_from_base
1481
1482#include "undef.F90"
1483#include "complex.F90"
1484#include "math_inc.F90"
1485
1486#include "undef.F90"
1487#include "real.F90"
1488#include "math_inc.F90"
1489
1490end module math_oct_m
1491
1492!! Local Variables:
1493!! mode: f90
1494!! coding: utf-8
1495!! End:
This is the common interface to a simple-minded polynomical interpolation procedure (simple use of th...
Definition: math.F90:187
double acos(double __x) __attribute__((__nothrow__
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
double floor(double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine zsymmetrize_matrix(n, a)
make a complex matrix Hermitian
Definition: math.F90:1465
complex(real64) pure function, public phi2(z)
Compute phi2(z) = (phi1(z)-1)/z = (exp(z) - z - 1)/z^2.
Definition: math.F90:974
integer pure function ilog2(xx)
Definition: math.F90:927
logical pure function, public even(n)
Returns if n is even.
Definition: math.F90:693
pure real(real64) function, dimension(dim, dim) ddiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
Definition: math.F90:1837
subroutine, public generalized_laguerre_polynomial(np, nn, mm, xx, cx)
Definition: math.F90:1357
complex(real64) pure function, public exponential(z)
Wrapper for exponential.
Definition: math.F90:943
subroutine dinterpolate_0(xa, ya, x, y)
Definition: math.F90:1915
subroutine, public ylmr_cmplx(xx, li, mi, ylm)
Computes spherical harmonics ylm at position (x, y, z)
Definition: math.F90:305
subroutine zinterpolate_0(xa, ya, x, y)
Definition: math.F90:1724
subroutine, public convert_to_base(num, base, converted)
convert an integer in base 10 to a different base
Definition: math.F90:1534
subroutine dinterpolate_1(xa, ya, x, y)
Definition: math.F90:1887
real(real64) pure function dlog2(xx)
Definition: math.F90:919
integer(int64) pure function pad88(size, blk)
Definition: math.F90:846
real(real64) pure function, public square_root(x)
Wrapper for sqrt.
Definition: math.F90:988
logical function, public member(n, a)
Considers a(1:ubound(a, 1)) as an integer set and checks if n is a member of it.
Definition: math.F90:649
integer pure function pad4(size, blk)
Definition: math.F90:877
subroutine dlower_triangular_to_hermitian(nn, aa)
Definition: math.F90:1422
pure complex(real64) function, dimension(1:3), public dzcross_product(a, b)
Definition: math.F90:1231
subroutine, public cartesian2hyperspherical(x, u)
Performs a transformation of an n-dimensional vector from Cartesian coordinates to hyperspherical coo...
Definition: math.F90:713
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
Definition: math.F90:890
subroutine zinterpolate_1(xa, ya, x, y)
Definition: math.F90:1696
pure real(real64) function, dimension(1:3), public dcross_product(a, b)
Definition: math.F90:1940
integer(int64) pure function pad48(size, blk)
Definition: math.F90:861
subroutine dsymmetrize_matrix(nn, aa)
symmetrize a real matrix
Definition: math.F90:1450
subroutine, public ylmr_real(xx, li, mi, ylm)
This is a Numerical Recipes-based subroutine computes real spherical harmonics ylm at position (x,...
Definition: math.F90:376
subroutine dupper_triangular_to_hermitian(nn, aa)
Definition: math.F90:1395
subroutine, public interpolation_coefficients(nn, xa, xx, cc)
Definition: math.F90:670
logical function, public is_symmetric(a, tol)
Check if a 2D array is symmetric.
Definition: math.F90:1491
subroutine, public generate_rotation_matrix(R, ff, tt)
Generates a rotation matrix R to rotate a vector f to t.
Definition: math.F90:1031
elemental logical function zis_close_scalar(x, y, rtol, atol)
Same as dis_close_scalar for complex numbers.
Definition: math.F90:238
logical pure function, public odd(n)
Returns if n is odd.
Definition: math.F90:703
real(real64) pure function, public ddelta(i, j)
Definition: math.F90:602
subroutine, public make_idx_set(n, out, length, in)
Construct out(1:length) = (/1, ..., n/) if in is not present, out(1:length) = in otherwise.
Definition: math.F90:620
logical function, public is_prime(n)
Definition: math.F90:996
complex(real64) pure function, public phi1(z)
Compute phi1(z) = (exp(z)-1)/z.
Definition: math.F90:955
pure complex(real64) function, dimension(1:3), public zcross_product(a, b)
Definition: math.F90:1749
subroutine, public weights(N, M, cc, side)
Compute the weights for finite-difference calculations:
Definition: math.F90:526
subroutine zinterpolate_2(xa, ya, x, y)
Definition: math.F90:1661
pure complex(real64) function, dimension(dim, dim) zdiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
Definition: math.F90:1646
subroutine zlower_triangular_to_hermitian(nn, aa)
Definition: math.F90:1435
logical pure function, public is_diagonal(dim, matrix, tol)
Returns true is the matrix of size dim x dim is diagonal, given a relative tolerance.
Definition: math.F90:1510
subroutine zupper_triangular_to_hermitian(nn, aa)
Definition: math.F90:1408
pure integer function, dimension(dim, dim) idiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the real and comp...
Definition: math.F90:254
integer(int64) pure function pad8(size, blk)
Definition: math.F90:869
subroutine dinterpolate_2(xa, ya, x, y)
Definition: math.F90:1852
subroutine, public numder_ridders(x, h, res, err, f)
Numerical derivative (Ridder`s algorithm).
Definition: math.F90:1175
integer(int64) pure function llog2(xx)
Definition: math.F90:935
recursive real(real64) function, public hermite(n, x)
Definition: math.F90:269
elemental logical function dis_close_scalar(x, y, rtol, atol)
Are and equal within a tolerance.
Definition: math.F90:225
pure complex(real64) function, dimension(1:3), public zdcross_product(a, b)
Definition: math.F90:1244
subroutine, public convert_from_base(converted, base, num)
convert an integer to base 10 from a different base
Definition: math.F90:1560
subroutine, public hyperspherical2cartesian(u, x)
Performs the inverse transformation of cartesian2hyperspherical.
Definition: math.F90:760
subroutine, public dzero_small_elements_matrix(nn, aa, tol)
Definition: math.F90:1480
subroutine, public hypersphere_grad_matrix(grad_matrix, r, x)
Gives the hyperspherical gradient matrix, which contains the derivatives of the Cartesian coordinates...
Definition: math.F90:806
recursive integer function, public factorial(n)
Definition: math.F90:289
static double f(double w, void *p)
int true(void)