TABLE OF CONTENTS

ABINIT/blocked_loop [ Functions ]

[ Top ] [ Functions ]

NAME

 blocked_loop

FUNCTION

  Helper function to implement blocked algorithms inside do loops i.e. algorithms
  operating on multiple items up to a maximum `batch_size`.

 Usage:

   batch_size = 4
   allocate(work(..., batch_size)

   do loop_index=1, loop_stop, batch_size
     ndat = blocked_loop(loop_index, loop_stop, batch_size)
     ! operate on ndat items in work
   end do

SOURCE

6425 integer pure function blocked_loop(loop_index, loop_stop, batch_size) result(ndat)
6426 
6427 !Arguments ----------------------------------------------
6428  integer,intent(in) :: loop_index, loop_stop, batch_size
6429 
6430 ! *********************************************************************
6431 
6432  ndat = merge(batch_size, loop_stop - loop_index + 1, loop_index + batch_size - 1 <= loop_stop)
6433 
6434 end function blocked_loop

ABINIT/bool2index [ Functions ]

[ Top ] [ Functions ]

NAME

 bool2index

FUNCTION

  Allocate and return array with the indices in the input boolean array `bool_list` that evaluates to .True.

SOURCE

6270 subroutine bool2index(bool_list, out_index)
6271 
6272 !Arguments ----------------------------------------------
6273 !scalars
6274  logical,intent(in) :: bool_list(:)
6275  integer,allocatable,intent(inout) :: out_index(:)
6276 
6277 !Local variables-------------------------------
6278  integer :: ii, cnt
6279 ! *********************************************************************
6280 
6281  cnt = count(bool_list)
6282  ABI_REMALLOC(out_index, (cnt))
6283  cnt = 0
6284  do ii=1,size(bool_list)
6285    if (bool_list(ii)) then
6286      cnt = cnt + 1
6287      out_index(cnt) = ii
6288    end if
6289  end do
6290 
6291 end subroutine bool2index

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

SOURCE

6120 function dotproduct(nv1,nv2,v1,v2)
6121 
6122 !Arguments ------------------------------------
6123 !scalars
6124  integer,intent(in) :: nv1,nv2
6125  real(dp) :: dotproduct
6126 !arrays
6127  real(dp),intent(in) :: v1(nv1,nv2),v2(nv1,nv2)
6128 
6129 !Local variables-------------------------------
6130 !scalars
6131  integer :: i,j
6132 
6133 ! *************************************************************************
6134  dotproduct=zero
6135  do j=1,nv2
6136   do i=1,nv1
6137    dotproduct=dotproduct+v1(i,j)*v2(i,j)
6138   end do
6139  end do
6140 
6141 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-2024 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 .

SOURCE

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

ABINIT/polynomial_regression [ Functions ]

[ Top ] [ Functions ]

NAME

  polynomial_regression

FUNCTION

  Perform a polynomial regression on incoming data points, the
  x-values of which are stored in array xvals and the y-values
  stored in array yvals. Returns a one dimensional array with
  fit coefficients (coeffs) and the unbiased RMS error of the
  fit as a scalar (RMSerr).

INPUTS

  degree = order of the polynomial
  npts = number of data points
  xvals(npts) = x-values of those data points
  yvals(npts) = y-values of those data points

OUTPUT

  coeffs(degree+1) = coefficients of the polynomial regression
  RMSerr = unbiased RMS error on the fit
            RMSerr=\sqrt{\frac{1}{npts-1}*
                      \sum_i^npts{(fitval-yvals(i))**2}}

SOURCE

6322 subroutine polynomial_regression(degree,npts,xvals,yvals,coeffs,RMSerr)
6323 
6324 !Arguments ------------------------------------
6325 
6326 !scalars
6327  integer                     :: degree,npts
6328  real(dp),intent(out)        :: RMSerr
6329 !arrays
6330  real(dp),intent(in)         :: xvals(1:npts),yvals(1:npts)
6331  real(dp),intent(out)        :: coeffs(degree+1)
6332 
6333 !Local variables-------------------------------
6334 !scalars
6335  integer                     :: ncoeffs,icoeff,ipoint,info
6336  real(dp)                    :: residual,fitval
6337 !arrays
6338  integer,allocatable         :: tmp(:)
6339  real(dp),allocatable        :: tmptwo(:)
6340  real(dp),allocatable        :: A(:,:),ATA(:,:)
6341 
6342 !####################################################################
6343 !#####################  Get Polynomial Fit  #########################
6344 
6345   ncoeffs=degree+1
6346 
6347   ABI_MALLOC(tmp,(ncoeffs))
6348   ABI_MALLOC(tmptwo,(ncoeffs))
6349   ABI_MALLOC(A,(npts,ncoeffs))
6350   ABI_MALLOC(ATA,(ncoeffs,ncoeffs))
6351 
6352   !Construct a polynomial for all input xvalues
6353   do icoeff=1,ncoeffs
6354     do ipoint=1,npts
6355        if (icoeff==1.and.xvals(ipoint)==0.0) then
6356           A(ipoint,icoeff) = 1.0
6357        else
6358           A(ipoint,icoeff) = xvals(ipoint)**(icoeff-1)
6359        end if
6360     end do
6361   end do
6362 
6363   !Get matrix product of transpose of A and A
6364   ATA = matmul(transpose(A),A)
6365 
6366   !Compute LU factorization of ATA
6367   call DGETRF(ncoeffs,ncoeffs,ATA,ncoeffs,tmp,info)
6368   ABI_CHECK(info == 0, sjoin('LAPACK DGETRF in polynomial regression returned:', itoa(info)))
6369 
6370   !Compute inverse of the LU factorized version of ATA
6371   call DGETRI(ncoeffs,ATA,ncoeffs,tmp,tmptwo,ncoeffs,info)
6372   ABI_CHECK(info == 0, sjoin('LAPACK DGETRI in polynomial regression returned:', itoa(info)))
6373 
6374   !Harvest polynomial coefficients
6375   coeffs = matmul(matmul(ATA,transpose(A)),yvals)
6376 
6377 !####################################################################
6378 !##############  RMS error on the polynomial fit  ###################
6379 
6380   residual=0.0d0
6381   do ipoint=1,npts
6382     fitval=0.0d0
6383     do icoeff=1,ncoeffs
6384       if (icoeff==1.and.xvals(ipoint)==0.0) then
6385         fitval=fitval+coeffs(icoeff)
6386       else
6387         fitval=fitval+coeffs(icoeff)*xvals(ipoint)**(icoeff-1)
6388       end if
6389     end do
6390     residual=residual+(fitval-yvals(ipoint))**2
6391   end do
6392   RMSerr=sqrt(residual/(real(npts-1,8)))
6393 
6394 !####################################################################
6395 !########################  Deallocations  ###########################
6396 
6397   ABI_FREE(A)
6398   ABI_FREE(ATA)
6399   ABI_FREE(tmp)
6400   ABI_FREE(tmptwo)
6401 
6402 end subroutine polynomial_regression

ABINIT/safe_div [ Functions ]

[ Top ] [ Functions ]

NAME

 safe_div

FUNCTION

  Subroutine safe_div performs "safe division", that is to prevent overflow,
  underflow, NaN, or infinity errors.  An alternate value is returned if the
  division cannot be performed. (bmy, 2/26/08)

  For more information, see the discussion on:
  http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/

  Taken by HM from:
  http://wiki.seas.harvard.edu/geos-chem/index.php/Floating_point_math_issues#Safe_floating-point_division

INPUTS

  n    : Numerator for the division
  d    : Divisor for the division
  altv : Alternate value to be returned if the division can't be done

OUTPUT

SOURCE

6243 elemental subroutine safe_div(n, d, altv, q)
6244 
6245 !Arguments ----------------------------------------------
6246 !scalars
6247  real(dp),intent(in) :: n, d, altv
6248  real(dp),intent(out) :: q
6249 
6250 ! *********************************************************************
6251 
6252  if ( exponent(n) - exponent(d) >= maxexponent(n) .or. d == zero) then
6253    q = altv
6254  else
6255    q = n / d
6256  endif
6257 
6258 end subroutine safe_div

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

313 pure function arth_int(start, step, nn)
314 
315 !Arguments ------------------------------------
316 !scalars
317  integer,intent(in) :: nn
318  integer,intent(in) :: start,step
319  integer :: arth_int(nn)
320 
321 !Local variables-------------------------------
322  integer :: ii
323 ! *********************************************************************
324 
325  select case (nn)
326 
327  case (1:)
328    arth_int(1)=start
329    do ii=2,nn
330     arth_int(ii)=arth_int(ii-1)+step
331    end do
332 
333  case (0)
334    return
335  end select
336 
337 end function arth_int

m_numeric_tools/arth_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  arth_rdp

FUNCTION

INPUTS

OUTPUT

SOURCE

354 pure function arth_rdp(start, step, nn)
355 
356 !Arguments ------------------------------------
357 !scalars
358  integer,intent(in) :: nn
359  real(dp),intent(in) :: start,step
360  real(dp) :: arth_rdp(nn)
361 
362 !Local variables-------------------------------
363  integer :: ii
364 ! *********************************************************************
365 
366  select case (nn)
367  case (1:)
368    arth_rdp(1)=start
369    do ii=2,nn
370     arth_rdp(ii)=arth_rdp(ii-1)+step
371    end do
372 
373  case (0)
374    return
375  end select
376 
377 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

1633 pure function bisect_int(AA,xx) result(loc)
1634 
1635 !Arguments ------------------------------------
1636 !scalars
1637  integer,intent(in) :: AA(:)
1638  integer,intent(in) :: xx
1639  integer :: loc
1640 
1641 !Local variables-------------------------------
1642  integer :: nn,jl,jm,ju
1643  logical :: ascnd
1644 ! *********************************************************************
1645 
1646  nn=SIZE(AA) ; ascnd=(AA(nn)>=AA(1))
1647 
1648  ! Initialize lower and upper limits
1649  jl=0 ; ju=nn+1
1650  do
1651   if (ju-jl<=1) EXIT
1652   jm=(ju+jl)/2  ! Compute a midpoint
1653   if (ascnd.EQV.(xx>=AA(jm))) then
1654    jl=jm ! Replace lower limit
1655   else
1656    ju=jm ! Replace upper limit
1657   end if
1658  end do
1659  !
1660  ! Set the output, being careful with the endpoints
1661  if (xx==AA(1)) then
1662   loc=1
1663  else if (xx==AA(nn)) then
1664   loc=nn-1
1665  else
1666   loc=jl
1667  end if
1668 
1669 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

1578 pure function bisect_rdp(AA, xx) result(loc)
1579 
1580 !Arguments ------------------------------------
1581 !scalars
1582  real(dp),intent(in) :: AA(:)
1583  real(dp),intent(in) :: xx
1584  integer :: loc
1585 
1586 !Local variables-------------------------------
1587  integer :: nn,jl,jm,ju
1588  logical :: ascnd
1589 ! *********************************************************************
1590 
1591  nn=SIZE(AA); ascnd=(AA(nn)>=AA(1))
1592  !
1593  ! Initialize lower and upper limits
1594  jl=0; ju=nn+1
1595  do
1596    if (ju-jl<=1) EXIT
1597    jm=(ju+jl)/2  ! Compute a midpoint,
1598    if (ascnd.EQV.(xx>=AA(jm))) then
1599      jl=jm ! Replace lower limit
1600    else
1601      ju=jm ! Replace upper limit
1602    end if
1603  end do
1604  !
1605  ! Set the output, being careful with the endpoints
1606  if (xx==AA(1)) then
1607    loc=1
1608  else if (xx==AA(nn)) then
1609    loc=nn-1
1610  else
1611    loc=jl
1612  end if
1613 
1614 end function bisect_rdp

m_numeric_tools/calculate_pade_a [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  calculate_pade_a

FUNCTION

INPUTS

OUTPUT

SOURCE

4132 subroutine calculate_pade_a(a, n, z, f)
4133 
4134 !Arguments ------------------------------------
4135 !scalars
4136  integer,intent(in) :: n
4137  complex(dpc),intent(in) :: z(n),f(n)
4138  complex(dpc),intent(out) :: a(n)
4139 
4140 !Local variables-------------------------------
4141 !scalars
4142  integer :: i,j
4143 !arrays
4144  complex(dpc) :: g(n,n)
4145 ! *************************************************************************
4146 
4147  g(1,1:n)=f(1:n)
4148 
4149  do i=2,n
4150    do j=i,n
4151      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)
4152      g(i,j)=(g(i-1,i-1)-g(i-1,j)) / ((z(j)-z(i-1))*g(i-1,j))
4153      !write(std_out,*) 'g_i(z_j)',i,j,g(i,j)
4154    end do
4155  end do
4156  do i=1,n
4157    a(i)=g(i,i)
4158  end do
4159  !write(std_out,*) 'a ',a(:)
4160 
4161 end subroutine calculate_pade_a

m_numeric_tools/cdp2rdp_0D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_0D

FUNCTION

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

INPUTS

  cc=the input complex number

OUTPUT

  rr(2=the real array

SOURCE

1181 pure function cdp2rdp_0D(cc) result(rr)
1182 
1183 !Arguments ------------------------------------
1184 !scalars
1185  complex(dpc),intent(in) :: cc
1186  real(dp) :: rr(2)
1187 
1188 ! *********************************************************************
1189 
1190  rr(1)=REAL (cc)
1191  rr(2)=AIMAG(cc)
1192 
1193 end function cdp2rdp_0D

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

SOURCE

1207 pure function cdp2rdp_1D(cc) result(rr)
1208 
1209 !Arguments ------------------------------------
1210 !scalars
1211  complex(dpc),intent(in) :: cc(:)
1212  real(dp) :: rr(2,SIZE(cc))
1213 
1214 ! *********************************************************************
1215 
1216  rr(1,:)=REAL (cc(:))
1217  rr(2,:)=AIMAG(cc(:))
1218 
1219 end function cdp2rdp_1D

m_numeric_tools/cdp2rdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_2D

FUNCTION

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

SOURCE

1233 pure function cdp2rdp_2D(cc) result(rr)
1234 
1235 !Arguments ------------------------------------
1236 !scalars
1237  complex(dpc),intent(in) :: cc(:,:)
1238  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2))
1239 ! *********************************************************************
1240 
1241  rr(1,:,:)=REAL (cc(:,:))
1242  rr(2,:,:)=AIMAG(cc(:,:))
1243 
1244 end function cdp2rdp_2D

m_numeric_tools/cdp2rdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_3D

FUNCTION

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

SOURCE

1258 pure function cdp2rdp_3D(cc) result(rr)
1259 
1260 !Arguments ------------------------------------
1261 !scalars
1262  complex(dpc),intent(in) :: cc(:,:,:)
1263  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3))
1264 
1265 ! *********************************************************************
1266 
1267  rr(1,:,:,:)=REAL (cc(:,:,:))
1268  rr(2,:,:,:)=AIMAG(cc(:,:,:))
1269 
1270 end function cdp2rdp_3D

m_numeric_tools/cdp2rdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_4D

FUNCTION

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

SOURCE

1284 pure function cdp2rdp_4D(cc) result(rr)
1285 
1286 !Arguments ------------------------------------
1287 !scalars
1288  complex(dpc),intent(in) :: cc(:,:,:,:)
1289  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4))
1290 ! *********************************************************************
1291 
1292  rr(1,:,:,:,:)=REAL (cc(:,:,:,:))
1293  rr(2,:,:,:,:)=AIMAG(cc(:,:,:,:))
1294 
1295 end function cdp2rdp_4D

m_numeric_tools/cdp2rdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_5D

FUNCTION

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

SOURCE

1309 pure function cdp2rdp_5D(cc) result(rr)
1310 
1311 !Arguments ------------------------------------
1312 !scalars
1313  complex(dpc),intent(in) :: cc(:,:,:,:,:)
1314  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5))
1315 
1316 ! *********************************************************************
1317 
1318  rr(1,:,:,:,:,:)=REAL (cc(:,:,:,:,:))
1319  rr(2,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:))
1320 
1321 end function cdp2rdp_5D

m_numeric_tools/cdp2rdp_6D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_6D

FUNCTION

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

SOURCE

1335 pure function cdp2rdp_6D(cc) result(rr)
1336 
1337 !Arguments ------------------------------------
1338 !scalars
1339  complex(dpc),intent(in) :: cc(:,:,:,:,:,:)
1340  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5),SIZE(cc,6))
1341 
1342 ! *********************************************************************
1343 
1344  rr(1,:,:,:,:,:,:)=REAL (cc(:,:,:,:,:,:))
1345  rr(2,:,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:,:))
1346 
1347 end function cdp2rdp_6D

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

SOURCE

