TABLE OF CONTENTS


ABINIT/dotproduct [ Functions ]

[ Top ] [ Functions ]

NAME

 dotproduct

FUNCTION

 scalar product of two vectors

INPUTS

 v1 and v2: two real(dp) vectors

OUTPUT

 scalar product of the two vectors

SIDE EFFECTS

WARNINGS

 vector size is not checked

NOTES

 I've benchmarked this to be speedier than the intrinsic dot_product even on
 big vectors. The point is that less check is performed.

 MG: FIXME: Well, optized blas1 is for sure better than what you wrote!
 Now I dont' have time to update ref files

PARENTS

 cgpr,brent

CHILDREN

SOURCE

7060 function dotproduct(nv1,nv2,v1,v2)
7061 
7062 
7063 !This section has been created automatically by the script Abilint (TD).
7064 !Do not modify the following lines by hand.
7065 #undef ABI_FUNC
7066 #define ABI_FUNC 'dotproduct'
7067 !End of the abilint section
7068 
7069  implicit none
7070 
7071 !Arguments ------------------------------------
7072 !scalars
7073  integer,intent(in) :: nv1,nv2
7074  real(dp) :: dotproduct
7075 !arrays
7076  real(dp),intent(in) :: v1(nv1,nv2),v2(nv1,nv2)
7077 
7078 !Local variables-------------------------------
7079 !scalars
7080  integer :: i,j
7081 
7082 ! *************************************************************************
7083  dotproduct=zero
7084  do j=1,nv2
7085   do i=1,nv1
7086    dotproduct=dotproduct+v1(i,j)*v2(i,j)
7087   end do
7088  end do
7089 end function dotproduct

ABINIT/m_numeric_tools [ Modules ]

[ Top ] [ Modules ]

NAME

  m_numeric_tools

FUNCTION

  This module contains basic tools for numeric computations.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG, GMR, MJV, XG, MVeithen, NH, FJ, MT, DCS, FrD, Olevano, Reining, Sottile, AL)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

 21 #if defined HAVE_CONFIG_H
 22 #include "config.h"
 23 #endif
 24 
 25 #include "abi_common.h"
 26 
 27 MODULE m_numeric_tools
 28 
 29  use defs_basis
 30  use m_errors
 31  use m_abicore
 32  use m_linalg_interfaces
 33 
 34  use m_fstrings,   only : itoa, sjoin
 35 
 36  implicit none
 37 
 38  private
 39 
 40  public :: arth                  ! Return an arithmetic progression
 41  public :: geop                  ! Return a geometric progression
 42  public :: reverse               ! Reverse a 1D array *IN PLACE*
 43  public :: set2unit              ! Set the matrix to be a unit matrix (if it is square)
 44  public :: get_trace             ! Calculate the trace of a square matrix
 45  public :: get_diag              ! Return the diagonal of a matrix as a vector
 46  public :: isdiagmat             ! True if matrix is diagonal
 47  public :: l2int                 ! convert logical data to int array
 48  public :: r2c,c2r               ! Transfer complex data stored in a real array to a complex array and vice versa
 49  public :: iseven                ! True if int is even
 50  public :: isinteger             ! True if all elements of rr differ from an integer by less than tol
 51  public :: is_zero               ! True if all elements of rr differ from zero by less than tol
 52  public :: isinside              ! True if float is inside an interval.
 53  public :: bisect                ! Given a monotonic array A and x find j such that A(j)>x>A(j+1) using bisection
 54  public :: imax_loc              ! Index of maxloc on an array returned as scalar instead of array-valued quantity
 55  public :: imin_loc              ! Index of minloc on an array returned as scalar instead of array-valued quantity
 56  public :: lfind                 ! Find the index of the first occurrence of .True. in a logical array.
 57  public :: list2blocks           ! Given a list of integers, find the number of contiguos groups of values.
 58  public :: mask2blocks           ! Find groups of .TRUE. elements in a logical mask.
 59  public :: linfit                ! Perform a linear fit, y=ax+b, of data
 60  public :: llsfit_svd            ! Linear least squares fit with SVD of an user-defined set of functions
 61  public :: polyn_interp          ! Polynomial interpolation with Nevilles"s algorithms, error estimate is reported
 62  public :: quadrature            ! Driver routine for performing quadratures in finite domains using different algorithms
 63  public :: cspint                ! Estimates the integral of a tabulated function.
 64  public :: ctrap                 ! Corrected trapezoidal integral on uniform grid of spacing hh.
 65  public :: coeffs_gausslegint    ! Compute the coefficients (supports and weights) for Gauss-Legendre integration.
 66  public :: simpson_cplx          ! Integrate a complex function via extended Simpson's rule.
 67  public :: hermitianize          ! Force a square matrix to be hermitian
 68  public :: mkherm                ! Make the complex array(2,ndim,ndim) hermitian, by adding half of it to its hermitian conjugate.
 69  public :: hermit                ! Rdefine diagonal elements of packed matrix to impose Hermiticity.
 70  public :: symmetrize            ! Force a square matrix to be symmetric
 71  public :: pack_matrix           ! Packs a matrix into hermitian format
 72  public :: print_arr             ! Print a vector/array
 73  public :: pade, dpade           ! Functions for Pade approximation (complex case)
 74  public :: newrap_step           ! Apply single step Newton-Raphson method to find root of a complex function
 75  public :: OPERATOR(.x.)         ! Cross product of two 3D vectors
 76  public :: l2norm                ! Return the length (ordinary L2 norm) of a vector
 77  public :: remove_copies         ! Find the subset of inequivalent items in a list.
 78  public :: denominator           ! Return the denominator of a rational number.
 79  public :: mincm                 ! Return the minimum common multiple of two integers.
 80  public :: continued_fract       ! Routine to calculate the continued fraction (see description).
 81  public :: cmplx_sphcart         ! Convert an array of cplx numbers from spherical to Cartesian coordinates or vice versa.
 82  public :: pfactorize            ! Factorize a number in terms of an user-specified set of prime factors.
 83  public :: isordered             ! Check the ordering of a sequence.
 84  public :: wrap2_zero_one        ! Transforms a real number in a reduced number in the interval [0,1[ where 1 is not included (tol12)
 85  public :: wrap2_pmhalf          ! Transforms a real number in areduced number in the interval ]-1/2,1/2] where -1/2 is not included (tol12)
 86  public :: interpol3d            ! Linear interpolation in 3D
 87  public :: interpol3d_indices    ! Computes the indices in a cube which are neighbors to the point to be interpolated in interpol3d
 88  public :: interpolate_denpot    ! Liner interpolation of scalar field e.g. density of potential
 89  public :: simpson_int           ! Simpson integral of a tabulated function. Returns arrays with integrated values
 90  public :: simpson               ! Simpson integral of a tabulated function. Returns scalar with the integral on the full mesh.
 91  public :: rhophi                ! Compute the phase and the module of a complex number.
 92  public :: smooth                ! Smooth data.
 93  public :: nderiv                ! Compute first or second derivative of input function y(x) on a regular grid.
 94  public :: central_finite_diff   ! Coefficients of the central differences, for several orders of accuracy.
 95  public :: uniformrandom         ! Returns a uniform random deviate between 0.0 and 1.0.
 96  public :: findmin               ! Compute the minimum of a function whose value and derivative are known at two points.
 97  public :: kramerskronig         ! check or apply the Kramers Kronig relation
 98  public :: invcb                 ! Compute a set of inverse cubic roots as fast as possible.
 99 
100  !MG FIXME: deprecated: just to avoid updating refs while refactoring.
101  public :: dotproduct
102 
103  interface arth
104    module procedure arth_int
105    module procedure arth_rdp
106  end interface arth
107 
108  interface reverse
109    module procedure arth_int
110    module procedure arth_rdp
111  end interface reverse
112 
113  interface set2unit
114    module procedure unit_matrix_int
115    module procedure unit_matrix_rdp
116    module procedure unit_matrix_cdp
117  end interface set2unit
118 
119  interface get_trace
120    module procedure get_trace_int
121    module procedure get_trace_rdp
122    module procedure get_trace_cdp
123  end interface get_trace
124 
125  interface get_diag
126    module procedure get_diag_int
127    module procedure get_diag_rdp
128    module procedure get_diag_cdp
129  end interface get_diag
130 
131  interface isdiagmat
132    module procedure isdiagmat_int
133    module procedure isdiagmat_rdp
134    !module procedure isdiagmat_cdp
135  end interface isdiagmat
136 
137  interface l2int
138    module procedure l2int_1D
139    module procedure l2int_2D
140    module procedure l2int_3D
141  end interface l2int
142 
143  interface r2c
144    module procedure rdp2cdp_1D
145    module procedure rdp2cdp_2D
146    module procedure rdp2cdp_3D
147    module procedure rdp2cdp_4D
148    module procedure rdp2cdp_5D
149    module procedure rdp2cdp_6D
150  end interface r2c
151 
152  interface c2r
153    module procedure cdp2rdp_1D
154    module procedure cdp2rdp_2D
155    module procedure cdp2rdp_3D
156    module procedure cdp2rdp_4D
157    module procedure cdp2rdp_5D
158  end interface c2r
159 
160  interface isinteger
161    module procedure is_integer_0d
162    module procedure is_integer_1d
163  end interface isinteger
164 
165  interface is_zero
166    module procedure is_zero_rdp_0d
167    module procedure is_zero_rdp_1d
168  end interface is_zero
169 
170  interface bisect
171    module procedure bisect_rdp
172    module procedure bisect_int
173  end interface bisect
174 
175  interface imax_loc
176    module procedure imax_loc_int
177    module procedure imax_loc_rdp
178  end interface imax_loc
179 
180  interface imin_loc
181    module procedure imin_loc_int
182    module procedure imin_loc_rdp
183  end interface imin_loc
184 
185  interface linfit
186    module procedure linfit_rdp
187    module procedure linfit_spc
188    module procedure linfit_dpc
189  end interface linfit
190 
191  interface hermitianize
192    module procedure hermitianize_spc
193    module procedure hermitianize_dpc
194  end interface hermitianize
195 
196  interface symmetrize
197    module procedure symmetrize_spc
198    module procedure symmetrize_dpc
199  end interface symmetrize
200 
201  interface print_arr  !TODO add prtm
202    module procedure print_arr1d_spc
203    module procedure print_arr1d_dpc
204    module procedure print_arr2d_spc
205    module procedure print_arr2d_dpc
206  end interface print_arr
207 
208  interface operator (.x.)
209    module procedure cross_product_int
210    module procedure cross_product_rdp
211  end interface
212 
213  interface l2norm
214    module procedure l2norm_rdp
215  end interface l2norm
216 
217  interface isordered
218    module procedure isordered_rdp
219  end interface isordered

m_numeric_tools/arth_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  arth_int

FUNCTION

  Returns an array of length nn containing an arithmetic progression whose
  starting value is start and whose step is step.

INPUTS

  start=initial point
  step=the increment
  nn=the number of points

OUTPUT

  arth(nn)=the progression

SOURCE

292 pure function arth_int(start,step,nn)
293 
294 
295 !This section has been created automatically by the script Abilint (TD).
296 !Do not modify the following lines by hand.
297 #undef ABI_FUNC
298 #define ABI_FUNC 'arth_int'
299 !End of the abilint section
300 
301  implicit none
302 
303 !Arguments ------------------------------------
304 !scalars
305  integer,intent(in) :: nn
306  integer,intent(in) :: start,step
307  integer :: arth_int(nn)
308 
309 !Local variables-------------------------------
310  integer :: ii
311 ! *********************************************************************
312 
313  select case (nn)
314 
315  case (1:)
316    arth_int(1)=start
317    do ii=2,nn
318     arth_int(ii)=arth_int(ii-1)+step
319    end do
320 
321  case (0)
322    RETURN
323  end select
324 
325 end function arth_int

m_numeric_tools/arth_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  arth_rdp

FUNCTION

INPUTS

OUTPUT

SOURCE

342 pure function arth_rdp(start,step,nn)
343 
344 
345 !This section has been created automatically by the script Abilint (TD).
346 !Do not modify the following lines by hand.
347 #undef ABI_FUNC
348 #define ABI_FUNC 'arth_rdp'
349 !End of the abilint section
350 
351  implicit none
352 
353 !Arguments ------------------------------------
354 !scalars
355  integer,intent(in) :: nn
356  real(dp),intent(in) :: start,step
357  real(dp) :: arth_rdp(nn)
358 
359 !Local variables-------------------------------
360  integer :: ii
361 ! *********************************************************************
362 
363  select case (nn)
364 
365  case (1:)
366   arth_rdp(1)=start
367   do ii=2,nn
368    arth_rdp(ii)=arth_rdp(ii-1)+step
369   end do
370 
371  case (0)
372   RETURN
373 
374  end select
375 
376 end function arth_rdp

m_numeric_tools/bisect_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  bisect_int

FUNCTION

  Given an array AA(1:N), and a value x, returns the index j such that AA(j)<=x<= AA(j + 1).
  AA must be monotonic, either increasing or decreasing. j=0 or
  j=N is returned to indicate that x is out of range.

INPUTS

OUTPUT

SOURCE

1865 pure function bisect_int(AA,xx) result(loc)
1866 
1867 
1868 !This section has been created automatically by the script Abilint (TD).
1869 !Do not modify the following lines by hand.
1870 #undef ABI_FUNC
1871 #define ABI_FUNC 'bisect_int'
1872 !End of the abilint section
1873 
1874  implicit none
1875 
1876 !Arguments ------------------------------------
1877 !scalars
1878  integer,intent(in) :: AA(:)
1879  integer,intent(in) :: xx
1880  integer :: loc
1881 
1882 !Local variables-------------------------------
1883  integer :: nn,jl,jm,ju
1884  logical :: ascnd
1885 ! *********************************************************************
1886 
1887  nn=SIZE(AA) ; ascnd=(AA(nn)>=AA(1))
1888  !
1889  ! === Initialize lower and upper limits ===
1890  jl=0 ; ju=nn+1
1891  do
1892   if (ju-jl<=1) EXIT
1893   jm=(ju+jl)/2  ! Compute a midpoint
1894   if (ascnd.EQV.(xx>=AA(jm))) then
1895    jl=jm ! Replace lower limit
1896   else
1897    ju=jm ! Replace upper limit
1898   end if
1899  end do
1900  !
1901  ! === Set the output, being careful with the endpoints ===
1902  if (xx==AA(1)) then
1903   loc=1
1904  else if (xx==AA(nn)) then
1905   loc=nn-1
1906  else
1907   loc=jl
1908  end if
1909 
1910 end function bisect_int

m_numeric_tools/bisect_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  bisect_rdp

FUNCTION

  Given an array AA(1:N), and a value x, returns the index j such that AA(j)<=x<= AA(j + 1).
  AA must be monotonic, either increasing or decreasing. j=0 or
  j=N is returned to indicate that x is out of range.

SOURCE

1801 pure function bisect_rdp(AA,xx) result(loc)
1802 
1803 
1804 !This section has been created automatically by the script Abilint (TD).
1805 !Do not modify the following lines by hand.
1806 #undef ABI_FUNC
1807 #define ABI_FUNC 'bisect_rdp'
1808 !End of the abilint section
1809 
1810  implicit none
1811 
1812 !Arguments ------------------------------------
1813 !scalars
1814  real(dp),intent(in) :: AA(:)
1815  real(dp),intent(in) :: xx
1816  integer :: loc
1817 
1818 !Local variables-------------------------------
1819  integer :: nn,jl,jm,ju
1820  logical :: ascnd
1821 ! *********************************************************************
1822 
1823  nn=SIZE(AA); ascnd=(AA(nn)>=AA(1))
1824  !
1825  ! === Initialize lower and upper limits ===
1826  jl=0; ju=nn+1
1827  do
1828    if (ju-jl<=1) EXIT
1829    jm=(ju+jl)/2  ! Compute a midpoint,
1830    if (ascnd.EQV.(xx>=AA(jm))) then
1831      jl=jm ! Replace lower limit
1832    else
1833      ju=jm ! Replace upper limit
1834    end if
1835  end do
1836  !
1837  ! === Set the output, being careful with the endpoints ===
1838  if (xx==AA(1)) then
1839    loc=1
1840  else if (xx==AA(nn)) then
1841    loc=nn-1
1842  else
1843    loc=jl
1844  end if
1845 
1846 end function bisect_rdp

m_numeric_tools/calculate_pade_a [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  calculate_pade_a

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_numeric_tools

CHILDREN

SOURCE

4739 subroutine calculate_pade_a(a,n,z,f)
4740 
4741 
4742 !This section has been created automatically by the script Abilint (TD).
4743 !Do not modify the following lines by hand.
4744 #undef ABI_FUNC
4745 #define ABI_FUNC 'calculate_pade_a'
4746 !End of the abilint section
4747 
4748  implicit none
4749 
4750 !Arguments ------------------------------------
4751 !scalars
4752  integer,intent(in) :: n
4753  complex(dpc),intent(in) :: z(n),f(n)
4754  complex(dpc),intent(out) :: a(n)
4755 
4756 !Local variables-------------------------------
4757 !scalars
4758  integer :: i,j
4759 !arrays
4760  complex(dpc) :: g(n,n)
4761 ! *************************************************************************
4762 
4763  g(1,1:n)=f(1:n)
4764 
4765  do i=2,n
4766    do j=i,n
4767      if (REAL(g(i-1,j))==zero.and.AIMAG(g(i-1,j))==zero) write(std_out,*) 'g_i(z_j)',i,j,g(i,j)
4768      g(i,j)=(g(i-1,i-1)-g(i-1,j)) / ((z(j)-z(i-1))*g(i-1,j))
4769      !write(std_out,*) 'g_i(z_j)',i,j,g(i,j)
4770    end do
4771  end do
4772  do i=1,n
4773    a(i)=g(i,i)
4774  end do
4775  !write(std_out,*) 'a ',a(:)
4776 
4777 end subroutine calculate_pade_a

m_numeric_tools/cdp2rdp_1D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_1D

FUNCTION

  Create a real array containing real and imaginary part starting from a complex array

INPUTS

  cc(:)=the input complex array

OUTPUT

  rr(2,:)=the real array

SOURCE

1385 pure function cdp2rdp_1D(cc) result(rr)
1386 
1387 
1388 !This section has been created automatically by the script Abilint (TD).
1389 !Do not modify the following lines by hand.
1390 #undef ABI_FUNC
1391 #define ABI_FUNC 'cdp2rdp_1D'
1392 !End of the abilint section
1393 
1394  implicit none
1395 
1396 !Arguments ------------------------------------
1397 !scalars
1398  complex(dpc),intent(in) :: cc(:)
1399  real(dp) :: rr(2,SIZE(cc))
1400 
1401 ! *********************************************************************
1402 
1403  rr(1,:)=REAL (cc(:))
1404  rr(2,:)=AIMAG(cc(:))
1405 
1406 end function cdp2rdp_1D

m_numeric_tools/cdp2rdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_2D

FUNCTION

INPUTS

OUTPUT

SOURCE

1423 pure function cdp2rdp_2D(cc) result(rr)
1424 
1425 
1426 !This section has been created automatically by the script Abilint (TD).
1427 !Do not modify the following lines by hand.
1428 #undef ABI_FUNC
1429 #define ABI_FUNC 'cdp2rdp_2D'
1430 !End of the abilint section
1431 
1432  implicit none
1433 
1434 !Arguments ------------------------------------
1435 !scalars
1436  complex(dpc),intent(in) :: cc(:,:)
1437  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2))
1438 ! *********************************************************************
1439 
1440  rr(1,:,:)=REAL (cc(:,:))
1441  rr(2,:,:)=AIMAG(cc(:,:))
1442 
1443 end function cdp2rdp_2D

m_numeric_tools/cdp2rdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_3D

FUNCTION

INPUTS

OUTPUT

SOURCE

1460 pure function cdp2rdp_3D(cc) result(rr)
1461 
1462 
1463 !This section has been created automatically by the script Abilint (TD).
1464 !Do not modify the following lines by hand.
1465 #undef ABI_FUNC
1466 #define ABI_FUNC 'cdp2rdp_3D'
1467 !End of the abilint section
1468 
1469  implicit none
1470 
1471 !Arguments ------------------------------------
1472 !scalars
1473  complex(dpc),intent(in) :: cc(:,:,:)
1474  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3))
1475 
1476 ! *********************************************************************
1477 
1478  rr(1,:,:,:)=REAL (cc(:,:,:))
1479  rr(2,:,:,:)=AIMAG(cc(:,:,:))
1480 
1481 end function cdp2rdp_3D

m_numeric_tools/cdp2rdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_4D

FUNCTION

INPUTS

OUTPUT

SOURCE

1498 pure function cdp2rdp_4D(cc) result(rr)
1499 
1500 
1501 !This section has been created automatically by the script Abilint (TD).
1502 !Do not modify the following lines by hand.
1503 #undef ABI_FUNC
1504 #define ABI_FUNC 'cdp2rdp_4D'
1505 !End of the abilint section
1506 
1507  implicit none
1508 
1509 !Arguments ------------------------------------
1510 !scalars
1511  complex(dpc),intent(in) :: cc(:,:,:,:)
1512  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4))
1513 ! *********************************************************************
1514 
1515  rr(1,:,:,:,:)=REAL (cc(:,:,:,:))
1516  rr(2,:,:,:,:)=AIMAG(cc(:,:,:,:))
1517 
1518 end function cdp2rdp_4D

m_numeric_tools/cdp2rdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_5D

FUNCTION

INPUTS

OUTPUT

SOURCE

1535 pure function cdp2rdp_5D(cc) result(rr)
1536 
1537 
1538 !This section has been created automatically by the script Abilint (TD).
1539 !Do not modify the following lines by hand.
1540 #undef ABI_FUNC
1541 #define ABI_FUNC 'cdp2rdp_5D'
1542 !End of the abilint section
1543 
1544  implicit none
1545 
1546 !Arguments ------------------------------------
1547 !scalars
1548  complex(dpc),intent(in) :: cc(:,:,:,:,:)
1549  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5))
1550 
1551 ! *********************************************************************
1552 
1553  rr(1,:,:,:,:,:)=REAL (cc(:,:,:,:,:))
1554  rr(2,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:))
1555 
1556 end function cdp2rdp_5D