5597 real(dp) function central_finite_diff(order, ipos, npts) result(fact)
5598 
5599 !Arguments ---------------------------------------------
5600 !scalars
5601  integer,intent(in) :: ipos,order,npts
5602 
5603 !Local variables ---------------------------------------
5604 !scalars
5605  real(dp),parameter :: empty=huge(one)
5606 ! 1st derivative.
5607  real(dp),parameter :: d1(9,4) = reshape([ &
5608   [-1/2._dp, 0._dp, 1/2._dp, empty, empty, empty, empty, empty, empty], &
5609   [ 1/12._dp, -2/3._dp, 0._dp, 2/3._dp, -1/12._dp, empty, empty, empty, empty], &
5610   [-1/60._dp, 3/20._dp, -3/4._dp, 0._dp, 3/4._dp, -3/20._dp, 1/60._dp, empty, empty], &
5611   [ 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])
5612 ! 2nd derivative.
5613  real(dp),parameter :: d2(9,4) = reshape([ &
5614    [ 1._dp, -2._dp, 1._dp, empty, empty, empty, empty, empty, empty], &
5615    [-1/12._dp, 4/3._dp, -5/2._dp, 4/3._dp, -1/12._dp, empty, empty, empty, empty], &
5616    [ 1/90._dp, -3/20._dp, 3/2._dp, -49/18._dp, 3/2._dp, -3/20._dp, 1/90._dp, empty, empty], &
5617    [-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])
5618 ! 3th derivative.
5619  real(dp),parameter :: d3(9,3) = reshape([ &
5620    [-1/2._dp, 1._dp, 0._dp, -1._dp, 1/2._dp, empty, empty, empty, empty], &
5621    [ 1/8._dp, -1._dp, 13/8._dp, 0._dp, -13/8._dp, 1._dp, -1/8._dp, empty, empty], &
5622    [ -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]], &
5623    [9,3])
5624 ! 4th derivative.
5625  real(dp),parameter :: d4(9,3) = reshape([ &
5626    [ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp, empty, empty, empty, empty], &
5627    [ -1/6._dp, 2._dp, -13/2._dp, 28/3._dp, -13/2._dp, 2._dp, -1/6._dp, empty, empty], &
5628    [ 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])
5629 ! 5th derivative.
5630  real(dp),parameter :: d5(7) = [ -1/2._dp, 2._dp, -5/2._dp, 0._dp, 5/2._dp, -2._dp, 1/2._dp]
5631 ! 6th derivative.
5632  real(dp),parameter :: d6(7) = [ 1._dp, -6._dp, 15._dp, -20._dp, 15._dp, -6._dp, 1._dp]
5633 ! *********************************************************************
5634 
5635  select case (order)
5636  case (1)
5637    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5638    fact = d1(ipos, npts/2)
5639  case (2)
5640    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5641    fact = d2(ipos, npts/2)
5642  case (3)
5643    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5644    fact = d3(ipos, npts/2)
5645  case (4)
5646    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5647    fact = d4(ipos, npts/2 - 1)
5648  case (5)
5649    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5650    fact = d5(ipos)
5651  case (6)
5652    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5653    fact = d6(ipos)
5654  case default
5655    ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
5656  end select
5657 
5658  if (fact == empty) then
5659    ABI_ERROR(sjoin("Invalid ipos:",itoa(ipos),"for order", itoa(order), "npts", itoa(npts)))
5660  end if
5661  return
5662 
5663 10 ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
5664 
5665 end function central_finite_diff

m_numeric_tools/check_vec_conjg [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 check_vec_conjg

FUNCTION

 Test whether two complex vectors `vec1` and `vec2` of size `nn` are conjugate of each other.
 within an optional tolerance abs_tol (default: 1e-6).
 Return max absolute difference for real and imag part in abs_diff(2) and exit status in ierr.

SOURCE

3742 integer function check_vec_conjg(nn, vec1, vec2, abs_diff, abs_tol) result(ierr)
3743 
3744 !Arguments ------------------------------------
3745  integer, intent(in) :: nn
3746  complex(dp), intent(in) :: vec1(nn), vec2(nn)
3747  real(dp),intent(out) :: abs_diff(2)
3748  real(dp),optional,intent(in) :: abs_tol
3749 
3750 !Local variables-------------------------------
3751  integer :: ii
3752  real(dp) :: my_abs_tol
3753 
3754  ! *************************************************************************
3755 
3756  my_abs_tol = tol6; if (present(abs_tol)) my_abs_tol = abs_tol
3757  ierr = 0
3758  abs_diff = zero
3759 
3760  do ii=1,nn
3761    abs_diff(1) = max(abs_diff(1), abs(real(vec1(ii)) - real(vec2(ii))))
3762    abs_diff(2) = max(abs_diff(2), abs(aimag(vec1(ii)) + aimag(vec2(ii))))
3763  end do
3764 
3765  if (any(abs_diff > abs_tol)) ierr = 1
3766 
3767 end function check_vec_conjg

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.

SOURCE

4580 subroutine cmplx_sphcart(carr, from, units)
4581 
4582 !Arguments ------------------------------------
4583 !scalars
4584  character(len=*),intent(in) :: from
4585  character(len=*),optional,intent(in) :: units
4586 !arrays
4587  complex(dpc),intent(inout) :: carr(:,:)
4588 
4589 !Local variables-------------------------------
4590 !scalars
4591  integer :: jj,ii
4592  real(dp) :: rho,theta,fact
4593  character(len=500) :: msg
4594 
4595 ! *************************************************************************
4596 
4597  select case (from(1:1))
4598 
4599  case ("S","s") ! Spherical --> Cartesian
4600 
4601    fact = one
4602    if (PRESENT(units)) then
4603      if (units(1:1) == "D" .or. units(1:1) == "d") fact = two_pi/360_dp
4604    end if
4605 
4606    do jj=1,SIZE(carr,DIM=2)
4607      do ii=1,SIZE(carr,DIM=1)
4608         rho  = DBLE(carr(ii,jj))
4609         theta= AIMAG(carr(ii,jj)) * fact
4610         carr(ii,jj) = CMPLX(rho*DCOS(theta), rho*DSIN(theta), kind=dpc)
4611      end do
4612    end do
4613 
4614  case ("C","c") ! Cartesian --> Spherical \theta = 2 arctan(y/(rho+x))
4615 
4616    fact = one
4617    if (PRESENT(units)) then
4618      if (units(1:1) == "D" .or. units(1:1) == "d") fact = 360_dp/two_pi
4619    end if
4620 
4621    do jj=1,SIZE(carr,DIM=2)
4622      do ii=1,SIZE(carr,DIM=1)
4623         rho  = SQRT(ABS(carr(ii,jj)))
4624         if (rho > tol16) then
4625           theta= two * ATAN( AIMAG(carr(ii,jj)) / (DBLE(carr(ii,jj)) + rho) )
4626         else
4627           theta= zero
4628         end if
4629         carr(ii,jj) = CMPLX(rho, theta*fact, kind=dpc)
4630      end do
4631    end do
4632 
4633  case default
4634    msg = " Wrong value for from: "//TRIM(from)
4635    ABI_BUG(msg)
4636  end select
4637 
4638 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

SOURCE

3097 subroutine coeffs_gausslegint(xmin,xmax,x,weights,n)
3098 
3099 !Arguments ------------------------------------
3100 !scalars
3101  integer,intent(in) :: n
3102  real(dp),intent(in) :: xmin,xmax
3103  real(dp),intent(out) :: x(n),weights(n)
3104 
3105 !Local variables ------------------------------
3106 !scalars
3107  integer :: i,j
3108  real(dp),parameter :: tol=1.d-13
3109  real(dp),parameter :: pi=4.d0*atan(1.d0)
3110  real(dp) :: z,z1,xmean,p1,p2,p3,pp,xl
3111 
3112 !************************************************************************
3113 
3114  xl=(xmax-xmin)*0.5d0
3115  xmean=(xmax+xmin)*0.5d0
3116 
3117  do i=1,(n+1)/2
3118   z=cos(pi*(i-0.25d0)/(n+0.5d0))
3119 
3120   do
3121     p1=1.d0
3122     p2=0.d0
3123 
3124     do j=1,n
3125      p3=p2
3126      p2=p1
3127      p1=((2.d0*j - 1.d0)*z*p2 - (j-1.d0)*p3)/j
3128     end do
3129 
3130     pp=n*(p2-z*p1)/(1.0d0-z**2)
3131     z1=z
3132     z=z1-p1/pp
3133 
3134     if(abs(z-z1) < tol) exit
3135   end do
3136 
3137   x(i)=xmean-xl*z
3138   x(n+1-i)=xmean+xl*z
3139   weights(i)=2.d0*xl/((1.d0-z**2)*pp**2)
3140   weights(n+1-i)=weights(i)
3141  end do
3142 
3143 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.

SOURCE

4476 subroutine continued_fract(nlev,term_type,aa,bb,nz,zpts,spectrum)
4477 
4478 !Arguments ------------------------------------
4479 !scalars
4480  integer,intent(in) :: nlev,term_type,nz
4481 !arrays
4482  real(dp),intent(in) ::  bb(nlev)
4483  complex(dpc),intent(in) :: aa(nlev)
4484  complex(dpc),intent(in) :: zpts(nz)
4485  complex(dpc),intent(out) :: spectrum(nz)
4486 
4487 !Local variables ------------------------------
4488 !scalars
4489  integer :: it
4490  real(dp) ::  bb_inf,bg,bu,swap
4491  complex(dpc) :: aa_inf
4492  character(len=500) :: msg
4493 !arrays
4494  complex(dpc),allocatable :: div(:),den(:)
4495 
4496 !************************************************************************
4497 
4498  ABI_MALLOC(div,(nz))
4499  ABI_MALLOC(den,(nz))
4500 
4501  select case (term_type)
4502  case (0) ! No terminator.
4503    div=czero
4504  case (-1,1)
4505    if (term_type==-1) then
4506      bb_inf=bb(nlev)
4507      aa_inf=aa(nlev)
4508    else
4509      bb_inf=SUM(bb)/nlev
4510      aa_inf=SUM(aa)/nlev
4511    end if
4512    ! Be careful with the sign of the SQRT.
4513    div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
4514  case (2)
4515    ABI_ERROR("To be tested")
4516    div = zero
4517    if (nlev>4) then
4518      bg=zero; bu=zero
4519      do it=1,nlev,2
4520        if (it+2<nlev) bg = bg + bb(it+2)
4521        bu = bu + bb(it)
4522      end do
4523      bg = bg/(nlev/2+MOD(nlev,2))
4524      bu = bg/((nlev+1)/2)
4525      !if (iseven(nlev)) then
4526      if (.not.iseven(nlev)) then
4527        swap = bg
4528        bg = bu
4529        bu = bg
4530      end if
4531      !write(std_out,*)nlev,bg,bu
4532      !Here be careful with the sign of SQRT
4533      do it=1,nz
4534        div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
4535 &        ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
4536      end do
4537    end if
4538 
4539  case default
4540    write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
4541    ABI_ERROR(msg)
4542  end select
4543 
4544  do it=nlev,2,-1
4545    den(:) = zpts(:) - aa(it) - div(:)
4546    div(:) = (bb(it-1)**2 )/ den(:)
4547  end do
4548 
4549  den = zpts(:) - aa(1) - div(:)
4550  div = one/den(:)
4551 
4552  spectrum = div
4553  ABI_FREE(div)
4554  ABI_FREE(den)
4555 
4556 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.

SOURCE

2943 subroutine cspint ( ftab, xtab, ntab, a, b, y, e, work, result )
2944 
2945 !Arguments ------------------------------------
2946 !scalars
2947   integer, intent(in) :: ntab
2948   real(dp), intent(in) :: a
2949   real(dp), intent(in) :: b
2950   real(dp), intent(inout) :: e(ntab)
2951   real(dp), intent(in) :: ftab(ntab)
2952   real(dp), intent(inout) :: work(ntab)
2953   real(dp), intent(in) :: xtab(ntab)
2954   real(dp), intent(inout) :: y(3,ntab)
2955   real(dp), intent(out) :: result
2956 
2957 !Local variables ------------------------------
2958 !scalars
2959   integer :: i
2960   integer :: j
2961   real(dp) :: r
2962   real(dp) :: s
2963   real(dp) :: term
2964   real(dp) :: u
2965 
2966 !************************************************************************
2967 
2968   if ( ntab < 3 ) then
2969     write(std_out,'(a)' ) ' '
2970     write(std_out,'(a)' ) 'CSPINT - Fatal error!'
2971     write(std_out,'(a,i6)' ) '  NTAB must be at least 3, but input NTAB = ',ntab
2972     ABI_ERROR("Aborting now")
2973   end if
2974 
2975   do i = 1, ntab-1
2976 
2977     if ( xtab(i+1) <= xtab(i) ) then
2978       write(std_out,'(a)' ) ' '
2979       write(std_out,'(a)' ) 'CSPINT - Fatal error!'
2980       write(std_out,'(a)' ) '  Nodes not in strict increasing order.'
2981       write(std_out,'(a,i6)' ) '  XTAB(I) <= XTAB(I-1) for I=',i
2982       write(std_out,'(a,g14.6)' ) '  XTAB(I) = ',xtab(i)
2983       write(std_out,'(a,g14.6)' ) '  XTAB(I-1) = ',xtab(i-1)
2984       ABI_ERROR("Aborting now")
2985     end if
2986 
2987   end do
2988 
2989   s = zero
2990   do i = 1, ntab-1
2991     r = ( ftab(i+1) - ftab(i) ) / ( xtab(i+1) - xtab(i) )
2992     y(2,i) = r - s
2993     s = r
2994   end do
2995 
2996   result = zero
2997   s = zero
2998   r = zero
2999   y(2,1) = zero
3000   y(2,ntab) = zero
3001 
3002   do i = 2, ntab-1
3003     y(2,i) = y(2,i) + r * y(2,i-1)
3004     work(i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s
3005     s = xtab(i+1) - xtab(i)
3006     r = s / work(i)
3007   end do
3008 
3009   do j = 2, ntab-1
3010     i = ntab+1-j
3011     y(2,i) = ( ( xtab(i+1) - xtab(i) ) * y(2,i+1) - y(2,i) ) / work(i)
3012   end do
3013 
3014   do i = 1, ntab-1
3015     s = xtab(i+1) - xtab(i)
3016     r = y(2,i+1) - y(2,i)
3017     y(3,i) = r / s
3018     y(2,i) = three * y(2,i)
3019     y(1,i) = ( ftab(i+1) - ftab(i) ) / s - ( y(2,i) + r ) * s
3020   end do
3021 
3022   e(1) = 0.0D+00
3023   do i = 1, ntab-1
3024     s = xtab(i+1)-xtab(i)
3025     term = ((( y(3,i) * quarter * s + y(2,i) * third ) * s &
3026       + y(1,i) * half ) * s + ftab(i) ) * s
3027     e(i+1) = e(i) + term
3028   end do
3029 !
3030 !  Determine where the endpoints A and B lie in the mesh of XTAB's.
3031 !
3032   r = a
3033   u = one
3034 
3035   do j = 1, 2
3036 !
3037 !  The endpoint is less than or equal to XTAB(1).
3038 !
3039     if ( r <= xtab(1) ) then
3040       result = result-u*((r-xtab(1))*y(1,1)*half +ftab(1))*(r-xtab(1))
3041 !
3042 !  The endpoint is greater than or equal to XTAB(NTAB).
3043 !
3044     else if ( xtab(ntab) <= r ) then
3045 
3046       result = result -u * ( e(ntab) + ( r - xtab(ntab) ) &
3047         * ( ftab(ntab) + half * ( ftab(ntab-1) &
3048         + ( xtab(ntab) - xtab(ntab-1) ) * y(1,ntab-1) ) &
3049         * ( r - xtab(ntab) )))
3050 !
3051 !  The endpoint is strictly between XTAB(1) and XTAB(NTAB).
3052 !
3053     else
3054 
3055       do i = 1, ntab-1
3056 
3057         if ( r <= xtab(i+1) ) then
3058           r = r-xtab(i)
3059           result = result-u*(e(i)+(((y(3,i)*quarter*r+y(2,i)*third)*r &
3060             +y(1,i)*half )*r+ftab(i))*r)
3061           go to 120
3062         end if
3063 
3064       end do
3065 
3066     end if
3067 
3068   120   continue
3069 
3070     u = -one
3071     r = b
3072 
3073   end do
3074 
3075 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

SOURCE

2791 subroutine ctrap(imax,ff,hh,ans)
2792 
2793 !Arguments ------------------------------------
2794 !scalars
2795  integer,intent(in) :: imax
2796  real(dp),intent(in) :: hh
2797  real(dp),intent(out) :: ans
2798 !arrays
2799  real(dp),intent(in) :: ff(imax)
2800 
2801 !Local variables-------------------------------
2802 !scalars
2803  integer :: ir,ir2
2804  real(dp) :: endpt,sum
2805 
2806 ! *************************************************************************
2807 
2808  if (imax>=10)then
2809 
2810 !  endpt=end point correction terms (low and high ends)
2811    endpt  = (23.75d0*(ff(1)+ff(imax  )) &
2812 &   + 95.10d0*(ff(2)+ff(imax-1)) &
2813 &   + 55.20d0*(ff(3)+ff(imax-2)) &
2814 &   + 79.30d0*(ff(4)+ff(imax-3)) &
2815 &   + 70.65d0*(ff(5)+ff(imax-4)))/ 72.d0
2816    ir2 = imax - 5
2817    sum=0.00d0
2818    if (ir2 > 5) then
2819      do ir=6,ir2
2820        sum = sum + ff(ir)
2821      end do
2822    end if
2823    ans = (sum + endpt ) * hh
2824 
2825  else if (imax>=8)then
2826    endpt  = (17.0d0*(ff(1)+ff(imax  )) &
2827 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2828 &   + 43.0d0*(ff(3)+ff(imax-2)) &
2829 &   + 49.0d0*(ff(4)+ff(imax-3)) )/ 48.d0
2830    sum=0.0d0
2831    if(imax==9)sum=ff(5)
2832    ans = (sum + endpt ) * hh
2833 
2834  else if (imax==7)then
2835    ans = (17.0d0*(ff(1)+ff(imax  )) &
2836 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2837 &   + 43.0d0*(ff(3)+ff(imax-2)) &
2838 &   + 50.0d0* ff(4)                )/ 48.d0  *hh
2839 
2840  else if (imax==6)then
2841    ans = (17.0d0*(ff(1)+ff(imax  )) &
2842 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2843 &   + 44.0d0*(ff(3)+ff(imax-2)) )/ 48.d0  *hh
2844 
2845  else if (imax==5)then
2846    ans = (     (ff(1)+ff(5)) &
2847 &   + four*(ff(2)+ff(4)) &
2848 &   + two * ff(3)         )/ three  *hh
2849 
2850  else if (imax==4)then
2851    ans = (three*(ff(1)+ff(4)) &
2852 &   + nine *(ff(2)+ff(3))  )/ eight  *hh
2853 
2854  else if (imax==3)then
2855    ans = (     (ff(1)+ff(3)) &
2856 &   + four* ff(2)         )/ three *hh
2857 
2858  else if (imax==2)then
2859    ans = (ff(1)+ff(2))/ two  *hh
2860 
2861  else if (imax==1)then
2862    ans = ff(1)*hh
2863 
2864  end if
2865 
2866 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

4375 integer function denominator(dd,ierr,tolerance)
4376 
4377 !Arguments ------------------------------------
4378 !scalars
4379  integer,intent(out) :: ierr
4380  real(dp),intent(in) :: dd
4381  real(dp),optional,intent(in) :: tolerance
4382 
4383 !Local variables ------------------------------
4384 !scalars
4385  integer,parameter :: largest_integer = HUGE(1)
4386  integer :: ii
4387  real(dp) :: my_tol
4388 
4389 !************************************************************************
4390 
4391  ii=1
4392  my_tol=0.0001 ; if (PRESENT(tolerance)) my_tol=ABS(tolerance)
4393  do
4394    if (ABS(dd*ii-NINT(dd*ii))<my_tol) then
4395      denominator=ii
4396      ierr=0
4397      RETURN
4398    end if
4399    ! Handle the case in which dd is not rational within my_tol.
4400    if (ii==largest_integer) then
4401      denominator=ii
4402      ierr=-1
4403      RETURN
4404    end if
4405    ii=ii+1
4406  end do
4407 
4408 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

4073 function dpade(n, z, f, zz)
4074 
4075 !Arguments ------------------------------------
4076 !scalars
4077  integer,intent(in) :: n
4078  complex(dpc),intent(in) :: zz
4079  complex(dpc) :: dpade
4080 !arrays
4081  complex(dpc),intent(in) :: z(n),f(n)
4082 
4083 !Local variables-------------------------------
4084 !scalars
4085  integer :: i
4086 !arrays
4087  complex(dpc) :: a(n)
4088  complex(dpc) :: Az(0:n), Bz(0:n)
4089  complex(dpc) :: dAz(0:n), dBz(0:n)
4090 ! *************************************************************************
4091 
4092  call calculate_pade_a(a, n, z, f)
4093 
4094  Az(0)=czero
4095  Az(1)=a(1)
4096  Bz(0)=cone
4097  Bz(1)=cone
4098  dAz(0)=czero
4099  dAz(1)=czero
4100  dBz(0)=czero
4101  dBz(1)=czero
4102 
4103  do i=1,n-1
4104    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4105    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4106    dAz(i+1)=dAz(i)+a(i+1)*Az(i-1)+(zz-z(i))*a(i+1)*dAz(i-1)
4107    dBz(i+1)=dBz(i)+a(i+1)*Bz(i-1)+(zz-z(i))*a(i+1)*dBz(i-1)
4108  end do
4109  !write(std_out,*) 'Bz(n)', Bz(n)
4110  !if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) 'Bz(n)',Bz(n)
4111  !pade_approx = Az(n) / Bz(n)
4112  dpade=dAz(n)/Bz(n) -Az(n)*dBz(n)/(Bz(n)*Bz(n))
4113  !write(std_out,*) 'pade_approx ', pade_approx
4114 
4115 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

SOURCE

5790 subroutine findmin(dedv_1,dedv_2,dedv_predict,&
5791 & d2edv2_1,d2edv2_2,d2edv2_predict,&
5792 & etotal_1,etotal_2,etotal_predict,&
5793 & lambda_1,lambda_2,lambda_predict,status)
5794 
5795 !Arguments ------------------------------------
5796 !scalars
5797  integer,intent(out) :: status
5798  real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2
5799  real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict
5800  real(dp),intent(out) :: etotal_predict,lambda_predict
5801 
5802 !Local variables-------------------------------
5803 !scalars
5804  real(dp) :: aa,bb,bbp,cc,ccp,d_lambda,dd
5805  real(dp) :: discr,ee,eep,lambda_shift,sum1,sum2,sum3,uu
5806  real(dp) :: uu3,vv,vv3
5807  character(len=500) :: msg
5808 
5809 ! *************************************************************************
5810 
5811 !DEBUG
5812 !write(std_out,*)' findmin : enter'
5813 !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2
5814 !ENDDEBUG
5815 
5816  status=0
5817  d_lambda=lambda_1-lambda_2
5818 
5819 !DEBUG
5820 !do choice=3,1,-1
5821 !ENDDEBUG
5822 
5823  if(abs(lambda_1-1.0_dp)>tol12 .or. abs(lambda_2)>tol12) then
5824    ABI_BUG('For choice=4, lambda_1 must be 1 and lambda_2 must be 0.')
5825  end if
5826 
5827 !Evaluate quartic interpolation
5828 !etotal = aa + bb * lambda + cc * lambda**2 + dd * lambda**3 + ee * lambda**4
5829 !Impose positive second derivative everywhere, with
5830 !one point where it vanishes :  3*dd**2=8*cc*ee
5831  aa=etotal_2
5832  bb=dedv_2
5833  sum1=etotal_1-aa-bb
5834  sum2=dedv_1-bb
5835  sum3=sum2-2.0_dp*sum1
5836 
5837 !Build the discriminant of the associated 2nd degree equation
5838  discr=sum2**2-3.0_dp*sum3**2
5839  if(discr<0.0_dp .or. sum2<0.0_dp)then
5840 
5841 ! jmb init
5842    d2edv2_2=0.0
5843    d2edv2_1=0.0
5844    d2edv2_predict=0.0
5845 
5846 !  Even if there is a problem, try to keep going ...
5847    ABI_WARNING('The 2nd degree equation has no positive root (choice=4).')
5848    status=2
5849    if(etotal_1<etotal_2)then
5850      write(msg, '(a,a,a)' )&
5851       'Will continue, since the new total energy is lower',ch10,&
5852       'than the old. Take a larger step in the same direction.'
5853      ABI_COMMENT(msg)
5854      lambda_predict=2.5_dp
5855    else
5856      write(msg, '(a,a,a,a,a)' )&
5857      'There is a problem, since the new total energy is larger',ch10,&
5858      'than the old (choice=4).',ch10,&
5859      'I take a point between the old and new, close to the old .'
5860      ABI_COMMENT(msg)
5861      lambda_predict=0.25_dp
5862    end if
5863 !  Mimick a zero-gradient lambda, in order to avoid spurious
5864 !  action of the inverse hessian (the next line would be a realistic estimation)
5865    dedv_predict=0.0_dp
5866 !  dedv_predict=dedv_2+lambda_predict*(dedv_1-dedv_2)
5867 !  Uses the energies, and the gradient at lambda_2
5868    etotal_predict=etotal_2+dedv_2*lambda_predict&
5869 &   +(etotal_1-etotal_2-dedv_2)*lambda_predict**2
5870 
5871  else
5872 
5873 !  Here, there is an acceptable solution to the 2nd degree equation
5874    discr=sqrt(discr)
5875 !  The root that gives the smallest ee corresponds to  -discr
5876 !  This is the one to be used: one aims at modelling the
5877 !  behaviour of the function as much as possible with the
5878 !  lowest orders of the polynomial, not the quartic term.
5879    ee=(sum2-discr)*0.5_dp
5880    dd=sum3-2.0_dp*ee
5881    cc=sum1-dd-ee
5882 
5883 !  DEBUG
5884 !  write(std_out,*)'aa,bb,cc,dd,ee',aa,bb,cc,dd,ee
5885 !  ENDDEBUG
5886 
5887 !  Now, must find the unique root of
5888 !  0 = bb + 2*cc * lambda + 3*dd * lambda^2 + 4*ee * lambda^3
5889 !  This root is unique because it was imposed that the second derivative
5890 !  of the quartic polynomial is everywhere positive.
5891 !  First, remove the quadratic term, by a shift of lambda
5892 !  lambdap=lambda-lambda_shift
5893 !  0 = bbp + ccp * lambdap + eep * lambdap^3
5894    eep=4.0_dp*ee
5895    lambda_shift=-dd/(4.0_dp*ee)
5896    ccp=2.0_dp*cc-12.0_dp*ee*lambda_shift**2
5897    bbp=bb+ccp*lambda_shift+eep*lambda_shift**3
5898 
5899 !  DEBUG
5900 !  write(std_out,*)'bbp,ccp,eep,lambda_shift',bbp,ccp,eep,lambda_shift
5901 !  ENDDEBUG
5902 
5903 !  The solution of a cubic polynomial equation is as follows :
5904    discr=(bbp/eep)**2+(4.0_dp/27.0_dp)*(ccp/eep)**3
5905 !  In the present case, discr will always be positive
5906    discr=sqrt(discr)
5907    uu3=0.5_dp*(-bbp/eep+discr) ; uu=sign((abs(uu3))**(1.0_dp/3.0_dp),uu3)
5908    vv3=0.5_dp*(-bbp/eep-discr) ; vv=sign((abs(vv3))**(1.0_dp/3.0_dp),vv3)
5909    lambda_predict=uu+vv
5910 
5911 !  Restore the shift
5912    lambda_predict=lambda_predict+lambda_shift
5913    etotal_predict=aa+bb*lambda_predict+cc*lambda_predict**2+&
5914 &   dd*lambda_predict**3+ee*lambda_predict**4
5915    dedv_predict=bb+2.0_dp*cc*lambda_predict+3.0_dp*dd*lambda_predict**2+&
5916 &   4.0_dp*ee*lambda_predict**3
5917    d2edv2_1=2*cc+6*dd*lambda_1+12*ee*lambda_1**2
5918    d2edv2_2=2*cc+6*dd*lambda_2+12*ee*lambda_2**2
5919    d2edv2_predict=2*cc+6*dd*lambda_predict+12*ee*lambda_predict**2
5920 
5921  end if
5922 
5923  write(msg, '(a,i3)' )'   line minimization, algorithm ',4
5924  call wrtout(std_out,msg,'COLL')
5925  write(msg, '(a,a)' )'                        lambda      etotal ','           dedv        d2edv2    '
5926  call wrtout(std_out,msg,'COLL')
5927  write(msg, '(a,es12.4,es18.10,2es12.4)' )'   old point         :',lambda_2,etotal_2,dedv_2,d2edv2_2
5928  call wrtout(std_out,msg,'COLL')
5929  write(msg, '(a,es12.4,es18.10,2es12.4)' )'   new point         :',lambda_1,etotal_1,dedv_1,d2edv2_1
5930  call wrtout(std_out,msg,'COLL')
5931  write(msg, '(a,es12.4,es18.10,2es12.4)' )'   predicted point   :',lambda_predict,etotal_predict,dedv_predict,d2edv2_predict
5932  call wrtout(std_out,msg,'COLL')
5933  write(msg, '(a)' ) ' '
5934  call wrtout(std_out,msg,'COLL')
5935 
5936 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

443 pure function geop(start,factor,nn) result(res)
444 
445 !Arguments ------------------------------------
446 !scalars
447  real(dp),intent(in) :: start,factor
448  integer,intent(in) :: nn
449  real(dp) :: res(nn)
450 
451 !Local variables-------------------------------
452  integer :: ii
453 ! *********************************************************************
454 
455  if (nn>0) res(1)=start
456  do ii=2,nn
457    res(ii)=res(ii-1)*factor
458  end do
459 
460 end function geop

m_numeric_tools/get_diag_cdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_cdp

FUNCTION

  Return the diagonal of a square matrix as a vector

SOURCE

804 function get_diag_cdp(cmat) result(cdiag)
805 
806 !Arguments ------------------------------------
807 !scalars
808  complex(dpc),intent(in) :: cmat(:,:)
809  complex(dpc) :: cdiag(SIZE(cmat,1))
810 
811 !Local variables-------------------------------
812  integer :: ii
813 ! *************************************************************************
814 
815  ABI_CHECK(SIZE(cmat,1) == SIZE(cmat,2), 'Matrix not square')
816 
817  do ii=1,SIZE(cmat,1)
818    cdiag(ii)=cmat(ii,ii)
819  end do
820 
821 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

SOURCE

742 function get_diag_int(mat) result(diag)
743 
744 !Arguments ------------------------------------
745 !scalars
746  integer,intent(in) :: mat(:,:)
747  integer :: diag(SIZE(mat,1))
748 
749 !Local variables-------------------------------
750  integer :: ii
751 ! *************************************************************************
752 
753  ii = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
754 
755  do ii=1,SIZE(mat,1)
756    diag(ii)=mat(ii,ii)
757  end do
758 
759 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

SOURCE

773 function get_diag_rdp(mat) result(diag)
774 
775 !Arguments ------------------------------------
776 !scalars
777  real(dp),intent(in) :: mat(:,:)
778  real(dp) :: diag(SIZE(mat,1))
779 
780 !Local variables-------------------------------
781  integer :: ii
782 ! *************************************************************************
783 
784  ABI_CHECK(SIZE(mat,1) == SIZE(mat,2), 'Matrix not square')
785 
786  do ii=1,SIZE(mat,1)
787    diag(ii) = mat(ii,ii)
788  end do
789 
790 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

714 pure function get_trace_cdp(matrix) result(trace)
715 
716 !Arguments ------------------------------------
717  complex(dpc) :: trace
718  complex(dpc),intent(in) :: matrix(:,:)
719 
720 !Local variables-------------------------------
721 !scalars
722  integer :: ii
723 ! *********************************************************************
724 
725  trace=czero
726  do ii=1,size(matrix,dim=1)
727     trace=trace+matrix(ii,ii)
728  end do
729 
730 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

644 pure function get_trace_int(matrix) result(trace)
645 
646 !Arguments ------------------------------------
647  integer :: trace
648  integer,intent(in) :: matrix(:,:)
649 
650 !Local variables-------------------------------
651 !scalars
652  integer :: ii
653 ! *********************************************************************
654 
655  trace=0
656  do ii=1,size(matrix,dim=1)
657    trace=trace+matrix(ii,ii)
658  end do
659 
660 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

680 pure function get_trace_rdp(matrix) result(trace)
681 
682 !Arguments ------------------------------------
683  real(dp) :: trace
684  real(dp),intent(in) :: matrix(:,:)
685 
686 !Local variables-------------------------------
687 !scalars
688  integer :: ii
689 ! *********************************************************************
690 
691  trace=zero
692  do ii=1,size(matrix,dim=1)
693     trace=trace+matrix(ii,ii)
694  end do
695 
696 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

SOURCE

3456 subroutine hermit(chmin, chmout, ierr, ndim)
3457 
3458 !Arguments ------------------------------------
3459 !scalars
3460  integer,intent(in) :: ndim
3461  integer,intent(out) :: ierr
3462 !arrays
3463  real(dp),intent(inout) :: chmin(ndim*ndim+ndim)
3464  real(dp),intent(inout) :: chmout(ndim*ndim+ndim)
3465 
3466 !Local variables-------------------------------
3467 !scalars
3468  integer,save :: mmesgs=20,nmesgs=0
3469  integer :: idim,merrors,nerrors
3470  real(dp),parameter :: eps=epsilon(0.0d0)
3471  real(dp) :: ch_im,ch_re,moduls,tol
3472  character(len=500) :: msg
3473 
3474 ! *************************************************************************
3475 
3476  tol=4096.0d0*eps
3477 
3478  ierr=0
3479  merrors=0
3480 
3481 !Copy matrix into possibly new location
3482  chmout(:)=chmin(:)
3483 
3484 !Loop over diagonal elements of matrix (off-diag not altered)
3485  do idim=1,ndim
3486 
3487    ch_im=chmout(idim*idim+idim  )
3488    ch_re=chmout(idim*idim+idim-1)
3489 
3490 !  check for large absolute Im part and print warning when
3491 !  larger than (some factor)*(machine precision)
3492    nerrors=0
3493    if( abs(ch_im) > tol .and. abs(ch_im) > tol8*abs(ch_re)) nerrors=2
3494    if( abs(ch_im) > tol .or. abs(ch_im) > tol8*abs(ch_re)) nerrors=1
3495 
3496    if( (abs(ch_im) > tol .and. nmesgs<mmesgs) .or. nerrors==2)then
3497      write(msg, '(3a,i0,a,es20.12,a,es20.12,a)' )&
3498      ' Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
3499      ' for component:',idim,' Im part is: ',ch_im,', Re part is: ',ch_re,'.'
3500      call wrtout(std_out,msg)
3501      nmesgs=nmesgs+1
3502    end if
3503 
3504    if( ( abs(ch_im) > tol8*abs(ch_re) .and. nmesgs<mmesgs) .or. nerrors==2)then
3505      write(msg, '(3a,i0,a,es20.12,a,es20.12,a)' )&
3506      ' Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
3507      ' for component',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
3508      call wrtout(std_out,msg)
3509      nmesgs=nmesgs+1
3510    end if
3511 
3512 !  compute modulus $= (\Re^2+\Im^2)^{1/2}$
3513    moduls=sqrt(ch_re**2+ch_im**2)
3514 
3515 !  set Re part to modulus with sign of original Re part
3516    chmout(idim*idim+idim-1)=sign(moduls,ch_re)
3517 
3518 !  set Im part to 0
3519    chmout(idim*idim+idim)=zero
3520 
3521    merrors=max(merrors,nerrors)
3522 
3523  end do
3524 
3525  if(merrors==2)then
3526    ierr=1
3527    write(msg, '(3a)' )&
3528     'Imaginary part(s) of diagonal Hermitian matrix element(s) is too large.',ch10,&
3529     'See previous messages.'
3530    ABI_BUG(msg)
3531  end if
3532 
3533 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

SOURCE

3321 subroutine hermitianize_dpc(mat, uplo)
3322 
3323 !Arguments ------------------------------------
3324 !scalars
3325  character(len=*),intent(in) :: uplo
3326 !arrays
3327  complex(dpc),intent(inout) :: mat(:,:)
3328 
3329 !Local variables-------------------------------
3330  integer :: nn,ii,jj
3331  complex(dpc),allocatable :: tmp(:)
3332 ! *************************************************************************
3333 
3334  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3335 
3336  select case (uplo(1:1))
3337 
3338  case ("A","a")
3339    ! Full matrix has been calculated.
3340    ABI_MALLOC(tmp,(nn))
3341    do ii=1,nn
3342      do jj=ii,nn
3343        tmp(jj)=half*(mat(ii,jj)+DCONJG(mat(jj,ii)))
3344      end do
3345      mat(ii,ii:nn)=tmp(ii:nn)
3346      mat(ii:nn,ii)=DCONJG(tmp(ii:nn))
3347    end do
3348    ABI_FREE(tmp)
3349 
3350  case ("U","u")
3351   ! Only the upper triangle is used.
3352    do jj=1,nn
3353      do ii=1,jj
3354        if (ii/=jj) then
3355          mat(jj,ii) = DCONJG(mat(ii,jj))
3356        else
3357          mat(ii,ii) = CMPLX(DBLE(mat(ii,ii)),zero, kind=dpc)
3358        end if
3359      end do
3360    end do
3361 
3362  case ("L","l")
3363   ! Only the lower triangle is used.
3364   do jj=1,nn
3365     do ii=1,jj
3366       if (ii/=jj) then
3367         mat(ii,jj) = DCONJG(mat(jj,ii))
3368       else
3369         mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),zero, kind=dpc)
3370       end if
3371     end do
3372   end do
3373 
3374  case default
3375    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3376  end select
3377 
3378 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

SOURCE

3236 subroutine hermitianize_spc(mat,uplo)
3237 
3238 !Arguments ------------------------------------
3239 !scalars
3240  character(len=*),intent(in) :: uplo
3241 !arrays
3242  complex(spc),intent(inout) :: mat(:,:)
3243 
3244 !Local variables-------------------------------
3245  integer :: nn,ii,jj
3246  complex(spc),allocatable :: tmp(:)
3247 ! *************************************************************************
3248 
3249  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3250 
3251  select case (uplo(1:1))
3252 
3253  case ("A","a")
3254    ! Full matrix has been calculated.
3255    ABI_MALLOC(tmp,(nn))
3256    do ii=1,nn
3257      do jj=ii,nn
3258        ! reference half constant is dp not sp
3259        tmp(jj)=real(half)*(mat(ii,jj)+CONJG(mat(jj,ii)))
3260      end do
3261      mat(ii,ii:nn)=tmp(ii:nn)
3262      mat(ii:nn,ii)=CONJG(tmp(ii:nn))
3263    end do
3264    ABI_FREE(tmp)
3265 
3266  case ("U","u")
3267    ! Only the upper triangle is used.
3268    do jj=1,nn
3269      do ii=1,jj
3270        if (ii/=jj) then
3271          mat(jj,ii) = CONJG(mat(ii,jj))
3272        else
3273          mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3274        end if
3275      end do
3276    end do
3277 
3278  case ("L","l")
3279   ! Only the lower triangle is used.
3280   do jj=1,nn
3281     do ii=1,jj
3282       if (ii/=jj) then
3283         mat(ii,jj) = CONJG(mat(jj,ii))
3284       else
3285         mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3286       end if
3287     end do
3288   end do
3289 
3290  case default
3291    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3292  end select
3293 
3294 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

1683 pure function imax_loc_int(iarr,mask)
1684 
1685 !Arguments ------------------------------------
1686 !scalars
1687  integer :: imax_loc_int
1688 !arrays
1689  integer,intent(in) :: iarr(:)
1690  logical,optional,intent(in) :: mask(:)
1691 
1692 !Local variables-------------------------------
1693  integer :: imax(1)
1694 ! *************************************************************************
1695 
1696  if (PRESENT(mask)) then
1697   imax=MAXLOC(iarr,MASK=mask)
1698  else
1699   imax=MAXLOC(iarr)
1700  end if
1701  imax_loc_int=imax(1)
1702 
1703 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

1719 pure function imax_loc_rdp(arr,mask)
1720 
1721 !Arguments ------------------------------------
1722 !scalars
1723  integer :: imax_loc_rdp
1724 !arrays
1725  real(dp),intent(in) :: arr(:)
1726  logical,optional,intent(in) :: mask(:)
1727 
1728 !Local variables-------------------------------
1729  integer :: imax(1)
1730 ! *************************************************************************
1731 
1732  if (PRESENT(mask)) then
1733   imax=MAXLOC(arr,MASK=mask)
1734  else
1735   imax=MAXLOC(arr)
1736  end if
1737  imax_loc_rdp=imax(1)
1738 
1739 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

1753 pure function imin_loc_int(arr, mask)
1754 
1755 !Arguments ------------------------------------
1756 !scalars
1757  integer :: imin_loc_int
1758 !arrays
1759  integer,intent(in) :: arr(:)
1760  logical,optional,intent(in) :: mask(:)
1761 
1762 !Local variables-------------------------------
1763  integer :: imin(1)
1764 ! *************************************************************************
1765 
1766  if (PRESENT(mask)) then
1767   imin=MINLOC(arr,MASK=mask)
1768  else
1769   imin=MINLOC(arr)
1770  end if
1771  imin_loc_int=imin(1)
1772 
1773 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

1790 pure function imin_loc_rdp(arr,mask)
1791 
1792 !Arguments ------------------------------------
1793 !scalars
1794  integer :: imin_loc_rdp
1795 !arrays
1796  real(dp),intent(in) :: arr(:)
1797  logical,optional,intent(in) :: mask(:)
1798 
1799 !Local variables-------------------------------
1800  integer :: imin(1)
1801 ! *************************************************************************
1802 
1803  if (PRESENT(mask)) then
1804   imin=MINLOC(arr,MASK=mask)
1805  else
1806   imin=MINLOC(arr)
1807  end if
1808 
1809  imin_loc_rdp=imin(1)
1810 
1811 end function imin_loc_rdp

m_numeric_tools/inrange_dp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  inrange_dp

FUNCTION

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

SOURCE

1555 pure logical function inrange_dp(xval, win)
1556 
1557 !Arguments ------------------------------------
1558 !scalars
1559  real(dp),intent(in) :: xval, win(2)
1560 ! *************************************************************************
1561 
1562  inrange_dp = (xval >= win(1) .and. xval <= win(2))
1563 
1564 end function inrange_dp

m_numeric_tools/inrange_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  inrange_int

FUNCTION

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

SOURCE

1532 pure logical function inrange_int(xval, win)
1533 
1534 !Arguments ------------------------------------
1535 !scalars
1536  integer,intent(in) :: xval,win(2)
1537 ! *************************************************************************
1538 
1539  inrange_int = (xval >= win(1) .and. xval <= win(2))
1540 
1541 end function inrange_int

m_numeric_tools/interpol3d_0d [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 interpol3d_0d

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

SOURCE

4929 pure function interpol3d_0d(r, nr1, nr2, nr3, grid) result(res)
4930 
4931 !Arguments-------------------------------------------------------------
4932 !scalars
4933  integer,intent(in) :: nr1, nr2, nr3
4934  real(dp) :: res
4935 !arrays
4936  real(dp),intent(in) :: grid(nr1,nr2,nr3),r(3)
4937 
4938 !Local variables--------------------------------------------------------
4939 !scalars
4940  integer :: ir1,ir2,ir3,pr1,pr2,pr3
4941  real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8
4942  real(dp) :: x1,x2,x3
4943 
4944 ! *************************************************************************
4945 
4946  call interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3, pr1,pr2,pr3)
4947 
4948 !weight
4949  x1=one+r(1)*nr1-real(ir1)
4950  x2=one+r(2)*nr2-real(ir2)
4951  x3=one+r(3)*nr3-real(ir3)
4952 
4953 !calculation of the density value
4954  res1=grid(ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
4955  res2=grid(pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
4956  res3=grid(ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
4957  res4=grid(ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
4958  res5=grid(pr1, pr2, ir3) * x1*x2*(one-x3)
4959  res6=grid(ir1, pr2, pr3) * (one-x1)*x2*x3
4960  res7=grid(pr1, ir2, pr3) * x1*(one-x2)*x3
4961  res8=grid(pr1, pr2, pr3) * x1*x2*x3
4962  res=res1+res2+res3+res4+res5+res6+res7+res8
4963 
4964 end function interpol3d_0d

m_numeric_tools/interpol3d_1d [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 interpol3d_1d

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(nd,nr1,nr2,nr3)=grid matrix

OUTPUT

 res(nd)=Interpolated value

SOURCE

4990 pure function interpol3d_1d(r, nr1, nr2, nr3, grid, nd) result(res)
4991 
4992 !Arguments-------------------------------------------------------------
4993 !scalars
4994  integer,intent(in) :: nr1, nr2, nr3, nd
4995  real(dp) :: res(nd)
4996 !arrays
4997  real(dp),intent(in) :: grid(nd,nr1,nr2,nr3),r(3)
4998 
4999 !Local variables--------------------------------------------------------
5000 !scalars
5001  integer :: id,ir1,ir2,ir3,pr1,pr2,pr3
5002  real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8
5003  real(dp) :: x1,x2,x3
5004 
5005 ! *************************************************************************
5006 
5007  call interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3, pr1,pr2,pr3)
5008 
5009 !weight
5010  x1=one+r(1)*nr1-real(ir1)
5011  x2=one+r(2)*nr2-real(ir2)
5012  x3=one+r(3)*nr3-real(ir3)
5013 
5014 !calculation of the density value
5015  do id=1,nd
5016    res1=grid(id,ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
5017    res2=grid(id,pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
5018    res3=grid(id,ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
5019    res4=grid(id,ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
5020    res5=grid(id,pr1, pr2, ir3) * x1*x2*(one-x3)
5021    res6=grid(id,ir1, pr2, pr3) * (one-x1)*x2*x3
5022    res7=grid(id,pr1, ir2, pr3) * x1*(one-x2)*x3
5023    res8=grid(id,pr1, pr2, pr3) * x1*x2*x3
5024    res(id)=res1+res2+res3+res4+res5+res6+res7+res8
5025  enddo
5026 
5027 end function interpol3d_1d

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

SOURCE

5052 pure subroutine interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3,pr1,pr2,pr3)
5053 
5054 !Arguments-------------------------------------------------------------
5055 !scalars
5056  integer,intent(in) :: nr1,nr2,nr3
5057  integer,intent(out) :: ir1,ir2,ir3
5058  integer,intent(out) :: pr1,pr2,pr3
5059 !arrays
5060  real(dp),intent(in) :: r(3)
5061 
5062 !Local variables-------------------------------
5063  real(dp) :: d1,d2,d3
5064 
5065 ! *************************************************************************
5066 
5067 !grid density
5068  d1=one/nr1
5069  d2=one/nr2
5070  d3=one/nr3
5071 
5072 !lower left
5073  ir1=int(r(1)/d1)+1
5074  ir2=int(r(2)/d2)+1
5075  ir3=int(r(3)/d3)+1
5076 
5077 !upper right
5078  pr1=mod(ir1+1,nr1)
5079  pr2=mod(ir2+1,nr2)
5080  pr3=mod(ir3+1,nr3)
5081 
5082  if(ir1==0) ir1=nr1
5083  if(ir2==0) ir2=nr2
5084  if(ir3==0) ir3=nr3
5085 
5086  if(ir1>nr1) ir1=ir1-nr1
5087  if(ir2>nr2) ir2=ir2-nr2
5088  if(ir3>nr3) ir3=ir3-nr3
5089 
5090  if(pr1==0) pr1=nr1
5091  if(pr2==0) pr2=nr2
5092  if(pr3==0) pr3=nr3
5093 
5094 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.

SOURCE

5119 subroutine interpolate_denpot(cplex, in_ngfft, nspden, in_rhor, out_ngfft, out_rhor)
5120 
5121 !Arguments-------------------------------------------------------------
5122 !scalars
5123  integer,intent(in) :: cplex,nspden
5124 !arrays
5125  integer,intent(in) :: in_ngfft(3), out_ngfft(3)
5126  real(dp),intent(in) :: in_rhor(cplex, product(in_ngfft), nspden)
5127  real(dp),intent(out) :: out_rhor(cplex, product(out_ngfft), nspden)
5128 
5129 !Local variables--------------------------------------------------------
5130 !scalars
5131  integer :: ispden, ir1, ir2, ir3, ifft
5132  real(dp) :: rr(3)
5133 
5134 ! *************************************************************************
5135 
5136  ! Linear interpolation.
5137  do ispden=1,nspden
5138    do ir3=0,out_ngfft(3)-1
5139      rr(3) = DBLE(ir3)/out_ngfft(3)
5140      do ir2=0,out_ngfft(2)-1
5141        rr(2) = DBLE(ir2)/out_ngfft(2)
5142        do ir1=0,out_ngfft(1)-1
5143          rr(1) = DBLE(ir1)/out_ngfft(1)
5144          ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5145          out_rhor(1:cplex, ifft, ispden) = interpol3d_1d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_rhor(:, :, ispden),cplex)
5146        end do
5147      end do
5148    end do
5149  end do
5150 
5151 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

SOURCE

6161 subroutine invcb(rhoarr,rspts,npts)
6162 
6163 !Arguments ------------------------------------
6164 !scalars
6165  integer,intent(in) :: npts
6166 !arrays
6167  real(dp),intent(in) :: rhoarr(npts)
6168  real(dp),intent(out) :: rspts(npts)
6169 
6170 !Local variables-------------------------------
6171 !scalars
6172  integer :: ii,ipts
6173  real(dp),parameter :: c2_27=2.0e0_dp/27.0e0_dp,c5_9=5.0e0_dp/9.0e0_dp
6174  real(dp),parameter :: c8_9=8.0e0_dp/9.0e0_dp,m1thrd=-third
6175  real(dp) :: del,prod,rho,rhom1,rhomtrd
6176  logical :: test
6177 !character(len=500) :: message
6178 
6179 ! *************************************************************************
6180 
6181 !Loop over points : here, brute force algorithm
6182 !do ipts=1,npts
6183 !rspts(ipts)=sign( (abs(rhoarr(ipts)))**m1thrd,rhoarr(ipts))
6184 !end do
6185 !
6186 
6187  rhomtrd=sign( (abs(rhoarr(1)))**m1thrd, rhoarr(1) )
6188  rhom1=one/rhoarr(1)
6189  rspts(1)=rhomtrd
6190  do ipts=2,npts
6191    rho=rhoarr(ipts)
6192    prod=rho*rhom1
6193 !  If the previous point is too far ...
6194    if(prod < 0.01_dp .or. prod > 10._dp )then
6195      rhomtrd=sign( (abs(rho))**m1thrd , rho )
6196      rhom1=one/rho
6197    else
6198      del=prod-one
6199      do ii=1,5
6200 !      Choose one of the two next lines, the last one is more accurate
6201 !      rhomtrd=((one+third*del)/(one+two_thirds*del))*rhomtrd
6202        rhomtrd=((one+c5_9*del)/(one+del*(c8_9+c2_27*del)))*rhomtrd
6203        rhom1=rhomtrd*rhomtrd*rhomtrd
6204        del=rho*rhom1-one
6205 !      write(std_out,*)rhomtrd,del
6206        test = del*del < 1.0e-24_dp
6207        if(test) exit
6208      end do
6209      if( .not. test) then
6210        rhomtrd=sign( (abs(rho))**m1thrd , rho )
6211      end if
6212    end if
6213    rspts(ipts)=rhomtrd
6214  end do
6215 
6216 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

1414 pure function is_integer_0d(rr,tol) result(ans)
1415 
1416 !Arguments ------------------------------------
1417 !scalars
1418  real(dp),intent(in) :: tol
1419  logical :: ans
1420 !arrays
1421  real(dp),intent(in) :: rr
1422 
1423 ! *************************************************************************
1424 
1425  ans=(ABS(rr-NINT(rr))<tol)
1426 
1427 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

1444 pure function is_integer_1d(rr,tol) result(ans)
1445 
1446 !Arguments ------------------------------------
1447 !scalars
1448  real(dp),intent(in) :: tol
1449  logical :: ans
1450 !arrays
1451  real(dp),intent(in) :: rr(:)
1452 
1453 ! *************************************************************************
1454 
1455  ans=ALL((ABS(rr-NINT(rr))<tol))
1456 
1457 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

SOURCE

1477 function is_zero_rdp_0d(rr,tol) result(ans)
1478 
1479 !Arguments ------------------------------------
1480 !scalars
1481  real(dp),intent(in) :: tol
1482  logical :: ans
1483 !arrays
1484  real(dp),intent(in) :: rr
1485 ! *************************************************************************
1486 
1487  ans=(ABS(rr)<tol)
1488 
1489 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

1506 function is_zero_rdp_1d(rr,tol) result(ans)
1507 
1508 !Arguments ------------------------------------
1509 !scalars
1510  real(dp),intent(in) :: tol
1511  logical :: ans
1512 !arrays
1513  real(dp),intent(in) :: rr(:)
1514 ! *************************************************************************
1515 
1516  ans=ALL(ABS(rr(:))<tol)
1517 
1518 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

835 pure logical function isdiagmat_int(mat) result(ans)
836 
837 !Arguments ------------------------------------
838 !scalars
839  integer,intent(in) :: mat(:,:)
840 
841 !Local variables-------------------------------
842  integer :: ii,jj
843 ! *************************************************************************
844 
845  ans = .True.
846  do jj=1,size(mat,dim=2)
847    do ii=1,size(mat,dim=1)
848      if (ii == jj) cycle
849      if (mat(ii,jj) /= 0) then
850        ans = .False.; return
851      end if
852    end do
853  end do
854 
855 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

869 pure logical function isdiagmat_rdp(mat, atol) result(ans)
870 
871 !Arguments ------------------------------------
872 !scalars
873  real(dp),intent(in) :: mat(:,:)
874  real(dp),optional,intent(in) :: atol
875 
876 !Local variables-------------------------------
877  integer :: ii,jj
878  real(dp) :: my_atol
879 ! *************************************************************************
880 
881  my_atol = tol12; if (present(atol)) my_atol = atol
882 
883  ans = .True.
884  do jj=1,size(mat,dim=2)
885    do ii=1,size(mat,dim=1)
886      if (ii == jj) cycle
887      if (abs(mat(ii,jj)) > my_atol) then
888        ans = .False.; return
889      end if
890    end do
891  end do
892 
893 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

SOURCE

1386 elemental function iseven(nn)
1387 
1388 !Arguments ------------------------------------
1389 !scalars
1390  integer,intent(in) :: nn
1391  logical :: iseven
1392 ! *********************************************************************
1393 
1394  iseven = ((nn / 2) * 2 == nn)
1395 
1396 end function iseven

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.

SOURCE

4717 function isordered_rdp(nn,arr,direction,tol) result(isord)
4718 
4719 !Arguments ------------------------------------
4720 !scalars
4721  integer,intent(in) :: nn
4722  real(dp),intent(in) :: tol
4723  logical :: isord
4724  character(len=*),intent(in) :: direction
4725 !arrays
4726  real(dp),intent(in) :: arr(nn)
4727 
4728 !Local variables ------------------------------
4729 !scalars
4730  integer :: ii
4731  real(dp) :: prev
4732  character(len=500) :: msg
4733 
4734 ! *************************************************************************
4735 
4736  prev = arr(1); isord =.TRUE.
4737 
4738  SELECT CASE (direction(1:1))
4739  CASE(">")
4740  ii=2;
4741  do while (ii<=nn .and. isord)
4742    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) >= prev)
4743    prev = arr(ii)
4744    ii = ii +1
4745  end do
4746 
4747  CASE("<")
4748  ii=2;
4749  do while (ii<=nn .and. isord)
4750    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) <= prev)
4751    prev = arr(ii)
4752    ii = ii +1
4753  end do
4754 
4755  CASE DEFAULT
4756    msg = "Wrong direction: "//TRIM(direction)
4757    ABI_ERROR(msg)
4758  END SELECT
4759 
4760 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

SOURCE

5968 subroutine kramerskronig(nomega,omega,eps,method,only_check)
5969 
5970 !Arguments ------------------------------------
5971 !scalars
5972  integer,intent(in) :: method,nomega,only_check
5973 !arrays
5974  real(dp),intent(in) :: omega(nomega)
5975  complex(dpc),intent(inout) :: eps(nomega)
5976 
5977 !Local variables-------------------------------
5978 !scalars
5979  integer,save :: enough=0
5980  integer :: ii,ip
5981  real(dp) :: acc,domega,eav,kkdif,kkrms,ww,wwp
5982  character(len=500) :: msg
5983 !arrays
5984  real(dp) :: e1kk(nomega),intkk(nomega),kk(nomega)
5985 
5986 ! *************************************************************************
5987 
5988 !Check whether the frequency grid is linear or not
5989  domega = (omega(nomega) - omega(1)) / (nomega-1)
5990  do ii=2,nomega
5991    if (ABS(domega-(omega(ii)-omega(ii-1))) > 0.001) then
5992      if (only_check/=1) then
5993        ABI_WARNING("Check cannot be performed since the frequency step is not constant")
5994        RETURN
5995      else
5996        ABI_ERROR('Cannot perform integration since frequency step is not constant')
5997      end if
5998    end if
5999  end do
6000 
6001 !Check whether omega(1) is small or not
6002  if (omega(1) > 0.1/Ha_eV) then
6003    if (only_check/=1) then
6004      ABI_WARNING('Check cannot be performed since first frequency on the grid > 0.1 eV')
6005      RETURN
6006    else
6007      ABI_ERROR('Cannot perform integration since first frequency on the grid > 0.1 eV')
6008    end if
6009  end if
6010 
6011 !If eps(nomega) is not 0 warn
6012  if (AIMAG(eps(nomega)) > 0.1 .and. enough<50) then
6013    enough=enough+1
6014    write(msg,'(a,f8.4,3a,f8.2,2a)')&
6015 &   'Im epsilon for omega = ',omega(nomega)*Ha_eV,' eV',ch10,&
6016 &   'is not yet zero, epsilon_2 = ',AIMAG(eps(nomega)),ch10,&
6017 &   'Kramers Kronig could give wrong results'
6018    ABI_WARNING(msg)
6019    if (enough==50) then
6020      write(msg,'(3a)')' sufficient number of WARNINGS-',ch10,' stop writing '
6021      call wrtout(std_out,msg,'COLL')
6022    end if
6023  end if
6024 
6025 
6026 !Perform Kramers-Kronig using naive integration
6027  select case (method)
6028  case (0)
6029 
6030    do ii=1,nomega
6031      ww = omega(ii)
6032      acc = 0.0_dp
6033      do ip=1,nomega
6034        if (ip == ii) CYCLE
6035        wwp = omega(ip)
6036        acc = acc + wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6037      end do
6038      e1kk(ii) = one + two/pi*domega* acc
6039    end do
6040 
6041 !    Perform Kramers-Kronig using Simpson integration
6042 !    Simpson O(1/N^4), from NumRec in C p 134  NumRec in Fortran p 128
6043  case (1)
6044 
6045    kk=zero
6046 
6047    do ii=1,nomega
6048      ww=omega(ii)
6049      do ip=1,nomega
6050        if (ip == ii) CYCLE
6051        wwp = omega(ip)
6052        kk(ip) = wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6053      end do
6054      call simpson_int(nomega,domega,kk,intkk)
6055      e1kk(ii) = one + two/pi * intkk(nomega)
6056    end do
6057 
6058  case default
6059    write(msg,'(a,i0)')' Wrong value for method ',method
6060    ABI_BUG(msg)
6061  end select
6062 
6063 !at this point real part is in e1kk, need to put it into eps
6064  do ii=1,nomega
6065    eps(ii)=CMPLX(e1kk(ii),AIMAG(eps(ii)), kind=dpc)
6066  end do
6067 
6068 !Verify Kramers-Kronig
6069  eav   = zero
6070  kkdif = zero
6071  kkrms = zero
6072 
6073  do ii=1,nomega
6074    kkdif = kkdif + ABS(REAL(eps(ii)) - e1kk(ii))
6075    kkrms = kkrms + (REAL(eps(ii)) - e1kk(ii))*(REAL(eps(ii)) - e1kk(ii))
6076    eav = eav + ABS(REAL(eps(ii)))
6077  end do
6078 
6079  eav = eav/nomega
6080  kkdif = (kkdif/nomega) / eav
6081  kkrms = (kkrms/nomega) / (eav*eav)
6082 
6083  kk = ABS(REAL(eps(1)) - e1kk(1)) / REAL(eps(1))
6084 
6085 !Write data
6086  write(msg,'(a,f7.2,a)')' Kramers-Kronig transform is verified within ',MAXVAL(kk)*100,"%"
6087  call wrtout(std_out,msg,'COLL')
6088 
6089 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

910 pure function l2int_1D(larr) result(int_arr)
911 
912 !Arguments ------------------------------------
913 !scalars
914  logical,intent(in) :: larr(:)
915  integer :: int_arr(size(larr))
916 
917 ! *********************************************************************
918 
919  where (larr)
920    int_arr = 1
921  elsewhere
922    int_arr = 0
923  end where
924 
925 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

942 pure function l2int_2D(larr) result(int_arr)
943 
944 !Arguments ------------------------------------
945 !scalars
946  logical,intent(in) :: larr(:,:)
947  integer :: int_arr(size(larr,1), size(larr,2))
948 
949 ! *********************************************************************
950 
951  where (larr)
952    int_arr = 1
953  elsewhere
954    int_arr = 0
955  end where
956 
957 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

974 pure function l2int_3D(larr) result(int_arr)
975 
976 !Arguments ------------------------------------
977 !scalars
978  logical,intent(in) :: larr(:,:,:)
979  integer :: int_arr(size(larr,1), size(larr,2), size(larr,3))
980 
981 ! *********************************************************************
982 
983  where (larr)
984    int_arr = 1
985  elsewhere
986    int_arr = 0
987  end where
988 
989 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

SOURCE

1830 integer pure function lfind(mask, back)
1831 
1832 !Arguments ------------------------------------
1833 !scalars
1834  logical,intent(in) :: mask(:)
1835  logical,optional,intent(in) :: back
1836 !arrays
1837 
1838 !Local variables-------------------------------
1839 !scalars
1840  integer :: ii,nitems
1841  logical :: do_back
1842 
1843 !************************************************************************
1844 
1845  do_back = .False.; if (present(back)) do_back = back
1846  lfind = -1; nitems = size(mask); if (nitems == 0) return
1847 
1848  if (do_back) then
1849    ! Backward search
1850    do ii=nitems,1,-1
1851      if (mask(ii)) then
1852        lfind = ii; return
1853      end if
1854    end do
1855  else
1856    ! Forward search.
1857    do ii=1,nitems
1858      if (mask(ii)) then
1859        lfind = ii; return
1860      end if
1861    end do
1862  end if
1863 
1864 end function lfind

m_numeric_tools/linfit_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linfit_dpc

FUNCTION

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

INPUTS

OUTPUT

SOURCE

2150 function linfit_dpc(nn,xx,zz,aa,bb) result(res)
2151 
2152 !Arguments ------------------------------------
2153 !scalars
2154  integer,intent(in) :: nn
2155  real(dp) :: res
2156  real(dp),intent(in) :: xx(nn)
2157  complex(dpc),intent(in) :: zz(nn)
2158  complex(dpc),intent(out) :: aa,bb
2159 !arrays
2160 
2161 !Local variables-------------------------------
2162 !scalars
2163  integer :: ii
2164  real(dp) :: sx,sx2,msrt
2165  complex(dpc) :: sz,sxz
2166 ! *************************************************************************
2167 
2168  sx=zero  ; sx2=zero ; msrt=zero
2169  sz=czero ; sxz=czero
2170  do ii=1,nn
2171   sx=sx+xx(ii)
2172   sz=sz+zz(ii)
2173   sxz=sxz+xx(ii)*zz(ii)
2174   sx2=sx2+xx(ii)*xx(ii)
2175  end do
2176 
2177  aa=(nn*sxz-sx*sz)/(nn*sx2-sx*sx)
2178  bb=sz/nn-sx*aa/nn
2179 
2180  do ii=1,nn
2181   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2182  end do
2183  msrt=SQRT(msrt) ; res=msrt
2184 
2185 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
  res=root mean square of differences between data and fit

SOURCE

2042 function linfit_rdp(nn,xx,yy,aa,bb) result(res)
2043 
2044 !Arguments ------------------------------------
2045 !scalars
2046  integer,intent(in) :: nn
2047  real(dp) :: res
2048  real(dp),intent(out) :: aa,bb
2049 !arrays
2050  real(dp),intent(in) :: xx(nn),yy(nn)
2051 
2052 !Local variables-------------------------------
2053 !scalars
2054  integer :: ii
2055  real(dp) :: msrt,sx2,sx,sxy,sy,tx,ty
2056 ! *************************************************************************
2057 
2058  sx=zero ; sy=zero ; sxy=zero ; sx2=zero
2059  do ii=1,nn
2060   tx=xx(ii)
2061   ty=yy(ii)
2062   sx=sx+tx
2063   sy=sy+ty
2064   sxy=sxy+tx*ty
2065   sx2=sx2+tx*tx
2066  end do
2067 
2068  aa=(nn*sxy-sx*sy)/(nn*sx2-sx*sx)
2069  bb=sy/nn-sx*aa/nn
2070 
2071  msrt=zero
2072  do ii=1,nn
2073   tx=xx(ii)
2074   ty=yy(ii)
2075   msrt=msrt+(ty-aa*tx-bb)**2
2076  end do
2077  msrt=SQRT(msrt/nn) ; res=msrt
2078 
2079 end function linfit_rdp

m_numeric_tools/linfit_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linfit_spc

FUNCTION

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

INPUTS

OUTPUT

SOURCE

2097 function linfit_spc(nn,xx,zz,aa,bb) result(res)
2098 
2099 !Arguments ------------------------------------
2100 !scalars
2101  integer,intent(in) :: nn
2102  real(dp) :: res
2103  real(dp),intent(in) :: xx(nn)
2104  complex(spc),intent(in) :: zz(nn)
2105  complex(spc),intent(out) :: aa,bb
2106 !arrays
2107 
2108 !Local variables-------------------------------
2109 !scalars
2110  integer :: ii
2111  real(dp) :: sx,sx2,msrt
2112  complex(dpc) :: sz,sxz
2113 ! *************************************************************************
2114 
2115  sx=zero ; sx2=zero ; msrt=zero
2116  sz=czero ; sxz=czero
2117  do ii=1,nn
2118   sx=sx+xx(ii)
2119   sz=sz+zz(ii)
2120   sxz=sxz+xx(ii)*zz(ii)
2121   sx2=sx2+xx(ii)*xx(ii)
2122  end do
2123 
2124  aa=CMPLX((nn*sxz-sx*sz)/(nn*sx2-sx*sx), kind=spc)
2125  bb=CMPLX(sz/nn-sx*aa/nn, kind=spc)
2126 
2127  do ii=1,nn
2128   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2129  end do
2130  msrt=SQRT(msrt) ; res=msrt
2131 
2132 end function linfit_spc

m_numeric_tools/linspace [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linspace

FUNCTION

INPUTS

OUTPUT

SOURCE

394 pure function linspace(start, stop, nn)
395 
396 !Arguments ------------------------------------
397 !scalars
398  integer,intent(in) :: nn
399  real(dp),intent(in) :: start, stop
400  real(dp) :: linspace(nn)
401 
402 !Local variables-------------------------------
403  real(dp) :: length
404  integer :: ii
405 ! *********************************************************************
406 
407  select case (nn)
408  case (1:)
409    length = stop - start
410    do ii=1,nn
411      linspace(ii) = start + length * (ii-1) / (nn-1)
412    end do
413 
414  case (0)
415    return
416  end select
417 
418 end function linspace

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

SOURCE

1891 subroutine list2blocks(list,nblocks,blocks)
1892 
1893 !Arguments ------------------------------------
1894 !scalars
1895  integer,intent(out) :: nblocks
1896  integer,intent(in) :: list(:)
1897 !arrays
1898  integer,intent(out),allocatable :: blocks(:,:)
1899 
1900 !Local variables-------------------------------
1901 !scalars
1902  integer :: ii,nitems
1903 !arrays
1904  integer :: work(2,size(list))
1905 
1906 !************************************************************************
1907 
1908  nitems = size(list)
1909 
1910  ! Handle nitems == 1 case
1911  if (nitems == 1) then
1912    ABI_MALLOC(blocks, (2,1))
1913    blocks = 1
1914    return
1915  end if
1916 
1917  nblocks = 1; work(1,1) = 1
1918 
1919  do ii=2,nitems
1920    if (list(ii) /= (list(ii-1) + 1)) then
1921      work(2,nblocks) = ii - 1
1922      nblocks = nblocks + 1
1923      work(1,nblocks) = ii
1924    end if
1925  end do
1926 
1927  work(2,nblocks) = nitems
1928 
1929  ABI_MALLOC(blocks, (2,nblocks))
1930  blocks = work(:,1:nblocks)
1931 
1932 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

SOURCE

2208 subroutine llsfit_svd(xx,yy,sigma,nfuncs,funcs,chisq,par,var,cov,info)
2209 
2210 !Arguments ------------------------------------
2211 !scalars
2212  integer,intent(in) :: nfuncs
2213  integer,intent(out) :: info
2214  real(dp),intent(out) :: chisq
2215 !arrays
2216  real(dp),intent(in) :: xx(:),yy(:),sigma(:)
2217  real(dp),intent(out) :: par(:),var(:),cov(:,:)
2218 
2219  interface
2220   function funcs(xx,nf)
2221   use defs_basis
2222   implicit none
2223   real(dp),intent(in) :: xx
2224   integer,intent(in) :: nf
2225   real(dp) :: funcs(nf)
2226   end function funcs
2227  end interface
2228 
2229 !Local variables-------------------------------
2230  integer,parameter :: PAD_=50
2231  integer :: ii,npts,lwork
2232  real(dp),parameter :: TOL_=1.0e-5_dp
2233 !arrays
2234  real(dp),dimension(SIZE(xx)) :: bb,sigm1
2235  real(dp),dimension(SIZE(xx),nfuncs) :: dmat,dmat_save
2236  real(dp) :: tmp(nfuncs)
2237  real(dp),allocatable :: work(:),Vt(:,:),U(:,:),S(:)
2238 ! *************************************************************************
2239 
2240  npts = assert_eq(SIZE(xx),SIZE(yy),SIZE(sigma),'Wrong size in xx,yy,sigma', __FILE__, __LINE__)
2241  call assert((npts>=nfuncs),'No. of functions must greater than no. of points', __FILE__, __LINE__)
2242  ii = assert_eq(nfuncs,SIZE(cov,1),SIZE(cov,2),SIZE(var),'Wrong size in covariance', __FILE__, __LINE__)
2243 
2244  !
2245  ! === Calculate design matrix and b vector ===
2246  ! * dmat_ij=f_j(x_i)/sigma_i, b_i=y_i/sigma_i
2247  sigm1(:)=one/sigma(:) ; bb(:)=yy(:)*sigm1(:)
2248  do ii=1,npts
2249   dmat_save(ii,:)=funcs(xx(ii),nfuncs)
2250  end do
2251  dmat=dmat_save*SPREAD(sigm1,DIM=2,ncopies=nfuncs)
2252  dmat_save(:,:)=dmat(:,:)
2253  !
2254  ! === Singular value decomposition ===
2255  lwork=MAX(3*MIN(npts,nfuncs)+MAX(npts,nfuncs),5*MIN(npts,nfuncs)-4)+PAD_
2256  ABI_MALLOC(work,(lwork))
2257  ABI_MALLOC(U,(npts,npts))
2258  ABI_MALLOC(S,(nfuncs))
2259  ABI_MALLOC(Vt,(nfuncs,nfuncs))
2260 
2261  call DGESVD('A','A',npts,nfuncs,dmat,npts,S,U,npts,Vt,nfuncs,work,lwork,info)
2262  ABI_FREE(work)
2263  GOTO 10
2264  !
2265  ! === Set to zero small singular values according to TOL_ and find coefficients ===
2266  WHERE (S>TOL_*MAXVAL(S))
2267   tmp=MATMUL(bb,U)/S
2268  ELSEWHERE
2269   S  =zero
2270   tmp=zero
2271  END WHERE
2272  par(:)=MATMUL(tmp,Vt)
2273  !
2274  ! === Evaluate chi-square ===
2275  chisq=l2norm(MATMUL(dmat_save,par)-bb)**2
2276  !
2277  ! === Calculate covariance and variance ===
2278  ! C_jk = V_ji V_ki / S_i^2
2279  WHERE (S/=zero) S=one/(S*S)
2280 
2281  ! check this but should be correct
2282  cov(:,:)=Vt*SPREAD(S,DIM=2,ncopies=nfuncs)
2283  cov(:,:)=MATMUL(TRANSPOSE(Vt),cov)
2284  var(:)=SQRT(get_diag(cov))
2285 
2286 10 continue
2287  ABI_FREE(U)
2288  ABI_FREE(S)
2289  ABI_FREE(Vt)
2290 
2291 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

SOURCE

1958 subroutine mask2blocks(mask,nblocks,blocks)
1959 
1960 !Arguments ------------------------------------
1961 !scalars
1962  integer,intent(out) :: nblocks
1963  logical,intent(in) :: mask(:)
1964 !arrays
1965  integer,allocatable :: blocks(:,:)
1966 
1967 !Local variables-------------------------------
1968 !scalars
1969  integer :: ii,nitems,start
1970  logical :: inblock
1971 !arrays
1972  integer :: work(2,SIZE(mask))
1973 
1974 !************************************************************************
1975 
1976  ! Find first element.
1977  nitems = size(mask); start = 0
1978  do ii=1,nitems
1979    if (mask(ii)) then
1980      start = ii
1981      exit
1982    end if
1983  end do
1984 
1985  ! Handle no true element or just one.
1986  if (start == 0) then
1987    nblocks = 0
1988    ABI_MALLOC(blocks, (0,0))
1989    return
1990  end if
1991  if (start /= 0 .and. nitems == 1) then
1992    nblocks = 1
1993    ABI_MALLOC(blocks, (2,1))
1994    blocks(:,1) = [1,1]
1995  end if
1996 
1997  nblocks = 1; work(1,1) = start; inblock = .True.
1998 
1999  do ii=start+1,nitems
2000    if (.not.mask(ii)) then
2001      if (inblock) then
2002        inblock = .False.
2003        work(2,nblocks) = ii - 1
2004      end if
2005    else
2006      if (.not. inblock) then
2007        inblock = .True.
2008        nblocks = nblocks + 1
2009        work(1,nblocks) = ii
2010      end if
2011    end if
2012  end do
2013 
2014  if (mask(nitems) .and. inblock) work(2,nblocks) = nitems
2015 
2016  ABI_MALLOC(blocks, (2,nblocks))
2017  blocks = work(:,1:nblocks)
2018 
2019 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.

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

2493  recursive subroutine midpoint_(func,nn,xmin,xmax,quad)
2494 
2495 !Arguments ------------------------------------
2496 !scalars
2497  integer,intent(in) :: nn
2498  !real(dp),external :: func
2499  real(dp),intent(in) :: xmin,xmax
2500  real(dp),intent(inout) :: quad
2501 
2502  interface
2503    function func(x)
2504      use defs_basis
2505      real(dp),intent(in) :: x
2506      real(dp) :: func
2507    end function func
2508  end interface
2509 
2510  !interface
2511  ! function func(x)
2512  !  use defs_basis
2513  !  real(dp),intent(in) :: x(:)
2514  !  real(dp) :: func(SIZE(x))
2515  ! end function func
2516  !end interface
2517 
2518 !Local variables-------------------------------
2519 !scalars
2520  integer  :: npt,ix
2521  real(dp) :: space
2522  character(len=500) :: msg
2523 !arrays
2524  real(dp),allocatable :: xx(:)
2525 
2526 !************************************************************************
2527 
2528  select case (nn)
2529 
2530  case (1)
2531    ! === Initial crude estimate done at the middle of the interval
2532    !quad=(xmax-xmin)*SUM(func((/half*(xmin+xmax)/))) !PARALLEL version
2533    quad=(xmax-xmin)*func(half*(xmin+xmax))
2534 
2535  case (2:)
2536    ! === Add npt interior points, they alternate in spacing between space and 2*space ===
2537    ABI_MALLOC(xx,(2*3**(nn-2)))
2538    npt=3**(nn-2) ; space=(xmax-xmin)/(three*npt)
2539    xx(1:2*npt-1:2)=arth(xmin+half*space,three*space,npt)
2540    xx(2:2*npt:2)=xx(1:2*npt-1:2)+two*space
2541    ! === The new sum is combined with the old integral to give a refined integral ===
2542    !quad=quad/three+space*SUM(func(xx))  !PARALLEL version
2543    quad=quad/three
2544    do ix=1,SIZE(xx)
2545      quad=quad+space*func(xx(ix))
2546    end do
2547    ABI_FREE(xx)
2548 
2549  case (:0)
2550    write(msg,'(a,i3)')' wrong value for nn ',nn
2551    ABI_BUG('Wrong value for nn')
2552  end select
2553 
2554 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

4422 integer function mincm(ii,jj)
4423 
4424 !Arguments ------------------------------------
4425 !scalars
4426  integer,intent(in) :: ii,jj
4427 
4428 !************************************************************************
4429 
4430  if (ii==0.or.jj==0) then
4431    ABI_BUG('ii==0 or jj==0')
4432  end if
4433 
4434  mincm=MAX(ii,jj)
4435  do
4436    if ( ((mincm/ii)*ii)==mincm .and. ((mincm/jj)*jj)==mincm ) RETURN
4437    mincm=mincm+1
4438  end do
4439 
4440 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

SOURCE

3400 pure subroutine mkherm(array,ndim)
3401 
3402 !Arguments -------------------------------
3403 !scalars
3404  integer,intent(in) :: ndim
3405 !arrays
3406  real(dp),intent(inout) :: array(2,ndim,ndim)
3407 
3408 !Local variables -------------------------
3409 !scalars
3410  integer :: i1,i2
3411 
3412 ! *********************************************************************
3413 
3414  do i1=1,ndim
3415    do i2=1,i1
3416      array(1,i1,i2)=(array(1,i1,i2)+array(1,i2,i1))*half
3417      array(2,i1,i2)=(array(2,i1,i2)-array(2,i2,i1))*half
3418      array(1,i2,i1)=array(1,i1,i2)
3419      array(2,i2,i1)=-array(2,i1,i2)
3420    end do
3421  end do
3422 
3423 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

SOURCE

5509 subroutine nderiv(hh, yy, zz, ndim, norder)
5510 
5511 !Arguments ---------------------------------------------
5512 !scalars
5513  integer,intent(in) :: ndim,norder
5514  real(dp),intent(in) :: hh
5515 !arrays
5516  real(dp),intent(in) :: yy(ndim)
5517  real(dp),intent(out) :: zz(ndim)
5518 
5519 !Local variables ---------------------------------------
5520 !scalars
5521  integer :: ier,ii
5522  real(dp) :: aa,bb,cc,h1,y1
5523 
5524 ! *********************************************************************
5525 
5526 !Initialization (common to 1st and 2nd derivative)
5527  h1=one/(12.d0*hh)
5528  y1=yy(ndim-4)
5529 
5530 !FIRST DERIVATIVE
5531 !================
5532  if (norder==1) then
5533 
5534 !  Prepare differentiation loop
5535    bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
5536    cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
5537 !  Start differentiation loop
5538    do ii=5,ndim
5539      aa=bb;bb=cc
5540      cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
5541      zz(ii-4)=aa
5542    end do
5543 !  Normal exit
5544    ier=0
5545    aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
5546    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))
5547    zz(ndim-1)=aa
5548    zz(ndim-2)=cc
5549    zz(ndim-3)=bb
5550 
5551 !  SECOND DERIVATIVE
5552 !  =================
5553  else
5554    h1=h1/hh
5555 !  Prepare differentiation loop
5556    bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
5557    cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
5558 !  Start differentiation loop
5559    do ii=5,ndim
5560      aa=bb;bb=cc
5561      cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
5562      zz(ii-4)=aa
5563    end do
5564 !  Normal exit
5565    ier=0
5566    aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
5567    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))
5568    zz(ndim-1)=aa
5569    zz(ndim-2)=cc
5570    zz(ndim-3)=bb
5571 
5572  end if !norder
5573 
5574 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

4176 complex(dp) function newrap_step(z, f, df)
4177 
4178 !Arguments ------------------------------------
4179 !scalars
4180  complex(dpc),intent(in) :: z,f,df
4181 
4182 !Local variables-------------------------------
4183  real(dp) :: dfm2
4184 ! *************************************************************************
4185 
4186  dfm2=ABS(df)*ABS(df)
4187 
4188  newrap_step = z - (f*CONJG(df))/dfm2
4189  !& z-one/(ABS(df)*ABS(df)) * CMPLX( REAL(f)*REAL(df)+AIMAG(f)*AIMAG(df), -REAL(f)*AIMAG(df)+AIMAG(f)*EAL(df) )
4190 
4191 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: 2 if matrix is complex, 1 for real matrix.
 mat_in(cplx, N*N)= matrix to be packed

OUTPUT

 mat_out(cplx*N*N+1/2)= packed matrix (upper triangle)

SOURCE

3702 subroutine pack_matrix(mat_in, mat_out, N, cplx)
3703 
3704 !Arguments ------------------------------------
3705 !scalars
3706  integer, intent(in) :: N, cplx
3707  real(dp), intent(in) :: mat_in(cplx, N*N)
3708  real(dp), intent(out) :: mat_out(cplx*N*(N+1)/2)
3709 
3710 !Local variables-------------------------------
3711  integer :: isubh, i, j
3712 
3713 ! *************************************************************************
3714 
3715  isubh = 1
3716  do j=1,N
3717    do i=1,j
3718      mat_out(isubh)    = mat_in(1, (j-1)*N+i)
3719      ! bad for vectorization, but it's not performance critical, so ...
3720      if(cplx == 2) then
3721        mat_out(isubh+1)  = mat_in(2, (j-1)*N+i)
3722      end if
3723      isubh=isubh+cplx
3724    end do
3725  end do
3726 
3727 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

4026 function pade(n, z, f, zz)
4027 
4028 !Arguments ------------------------------------
4029 !scalars
4030  integer,intent(in) :: n
4031  complex(dpc),intent(in) :: zz
4032  complex(dpc) :: pade
4033 !arrays
4034  complex(dpc),intent(in) :: z(n), f(n)
4035 
4036 !Local variables-------------------------------
4037 !scalars
4038  complex(dpc) :: a(n)
4039  complex(dpc) :: Az(0:n), Bz(0:n)
4040  integer :: i
4041 ! *************************************************************************
4042 
4043  call calculate_pade_a(a, n, z, f)
4044 
4045  Az(0)=czero
4046  Az(1)=a(1)
4047  Bz(0)=cone
4048  Bz(1)=cone
4049 
4050  do i=1,n-1
4051    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4052    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4053  end do
4054  !write(std_out,*) 'Bz(n)',Bz(n)
4055  !if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) ' Bz(n) ',Bz(n)
4056  pade=Az(n)/Bz(n)
4057  !write(std_out,*) 'pade_approx ', pade_approx
4058 
4059 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.

SOURCE

4662 subroutine pfactorize(nn,nfactors,pfactors,powers)
4663 
4664 !Arguments ------------------------------------
4665 !scalars
4666  integer,intent(in) :: nn,nfactors
4667  integer,intent(in) :: pfactors(nfactors)
4668  integer,intent(out) :: powers (nfactors+1)
4669 
4670 !Local variables ------------------------------
4671 !scalars
4672  integer :: tnn,ifc,fact,ipow,maxpwr
4673 
4674 ! *************************************************************************
4675 
4676  powers=0; tnn=nn
4677 
4678  fact_loop: do ifc=1,nfactors
4679    fact = pfactors (ifc)
4680    maxpwr = NINT ( LOG(DBLE(tnn))/LOG(DBLE(fact) ) ) + 1
4681    do ipow=1,maxpwr
4682      if (tnn==1) EXIT fact_loop
4683      if ( MOD(tnn,fact)==0 ) then
4684        tnn=tnn/fact
4685        powers(ifc)=powers(ifc) + 1
4686      end if
4687    end do
4688  end do fact_loop
4689 
4690  if ( nn /= tnn * PRODUCT( pfactors**powers(1:nfactors)) ) then
4691    ABI_BUG('nn/=tnn!')
4692  end if
4693 
4694  powers(nfactors+1) = tnn
4695 
4696 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