m_numeric_tools/central_finite_diff [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 central_finite_diff

FUNCTION

 Coefficients of the central differences, for several orders of accuracy.
 See: https://en.wikipedia.org/wiki/Finite_difference_coefficient

INPUTS

  order=Derivative order.
  ipos=Index of the point must be in [1,npts]
  npts=Number of points used in finite difference, origin at npts/2 + 1

OUTPUT

  coeffient for central finite difference

PARENTS

CHILDREN

SOURCE

6474 real(dp) function central_finite_diff(order, ipos, npts) result(fact)
6475 
6476 
6477 !This section has been created automatically by the script Abilint (TD).
6478 !Do not modify the following lines by hand.
6479 #undef ABI_FUNC
6480 #define ABI_FUNC 'central_finite_diff'
6481 !End of the abilint section
6482 
6483  implicit none
6484 
6485 !Arguments ---------------------------------------------
6486 !scalars
6487  integer,intent(in) :: ipos,order,npts
6488 
6489 !Local variables ---------------------------------------
6490 !scalars
6491  real(dp),parameter :: empty=huge(one)
6492 ! 1st derivative.
6493  real(dp),parameter :: d1(9,4) = reshape([ &
6494   [-1/2._dp, 0._dp, 1/2._dp, empty, empty, empty, empty, empty, empty], &
6495   [ 1/12._dp, -2/3._dp, 0._dp, 2/3._dp, -1/12._dp, empty, empty, empty, empty], &
6496   [-1/60._dp, 3/20._dp, -3/4._dp, 0._dp, 3/4._dp, -3/20._dp, 1/60._dp, empty, empty], &
6497   [ 1/280._dp, -4/105._dp, 1/5._dp, -4/5._dp, 0._dp, 4/5._dp, -1/5._dp, 4/105._dp, -1/280._dp]], [9,4])
6498 ! 2nd derivative.
6499  real(dp),parameter :: d2(9,4) = reshape([ &
6500    [ 1._dp, -2._dp, 1._dp, empty, empty, empty, empty, empty, empty], &
6501    [-1/12._dp, 4/3._dp, -5/2._dp, 4/3._dp, -1/12._dp, empty, empty, empty, empty], &
6502    [ 1/90._dp, -3/20._dp, 3/2._dp, -49/18._dp, 3/2._dp, -3/20._dp, 1/90._dp, empty, empty], &
6503    [-1/560._dp, 8/315._dp, -1/5._dp, 8/5._dp, -205/72._dp, 8/5._dp, -1/5._dp, 8/315._dp, -1/560._dp]], [9,4])
6504 ! 3th derivative.
6505  real(dp),parameter :: d3(9,3) = reshape([ &
6506    [-1/2._dp, 1._dp, 0._dp, -1._dp, 1/2._dp, empty, empty, empty, empty], &
6507    [ 1/8._dp, -1._dp, 13/8._dp, 0._dp, -13/8._dp, 1._dp, -1/8._dp, empty, empty], &
6508    [ -7/240._dp, 3/10._dp, -169/120._dp, 61/30._dp, 0._dp, -61/30._dp, 169/120._dp, -3/10._dp, 7/240._dp]], &
6509    [9,3])
6510 ! 4th derivative.
6511  real(dp),parameter :: d4(9,3) = reshape([ &
6512    [ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp, empty, empty, empty, empty], &
6513    [ -1/6._dp, 2._dp, -13/2._dp, 28/3._dp, -13/2._dp, 2._dp, -1/6._dp, empty, empty], &
6514    [ 7/240._dp, -2/5._dp, 169/60._dp, -122/15._dp, 91/8._dp, -122/15._dp, 169/60._dp, -2/5._dp, 7/240._dp]], [9,3])
6515 ! 5th derivative.
6516  real(dp),parameter :: d5(7) = [ -1/2._dp, 2._dp, -5/2._dp, 0._dp, 5/2._dp, -2._dp, 1/2._dp]
6517 ! 6th derivative.
6518  real(dp),parameter :: d6(7) = [ 1._dp, -6._dp, 15._dp, -20._dp, 15._dp, -6._dp, 1._dp]
6519 ! *********************************************************************
6520 
6521  select case (order)
6522  case (1)
6523    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
6524    fact = d1(ipos, npts/2)
6525  case (2)
6526    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
6527    fact = d2(ipos, npts/2)
6528  case (3)
6529    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
6530    fact = d3(ipos, npts/2)
6531  case (4)
6532    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
6533    fact = d4(ipos, npts/2 - 1)
6534  case (5)
6535    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
6536    fact = d5(ipos)
6537  case (6)
6538    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
6539    fact = d6(ipos)
6540  case default
6541    MSG_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
6542  end select
6543 
6544  if (fact == empty) then
6545    MSG_ERROR(sjoin("Invalid ipos:",itoa(ipos),"for order", itoa(order), "npts", itoa(npts)))
6546  end if
6547  return
6548 
6549 10 MSG_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
6550 
6551 end function central_finite_diff

m_numeric_tools/cmplx_sphcart [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 cmplx_sphcart

FUNCTION

 Convert an array of complex values stored in spherical coordinates
 to Cartesian coordinates with real and imaginary part or vice versa.

INPUTS

  from=Option specifying the format used to store the complex values. See below.
  [units]=Option to specify if angles are given in  "Radians" (default) or "Degrees".

SIDE EFFECTS

  carr(:,:):
    input:  array with complex values in Cartesian form if from="C" or spherical form if from="S"
    output: array with values converted to the new representation.

PARENTS

      m_ptgroups

CHILDREN

SOURCE

5286 subroutine cmplx_sphcart(carr, from, units)
5287 
5288 
5289 !This section has been created automatically by the script Abilint (TD).
5290 !Do not modify the following lines by hand.
5291 #undef ABI_FUNC
5292 #define ABI_FUNC 'cmplx_sphcart'
5293 !End of the abilint section
5294 
5295  implicit none
5296 
5297 !Arguments ------------------------------------
5298 !scalars
5299  character(len=*),intent(in) :: from
5300  character(len=*),optional,intent(in) :: units
5301 !arrays
5302  complex(dpc),intent(inout) :: carr(:,:)
5303 
5304 !Local variables-------------------------------
5305 !scalars
5306  integer :: jj,ii
5307  real(dp) :: rho,theta,fact
5308  character(len=500) :: msg
5309 
5310 ! *************************************************************************
5311 
5312  select case (from(1:1))
5313 
5314  case ("S","s") ! Spherical --> Cartesian
5315 
5316    fact = one
5317    if (PRESENT(units)) then
5318      if (units(1:1) == "D" .or. units(1:1) == "d") fact = two_pi/360_dp
5319    end if
5320 
5321    do jj=1,SIZE(carr,DIM=2)
5322      do ii=1,SIZE(carr,DIM=1)
5323         rho  = DBLE(carr(ii,jj))
5324         theta= AIMAG(carr(ii,jj)) * fact
5325         carr(ii,jj) = DCMPLX(rho*DCOS(theta), rho*DSIN(theta))
5326      end do
5327    end do
5328 
5329  case ("C","c") ! Cartesian --> Spherical \theta = 2 arctan(y/(rho+x))
5330 
5331    fact = one
5332    if (PRESENT(units)) then
5333      if (units(1:1) == "D" .or. units(1:1) == "d") fact = 360_dp/two_pi
5334    end if
5335 
5336    do jj=1,SIZE(carr,DIM=2)
5337      do ii=1,SIZE(carr,DIM=1)
5338         rho  = SQRT(ABS(carr(ii,jj)))
5339         if (rho > tol16) then
5340           theta= two * ATAN( AIMAG(carr(ii,jj)) / (DBLE(carr(ii,jj)) + rho) )
5341         else
5342           theta= zero
5343         end if
5344         carr(ii,jj) = DCMPLX(rho, theta*fact)
5345      end do
5346    end do
5347 
5348  case default
5349    msg = " Wrong value for from: "//TRIM(from)
5350    MSG_BUG(msg)
5351  end select
5352 
5353 end subroutine cmplx_sphcart

m_numeric_tools/coeffs_gausslegint [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  coeffs_gausslegint

FUNCTION

 Compute the coefficients (supports and weights) for Gauss-Legendre integration.
 Inspired by a routine due to G. Rybicki.

INPUTS

 xmin=lower bound of integration
 xmax=upper bound of integration
 n=order of integration

OUTPUT

 x(n)=array of support points
 weights(n)=array of integration weights

PARENTS

      calc_rpa_functional,calc_sigc_me,integrho,integvol,m_numeric_tools
      screening,surf

CHILDREN

SOURCE

3548 subroutine coeffs_gausslegint(xmin,xmax,x,weights,n)
3549 
3550 
3551 !This section has been created automatically by the script Abilint (TD).
3552 !Do not modify the following lines by hand.
3553 #undef ABI_FUNC
3554 #define ABI_FUNC 'coeffs_gausslegint'
3555 !End of the abilint section
3556 
3557  implicit none
3558 
3559 !Arguments ------------------------------------
3560 !scalars
3561  integer,intent(in) :: n
3562  real(dp),intent(in) :: xmin,xmax
3563  real(dp),intent(out) :: x(n),weights(n)
3564 
3565 !Local variables ------------------------------
3566 !scalars
3567  integer :: i,j
3568  real(dp),parameter :: tol=1.d-13
3569  real(dp),parameter :: pi=4.d0*atan(1.d0)
3570  real(dp) :: z,z1,xmean,p1,p2,p3,pp,xl
3571 
3572 !************************************************************************
3573 
3574  xl=(xmax-xmin)*0.5d0
3575  xmean=(xmax+xmin)*0.5d0
3576 
3577  do i=1,(n+1)/2
3578   z=cos(pi*(i-0.25d0)/(n+0.5d0))
3579 
3580   do
3581     p1=1.d0
3582     p2=0.d0
3583 
3584     do j=1,n
3585      p3=p2
3586      p2=p1
3587      p1=((2.d0*j - 1.d0)*z*p2 - (j-1.d0)*p3)/j
3588     end do
3589 
3590     pp=n*(p2-z*p1)/(1.0d0-z**2)
3591     z1=z
3592     z=z1-p1/pp
3593 
3594     if(abs(z-z1) < tol) exit
3595   end do
3596 
3597   x(i)=xmean-xl*z
3598   x(n+1-i)=xmean+xl*z
3599   weights(i)=2.d0*xl/((1.d0-z**2)*pp**2)
3600   weights(n+1-i)=weights(i)
3601  end do
3602 
3603 end subroutine coeffs_gausslegint

m_numeric_tools/continued_fract [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  continued_fract

FUNCTION

  This routine calculates the continued fraction:

                        1
 f(z) =  _______________________________
           z - a1 -        b1^2
                   _____________________
                     z - a2 -    b2^2
                             ___________
                                z -a3 -    ........

INPUTS

  nlev=Number of "levels" in the continued fraction.
  term_type=Type of the terminator.
    0 --> No terminator.
   -1 --> Assume constant coefficients for a_i and b_i for i>nlev with a_inf = a(nlev) and b_inf = b(nleb)
    1 --> Same as above but a_inf and b_inf are obtained by averaging over the nlev values.
  aa(nlev)=Set of a_i coefficients.
  bb(nlev)=Set of b_i coefficients.
  nz=Number of points on the z-mesh.
  zpts(nz)=z-mesh.

OUTPUT

  spectrum(nz)=Contains f(z) on the input mesh.

PARENTS

      bsepostproc,m_haydock

CHILDREN

SOURCE

5168 subroutine continued_fract(nlev,term_type,aa,bb,nz,zpts,spectrum)
5169 
5170 
5171 !This section has been created automatically by the script Abilint (TD).
5172 !Do not modify the following lines by hand.
5173 #undef ABI_FUNC
5174 #define ABI_FUNC 'continued_fract'
5175 !End of the abilint section
5176 
5177  implicit none
5178 
5179 !Arguments ------------------------------------
5180 !scalars
5181  integer,intent(in) :: nlev,term_type,nz
5182 !arrays
5183  real(dp),intent(in) ::  bb(nlev)
5184  complex(dpc),intent(in) :: aa(nlev)
5185  complex(dpc),intent(in) :: zpts(nz)
5186  complex(dpc),intent(out) :: spectrum(nz)
5187 
5188 !Local variables ------------------------------
5189 !scalars
5190  integer :: it
5191  real(dp) ::  bb_inf,bg,bu,swap
5192  complex(dpc) :: aa_inf
5193  character(len=500) :: msg
5194 !arrays
5195  complex(dpc),allocatable :: div(:),den(:)
5196 
5197 !************************************************************************
5198 
5199  ABI_MALLOC(div,(nz))
5200  ABI_MALLOC(den,(nz))
5201 
5202  select case (term_type)
5203  case (0) ! No terminator.
5204    div=czero
5205  case (-1,1)
5206    if (term_type==-1) then
5207      bb_inf=bb(nlev)
5208      aa_inf=aa(nlev)
5209    else
5210      bb_inf=SUM(bb)/nlev
5211      aa_inf=SUM(aa)/nlev
5212    end if
5213    ! Be careful with the sign of the SQRT.
5214    div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
5215  case (2)
5216    MSG_ERROR("To be tested")
5217    div = zero
5218    if (nlev>4) then
5219      bg=zero; bu=zero
5220      do it=1,nlev,2
5221        if (it+2<nlev) bg = bg + bb(it+2)
5222        bu = bu + bb(it)
5223      end do
5224      bg = bg/(nlev/2+MOD(nlev,2))
5225      bu = bg/((nlev+1)/2)
5226      !if (iseven(nlev)) then
5227      if (.not.iseven(nlev)) then
5228        swap = bg
5229        bg = bu
5230        bu = bg
5231      end if
5232      !write(std_out,*)nlev,bg,bu
5233      !Here be careful with the sign of SQRT
5234      do it=1,nz
5235        div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
5236 &        ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
5237      end do
5238    end if
5239 
5240  case default
5241    write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
5242    MSG_ERROR(msg)
5243  end select
5244 
5245  do it=nlev,2,-1
5246    den(:) = zpts(:) - aa(it) - div(:)
5247    div(:) = (bb(it-1)**2 )/ den(:)
5248  end do
5249 
5250  den = zpts(:) - aa(1) - div(:)
5251  div = one/den(:)
5252 
5253  spectrum = div
5254  ABI_FREE(div)
5255  ABI_FREE(den)
5256 
5257 end subroutine continued_fract

m_numeric_tools/cross_product_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cross_product_int

FUNCTION

  Return the cross product of two vectors with integer components.


m_numeric_tools/cross_product_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cross_product_rdp

FUNCTION

  Return the cross product of two vectors with real double precision components.


m_numeric_tools/cspint [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cspint

FUNCTION

  Estimates the integral of a tabulated function.

INPUTS

OUTPUT

NOTES

    The routine is given the value of a function F(X) at a set of
    nodes XTAB, and estimates

      Integral ( A <= X <= B ) F(X) DX

    by computing the cubic natural spline S(X) that interpolates
    F(X) at the nodes, and then computing

      Integral ( A <= X <= B ) S(X) DX

    exactly.

    Other output from the program includes the definite integral
    from X(1) to X(I) of S(X), and the coefficients necessary for
    the user to evaluate the spline S(X) at any point.

  Modified:

    30 October 2000

  Reference:

    Philip Davis and Philip Rabinowitz,
    Methods of Numerical Integration,
    Blaisdell Publishing, 1967.

  Parameters:

    Input, real (dp) FTAB(NTAB), contains the tabulated values of
    the function, FTAB(I) = F(XTAB(I)).

    Input, real (dp) XTAB(NTAB), contains the points at which the
    function was evaluated.  The XTAB's must be distinct and
    in ascending order.

    Input, integer NTAB, the number of entries in FTAB and
    XTAB.  NTAB must be at least 3.

    Input, real (dp) A, lower limit of integration.

    Input, real (dp) B, upper limit of integration.

    Output, real (dp) Y(3,NTAB), will contain the coefficients
    of the interpolating natural spline over each subinterval.

    For XTAB(I) <= X <= XTAB(I+1),

      S(X) = FTAB(I) + Y(1,I)*(X-XTAB(I))
                   + Y(2,I)*(X-XTAB(I))**2
                   + Y(3,I)*(X-XTAB(I))**3

    Output, real (dp) E(NTAB), E(I) = the definite integral from
    XTAB(1) to XTAB(I) of S(X).

    Workspace, real (dp) WORK(NTAB).

    Output, real (dp) RESULT, the estimated value of the integral.

PARENTS

      m_xc_vdw,mrgscr

CHILDREN

SOURCE

3379 subroutine cspint ( ftab, xtab, ntab, a, b, y, e, work, result )
3380 
3381 
3382 !This section has been created automatically by the script Abilint (TD).
3383 !Do not modify the following lines by hand.
3384 #undef ABI_FUNC
3385 #define ABI_FUNC 'cspint'
3386 !End of the abilint section
3387 
3388   implicit none
3389 
3390 !Arguments ------------------------------------
3391 !scalars
3392   integer, intent(in) :: ntab
3393   real(dp), intent(in) :: a
3394   real(dp), intent(in) :: b
3395   real(dp), intent(inout) :: e(ntab)
3396   real(dp), intent(in) :: ftab(ntab)
3397   real(dp), intent(inout) :: work(ntab)
3398   real(dp), intent(in) :: xtab(ntab)
3399   real(dp), intent(inout) :: y(3,ntab)
3400   real(dp), intent(out) :: result
3401 
3402 !Local variables ------------------------------
3403 !scalars
3404   integer :: i
3405   integer :: j
3406   real(dp) :: r
3407   real(dp) :: s
3408   real(dp) :: term
3409   real(dp) :: u
3410 
3411 !************************************************************************
3412 
3413   if ( ntab < 3 ) then
3414     write(std_out,'(a)' ) ' '
3415     write(std_out,'(a)' ) 'CSPINT - Fatal error!'
3416     write(std_out,'(a,i6)' ) '  NTAB must be at least 3, but input NTAB = ',ntab
3417     MSG_ERROR("Aborting now")
3418   end if
3419 
3420   do i = 1, ntab-1
3421 
3422     if ( xtab(i+1) <= xtab(i) ) then
3423       write(std_out,'(a)' ) ' '
3424       write(std_out,'(a)' ) 'CSPINT - Fatal error!'
3425       write(std_out,'(a)' ) '  Nodes not in strict increasing order.'
3426       write(std_out,'(a,i6)' ) '  XTAB(I) <= XTAB(I-1) for I=',i
3427       write(std_out,'(a,g14.6)' ) '  XTAB(I) = ',xtab(i)
3428       write(std_out,'(a,g14.6)' ) '  XTAB(I-1) = ',xtab(i-1)
3429       MSG_ERROR("Aborting now")
3430     end if
3431 
3432   end do
3433 
3434   s = zero
3435   do i = 1, ntab-1
3436     r = ( ftab(i+1) - ftab(i) ) / ( xtab(i+1) - xtab(i) )
3437     y(2,i) = r - s
3438     s = r
3439   end do
3440 
3441   result = zero
3442   s = zero
3443   r = zero
3444   y(2,1) = zero
3445   y(2,ntab) = zero
3446 
3447   do i = 2, ntab-1
3448     y(2,i) = y(2,i) + r * y(2,i-1)
3449     work(i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s
3450     s = xtab(i+1) - xtab(i)
3451     r = s / work(i)
3452   end do
3453 
3454   do j = 2, ntab-1
3455     i = ntab+1-j
3456     y(2,i) = ( ( xtab(i+1) - xtab(i) ) * y(2,i+1) - y(2,i) ) / work(i)
3457   end do
3458 
3459   do i = 1, ntab-1
3460     s = xtab(i+1) - xtab(i)
3461     r = y(2,i+1) - y(2,i)
3462     y(3,i) = r / s
3463     y(2,i) = three * y(2,i)
3464     y(1,i) = ( ftab(i+1) - ftab(i) ) / s - ( y(2,i) + r ) * s
3465   end do
3466 
3467   e(1) = 0.0D+00
3468   do i = 1, ntab-1
3469     s = xtab(i+1)-xtab(i)
3470     term = ((( y(3,i) * quarter * s + y(2,i) * third ) * s &
3471       + y(1,i) * half ) * s + ftab(i) ) * s
3472     e(i+1) = e(i) + term
3473   end do
3474 !
3475 !  Determine where the endpoints A and B lie in the mesh of XTAB's.
3476 !
3477   r = a
3478   u = one
3479 
3480   do j = 1, 2
3481 !
3482 !  The endpoint is less than or equal to XTAB(1).
3483 !
3484     if ( r <= xtab(1) ) then
3485       result = result-u*((r-xtab(1))*y(1,1)*half +ftab(1))*(r-xtab(1))
3486 !
3487 !  The endpoint is greater than or equal to XTAB(NTAB).
3488 !
3489     else if ( xtab(ntab) <= r ) then
3490 
3491       result = result -u * ( e(ntab) + ( r - xtab(ntab) ) &
3492         * ( ftab(ntab) + half * ( ftab(ntab-1) &
3493         + ( xtab(ntab) - xtab(ntab-1) ) * y(1,ntab-1) ) &
3494         * ( r - xtab(ntab) )))
3495 !
3496 !  The endpoint is strictly between XTAB(1) and XTAB(NTAB).
3497 !
3498     else
3499 
3500       do i = 1, ntab-1
3501 
3502         if ( r <= xtab(i+1) ) then
3503           r = r-xtab(i)
3504           result = result-u*(e(i)+(((y(3,i)*quarter*r+y(2,i)*third)*r &
3505             +y(1,i)*half )*r+ftab(i))*r)
3506           go to 120
3507         end if
3508 
3509       end do
3510 
3511     end if
3512 
3513   120   continue
3514 
3515     u = -one
3516     r = b
3517 
3518   end do
3519 
3520 end subroutine cspint

m_numeric_tools/ctrap [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 ctrap

FUNCTION

 Do corrected trapezoidal integral on uniform grid of spacing hh.

INPUTS

  imax=highest index of grid=grid point number of upper limit
  ff(imax)=integrand values
  hh=spacing between x points

OUTPUT

  ans=resulting integral by corrected trapezoid

NOTES

PARENTS

      psden,psp11lo,psp11nl,psp5lo,psp5nl,psp8lo,psp8nl,vhtnzc

CHILDREN

SOURCE

3213 subroutine ctrap(imax,ff,hh,ans)
3214 
3215 
3216 !This section has been created automatically by the script Abilint (TD).
3217 !Do not modify the following lines by hand.
3218 #undef ABI_FUNC
3219 #define ABI_FUNC 'ctrap'
3220 !End of the abilint section
3221 
3222  implicit none
3223 
3224 !Arguments ------------------------------------
3225 !scalars
3226  integer,intent(in) :: imax
3227  real(dp),intent(in) :: hh
3228  real(dp),intent(out) :: ans
3229 !arrays
3230  real(dp),intent(in) :: ff(imax)
3231 
3232 !Local variables-------------------------------
3233 !scalars
3234  integer :: ir,ir2
3235  real(dp) :: endpt,sum
3236 
3237 ! *************************************************************************
3238 
3239  if (imax>=10)then
3240 
3241 !  endpt=end point correction terms (low and high ends)
3242    endpt  = (23.75d0*(ff(1)+ff(imax  )) &
3243 &   + 95.10d0*(ff(2)+ff(imax-1)) &
3244 &   + 55.20d0*(ff(3)+ff(imax-2)) &
3245 &   + 79.30d0*(ff(4)+ff(imax-3)) &
3246 &   + 70.65d0*(ff(5)+ff(imax-4)))/ 72.d0
3247    ir2 = imax - 5
3248    sum=0.00d0
3249    if (ir2 > 5) then
3250      do ir=6,ir2
3251        sum = sum + ff(ir)
3252      end do
3253    end if
3254    ans = (sum + endpt ) * hh
3255 
3256  else if (imax>=8)then
3257    endpt  = (17.0d0*(ff(1)+ff(imax  )) &
3258 &   + 59.0d0*(ff(2)+ff(imax-1)) &
3259 &   + 43.0d0*(ff(3)+ff(imax-2)) &
3260 &   + 49.0d0*(ff(4)+ff(imax-3)) )/ 48.d0
3261    sum=0.0d0
3262    if(imax==9)sum=ff(5)
3263    ans = (sum + endpt ) * hh
3264 
3265  else if (imax==7)then
3266    ans = (17.0d0*(ff(1)+ff(imax  )) &
3267 &   + 59.0d0*(ff(2)+ff(imax-1)) &
3268 &   + 43.0d0*(ff(3)+ff(imax-2)) &
3269 &   + 50.0d0* ff(4)                )/ 48.d0  *hh
3270 
3271  else if (imax==6)then
3272    ans = (17.0d0*(ff(1)+ff(imax  )) &
3273 &   + 59.0d0*(ff(2)+ff(imax-1)) &
3274 &   + 44.0d0*(ff(3)+ff(imax-2)) )/ 48.d0  *hh
3275 
3276  else if (imax==5)then
3277    ans = (     (ff(1)+ff(5)) &
3278 &   + four*(ff(2)+ff(4)) &
3279 &   + two * ff(3)         )/ three  *hh
3280 
3281  else if (imax==4)then
3282    ans = (three*(ff(1)+ff(4)) &
3283 &   + nine *(ff(2)+ff(3))  )/ eight  *hh
3284 
3285  else if (imax==3)then
3286    ans = (     (ff(1)+ff(3)) &
3287 &   + four* ff(2)         )/ three *hh
3288 
3289  else if (imax==2)then
3290    ans = (ff(1)+ff(2))/ two  *hh
3291 
3292  else if (imax==1)then
3293    ans = ff(1)*hh
3294 
3295  end if
3296 
3297 end subroutine ctrap

m_numeric_tools/denominator [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

   denominator

FUNCTION

  Return the denominator of the rational number dd, sign is not considered.

INPUTS

  dd=The rational number
  tolerance=Absolute tolerance

OUTPUT

  ierr=If /=0  the input number is not rational within the given tolerance.

SOURCE

5044 integer function denominator(dd,ierr,tolerance)
5045 
5046 
5047 !This section has been created automatically by the script Abilint (TD).
5048 !Do not modify the following lines by hand.
5049 #undef ABI_FUNC
5050 #define ABI_FUNC 'denominator'
5051 !End of the abilint section
5052 
5053  implicit none
5054 
5055 !Arguments ------------------------------------
5056 !scalars
5057  integer,intent(out) :: ierr
5058  real(dp),intent(in) :: dd
5059  real(dp),optional,intent(in) :: tolerance
5060 
5061 !Local variables ------------------------------
5062 !scalars
5063  integer,parameter :: largest_integer = HUGE(1)
5064  integer :: ii
5065  real(dp) :: my_tol
5066 
5067 !************************************************************************
5068 
5069  ii=1
5070  my_tol=0.0001 ; if (PRESENT(tolerance)) my_tol=ABS(tolerance)
5071  do
5072    if (ABS(dd*ii-NINT(dd*ii))<my_tol) then
5073      denominator=ii
5074      ierr=0
5075      RETURN
5076    end if
5077    ! Handle the case in which dd is not rational within my_tol.
5078    if (ii==largest_integer) then
5079      denominator=ii
5080      ierr=-1
5081      RETURN
5082    end if
5083    ii=ii+1
5084  end do
5085 
5086 end function denominator

m_numeric_tools/dpade [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  dpade

FUNCTION

  Calculate the derivative of the pade approximant in zz of the function f calculated at the n points z

SOURCE

4666 function dpade(n,z,f,zz)
4667 
4668 
4669 !This section has been created automatically by the script Abilint (TD).
4670 !Do not modify the following lines by hand.
4671 #undef ABI_FUNC
4672 #define ABI_FUNC 'dpade'
4673 !End of the abilint section
4674 
4675  implicit none
4676 
4677 !Arguments ------------------------------------
4678 !scalars
4679  integer,intent(in) :: n
4680  complex(dpc),intent(in) :: zz
4681  complex(dpc) :: dpade
4682 !arrays
4683  complex(dpc),intent(in) :: z(n),f(n)
4684 
4685 !Local variables-------------------------------
4686 !scalars
4687  integer :: i
4688 !arrays
4689  complex(dpc) :: a(n)
4690  complex(dpc) :: Az(0:n), Bz(0:n)
4691  complex(dpc) :: dAz(0:n), dBz(0:n)
4692 ! *************************************************************************
4693 
4694  call calculate_pade_a(a,n,z,f)
4695 
4696  Az(0)=czero
4697  Az(1)=a(1)
4698  Bz(0)=cone
4699  Bz(1)=cone
4700  dAz(0)=czero
4701  dAz(1)=czero
4702  dBz(0)=czero
4703  dBz(1)=czero
4704 
4705  do i=1,n-1
4706    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4707    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4708    dAz(i+1)=dAz(i)+a(i+1)*Az(i-1)+(zz-z(i))*a(i+1)*dAz(i-1)
4709    dBz(i+1)=dBz(i)+a(i+1)*Bz(i-1)+(zz-z(i))*a(i+1)*dBz(i-1)
4710  end do
4711  !write(std_out,*) 'Bz(n)', Bz(n)
4712  if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) 'Bz(n)',Bz(n)
4713  !pade_approx = Az(n) / Bz(n)
4714  dpade=dAz(n)/Bz(n) -Az(n)*dBz(n)/(Bz(n)*Bz(n))
4715  !write(std_out,*) 'pade_approx ', pade_approx
4716 
4717 end function dpade

m_numeric_tools/findmin [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 findmin

FUNCTION

 Compute the minimum of a function whose value and derivative are known at two points.
 Also deduce different quantities at this predicted point, and at the two other points
 It uses a quartic interpolation, with the supplementary
 condition that the second derivative vanishes at one and
 only one point. See Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].
 For this option, lambda_1 must be 1 (new point),
 and lambda_2 must be 0 (old point).
 Also, if the derivative at the new point is more negative
 than the derivative at the old point, the predicted
 point cannot correspond to a minimum, but will be lambda=2.5_dp,
 if the energy of the second point is lower than the energy
 of the first point.

INPUTS

 etotal_1=first value of the function
 etotal_2=second value of the function
 dedv_1=first value of the derivative
 dedv_2=second value of the derivative
 lambda_1=first value of the argument
 lambda_2=second value of the argument

OUTPUT

 dedv_predict=predicted value of the derivative (usually zero,
  except if choice=4, if it happens that a minimum cannot be located,
  and a trial step is taken)
 d2edv2_predict=predicted value of the second derivative (not if choice=4)
 d2edv2_1=first value of the second derivative (not if choice=4)
 d2edv2_2=second value of the second derivative (not if choice=4)
 etotal_predict=predicted value of the function
 lambda_predict=predicted value of the argument
 status= 0 if everything went normally ;
         1 if negative second derivative
         2 if some other problem

PARENTS

      m_bfgs

CHILDREN

      wrtout

SOURCE

6695 subroutine findmin(dedv_1,dedv_2,dedv_predict,&
6696 & d2edv2_1,d2edv2_2,d2edv2_predict,&
6697 & etotal_1,etotal_2,etotal_predict,&
6698 & lambda_1,lambda_2,lambda_predict,status)
6699 
6700 
6701 !This section has been created automatically by the script Abilint (TD).
6702 !Do not modify the following lines by hand.
6703 #undef ABI_FUNC
6704 #define ABI_FUNC 'findmin'
6705 !End of the abilint section
6706 
6707  implicit none
6708 
6709 !Arguments ------------------------------------
6710 !scalars
6711  integer,intent(out) :: status
6712  real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2
6713  real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict
6714  real(dp),intent(out) :: etotal_predict,lambda_predict
6715 
6716 !Local variables-------------------------------
6717 !scalars
6718  real(dp) :: aa,bb,bbp,cc,ccp,d_lambda,dd
6719  real(dp) :: discr,ee,eep,lambda_shift,sum1,sum2,sum3,uu
6720  real(dp) :: uu3,vv,vv3
6721  character(len=500) :: message
6722 
6723 ! *************************************************************************
6724 
6725 !DEBUG
6726 !write(std_out,*)' findmin : enter'
6727 !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2
6728 !ENDDEBUG
6729 
6730  status=0
6731  d_lambda=lambda_1-lambda_2
6732 
6733 !DEBUG
6734 !do choice=3,1,-1
6735 !ENDDEBUG
6736 
6737  if(abs(lambda_1-1.0_dp)>tol12 .or. abs(lambda_2)>tol12) then
6738    message = '  For choice=4, lambda_1 must be 1 and lambda_2 must be 0.'
6739    MSG_BUG(message)
6740  end if
6741 
6742 !Evaluate quartic interpolation
6743 !etotal = aa + bb * lambda + cc * lambda**2 + dd * lambda**3 + ee * lambda**4
6744 !Impose positive second derivative everywhere, with
6745 !one point where it vanishes :  3*dd**2=8*cc*ee
6746  aa=etotal_2
6747  bb=dedv_2
6748  sum1=etotal_1-aa-bb
6749  sum2=dedv_1-bb
6750  sum3=sum2-2.0_dp*sum1
6751 
6752 !Build the discriminant of the associated 2nd degree equation
6753  discr=sum2**2-3.0_dp*sum3**2
6754  if(discr<0.0_dp .or. sum2<0.0_dp)then
6755 
6756 ! jmb init
6757    d2edv2_2=0.0
6758    d2edv2_1=0.0
6759    d2edv2_predict=0.0
6760 
6761 !  Even if there is a problem, try to keep going ...
6762    message = 'The 2nd degree equation has no positive root (choice=4).'
6763    MSG_WARNING(message)
6764    status=2
6765    if(etotal_1<etotal_2)then
6766      write(message, '(a,a,a)' )&
6767 &     'Will continue, since the new total energy is lower',ch10,&
6768 &     'than the old. Take a larger step in the same direction.'
6769      MSG_COMMENT(message)
6770      lambda_predict=2.5_dp
6771    else
6772      write(message, '(a,a,a,a,a)' )&
6773 &     'There is a problem, since the new total energy is larger',ch10,&
6774 &     'than the old (choice=4).',ch10,&
6775 &     'I take a point between the old and new, close to the old .'
6776      MSG_COMMENT(message)
6777      lambda_predict=0.25_dp
6778    end if
6779 !  Mimick a zero-gradient lambda, in order to avoid spurious
6780 !  action of the inverse hessian (the next line would be a realistic estimation)
6781    dedv_predict=0.0_dp
6782 !  dedv_predict=dedv_2+lambda_predict*(dedv_1-dedv_2)
6783 !  Uses the energies, and the gradient at lambda_2
6784    etotal_predict=etotal_2+dedv_2*lambda_predict&
6785 &   +(etotal_1-etotal_2-dedv_2)*lambda_predict**2
6786 
6787  else
6788 
6789 !  Here, there is an acceptable solution to the 2nd degree equation
6790    discr=sqrt(discr)
6791 !  The root that gives the smallest ee corresponds to  -discr
6792 !  This is the one to be used: one aims at modelling the
6793 !  behaviour of the function as much as possible with the
6794 !  lowest orders of the polynomial, not the quartic term.
6795    ee=(sum2-discr)*0.5_dp
6796    dd=sum3-2.0_dp*ee
6797    cc=sum1-dd-ee
6798 
6799 !  DEBUG
6800 !  write(std_out,*)'aa,bb,cc,dd,ee',aa,bb,cc,dd,ee
6801 !  ENDDEBUG
6802 
6803 !  Now, must find the unique root of
6804 !$0 = bb + 2*cc * lambda + 3*dd * lambda^2 + 4*ee * lambda^3$
6805 !  This root is unique because it was imposed that the second derivative
6806 !  of the quartic polynomial is everywhere positive.
6807 !  First, remove the quadratic term, by a shift of lambda
6808 !  lambdap=lambda-lambda_shift
6809 !$0 = bbp + ccp * lambdap + eep * lambdap^3$
6810    eep=4.0_dp*ee
6811    lambda_shift=-dd/(4.0_dp*ee)
6812    ccp=2.0_dp*cc-12.0_dp*ee*lambda_shift**2
6813    bbp=bb+ccp*lambda_shift+eep*lambda_shift**3
6814 
6815 !  DEBUG
6816 !  write(std_out,*)'bbp,ccp,eep,lambda_shift',bbp,ccp,eep,lambda_shift
6817 !  ENDDEBUG
6818 
6819 !  The solution of a cubic polynomial equation is as follows :
6820    discr=(bbp/eep)**2+(4.0_dp/27.0_dp)*(ccp/eep)**3
6821 !  In the present case, discr will always be positive
6822    discr=sqrt(discr)
6823    uu3=0.5_dp*(-bbp/eep+discr) ; uu=sign((abs(uu3))**(1.0_dp/3.0_dp),uu3)
6824    vv3=0.5_dp*(-bbp/eep-discr) ; vv=sign((abs(vv3))**(1.0_dp/3.0_dp),vv3)
6825    lambda_predict=uu+vv
6826 
6827 !  Restore the shift
6828    lambda_predict=lambda_predict+lambda_shift
6829    etotal_predict=aa+bb*lambda_predict+cc*lambda_predict**2+&
6830 &   dd*lambda_predict**3+ee*lambda_predict**4
6831    dedv_predict=bb+2.0_dp*cc*lambda_predict+3.0_dp*dd*lambda_predict**2+&
6832 &   4.0_dp*ee*lambda_predict**3
6833    d2edv2_1=2*cc+6*dd*lambda_1+12*ee*lambda_1**2
6834    d2edv2_2=2*cc+6*dd*lambda_2+12*ee*lambda_2**2
6835    d2edv2_predict=2*cc+6*dd*lambda_predict+12*ee*lambda_predict**2
6836 
6837  end if
6838 
6839  write(message, '(a,i3)' )'   line minimization, algorithm ',4
6840  call wrtout(std_out,message,'COLL')
6841  write(message, '(a,a)' )'                        lambda      etotal ','           dedv        d2edv2    '
6842  call wrtout(std_out,message,'COLL')
6843  write(message, '(a,es12.4,es18.10,2es12.4)' )'   old point         :',lambda_2,etotal_2,dedv_2,d2edv2_2
6844  call wrtout(std_out,message,'COLL')
6845  write(message, '(a,es12.4,es18.10,2es12.4)' )'   new point         :',lambda_1,etotal_1,dedv_1,d2edv2_1
6846  call wrtout(std_out,message,'COLL')
6847  write(message, '(a,es12.4,es18.10,2es12.4)' )'   predicted point   :',lambda_predict,etotal_predict,dedv_predict,d2edv2_predict
6848  call wrtout(std_out,message,'COLL')
6849  write(message, '(a)' ) ' '
6850  call wrtout(std_out,message,'COLL')
6851 
6852 end subroutine findmin

m_numeric_tools/geop [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  geop

FUNCTION

  Returns an array of length nn containing a geometric progression whose
  starting value is start and whose factor is factor!

INPUTS

  start=initial point
  factor=the factor of the geometric progression
  nn=the number of points

OUTPUT

  geop(nn)=the progression

SOURCE

400 pure function geop(start,factor,nn) result(res)
401 
402 
403 !This section has been created automatically by the script Abilint (TD).
404 !Do not modify the following lines by hand.
405 #undef ABI_FUNC
406 #define ABI_FUNC 'geop'
407 !End of the abilint section
408 
409  implicit none
410 
411 !Arguments ------------------------------------
412 !scalars
413  real(dp),intent(in) :: start,factor
414  integer,intent(in) :: nn
415  real(dp) :: res(nn)
416 
417 !Local variables-------------------------------
418  integer :: ii
419 ! *********************************************************************
420 
421  if (nn>0) res(1)=start
422  do ii=2,nn
423    res(ii)=res(ii-1)*factor
424  end do
425 
426 end function geop

m_numeric_tools/get_diag_cdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_cdp

FUNCTION

INPUTS

OUTPUT

SOURCE

896 function get_diag_cdp(cmat) result(cdiag)
897 
898 !Arguments ------------------------------------
899 !scalars
900 
901 !This section has been created automatically by the script Abilint (TD).
902 !Do not modify the following lines by hand.
903 #undef ABI_FUNC
904 #define ABI_FUNC 'get_diag_cdp'
905 !End of the abilint section
906 
907  complex(dpc),intent(in) :: cmat(:,:)
908  complex(dpc) :: cdiag(SIZE(cmat,1))
909 
910 !Local variables-------------------------------
911  integer :: ii
912 ! *************************************************************************
913 
914  ii=assert_eq(SIZE(cmat,1),SIZE(cmat,2),'Matrix not square',&
915 & __FILE__,__LINE__)
916 
917  do ii=1,SIZE(cmat,1)
918   cdiag(ii)=cmat(ii,ii)
919  end do
920 
921 end function get_diag_cdp

m_numeric_tools/get_diag_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_int

FUNCTION

  Return the diagonal of a square matrix as a vector

INPUTS

  matrix(:,:)

OUTPUT

  diag(:)=the diagonal

SOURCE

806 function get_diag_int(mat) result(diag)
807 
808 
809 !This section has been created automatically by the script Abilint (TD).
810 !Do not modify the following lines by hand.
811 #undef ABI_FUNC
812 #define ABI_FUNC 'get_diag_int'
813 !End of the abilint section
814 
815  implicit none
816 
817 !Arguments ------------------------------------
818 !scalars
819  integer,intent(in) :: mat(:,:)
820  integer :: diag(SIZE(mat,1))
821 
822 !Local variables-------------------------------
823  integer :: ii
824 ! *************************************************************************
825 
826  ii=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
827 & __FILE__,__LINE__)
828 
829  do ii=1,SIZE(mat,1)
830   diag(ii)=mat(ii,ii)
831  end do
832 
833 end function get_diag_int

m_numeric_tools/get_diag_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_rdp

FUNCTION

  Return the diagonal of a square matrix as a vector

INPUTS

  matrix(:,:)

OUTPUT

  diag(:)=the diagonal

SOURCE

852 function get_diag_rdp(mat) result(diag)
853 
854 
855 !This section has been created automatically by the script Abilint (TD).
856 !Do not modify the following lines by hand.
857 #undef ABI_FUNC
858 #define ABI_FUNC 'get_diag_rdp'
859 !End of the abilint section
860 
861  implicit none
862 
863 !Arguments ------------------------------------
864 !scalars
865  real(dp),intent(in) :: mat(:,:)
866  real(dp) :: diag(SIZE(mat,1))
867 
868 !Local variables-------------------------------
869  integer :: ii
870 ! *************************************************************************
871 
872  ii=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
873 & __FILE__,__LINE__)
874 
875  do ii=1,SIZE(mat,1)
876   diag(ii)=mat(ii,ii)
877  end do
878 
879 end function get_diag_rdp

m_numeric_tools/get_trace_cdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_trace_cdp

FUNCTION

  Calculate the trace of a square matrix (complex(dpc) version)

INPUTS

OUTPUT

SOURCE

763 pure function get_trace_cdp(matrix) result(trace)
764 
765 
766 !This section has been created automatically by the script Abilint (TD).
767 !Do not modify the following lines by hand.
768 #undef ABI_FUNC
769 #define ABI_FUNC 'get_trace_cdp'
770 !End of the abilint section
771 
772  implicit none
773 
774 !Arguments ------------------------------------
775  complex(dpc) :: trace
776  complex(dpc),intent(in) :: matrix(:,:)
777 
778 !Local variables-------------------------------
779 !scalars
780  integer :: ii
781 ! *********************************************************************
782 
783  trace=czero
784  do ii=1,size(matrix,dim=1)
785     trace=trace+matrix(ii,ii)
786  end do
787 
788 end function get_trace_cdp

m_numeric_tools/get_trace_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_trace_int

FUNCTION

  Calculate the trace of a square matrix

INPUTS

  matrix(:,:)

OUTPUT

  trace=the trace

SOURCE

675 pure function get_trace_int(matrix) result(trace)
676 
677 
678 !This section has been created automatically by the script Abilint (TD).
679 !Do not modify the following lines by hand.
680 #undef ABI_FUNC
681 #define ABI_FUNC 'get_trace_int'
682 !End of the abilint section
683 
684  implicit none
685 
686 !Arguments ------------------------------------
687  integer :: trace
688  integer,intent(in) :: matrix(:,:)
689 
690 !Local variables-------------------------------
691 !scalars
692  integer :: ii
693 ! *********************************************************************
694 
695  trace=0
696  do ii=1,size(matrix,dim=1)
697    trace=trace+matrix(ii,ii)
698  end do
699 
700 end function get_trace_int

m_numeric_tools/get_trace_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_trace_int

FUNCTION

  Calculate the trace of a square matrix (real(dp) version)

INPUTS

  matrix(:,:)

OUTPUT

  trace=the trace

SOURCE

720 pure function get_trace_rdp(matrix) result(trace)
721 
722 
723 !This section has been created automatically by the script Abilint (TD).
724 !Do not modify the following lines by hand.
725 #undef ABI_FUNC
726 #define ABI_FUNC 'get_trace_rdp'
727 !End of the abilint section
728 
729  implicit none
730 
731 !Arguments ------------------------------------
732  real(dp) :: trace
733  real(dp),intent(in) :: matrix(:,:)
734 
735 !Local variables-------------------------------
736 !scalars
737  integer :: ii
738 ! *********************************************************************
739 
740  trace=zero
741  do ii=1,size(matrix,dim=1)
742     trace=trace+matrix(ii,ii)
743  end do
744 
745 end function get_trace_rdp

m_numeric_tools/hermit [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 hermit

FUNCTION

 Take a matrix in hermitian storage mode (lower triangle stored)
 and redefine diagonal elements to impose Hermiticity
 (diagonal terms have to be real).
 If abs(Im(H(i,i)))>4096*machine precision, print error warning.
 (Typical 64 bit machine precision is 2^-52 or 2.22e-16)

INPUTS

  chmin(n*n+n)=complex hermitian matrix with numerical noise possibly
   rendering Im(diagonal elements) approximately 1e-15 or so
  ndim=size of complex hermitian matrix

OUTPUT

  chmout(n*n+n)=redefined matrix with strictly real diagonal elements.
   May be same storage location as chmin.
  ierr=0 if no problem, 1 if the imaginary part of some element
   too large (at present, stop in this case).

TODO

  Name is misleading, perhaps hermit_force_diago?
  Interface allows aliasing

PARENTS

      extrapwf,subdiago

CHILDREN

      wrtout

SOURCE

3974 subroutine hermit(chmin,chmout,ierr,ndim)
3975 
3976 
3977 !This section has been created automatically by the script Abilint (TD).
3978 !Do not modify the following lines by hand.
3979 #undef ABI_FUNC
3980 #define ABI_FUNC 'hermit'
3981 !End of the abilint section
3982 
3983  implicit none
3984 
3985 !Arguments ------------------------------------
3986 !scalars
3987  integer,intent(in) :: ndim
3988  integer,intent(out) :: ierr
3989 !arrays
3990  real(dp),intent(inout) :: chmin(ndim*ndim+ndim)
3991  real(dp),intent(inout) :: chmout(ndim*ndim+ndim)
3992 
3993 !Local variables-------------------------------
3994 !scalars
3995  integer,save :: mmesgs=20,nmesgs=0
3996  integer :: idim,merrors,nerrors
3997  real(dp),parameter :: eps=epsilon(0.0d0)
3998  real(dp) :: ch_im,ch_re,moduls,tol
3999  character(len=500) :: message
4000 
4001 ! *************************************************************************
4002 
4003  tol=4096.0d0*eps
4004 
4005  ierr=0
4006  merrors=0
4007 
4008 !Copy matrix into possibly new location
4009  chmout(:)=chmin(:)
4010 
4011 !Loop over diagonal elements of matrix (off-diag not altered)
4012  do idim=1,ndim
4013 
4014    ch_im=chmout(idim*idim+idim  )
4015    ch_re=chmout(idim*idim+idim-1)
4016 
4017 !  check for large absolute Im part and print warning when
4018 !  larger than (some factor)*(machine precision)
4019    nerrors=0
4020    if( abs(ch_im) > tol .and. abs(ch_im) > tol8*abs(ch_re)) nerrors=2
4021    if( abs(ch_im) > tol .or. abs(ch_im) > tol8*abs(ch_re)) nerrors=1
4022 
4023    if( (abs(ch_im) > tol .and. nmesgs<mmesgs) .or. nerrors==2)then
4024      write(message, '(3a,i0,a,es20.12,a,es20.12,a)' )&
4025 &     'Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
4026 &     'for component',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
4027      call wrtout(std_out,message,'PERS')
4028      nmesgs=nmesgs+1
4029    end if
4030 
4031    if( ( abs(ch_im) > tol8*abs(ch_re) .and. nmesgs<mmesgs) .or. nerrors==2)then
4032      write(message, '(3a,i0,a,es20.12,a,es20.12,a)' )&
4033 &     'Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
4034 &     'for component',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
4035      call wrtout(std_out,message,'PERS')
4036      nmesgs=nmesgs+1
4037    end if
4038 
4039 !  compute modulus $= (\Re^2+\Im^2)^{1/2}$
4040    moduls=sqrt(ch_re**2+ch_im**2)
4041 
4042 !  set Re part to modulus with sign of original Re part
4043    chmout(idim*idim+idim-1)=sign(moduls,ch_re)
4044 
4045 !  set Im part to 0
4046    chmout(idim*idim+idim)=zero
4047 
4048    merrors=max(merrors,nerrors)
4049 
4050  end do
4051 
4052  if(merrors==2)then
4053    ierr=1
4054    write(message, '(3a)' )&
4055 &   'Imaginary part(s) of diagonal Hermitian matrix element(s) is too large.',ch10,&
4056 &   'See previous messages.'
4057    MSG_BUG(message)
4058  end if
4059 
4060 end subroutine hermit

m_numeric_tools/hermitianize_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  hermitianize_dpc

FUNCTION

  Force a square matrix to be hermitian

INPUTS

  uplo=String describing which part of the matrix has been calculated.
    Only the first character is tested (no case sensitive). Possible values are:
    "All"= Full matrix is supplied in input
    "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
    "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.

OUTPUT

  (see side effects)

SIDE EFFECTS

  mat(:,:)=complex input matrix, hermitianized in output

PARENTS

CHILDREN

SOURCE

3810 subroutine hermitianize_dpc(mat,uplo)
3811 
3812 
3813 !This section has been created automatically by the script Abilint (TD).
3814 !Do not modify the following lines by hand.
3815 #undef ABI_FUNC
3816 #define ABI_FUNC 'hermitianize_dpc'
3817 !End of the abilint section
3818 
3819  implicit none
3820 
3821 !Arguments ------------------------------------
3822 !scalars
3823  character(len=*),intent(in) :: uplo
3824 !arrays
3825  complex(dpc),intent(inout) :: mat(:,:)
3826 
3827 !Local variables-------------------------------
3828 !scalars
3829  integer :: nn,ii,jj
3830 !arrays
3831  complex(dpc),allocatable :: tmp(:)
3832 ! *************************************************************************
3833 
3834  nn=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
3835 & __FILE__,__LINE__)
3836 
3837  select case (uplo(1:1))
3838 
3839  case ("A","a") ! Full matrix has been calculated.
3840    ABI_MALLOC(tmp,(nn))
3841    do ii=1,nn
3842      do jj=ii,nn
3843        tmp(jj)=half*(mat(ii,jj)+DCONJG(mat(jj,ii)))
3844      end do
3845      mat(ii,ii:nn)=tmp(ii:nn)
3846      mat(ii:nn,ii)=DCONJG(tmp(ii:nn))
3847    end do
3848    ABI_FREE(tmp)
3849 
3850  case ("U","u") ! Only the upper triangle is used.
3851    do jj=1,nn
3852      do ii=1,jj
3853        if (ii/=jj) then
3854          mat(jj,ii) = DCONJG(mat(ii,jj))
3855        else
3856          mat(ii,ii) = DCMPLX(DBLE(mat(ii,ii)),zero)
3857        end if
3858      end do
3859    end do
3860 
3861  case ("L","l") ! Only the lower triangle is used.
3862   do jj=1,nn
3863     do ii=1,jj
3864       if (ii/=jj) then
3865         mat(ii,jj) = DCONJG(mat(jj,ii))
3866       else
3867         mat(ii,ii) = DCMPLX(REAL(mat(ii,ii)),zero)
3868       end if
3869     end do
3870   end do
3871 
3872  case default
3873    MSG_ERROR("Wrong uplo"//TRIM(uplo))
3874  end select
3875 
3876 end subroutine hermitianize_dpc

m_numeric_tools/hermitianize_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  hermitianize_spc

FUNCTION

  Force a square matrix to be hermitian

INPUTS

  uplo=String describing which part of the matrix has been calculated.
    Only the first character is tested (no case sensitive). Possible values are:
    "All"= Full matrix is supplied in input
    "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
    "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.

OUTPUT

  (see side effects)

SIDE EFFECTS

  mat(:,:)=complex input matrix, hermitianized at output

PARENTS

CHILDREN

SOURCE

3713 subroutine hermitianize_spc(mat,uplo)
3714 
3715 
3716 !This section has been created automatically by the script Abilint (TD).
3717 !Do not modify the following lines by hand.
3718 #undef ABI_FUNC
3719 #define ABI_FUNC 'hermitianize_spc'
3720 !End of the abilint section
3721 
3722  implicit none
3723 
3724 !Arguments ------------------------------------
3725 !scalars
3726  character(len=*),intent(in) :: uplo
3727 !arrays
3728  complex(spc),intent(inout) :: mat(:,:)
3729 
3730 !Local variables-------------------------------
3731 !scalars
3732  integer :: nn,ii,jj
3733 !arrays
3734  complex(spc),allocatable :: tmp(:)
3735 ! *************************************************************************
3736 
3737  nn=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
3738 & __FILE__,__LINE__)
3739 
3740  select case (uplo(1:1))
3741 
3742  case ("A","a") ! Full matrix has been calculated.
3743    ABI_MALLOC(tmp,(nn))
3744    do ii=1,nn
3745      do jj=ii,nn
3746        tmp(jj)=half*(mat(ii,jj)+CONJG(mat(jj,ii)))
3747      end do
3748      mat(ii,ii:nn)=tmp(ii:nn)
3749      mat(ii:nn,ii)=CONJG(tmp(ii:nn))
3750    end do
3751    ABI_FREE(tmp)
3752 
3753  case ("U","u") ! Only the upper triangle is used.
3754    do jj=1,nn
3755      do ii=1,jj
3756        if (ii/=jj) then
3757          mat(jj,ii) = CONJG(mat(ii,jj))
3758        else
3759          mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3760        end if
3761      end do
3762    end do
3763 
3764  case ("L","l") ! Only the lower triangle is used.
3765   do jj=1,nn
3766     do ii=1,jj
3767       if (ii/=jj) then
3768         mat(ii,jj) = CONJG(mat(jj,ii))
3769       else
3770         mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3771       end if
3772     end do
3773   end do
3774 
3775  case default
3776    MSG_ERROR("Wrong uplo"//TRIM(uplo))
3777  end select
3778 
3779 end subroutine hermitianize_spc

m_numeric_tools/imax_loc_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  imax_loc_int

FUNCTION

  Index of maxloc on an array returned as scalar instead of array-valued

SOURCE

1924 pure function imax_loc_int(iarr,mask)
1925 
1926 
1927 !This section has been created automatically by the script Abilint (TD).
1928 !Do not modify the following lines by hand.
1929 #undef ABI_FUNC
1930 #define ABI_FUNC 'imax_loc_int'
1931 !End of the abilint section
1932 
1933  implicit none
1934 
1935 !Arguments ------------------------------------
1936 !scalars
1937  integer :: imax_loc_int
1938 !arrays
1939  integer,intent(in) :: iarr(:)
1940  logical,optional,intent(in) :: mask(:)
1941 
1942 !Local variables-------------------------------
1943  integer :: imax(1)
1944 ! *************************************************************************
1945 
1946  if (PRESENT(mask)) then
1947   imax=MAXLOC(iarr,MASK=mask)
1948  else
1949   imax=MAXLOC(iarr)
1950  end if
1951  imax_loc_int=imax(1)
1952 
1953 end function imax_loc_int

m_numeric_tools/imax_loc_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  imax_loc_rdp

FUNCTION

INPUTS

OUTPUT

SOURCE

1969 pure function imax_loc_rdp(arr,mask)
1970 
1971 
1972 !This section has been created automatically by the script Abilint (TD).
1973 !Do not modify the following lines by hand.
1974 #undef ABI_FUNC
1975 #define ABI_FUNC 'imax_loc_rdp'
1976 !End of the abilint section
1977 
1978  implicit none
1979 
1980 !Arguments ------------------------------------
1981 !scalars
1982  integer :: imax_loc_rdp
1983 !arrays
1984  real(dp),intent(in) :: arr(:)
1985  logical,optional,intent(in) :: mask(:)
1986 
1987 !Local variables-------------------------------
1988  integer :: imax(1)
1989 ! *************************************************************************
1990 
1991  if (PRESENT(mask)) then
1992   imax=MAXLOC(arr,MASK=mask)
1993  else
1994   imax=MAXLOC(arr)
1995  end if
1996  imax_loc_rdp=imax(1)
1997 
1998 end function imax_loc_rdp

m_numeric_tools/imin_loc_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  imin_loc_int

FUNCTION

  Index of minloc on an array returned as scalar instead of array-valued

SOURCE

2012 pure function imin_loc_int(arr,mask)
2013 
2014 
2015 !This section has been created automatically by the script Abilint (TD).
2016 !Do not modify the following lines by hand.
2017 #undef ABI_FUNC
2018 #define ABI_FUNC 'imin_loc_int'
2019 !End of the abilint section
2020 
2021  implicit none
2022 
2023 !Arguments ------------------------------------
2024 !scalars
2025  integer :: imin_loc_int
2026 !arrays
2027  integer,intent(in) :: arr(:)
2028  logical,optional,intent(in) :: mask(:)
2029 
2030 !Local variables-------------------------------
2031  integer :: imin(1)
2032 ! *************************************************************************
2033 
2034  if (PRESENT(mask)) then
2035   imin=MINLOC(arr,MASK=mask)
2036  else
2037   imin=MINLOC(arr)
2038  end if
2039  imin_loc_int=imin(1)
2040 
2041 end function imin_loc_int

m_numeric_tools/imin_loc_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  imin_loc_rdp

FUNCTION

INPUTS

OUTPUT

SOURCE

2058 pure function imin_loc_rdp(arr,mask)
2059 
2060 
2061 !This section has been created automatically by the script Abilint (TD).
2062 !Do not modify the following lines by hand.
2063 #undef ABI_FUNC
2064 #define ABI_FUNC 'imin_loc_rdp'
2065 !End of the abilint section
2066 
2067  implicit none
2068 
2069 !Arguments ------------------------------------
2070 !scalars
2071  integer :: imin_loc_rdp
2072 !arrays
2073  real(dp),intent(in) :: arr(:)
2074  logical,optional,intent(in) :: mask(:)
2075 
2076 !Local variables-------------------------------
2077  integer :: imin(1)
2078 ! *************************************************************************
2079 
2080  if (PRESENT(mask)) then
2081   imin=MINLOC(arr,MASK=mask)
2082  else
2083   imin=MINLOC(arr)
2084  end if
2085 
2086  imin_loc_rdp=imin(1)
2087 
2088 end function imin_loc_rdp

m_numeric_tools/interpol3d [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 interpol3d

FUNCTION

 Computes the value at any point r by linear interpolation
 inside the eight vertices of the surrounding cube
 r is presumed to be normalized, in a unit cube for the full grid

INPUTS

 r(3)=point coordinate
 nr1=grid size along x
 nr2=grid size along y
 nr3=grid size along z
 grid(nr1,nr2,nr3)=grid matrix

OUTPUT

 res=Interpolated value

PARENTS

      integrate_gamma_alt,lin_interpq_gam,lineint,m_nesting,m_qparticles
      planeint,pointint,volumeint

CHILDREN

SOURCE

5725 pure function interpol3d(r,nr1,nr2,nr3,grid) result(res)
5726 
5727 
5728 !This section has been created automatically by the script Abilint (TD).
5729 !Do not modify the following lines by hand.
5730 #undef ABI_FUNC
5731 #define ABI_FUNC 'interpol3d'
5732 !End of the abilint section
5733 
5734  implicit none
5735 
5736 !Arguments-------------------------------------------------------------
5737 !scalars
5738  integer,intent(in) :: nr1,nr2,nr3
5739  real(dp) :: res
5740 !arrays
5741  real(dp),intent(in) :: grid(nr1,nr2,nr3),r(3)
5742 
5743 !Local variables--------------------------------------------------------
5744 !scalars
5745  integer :: ir1,ir2,ir3,pr1,pr2,pr3
5746  real(dp) :: x1,x2,x3
5747 
5748 ! *************************************************************************
5749 
5750  call interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3, pr1,pr2,pr3)
5751 
5752 !weight
5753  x1=one+r(1)*nr1-real(ir1)
5754  x2=one+r(2)*nr2-real(ir2)
5755  x3=one+r(3)*nr3-real(ir3)
5756 
5757 !calculation of the density value
5758  res=zero
5759  res=res + grid(ir1,ir2,ir3)*(one-x1)*(one-x2)*(one-x3)
5760  res=res + grid(pr1,ir2,ir3)*x1*(one-x2)*(one-x3)
5761  res=res + grid(ir1,pr2,ir3)*(one-x1)*x2*(one-x3)
5762  res=res + grid(ir1,ir2,pr3)*(one-x1)*(one-x2)*x3
5763  res=res + grid(pr1,pr2,ir3)*x1*x2*(one-x3)
5764  res=res + grid(ir1,pr2,pr3)*(one-x1)*x2*x3
5765  res=res + grid(pr1,ir2,pr3)*x1*(one-x2)*x3
5766  res=res + grid(pr1,pr2,pr3)*x1*x2*x3
5767 
5768 end function interpol3d

m_numeric_tools/interpol3d_indices [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 interpol3d_indices

FUNCTION

 Computes the indices in a cube which are neighbors to the point to be
  interpolated in interpol3d

INPUTS

 r(3)=point coordinate
 nr1=grid size along x
 nr2=grid size along y
 nr3=grid size along z

OUTPUT

 ir1,ir2,ir3 = bottom left neighbor
 pr1,pr2,pr3 = top right neighbor

PARENTS

      interpol3d,k_neighbors

CHILDREN

SOURCE

5798 pure subroutine interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3,pr1,pr2,pr3)
5799 
5800 
5801 !This section has been created automatically by the script Abilint (TD).
5802 !Do not modify the following lines by hand.
5803 #undef ABI_FUNC
5804 #define ABI_FUNC 'interpol3d_indices'
5805 !End of the abilint section
5806 
5807  implicit none
5808 
5809 !Arguments-------------------------------------------------------------
5810 !scalars
5811  integer,intent(in) :: nr1,nr2,nr3
5812  integer,intent(out) :: ir1,ir2,ir3
5813  integer,intent(out) :: pr1,pr2,pr3
5814 !arrays
5815  real(dp),intent(in) :: r(3)
5816 
5817 !Local variables-------------------------------
5818  real(dp) :: d1,d2,d3
5819 
5820 ! *************************************************************************
5821 
5822 !grid density
5823  d1=one/nr1
5824  d2=one/nr2
5825  d3=one/nr3
5826 
5827 !lower left
5828  ir1=int(r(1)/d1)+1
5829  ir2=int(r(2)/d2)+1
5830  ir3=int(r(3)/d3)+1
5831 
5832 !upper right
5833  pr1=mod(ir1+1,nr1)
5834  pr2=mod(ir2+1,nr2)
5835  pr3=mod(ir3+1,nr3)
5836 
5837  if(ir1==0) ir1=nr1
5838  if(ir2==0) ir2=nr2
5839  if(ir3==0) ir3=nr3
5840 
5841  if(ir1>nr1) ir1=ir1-nr1
5842  if(ir2>nr2) ir2=ir2-nr2
5843  if(ir3>nr3) ir3=ir3-nr3
5844 
5845  if(pr1==0) pr1=nr1
5846  if(pr2==0) pr2=nr2
5847  if(pr3==0) pr3=nr3
5848 
5849 end subroutine interpol3d_indices

m_numeric_tools/interpolate_denpot [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 interpolate_denpot

FUNCTION

  Linear interpolation of density/potential given on the real space FFT mesh.
  Assumes array on full mesh i.e. no MPI-FFT.

INPUTS

  cplex=1 for real, 2 for complex data.
  in_ngfft(3)=Mesh divisions of input array
  nspden=Number of density components.
  in_rhor(cplex * in_nfftot * nspden)=Input array
  out_ngfft(3)=Mesh divisions of output array

OUTPUT

  outrhor(cplex * out_nfftot * nspden)=Output array with interpolated data.

PARENTS

      m_ioarr

CHILDREN

SOURCE

5879 subroutine interpolate_denpot(cplex, in_ngfft, nspden, in_rhor, out_ngfft, out_rhor)
5880 
5881 
5882 !This section has been created automatically by the script Abilint (TD).
5883 !Do not modify the following lines by hand.
5884 #undef ABI_FUNC
5885 #define ABI_FUNC 'interpolate_denpot'
5886 !End of the abilint section
5887 
5888  implicit none
5889 
5890 !Arguments-------------------------------------------------------------
5891 !scalars
5892  integer,intent(in) :: cplex,nspden
5893 !arrays
5894  integer,intent(in) :: in_ngfft(3), out_ngfft(3)
5895  real(dp),intent(in) :: in_rhor(cplex, product(in_ngfft), nspden)
5896  real(dp),intent(out) :: out_rhor(cplex, product(out_ngfft), nspden)
5897 
5898 !Local variables--------------------------------------------------------
5899 !scalars
5900  integer :: ispden, ir1, ir2, ir3, ifft
5901  real(dp) :: rr(3)
5902  real(dp),allocatable :: re(:,:),im(:,:)
5903 
5904 ! *************************************************************************
5905 
5906  if (cplex == 2) then
5907    ! copy slices for efficiency reasons (the best would be to have stride option in interpol3d)
5908    ABI_MALLOC(re, (product(in_ngfft), nspden))
5909    ABI_MALLOC(im, (product(in_ngfft), nspden))
5910    re = in_rhor(1, :, :)
5911    im = in_rhor(2, :, :)
5912  end if
5913 
5914  ! Linear interpolation.
5915  do ispden=1,nspden
5916    do ir3=0,out_ngfft(3)-1
5917      rr(3) = DBLE(ir3)/out_ngfft(3)
5918      do ir2=0,out_ngfft(2)-1
5919        rr(2) = DBLE(ir2)/out_ngfft(2)
5920        do ir1=0,out_ngfft(1)-1
5921          rr(1) = DBLE(ir1)/out_ngfft(1)
5922          ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5923          if (cplex == 1) then
5924            out_rhor(1, ifft, ispden) = interpol3d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_rhor(1, :, ispden))
5925          else
5926            out_rhor(1, ifft, ispden) = interpol3d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), re(:, ispden))
5927            out_rhor(2, ifft, ispden) = interpol3d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), im(:, ispden))
5928          end if
5929        end do
5930      end do
5931    end do
5932  end do
5933 
5934  if (cplex == 2) then
5935    ABI_FREE(re)
5936    ABI_FREE(im)
5937  end if
5938 
5939 end subroutine interpolate_denpot

m_numeric_tools/invcb [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 invcb

FUNCTION

 Compute a set of inverse cubic roots as fast as possible :
 rspts(:)=rhoarr(:)$^\frac{-1}{3}$

INPUTS

  npts=number of real space points on which density is provided
  rhoarr(npts)=input data

OUTPUT

  rspts(npts)=inverse cubic root of rhoarr

PARENTS

      drivexc,gammapositron,xchcth,xcpbe,xcpositron,xctfw

CHILDREN

SOURCE

7114 subroutine invcb(rhoarr,rspts,npts)
7115 
7116 
7117 !This section has been created automatically by the script Abilint (TD).
7118 !Do not modify the following lines by hand.
7119 #undef ABI_FUNC
7120 #define ABI_FUNC 'invcb'
7121 !End of the abilint section
7122 
7123  implicit none
7124 
7125 !Arguments ------------------------------------
7126 !scalars
7127  integer,intent(in) :: npts
7128 !arrays
7129  real(dp),intent(in) :: rhoarr(npts)
7130  real(dp),intent(out) :: rspts(npts)
7131 
7132 !Local variables-------------------------------
7133 !scalars
7134  integer :: ii,ipts
7135  real(dp),parameter :: c2_27=2.0e0_dp/27.0e0_dp,c5_9=5.0e0_dp/9.0e0_dp
7136  real(dp),parameter :: c8_9=8.0e0_dp/9.0e0_dp,m1thrd=-third
7137  real(dp) :: del,prod,rho,rhom1,rhomtrd
7138  logical :: test
7139 !character(len=500) :: message
7140 
7141 ! *************************************************************************
7142 
7143 !Loop over points : here, brute force algorithm
7144 !do ipts=1,npts
7145 !rspts(ipts)=sign( (abs(rhoarr(ipts)))**m1thrd,rhoarr(ipts))
7146 !end do
7147 !
7148 
7149  rhomtrd=sign( (abs(rhoarr(1)))**m1thrd, rhoarr(1) )
7150  rhom1=one/rhoarr(1)
7151  rspts(1)=rhomtrd
7152  do ipts=2,npts
7153    rho=rhoarr(ipts)
7154    prod=rho*rhom1
7155 !  If the previous point is too far ...
7156    if(prod < 0.01_dp .or. prod > 10._dp )then
7157      rhomtrd=sign( (abs(rho))**m1thrd , rho )
7158      rhom1=one/rho
7159    else
7160      del=prod-one
7161      do ii=1,5
7162 !      Choose one of the two next lines, the last one is more accurate
7163 !      rhomtrd=((one+third*del)/(one+two_thirds*del))*rhomtrd
7164        rhomtrd=((one+c5_9*del)/(one+del*(c8_9+c2_27*del)))*rhomtrd
7165        rhom1=rhomtrd*rhomtrd*rhomtrd
7166        del=rho*rhom1-one
7167 !      write(std_out,*)rhomtrd,del
7168        test = del*del < 1.0e-24_dp
7169        if(test) exit
7170      end do
7171      if( .not. test) then
7172        rhomtrd=sign( (abs(rho))**m1thrd , rho )
7173      end if
7174    end if
7175    rspts(ipts)=rhomtrd
7176  end do
7177 
7178 end subroutine invcb

m_numeric_tools/is_integer_0D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  is_integer_0D

FUNCTION

  Return .TRUE. if all elements differ from an integer by less that tol

INPUTS

  rr=the set of real values to be checked
  tol=tolerance on the difference between real and integer

SOURCE

1611 pure function is_integer_0d(rr,tol) result(ans)
1612 
1613 
1614 !This section has been created automatically by the script Abilint (TD).
1615 !Do not modify the following lines by hand.
1616 #undef ABI_FUNC
1617 #define ABI_FUNC 'is_integer_0d'
1618 !End of the abilint section
1619 
1620  implicit none
1621 
1622 !Arguments ------------------------------------
1623 !scalars
1624  real(dp),intent(in) :: tol
1625  logical :: ans
1626 !arrays
1627  real(dp),intent(in) :: rr
1628 
1629 ! *************************************************************************
1630 
1631  ans=(ABS(rr-NINT(rr))<tol)
1632 
1633 end function is_integer_0d

m_numeric_tools/is_integer_1D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  is_integer_1D

FUNCTION

INPUTS

OUTPUT

SOURCE

1650 pure function is_integer_1d(rr,tol) result(ans)
1651 
1652 
1653 !This section has been created automatically by the script Abilint (TD).
1654 !Do not modify the following lines by hand.
1655 #undef ABI_FUNC
1656 #define ABI_FUNC 'is_integer_1d'
1657 !End of the abilint section
1658 
1659  implicit none
1660 
1661 !Arguments ------------------------------------
1662 !scalars
1663  real(dp),intent(in) :: tol
1664  logical :: ans
1665 !arrays
1666  real(dp),intent(in) :: rr(:)
1667 
1668 ! *************************************************************************
1669 
1670  ans=ALL((ABS(rr-NINT(rr))<tol))
1671 
1672 end function is_integer_1d

m_numeric_tools/is_zero_rdp_0D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  is_zero_rdp_0D

FUNCTION

  Return .TRUE. if all elements differ from zero by less that tol

INPUTS

  rr=the set of real values to be checked
  tol=tolerance

OUTPUT

PARENTS

CHILDREN

SOURCE

1696 function is_zero_rdp_0d(rr,tol) result(ans)
1697 
1698 
1699 !This section has been created automatically by the script Abilint (TD).
1700 !Do not modify the following lines by hand.
1701 #undef ABI_FUNC
1702 #define ABI_FUNC 'is_zero_rdp_0d'
1703 !End of the abilint section
1704 
1705  implicit none
1706 
1707 !Arguments ------------------------------------
1708 !scalars
1709  real(dp),intent(in) :: tol
1710  logical :: ans
1711 !arrays
1712  real(dp),intent(in) :: rr
1713 ! *************************************************************************
1714 
1715  ans=(ABS(rr)<tol)
1716 
1717 end function is_zero_rdp_0d

m_numeric_tools/is_zero_rdp_1d [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  is_zero_rdp_1d

FUNCTION

INPUTS

OUTPUT

SOURCE

1734 function is_zero_rdp_1d(rr,tol) result(ans)
1735 
1736 
1737 !This section has been created automatically by the script Abilint (TD).
1738 !Do not modify the following lines by hand.
1739 #undef ABI_FUNC
1740 #define ABI_FUNC 'is_zero_rdp_1d'
1741 !End of the abilint section
1742 
1743  implicit none
1744 
1745 !Arguments ------------------------------------
1746 !scalars
1747  real(dp),intent(in) :: tol
1748  logical :: ans
1749 !arrays
1750  real(dp),intent(in) :: rr(:)
1751 ! *************************************************************************
1752 
1753  ans=ALL(ABS(rr(:))<tol)
1754 
1755 end function is_zero_rdp_1d

m_numeric_tools/isdiagmat_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  isdiagmat_int

FUNCTION

  True if matrix mat is diagonal

SOURCE

935 pure logical function isdiagmat_int(mat) result(ans)
936 
937 
938 !This section has been created automatically by the script Abilint (TD).
939 !Do not modify the following lines by hand.
940 #undef ABI_FUNC
941 #define ABI_FUNC 'isdiagmat_int'
942 !End of the abilint section
943 
944  implicit none
945 
946 !Arguments ------------------------------------
947 !scalars
948  integer,intent(in) :: mat(:,:)
949 
950 !Local variables-------------------------------
951  integer :: ii,jj
952 ! *************************************************************************
953 
954  ans = .True.
955  do jj=1,size(mat,dim=2)
956    do ii=1,size(mat,dim=1)
957      if (ii == jj) cycle
958      if (mat(ii,jj) /= 0) then
959        ans = .False.; return
960      end if
961    end do
962  end do
963 
964 end function isdiagmat_int

m_numeric_tools/isdiagmat_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  isdiagmat_rdp

FUNCTION

  True if matrix mat is diagonal withing the given absolute tolerance (default: tol12)

SOURCE

 978 pure logical function isdiagmat_rdp(mat, atol) result(ans)
 979 
 980 
 981 !This section has been created automatically by the script Abilint (TD).
 982 !Do not modify the following lines by hand.
 983 #undef ABI_FUNC
 984 #define ABI_FUNC 'isdiagmat_rdp'
 985 !End of the abilint section
 986 
 987  implicit none
 988 
 989 !This section has been created automatically by the script Abilint (TD).
 990 !Do not modify the following lines by hand.
 991 #undef ABI_FUNC
 992 #define ABI_FUNC 'isdiagmat_rdp'
 993 !End of the abilint section
 994 
 995 !Arguments ------------------------------------
 996 !scalars
 997  real(dp),intent(in) :: mat(:,:)
 998  real(dp),optional,intent(in) :: atol
 999 
1000 !Local variables-------------------------------
1001  integer :: ii,jj
1002  real(dp) :: my_atol
1003 ! *************************************************************************
1004 
1005  my_atol = tol12; if (present(atol)) my_atol = atol
1006 
1007  ans = .True.
1008  do jj=1,size(mat,dim=2)
1009    do ii=1,size(mat,dim=1)
1010      if (ii == jj) cycle
1011      if (abs(mat(ii,jj)) > my_atol) then
1012        ans = .False.; return
1013      end if
1014    end do
1015  end do
1016 
1017 end function isdiagmat_rdp

m_numeric_tools/iseven [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  iseven

FUNCTION

  Return .TRUE. if the given integer is even

PARENTS

CHILDREN

SOURCE

1574 elemental function iseven(nn)
1575 
1576 
1577 !This section has been created automatically by the script Abilint (TD).
1578 !Do not modify the following lines by hand.
1579 #undef ABI_FUNC
1580 #define ABI_FUNC 'iseven'
1581 !End of the abilint section
1582 
1583  implicit none
1584 
1585 !Arguments ------------------------------------
1586 !scalars
1587  integer,intent(in) :: nn
1588  logical :: iseven
1589 ! *********************************************************************
1590 
1591  iseven=.FALSE.; if ((nn/2)*2==nn) iseven=.TRUE.
1592 
1593 end function iseven

m_numeric_tools/isinside [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  isinside

FUNCTION

  True if float `xval` is inside the interval [win(1), win(2)]

SOURCE

1769 pure logical function isinside(xval, win)
1770 
1771 
1772 !This section has been created automatically by the script Abilint (TD).
1773 !Do not modify the following lines by hand.
1774 #undef ABI_FUNC
1775 #define ABI_FUNC 'isinside'
1776 !End of the abilint section
1777 
1778  implicit none
1779 
1780 !Arguments ------------------------------------
1781 !scalars
1782  real(dp),intent(in) :: xval,win(2)
1783 ! *************************************************************************
1784 
1785  isinside = (xval >= win(1) .and. xval <= win(2))
1786 
1787 end function isinside

m_numeric_tools/isordered [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  isordered

FUNCTION

  Return .TRUE. if values in array arr are ordered.
  Consider that two double precision numbers within tolerance tol are equal.

INPUTS

  nn=Size of arr.
  arr(nn)=The array with real values to be tested.
  direction= ">" for ascending numerical order.
             ">" for decreasing numerical order.

PARENTS

CHILDREN

SOURCE

5449 function isordered_rdp(nn,arr,direction,tol) result(isord)
5450 
5451 
5452 !This section has been created automatically by the script Abilint (TD).
5453 !Do not modify the following lines by hand.
5454 #undef ABI_FUNC
5455 #define ABI_FUNC 'isordered_rdp'
5456 !End of the abilint section
5457 
5458  implicit none
5459 
5460 !Arguments ------------------------------------
5461 !scalars
5462  integer,intent(in) :: nn
5463  real(dp),intent(in) :: tol
5464  logical :: isord
5465  character(len=*),intent(in) :: direction
5466 !arrays
5467  real(dp),intent(in) :: arr(nn)
5468 
5469 !Local variables ------------------------------
5470 !scalars
5471  integer :: ii
5472  real(dp) :: prev
5473  character(len=500) :: msg
5474 
5475 ! *************************************************************************
5476 
5477  prev = arr(1); isord =.TRUE.
5478 
5479  SELECT CASE (direction(1:1))
5480  CASE(">")
5481  ii=2;
5482  do while (ii<=nn .and. isord)
5483    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) >= prev)
5484    prev = arr(ii)
5485    ii = ii +1
5486  end do
5487 
5488  CASE("<")
5489  ii=2;
5490  do while (ii<=nn .and. isord)
5491    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) <= prev)
5492    prev = arr(ii)
5493    ii = ii +1
5494  end do
5495 
5496  CASE DEFAULT
5497    msg = "Wrong direction: "//TRIM(direction)
5498    MSG_ERROR(msg)
5499  END SELECT
5500 
5501 end function isordered_rdp

m_numeric_tools/kramerskronig [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 kramerskronig

FUNCTION

 check or apply the Kramers Kronig relation:
  Re \epsilon(\omega) = 1 + \frac{2}{\pi}
  \int_0^\infty d\omega' frac{\omega'}{\omega'^2 - \omega^2} Im \epsilon(\omega')

INPUTS

  nomega=number of real frequencies
  omega(nomega)= real frequencies
  eps(nomega)= function on the frequency grid (both real and imaginary part)
   real part can be used to check whether the K-K relation is satisfied or not
  method=method used to perform the integration
   0= naive integration
   1=simpson rule
  only_check= if /=0 the real part of eps is checked against the imaginary part,
                a final report in written but the array eps is not modified
              if ==0 the real part of eps is overwritten using the
              results obtained using the Kramers-Kronig relation

OUTPUT

NOTES

 Inspired to check_kramerskronig of the DP code

PARENTS

      linear_optics_paw

CHILDREN

      simpson_int,wrtout

SOURCE

6890 subroutine kramerskronig(nomega,omega,eps,method,only_check)
6891 
6892 
6893 !This section has been created automatically by the script Abilint (TD).
6894 !Do not modify the following lines by hand.
6895 #undef ABI_FUNC
6896 #define ABI_FUNC 'kramerskronig'
6897 !End of the abilint section
6898 
6899  implicit none
6900 
6901 !Arguments ------------------------------------
6902 !scalars
6903  integer,intent(in) :: method,nomega,only_check
6904 !arrays
6905  real(dp),intent(in) :: omega(nomega)
6906  complex,intent(inout) :: eps(nomega)
6907 
6908 !Local variables-------------------------------
6909 !scalars
6910  integer,save :: enough=0
6911  integer :: ii,ip
6912  real(dp) :: acc,domega,eav,kkdif,kkrms,ww,wwp
6913  character(len=500) :: msg
6914 !arrays
6915  real(dp) :: e1kk(nomega),intkk(nomega),kk(nomega)
6916 
6917 ! *************************************************************************
6918 
6919 !Check whether the frequency grid is linear or not
6920  domega = (omega(nomega) - omega(1)) / (nomega-1)
6921  do ii=2,nomega
6922    if (ABS(domega-(omega(ii)-omega(ii-1))) > 0.001) then
6923      if (only_check/=1) then
6924        msg="check cannot be performed since the frequency step is not constant"
6925        MSG_WARNING(msg)
6926        RETURN
6927      else
6928        msg=' Cannot perform integration since frequency step is not constant'
6929        MSG_ERROR(msg)
6930      end if
6931    end if
6932  end do
6933 
6934 !Check whether omega(1) is small or not
6935  if (omega(1) > 0.1/Ha_eV) then
6936    if (only_check/=1) then
6937      msg=' Check cannot be performed since first frequency on the grid > 0.1 eV'
6938      MSG_WARNING(msg)
6939      RETURN
6940    else
6941      msg=' Cannot perform integration since first frequency on the grid > 0.1 eV'
6942      MSG_ERROR(msg)
6943    end if
6944  end if
6945 
6946 !If eps(nomega) is not 0 warn
6947  if (AIMAG(eps(nomega)) > 0.1 .and. enough<50) then
6948    enough=enough+1
6949    write(msg,'(a,f8.4,3a,f8.2,2a)')&
6950 &   'Im epsilon for omega = ',omega(nomega)*Ha_eV,' eV',ch10,&
6951 &   'is not yet zero, epsilon_2 = ',AIMAG(eps(nomega)),ch10,&
6952 &   'Kramers Kronig could give wrong results'
6953    MSG_WARNING(msg)
6954    if (enough==50) then
6955      write(msg,'(3a)')' sufficient number of WARNINGS-',ch10,' stop writing '
6956      call wrtout(std_out,msg,'COLL')
6957    end if
6958  end if
6959 
6960 
6961 !Perform Kramers-Kronig using naive integration
6962  select case (method)
6963  case (0)
6964 
6965    do ii=1,nomega
6966      ww = omega(ii)
6967      acc = 0.0_dp
6968      do ip=1,nomega
6969        if (ip == ii) CYCLE
6970        wwp = omega(ip)
6971        acc = acc + wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6972      end do
6973      e1kk(ii) = one + two/pi*domega* acc
6974    end do
6975 
6976 !    Perform Kramers-Kronig using Simpson integration
6977 !    Simpson O(1/N^4), from NumRec in C p 134  NumRec in Fortran p 128
6978  case (1)
6979 
6980    kk=zero
6981 
6982    do ii=1,nomega
6983      ww=omega(ii)
6984      do ip=1,nomega
6985        if (ip == ii) CYCLE
6986        wwp = omega(ip)
6987        kk(ip) = wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6988      end do
6989      call simpson_int(nomega,domega,kk,intkk)
6990      e1kk(ii) = one + two/pi * intkk(nomega)
6991    end do
6992 
6993  case default
6994    write(msg,'(a,i0)')' Wrong value for method ',method
6995    MSG_BUG(msg)
6996  end select
6997 
6998 !at this point real part is in e1kk, need to put it into eps
6999  do ii=1,nomega
7000    eps(ii)=CMPLX(e1kk(ii),AIMAG(eps(ii)))
7001  end do
7002 
7003 !Verify Kramers-Kronig
7004  eav   = zero
7005  kkdif = zero
7006  kkrms = zero
7007 
7008  do ii=1,nomega
7009    kkdif = kkdif + ABS(REAL(eps(ii)) - e1kk(ii))
7010    kkrms = kkrms + (REAL(eps(ii)) - e1kk(ii))*(REAL(eps(ii)) - e1kk(ii))
7011    eav = eav + ABS(REAL(eps(ii)))
7012  end do
7013 
7014  eav = eav/nomega
7015  kkdif = (kkdif/nomega) / eav
7016  kkrms = (kkrms/nomega) / (eav*eav)
7017 
7018  kk = ABS(REAL(eps(1)) - e1kk(1)) / REAL(eps(1))
7019 
7020 !Write data
7021  write(msg,'(a,f7.2,a)')' Kramers-Kronig transform is verified within ',MAXVAL(kk)*100,"%"
7022  call wrtout(std_out,msg,'COLL')
7023 
7024 end subroutine kramerskronig

m_numeric_tools/l2int_1D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  l2int_1D

FUNCTION

  Convert a logical array into an int array (True --> 1, False --> 0)

INPUTS

  larr(:)=the input logical array

SOURCE

1034 pure function l2int_1D(larr) result(int_arr)
1035 
1036 
1037 !This section has been created automatically by the script Abilint (TD).
1038 !Do not modify the following lines by hand.
1039 #undef ABI_FUNC
1040 #define ABI_FUNC 'l2int_1D'
1041 !End of the abilint section
1042 
1043  implicit none
1044 
1045 !Arguments ------------------------------------
1046 !scalars
1047  logical,intent(in) :: larr(:)
1048  integer :: int_arr(size(larr))
1049 
1050 ! *********************************************************************
1051 
1052  where (larr)
1053    int_arr = 1
1054  elsewhere
1055    int_arr = 0
1056  end where
1057 
1058 end function l2int_1D

m_numeric_tools/l2int_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  l2int_2D

FUNCTION

  Convert a logical array into an int array (True --> 1, False --> 0)

INPUTS

  larr(:)=the input logical array

SOURCE

1075 pure function l2int_2D(larr) result(int_arr)
1076 
1077 
1078 !This section has been created automatically by the script Abilint (TD).
1079 !Do not modify the following lines by hand.
1080 #undef ABI_FUNC
1081 #define ABI_FUNC 'l2int_2D'
1082 !End of the abilint section
1083 
1084  implicit none
1085 
1086 !Arguments ------------------------------------
1087 !scalars
1088  logical,intent(in) :: larr(:,:)
1089  integer :: int_arr(size(larr,1), size(larr,2))
1090 
1091 ! *********************************************************************
1092 
1093  where (larr)
1094    int_arr = 1
1095  elsewhere
1096    int_arr = 0
1097  end where
1098 
1099 end function l2int_2D

m_numeric_tools/l2int_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  l2int_3D

FUNCTION

  Convert a logical array into an int array (True --> 1, False --> 0)

INPUTS

  larr(:)=the input logical array

SOURCE

1116 pure function l2int_3D(larr) result(int_arr)
1117 
1118 
1119 !This section has been created automatically by the script Abilint (TD).
1120 !Do not modify the following lines by hand.
1121 #undef ABI_FUNC
1122 #define ABI_FUNC 'l2int_3D'
1123 !End of the abilint section
1124 
1125  implicit none
1126 
1127 !Arguments ------------------------------------
1128 !scalars
1129  logical,intent(in) :: larr(:,:,:)
1130  integer :: int_arr(size(larr,1), size(larr,2), size(larr,3))
1131 
1132 ! *********************************************************************
1133 
1134  where (larr)
1135    int_arr = 1
1136  elsewhere
1137    int_arr = 0
1138  end where
1139 
1140 end function l2int_3D

m_numeric_tools/l2norm_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  l2norm_rdp

FUNCTION

  Return the length (ordinary L2 norm) of a vector.


m_numeric_tools/lfind [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  lfind

FUNCTION

  Find the index of the first occurrence of .True. in a logical array.
  Return -1 if not found. If back is True, the search starts from the
  last element of the array (default: False).

INPUTS

  mask(:)=Input logical mask

PARENTS

CHILDREN

SOURCE

2111 integer pure function lfind(mask, back)
2112 
2113 
2114 !This section has been created automatically by the script Abilint (TD).
2115 !Do not modify the following lines by hand.
2116 #undef ABI_FUNC
2117 #define ABI_FUNC 'lfind'
2118 !End of the abilint section
2119 
2120  implicit none
2121 
2122 !Arguments ------------------------------------
2123 !scalars
2124  logical,intent(in) :: mask(:)
2125  logical,optional,intent(in) :: back
2126 !arrays
2127 
2128 !Local variables-------------------------------
2129 !scalars
2130  integer :: ii,nitems
2131  logical :: do_back
2132 
2133 !************************************************************************
2134 
2135  do_back = .False.; if (present(back)) do_back = back
2136  lfind = -1; nitems = size(mask); if (nitems == 0) return
2137 
2138  if (do_back) then
2139    ! Backward search
2140    do ii=nitems,1,-1
2141      if (mask(ii)) then
2142        lfind = ii; return
2143      end if
2144    end do
2145  else
2146    ! Forward search.
2147    do ii=1,nitems
2148      if (mask(ii)) then
2149        lfind = ii; return
2150      end if
2151    end do
2152  end if
2153 
2154 end function lfind

m_numeric_tools/linfit_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linfit_dpc

FUNCTION

INPUTS

OUTPUT

SOURCE

2488 function linfit_dpc(nn,xx,zz,aa,bb) result(res)
2489 
2490 
2491 !This section has been created automatically by the script Abilint (TD).
2492 !Do not modify the following lines by hand.
2493 #undef ABI_FUNC
2494 #define ABI_FUNC 'linfit_dpc'
2495 !End of the abilint section
2496 
2497  implicit none
2498 
2499 !Arguments ------------------------------------
2500 !scalars
2501  integer,intent(in) :: nn
2502  real(dp) :: res
2503  real(dp),intent(in) :: xx(nn)
2504  complex(dpc),intent(in) :: zz(nn)
2505  complex(dpc),intent(out) :: aa,bb
2506 !arrays
2507 
2508 !Local variables-------------------------------
2509 !scalars
2510  integer :: ii
2511  real(dp) :: sx,sx2,msrt
2512  complex(dpc) :: sz,sxz
2513 ! *************************************************************************
2514 
2515  sx=zero  ; sx2=zero ; msrt=zero
2516  sz=czero ; sxz=czero
2517  do ii=1,nn
2518   sx=sx+xx(ii)
2519   sz=sz+zz(ii)
2520   sxz=sxz+xx(ii)*zz(ii)
2521   sx2=sx2+xx(ii)*xx(ii)
2522  end do
2523 
2524  aa=(nn*sxz-sx*sz)/(nn*sx2-sx*sx)
2525  bb=sz/nn-sx*aa/nn
2526 
2527  do ii=1,nn
2528   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2529  end do
2530  msrt=SQRT(msrt) ; res=msrt
2531 
2532 end function linfit_dpc

m_numeric_tools/linfit_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linfit_rdp

FUNCTION

  Perform a linear fit, y=ax+b, of data

INPUTS

  xx(nn)=xx coordinates
  yy(nn)=yy coordinates

OUTPUT

  aa=coefficient of linear term of fit
  bb=coefficient of constant term of fit
  function linfit=root mean square of differences between data and fit

SOURCE

2360 function linfit_rdp(nn,xx,yy,aa,bb) result(res)
2361 
2362 
2363 !This section has been created automatically by the script Abilint (TD).
2364 !Do not modify the following lines by hand.
2365 #undef ABI_FUNC
2366 #define ABI_FUNC 'linfit_rdp'
2367 !End of the abilint section
2368 
2369  implicit none
2370 
2371 !This section has been created automatically by the script Abilint (TD).
2372 !Do not modify the following lines by hand.
2373 #undef ABI_FUNC
2374 #define ABI_FUNC 'linfit_rdp'
2375 !End of the abilint section
2376 
2377 !Arguments ------------------------------------
2378 !scalars
2379  integer,intent(in) :: nn
2380  real(dp) :: res
2381  real(dp),intent(out) :: aa,bb
2382 !arrays
2383  real(dp),intent(in) :: xx(nn),yy(nn)
2384 
2385 !Local variables-------------------------------
2386 !scalars
2387  integer :: ii
2388  real(dp) :: msrt,sx2,sx,sxy,sy,tx,ty
2389 ! *************************************************************************
2390 
2391  sx=zero ; sy=zero ; sxy=zero ; sx2=zero
2392  do ii=1,nn
2393   tx=xx(ii)
2394   ty=yy(ii)
2395   sx=sx+tx
2396   sy=sy+ty
2397   sxy=sxy+tx*ty
2398   sx2=sx2+tx*tx
2399  end do
2400 
2401  aa=(nn*sxy-sx*sy)/(nn*sx2-sx*sx)
2402  bb=sy/nn-sx*aa/nn
2403 
2404  msrt=zero
2405  do ii=1,nn
2406   tx=xx(ii)
2407   ty=yy(ii)
2408   msrt=msrt+(ty-aa*tx-bb)**2
2409  end do
2410  msrt=SQRT(msrt/nn) ; res=msrt
2411 
2412 end function linfit_rdp

m_numeric_tools/linfit_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linfit_spc

FUNCTION

INPUTS

OUTPUT

SOURCE

2429 function linfit_spc(nn,xx,zz,aa,bb) result(res)
2430 
2431 !Arguments ------------------------------------
2432 !scalars
2433 
2434 !This section has been created automatically by the script Abilint (TD).
2435 !Do not modify the following lines by hand.
2436 #undef ABI_FUNC
2437 #define ABI_FUNC 'linfit_spc'
2438 !End of the abilint section
2439 
2440  integer,intent(in) :: nn
2441  real(dp) :: res
2442  real(dp),intent(in) :: xx(nn)
2443  complex(spc),intent(in) :: zz(nn)
2444  complex(spc),intent(out) :: aa,bb
2445 !arrays
2446 
2447 !Local variables-------------------------------
2448 !scalars
2449  integer :: ii
2450  real(dp) :: sx,sx2,msrt
2451  complex(dpc) :: sz,sxz
2452 ! *************************************************************************
2453 
2454  sx=zero ; sx2=zero ; msrt=zero
2455  sz=czero ; sxz=czero
2456  do ii=1,nn
2457   sx=sx+xx(ii)
2458   sz=sz+zz(ii)
2459   sxz=sxz+xx(ii)*zz(ii)
2460   sx2=sx2+xx(ii)*xx(ii)
2461  end do
2462 
2463  aa=(nn*sxz-sx*sz)/(nn*sx2-sx*sx)
2464  bb=sz/nn-sx*aa/nn
2465 
2466  do ii=1,nn
2467   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2468  end do
2469  msrt=SQRT(msrt) ; res=msrt
2470 
2471 end function linfit_spc

m_numeric_tools/list2blocks [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  list2blocks

FUNCTION

  Given a list of integers, find the number of contiguos groups of values.
  and returns the set of indices that can be used to loop over these groups
  Example list = [1,2,3,5,6] --> blocks = [[1,3], [4,5]]

INPUTS

  list(:)=List of integers

 OUTPUTS
 nblocks=Number of blocks
 blocks(2,nblocks)=
    allocatable array in input
    in output:
       blocks(1,i) gives the start of the i-th block
       blocks(2,i) gives the end of the i-th block

PARENTS

      m_wfd

CHILDREN

SOURCE

2186 subroutine list2blocks(list,nblocks,blocks)
2187 
2188 
2189 !This section has been created automatically by the script Abilint (TD).
2190 !Do not modify the following lines by hand.
2191 #undef ABI_FUNC
2192 #define ABI_FUNC 'list2blocks'
2193 !End of the abilint section
2194 
2195  implicit none
2196 
2197 !Arguments ------------------------------------
2198 !scalars
2199  integer,intent(out) :: nblocks
2200  integer,intent(in) :: list(:)
2201 !arrays
2202  integer,intent(out),allocatable :: blocks(:,:)
2203 
2204 !Local variables-------------------------------
2205 !scalars
2206  integer :: ii,nitems
2207 !arrays
2208  integer :: work(2,size(list))
2209 
2210 !************************************************************************
2211 
2212  nitems = size(list)
2213 
2214  ! Handle nitems == 1 case
2215  if (nitems == 1) then
2216    ABI_MALLOC(blocks, (2,1))
2217    blocks = 1
2218    return
2219  end if
2220 
2221  nblocks = 1; work(1,1) = 1
2222 
2223  do ii=2,nitems
2224    if (list(ii) /= (list(ii-1) + 1)) then
2225      work(2,nblocks) = ii - 1
2226      nblocks = nblocks + 1
2227      work(1,nblocks) = ii
2228    end if
2229  end do
2230 
2231  work(2,nblocks) = nitems
2232 
2233  ABI_MALLOC(blocks, (2,nblocks))
2234  blocks = work(:,1:nblocks)
2235 
2236 end subroutine list2blocks

m_numeric_tools/llsfit_svd [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  llsfit_svd

FUNCTION

  Given a set of N data points (x,y) with individual standard deviations sigma_i,
  use chi-square minimization to determine the M coefficients, par, of a function that
  depends linearly on nfuncs functions, i.e f(x) = \sum_i^{nfuncs} par_i * func_i(x).
  Solve the fitting equations using singular value decomposition of the design matrix as in Eq 14.3.17
  of Numerical Recipes. The program returns values for the M fit parameters par, and chi-square.
  The user supplies a subroutine funcs(x,nfuncs) that returns the M basis functions evaluated at xx.

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

2559 subroutine llsfit_svd(xx,yy,sigma,nfuncs,funcs,chisq,par,var,cov,info)
2560 
2561 
2562 !This section has been created automatically by the script Abilint (TD).
2563 !Do not modify the following lines by hand.
2564 #undef ABI_FUNC
2565 #define ABI_FUNC 'llsfit_svd'
2566 !End of the abilint section
2567 
2568  implicit none
2569 
2570 !Arguments ------------------------------------
2571 !scalars
2572  integer,intent(in) :: nfuncs
2573  integer,intent(out) :: info
2574  real(dp),intent(out) :: chisq
2575 !arrays
2576  real(dp),intent(in) :: xx(:),yy(:),sigma(:)
2577  real(dp),intent(out) :: par(:),var(:),cov(:,:)
2578 
2579  interface
2580   function funcs(xx,nf)
2581   use defs_basis
2582   implicit none
2583   real(dp),intent(in) :: xx
2584   integer,intent(in) :: nf
2585   real(dp) :: funcs(nf)
2586   end function funcs
2587  end interface
2588 
2589 !Local variables-------------------------------
2590  integer,parameter :: PAD_=50
2591  integer :: ii,npts,lwork
2592  real(dp),parameter :: TOL_=1.0e-5_dp
2593 !arrays
2594  real(dp),dimension(SIZE(xx)) :: bb,sigm1
2595  real(dp),dimension(SIZE(xx),nfuncs) :: dmat,dmat_save
2596  real(dp) :: tmp(nfuncs)
2597  real(dp),allocatable :: work(:),Vt(:,:),U(:,:),S(:)
2598 ! *************************************************************************
2599 
2600  npts=assert_eq(SIZE(xx),SIZE(yy),SIZE(sigma),'Wrong size in xx,yy,sigma',&
2601 & __FILE__,__LINE__)
2602  call assert((npts>=nfuncs),'No. of functions must greater than no. of points',&
2603 & __FILE__,__LINE__)
2604  ii=assert_eq(nfuncs,SIZE(cov,1),SIZE(cov,2),SIZE(var),'Wrong size in covariance',&
2605 & __FILE__,__LINE__)
2606 
2607  !
2608  ! === Calculate design matrix and b vector ===
2609  ! * dmat_ij=f_j(x_i)/sigma_i, b_i=y_i/sigma_i
2610  sigm1(:)=one/sigma(:) ; bb(:)=yy(:)*sigm1(:)
2611  do ii=1,npts
2612   dmat_save(ii,:)=funcs(xx(ii),nfuncs)
2613  end do
2614  dmat=dmat_save*SPREAD(sigm1,DIM=2,ncopies=nfuncs)
2615  dmat_save(:,:)=dmat(:,:)
2616  !
2617  ! === Singular value decomposition ===
2618  lwork=MAX(3*MIN(npts,nfuncs)+MAX(npts,nfuncs),5*MIN(npts,nfuncs)-4)+PAD_
2619  ABI_MALLOC(work,(lwork))
2620  ABI_MALLOC(U,(npts,npts))
2621  ABI_MALLOC(S,(nfuncs))
2622  ABI_MALLOC(Vt,(nfuncs,nfuncs))
2623 
2624  call DGESVD('A','A',npts,nfuncs,dmat,npts,S,U,npts,Vt,nfuncs,work,lwork,info)
2625  ABI_FREE(work)
2626  GOTO 10
2627  !
2628  ! === Set to zero small singular values according to TOL_ and find coefficients ===
2629  WHERE (S>TOL_*MAXVAL(S))
2630   tmp=MATMUL(bb,U)/S
2631  ELSEWHERE
2632   S  =zero
2633   tmp=zero
2634  END WHERE
2635  par(:)=MATMUL(tmp,Vt)
2636  !
2637  ! === Evaluate chi-square ===
2638  chisq=l2norm(MATMUL(dmat_save,par)-bb)**2
2639  !
2640  ! === Calculate covariance and variance ===
2641  ! C_jk = V_ji V_ki / S_i^2
2642  WHERE (S/=zero) S=one/(S*S)
2643 
2644  ! check this but should be correct
2645  cov(:,:)=Vt*SPREAD(S,DIM=2,ncopies=nfuncs)
2646  cov(:,:)=MATMUL(TRANSPOSE(Vt),cov)
2647  var(:)=SQRT(get_diag(cov))
2648 
2649 10 continue
2650  ABI_FREE(U)
2651  ABI_FREE(S)
2652  ABI_FREE(Vt)
2653 
2654 end subroutine llsfit_svd

m_numeric_tools/mask2blocks [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  mask2blocks

FUNCTION

  Give a logical mask, find the number of contiguos groups of .TRUE. values.
  and return the set of indices that can be used to loop over these groups

INPUTS

  mask(:)=Input logical mask

 OUTPUTS
 nblocks=Number of blocks

SIDE EFFECTS

 blocks(:,:)= Null pointer in input. blocks(2,nblocks) in output where
   blocks(1,i) gives the start of the i-th block
   blocks(2,i) gives the end of the i-th block

PARENTS

      m_wfk

CHILDREN

SOURCE

2267 subroutine mask2blocks(mask,nblocks,blocks)
2268 
2269 
2270 !This section has been created automatically by the script Abilint (TD).
2271 !Do not modify the following lines by hand.
2272 #undef ABI_FUNC
2273 #define ABI_FUNC 'mask2blocks'
2274 !End of the abilint section
2275 
2276  implicit none
2277 
2278 !Arguments ------------------------------------
2279 !scalars
2280  integer,intent(out) :: nblocks
2281  logical,intent(in) :: mask(:)
2282 !arrays
2283  integer,allocatable :: blocks(:,:)
2284 
2285 !Local variables-------------------------------
2286 !scalars
2287  integer :: ii,nitems,start
2288  logical :: inblock
2289 !arrays
2290  integer :: work(2,SIZE(mask))
2291 
2292 !************************************************************************
2293 
2294  ! Find first element.
2295  nitems = size(mask); start = 0
2296  do ii=1,nitems
2297    if (mask(ii)) then
2298      start = ii
2299      exit
2300    end if
2301  end do
2302 
2303  ! Handle no true element or just one.
2304  if (start == 0) then
2305    nblocks = 0
2306    ABI_MALLOC(blocks, (0,0))
2307    return
2308  end if
2309  if (start /= 0 .and. nitems == 1) then
2310    nblocks = 1
2311    ABI_MALLOC(blocks, (2,1))
2312    blocks(:,1) = [1,1]
2313  end if
2314 
2315  nblocks = 1; work(1,1) = start; inblock = .True.
2316 
2317  do ii=start+1,nitems
2318    if (.not.mask(ii)) then
2319      if (inblock) then
2320        inblock = .False.
2321        work(2,nblocks) = ii - 1
2322      end if
2323    else
2324      if (.not. inblock) then
2325        inblock = .True.
2326        nblocks = nblocks + 1
2327        work(1,nblocks) = ii
2328      end if
2329    end if
2330  end do
2331 
2332  if (mask(nitems) .and. inblock) work(2,nblocks) = nitems
2333 
2334  ABI_MALLOC(blocks, (2,nblocks))
2335  blocks = work(:,1:nblocks)
2336 
2337 end subroutine mask2blocks

m_numeric_tools/midpoint_ [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  midpoint_ (PRIVATE)

FUNCTION

  This routine computes the n-th stage of refinement of an extended midpoint rule.

INPUTS

  func(external)=the name of the function to be integrated
  xmin,xmax=the limits of integration
  nn=integer defining the refinement of the mesh, each call adds (2/3)*3n-1 additional
   interior points between xmin ans xmax

OUTPUT

  See SIDE EFFECTS

SIDE EFFECTS

  quad=the integral at the n-th stage.

PARENTS

CHILDREN

NOTES

  When called with nn=1, the routine returns as quad the crudest estimate of the integral
  Subsequent calls with nn=2,3,... (in that sequential order) will improve the accuracy of quad by adding
  (2/3)*3n-1 additional interior points. quad should not be modified between sequential calls.
  Subroutine is defined as recursive to allow multi-dimensional integrations

SOURCE

2888  recursive subroutine midpoint_(func,nn,xmin,xmax,quad)
2889 
2890 
2891 !This section has been created automatically by the script Abilint (TD).
2892 !Do not modify the following lines by hand.
2893 #undef ABI_FUNC
2894 #define ABI_FUNC 'midpoint_'
2895 !End of the abilint section
2896 
2897  implicit none
2898 
2899 !Arguments ------------------------------------
2900 !scalars
2901  integer,intent(in) :: nn
2902  !real(dp),external :: func
2903  real(dp),intent(in) :: xmin,xmax
2904  real(dp),intent(inout) :: quad
2905 
2906  interface
2907    function func(x)
2908      use defs_basis
2909      real(dp),intent(in) :: x
2910      real(dp) :: func
2911    end function func
2912  end interface
2913 
2914  !interface
2915  ! function func(x)
2916  !  use defs_basis
2917  !  real(dp),intent(in) :: x(:)
2918  !  real(dp) :: func(SIZE(x))
2919  ! end function func
2920  !end interface
2921 
2922 !Local variables-------------------------------
2923 !scalars
2924  integer  :: npt,ix
2925  real(dp) :: space
2926  character(len=500) :: msg
2927 !arrays
2928  real(dp),allocatable :: xx(:)
2929 
2930 !************************************************************************
2931 
2932  select case (nn)
2933 
2934  case (1)
2935    ! === Initial crude estimate done at the middle of the interval
2936    !quad=(xmax-xmin)*SUM(func((/half*(xmin+xmax)/))) !PARALLEL version
2937    quad=(xmax-xmin)*func(half*(xmin+xmax))
2938 
2939  case (2:)
2940    ! === Add npt interior points, they alternate in spacing between space and 2*space ===
2941    ABI_MALLOC(xx,(2*3**(nn-2)))
2942    npt=3**(nn-2) ; space=(xmax-xmin)/(three*npt)
2943    xx(1:2*npt-1:2)=arth(xmin+half*space,three*space,npt)
2944    xx(2:2*npt:2)=xx(1:2*npt-1:2)+two*space
2945    ! === The new sum is combined with the old integral to give a refined integral ===
2946    !quad=quad/three+space*SUM(func(xx))  !PARALLEL version
2947    quad=quad/three
2948    do ix=1,SIZE(xx)
2949      quad=quad+space*func(xx(ix))
2950    end do
2951    ABI_FREE(xx)
2952 
2953  case (:0)
2954    write(msg,'(a,i3)')' wrong value for nn ',nn
2955    MSG_BUG('Wrong value for nn')
2956  end select
2957 
2958 end subroutine midpoint_

m_numeric_tools/mincm [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  mincm

FUNCTION

  Return the minimum common multiple of ii and jj.

SOURCE

5100 integer function mincm(ii,jj)
5101 
5102 
5103 !This section has been created automatically by the script Abilint (TD).
5104 !Do not modify the following lines by hand.
5105 #undef ABI_FUNC
5106 #define ABI_FUNC 'mincm'
5107 !End of the abilint section
5108 
5109  implicit none
5110 
5111 !Arguments ------------------------------------
5112 !scalars
5113  integer,intent(in) :: ii,jj
5114 
5115 !************************************************************************
5116 
5117  if (ii==0.or.jj==0) then
5118    MSG_BUG('ii==0 or jj==0')
5119  end if
5120 
5121  mincm=MAX(ii,jj)
5122  do
5123    if ( ((mincm/ii)*ii)==mincm .and. ((mincm/jj)*jj)==mincm ) RETURN
5124    mincm=mincm+1
5125  end do
5126 
5127 end function mincm

m_numeric_tools/mkherm [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 mkherm

FUNCTION

 Make the complex array(ndim,ndim) hermitian,
 by adding half of it to its hermitian conjugate.

INPUTS

 ndim=dimension of the matrix
 array= complex matrix

SIDE EFFECTS

 array= hermitian matrix made by adding half of array to its hermitian conjugate

PARENTS

      anaddb,dfpt_phfrq,m_ddb,m_phgamma

CHILDREN

SOURCE

3903 pure subroutine mkherm(array,ndim)
3904 
3905 
3906 !This section has been created automatically by the script Abilint (TD).
3907 !Do not modify the following lines by hand.
3908 #undef ABI_FUNC
3909 #define ABI_FUNC 'mkherm'
3910 !End of the abilint section
3911 
3912  implicit none
3913 
3914 !Arguments -------------------------------
3915 !scalars
3916  integer,intent(in) :: ndim
3917 !arrays
3918  real(dp),intent(inout) :: array(2,ndim,ndim)
3919 
3920 !Local variables -------------------------
3921 !scalars
3922  integer :: i1,i2
3923 
3924 ! *********************************************************************
3925 
3926  do i1=1,ndim
3927    do i2=1,i1
3928      array(1,i1,i2)=(array(1,i1,i2)+array(1,i2,i1))*half
3929      array(2,i1,i2)=(array(2,i1,i2)-array(2,i2,i1))*half
3930      array(1,i2,i1)=array(1,i1,i2)
3931      array(2,i2,i1)=-array(2,i1,i2)
3932    end do
3933  end do
3934 
3935 end subroutine mkherm

m_numeric_tools/nderiv [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 nderiv

FUNCTION

 Given an input function y(x) on a regular grid,
 compute its first or second derivative.

INPUTS

  hh= radial step
  ndim= radial mesh size
  yy(ndim)= input function
  norder= order of derivation (1 or 2)

OUTPUT

  zz(ndim)= first or second derivative of y

PARENTS

      upf2abinit

CHILDREN

SOURCE

6373 subroutine nderiv(hh,yy,zz,ndim,norder)
6374 
6375 
6376 !This section has been created automatically by the script Abilint (TD).
6377 !Do not modify the following lines by hand.
6378 #undef ABI_FUNC
6379 #define ABI_FUNC 'nderiv'
6380 !End of the abilint section
6381 
6382  implicit none
6383 
6384 !Arguments ---------------------------------------------
6385 !scalars
6386  integer,intent(in) :: ndim,norder
6387  real(dp),intent(in) :: hh
6388 !arrays
6389  real(dp),intent(in) :: yy(ndim)
6390  real(dp),intent(out) :: zz(ndim)
6391 
6392 !Local variables ---------------------------------------
6393 !scalars
6394  integer :: ier,ii
6395  real(dp) :: aa,bb,cc,h1,y1
6396 
6397 ! *********************************************************************
6398 
6399 !Initialization (common to 1st and 2nd derivative)
6400  h1=one/(12.d0*hh)
6401  y1=yy(ndim-4)
6402 
6403 !FIRST DERIVATIVE
6404 !================
6405  if (norder==1) then
6406 
6407 !  Prepare differentiation loop
6408    bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
6409    cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
6410 !  Start differentiation loop
6411    do ii=5,ndim
6412      aa=bb;bb=cc
6413      cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
6414      zz(ii-4)=aa
6415    end do
6416 !  Normal exit
6417    ier=0
6418    aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
6419    zz(ndim)=h1*(3.d0*y1-16.d0*yy(ndim-3)+36.d0*yy(ndim-2) -48.d0*yy(ndim-1)+25.d0*yy(ndim))
6420    zz(ndim-1)=aa
6421    zz(ndim-2)=cc
6422    zz(ndim-3)=bb
6423 
6424 !  SECOND DERIVATIVE
6425 !  =================
6426  else
6427    h1=h1/hh
6428 !  Prepare differentiation loop
6429    bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
6430    cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
6431 !  Start differentiation loop
6432    do ii=5,ndim
6433      aa=bb;bb=cc
6434      cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
6435      zz(ii-4)=aa
6436    end do
6437 !  Normal exit
6438    ier=0
6439    aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
6440    zz(ndim)=h1*(11.d0*y1-56.d0*yy(ndim-3)+114.d0*yy(ndim-2) -104.d0*yy(ndim-1)+35.d0*yy(ndim))
6441    zz(ndim-1)=aa
6442    zz(ndim-2)=cc
6443    zz(ndim-3)=bb
6444 
6445  end if !norder
6446 
6447 end subroutine nderiv

m_numeric_tools/newrap_step [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  newrap_step

FUNCTION

  Apply single step newton-raphson method to find the root of a complex function
   z_k+1=z_k-f(z_k)/(df/dz(z_k))

SOURCE

4792 function newrap_step(z,f,df)
4793 
4794 
4795 !This section has been created automatically by the script Abilint (TD).
4796 !Do not modify the following lines by hand.
4797 #undef ABI_FUNC
4798 #define ABI_FUNC 'newrap_step'
4799 !End of the abilint section
4800 
4801  implicit none
4802 
4803 !Arguments ------------------------------------
4804 !scalars
4805  complex(dpc),intent(in) :: z,f,df
4806  complex(dpc) :: newrap_step
4807 
4808 !Local variables-------------------------------
4809 !scalars
4810  real(dp) :: dfm2
4811 ! *************************************************************************
4812 
4813  dfm2=ABS(df)*ABS(df)
4814 
4815  newrap_step= z - (f*CONJG(df))/dfm2
4816  !& z-one/(ABS(df)*ABS(df)) * CMPLX( REAL(f)*REAL(df)+AIMAG(f)*AIMAG(df), -REAL(f)*AIMAG(df)+AIMAG(f)*EAL(df) )
4817 
4818 end function newrap_step

m_numeric_tools/pack_matrix [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 pack_matrix

FUNCTION

 Packs a matrix into hermitian format

INPUTS

 N: size of matrix
 cplx: is the matrix complex
 mat_in(2, N*N)= matrix to be packed

OUTPUT

 mat_out(N*N+1)= packed matrix

PARENTS

      rayleigh_ritz

CHILDREN

SOURCE

4263 subroutine pack_matrix(mat_in, mat_out, N, cplx)
4264 
4265 
4266 !This section has been created automatically by the script Abilint (TD).
4267 !Do not modify the following lines by hand.
4268 #undef ABI_FUNC
4269 #define ABI_FUNC 'pack_matrix'
4270 !End of the abilint section
4271 
4272  implicit none
4273 
4274 !This section has been created automatically by the script Abilint (TD).
4275 !Do not modify the following lines by hand.
4276 #undef ABI_FUNC
4277 #define ABI_FUNC 'pack_matrix'
4278 !End of the abilint section
4279 
4280  integer, intent(in) :: N, cplx
4281  real(dp), intent(in) :: mat_in(cplx, N*N)
4282  real(dp), intent(out) :: mat_out(cplx*N*(N+1)/2)
4283  integer :: isubh, i, j
4284 
4285  ! *************************************************************************
4286 
4287  isubh = 1
4288  do j=1,N
4289    do i=1,j
4290      mat_out(isubh)    = mat_in(1, (j-1)*N+i)
4291      ! bad for vectorization, but it's not performance critical, so ...
4292      if(cplx == 2) then
4293        mat_out(isubh+1)  = mat_in(2, (j-1)*N+i)
4294      end if
4295      isubh=isubh+cplx
4296    end do
4297  end do
4298 
4299 end subroutine pack_matrix

m_numeric_tools/pade [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  pade

FUNCTION

  Calculate the pade approximant in zz of the function f calculated at the n points z

SOURCE

4610 function pade(n,z,f,zz)
4611 
4612 
4613 !This section has been created automatically by the script Abilint (TD).
4614 !Do not modify the following lines by hand.
4615 #undef ABI_FUNC
4616 #define ABI_FUNC 'pade'
4617 !End of the abilint section
4618 
4619  implicit none
4620 
4621 !Arguments ------------------------------------
4622 !scalars
4623  integer,intent(in) :: n
4624  complex(dpc),intent(in) :: zz
4625  complex(dpc) :: pade
4626 !arrays
4627  complex(dpc),intent(in) :: z(n),f(n)
4628 
4629 !Local variables-------------------------------
4630 !scalars
4631  complex(dpc) :: a(n)
4632  complex(dpc) :: Az(0:n), Bz(0:n)
4633  integer :: i
4634 ! *************************************************************************
4635 
4636  call calculate_pade_a(a,n,z,f)
4637 
4638  Az(0)=czero
4639  Az(1)=a(1)
4640  Bz(0)=cone
4641  Bz(1)=cone
4642 
4643  do i=1,n-1
4644    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4645    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4646  end do
4647  !write(std_out,*) 'Bz(n)',Bz(n)
4648  if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) ' Bz(n) ',Bz(n)
4649  pade=Az(n)/Bz(n)
4650  !write(std_out,*) 'pade_approx ', pade_approx
4651 
4652 end function pade

m_numeric_tools/pfactorize [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  pfactorize

FUNCTION

  Factorize a number in terms of an user-specified set of prime factors
  nn = alpha * Prod_i p^i   1)

INPUTS

  nn=The number to be factorized.
  nfactors=The number of factors
  pfactors(nfactors)=The list of prime number e.g. (/ 2, 3, 5, 7, 11 /)

OUTPUT

  powers(nfactors+1)=
   The first nfactors entries are the powers i in Eq.1. powers(nfactors+1) is alpha.

PARENTS

CHILDREN

SOURCE

5381 subroutine pfactorize(nn,nfactors,pfactors,powers)
5382 
5383 
5384 !This section has been created automatically by the script Abilint (TD).
5385 !Do not modify the following lines by hand.
5386 #undef ABI_FUNC
5387 #define ABI_FUNC 'pfactorize'
5388 !End of the abilint section
5389 
5390  implicit none
5391 
5392 !Arguments ------------------------------------
5393 !scalars
5394  integer,intent(in) :: nn,nfactors
5395  integer,intent(in) :: pfactors(nfactors)
5396  integer,intent(out) :: powers (nfactors+1)
5397 
5398 !Local variables ------------------------------
5399 !scalars
5400  integer :: tnn,ifc,fact,ipow,maxpwr
5401 
5402 ! *************************************************************************
5403 
5404  powers=0; tnn=nn
5405 
5406  fact_loop: do ifc=1,nfactors
5407    fact = pfactors (ifc)
5408    maxpwr = NINT ( LOG(DBLE(tnn))/LOG(DBLE(fact) ) ) + 1
5409    do ipow=1,maxpwr
5410      if (tnn==1) EXIT fact_loop
5411      if ( MOD(tnn,fact)==0 ) then
5412        tnn=tnn/fact
5413        powers(ifc)=powers(ifc) + 1
5414      end if
5415    end do
5416  end do fact_loop
5417 
5418  if ( nn /= tnn * PRODUCT( pfactors**powers(1:nfactors)) ) then
5419    MSG_BUG('nn/=tnn!')
5420  end if
5421 
5422  powers(nfactors+1) = tnn
5423 
5424 end subroutine pfactorize

m_numeric_tools/polyn_interp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  polyn_interp

FUNCTION

  Given arrays xa and ya of length N, and given a value x, return a value y, and an error estimate dy.
  If P(x) is the polynomial of degree N-1 such that P(xai)=yai, i=1,...,N, then the returned value y=P(x).

INPUTS

  xa(:)=abscissas in ascending order
  ya(:)=ordinates
  x=the point where the set of data has to be interpolated

OUTPUT

  y=the interpolated value
  dy=error estimate

NOTES

  Based on the polint routine reported in Numerical Recipies

PARENTS

      m_numeric_tools

CHILDREN

SOURCE

2686 subroutine polyn_interp(xa,ya,x,y,dy)
2687 
2688 
2689 !This section has been created automatically by the script Abilint (TD).
2690 !Do not modify the following lines by hand.
2691 #undef ABI_FUNC
2692 #define ABI_FUNC 'polyn_interp'
2693 !End of the abilint section
2694 
2695  implicit none
2696 
2697 !Arguments ------------------------------------
2698 !scalars
2699  real(dp),intent(in) :: xa(:),ya(:)
2700  real(dp),intent(in) :: x
2701  real(dp),intent(out) :: y,dy
2702 !Local variables-------------------------------
2703 !scalars
2704  integer :: m,n,ns
2705 !arrays
2706  real(dp),dimension(SIZE(xa)) :: c,d,den,ho
2707 ! *************************************************************************
2708 
2709  n=assert_eq(SIZE(xa),SIZE(ya),'Different size in xa and ya',&
2710 & __FILE__,__LINE__)
2711 
2712  ! === Initialize the tables of c and d ===
2713  c(:)=ya(:) ; d(:)=ya(:) ; ho(:)=xa(:)-x
2714  ! === Find closest table entry and initial approximation to y ===
2715  ns=imin_loc(ABS(x-xa)) ; y=ya(ns)
2716  ns=ns-1
2717  !
2718  ! === For each column of the tableau loop over current c and d and up-date them ===
2719  do m=1,n-1
2720   den(1:n-m)=ho(1:n-m)-ho(1+m:n)
2721   if (ANY(den(1:n-m)==zero)) then
2722    MSG_ERROR('Two input xa are identical')
2723   end if
2724 
2725   den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
2726   d(1:n-m)=ho(1+m:n)*den(1:n-m) ! Update c and d
2727   c(1:n-m)=ho(1:n-m)*den(1:n-m)
2728 
2729   if (2*ns<n-m) then  ! Now decide which correction, c or d, we want to add to the
2730    dy=c(ns+1)         ! accumulating value of y, The last dy added is the error indication.
2731   else
2732    dy=d(ns)
2733    ns=ns-1
2734   end if
2735 
2736   y=y+dy
2737  end do
2738 
2739 end subroutine polyn_interp

m_numeric_tools/print_arr1d_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr1d_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

4400 subroutine print_arr1d_dpc(arr,max_r,unit,mode_paral)
4401 
4402 
4403 !This section has been created automatically by the script Abilint (TD).
4404 !Do not modify the following lines by hand.
4405 #undef ABI_FUNC
4406 #define ABI_FUNC 'print_arr1d_dpc'
4407 !End of the abilint section
4408 
4409  implicit none
4410 
4411 !Arguments ------------------------------------
4412 !scalars
4413  integer,optional,intent(in) :: unit,max_r
4414  character(len=4),optional,intent(in) :: mode_paral
4415 !arrays
4416  complex(dpc),intent(in) :: arr(:)
4417 
4418 !Local variables-------------------------------
4419 !scalars
4420  integer :: unt,ii,nr,mr
4421  character(len=4) :: mode
4422  character(len=500) :: msg
4423  character(len=100) :: fmth,fmt1
4424 ! *************************************************************************
4425 
4426  unt=std_out ; if (PRESENT(unit      )) unt=unit
4427  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
4428  mr=15       ; if (PRESENT(max_r     )) mr=max_r
4429 
4430  if (mode/='COLL'.and.mode/='PERS') then
4431   write(msg,'(2a)')' Wrong value of mode_paral ',mode
4432   MSG_BUG(msg)
4433  end if
4434  !
4435  ! === Print out matrix ===
4436  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
4437 
4438  write(fmth,*)'(6x,',mr,'(i2,6x))'
4439  write(fmt1,*)'(3x,',mr,'f8.3)'
4440 
4441  write(msg,fmth)(ii,ii=1,mr)
4442  call wrtout(unt,msg,mode) !header
4443  write(msg,fmt1)REAL (arr(1:mr))
4444  call wrtout(unt,msg,mode) !real part
4445  write(msg,fmt1)AIMAG(arr(1:mr))
4446  call wrtout(unt,msg,mode) !imag part
4447 
4448 end subroutine print_arr1d_dpc

m_numeric_tools/print_arr1d_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr1d_spc

FUNCTION

 Print an array using a nice (?) format

INPUTS

  arr(:)=vector/matrix to be printed
  mode_paral(optional)=parallel mode, DEFAULT is "COLL"
   "COLL" if all procs are calling the routine with the same message to be written only once
   "PERS" if the procs are calling the routine with different mesgs each to be written,
          or if one proc is calling the routine
  unit(optional)=the unit number of the file, DEFAULT=std_out
  max_r,max_c(optional)=Max number of rows and columns to be printed
   (DEFAULT is 9, output format assumes to be less that 99, but there might be
    problems with wrtout if message size exceeds 500 thus max number of elements should be ~60)

OUTPUT

  (only printing)

PARENTS

CHILDREN

SOURCE

4331 subroutine print_arr1d_spc(arr,max_r,unit,mode_paral)
4332 
4333 
4334 !This section has been created automatically by the script Abilint (TD).
4335 !Do not modify the following lines by hand.
4336 #undef ABI_FUNC
4337 #define ABI_FUNC 'print_arr1d_spc'
4338 !End of the abilint section
4339 
4340  implicit none
4341 
4342 !Arguments ------------------------------------
4343 !scalars
4344  integer,optional,intent(in) :: unit,max_r
4345  character(len=4),optional,intent(in) :: mode_paral
4346 !arrays
4347  complex(spc),intent(in) :: arr(:)
4348 
4349 !Local variables-------------------------------
4350 !scalars
4351  integer :: unt,ii,nr,mr
4352  character(len=4) :: mode
4353  character(len=500) :: msg
4354  character(len=100) :: fmth,fmt1
4355 ! *************************************************************************
4356 
4357  unt=std_out ; if (PRESENT(unit      )) unt=unit
4358  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
4359  mr=15       ; if (PRESENT(max_r     )) mr=max_r
4360 
4361  if (mode/='COLL'.and.mode/='PERS') then
4362   write(msg,'(2a)')' Wrong value of mode_paral ',mode
4363   MSG_BUG(msg)
4364  end if
4365  !
4366  ! === Print out matrix ===
4367  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
4368 
4369  write(fmth,*)'(6x,',mr,'(i2,6x))'
4370  write(fmt1,*)'(3x,',mr,'f8.3)'
4371 
4372  write(msg,fmth)(ii,ii=1,mr)
4373  call wrtout(unt,msg,mode) !header
4374  write(msg,fmt1)REAL (arr(1:mr))
4375  call wrtout(unt,msg,mode) !real part
4376  write(msg,fmt1)AIMAG(arr(1:mr))
4377  call wrtout(unt,msg,mode) !imag part
4378 
4379 end subroutine print_arr1d_spc

m_numeric_tools/print_arr2d_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr2d_dpc

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

4543 subroutine print_arr2d_dpc(arr,max_r,max_c,unit,mode_paral)
4544 
4545 
4546 !This section has been created automatically by the script Abilint (TD).
4547 !Do not modify the following lines by hand.
4548 #undef ABI_FUNC
4549 #define ABI_FUNC 'print_arr2d_dpc'
4550 !End of the abilint section
4551 
4552  implicit none
4553 
4554 !Arguments ------------------------------------
4555 !scalars
4556  integer,optional,intent(in) :: unit,max_r,max_c
4557  character(len=4),optional,intent(in) :: mode_paral
4558 !arrays
4559  complex(dpc),intent(in) :: arr(:,:)
4560 
4561 !Local variables-------------------------------
4562 !scalars
4563  integer :: unt,ii,jj,nc,nr,mc,mr
4564  character(len=4) :: mode
4565  character(len=500) :: msg
4566  character(len=100) :: fmth,fmt1,fmt2
4567 ! *************************************************************************
4568 
4569  unt =std_out; if (PRESENT(unit      )) unt =unit
4570  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
4571  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
4572  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
4573 
4574  if (mode/='COLL'.and.mode/='PERS') then
4575    write(msg,'(2a)')'Wrong value of mode_paral ',mode
4576    MSG_BUG(msg)
4577  end if
4578  !
4579  ! === Print out matrix ===
4580  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
4581  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
4582 
4583  write(fmth,*)'(6x,',mc,'(i2,6x))'
4584  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
4585  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
4586 
4587  write(msg,fmth)(jj,jj=1,mc)
4588  call wrtout(unt,msg,mode) !header
4589  do ii=1,mr
4590    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
4591    call wrtout(unt,msg,mode) !real part
4592    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
4593    call wrtout(unt,msg,mode) !imag part
4594  end do
4595 
4596 end subroutine print_arr2d_dpc

m_numeric_tools/print_arr2d_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr2d_spc

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

4469 subroutine print_arr2d_spc(arr,max_r,max_c,unit,mode_paral)
4470 
4471 
4472 !This section has been created automatically by the script Abilint (TD).
4473 !Do not modify the following lines by hand.
4474 #undef ABI_FUNC
4475 #define ABI_FUNC 'print_arr2d_spc'
4476 !End of the abilint section
4477 
4478  implicit none
4479 
4480 !Arguments ------------------------------------
4481 !scalars
4482  integer,optional,intent(in) :: unit,max_r,max_c
4483  character(len=4),optional,intent(in) :: mode_paral
4484 !arrays
4485  complex(spc),intent(in) :: arr(:,:)
4486 
4487 !Local variables-------------------------------
4488 !scalars
4489  integer :: unt,ii,jj,nc,nr,mc,mr
4490  character(len=4) :: mode
4491  character(len=500) :: msg
4492  character(len=100) :: fmth,fmt1,fmt2
4493 ! *************************************************************************
4494 
4495  unt =std_out; if (PRESENT(unit      )) unt =unit
4496  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
4497  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
4498  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
4499 
4500  if (mode/='COLL'.and.mode/='PERS') then
4501    write(msg,'(2a)')'Wrong value of mode_paral ',mode
4502    MSG_BUG(msg)
4503  end if
4504  !
4505  ! === Print out matrix ===
4506  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
4507  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
4508 
4509  write(fmth,*)'(6x,',mc,'(i2,6x))'
4510  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
4511  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
4512 
4513  write(msg,fmth)(jj,jj=1,mc)
4514  call wrtout(unt,msg,mode) !header
4515  do ii=1,mr
4516    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
4517    call wrtout(unt,msg,mode) !real part
4518    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
4519    call wrtout(unt,msg,mode) !imag part
4520  end do
4521 
4522 end subroutine print_arr2d_spc

m_numeric_tools/quadrature [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  quadrature

FUNCTION

  Driver routine to perform quadratures in finite domains using different techniques.
  The routine improves the resolution of the grid until a given accuracy is reached

INPUTS

  func(external)=the name of the function to be integrated
  xmin,xmax=the limits of integration
  npts=Initial number of points, only for Gauss-Legendre. At each step this number is doubled
  accuracy=fractional accuracy required
  ntrial=Max number of attempts
  qopt=integer flag defining the algorithm for the quadrature:
    1 for Trapezoidal rule, closed, O(1/N^2)
    2 for Simpson based on trapezoidal,closed, O(1/N^4)
    3 for Midpoint rule, open, O(1/N^2)
    4 for midpoint rule with cancellation of leading error, open, O(1/N^4)
    5 for Romberg integration (closed form) and extrapolation for h-->0 (order 10 is hard-coded)
    6 for Romberg integration with midpoint rule and extrapolation for h-->0 (order 10 is hard-coded)
    7 for Gauss-Legendre

OUTPUT

  quad=the integral
  ierr=0 if quadrature converged.

PARENTS

CHILDREN

SOURCE

2996 recursive subroutine quadrature(func,xmin,xmax,qopt,quad,ierr,ntrial,accuracy,npts)
2997 
2998 
2999 !This section has been created automatically by the script Abilint (TD).
3000 !Do not modify the following lines by hand.
3001 #undef ABI_FUNC
3002 #define ABI_FUNC 'quadrature'
3003 !End of the abilint section
3004 
3005  implicit none
3006 
3007 !Arguments ------------------------------------
3008 !scalars
3009  integer,intent(in) :: qopt
3010  integer,intent(out) :: ierr
3011  integer,optional,intent(in) :: ntrial,npts
3012  real(dp),intent(in) :: xmin,xmax
3013  real(dp),optional,intent(in) :: accuracy
3014  real(dp),intent(out) :: quad
3015 
3016  interface
3017    function func(x)
3018      use defs_basis
3019      real(dp),intent(in) :: x
3020      real(dp) :: func
3021    end function func
3022  end interface
3023 
3024  !interface
3025  ! function func(x)
3026  !  use defs_basis
3027  !  real(dp),intent(in) :: x(:)
3028  !  real(dp) :: func(SIZE(x))
3029  ! end function func
3030  !end interface
3031 
3032 !Local variables-------------------------------
3033 !scalars
3034  integer :: K,KM,NT,NX,NX0,it,ix
3035  real(dp) :: EPS,old_st,st,old_quad,dqromb
3036  real(dp) :: TOL
3037  character(len=500) :: msg
3038 !arrays
3039  real(dp),allocatable :: h(:),s(:)
3040  real(dp),allocatable :: wx(:),xx(:)
3041 ! *************************************************************************
3042 
3043  ierr = 0
3044  TOL  =tol12
3045  EPS  =tol6  ; if (PRESENT(accuracy)) EPS=accuracy
3046  NT   =20    ; if (PRESENT(ntrial  )) NT=ntrial
3047  quad =zero
3048 
3049  select case (qopt)
3050 
3051  case (1)
3052    ! === Trapezoidal, closed form, O(1/N^2)
3053    do it=1,NT
3054      call trapezoidal_(func,it,xmin,xmax,quad)
3055      if (it>5) then ! Avoid spurious early convergence
3056        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
3057      end if
3058      old_quad=quad
3059    end do
3060 
3061  case (2)
3062   ! === Extended Simpson rule based on trapezoidal O(1/N^4) ===
3063    do it=1,NT
3064      call trapezoidal_(func,it,xmin,xmax,st)
3065      if (it==1) then
3066        quad=st
3067      else
3068        quad=(four*st-old_st)/three
3069      end if
3070      if (it>5) then ! Avoid spurious early convergence
3071        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
3072      end if
3073      old_quad=quad
3074      old_st=st
3075    end do
3076 
3077  case (3)
3078   ! === Midpoint rule, open form, O(1/N^2) ===
3079   do it=1,NT
3080     call midpoint_(func,it,xmin,xmax,quad)
3081     if (it>4) then ! Avoid spurious early convergence
3082       if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
3083     end if
3084     old_quad=quad
3085   end do
3086 
3087  case (4)
3088    ! === Midpoint rule with cancellation of leading 1/N^2 term, open form, O(1/N^4) ===
3089    do it=1,NT
3090      call midpoint_(func,it,xmin,xmax,st)
3091      if (it==1) then
3092        quad=st
3093      else
3094        quad=(nine*st-old_st)/eight
3095      end if
3096      if (it>4) then ! Avoid spurious early convergence
3097        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
3098      end if
3099      old_quad=quad
3100      old_st=st
3101    end do
3102 
3103  case (5)
3104    ! === Romberg Integration, closed form ===
3105    K=5 ; KM=K-1 ! Order 10
3106    ABI_MALLOC(h,(NT+1))
3107    ABI_MALLOC(s,(NT+1))
3108    h=zero
3109    s=zero
3110    h(1)=one
3111    do it=1,NT
3112      call trapezoidal_(func,it,xmin,xmax,s(it))
3113      !write(std_out,*) ' romberg-trap at ',ncall,it,s(it)
3114      if (it>=K) then
3115        call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
3116        if (ABS(dqromb)<EPS*ABS(quad)) then
3117          ABI_FREE(h)
3118          ABI_FREE(s)
3119          RETURN
3120        end if
3121      end if
3122      s(it+1)=s(it)
3123      h(it+1)=quarter*h(it) ! Quarter makes the extrapolation a polynomial in h^2,
3124    end do                 ! This is required to use the Euler-Maclaurin formula
3125    ABI_FREE(h)
3126    ABI_FREE(s)
3127 
3128  case (6)
3129    ! === Romberg Integration, closed form ===
3130    K=5 ; KM=K-1 ! Order 10
3131    ABI_MALLOC(h,(NT+1))
3132    ABI_MALLOC(s,(NT+1))
3133    h=zero
3134    s=zero
3135    h(1)=one
3136    do it=1,NT
3137      call midpoint_(func,it,xmin,xmax,s(it))
3138      if (it>=K) then
3139        call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
3140        !write(std_out,*) quad,dqromb
3141        if (ABS(dqromb)<EPS*ABS(quad)) then
3142          ABI_FREE(h)
3143          ABI_FREE(s)
3144          RETURN
3145        end if
3146      end if
3147      s(it+1)=s(it)
3148      h(it+1)=ninth*h(it) ! factor is due to step tripling in midpoint and even error series
3149    end do
3150    ABI_FREE(h)
3151    ABI_FREE(s)
3152 
3153  case (7)
3154    ! === Gauss-Legendre ===
3155    NX0=5 ; if (PRESENT(npts)) NX0=npts
3156    NX=NX0
3157    do it=1,NT
3158      ABI_MALLOC(wx,(NX))
3159      ABI_MALLOC(xx,(NX))
3160      call coeffs_gausslegint(xmin,xmax,xx,wx,NX)
3161      quad=zero
3162      do ix=1,NX
3163        quad=quad+wx(ix)*func(xx(ix))
3164      end do
3165      ABI_FREE(wx)
3166      ABI_FREE(xx)
3167      if (it>1) then
3168        !write(std_out,*) quad
3169        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
3170      end if
3171      old_quad=quad
3172      NX=NX+NX0
3173      !NX=2*NX
3174    end do
3175 
3176  case default
3177    write(msg,'(a,i3)')'Wrong value for qopt',qopt
3178    MSG_BUG(msg)
3179  end select
3180 
3181  write(msg,'(a,i0,2(a,es14.6))')&
3182 &  "Results are not converged within the given accuracy. ntrial= ",NT,"; EPS= ",EPS,"; TOL= ",TOL
3183  MSG_WARNING(msg)
3184  ierr = -1
3185 
3186 end subroutine quadrature

m_numeric_tools/rdp2cdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_2D

FUNCTION

INPUTS

OUTPUT

SOURCE

1197 pure function rdp2cdp_2D(rr) result(cc)
1198 
1199 
1200 !This section has been created automatically by the script Abilint (TD).
1201 !Do not modify the following lines by hand.
1202 #undef ABI_FUNC
1203 #define ABI_FUNC 'rdp2cdp_2D'
1204 !End of the abilint section
1205 
1206  implicit none
1207 
1208 !Arguments ------------------------------------
1209 !scalars
1210  real(dp),intent(in) :: rr(:,:,:)
1211  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3))
1212 
1213 ! *********************************************************************
1214 
1215  cc(:,:)=CMPLX(rr(1,:,:),rr(2,:,:), kind=dpc)
1216 
1217 end function rdp2cdp_2D

m_numeric_tools/rdp2cdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_3D

FUNCTION

INPUTS

OUTPUT

SOURCE

1234 pure function rdp2cdp_3D(rr) result(cc)
1235 
1236 
1237 !This section has been created automatically by the script Abilint (TD).
1238 !Do not modify the following lines by hand.
1239 #undef ABI_FUNC
1240 #define ABI_FUNC 'rdp2cdp_3D'
1241 !End of the abilint section
1242 
1243  implicit none
1244 
1245 !Arguments ------------------------------------
1246 !scalars
1247  real(dp),intent(in) :: rr(:,:,:,:)
1248  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4))
1249 
1250 ! *********************************************************************
1251 
1252  cc(:,:,:)=CMPLX(rr(1,:,:,:),rr(2,:,:,:), kind=dpc)
1253 
1254 end function rdp2cdp_3D

m_numeric_tools/rdp2cdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_4D

FUNCTION

INPUTS

OUTPUT

SOURCE

1271 pure function rdp2cdp_4D(rr) result(cc)
1272 
1273 
1274 !This section has been created automatically by the script Abilint (TD).
1275 !Do not modify the following lines by hand.
1276 #undef ABI_FUNC
1277 #define ABI_FUNC 'rdp2cdp_4D'
1278 !End of the abilint section
1279 
1280  implicit none
1281 
1282 !Arguments ------------------------------------
1283 !scalars
1284  real(dp),intent(in) :: rr(:,:,:,:,:)
1285  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5))
1286 
1287 ! *********************************************************************
1288 
1289  cc(:,:,:,:)=CMPLX(rr(1,:,:,:,:),rr(2,:,:,:,:), kind=dpc)
1290 
1291 end function rdp2cdp_4D

m_numeric_tools/rdp2cdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_5D

FUNCTION

INPUTS

OUTPUT

SOURCE

1308 pure function rdp2cdp_5D(rr) result(cc)
1309 
1310 
1311 !This section has been created automatically by the script Abilint (TD).
1312 !Do not modify the following lines by hand.
1313 #undef ABI_FUNC
1314 #define ABI_FUNC 'rdp2cdp_5D'
1315 !End of the abilint section
1316 
1317  implicit none
1318 
1319 !Arguments ------------------------------------
1320 !scalars
1321  real(dp),intent(in) :: rr(:,:,:,:,:,:)
1322  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6))
1323 
1324 ! *********************************************************************
1325 
1326  cc(:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:),rr(2,:,:,:,:,:), kind=dpc)
1327 
1328 end function rdp2cdp_5D

m_numeric_tools/rdp2cdp_6D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_6D

FUNCTION

INPUTS

OUTPUT

SOURCE

1345 pure function rdp2cdp_6D(rr) result(cc)
1346 
1347 
1348 !This section has been created automatically by the script Abilint (TD).
1349 !Do not modify the following lines by hand.
1350 #undef ABI_FUNC
1351 #define ABI_FUNC 'rdp2cdp_6D'
1352 !End of the abilint section
1353 
1354  implicit none
1355 
1356 !Arguments ------------------------------------
1357 !scalars
1358  real(dp),intent(in) :: rr(:,:,:,:,:,:,:)
1359  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6),SIZE(rr,7))
1360 
1361 ! *********************************************************************
1362 
1363  cc(:,:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:,:),rr(2,:,:,:,:,:,:), kind=dpc)
1364 
1365 end function rdp2cdp_6D

m_numeric_tools/remove_copies [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  remove_copies

FUNCTION

  Given an initial set of elements, set_in, return the subset of inequivalent items
  packed in the firt n_out positions of set_in. Use the logical function is_equal
  to define whether two items are equivalent.

INPUTS

 n_in=Initial number of elements.
 is_equal=logical function used to discern if two items are equal.

OUTPUT

 n_out=Number of inequivalent items found.

SIDE EFFECTS

 set_in(3,n_in)=
  In input the initial set of n_in elements
  In output set_in(3,1:n_out) contains the inequivalent elements found.

NOTES

 The routines only deals with arrays of 3D-vectors although generalizing the
 algorithm to nD-space is straightforward.

PARENTS

      m_io_screening

CHILDREN

SOURCE

4950 subroutine remove_copies(n_in,set_in,n_out,is_equal)
4951 
4952 
4953 !This section has been created automatically by the script Abilint (TD).
4954 !Do not modify the following lines by hand.
4955 #undef ABI_FUNC
4956 #define ABI_FUNC 'remove_copies'
4957 !End of the abilint section
4958 
4959  implicit none
4960 
4961 !Arguments ------------------------------------
4962 !scalars
4963  integer,intent(in) :: n_in
4964  integer,intent(out) :: n_out
4965 
4966 !arrays
4967  real(dp),target,intent(inout) :: set_in(3,n_in)
4968 
4969  interface
4970    function is_equal(k1,k2)
4971      use defs_basis
4972      real(dp),intent(in) :: k1(3),k2(3)
4973      logical :: is_equal
4974    end function is_equal
4975  end interface
4976 
4977 !Local variables-------------------------------
4978 !scalars
4979  integer :: ii,jj
4980  logical :: isnew
4981 !arrays
4982  type rdp1d_pt
4983   integer :: idx
4984   real(dp),pointer :: rpt(:)
4985  end type rdp1d_pt
4986  type(rdp1d_pt),allocatable :: Ap(:)
4987 
4988 ! *************************************************************************
4989 
4990  ABI_DATATYPE_ALLOCATE(Ap,(n_in))
4991  Ap(1)%idx = 1
4992  Ap(1)%rpt => set_in(:,1)
4993 
4994  n_out=1
4995  do ii=2,n_in
4996 
4997    isnew=.TRUE.
4998    do jj=1,n_out
4999      if (is_equal(set_in(:,ii),Ap(jj)%rpt(:))) then
5000        isnew=.FALSE.
5001        exit
5002      end if
5003    end do
5004 
5005    if (isnew) then
5006      n_out=n_out+1
5007      Ap(n_out)%rpt => set_in(:,ii)
5008      Ap(n_out)%idx = ii
5009    end if
5010  end do
5011 
5012  ! The n_out inequivalent items are packed first.
5013  if (n_out/=n_in) then
5014    do ii=1,n_out
5015      jj=Ap(ii)%idx
5016      set_in(:,ii) = set_in(:,jj)
5017      !write(std_out,*) Ap(ii)%idx,Ap(ii)%rpt(:)
5018    end do
5019  end if
5020 
5021  ABI_DATATYPE_DEALLOCATE(Ap)
5022 
5023 end subroutine remove_copies

m_numeric_tools/reverse_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  reverse_int

FUNCTION

   Reverse a 1D array *IN PLACE*. Target: INT arrays

PARENTS

CHILDREN

SOURCE

444 subroutine reverse_int(arr)
445 
446 
447 !This section has been created automatically by the script Abilint (TD).
448 !Do not modify the following lines by hand.
449 #undef ABI_FUNC
450 #define ABI_FUNC 'reverse_int'
451 !End of the abilint section
452 
453  implicit none
454 
455 !Arguments ------------------------------------
456 !scalars
457  integer,intent(inout) :: arr(:)
458 !arrays
459  integer :: ii,nn,swap
460 ! *************************************************************************
461 
462  nn = SIZE(arr)
463  if (nn <= 1) return
464 
465  do ii=1,nn/2
466    swap = arr(ii)
467    arr(ii) = arr(nn-ii+1)
468    arr(nn-ii+1) = swap
469  end do
470 
471 end subroutine reverse_int

m_numeric_tools/reverse_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  reverse_rdp

FUNCTION

   Reverse a 1D array *IN PLACE*. Target: DP arrays

PARENTS

CHILDREN

SOURCE

489 subroutine reverse_rdp(arr)
490 
491 
492 !This section has been created automatically by the script Abilint (TD).
493 !Do not modify the following lines by hand.
494 #undef ABI_FUNC
495 #define ABI_FUNC 'reverse_rdp'
496 !End of the abilint section
497 
498  implicit none
499 
500 !Arguments ------------------------------------
501 !scalars
502  real(dp),intent(inout) :: arr(:)
503 !arrays
504  integer :: ii,nn
505  real(dp) :: swap
506 ! *************************************************************************
507 
508  nn = SIZE(arr)
509  if (nn <= 1) return
510 
511  do ii=1,nn/2
512    swap = arr(ii)
513    arr(ii) = arr(nn-ii+1)
514    arr(nn-ii+1) = swap
515  end do
516 
517 end subroutine reverse_rdp

m_numeric_tools/rhophi [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 rhophi

FUNCTION

 Compute the phase and the module of a complex number.
 The phase angle is fold into the interval [-pi,pi]

INPUTS

  cx(2) = complex number

OUTPUT

  phi = phase of cx fold into [-pi,pi]
  rho = modul of cx

PARENTS

      berryphase_new,etheta,linemin

SOURCE

6099 pure subroutine rhophi(cx,phi,rho)
6100 
6101 
6102 !This section has been created automatically by the script Abilint (TD).
6103 !Do not modify the following lines by hand.
6104 #undef ABI_FUNC
6105 #define ABI_FUNC 'rhophi'
6106 !End of the abilint section
6107 
6108  implicit none
6109 
6110 !Arguments ------------------------------------
6111 !scalars
6112  real(dp),intent(out) :: phi,rho
6113 !arrays
6114  real(dp),intent(in) :: cx(2)
6115 
6116 
6117 ! ***********************************************************************
6118 
6119  rho = sqrt(cx(1)*cx(1) + cx(2)*cx(2))
6120 
6121  if (abs(cx(1)) > tol8) then
6122 
6123    phi = atan(cx(2)/cx(1))
6124 
6125 !  phi is an element of [-pi,pi]
6126    if (cx(1) < zero) then
6127      if (phi < zero) then
6128        phi = phi + pi
6129      else
6130        phi = phi - pi
6131      end if
6132    end if
6133 
6134  else
6135 
6136    if (cx(2) > tol8) then
6137      phi = pi*half
6138    else if (cx(2) < tol8) then
6139      phi = -0.5_dp*pi
6140    else
6141      phi = 0
6142    end if
6143 
6144  end if
6145 
6146 end subroutine rhophi

m_numeric_tools/simpson [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 simpson

FUNCTION

   Simpson integral of input function

INPUTS

  step = space between integral arguments
  values(npts)=integrand function.

OUTPUT

  integral of values on the full mesh.

PARENTS

CHILDREN

SOURCE

6047 function simpson(step,values) result(res)
6048 
6049 
6050 !This section has been created automatically by the script Abilint (TD).
6051 !Do not modify the following lines by hand.
6052 #undef ABI_FUNC
6053 #define ABI_FUNC 'simpson'
6054 !End of the abilint section
6055 
6056  implicit none
6057 
6058 !Arguments ------------------------------------
6059 !scalars
6060  real(dp),intent(in) :: step
6061  real(dp) :: res
6062 !arrays
6063  real(dp),intent(in) :: values(:)
6064 
6065 !Local variables -------------------------
6066 !scalars
6067  real(dp) :: int_values(size(values))
6068 
6069 ! *********************************************************************
6070 
6071  call simpson_int(size(values),step,values,int_values)
6072  res = int_values(size(values))
6073 
6074 end function simpson

m_numeric_tools/simpson_cplx [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  simpson_cplx

FUNCTION

  Integrate a complex function using extended Simpson's rule.

INPUTS

  npts=Number of points.
  step=Step of the mesh.
  ff(npts)=Values of the integrand.

OUTPUT

  simpson_cplx=Integral of ff.

NOTES

  If npts is odd, the integration is done with the extended Simpson's rule (error = O^(step^4))
  If npts is even, the last 4 four points are integrated separately via Simpson's 3/8 rule. Error = O(step^5)
  while the first npts-3 points are integrared with the extended Simpson's rule.

PARENTS

CHILDREN

SOURCE

3634 function simpson_cplx(npts,step,ff)
3635 
3636 
3637 !This section has been created automatically by the script Abilint (TD).
3638 !Do not modify the following lines by hand.
3639 #undef ABI_FUNC
3640 #define ABI_FUNC 'simpson_cplx'
3641 !End of the abilint section
3642 
3643  implicit none
3644 
3645 !Arguments ------------------------------------
3646 !scalars
3647  integer,intent(in) :: npts
3648  real(dp),intent(in)  :: step
3649  complex(dpc),intent(in) :: ff(npts)
3650  complex(dpc) :: simpson_cplx
3651 
3652 !Local variables ------------------------------
3653 !scalars
3654  integer :: ii,my_n
3655  complex(dpc) :: sum_even, sum_odd
3656 
3657 !************************************************************************
3658 
3659  my_n=npts; if ((npts/2)*2 == npts) my_n=npts-3
3660 
3661  if (my_n<2) then
3662    MSG_ERROR("Too few points")
3663  end if
3664 
3665  sum_odd=czero
3666  do ii=2,my_n-1,2
3667    sum_odd = sum_odd + ff(ii)
3668  end do
3669 
3670  sum_even=zero
3671  do ii=3,my_n-2,2
3672    sum_even = sum_even + ff(ii)
3673  end do
3674 
3675  ! Eq 25.4.6 Abramowitz. Error is O(step^4)
3676  simpson_cplx = step/three * (ff(1) + four*sum_odd + two*sum_even + ff(my_n))
3677 
3678  if (my_n/=npts) then ! Simpson's 3/8 rule. Eq 25.4.13 Abramowitz. Error is O(step^5)
3679   simpson_cplx = simpson_cplx + three*step/eight * (ff(npts-3) + 3*ff(npts-2) + 3*ff(npts-1) + ff(npts))
3680  end if
3681 
3682 end function simpson_cplx

m_numeric_tools/simpson_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 simpson_int

FUNCTION

   Simpson integral of input function

INPUTS

  npts=max number of points on grid for integral
  step = space between integral arguments
  values(npts)=integrand function.

OUTPUT

  int_values(npts)=integral of values.

PARENTS

      eliashberg_1d,evdw_wannier,kramerskronig,m_ebands,m_exc_spectra
      m_numeric_tools,m_phgamma,m_phonons,m_xc_vdw,mka2f,mka2fQgrid,mka2f_tr
      mka2f_tr_lova,mlwfovlp_projpaw,mlwfovlp_radial,outscfcv,radsintr

CHILDREN

SOURCE

5968 subroutine simpson_int(npts,step,values,int_values)
5969 
5970 
5971 !This section has been created automatically by the script Abilint (TD).
5972 !Do not modify the following lines by hand.
5973 #undef ABI_FUNC
5974 #define ABI_FUNC 'simpson_int'
5975 !End of the abilint section
5976 
5977  implicit none
5978 
5979 !Arguments ------------------------------------
5980 !scalars
5981  integer,intent(in) :: npts
5982  real(dp),intent(in) :: step
5983 !arrays
5984  real(dp),intent(in) :: values(npts)
5985  real(dp),intent(out) :: int_values(npts)
5986 
5987 !Local variables -------------------------
5988 !scalars
5989  integer :: ii
5990  real(dp),parameter :: coef1 = 0.375_dp                          !9.0_dp  / 24.0_dp
5991  real(dp),parameter :: coef2 = 1.166666666666666666666666667_dp  !28.0_dp / 24.0_dp
5992  real(dp),parameter :: coef3 = 0.958333333333333333333333333_dp  !23.0_dp / 24.0_dp
5993  character(len=500) :: msg
5994 
5995 ! *********************************************************************
5996 
5997  if (npts < 6) then
5998    write(msg,"(a,i0)")"Number of points in integrand function must be >=6 while it is: ",npts
5999    MSG_ERROR(msg)
6000  end if
6001 
6002 !-----------------------------------------------------------------
6003 !Simpson integral of input function
6004 !-----------------------------------------------------------------
6005 
6006 !first point is 0: don t store it
6007 !do integration equivalent to Simpson O(1/N^4) from NumRec in C p 134  NumRec in Fortran p 128
6008  int_values(1) =               coef1*values(1)
6009  int_values(2) = int_values(1) + coef2*values(2)
6010  int_values(3) = int_values(2) + coef3*values(3)
6011 
6012  do ii=4,npts-3
6013    int_values(ii) = int_values(ii-1) + values(ii)
6014  end do
6015 
6016  int_values(npts-2) = int_values(npts-3) + coef3*values(npts-2)
6017  int_values(npts-1) = int_values(npts-2) + coef2*values(npts-1)
6018  int_values(npts  ) = int_values(npts-1) + coef1*values(npts  )
6019 
6020  int_values(:) = int_values(:) * step
6021 
6022 end subroutine simpson_int

m_numeric_tools/smooth [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  smooth data.

FUNCTION

  smooth

INPUTS

  mesh=Number of points.
  it=Number of iterations.

SIDE EFFECTS

  a(mesh)=Input values, smoothed in output

PARENTS

      psp6cc,upf2abinit

CHILDREN

SOURCE

6298 subroutine smooth(a,mesh,it)
6299 
6300 
6301 !This section has been created automatically by the script Abilint (TD).
6302 !Do not modify the following lines by hand.
6303 #undef ABI_FUNC
6304 #define ABI_FUNC 'smooth'
6305 !End of the abilint section
6306 
6307  implicit none
6308 
6309 !Arguments ------------------------------------
6310 !scalars
6311  integer, intent(in) :: it,mesh
6312  real(dp), intent(inout) :: a(mesh)
6313 !Local variables-------------------------------
6314  integer :: i,k
6315  real(dp) :: asm(mesh)
6316 ! *********************************************************************
6317 
6318  do k=1,it
6319    asm(1)=1.0d0/3.0d0*(a(1)+a(2)+a(3))
6320    asm(2)=0.25d0*(a(1)+a(2)+a(3)+a(4))
6321    asm(3)=0.2d0*(a(1)+a(2)+a(3)+a(4)+a(5))
6322    asm(4)=0.2d0*(a(2)+a(3)+a(4)+a(5)+a(6))
6323    asm(5)=0.2d0*(a(3)+a(4)+a(5)+a(6)+a(7))
6324    asm(mesh-4)=0.2d0*(a(mesh-2)+a(mesh-3)+a(mesh-4)+&
6325 &                   a(mesh-5)+a(mesh-6))
6326    asm(mesh-3)=0.2d0*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
6327 &                   a(mesh-4)+a(mesh-5))
6328    asm(mesh-2)=0.2d0*(a(mesh)+a(mesh-1)+a(mesh-2)+&
6329 &                   a(mesh-3)+a(mesh-4))
6330    asm(mesh-1)=0.25d0*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
6331    asm(mesh)=1.0d0/3.0d0*(a(mesh)+a(mesh-1)+a(mesh-2))
6332 
6333    do i=6,mesh-5
6334      asm(i)=0.1d0*a(i)+0.1d0*(a(i+1)+a(i-1))+&
6335 &             0.1d0*(a(i+2)+a(i-2))+&
6336 &             0.1d0*(a(i+3)+a(i-3))+&
6337 &             0.1d0*(a(i+4)+a(i-4))+&
6338 &             0.05d0*(a(i+5)+a(i-5))
6339    end do
6340 
6341    do i=1,mesh
6342      a(i)=asm(i)
6343    end do
6344  end do
6345 
6346 end subroutine smooth

m_numeric_tools/stats_eval [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  stats_eval

FUNCTION

  Helper function used to calculate the statistical parameters of a data set.

 INPUT
  arr(:)=Array with the values.

OUTPUT

  stats<stats_t>=Data type storing the parameters of the data set.

PARENTS

      m_shirley

CHILDREN

SOURCE

5527 pure function stats_eval(arr) result(stats)
5528 
5529 
5530 !This section has been created automatically by the script Abilint (TD).
5531 !Do not modify the following lines by hand.
5532 #undef ABI_FUNC
5533 #define ABI_FUNC 'stats_eval'
5534 !End of the abilint section
5535 
5536  implicit none
5537 
5538 !Arguments ------------------------------------
5539 !scalars
5540  type(stats_t) :: stats
5541 !arrays
5542  real(dp),intent(in) :: arr(:)
5543 
5544 !Local variables ------------------------------
5545 !scalars
5546  integer :: ii,nn
5547  real(dp) :: xx,x2_sum
5548 
5549 ! *************************************************************************
5550 
5551  !@stats_t
5552  stats%min   = +HUGE(one)
5553  stats%max   = -HUGE(one)
5554  stats%mean  = zero
5555 
5556  nn = SIZE(arr)
5557 
5558  do ii=1,nn
5559    xx = arr(ii)
5560    stats%max  = MAX(stats%max, xx)
5561    stats%min  = MIN(stats%min, xx)
5562    stats%mean = stats%mean + xx
5563  end do
5564 
5565  stats%mean  = stats%mean/nn
5566 
5567  ! Two-pass algorithm for the variance (more stable than the single-pass one).
5568  x2_sum = zero
5569  do ii=1,nn
5570    xx     = arr(ii)
5571    x2_sum = x2_sum + (xx - stats%mean)*(xx - stats%mean)
5572  end do
5573 
5574  if (nn>1) then
5575    stats%stdev  = x2_sum/(nn-1)
5576    stats%stdev = SQRT(ABS(stats%stdev))
5577  else
5578    stats%stdev = zero
5579  end if
5580 
5581 end function stats_eval

m_numeric_tools/stats_t [ Types ]

[ Top ] [ m_numeric_tools ] [ Types ]

NAME

 stats_t

FUNCTION

  Statistical parameters of a data distribution.

SOURCE

233  type,public :: stats_t
234    real(dp) :: mean
235    real(dp) :: stdev
236    real(dp) :: min
237    real(dp) :: max
238  end type stats_t
239 
240  public :: stats_eval  ! Calculate the statistical parameters of a data distribution.

m_numeric_tools/symmetrize_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  symmetrize_dpc

FUNCTION

  Force a square matrix to be symmetric.

INPUTS

  uplo=String describing which part of the matrix has been calculated.
    Only the first character is tested (no case sensitive). Possible values are:
    "All"= Full matrix is supplied in input
    "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
    "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.

OUTPUT

  (see side effects)

SIDE EFFECTS

  mat(:,:)=complex input matrix, symmetrized in output

PARENTS

CHILDREN

SOURCE

4180 subroutine symmetrize_dpc(mat,uplo)
4181 
4182 
4183 !This section has been created automatically by the script Abilint (TD).
4184 !Do not modify the following lines by hand.
4185 #undef ABI_FUNC
4186 #define ABI_FUNC 'symmetrize_dpc'
4187 !End of the abilint section
4188 
4189  implicit none
4190 
4191 !Arguments ------------------------------------
4192 !scalars
4193  character(len=*),intent(in) :: uplo
4194 !arrays
4195  complex(dpc),intent(inout) :: mat(:,:)
4196 
4197 !Local variables-------------------------------
4198 !scalars
4199  integer :: nn,ii,jj
4200 !arrays
4201  complex(dpc),allocatable :: tmp(:)
4202 ! *************************************************************************
4203 
4204  nn=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
4205 & __FILE__,__LINE__)
4206 
4207  select case (uplo(1:1))
4208 
4209  case ("A","a") ! Full matrix has been calculated.
4210    ABI_MALLOC(tmp,(nn))
4211    do ii=1,nn
4212      do jj=ii,nn
4213        tmp(jj)=half*(mat(ii,jj)+mat(jj,ii))
4214      end do
4215      mat(ii,ii:nn)=tmp(ii:nn)
4216      mat(ii:nn,ii)=tmp(ii:nn)
4217    end do
4218    ABI_FREE(tmp)
4219 
4220  case ("U","u") ! Only the upper triangle is used.
4221    do jj=1,nn
4222      do ii=1,jj-1
4223        mat(jj,ii) = mat(ii,jj)
4224      end do
4225    end do
4226 
4227  case ("L","l") ! Only the lower triangle is used.
4228   do jj=1,nn
4229     do ii=1,jj-1
4230       mat(ii,jj) = mat(jj,ii)
4231     end do
4232   end do
4233 
4234  case default
4235    MSG_ERROR("Wrong uplo"//TRIM(uplo))
4236  end select
4237 
4238 end subroutine symmetrize_dpc

m_numeric_tools/symmetrize_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  symmetrize_spc

FUNCTION

  Force a square matrix to be symmetric.

INPUTS

  uplo=String describing which part of the matrix has been calculated.
    Only the first character is tested (no case sensitive). Possible values are:
    "All"= Full matrix is supplied in input
    "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
    "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.

OUTPUT

  (see side effects)

SIDE EFFECTS

  mat(:,:)=complex input matrix, symmetrized at output

PARENTS

CHILDREN

SOURCE

4091 subroutine symmetrize_spc(mat,uplo)
4092 
4093 
4094 !This section has been created automatically by the script Abilint (TD).
4095 !Do not modify the following lines by hand.
4096 #undef ABI_FUNC
4097 #define ABI_FUNC 'symmetrize_spc'
4098 !End of the abilint section
4099 
4100  implicit none
4101 
4102 !Arguments ------------------------------------
4103 !scalars
4104  character(len=*),intent(in) :: uplo
4105 !arrays
4106  complex(spc),intent(inout) :: mat(:,:)
4107 
4108 !Local variables-------------------------------
4109 !scalars
4110  integer :: nn,ii,jj
4111 !arrays
4112  complex(spc),allocatable :: tmp(:)
4113 ! *************************************************************************
4114 
4115  nn=assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',&
4116 & __FILE__,__LINE__)
4117 
4118  select case (uplo(1:1))
4119 
4120  case ("A","a") ! Full matrix has been calculated.
4121    ABI_MALLOC(tmp,(nn))
4122    do ii=1,nn
4123      do jj=ii,nn
4124        tmp(jj)=half*(mat(ii,jj)+mat(jj,ii))
4125      end do
4126      mat(ii,ii:nn)=tmp(ii:nn)
4127      mat(ii:nn,ii)=tmp(ii:nn)
4128    end do
4129    ABI_FREE(tmp)
4130 
4131  case ("U","u") ! Only the upper triangle is used.
4132    do jj=1,nn
4133      do ii=1,jj-1
4134        mat(jj,ii) = mat(ii,jj)
4135      end do
4136    end do
4137 
4138  case ("L","l") ! Only the lower triangle is used.
4139   do jj=1,nn
4140     do ii=1,jj-1
4141       mat(ii,jj) = mat(jj,ii)
4142     end do
4143   end do
4144 
4145  case default
4146    MSG_ERROR("Wrong uplo"//TRIM(uplo))
4147  end select
4148 
4149 end subroutine symmetrize_spc

m_numeric_tools/trapezoidal_ [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  trapezoidal_ (PRIVATE)

FUNCTION

  Compute the n-th stage of refinement of an extended trapezoidal rule
  adding 2^(n-2) additional interior point in the finite range of integration

INPUTS

  func(external)=the name of the function to be integrated
  xmin,xmax=the limits of integration
  nn=integer defining the refinement of the mesh, each call adds 2^(n-2) additional interior points

OUTPUT

  See SIDE EFFECTS

SIDE EFFECTS

  quad=the integral at the n-th stage.

PARENTS

CHILDREN

NOTES

  When called with nn=1, the routine returns the crudest estimate of the integral
  Subsequent calls with nn=2,3,... (in that sequential order) will improve the accuracy
  by adding 2^(n-2) additional interior points. Note that quad should not be modified between sequential calls.
  Subroutine is defined as recursive to allow multi-dimensional integrations

SOURCE

2775 recursive subroutine trapezoidal_(func,nn,xmin,xmax,quad)
2776 
2777 
2778 !This section has been created automatically by the script Abilint (TD).
2779 !Do not modify the following lines by hand.
2780 #undef ABI_FUNC
2781 #define ABI_FUNC 'trapezoidal_'
2782 !End of the abilint section
2783 
2784  implicit none
2785 
2786 !Arguments ------------------------------------
2787 !scalars
2788  integer,intent(in) :: nn
2789  !real(dp),external :: func
2790  real(dp),intent(in) :: xmin,xmax
2791  real(dp),intent(inout) :: quad
2792 
2793  interface
2794    function func(x)
2795      use defs_basis
2796      real(dp),intent(in) :: x
2797      real(dp) :: func
2798    end function func
2799  end interface
2800 
2801  !interface
2802  ! function func(x)
2803  !  use defs_basis
2804  !  real(dp),intent(in) :: x(:)
2805  !  real(dp) :: func(SIZE(x))
2806  ! end function func
2807  !end interface
2808 
2809 !Local variables-------------------------------
2810 !scalars
2811  integer :: npt,ix
2812  real(dp) :: space,new,yy
2813  character(len=500) :: msg
2814 !arrays
2815  !real(dp),allocatable :: xx(:)
2816 !************************************************************************
2817 
2818  select case (nn)
2819 
2820  case (1)
2821    ! === Initial crude estimate (xmax-xmin)(f1+f2)/2 ===
2822    !quad=half*(xmax-xmin)*SUM(func((/xmin,xmax/)))
2823    quad=half*(xmax-xmin)*(func(xmin)+func(xmax))
2824 
2825  case (2:)
2826    ! === Add npt interior points of spacing space ===
2827    npt=2**(nn-2) ; space=(xmax-xmin)/npt
2828    ! === The new sum is combined with the old integral to give a refined integral ===
2829    !new=SUM(func(arth(xmin+half*space,space,npt))) !PARALLEL version
2830    !allocate(xx(npt))
2831    !xx(:)=arth(xmin+half*space,space,npt)
2832    !xx(1)=xmin+half*space
2833    !do ii=2,nn
2834    ! xx(ii)=xx(ii-1)+space
2835    !end do
2836    new=zero
2837    yy=xmin+half*space
2838    do ix=1,npt
2839     !new=new+func(xx(ix))
2840     new=new+func(yy)
2841     yy=yy+space
2842    end do
2843    !deallocate(xx)
2844    quad=half*(quad+space*new)
2845    !write(std_out,*) 'trapezoidal',quad
2846 
2847  case (:0)
2848    write(msg,'(a,i3)')'Wrong value for nn ',nn
2849    MSG_BUG(msg)
2850  end select
2851 
2852 end subroutine trapezoidal_

m_numeric_tools/uniformrandom [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  uniformrandom

FUNCTION

 Returns a uniform random deviate between 0.0 and 1.0.
 Set seed to any value < 0 to initialize or reinitialize sequence.
 Parameters are chosen from integer overflow=2**23 (conservative).
 For some documentation, see Numerical Recipes, 1986, p196.

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

6576 function uniformrandom(seed)
6577 
6578 
6579 !This section has been created automatically by the script Abilint (TD).
6580 !Do not modify the following lines by hand.
6581 #undef ABI_FUNC
6582 #define ABI_FUNC 'uniformrandom'
6583 !End of the abilint section
6584 
6585  implicit none
6586 
6587 !Arguments ------------------------------------
6588 !scalars
6589  real(dp) :: uniformrandom
6590  integer,intent(inout) :: seed
6591 
6592 !Local variables ---------------------------------------
6593  integer, parameter :: im1=11979,ia1= 430,ic1=2531
6594  integer, parameter :: im2= 6655,ia2= 936,ic2=1399
6595  integer, parameter :: im3= 6075,ia3=1366,ic3=1283
6596  integer, save :: init=0
6597  integer, save :: ii1,ii2,ii3
6598  integer :: kk
6599  real(dp) :: im1inv,im2inv
6600  real(dp), save :: table(97)
6601  character(len=500) :: message
6602 
6603 ! *********************************************************************
6604 
6605  im1inv=1.0d0/im1 ; im2inv=1.0d0/im2
6606 
6607 !Initialize on first call or when seed<0:
6608  if (seed<0.or.init==0) then
6609    seed=-abs(seed)
6610 
6611 !  First generator
6612    ii1=mod(ic1-seed,im1)
6613    ii1=mod(ia1*ii1+ic1,im1)
6614 !  Second generator
6615    ii2=mod(ii1,im2)
6616    ii1=mod(ia1*ii1+ic1,im1)
6617 !  Third generator
6618    ii3=mod(ii1,im3)
6619 
6620 !  Fill table
6621    do kk=1,97
6622      ii1=mod(ia1*ii1+ic1,im1)
6623      ii2=mod(ia2*ii2+ic2,im2)
6624      table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
6625    enddo
6626 
6627    init=1 ; seed=1
6628  end if
6629 
6630 !Third generator gives index
6631  ii3=mod(ia3*ii3+ic3,im3)
6632  kk=1+(97*ii3)/im3
6633  if (kk<1.or.kk>97) then
6634    write(message,'(a,2i0,a)' ) ' trouble in uniformrandom; ii3,kk=',ii3,kk,' =>stop'
6635    MSG_ERROR(message)
6636  end if
6637  uniformrandom=table(kk)
6638 
6639 !Replace old value, based on generators 1 and 2
6640  ii1=mod(ia1*ii1+ic1,im1)
6641  ii2=mod(ia2*ii2+ic2,im2)
6642  table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
6643 
6644 end function uniformrandom

m_numeric_tools/unit_matrix_cdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  unit_matrix_cdp

FUNCTION

  Set the matrix matrix to be a unit matrix (if it is square).

SIDE EFFECTS

  matrix(:,:)=set to unit on exit

PARENTS

CHILDREN

SOURCE

630 pure subroutine unit_matrix_cdp(matrix)
631 
632 
633 !This section has been created automatically by the script Abilint (TD).
634 !Do not modify the following lines by hand.
635 #undef ABI_FUNC
636 #define ABI_FUNC 'unit_matrix_cdp'
637 !End of the abilint section
638 
639  implicit none
640 
641 !Arguments ------------------------------------
642  complex(dpc),intent(inout) :: matrix(:,:)
643 
644 !Local variables-------------------------------
645 !scalars
646  integer :: ii,nn
647 ! *********************************************************************
648 
649  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
650  matrix=czero
651  do ii=1,nn
652    matrix(ii,ii)=cone
653  end do
654 
655 end subroutine unit_matrix_cdp

m_numeric_tools/unit_matrix_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  unit_matrix_int

FUNCTION

  Set the matrix matrix to be a unit matrix (if it is square).

SIDE EFFECTS

  matrix(:,:)=set to unit on exit

PARENTS

CHILDREN

SOURCE

538 pure subroutine unit_matrix_int(matrix)
539 
540 
541 !This section has been created automatically by the script Abilint (TD).
542 !Do not modify the following lines by hand.
543 #undef ABI_FUNC
544 #define ABI_FUNC 'unit_matrix_int'
545 !End of the abilint section
546 
547  implicit none
548 
549 !Arguments ------------------------------------
550  integer,intent(inout) :: matrix(:,:)
551 
552 !Local variables-------------------------------
553 !scalars
554  integer :: ii,nn
555 ! *********************************************************************
556 
557  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
558  matrix(:,:)=0
559  do ii=1,nn
560   matrix(ii,ii)=1
561  end do
562 
563 end subroutine unit_matrix_int

m_numeric_tools/unit_matrix_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  unit_matrix_rdp

FUNCTION

  Set the matrix matrix to be a unit matrix (if it is square).

SIDE EFFECTS

  matrix(:,:)=set to unit on exit

PARENTS

CHILDREN

SOURCE

584 pure subroutine unit_matrix_rdp(matrix)
585 
586 
587 !This section has been created automatically by the script Abilint (TD).
588 !Do not modify the following lines by hand.
589 #undef ABI_FUNC
590 #define ABI_FUNC 'unit_matrix_rdp'
591 !End of the abilint section
592 
593  implicit none
594 
595 !Arguments ------------------------------------
596  real(dp),intent(inout) :: matrix(:,:)
597 
598 !Local variables-------------------------------
599 !scalars
600  integer :: ii,nn
601 ! *********************************************************************
602 
603  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
604  matrix(:,:)=zero
605  do ii=1,nn
606    matrix(ii,ii)=one
607  end do
608 
609 end subroutine unit_matrix_rdp

m_numeric_tools/vdiff_eval [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 vdiff_eval

FUNCTION

  Estimate the "distance" between two functions tabulated on a homogeneous grid.
  See vdiff_t

INPUTS

  cplex=1 if f1 and f2 are real, 2 for complex.
  nr=Number of points in the mesh.
  f1(cplex,nr), f2(cplex,nr)=Vectors with values

OUTPUT

  vdiff_t object

PARENTS

SOURCE

6171 pure function vdiff_eval(cplex,nr,f1,f2,volume) result(vd)
6172 
6173 
6174 !This section has been created automatically by the script Abilint (TD).
6175 !Do not modify the following lines by hand.
6176 #undef ABI_FUNC
6177 #define ABI_FUNC 'vdiff_eval'
6178 !End of the abilint section
6179 
6180  implicit none
6181 
6182 !Arguments ------------------------------------
6183 !scalars
6184  integer,intent(in) :: cplex,nr
6185  real(dp),intent(in) :: volume
6186  type(vdiff_t) :: vd
6187 !arrays
6188  real(dp),intent(in) :: f1(cplex,nr),f2(cplex,nr)
6189 
6190 !Local variables-------------------------------
6191 !scalars
6192  integer :: ir
6193  real(dp) :: num,den,dr
6194  type(stats_t) :: stats
6195 !arrays
6196  real(dp) :: abs_diff(nr)
6197 ! *********************************************************************
6198 
6199  dr = volume / nr
6200 
6201  if (cplex == 1) then
6202    abs_diff = abs(f1(1,:) - f2(1,:))
6203    num = sum(abs_diff)
6204    den = sum(abs(f2(1,:)))
6205 
6206  else if (cplex == 2) then
6207    do ir=1,nr
6208      abs_diff(ir) = sqrt((f1(1,ir) - f2(1,ir))**2 + (f1(2,ir) - f2(2,ir))**2)
6209    end do
6210    num = sum(abs_diff)
6211    den = zero
6212    do ir=1,nr
6213      den = den + sqrt(f2(1,ir)**2 + f2(2,ir)**2)
6214    end do
6215  end if
6216 
6217  vd%int_adiff = num * dr
6218  vd%l1_rerr = num / den
6219 
6220  stats = stats_eval(abs_diff)
6221  vd%mean_adiff = stats%mean
6222  vd%stdev_adiff = stats%stdev
6223  vd%min_adiff = stats%min
6224  vd%max_adiff = stats%max
6225 
6226 end function vdiff_eval

m_numeric_tools/vdiff_print [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 vdiff_print

FUNCTION

  Print vdiff_t to unit

PARENTS

      m_dvdb

CHILDREN

SOURCE

6245 subroutine vdiff_print(vd, unit)
6246 
6247 
6248 !This section has been created automatically by the script Abilint (TD).
6249 !Do not modify the following lines by hand.
6250 #undef ABI_FUNC
6251 #define ABI_FUNC 'vdiff_print'
6252 !End of the abilint section
6253 
6254  implicit none
6255 
6256 !Arguments ------------------------------------
6257 !scalars
6258  integer,optional,intent(in) :: unit
6259  type(vdiff_t),intent(in) :: vd
6260 
6261 !Local variables-------------------------------
6262 !scalars
6263  integer :: unt
6264 ! *********************************************************************
6265 
6266  unt = std_out; if (present(unit)) unt = unit
6267  write(unt,"(a,es10.3)")"  L1_rerr: ", vd%l1_rerr
6268  write(unt,"(a,es10.3)")"  Integral |f1-f2| dr: ", vd%int_adiff
6269  write(unt,"(a,es10.3)")"  min {|f1-f2|}: ", vd%min_adiff
6270  write(unt,"(a,es10.3)")"  Max {|f1-f2|}: ", vd%max_adiff
6271  write(unt,"(a,es10.3)")"  mean {|f1-f2|}: ", vd%mean_adiff
6272  write(unt,"(a,es10.3)")"  stdev {|f1-f2|}: ", vd%stdev_adiff
6273 
6274 end subroutine vdiff_print

m_numeric_tools/vdiff_t [ Types ]

[ Top ] [ m_numeric_tools ] [ Types ]

NAME

 vdiff_t

FUNCTION

  Estimate the "distance" between two functions tabulated on a homogeneous grid.
  Use `vidff` function to construct the object.

SOURCE

255  type,public :: vdiff_t
256 
257    real(dp) :: int_adiff      ! \int |f1-f2| dr
258    real(dp) :: mean_adiff     ! Mean {|f1-f2|}
259    real(dp) :: stdev_adiff    ! Standard deviation of {|f1-f2|}
260    real(dp) :: min_adiff      ! Min {|f1-f2|}
261    real(dp) :: max_adiff      ! Max {|f1-f2|}
262    real(dp) :: l1_rerr        ! (\int |f1-f2| dr) / (\int |f2| dr)
263 
264  end type vdiff_t
265 
266  public :: vdiff_eval         ! Estimate the "distance" between two functions tabulated on a homogeneous grid.
267  public :: vdiff_print        ! Print vdiff_t to formatted file.

m_numeric_tools/wrap2_pmhalf [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 wrap2_pmhalf

FUNCTION

 Transforms a real number (num) in its corresponding reduced number
 (red) in the interval ]-1/2,1/2] where -1/2 is not included (tol12)
 num=red+shift

INPUTS

 num=real number

OUTPUT

 red=reduced number of num in the interval ]-1/2,1/2] where -1/2 is not included
 shift=num-red

PARENTS

      canat9,dist2,elphon,ep_setupqpt,get_full_kgrid,interpolate_gkk
      m_bz_mesh,m_cprj_bspline,m_gamma,m_io_gkk,m_optic_tools,m_shirley
      mkfskgrid,mkfsqgrid,mkph_linwid,mkphbs,order_fs_kpts,printvtk,read_gkk
      smpbz,littlegroup_q

CHILDREN

SOURCE

5667 elemental subroutine wrap2_pmhalf(num,red,shift)
5668 
5669 
5670 !This section has been created automatically by the script Abilint (TD).
5671 !Do not modify the following lines by hand.
5672 #undef ABI_FUNC
5673 #define ABI_FUNC 'wrap2_pmhalf'
5674 !End of the abilint section
5675 
5676  implicit none
5677 
5678 !Arguments -------------------------------
5679 !scalars
5680  real(dp),intent(in) :: num
5681  real(dp),intent(out) :: red,shift
5682 
5683 ! *********************************************************************
5684 
5685  if (num>zero) then
5686    red=mod((num+half-tol12),one)-half+tol12
5687  else
5688    red=-mod(-(num-half-tol12),one)+half+tol12
5689  end if
5690  if(abs(red)<tol12)red=0.0d0
5691  shift=num-red
5692 
5693 end subroutine wrap2_pmhalf

m_numeric_tools/wrap2_zero_one [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 wrap2_zero_one

FUNCTION

 Transforms a real number (num) in its corresponding reduced number
 (red) in the interval [0,1[ where 1 is not included (tol12)
 num=red+shift

INPUTS

  num=real number

OUTPUT

 red=reduced number of num in the interval [0,1[ where 1 is not included
 shift=num-red

PARENTS

      exc_plot,k_neighbors,lin_interpq_gam,m_nesting,m_paw_pwaves_lmn
      pawmkaewf

CHILDREN

SOURCE

5610 elemental subroutine wrap2_zero_one(num,red,shift)
5611 
5612 
5613 !This section has been created automatically by the script Abilint (TD).
5614 !Do not modify the following lines by hand.
5615 #undef ABI_FUNC
5616 #define ABI_FUNC 'wrap2_zero_one'
5617 !End of the abilint section
5618 
5619  implicit none
5620 
5621 !Arguments ------------------------------------
5622 !scalars
5623  real(dp),intent(in) :: num
5624  real(dp),intent(out) :: red,shift
5625 
5626 ! *************************************************************************
5627 
5628  if (num>zero) then
5629    red=mod((num+tol12),one)-tol12
5630  else
5631    red=-mod(-(num-one+tol12),one)+one-tol12
5632  end if
5633  if(abs(red)<tol12)red=0.0_dp
5634  shift=num-red
5635 
5636 end subroutine wrap2_zero_one