SOURCE

2318 subroutine polyn_interp(xa,ya,x,y,dy)
2319 
2320 !Arguments ------------------------------------
2321 !scalars
2322  real(dp),intent(in) :: xa(:),ya(:)
2323  real(dp),intent(in) :: x
2324  real(dp),intent(out) :: y,dy
2325 !Local variables-------------------------------
2326 !scalars
2327  integer :: m,n,ns
2328 !arrays
2329  real(dp),dimension(SIZE(xa)) :: c,d,den,ho
2330 ! *************************************************************************
2331 
2332  n = assert_eq(SIZE(xa),SIZE(ya),'Different size in xa and ya',__FILE__,__LINE__)
2333 
2334  ! === Initialize the tables of c and d ===
2335  c(:)=ya(:) ; d(:)=ya(:) ; ho(:)=xa(:)-x
2336  ! === Find closest table entry and initial approximation to y ===
2337  ns=imin_loc(ABS(x-xa)) ; y=ya(ns)
2338  ns=ns-1
2339  !
2340  ! === For each column of the tableau loop over current c and d and up-date them ===
2341  do m=1,n-1
2342   den(1:n-m)=ho(1:n-m)-ho(1+m:n)
2343   if (ANY(den(1:n-m)==zero)) then
2344    ABI_ERROR('Two input xa are identical')
2345   end if
2346 
2347   den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
2348   d(1:n-m)=ho(1+m:n)*den(1:n-m) ! Update c and d
2349   c(1:n-m)=ho(1:n-m)*den(1:n-m)
2350 
2351   if (2*ns<n-m) then  ! Now decide which correction, c or d, we want to add to the
2352    dy=c(ns+1)         ! accumulating value of y, The last dy added is the error indication.
2353   else
2354    dy=d(ns)
2355    ns=ns-1
2356   end if
2357 
2358   y=y+dy
2359  end do
2360 
2361 end subroutine polyn_interp

m_numeric_tools/print_arr1d_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr1d_dpc

FUNCTION

INPUTS

OUTPUT

SOURCE

3851 subroutine print_arr1d_dpc(arr, max_r, unit, mode_paral)
3852 
3853 !Arguments ------------------------------------
3854 !scalars
3855  integer,optional,intent(in) :: unit, max_r
3856  character(len=4),optional,intent(in) :: mode_paral
3857 !arrays
3858  complex(dpc),intent(in) :: arr(:)
3859 
3860 !Local variables-------------------------------
3861 !scalars
3862  integer :: unt,ii,nr,mr
3863  character(len=4) :: mode
3864  character(len=500) :: msg
3865  character(len=100) :: fmth,fmt1
3866 ! *************************************************************************
3867 
3868  unt=std_out ; if (PRESENT(unit      )) unt=unit
3869  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3870  mr=15       ; if (PRESENT(max_r     )) mr=max_r
3871 
3872  if (mode/='COLL'.and.mode/='PERS') then
3873   write(msg,'(2a)')' Wrong value of mode_paral ',mode
3874   ABI_BUG(msg)
3875  end if
3876 
3877  ! Print out matrix.
3878  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
3879 
3880  write(fmth,*)'(6x,',mr,'(i2,6x))'
3881  write(fmt1,*)'(3x,',mr,'f8.3)'
3882 
3883  write(msg,fmth)(ii,ii=1,mr)
3884  call wrtout(unt,msg,mode) ! header
3885  write(msg,fmt1)REAL (arr(1:mr))
3886  call wrtout(unt,msg,mode) !real part
3887  write(msg,fmt1)AIMAG(arr(1:mr))
3888  call wrtout(unt,msg,mode) !imag part
3889 
3890 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)

SOURCE

3795 subroutine print_arr1d_spc(arr,max_r,unit,mode_paral)
3796 
3797 !Arguments ------------------------------------
3798 !scalars
3799  integer,optional,intent(in) :: unit,max_r
3800  character(len=4),optional,intent(in) :: mode_paral
3801 !arrays
3802  complex(spc),intent(in) :: arr(:)
3803 
3804 !Local variables-------------------------------
3805 !scalars
3806  integer :: unt,ii,nr,mr
3807  character(len=4) :: mode
3808  character(len=500) :: msg
3809  character(len=100) :: fmth,fmt1
3810 ! *************************************************************************
3811 
3812  unt=std_out ; if (PRESENT(unit      )) unt=unit
3813  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3814  mr=15       ; if (PRESENT(max_r     )) mr=max_r
3815 
3816  if (mode/='COLL'.and.mode/='PERS') then
3817   write(msg,'(2a)')' Wrong value of mode_paral ',mode
3818   ABI_BUG(msg)
3819  end if
3820  !
3821  ! === Print out matrix ===
3822  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
3823 
3824  write(fmth,*)'(6x,',mr,'(i2,6x))'
3825  write(fmt1,*)'(3x,',mr,'f8.3)'
3826 
3827  write(msg,fmth)(ii,ii=1,mr)
3828  call wrtout(unt,msg,mode) !header
3829  write(msg,fmt1)REAL (arr(1:mr))
3830  call wrtout(unt,msg,mode) !real part
3831  write(msg,fmt1)AIMAG(arr(1:mr))
3832  call wrtout(unt,msg,mode) !imag part
3833 
3834 end subroutine print_arr1d_spc

m_numeric_tools/print_arr2d_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr2d_dpc

FUNCTION

INPUTS

OUTPUT

SOURCE

3968 subroutine print_arr2d_dpc(arr,max_r,max_c,unit,mode_paral)
3969 
3970 !Arguments ------------------------------------
3971 !scalars
3972  integer,optional,intent(in) :: unit,max_r,max_c
3973  character(len=4),optional,intent(in) :: mode_paral
3974 !arrays
3975  complex(dpc),intent(in) :: arr(:,:)
3976 
3977 !Local variables-------------------------------
3978 !scalars
3979  integer :: unt,ii,jj,nc,nr,mc,mr
3980  character(len=4) :: mode
3981  character(len=500) :: msg
3982  character(len=100) :: fmth,fmt1,fmt2
3983 ! *************************************************************************
3984 
3985  unt =std_out; if (PRESENT(unit      )) unt =unit
3986  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3987  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
3988  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
3989 
3990  if (mode/='COLL'.and.mode/='PERS') then
3991    write(msg,'(2a)')'Wrong value of mode_paral ',mode
3992    ABI_BUG(msg)
3993  end if
3994  !
3995  ! === Print out matrix ===
3996  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3997  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
3998 
3999  write(fmth,*)'(6x,',mc,'(i2,6x))'
4000  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
4001  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
4002 
4003  write(msg,fmth)(jj,jj=1,mc)
4004  call wrtout(unt, msg, mode) ! header
4005  do ii=1,mr
4006    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
4007    call wrtout(unt,msg,mode) ! real part
4008    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
4009    call wrtout(unt,msg,mode) ! imag part
4010  end do
4011 
4012 end subroutine print_arr2d_dpc

m_numeric_tools/print_arr2d_spc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr2d_spc

FUNCTION

INPUTS

OUTPUT

SOURCE

3907 subroutine print_arr2d_spc(arr, max_r, max_c, unit, mode_paral)
3908 
3909 !Arguments ------------------------------------
3910 !scalars
3911  integer,optional,intent(in) :: unit, max_r, max_c
3912  character(len=4),optional,intent(in) :: mode_paral
3913 !arrays
3914  complex(spc),intent(in) :: arr(:,:)
3915 
3916 !Local variables-------------------------------
3917 !scalars
3918  integer :: unt,ii,jj,nc,nr,mc,mr
3919  character(len=4) :: mode
3920  character(len=500) :: msg
3921  character(len=100) :: fmth,fmt1,fmt2
3922 ! *************************************************************************
3923 
3924  unt =std_out; if (PRESENT(unit      )) unt =unit
3925  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3926  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
3927  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
3928 
3929  if (mode/='COLL'.and.mode/='PERS') then
3930    write(msg,'(2a)')'Wrong value of mode_paral ',mode
3931    ABI_BUG(msg)
3932  end if
3933  !
3934  ! === Print out matrix ===
3935  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3936  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
3937 
3938  write(fmth,*)'(6x,',mc,'(i2,6x))'
3939  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
3940  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
3941 
3942  write(msg,fmth)(jj,jj=1,mc)
3943  call wrtout(unt,msg,mode) !header
3944  do ii=1,mr
3945    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
3946    call wrtout(unt,msg,mode) !real part
3947    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
3948    call wrtout(unt,msg,mode) !imag part
3949  end do
3950 
3951 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 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.

SOURCE

2588 recursive subroutine quadrature(func,xmin,xmax,qopt,quad,ierr,ntrial,accuracy,npts)
2589 
2590 !Arguments ------------------------------------
2591 !scalars
2592  integer,intent(in) :: qopt
2593  integer,intent(out) :: ierr
2594  integer,optional,intent(in) :: ntrial,npts
2595  real(dp),intent(in) :: xmin,xmax
2596  real(dp),optional,intent(in) :: accuracy
2597  real(dp),intent(out) :: quad
2598 
2599  interface
2600    function func(x)
2601      use defs_basis
2602      real(dp),intent(in) :: x
2603      real(dp) :: func
2604    end function func
2605  end interface
2606 
2607  !interface
2608  ! function func(x)
2609  !  use defs_basis
2610  !  real(dp),intent(in) :: x(:)
2611  !  real(dp) :: func(SIZE(x))
2612  ! end function func
2613  !end interface
2614 
2615 !Local variables-------------------------------
2616 !scalars
2617  integer :: K,KM,NT,NX,NX0,it,ix
2618  real(dp) :: EPS,old_st,st,old_quad,dqromb
2619  real(dp) :: TOL
2620  character(len=500) :: msg
2621 !arrays
2622  real(dp),allocatable :: h(:),s(:)
2623  real(dp),allocatable :: wx(:),xx(:)
2624 ! *************************************************************************
2625 
2626  ierr = 0
2627  TOL  =tol12
2628  EPS  =tol6  ; if (PRESENT(accuracy)) EPS=accuracy
2629  NT   =20    ; if (PRESENT(ntrial  )) NT=ntrial
2630  quad =zero
2631 
2632  select case (qopt)
2633 
2634  case (1)
2635    ! === Trapezoidal, closed form, O(1/N^2)
2636    do it=1,NT
2637      call trapezoidal_(func,it,xmin,xmax,quad)
2638      if (it>5) then ! Avoid spurious early convergence
2639        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2640      end if
2641      old_quad=quad
2642    end do
2643 
2644  case (2)
2645   ! === Extended Simpson rule based on trapezoidal O(1/N^4) ===
2646    do it=1,NT
2647      call trapezoidal_(func,it,xmin,xmax,st)
2648      if (it==1) then
2649        quad=st
2650      else
2651        quad=(four*st-old_st)/three
2652      end if
2653      if (it>5) then ! Avoid spurious early convergence
2654        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2655      end if
2656      old_quad=quad
2657      old_st=st
2658    end do
2659 
2660  case (3)
2661   ! === Midpoint rule, open form, O(1/N^2) ===
2662   do it=1,NT
2663     call midpoint_(func,it,xmin,xmax,quad)
2664     if (it>4) then ! Avoid spurious early convergence
2665       if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2666     end if
2667     old_quad=quad
2668   end do
2669 
2670  case (4)
2671    ! === Midpoint rule with cancellation of leading 1/N^2 term, open form, O(1/N^4) ===
2672    do it=1,NT
2673      call midpoint_(func,it,xmin,xmax,st)
2674      if (it==1) then
2675        quad=st
2676      else
2677        quad=(nine*st-old_st)/eight
2678      end if
2679      if (it>4) then ! Avoid spurious early convergence
2680        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2681      end if
2682      old_quad=quad
2683      old_st=st
2684    end do
2685 
2686  case (5)
2687    ! === Romberg Integration, closed form ===
2688    K=5 ; KM=K-1 ! Order 10
2689    ABI_MALLOC(h,(NT+1))
2690    ABI_MALLOC(s,(NT+1))
2691    h=zero
2692    s=zero
2693    h(1)=one
2694    do it=1,NT
2695      call trapezoidal_(func,it,xmin,xmax,s(it))
2696      !write(std_out,*) ' romberg-trap at ',ncall,it,s(it)
2697      if (it>=K) then
2698        call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
2699        if (ABS(dqromb)<EPS*ABS(quad)) then
2700          ABI_FREE(h)
2701          ABI_FREE(s)
2702          RETURN
2703        end if
2704      end if
2705      s(it+1)=s(it)
2706      h(it+1)=quarter*h(it) ! Quarter makes the extrapolation a polynomial in h^2,
2707    end do                 ! This is required to use the Euler-Maclaurin formula
2708    ABI_FREE(h)
2709    ABI_FREE(s)
2710 
2711  case (6)
2712    ! === Romberg Integration, closed form ===
2713    K=5 ; KM=K-1 ! Order 10
2714    ABI_MALLOC(h,(NT+1))
2715    ABI_MALLOC(s,(NT+1))
2716    h=zero
2717    s=zero
2718    h(1)=one
2719    do it=1,NT
2720      call midpoint_(func,it,xmin,xmax,s(it))
2721      if (it>=K) then
2722        call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
2723        !write(std_out,*) quad,dqromb
2724        if (ABS(dqromb)<EPS*ABS(quad)) then
2725          ABI_FREE(h)
2726          ABI_FREE(s)
2727          RETURN
2728        end if
2729      end if
2730      s(it+1)=s(it)
2731      h(it+1)=ninth*h(it) ! factor is due to step tripling in midpoint and even error series
2732    end do
2733    ABI_FREE(h)
2734    ABI_FREE(s)
2735 
2736  case (7)
2737    ! === Gauss-Legendre ===
2738    NX0=5 ; if (PRESENT(npts)) NX0=npts
2739    NX=NX0
2740    do it=1,NT
2741      ABI_MALLOC(wx,(NX))
2742      ABI_MALLOC(xx,(NX))
2743      call coeffs_gausslegint(xmin,xmax,xx,wx,NX)
2744      quad=zero
2745      do ix=1,NX
2746        quad=quad+wx(ix)*func(xx(ix))
2747      end do
2748      ABI_FREE(wx)
2749      ABI_FREE(xx)
2750      if (it>1) then
2751        !write(std_out,*) quad
2752        if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2753      end if
2754      old_quad=quad
2755      NX=NX+NX0
2756      !NX=2*NX
2757    end do
2758 
2759  case default
2760    write(msg,'(a,i3)')'Wrong value for qopt',qopt
2761    ABI_BUG(msg)
2762  end select
2763 
2764  write(msg,'(a,i0,2(a,es14.6))')&
2765   "Results are not converged within the given accuracy. ntrial= ",NT,"; EPS= ",EPS,"; TOL= ",TOL
2766  ABI_WARNING(msg)
2767  ierr = -1
2768 
2769 end subroutine quadrature

m_numeric_tools/rdp2cdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_2D

FUNCTION

INPUTS

OUTPUT

SOURCE

1037 pure function rdp2cdp_2D(rr) result(cc)
1038 
1039 !Arguments ------------------------------------
1040 !scalars
1041  real(dp),intent(in) :: rr(:,:,:)
1042  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3))
1043 
1044 ! *********************************************************************
1045 
1046  cc(:,:)=CMPLX(rr(1,:,:),rr(2,:,:), kind=dpc)
1047 
1048 end function rdp2cdp_2D

m_numeric_tools/rdp2cdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_3D

FUNCTION

INPUTS

OUTPUT

SOURCE

1065 pure function rdp2cdp_3D(rr) result(cc)
1066 
1067 !Arguments ------------------------------------
1068 !scalars
1069  real(dp),intent(in) :: rr(:,:,:,:)
1070  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4))
1071 
1072 ! *********************************************************************
1073 
1074  cc(:,:,:)=CMPLX(rr(1,:,:,:),rr(2,:,:,:), kind=dpc)
1075 
1076 end function rdp2cdp_3D

m_numeric_tools/rdp2cdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_4D

FUNCTION

INPUTS

OUTPUT

SOURCE

1093 pure function rdp2cdp_4D(rr) result(cc)
1094 
1095 !Arguments ------------------------------------
1096 !scalars
1097  real(dp),intent(in) :: rr(:,:,:,:,:)
1098  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5))
1099 
1100 ! *********************************************************************
1101 
1102  cc(:,:,:,:)=CMPLX(rr(1,:,:,:,:),rr(2,:,:,:,:), kind=dpc)
1103 
1104 end function rdp2cdp_4D

m_numeric_tools/rdp2cdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_5D

FUNCTION

INPUTS

OUTPUT

SOURCE

1121 pure function rdp2cdp_5D(rr) result(cc)
1122 
1123 !Arguments ------------------------------------
1124 !scalars
1125  real(dp),intent(in) :: rr(:,:,:,:,:,:)
1126  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6))
1127 
1128 ! *********************************************************************
1129 
1130  cc(:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:),rr(2,:,:,:,:,:), kind=dpc)
1131 
1132 end function rdp2cdp_5D

m_numeric_tools/rdp2cdp_6D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_6D

FUNCTION

INPUTS

OUTPUT

SOURCE

1149 pure function rdp2cdp_6D(rr) result(cc)
1150 
1151 !Arguments ------------------------------------
1152 !scalars
1153  real(dp),intent(in) :: rr(:,:,:,:,:,:,:)
1154  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6),SIZE(rr,7))
1155 
1156 ! *********************************************************************
1157 
1158  cc(:,:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:,:),rr(2,:,:,:,:,:,:), kind=dpc)
1159 
1160 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 first 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.

SOURCE

4291 subroutine remove_copies(n_in, set_in, n_out, is_equal)
4292 
4293 !Arguments ------------------------------------
4294 !scalars
4295  integer,intent(in) :: n_in
4296  integer,intent(out) :: n_out
4297 !arrays
4298  real(dp),target,intent(inout) :: set_in(3,n_in)
4299 
4300  interface
4301    function is_equal(k1,k2)
4302      use defs_basis
4303      real(dp),intent(in) :: k1(3),k2(3)
4304      logical :: is_equal
4305    end function is_equal
4306  end interface
4307 
4308 !Local variables-------------------------------
4309 !scalars
4310  integer :: ii,jj
4311  logical :: isnew
4312 !arrays
4313  type rdp1d_pt
4314   integer :: idx
4315   real(dp),pointer :: rpt(:)
4316  end type rdp1d_pt
4317  type(rdp1d_pt),allocatable :: Ap(:)
4318 
4319 ! *************************************************************************
4320 
4321  ABI_MALLOC(Ap,(n_in))
4322  Ap(1)%idx = 1
4323  Ap(1)%rpt => set_in(:,1)
4324 
4325  n_out=1
4326  do ii=2,n_in
4327 
4328    isnew=.TRUE.
4329    do jj=1,n_out
4330      if (is_equal(set_in(:,ii),Ap(jj)%rpt(:))) then
4331        isnew=.FALSE.
4332        exit
4333      end if
4334    end do
4335 
4336    if (isnew) then
4337      n_out=n_out+1
4338      Ap(n_out)%rpt => set_in(:,ii)
4339      Ap(n_out)%idx = ii
4340    end if
4341  end do
4342 
4343  ! The n_out inequivalent items are packed first.
4344  if (n_out/=n_in) then
4345    do ii=1,n_out
4346      jj=Ap(ii)%idx
4347      set_in(:,ii) = set_in(:,jj)
4348      !write(std_out,*) Ap(ii)%idx,Ap(ii)%rpt(:)
4349    end do
4350  end if
4351 
4352  ABI_FREE(Ap)
4353 
4354 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

SOURCE

474 subroutine reverse_int(arr)
475 
476 !Arguments ------------------------------------
477 !scalars
478  integer,intent(inout) :: arr(:)
479 !arrays
480  integer :: ii,nn,swap
481 ! *************************************************************************
482 
483  nn = SIZE(arr)
484  if (nn <= 1) return
485 
486  do ii=1,nn/2
487    swap = arr(ii)
488    arr(ii) = arr(nn-ii+1)
489    arr(nn-ii+1) = swap
490  end do
491 
492 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

SOURCE

506 subroutine reverse_rdp(arr)
507 
508 !Arguments ------------------------------------
509 !scalars
510  real(dp),intent(inout) :: arr(:)
511 !arrays
512  integer :: ii,nn
513  real(dp) :: swap
514 ! *************************************************************************
515 
516  nn = SIZE(arr)
517  if (nn <= 1) return
518 
519  do ii=1,nn/2
520    swap = arr(ii)
521    arr(ii) = arr(nn-ii+1)
522    arr(nn-ii+1) = swap
523  end do
524 
525 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 = module of cx

SOURCE

5279 pure subroutine rhophi(cx, phi, rho)
5280 
5281 !Arguments ------------------------------------
5282 !scalars
5283  real(dp),intent(out) :: phi,rho
5284 !arrays
5285  real(dp),intent(in) :: cx(2)
5286 
5287 ! ***********************************************************************
5288 
5289  rho = sqrt(cx(1)*cx(1) + cx(2)*cx(2))
5290 
5291  if (abs(cx(1)) > tol8) then
5292 
5293    phi = atan(cx(2)/cx(1))
5294 
5295    ! phi is an element of [-pi,pi]
5296    if (cx(1) < zero) then
5297      if (phi < zero) then
5298        phi = phi + pi
5299      else
5300        phi = phi - pi
5301      end if
5302    end if
5303 
5304  else
5305 
5306    if (cx(2) > tol8) then
5307      phi = pi*half
5308    else if (cx(2) < tol8) then
5309      phi = -pi*half
5310    else
5311      phi = zero
5312    end if
5313 
5314  end if
5315 
5316 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.

SOURCE

5239 function simpson(step, values) result(res)
5240 
5241 !Arguments ------------------------------------
5242 !scalars
5243  real(dp),intent(in) :: step
5244  real(dp) :: res
5245 !arrays
5246  real(dp),intent(in) :: values(:)
5247 
5248 !Local variables -------------------------
5249 !scalars
5250  real(dp) :: int_values(size(values))
5251 
5252 ! *********************************************************************
5253 
5254  call simpson_int(size(values),step,values,int_values)
5255  res = int_values(size(values))
5256 
5257 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.

SOURCE

3170 function simpson_cplx(npts,step,ff)
3171 
3172 !Arguments ------------------------------------
3173 !scalars
3174  integer,intent(in) :: npts
3175  real(dp),intent(in)  :: step
3176  complex(dpc),intent(in) :: ff(npts)
3177  complex(dpc) :: simpson_cplx
3178 
3179 !Local variables ------------------------------
3180 !scalars
3181  integer :: ii,my_n
3182  complex(dpc) :: sum_even, sum_odd
3183 
3184 !************************************************************************
3185 
3186  my_n=npts; if ((npts/2)*2 == npts) my_n=npts-3
3187 
3188  if (my_n<2) then
3189    ABI_ERROR("Too few points")
3190  end if
3191 
3192  sum_odd=czero
3193  do ii=2,my_n-1,2
3194    sum_odd = sum_odd + ff(ii)
3195  end do
3196 
3197  sum_even=zero
3198  do ii=3,my_n-2,2
3199    sum_even = sum_even + ff(ii)
3200  end do
3201 
3202  ! Eq 25.4.6 Abramowitz. Error is O(step^4)
3203  simpson_cplx = step/three * (ff(1) + four*sum_odd + two*sum_even + ff(my_n))
3204 
3205  if (my_n/=npts) then ! Simpson's 3/8 rule. Eq 25.4.13 Abramowitz. Error is O(step^5)
3206   simpson_cplx = simpson_cplx + three*step/eight * (ff(npts-3) + 3*ff(npts-2) + 3*ff(npts-1) + ff(npts))
3207  end if
3208 
3209 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.

SOURCE

5173 subroutine simpson_int(npts, step, values, int_values)
5174 
5175 !Arguments ------------------------------------
5176 !scalars
5177  integer,intent(in) :: npts
5178  real(dp),intent(in) :: step
5179 !arrays
5180  real(dp),intent(in) :: values(npts)
5181  real(dp),intent(out) :: int_values(npts)
5182 
5183 !Local variables -------------------------
5184 !scalars
5185  integer :: ii
5186  real(dp),parameter :: coef1 = 0.375_dp                          !9.0_dp  / 24.0_dp
5187  real(dp),parameter :: coef2 = 1.166666666666666666666666667_dp  !28.0_dp / 24.0_dp
5188  real(dp),parameter :: coef3 = 0.958333333333333333333333333_dp  !23.0_dp / 24.0_dp
5189  character(len=500) :: msg
5190 
5191 ! *********************************************************************
5192 
5193  if (npts < 6) then
5194    write(msg,"(a,i0)")"Number of points in integrand function must be >=6 while it is: ",npts
5195    ABI_ERROR(msg)
5196  end if
5197 
5198 !-----------------------------------------------------------------
5199 !Simpson integral of input function
5200 !-----------------------------------------------------------------
5201 
5202 !first point is 0: don t store it
5203 !do integration equivalent to Simpson O(1/N^4) from NumRec in C p 134  NumRec in Fortran p 128
5204  int_values(1) =               coef1*values(1)
5205  int_values(2) = int_values(1) + coef2*values(2)
5206  int_values(3) = int_values(2) + coef3*values(3)
5207 
5208  do ii=4,npts-3
5209    int_values(ii) = int_values(ii-1) + values(ii)
5210  end do
5211 
5212  int_values(npts-2) = int_values(npts-3) + coef3*values(npts-2)
5213  int_values(npts-1) = int_values(npts-2) + coef2*values(npts-1)
5214  int_values(npts  ) = int_values(npts-1) + coef1*values(npts  )
5215 
5216  int_values(:) = int_values(:) * step
5217 
5218 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. <= 0 to return a unchanged.

SIDE EFFECTS

  a(mesh)=Input values, smoothed in output

SOURCE

5448 subroutine smooth(a, mesh, it)
5449 
5450 !Arguments ------------------------------------
5451 !scalars
5452  integer, intent(in) :: it,mesh
5453  real(dp), intent(inout) :: a(mesh)
5454 !Local variables-------------------------------
5455  integer :: i,k
5456  real(dp) :: asm(mesh)
5457 ! *********************************************************************
5458 
5459  do k=1,it
5460    asm(1)=1.0d0/3.0d0*(a(1)+a(2)+a(3))
5461    asm(2)=0.25d0*(a(1)+a(2)+a(3)+a(4))
5462    asm(3)=0.2d0*(a(1)+a(2)+a(3)+a(4)+a(5))
5463    asm(4)=0.2d0*(a(2)+a(3)+a(4)+a(5)+a(6))
5464    asm(5)=0.2d0*(a(3)+a(4)+a(5)+a(6)+a(7))
5465    asm(mesh-4)=0.2d0*(a(mesh-2)+a(mesh-3)+a(mesh-4)+&
5466 &                   a(mesh-5)+a(mesh-6))
5467    asm(mesh-3)=0.2d0*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
5468 &                   a(mesh-4)+a(mesh-5))
5469    asm(mesh-2)=0.2d0*(a(mesh)+a(mesh-1)+a(mesh-2)+&
5470 &                   a(mesh-3)+a(mesh-4))
5471    asm(mesh-1)=0.25d0*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
5472    asm(mesh)=1.0d0/3.0d0*(a(mesh)+a(mesh-1)+a(mesh-2))
5473 
5474    do i=6,mesh-5
5475      asm(i)=0.1d0*a(i)+0.1d0*(a(i+1)+a(i-1))+&
5476 &             0.1d0*(a(i+2)+a(i-2))+&
5477 &             0.1d0*(a(i+3)+a(i-3))+&
5478 &             0.1d0*(a(i+4)+a(i-4))+&
5479 &             0.05d0*(a(i+5)+a(i-5))
5480    end do
5481 
5482    do i=1,mesh
5483      a(i)=asm(i)
5484    end do
5485  end do
5486 
5487 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 dataset.

 INPUT
  arr(:)=Array with the values.

OUTPUT

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

SOURCE

4780 pure function stats_eval(arr) result(stats)
4781 
4782 !Arguments ------------------------------------
4783 !scalars
4784  type(stats_t) :: stats
4785 !arrays
4786  real(dp),intent(in) :: arr(:)
4787 
4788 !Local variables ------------------------------
4789 !scalars
4790  integer :: ii,nn
4791  real(dp) :: xx,x2_sum
4792 
4793 ! *************************************************************************
4794 
4795  stats%min   = +HUGE(one)
4796  stats%max   = -HUGE(one)
4797  stats%mean  = zero
4798 
4799  nn = SIZE(arr)
4800  do ii=1,nn
4801    xx = arr(ii)
4802    stats%max  = MAX(stats%max, xx)
4803    stats%min  = MIN(stats%min, xx)
4804    stats%mean = stats%mean + xx
4805  end do
4806 
4807  stats%mean  = stats%mean/nn
4808 
4809  ! Two-pass algorithm for the variance (more stable than the single-pass one).
4810  x2_sum = zero
4811  do ii=1,nn
4812    xx     = arr(ii)
4813    x2_sum = x2_sum + (xx - stats%mean)*(xx - stats%mean)
4814  end do
4815 
4816  if (nn>1) then
4817    stats%stdev  = x2_sum/(nn-1)
4818    stats%stdev = SQRT(ABS(stats%stdev))
4819  else
4820    stats%stdev = zero
4821  end if
4822 
4823 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

254  type,public :: stats_t
255    real(dp) :: mean
256    real(dp) :: stdev
257    real(dp) :: min
258    real(dp) :: max
259  end type stats_t
260 
261  public :: stats_eval  ! Calculate 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

SOURCE

3635 subroutine symmetrize_dpc(mat, uplo)
3636 
3637 !Arguments ------------------------------------
3638 !scalars
3639  character(len=*),intent(in) :: uplo
3640 !arrays
3641  complex(dpc),intent(inout) :: mat(:,:)
3642 
3643 !Local variables-------------------------------
3644 !scalars
3645  integer :: nn,ii,jj
3646 !arrays
3647  complex(dpc),allocatable :: tmp(:)
3648 ! *************************************************************************
3649 
3650  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3651 
3652  select case (uplo(1:1))
3653  case ("A","a") ! Full matrix has been calculated.
3654    ABI_MALLOC(tmp,(nn))
3655    do ii=1,nn
3656      do jj=ii,nn
3657        tmp(jj)=half*(mat(ii,jj)+mat(jj,ii))
3658      end do
3659      mat(ii,ii:nn)=tmp(ii:nn)
3660      mat(ii:nn,ii)=tmp(ii:nn)
3661    end do
3662    ABI_FREE(tmp)
3663 
3664  case ("U","u") ! Only the upper triangle is used.
3665    do jj=1,nn
3666      do ii=1,jj-1
3667        mat(jj,ii) = mat(ii,jj)
3668      end do
3669    end do
3670 
3671  case ("L","l") ! Only the lower triangle is used.
3672   do jj=1,nn
3673     do ii=1,jj-1
3674       mat(ii,jj) = mat(jj,ii)
3675     end do
3676   end do
3677 
3678  case default
3679    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3680  end select
3681 
3682 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

SOURCE

3560 subroutine symmetrize_spc(mat,uplo)
3561 
3562 !Arguments ------------------------------------
3563 !scalars
3564  character(len=*),intent(in) :: uplo
3565 !arrays
3566  complex(spc),intent(inout) :: mat(:,:)
3567 
3568 !Local variables-------------------------------
3569 !scalars
3570  integer :: nn,ii,jj
3571 !arrays
3572  complex(spc),allocatable :: tmp(:)
3573 ! *************************************************************************
3574 
3575  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3576 
3577  select case (uplo(1:1))
3578 
3579  case ("A","a") ! Full matrix has been calculated.
3580    ABI_MALLOC(tmp,(nn))
3581    do ii=1,nn
3582      do jj=ii,nn
3583        tmp(jj)=REAL(half)*(mat(ii,jj)+mat(jj,ii))
3584      end do
3585      mat(ii,ii:nn)=tmp(ii:nn)
3586      mat(ii:nn,ii)=tmp(ii:nn)
3587    end do
3588    ABI_FREE(tmp)
3589 
3590  case ("U","u") ! Only the upper triangle is used.
3591    do jj=1,nn
3592      do ii=1,jj-1
3593        mat(jj,ii) = mat(ii,jj)
3594      end do
3595    end do
3596 
3597  case ("L","l") ! Only the lower triangle is used.
3598   do jj=1,nn
3599     do ii=1,jj-1
3600       mat(ii,jj) = mat(jj,ii)
3601     end do
3602   end do
3603 
3604  case default
3605    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3606  end select
3607 
3608 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.

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

2393 recursive subroutine trapezoidal_(func,nn,xmin,xmax,quad)
2394 
2395 !Arguments ------------------------------------
2396 !scalars
2397  integer,intent(in) :: nn
2398  !real(dp),external :: func
2399  real(dp),intent(in) :: xmin,xmax
2400  real(dp),intent(inout) :: quad
2401 
2402  interface
2403    function func(x)
2404      use defs_basis
2405      real(dp),intent(in) :: x
2406      real(dp) :: func
2407    end function func
2408  end interface
2409 
2410  !interface
2411  ! function func(x)
2412  !  use defs_basis
2413  !  real(dp),intent(in) :: x(:)
2414  !  real(dp) :: func(SIZE(x))
2415  ! end function func
2416  !end interface
2417 
2418 !Local variables-------------------------------
2419 !scalars
2420  integer :: npt,ix
2421  real(dp) :: space,new,yy
2422  character(len=500) :: msg
2423 !arrays
2424  !real(dp),allocatable :: xx(:)
2425 !************************************************************************
2426 
2427  select case (nn)
2428 
2429  case (1)
2430    ! === Initial crude estimate (xmax-xmin)(f1+f2)/2 ===
2431    !quad=half*(xmax-xmin)*SUM(func((/xmin,xmax/)))
2432    quad=half*(xmax-xmin)*(func(xmin)+func(xmax))
2433 
2434  case (2:)
2435    ! === Add npt interior points of spacing space ===
2436    npt=2**(nn-2) ; space=(xmax-xmin)/npt
2437    ! === The new sum is combined with the old integral to give a refined integral ===
2438    !new=SUM(func(arth(xmin+half*space,space,npt))) !PARALLEL version
2439    !allocate(xx(npt))
2440    !xx(:)=arth(xmin+half*space,space,npt)
2441    !xx(1)=xmin+half*space
2442    !do ii=2,nn
2443    ! xx(ii)=xx(ii-1)+space
2444    !end do
2445    new=zero
2446    yy=xmin+half*space
2447    do ix=1,npt
2448     !new=new+func(xx(ix))
2449     new=new+func(yy)
2450     yy=yy+space
2451    end do
2452    !deallocate(xx)
2453    quad=half*(quad+space*new)
2454    !write(std_out,*) 'trapezoidal',quad
2455 
2456  case (:0)
2457    write(msg,'(a,i3)')'Wrong value for nn ',nn
2458    ABI_BUG(msg)
2459  end select
2460 
2461 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

SOURCE

5686 function uniformrandom(seed)
5687 
5688 !Arguments ------------------------------------
5689 !scalars
5690  real(dp) :: uniformrandom
5691  integer,intent(inout) :: seed
5692 
5693 !Local variables ---------------------------------------
5694  integer, parameter :: im1=11979,ia1= 430,ic1=2531
5695  integer, parameter :: im2= 6655,ia2= 936,ic2=1399
5696  integer, parameter :: im3= 6075,ia3=1366,ic3=1283
5697  integer, save :: init=0
5698  integer, save :: ii1,ii2,ii3
5699  integer :: kk
5700  real(dp) :: im1inv,im2inv
5701  real(dp), save :: table(97)
5702  character(len=500) :: msg
5703 
5704 ! *********************************************************************
5705 
5706  im1inv=1.0d0/im1 ; im2inv=1.0d0/im2
5707 
5708 !Initialize on first call or when seed<0:
5709  if (seed<0.or.init==0) then
5710    seed=-abs(seed)
5711 
5712 !  First generator
5713    ii1=mod(ic1-seed,im1)
5714    ii1=mod(ia1*ii1+ic1,im1)
5715 !  Second generator
5716    ii2=mod(ii1,im2)
5717    ii1=mod(ia1*ii1+ic1,im1)
5718 !  Third generator
5719    ii3=mod(ii1,im3)
5720 
5721 !  Fill table
5722    do kk=1,97
5723      ii1=mod(ia1*ii1+ic1,im1)
5724      ii2=mod(ia2*ii2+ic2,im2)
5725      table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
5726    enddo
5727 
5728    init=1 ; seed=1
5729  end if
5730 
5731 !Third generator gives index
5732  ii3=mod(ia3*ii3+ic3,im3)
5733  kk=1+(97*ii3)/im3
5734  if (kk<1.or.kk>97) then
5735    write(msg,'(a,2i0,a)' ) ' trouble in uniformrandom; ii3,kk=',ii3,kk,' =>stop'
5736    ABI_ERROR(msg)
5737  end if
5738  uniformrandom=table(kk)
5739 
5740 !Replace old value, based on generators 1 and 2
5741  ii1=mod(ia1*ii1+ic1,im1)
5742  ii2=mod(ia2*ii2+ic2,im2)
5743  table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
5744 
5745 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

SOURCE

608 pure subroutine unit_matrix_cdp(matrix)
609 
610 !Arguments ------------------------------------
611  complex(dpc),intent(inout) :: matrix(:,:)
612 
613 !Local variables-------------------------------
614 !scalars
615  integer :: ii,nn
616 ! *********************************************************************
617 
618  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
619  matrix=czero
620  do ii=1,nn
621    matrix(ii,ii)=cone
622  end do
623 
624 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

SOURCE

542 pure subroutine unit_matrix_int(matrix)
543 
544 !Arguments ------------------------------------
545  integer,intent(inout) :: matrix(:,:)
546 
547 !Local variables-------------------------------
548 !scalars
549  integer :: ii,nn
550 ! *********************************************************************
551 
552  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
553  matrix(:,:)=0
554  do ii=1,nn
555   matrix(ii,ii)=1
556  end do
557 
558 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

SOURCE

575 pure subroutine unit_matrix_rdp(matrix)
576 
577 !Arguments ------------------------------------
578  real(dp),intent(inout) :: matrix(:,:)
579 
580 !Local variables-------------------------------
581 !scalars
582  integer :: ii,nn
583 ! *********************************************************************
584 
585  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
586  matrix(:,:)=zero
587  do ii=1,nn
588    matrix(ii,ii)=one
589  end do
590 
591 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
  [vd_max]= Compute max value of the different entries.

OUTPUT

  vdiff_t object

SOURCE

5340 type(vdiff_t) function vdiff_eval(cplex, nr, f1, f2, volume, vd_max) result(vd)
5341 
5342 !Arguments ------------------------------------
5343 !scalars
5344  integer,intent(in) :: cplex,nr
5345  real(dp),intent(in) :: volume
5346  type(vdiff_t),optional,intent(inout) :: vd_max
5347 !arrays
5348  real(dp),intent(in) :: f1(cplex,nr),f2(cplex,nr)
5349 
5350 !Local variables-------------------------------
5351 !scalars
5352  integer :: ir
5353  real(dp) :: num,den,dr
5354  type(stats_t) :: stats
5355 !arrays
5356  real(dp) :: abs_diff(nr)
5357 ! *********************************************************************
5358 
5359  dr = volume / nr
5360 
5361  if (cplex == 1) then
5362    abs_diff = abs(f1(1,:) - f2(1,:))
5363    num = sum(abs_diff)
5364    den = sum(abs(f2(1,:)))
5365 
5366  else if (cplex == 2) then
5367    do ir=1,nr
5368      abs_diff(ir) = sqrt((f1(1,ir) - f2(1,ir))**2 + (f1(2,ir) - f2(2,ir))**2)
5369    end do
5370    num = sum(abs_diff)
5371    den = zero
5372    do ir=1,nr
5373      den = den + sqrt(f2(1,ir)**2 + f2(2,ir)**2)
5374    end do
5375  end if
5376 
5377  vd%int_adiff = num * dr
5378  call safe_div(num,den,zero,vd%l1_rerr)
5379 
5380  stats = stats_eval(abs_diff)
5381  vd%mean_adiff = stats%mean
5382  vd%stdev_adiff = stats%stdev
5383  vd%min_adiff = stats%min
5384  vd%max_adiff = stats%max
5385 
5386  if (present(vd_max)) then
5387    vd_max%int_adiff =   max(vd_max%int_adiff, vd%int_adiff)
5388    vd_max%mean_adiff =  max(vd_max%mean_adiff, vd%mean_adiff)
5389    vd_max%stdev_adiff = max(vd_max%stdev_adiff, vd%stdev_adiff)
5390    vd_max%min_adiff =   max(vd_max%min_adiff, vd%min_adiff)
5391    vd_max%max_adiff =   max(vd_max%max_adiff, vd%max_adiff)
5392    vd_max%l1_rerr =     max(vd_max%l1_rerr, vd%L1_rerr)
5393  end if
5394 
5395 end function vdiff_eval

m_numeric_tools/vdiff_print [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 vdiff_print

FUNCTION

  Print vdiff_t to unit

SOURCE

5409 subroutine vdiff_print(vd, unit)
5410 
5411 !Arguments ------------------------------------
5412 !scalars
5413  integer,optional,intent(in) :: unit
5414  type(vdiff_t),intent(in) :: vd
5415 
5416 !Local variables-------------------------------
5417 !scalars
5418  integer :: unt
5419 ! *********************************************************************
5420 
5421  unt = std_out; if (present(unit)) unt = unit
5422  write(unt,"(a,es10.3,a)")"  L1_rerr: ", vd%l1_rerr, ","
5423  write(unt,"(a,es10.3,a)")"  'Integral |f1-f2|dr': ", vd%int_adiff, ","
5424  write(unt,"(a,es10.3,a)")"  'min {|f1-f2|}': ", vd%min_adiff, ","
5425  write(unt,"(a,es10.3,a)")"  'Max {|f1-f2|}': ", vd%max_adiff, ","
5426  write(unt,"(a,es10.3,a)")"  'mean {|f1-f2|}': ", vd%mean_adiff, ","
5427  write(unt,"(a,es10.3,a)")"  'stdev {|f1-f2|}': ", vd%stdev_adiff, ","
5428 
5429 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

276  type,public :: vdiff_t
277 
278    real(dp) :: int_adiff = zero     ! \int |f1-f2| dr
279    real(dp) :: mean_adiff = zero    ! Mean {|f1-f2|}
280    real(dp) :: stdev_adiff = zero   ! Standard deviation of {|f1-f2|}
281    real(dp) :: min_adiff = zero     ! Min {|f1-f2|}
282    real(dp) :: max_adiff = zero     ! Max {|f1-f2|}
283    real(dp) :: l1_rerr = zero       ! (\int |f1-f2| dr) / (\int |f2| dr)
284 
285  end type vdiff_t
286 
287  public :: vdiff_eval         ! Estimate the "distance" between two functions tabulated on a homogeneous grid.
288  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

SOURCE

4886 elemental subroutine wrap2_pmhalf(num,red,shift)
4887 
4888 !Arguments -------------------------------
4889 !scalars
4890  real(dp),intent(in) :: num
4891  real(dp),intent(out) :: red,shift
4892 
4893 ! *********************************************************************
4894 
4895  if (num>zero) then
4896    red=mod((num+half-tol12),one)-half+tol12
4897  else
4898    red=-mod(-(num-half-tol12),one)+half+tol12
4899  end if
4900  if(abs(red)<tol12)red=0.0d0
4901  shift=num-red
4902 
4903 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

SOURCE

4846 elemental subroutine wrap2_zero_one(num, red, shift)
4847 
4848 !Arguments ------------------------------------
4849 !scalars
4850  real(dp),intent(in) :: num
4851  real(dp),intent(out) :: red,shift
4852 
4853 ! *************************************************************************
4854 
4855  if (num>zero) then
4856    red=mod((num+tol12),one)-tol12
4857  else
4858    red=-mod(-(num-one+tol12),one)+one-tol12
4859  end if
4860  if(abs(red)<tol12)red=0.0_dp
4861  shift=num-red
4862 
4863 end subroutine wrap2_zero_one