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

6403 integer pure function blocked_loop(loop_index, loop_stop, batch_size) result(ndat)
6404 
6405 !Arguments ----------------------------------------------
6406  integer,intent(in) :: loop_index, loop_stop, batch_size
6407 
6408 ! *********************************************************************
6409 
6410  ndat = merge(batch_size, loop_stop - loop_index + 1, loop_index + batch_size - 1 <= loop_stop)
6411 
6412 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

6248 subroutine bool2index(bool_list, out_index)
6249 
6250 !Arguments ----------------------------------------------
6251 !scalars
6252  logical,intent(in) :: bool_list(:)
6253  integer,allocatable,intent(inout) :: out_index(:)
6254 
6255 !Local variables-------------------------------
6256  integer :: ii, cnt
6257 ! *********************************************************************
6258 
6259  cnt = count(bool_list)
6260  ABI_REMALLOC(out_index, (cnt))
6261  cnt = 0
6262  do ii=1,size(bool_list)
6263    if (bool_list(ii)) then
6264      cnt = cnt + 1
6265      out_index(cnt) = ii
6266    end if
6267  end do
6268 
6269 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

6098 function dotproduct(nv1,nv2,v1,v2)
6099 
6100 !Arguments ------------------------------------
6101 !scalars
6102  integer,intent(in) :: nv1,nv2
6103  real(dp) :: dotproduct
6104 !arrays
6105  real(dp),intent(in) :: v1(nv1,nv2),v2(nv1,nv2)
6106 
6107 !Local variables-------------------------------
6108 !scalars
6109  integer :: i,j
6110 
6111 ! *************************************************************************
6112  dotproduct=zero
6113  do j=1,nv2
6114   do i=1,nv1
6115    dotproduct=dotproduct+v1(i,j)*v2(i,j)
6116   end do
6117  end do
6118 
6119 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-2022 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  end interface c2r
177 
178  interface isinteger
179    module procedure is_integer_0d
180    module procedure is_integer_1d
181  end interface isinteger
182 
183  interface is_zero
184    module procedure is_zero_rdp_0d
185    module procedure is_zero_rdp_1d
186  end interface is_zero
187 
188  interface bisect
189    module procedure bisect_rdp
190    module procedure bisect_int
191  end interface bisect
192 
193  interface imax_loc
194    module procedure imax_loc_int
195    module procedure imax_loc_rdp
196  end interface imax_loc
197 
198  interface imin_loc
199    module procedure imin_loc_int
200    module procedure imin_loc_rdp
201  end interface imin_loc
202 
203  interface linfit
204    module procedure linfit_rdp
205    module procedure linfit_spc
206    module procedure linfit_dpc
207  end interface linfit
208 
209  interface hermitianize
210    module procedure hermitianize_spc
211    module procedure hermitianize_dpc
212  end interface hermitianize
213 
214  interface symmetrize
215    module procedure symmetrize_spc
216    module procedure symmetrize_dpc
217  end interface symmetrize
218 
219  interface print_arr  !TODO add prtm
220    module procedure print_arr1d_spc
221    module procedure print_arr1d_dpc
222    module procedure print_arr2d_spc
223    module procedure print_arr2d_dpc
224  end interface print_arr
225 
226  interface operator (.x.)
227    module procedure cross_product_int
228    module procedure cross_product_rdp
229  end interface
230 
231  interface l2norm
232    module procedure l2norm_rdp
233  end interface l2norm
234 
235  interface isordered
236    module procedure isordered_rdp
237  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

6300 subroutine polynomial_regression(degree,npts,xvals,yvals,coeffs,RMSerr)
6301 
6302 !Arguments ------------------------------------
6303 
6304 !scalars
6305  integer                     :: degree,npts
6306  real(dp),intent(out)        :: RMSerr
6307 !arrays
6308  real(dp),intent(in)         :: xvals(1:npts),yvals(1:npts)
6309  real(dp),intent(out)        :: coeffs(degree+1)
6310 
6311 !Local variables-------------------------------
6312 !scalars
6313  integer                     :: ncoeffs,icoeff,ipoint,info
6314  real(dp)                    :: residual,fitval
6315 !arrays
6316  integer,allocatable         :: tmp(:)
6317  real(dp),allocatable        :: tmptwo(:)
6318  real(dp),allocatable        :: A(:,:),ATA(:,:)
6319 
6320 !####################################################################
6321 !#####################  Get Polynomial Fit  #########################
6322 
6323   ncoeffs=degree+1
6324 
6325   ABI_MALLOC(tmp,(ncoeffs))
6326   ABI_MALLOC(tmptwo,(ncoeffs))
6327   ABI_MALLOC(A,(npts,ncoeffs))
6328   ABI_MALLOC(ATA,(ncoeffs,ncoeffs))
6329 
6330   !Construct a polynomial for all input xvalues
6331   do icoeff=1,ncoeffs
6332     do ipoint=1,npts
6333        if (icoeff==1.and.xvals(ipoint)==0.0) then
6334           A(ipoint,icoeff) = 1.0
6335        else
6336           A(ipoint,icoeff) = xvals(ipoint)**(icoeff-1)
6337        end if
6338     end do
6339   end do
6340 
6341   !Get matrix product of transpose of A and A
6342   ATA = matmul(transpose(A),A)
6343 
6344   !Compute LU factorization of ATA
6345   call DGETRF(ncoeffs,ncoeffs,ATA,ncoeffs,tmp,info)
6346   ABI_CHECK(info == 0, sjoin('LAPACK DGETRF in polynomial regression returned:', itoa(info)))
6347 
6348   !Compute inverse of the LU factorized version of ATA
6349   call DGETRI(ncoeffs,ATA,ncoeffs,tmp,tmptwo,ncoeffs,info)
6350   ABI_CHECK(info == 0, sjoin('LAPACK DGETRI in polynomial regression returned:', itoa(info)))
6351 
6352   !Harvest polynomial coefficients
6353   coeffs = matmul(matmul(ATA,transpose(A)),yvals)
6354 
6355 !####################################################################
6356 !##############  RMS error on the polynomial fit  ###################
6357 
6358   residual=0.0d0
6359   do ipoint=1,npts
6360     fitval=0.0d0
6361     do icoeff=1,ncoeffs
6362       if (icoeff==1.and.xvals(ipoint)==0.0) then
6363         fitval=fitval+coeffs(icoeff)
6364       else
6365         fitval=fitval+coeffs(icoeff)*xvals(ipoint)**(icoeff-1)
6366       end if
6367     end do
6368     residual=residual+(fitval-yvals(ipoint))**2
6369   end do
6370   RMSerr=sqrt(residual/(real(npts-1,8)))
6371 
6372 !####################################################################
6373 !########################  Deallocations  ###########################
6374 
6375   ABI_FREE(A)
6376   ABI_FREE(ATA)
6377   ABI_FREE(tmp)
6378   ABI_FREE(tmptwo)
6379 
6380 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

6221 elemental subroutine safe_div(n, d, altv, q)
6222 
6223 !Arguments ----------------------------------------------
6224 !scalars
6225  real(dp),intent(in) :: n, d, altv
6226  real(dp),intent(out) :: q
6227 
6228 ! *********************************************************************
6229 
6230  if ( exponent(n) - exponent(d) >= maxexponent(n) .or. d == zero) then
6231    q = altv
6232  else
6233    q = n / d
6234  endif
6235 
6236 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

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

m_numeric_tools/arth_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  arth_rdp

FUNCTION

INPUTS

OUTPUT

SOURCE

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

1613 pure function bisect_int(AA,xx) result(loc)
1614 
1615 !Arguments ------------------------------------
1616 !scalars
1617  integer,intent(in) :: AA(:)
1618  integer,intent(in) :: xx
1619  integer :: loc
1620 
1621 !Local variables-------------------------------
1622  integer :: nn,jl,jm,ju
1623  logical :: ascnd
1624 ! *********************************************************************
1625 
1626  nn=SIZE(AA) ; ascnd=(AA(nn)>=AA(1))
1627 
1628  ! Initialize lower and upper limits
1629  jl=0 ; ju=nn+1
1630  do
1631   if (ju-jl<=1) EXIT
1632   jm=(ju+jl)/2  ! Compute a midpoint
1633   if (ascnd.EQV.(xx>=AA(jm))) then
1634    jl=jm ! Replace lower limit
1635   else
1636    ju=jm ! Replace upper limit
1637   end if
1638  end do
1639  !
1640  ! Set the output, being careful with the endpoints
1641  if (xx==AA(1)) then
1642   loc=1
1643  else if (xx==AA(nn)) then
1644   loc=nn-1
1645  else
1646   loc=jl
1647  end if
1648 
1649 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

1558 pure function bisect_rdp(AA, xx) result(loc)
1559 
1560 !Arguments ------------------------------------
1561 !scalars
1562  real(dp),intent(in) :: AA(:)
1563  real(dp),intent(in) :: xx
1564  integer :: loc
1565 
1566 !Local variables-------------------------------
1567  integer :: nn,jl,jm,ju
1568  logical :: ascnd
1569 ! *********************************************************************
1570 
1571  nn=SIZE(AA); ascnd=(AA(nn)>=AA(1))
1572  !
1573  ! Initialize lower and upper limits
1574  jl=0; ju=nn+1
1575  do
1576    if (ju-jl<=1) EXIT
1577    jm=(ju+jl)/2  ! Compute a midpoint,
1578    if (ascnd.EQV.(xx>=AA(jm))) then
1579      jl=jm ! Replace lower limit
1580    else
1581      ju=jm ! Replace upper limit
1582    end if
1583  end do
1584  !
1585  ! Set the output, being careful with the endpoints
1586  if (xx==AA(1)) then
1587    loc=1
1588  else if (xx==AA(nn)) then
1589    loc=nn-1
1590  else
1591    loc=jl
1592  end if
1593 
1594 end function bisect_rdp

m_numeric_tools/calculate_pade_a [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  calculate_pade_a

FUNCTION

INPUTS

OUTPUT

SOURCE

4110 subroutine calculate_pade_a(a, n, z, f)
4111 
4112 !Arguments ------------------------------------
4113 !scalars
4114  integer,intent(in) :: n
4115  complex(dpc),intent(in) :: z(n),f(n)
4116  complex(dpc),intent(out) :: a(n)
4117 
4118 !Local variables-------------------------------
4119 !scalars
4120  integer :: i,j
4121 !arrays
4122  complex(dpc) :: g(n,n)
4123 ! *************************************************************************
4124 
4125  g(1,1:n)=f(1:n)
4126 
4127  do i=2,n
4128    do j=i,n
4129      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)
4130      g(i,j)=(g(i-1,i-1)-g(i-1,j)) / ((z(j)-z(i-1))*g(i-1,j))
4131      !write(std_out,*) 'g_i(z_j)',i,j,g(i,j)
4132    end do
4133  end do
4134  do i=1,n
4135    a(i)=g(i,i)
4136  end do
4137  !write(std_out,*) 'a ',a(:)
4138 
4139 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

1194 pure function cdp2rdp_0D(cc) result(rr)
1195 
1196 !Arguments ------------------------------------
1197 !scalars
1198  complex(dpc),intent(in) :: cc
1199  real(dp) :: rr(2)
1200 
1201 ! *********************************************************************
1202 
1203  rr(1)=REAL (cc)
1204  rr(2)=AIMAG(cc)
1205 
1206 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

INPUTS

  cc(:)=the input complex array

OUTPUT

  rr(2,:)=the real array

SOURCE

1226 pure function cdp2rdp_1D(cc) result(rr)
1227 
1228 !Arguments ------------------------------------
1229 !scalars
1230  complex(dpc),intent(in) :: cc(:)
1231  real(dp) :: rr(2,SIZE(cc))
1232 
1233 ! *********************************************************************
1234 
1235  rr(1,:)=REAL (cc(:))
1236  rr(2,:)=AIMAG(cc(:))
1237 
1238 end function cdp2rdp_1D

m_numeric_tools/cdp2rdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_2D

FUNCTION

INPUTS

OUTPUT

SOURCE

1255 pure function cdp2rdp_2D(cc) result(rr)
1256 
1257 !Arguments ------------------------------------
1258 !scalars
1259  complex(dpc),intent(in) :: cc(:,:)
1260  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2))
1261 ! *********************************************************************
1262 
1263  rr(1,:,:)=REAL (cc(:,:))
1264  rr(2,:,:)=AIMAG(cc(:,:))
1265 
1266 end function cdp2rdp_2D

m_numeric_tools/cdp2rdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_3D

FUNCTION

INPUTS

OUTPUT

SOURCE

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

m_numeric_tools/cdp2rdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_4D

FUNCTION

INPUTS

OUTPUT

SOURCE

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

m_numeric_tools/cdp2rdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  cdp2rdp_5D

FUNCTION

INPUTS

OUTPUT

SOURCE

1340 pure function cdp2rdp_5D(cc) result(rr)
1341 
1342 !Arguments ------------------------------------
1343 !scalars
1344  complex(dpc),intent(in) :: cc(:,:,:,:,:)
1345  real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5))
1346 
1347 ! *********************************************************************
1348 
1349  rr(1,:,:,:,:,:)=REAL (cc(:,:,:,:,:))
1350  rr(2,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:))
1351 
1352 end function cdp2rdp_5D

m_numeric_tools/central_finite_diff [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

 central_finite_diff

FUNCTION

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

INPUTS

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

OUTPUT

  coeffient for central finite difference

SOURCE

5575 real(dp) function central_finite_diff(order, ipos, npts) result(fact)
5576 
5577 !Arguments ---------------------------------------------
5578 !scalars
5579  integer,intent(in) :: ipos,order,npts
5580 
5581 !Local variables ---------------------------------------
5582 !scalars
5583  real(dp),parameter :: empty=huge(one)
5584 ! 1st derivative.
5585  real(dp),parameter :: d1(9,4) = reshape([ &
5586   [-1/2._dp, 0._dp, 1/2._dp, empty, empty, empty, empty, empty, empty], &
5587   [ 1/12._dp, -2/3._dp, 0._dp, 2/3._dp, -1/12._dp, empty, empty, empty, empty], &
5588   [-1/60._dp, 3/20._dp, -3/4._dp, 0._dp, 3/4._dp, -3/20._dp, 1/60._dp, empty, empty], &
5589   [ 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])
5590 ! 2nd derivative.
5591  real(dp),parameter :: d2(9,4) = reshape([ &
5592    [ 1._dp, -2._dp, 1._dp, empty, empty, empty, empty, empty, empty], &
5593    [-1/12._dp, 4/3._dp, -5/2._dp, 4/3._dp, -1/12._dp, empty, empty, empty, empty], &
5594    [ 1/90._dp, -3/20._dp, 3/2._dp, -49/18._dp, 3/2._dp, -3/20._dp, 1/90._dp, empty, empty], &
5595    [-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])
5596 ! 3th derivative.
5597  real(dp),parameter :: d3(9,3) = reshape([ &
5598    [-1/2._dp, 1._dp, 0._dp, -1._dp, 1/2._dp, empty, empty, empty, empty], &
5599    [ 1/8._dp, -1._dp, 13/8._dp, 0._dp, -13/8._dp, 1._dp, -1/8._dp, empty, empty], &
5600    [ -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]], &
5601    [9,3])
5602 ! 4th derivative.
5603  real(dp),parameter :: d4(9,3) = reshape([ &
5604    [ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp, empty, empty, empty, empty], &
5605    [ -1/6._dp, 2._dp, -13/2._dp, 28/3._dp, -13/2._dp, 2._dp, -1/6._dp, empty, empty], &
5606    [ 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])
5607 ! 5th derivative.
5608  real(dp),parameter :: d5(7) = [ -1/2._dp, 2._dp, -5/2._dp, 0._dp, 5/2._dp, -2._dp, 1/2._dp]
5609 ! 6th derivative.
5610  real(dp),parameter :: d6(7) = [ 1._dp, -6._dp, 15._dp, -20._dp, 15._dp, -6._dp, 1._dp]
5611 ! *********************************************************************
5612 
5613  select case (order)
5614  case (1)
5615    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5616    fact = d1(ipos, npts/2)
5617  case (2)
5618    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5619    fact = d2(ipos, npts/2)
5620  case (3)
5621    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5622    fact = d3(ipos, npts/2)
5623  case (4)
5624    if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5625    fact = d4(ipos, npts/2 - 1)
5626  case (5)
5627    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5628    fact = d5(ipos)
5629  case (6)
5630    if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5631    fact = d6(ipos)
5632  case default
5633    ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
5634  end select
5635 
5636  if (fact == empty) then
5637    ABI_ERROR(sjoin("Invalid ipos:",itoa(ipos),"for order", itoa(order), "npts", itoa(npts)))
5638  end if
5639  return
5640 
5641 10 ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
5642 
5643 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

3720 integer function check_vec_conjg(nn, vec1, vec2, abs_diff, abs_tol) result(ierr)
3721 
3722 !Arguments ------------------------------------
3723  integer, intent(in) :: nn
3724  complex(dp), intent(in) :: vec1(nn), vec2(nn)
3725  real(dp),intent(out) :: abs_diff(2)
3726  real(dp),optional,intent(in) :: abs_tol
3727 
3728 !Local variables-------------------------------
3729  integer :: ii
3730  real(dp) :: my_abs_tol
3731 
3732  ! *************************************************************************
3733 
3734  my_abs_tol = tol6; if (present(abs_tol)) my_abs_tol = abs_tol
3735  ierr = 0
3736  abs_diff = zero
3737 
3738  do ii=1,nn
3739    abs_diff(1) = max(abs_diff(1), abs(real(vec1(ii)) - real(vec2(ii))))
3740    abs_diff(2) = max(abs_diff(2), abs(aimag(vec1(ii)) + aimag(vec2(ii))))
3741  end do
3742 
3743  if (any(abs_diff > abs_tol)) ierr = 1
3744 
3745 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

4558 subroutine cmplx_sphcart(carr, from, units)
4559 
4560 !Arguments ------------------------------------
4561 !scalars
4562  character(len=*),intent(in) :: from
4563  character(len=*),optional,intent(in) :: units
4564 !arrays
4565  complex(dpc),intent(inout) :: carr(:,:)
4566 
4567 !Local variables-------------------------------
4568 !scalars
4569  integer :: jj,ii
4570  real(dp) :: rho,theta,fact
4571  character(len=500) :: msg
4572 
4573 ! *************************************************************************
4574 
4575  select case (from(1:1))
4576 
4577  case ("S","s") ! Spherical --> Cartesian
4578 
4579    fact = one
4580    if (PRESENT(units)) then
4581      if (units(1:1) == "D" .or. units(1:1) == "d") fact = two_pi/360_dp
4582    end if
4583 
4584    do jj=1,SIZE(carr,DIM=2)
4585      do ii=1,SIZE(carr,DIM=1)
4586         rho  = DBLE(carr(ii,jj))
4587         theta= AIMAG(carr(ii,jj)) * fact
4588         carr(ii,jj) = CMPLX(rho*DCOS(theta), rho*DSIN(theta), kind=dpc)
4589      end do
4590    end do
4591 
4592  case ("C","c") ! Cartesian --> Spherical \theta = 2 arctan(y/(rho+x))
4593 
4594    fact = one
4595    if (PRESENT(units)) then
4596      if (units(1:1) == "D" .or. units(1:1) == "d") fact = 360_dp/two_pi
4597    end if
4598 
4599    do jj=1,SIZE(carr,DIM=2)
4600      do ii=1,SIZE(carr,DIM=1)
4601         rho  = SQRT(ABS(carr(ii,jj)))
4602         if (rho > tol16) then
4603           theta= two * ATAN( AIMAG(carr(ii,jj)) / (DBLE(carr(ii,jj)) + rho) )
4604         else
4605           theta= zero
4606         end if
4607         carr(ii,jj) = CMPLX(rho, theta*fact, kind=dpc)
4608      end do
4609    end do
4610 
4611  case default
4612    msg = " Wrong value for from: "//TRIM(from)
4613    ABI_BUG(msg)
4614  end select
4615 
4616 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

3077 subroutine coeffs_gausslegint(xmin,xmax,x,weights,n)
3078 
3079 !Arguments ------------------------------------
3080 !scalars
3081  integer,intent(in) :: n
3082  real(dp),intent(in) :: xmin,xmax
3083  real(dp),intent(out) :: x(n),weights(n)
3084 
3085 !Local variables ------------------------------
3086 !scalars
3087  integer :: i,j
3088  real(dp),parameter :: tol=1.d-13
3089  real(dp),parameter :: pi=4.d0*atan(1.d0)
3090  real(dp) :: z,z1,xmean,p1,p2,p3,pp,xl
3091 
3092 !************************************************************************
3093 
3094  xl=(xmax-xmin)*0.5d0
3095  xmean=(xmax+xmin)*0.5d0
3096 
3097  do i=1,(n+1)/2
3098   z=cos(pi*(i-0.25d0)/(n+0.5d0))
3099 
3100   do
3101     p1=1.d0
3102     p2=0.d0
3103 
3104     do j=1,n
3105      p3=p2
3106      p2=p1
3107      p1=((2.d0*j - 1.d0)*z*p2 - (j-1.d0)*p3)/j
3108     end do
3109 
3110     pp=n*(p2-z*p1)/(1.0d0-z**2)
3111     z1=z
3112     z=z1-p1/pp
3113 
3114     if(abs(z-z1) < tol) exit
3115   end do
3116 
3117   x(i)=xmean-xl*z
3118   x(n+1-i)=xmean+xl*z
3119   weights(i)=2.d0*xl/((1.d0-z**2)*pp**2)
3120   weights(n+1-i)=weights(i)
3121  end do
3122 
3123 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

4454 subroutine continued_fract(nlev,term_type,aa,bb,nz,zpts,spectrum)
4455 
4456 !Arguments ------------------------------------
4457 !scalars
4458  integer,intent(in) :: nlev,term_type,nz
4459 !arrays
4460  real(dp),intent(in) ::  bb(nlev)
4461  complex(dpc),intent(in) :: aa(nlev)
4462  complex(dpc),intent(in) :: zpts(nz)
4463  complex(dpc),intent(out) :: spectrum(nz)
4464 
4465 !Local variables ------------------------------
4466 !scalars
4467  integer :: it
4468  real(dp) ::  bb_inf,bg,bu,swap
4469  complex(dpc) :: aa_inf
4470  character(len=500) :: msg
4471 !arrays
4472  complex(dpc),allocatable :: div(:),den(:)
4473 
4474 !************************************************************************
4475 
4476  ABI_MALLOC(div,(nz))
4477  ABI_MALLOC(den,(nz))
4478 
4479  select case (term_type)
4480  case (0) ! No terminator.
4481    div=czero
4482  case (-1,1)
4483    if (term_type==-1) then
4484      bb_inf=bb(nlev)
4485      aa_inf=aa(nlev)
4486    else
4487      bb_inf=SUM(bb)/nlev
4488      aa_inf=SUM(aa)/nlev
4489    end if
4490    ! Be careful with the sign of the SQRT.
4491    div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
4492  case (2)
4493    ABI_ERROR("To be tested")
4494    div = zero
4495    if (nlev>4) then
4496      bg=zero; bu=zero
4497      do it=1,nlev,2
4498        if (it+2<nlev) bg = bg + bb(it+2)
4499        bu = bu + bb(it)
4500      end do
4501      bg = bg/(nlev/2+MOD(nlev,2))
4502      bu = bg/((nlev+1)/2)
4503      !if (iseven(nlev)) then
4504      if (.not.iseven(nlev)) then
4505        swap = bg
4506        bg = bu
4507        bu = bg
4508      end if
4509      !write(std_out,*)nlev,bg,bu
4510      !Here be careful with the sign of SQRT
4511      do it=1,nz
4512        div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
4513 &        ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
4514      end do
4515    end if
4516 
4517  case default
4518    write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
4519    ABI_ERROR(msg)
4520  end select
4521 
4522  do it=nlev,2,-1
4523    den(:) = zpts(:) - aa(it) - div(:)
4524    div(:) = (bb(it-1)**2 )/ den(:)
4525  end do
4526 
4527  den = zpts(:) - aa(1) - div(:)
4528  div = one/den(:)
4529 
4530  spectrum = div
4531  ABI_FREE(div)
4532  ABI_FREE(den)
4533 
4534 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

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

2771 subroutine ctrap(imax,ff,hh,ans)
2772 
2773 !Arguments ------------------------------------
2774 !scalars
2775  integer,intent(in) :: imax
2776  real(dp),intent(in) :: hh
2777  real(dp),intent(out) :: ans
2778 !arrays
2779  real(dp),intent(in) :: ff(imax)
2780 
2781 !Local variables-------------------------------
2782 !scalars
2783  integer :: ir,ir2
2784  real(dp) :: endpt,sum
2785 
2786 ! *************************************************************************
2787 
2788  if (imax>=10)then
2789 
2790 !  endpt=end point correction terms (low and high ends)
2791    endpt  = (23.75d0*(ff(1)+ff(imax  )) &
2792 &   + 95.10d0*(ff(2)+ff(imax-1)) &
2793 &   + 55.20d0*(ff(3)+ff(imax-2)) &
2794 &   + 79.30d0*(ff(4)+ff(imax-3)) &
2795 &   + 70.65d0*(ff(5)+ff(imax-4)))/ 72.d0
2796    ir2 = imax - 5
2797    sum=0.00d0
2798    if (ir2 > 5) then
2799      do ir=6,ir2
2800        sum = sum + ff(ir)
2801      end do
2802    end if
2803    ans = (sum + endpt ) * hh
2804 
2805  else if (imax>=8)then
2806    endpt  = (17.0d0*(ff(1)+ff(imax  )) &
2807 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2808 &   + 43.0d0*(ff(3)+ff(imax-2)) &
2809 &   + 49.0d0*(ff(4)+ff(imax-3)) )/ 48.d0
2810    sum=0.0d0
2811    if(imax==9)sum=ff(5)
2812    ans = (sum + endpt ) * hh
2813 
2814  else if (imax==7)then
2815    ans = (17.0d0*(ff(1)+ff(imax  )) &
2816 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2817 &   + 43.0d0*(ff(3)+ff(imax-2)) &
2818 &   + 50.0d0* ff(4)                )/ 48.d0  *hh
2819 
2820  else if (imax==6)then
2821    ans = (17.0d0*(ff(1)+ff(imax  )) &
2822 &   + 59.0d0*(ff(2)+ff(imax-1)) &
2823 &   + 44.0d0*(ff(3)+ff(imax-2)) )/ 48.d0  *hh
2824 
2825  else if (imax==5)then
2826    ans = (     (ff(1)+ff(5)) &
2827 &   + four*(ff(2)+ff(4)) &
2828 &   + two * ff(3)         )/ three  *hh
2829 
2830  else if (imax==4)then
2831    ans = (three*(ff(1)+ff(4)) &
2832 &   + nine *(ff(2)+ff(3))  )/ eight  *hh
2833 
2834  else if (imax==3)then
2835    ans = (     (ff(1)+ff(3)) &
2836 &   + four* ff(2)         )/ three *hh
2837 
2838  else if (imax==2)then
2839    ans = (ff(1)+ff(2))/ two  *hh
2840 
2841  else if (imax==1)then
2842    ans = ff(1)*hh
2843 
2844  end if
2845 
2846 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

4353 integer function denominator(dd,ierr,tolerance)
4354 
4355 !Arguments ------------------------------------
4356 !scalars
4357  integer,intent(out) :: ierr
4358  real(dp),intent(in) :: dd
4359  real(dp),optional,intent(in) :: tolerance
4360 
4361 !Local variables ------------------------------
4362 !scalars
4363  integer,parameter :: largest_integer = HUGE(1)
4364  integer :: ii
4365  real(dp) :: my_tol
4366 
4367 !************************************************************************
4368 
4369  ii=1
4370  my_tol=0.0001 ; if (PRESENT(tolerance)) my_tol=ABS(tolerance)
4371  do
4372    if (ABS(dd*ii-NINT(dd*ii))<my_tol) then
4373      denominator=ii
4374      ierr=0
4375      RETURN
4376    end if
4377    ! Handle the case in which dd is not rational within my_tol.
4378    if (ii==largest_integer) then
4379      denominator=ii
4380      ierr=-1
4381      RETURN
4382    end if
4383    ii=ii+1
4384  end do
4385 
4386 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

4051 function dpade(n, z, f, zz)
4052 
4053 !Arguments ------------------------------------
4054 !scalars
4055  integer,intent(in) :: n
4056  complex(dpc),intent(in) :: zz
4057  complex(dpc) :: dpade
4058 !arrays
4059  complex(dpc),intent(in) :: z(n),f(n)
4060 
4061 !Local variables-------------------------------
4062 !scalars
4063  integer :: i
4064 !arrays
4065  complex(dpc) :: a(n)
4066  complex(dpc) :: Az(0:n), Bz(0:n)
4067  complex(dpc) :: dAz(0:n), dBz(0:n)
4068 ! *************************************************************************
4069 
4070  call calculate_pade_a(a, n, z, f)
4071 
4072  Az(0)=czero
4073  Az(1)=a(1)
4074  Bz(0)=cone
4075  Bz(1)=cone
4076  dAz(0)=czero
4077  dAz(1)=czero
4078  dBz(0)=czero
4079  dBz(1)=czero
4080 
4081  do i=1,n-1
4082    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4083    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4084    dAz(i+1)=dAz(i)+a(i+1)*Az(i-1)+(zz-z(i))*a(i+1)*dAz(i-1)
4085    dBz(i+1)=dBz(i)+a(i+1)*Bz(i-1)+(zz-z(i))*a(i+1)*dBz(i-1)
4086  end do
4087  !write(std_out,*) 'Bz(n)', Bz(n)
4088  !if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) 'Bz(n)',Bz(n)
4089  !pade_approx = Az(n) / Bz(n)
4090  dpade=dAz(n)/Bz(n) -Az(n)*dBz(n)/(Bz(n)*Bz(n))
4091  !write(std_out,*) 'pade_approx ', pade_approx
4092 
4093 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

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

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

m_numeric_tools/get_diag_cdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_cdp

FUNCTION

INPUTS

OUTPUT

SOURCE

817 function get_diag_cdp(cmat) result(cdiag)
818 
819 !Arguments ------------------------------------
820 !scalars
821  complex(dpc),intent(in) :: cmat(:,:)
822  complex(dpc) :: cdiag(SIZE(cmat,1))
823 
824 !Local variables-------------------------------
825  integer :: ii
826 ! *************************************************************************
827 
828  ABI_CHECK(SIZE(cmat,1) == SIZE(cmat,2), 'Matrix not square')
829 
830  do ii=1,SIZE(cmat,1)
831    cdiag(ii)=cmat(ii,ii)
832  end do
833 
834 end function get_diag_cdp

m_numeric_tools/get_diag_int [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_int

FUNCTION

  Return the diagonal of a square matrix as a vector

INPUTS

  matrix(:,:)

OUTPUT

  diag(:)=the diagonal

SOURCE

746 function get_diag_int(mat) result(diag)
747 
748 !Arguments ------------------------------------
749 !scalars
750  integer,intent(in) :: mat(:,:)
751  integer :: diag(SIZE(mat,1))
752 
753 !Local variables-------------------------------
754  integer :: ii
755 ! *************************************************************************
756 
757  ii = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
758 
759  do ii=1,SIZE(mat,1)
760    diag(ii)=mat(ii,ii)
761  end do
762 
763 end function get_diag_int

m_numeric_tools/get_diag_rdp [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  get_diag_rdp

FUNCTION

  Return the diagonal of a square matrix as a vector

INPUTS

  matrix(:,:)

OUTPUT

  diag(:)=the diagonal

SOURCE

783 function get_diag_rdp(mat) result(diag)
784 
785 !Arguments ------------------------------------
786 !scalars
787  real(dp),intent(in) :: mat(:,:)
788  real(dp) :: diag(SIZE(mat,1))
789 
790 !Local variables-------------------------------
791  integer :: ii
792 ! *************************************************************************
793 
794  ABI_CHECK(SIZE(mat,1) == SIZE(mat,2), 'Matrix not square')
795 
796  do ii=1,SIZE(mat,1)
797    diag(ii) = mat(ii,ii)
798  end do
799 
800 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

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

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

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

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

3300 subroutine hermitianize_dpc(mat,uplo)
3301 
3302 !Arguments ------------------------------------
3303 !scalars
3304  character(len=*),intent(in) :: uplo
3305 !arrays
3306  complex(dpc),intent(inout) :: mat(:,:)
3307 
3308 !Local variables-------------------------------
3309 !scalars
3310  integer :: nn,ii,jj
3311 !arrays
3312  complex(dpc),allocatable :: tmp(:)
3313 ! *************************************************************************
3314 
3315  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3316 
3317  select case (uplo(1:1))
3318 
3319  case ("A","a") ! Full matrix has been calculated.
3320    ABI_MALLOC(tmp,(nn))
3321    do ii=1,nn
3322      do jj=ii,nn
3323        tmp(jj)=half*(mat(ii,jj)+DCONJG(mat(jj,ii)))
3324      end do
3325      mat(ii,ii:nn)=tmp(ii:nn)
3326      mat(ii:nn,ii)=DCONJG(tmp(ii:nn))
3327    end do
3328    ABI_FREE(tmp)
3329 
3330  case ("U","u") ! Only the upper triangle is used.
3331    do jj=1,nn
3332      do ii=1,jj
3333        if (ii/=jj) then
3334          mat(jj,ii) = DCONJG(mat(ii,jj))
3335        else
3336          mat(ii,ii) = CMPLX(DBLE(mat(ii,ii)),zero, kind=dpc)
3337        end if
3338      end do
3339    end do
3340 
3341  case ("L","l") ! Only the lower triangle is used.
3342   do jj=1,nn
3343     do ii=1,jj
3344       if (ii/=jj) then
3345         mat(ii,jj) = DCONJG(mat(jj,ii))
3346       else
3347         mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),zero, kind=dpc)
3348       end if
3349     end do
3350   end do
3351 
3352  case default
3353    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3354  end select
3355 
3356 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

3216 subroutine hermitianize_spc(mat,uplo)
3217 
3218 !Arguments ------------------------------------
3219 !scalars
3220  character(len=*),intent(in) :: uplo
3221 !arrays
3222  complex(spc),intent(inout) :: mat(:,:)
3223 
3224 !Local variables-------------------------------
3225 !scalars
3226  integer :: nn,ii,jj
3227 !arrays
3228  complex(spc),allocatable :: tmp(:)
3229 ! *************************************************************************
3230 
3231  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3232 
3233  select case (uplo(1:1))
3234 
3235  case ("A","a") ! Full matrix has been calculated.
3236    ABI_MALLOC(tmp,(nn))
3237    do ii=1,nn
3238      do jj=ii,nn
3239        ! reference half constant is dp not sp
3240        tmp(jj)=real(half)*(mat(ii,jj)+CONJG(mat(jj,ii)))
3241      end do
3242      mat(ii,ii:nn)=tmp(ii:nn)
3243      mat(ii:nn,ii)=CONJG(tmp(ii:nn))
3244    end do
3245    ABI_FREE(tmp)
3246 
3247  case ("U","u") ! Only the upper triangle is used.
3248    do jj=1,nn
3249      do ii=1,jj
3250        if (ii/=jj) then
3251          mat(jj,ii) = CONJG(mat(ii,jj))
3252        else
3253          mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3254        end if
3255      end do
3256    end do
3257 
3258  case ("L","l") ! Only the lower triangle is used.
3259   do jj=1,nn
3260     do ii=1,jj
3261       if (ii/=jj) then
3262         mat(ii,jj) = CONJG(mat(jj,ii))
3263       else
3264         mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3265       end if
3266     end do
3267   end do
3268 
3269  case default
3270    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3271  end select
3272 
3273 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

1663 pure function imax_loc_int(iarr,mask)
1664 
1665 !Arguments ------------------------------------
1666 !scalars
1667  integer :: imax_loc_int
1668 !arrays
1669  integer,intent(in) :: iarr(:)
1670  logical,optional,intent(in) :: mask(:)
1671 
1672 !Local variables-------------------------------
1673  integer :: imax(1)
1674 ! *************************************************************************
1675 
1676  if (PRESENT(mask)) then
1677   imax=MAXLOC(iarr,MASK=mask)
1678  else
1679   imax=MAXLOC(iarr)
1680  end if
1681  imax_loc_int=imax(1)
1682 
1683 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

1699 pure function imax_loc_rdp(arr,mask)
1700 
1701 !Arguments ------------------------------------
1702 !scalars
1703  integer :: imax_loc_rdp
1704 !arrays
1705  real(dp),intent(in) :: arr(:)
1706  logical,optional,intent(in) :: mask(:)
1707 
1708 !Local variables-------------------------------
1709  integer :: imax(1)
1710 ! *************************************************************************
1711 
1712  if (PRESENT(mask)) then
1713   imax=MAXLOC(arr,MASK=mask)
1714  else
1715   imax=MAXLOC(arr)
1716  end if
1717  imax_loc_rdp=imax(1)
1718 
1719 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

1733 pure function imin_loc_int(arr, mask)
1734 
1735 !Arguments ------------------------------------
1736 !scalars
1737  integer :: imin_loc_int
1738 !arrays
1739  integer,intent(in) :: arr(:)
1740  logical,optional,intent(in) :: mask(:)
1741 
1742 !Local variables-------------------------------
1743  integer :: imin(1)
1744 ! *************************************************************************
1745 
1746  if (PRESENT(mask)) then
1747   imin=MINLOC(arr,MASK=mask)
1748  else
1749   imin=MINLOC(arr)
1750  end if
1751  imin_loc_int=imin(1)
1752 
1753 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

1770 pure function imin_loc_rdp(arr,mask)
1771 
1772 !Arguments ------------------------------------
1773 !scalars
1774  integer :: imin_loc_rdp
1775 !arrays
1776  real(dp),intent(in) :: arr(:)
1777  logical,optional,intent(in) :: mask(:)
1778 
1779 !Local variables-------------------------------
1780  integer :: imin(1)
1781 ! *************************************************************************
1782 
1783  if (PRESENT(mask)) then
1784   imin=MINLOC(arr,MASK=mask)
1785  else
1786   imin=MINLOC(arr)
1787  end if
1788 
1789  imin_loc_rdp=imin(1)
1790 
1791 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

1535 pure logical function inrange_dp(xval, win)
1536 
1537 !Arguments ------------------------------------
1538 !scalars
1539  real(dp),intent(in) :: xval, win(2)
1540 ! *************************************************************************
1541 
1542  inrange_dp = (xval >= win(1) .and. xval <= win(2))
1543 
1544 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

1512 pure logical function inrange_int(xval, win)
1513 
1514 !Arguments ------------------------------------
1515 !scalars
1516  integer,intent(in) :: xval,win(2)
1517 ! *************************************************************************
1518 
1519  inrange_int = (xval >= win(1) .and. xval <= win(2))
1520 
1521 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

4907 pure function interpol3d_0d(r, nr1, nr2, nr3, grid) result(res)
4908 
4909 !Arguments-------------------------------------------------------------
4910 !scalars
4911  integer,intent(in) :: nr1, nr2, nr3
4912  real(dp) :: res
4913 !arrays
4914  real(dp),intent(in) :: grid(nr1,nr2,nr3),r(3)
4915 
4916 !Local variables--------------------------------------------------------
4917 !scalars
4918  integer :: ir1,ir2,ir3,pr1,pr2,pr3
4919  real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8
4920  real(dp) :: x1,x2,x3
4921 
4922 ! *************************************************************************
4923 
4924  call interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3, pr1,pr2,pr3)
4925 
4926 !weight
4927  x1=one+r(1)*nr1-real(ir1)
4928  x2=one+r(2)*nr2-real(ir2)
4929  x3=one+r(3)*nr3-real(ir3)
4930 
4931 !calculation of the density value
4932  res1=grid(ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
4933  res2=grid(pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
4934  res3=grid(ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
4935  res4=grid(ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
4936  res5=grid(pr1, pr2, ir3) * x1*x2*(one-x3)
4937  res6=grid(ir1, pr2, pr3) * (one-x1)*x2*x3
4938  res7=grid(pr1, ir2, pr3) * x1*(one-x2)*x3
4939  res8=grid(pr1, pr2, pr3) * x1*x2*x3
4940  res=res1+res2+res3+res4+res5+res6+res7+res8
4941 
4942 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

4968 pure function interpol3d_1d(r, nr1, nr2, nr3, grid, nd) result(res)
4969 
4970 !Arguments-------------------------------------------------------------
4971 !scalars
4972  integer,intent(in) :: nr1, nr2, nr3, nd
4973  real(dp) :: res(nd)
4974 !arrays
4975  real(dp),intent(in) :: grid(nd,nr1,nr2,nr3),r(3)
4976 
4977 !Local variables--------------------------------------------------------
4978 !scalars
4979  integer :: id,ir1,ir2,ir3,pr1,pr2,pr3
4980  real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8
4981  real(dp) :: x1,x2,x3
4982 
4983 ! *************************************************************************
4984 
4985  call interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3, pr1,pr2,pr3)
4986 
4987 !weight
4988  x1=one+r(1)*nr1-real(ir1)
4989  x2=one+r(2)*nr2-real(ir2)
4990  x3=one+r(3)*nr3-real(ir3)
4991 
4992 !calculation of the density value
4993  do id=1,nd
4994    res1=grid(id,ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
4995    res2=grid(id,pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
4996    res3=grid(id,ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
4997    res4=grid(id,ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
4998    res5=grid(id,pr1, pr2, ir3) * x1*x2*(one-x3)
4999    res6=grid(id,ir1, pr2, pr3) * (one-x1)*x2*x3
5000    res7=grid(id,pr1, ir2, pr3) * x1*(one-x2)*x3
5001    res8=grid(id,pr1, pr2, pr3) * x1*x2*x3
5002    res(id)=res1+res2+res3+res4+res5+res6+res7+res8
5003  enddo
5004 
5005 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

5030 pure subroutine interpol3d_indices (r,nr1,nr2,nr3,ir1,ir2,ir3,pr1,pr2,pr3)
5031 
5032 !Arguments-------------------------------------------------------------
5033 !scalars
5034  integer,intent(in) :: nr1,nr2,nr3
5035  integer,intent(out) :: ir1,ir2,ir3
5036  integer,intent(out) :: pr1,pr2,pr3
5037 !arrays
5038  real(dp),intent(in) :: r(3)
5039 
5040 !Local variables-------------------------------
5041  real(dp) :: d1,d2,d3
5042 
5043 ! *************************************************************************
5044 
5045 !grid density
5046  d1=one/nr1
5047  d2=one/nr2
5048  d3=one/nr3
5049 
5050 !lower left
5051  ir1=int(r(1)/d1)+1
5052  ir2=int(r(2)/d2)+1
5053  ir3=int(r(3)/d3)+1
5054 
5055 !upper right
5056  pr1=mod(ir1+1,nr1)
5057  pr2=mod(ir2+1,nr2)
5058  pr3=mod(ir3+1,nr3)
5059 
5060  if(ir1==0) ir1=nr1
5061  if(ir2==0) ir2=nr2
5062  if(ir3==0) ir3=nr3
5063 
5064  if(ir1>nr1) ir1=ir1-nr1
5065  if(ir2>nr2) ir2=ir2-nr2
5066  if(ir3>nr3) ir3=ir3-nr3
5067 
5068  if(pr1==0) pr1=nr1
5069  if(pr2==0) pr2=nr2
5070  if(pr3==0) pr3=nr3
5071 
5072 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

5097 subroutine interpolate_denpot(cplex, in_ngfft, nspden, in_rhor, out_ngfft, out_rhor)
5098 
5099 !Arguments-------------------------------------------------------------
5100 !scalars
5101  integer,intent(in) :: cplex,nspden
5102 !arrays
5103  integer,intent(in) :: in_ngfft(3), out_ngfft(3)
5104  real(dp),intent(in) :: in_rhor(cplex, product(in_ngfft), nspden)
5105  real(dp),intent(out) :: out_rhor(cplex, product(out_ngfft), nspden)
5106 
5107 !Local variables--------------------------------------------------------
5108 !scalars
5109  integer :: ispden, ir1, ir2, ir3, ifft
5110  real(dp) :: rr(3)
5111 
5112 ! *************************************************************************
5113 
5114  ! Linear interpolation.
5115  do ispden=1,nspden
5116    do ir3=0,out_ngfft(3)-1
5117      rr(3) = DBLE(ir3)/out_ngfft(3)
5118      do ir2=0,out_ngfft(2)-1
5119        rr(2) = DBLE(ir2)/out_ngfft(2)
5120        do ir1=0,out_ngfft(1)-1
5121          rr(1) = DBLE(ir1)/out_ngfft(1)
5122          ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5123          out_rhor(1:cplex, ifft, ispden) = interpol3d_1d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_rhor(:, :, ispden),cplex)
5124        end do
5125      end do
5126    end do
5127  end do
5128 
5129 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

6139 subroutine invcb(rhoarr,rspts,npts)
6140 
6141 !Arguments ------------------------------------
6142 !scalars
6143  integer,intent(in) :: npts
6144 !arrays
6145  real(dp),intent(in) :: rhoarr(npts)
6146  real(dp),intent(out) :: rspts(npts)
6147 
6148 !Local variables-------------------------------
6149 !scalars
6150  integer :: ii,ipts
6151  real(dp),parameter :: c2_27=2.0e0_dp/27.0e0_dp,c5_9=5.0e0_dp/9.0e0_dp
6152  real(dp),parameter :: c8_9=8.0e0_dp/9.0e0_dp,m1thrd=-third
6153  real(dp) :: del,prod,rho,rhom1,rhomtrd
6154  logical :: test
6155 !character(len=500) :: message
6156 
6157 ! *************************************************************************
6158 
6159 !Loop over points : here, brute force algorithm
6160 !do ipts=1,npts
6161 !rspts(ipts)=sign( (abs(rhoarr(ipts)))**m1thrd,rhoarr(ipts))
6162 !end do
6163 !
6164 
6165  rhomtrd=sign( (abs(rhoarr(1)))**m1thrd, rhoarr(1) )
6166  rhom1=one/rhoarr(1)
6167  rspts(1)=rhomtrd
6168  do ipts=2,npts
6169    rho=rhoarr(ipts)
6170    prod=rho*rhom1
6171 !  If the previous point is too far ...
6172    if(prod < 0.01_dp .or. prod > 10._dp )then
6173      rhomtrd=sign( (abs(rho))**m1thrd , rho )
6174      rhom1=one/rho
6175    else
6176      del=prod-one
6177      do ii=1,5
6178 !      Choose one of the two next lines, the last one is more accurate
6179 !      rhomtrd=((one+third*del)/(one+two_thirds*del))*rhomtrd
6180        rhomtrd=((one+c5_9*del)/(one+del*(c8_9+c2_27*del)))*rhomtrd
6181        rhom1=rhomtrd*rhomtrd*rhomtrd
6182        del=rho*rhom1-one
6183 !      write(std_out,*)rhomtrd,del
6184        test = del*del < 1.0e-24_dp
6185        if(test) exit
6186      end do
6187      if( .not. test) then
6188        rhomtrd=sign( (abs(rho))**m1thrd , rho )
6189      end if
6190    end if
6191    rspts(ipts)=rhomtrd
6192  end do
6193 
6194 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

1394 pure function is_integer_0d(rr,tol) result(ans)
1395 
1396 !Arguments ------------------------------------
1397 !scalars
1398  real(dp),intent(in) :: tol
1399  logical :: ans
1400 !arrays
1401  real(dp),intent(in) :: rr
1402 
1403 ! *************************************************************************
1404 
1405  ans=(ABS(rr-NINT(rr))<tol)
1406 
1407 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

1424 pure function is_integer_1d(rr,tol) result(ans)
1425 
1426 !Arguments ------------------------------------
1427 !scalars
1428  real(dp),intent(in) :: tol
1429  logical :: ans
1430 !arrays
1431  real(dp),intent(in) :: rr(:)
1432 
1433 ! *************************************************************************
1434 
1435  ans=ALL((ABS(rr-NINT(rr))<tol))
1436 
1437 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

1457 function is_zero_rdp_0d(rr,tol) result(ans)
1458 
1459 !Arguments ------------------------------------
1460 !scalars
1461  real(dp),intent(in) :: tol
1462  logical :: ans
1463 !arrays
1464  real(dp),intent(in) :: rr
1465 ! *************************************************************************
1466 
1467  ans=(ABS(rr)<tol)
1468 
1469 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

1486 function is_zero_rdp_1d(rr,tol) result(ans)
1487 
1488 !Arguments ------------------------------------
1489 !scalars
1490  real(dp),intent(in) :: tol
1491  logical :: ans
1492 !arrays
1493  real(dp),intent(in) :: rr(:)
1494 ! *************************************************************************
1495 
1496  ans=ALL(ABS(rr(:))<tol)
1497 
1498 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

848 pure logical function isdiagmat_int(mat) result(ans)
849 
850 !Arguments ------------------------------------
851 !scalars
852  integer,intent(in) :: mat(:,:)
853 
854 !Local variables-------------------------------
855  integer :: ii,jj
856 ! *************************************************************************
857 
858  ans = .True.
859  do jj=1,size(mat,dim=2)
860    do ii=1,size(mat,dim=1)
861      if (ii == jj) cycle
862      if (mat(ii,jj) /= 0) then
863        ans = .False.; return
864      end if
865    end do
866  end do
867 
868 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

882 pure logical function isdiagmat_rdp(mat, atol) result(ans)
883 
884 !Arguments ------------------------------------
885 !scalars
886  real(dp),intent(in) :: mat(:,:)
887  real(dp),optional,intent(in) :: atol
888 
889 !Local variables-------------------------------
890  integer :: ii,jj
891  real(dp) :: my_atol
892 ! *************************************************************************
893 
894  my_atol = tol12; if (present(atol)) my_atol = atol
895 
896  ans = .True.
897  do jj=1,size(mat,dim=2)
898    do ii=1,size(mat,dim=1)
899      if (ii == jj) cycle
900      if (abs(mat(ii,jj)) > my_atol) then
901        ans = .False.; return
902      end if
903    end do
904  end do
905 
906 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

1366 elemental function iseven(nn)
1367 
1368 !Arguments ------------------------------------
1369 !scalars
1370  integer,intent(in) :: nn
1371  logical :: iseven
1372 ! *********************************************************************
1373 
1374  iseven = ((nn / 2) * 2 == nn)
1375 
1376 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

4695 function isordered_rdp(nn,arr,direction,tol) result(isord)
4696 
4697 !Arguments ------------------------------------
4698 !scalars
4699  integer,intent(in) :: nn
4700  real(dp),intent(in) :: tol
4701  logical :: isord
4702  character(len=*),intent(in) :: direction
4703 !arrays
4704  real(dp),intent(in) :: arr(nn)
4705 
4706 !Local variables ------------------------------
4707 !scalars
4708  integer :: ii
4709  real(dp) :: prev
4710  character(len=500) :: msg
4711 
4712 ! *************************************************************************
4713 
4714  prev = arr(1); isord =.TRUE.
4715 
4716  SELECT CASE (direction(1:1))
4717  CASE(">")
4718  ii=2;
4719  do while (ii<=nn .and. isord)
4720    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) >= prev)
4721    prev = arr(ii)
4722    ii = ii +1
4723  end do
4724 
4725  CASE("<")
4726  ii=2;
4727  do while (ii<=nn .and. isord)
4728    if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) <= prev)
4729    prev = arr(ii)
4730    ii = ii +1
4731  end do
4732 
4733  CASE DEFAULT
4734    msg = "Wrong direction: "//TRIM(direction)
4735    ABI_ERROR(msg)
4736  END SELECT
4737 
4738 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

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

923 pure function l2int_1D(larr) result(int_arr)
924 
925 !Arguments ------------------------------------
926 !scalars
927  logical,intent(in) :: larr(:)
928  integer :: int_arr(size(larr))
929 
930 ! *********************************************************************
931 
932  where (larr)
933    int_arr = 1
934  elsewhere
935    int_arr = 0
936  end where
937 
938 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

955 pure function l2int_2D(larr) result(int_arr)
956 
957 !Arguments ------------------------------------
958 !scalars
959  logical,intent(in) :: larr(:,:)
960  integer :: int_arr(size(larr,1), size(larr,2))
961 
962 ! *********************************************************************
963 
964  where (larr)
965    int_arr = 1
966  elsewhere
967    int_arr = 0
968  end where
969 
970 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

 987 pure function l2int_3D(larr) result(int_arr)
 988 
 989 !Arguments ------------------------------------
 990 !scalars
 991  logical,intent(in) :: larr(:,:,:)
 992  integer :: int_arr(size(larr,1), size(larr,2), size(larr,3))
 993 
 994 ! *********************************************************************
 995 
 996  where (larr)
 997    int_arr = 1
 998  elsewhere
 999    int_arr = 0
1000  end where
1001 
1002 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

1810 integer pure function lfind(mask, back)
1811 
1812 !Arguments ------------------------------------
1813 !scalars
1814  logical,intent(in) :: mask(:)
1815  logical,optional,intent(in) :: back
1816 !arrays
1817 
1818 !Local variables-------------------------------
1819 !scalars
1820  integer :: ii,nitems
1821  logical :: do_back
1822 
1823 !************************************************************************
1824 
1825  do_back = .False.; if (present(back)) do_back = back
1826  lfind = -1; nitems = size(mask); if (nitems == 0) return
1827 
1828  if (do_back) then
1829    ! Backward search
1830    do ii=nitems,1,-1
1831      if (mask(ii)) then
1832        lfind = ii; return
1833      end if
1834    end do
1835  else
1836    ! Forward search.
1837    do ii=1,nitems
1838      if (mask(ii)) then
1839        lfind = ii; return
1840      end if
1841    end do
1842  end if
1843 
1844 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

2130 function linfit_dpc(nn,xx,zz,aa,bb) result(res)
2131 
2132 !Arguments ------------------------------------
2133 !scalars
2134  integer,intent(in) :: nn
2135  real(dp) :: res
2136  real(dp),intent(in) :: xx(nn)
2137  complex(dpc),intent(in) :: zz(nn)
2138  complex(dpc),intent(out) :: aa,bb
2139 !arrays
2140 
2141 !Local variables-------------------------------
2142 !scalars
2143  integer :: ii
2144  real(dp) :: sx,sx2,msrt
2145  complex(dpc) :: sz,sxz
2146 ! *************************************************************************
2147 
2148  sx=zero  ; sx2=zero ; msrt=zero
2149  sz=czero ; sxz=czero
2150  do ii=1,nn
2151   sx=sx+xx(ii)
2152   sz=sz+zz(ii)
2153   sxz=sxz+xx(ii)*zz(ii)
2154   sx2=sx2+xx(ii)*xx(ii)
2155  end do
2156 
2157  aa=(nn*sxz-sx*sz)/(nn*sx2-sx*sx)
2158  bb=sz/nn-sx*aa/nn
2159 
2160  do ii=1,nn
2161   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2162  end do
2163  msrt=SQRT(msrt) ; res=msrt
2164 
2165 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

2022 function linfit_rdp(nn,xx,yy,aa,bb) result(res)
2023 
2024 !Arguments ------------------------------------
2025 !scalars
2026  integer,intent(in) :: nn
2027  real(dp) :: res
2028  real(dp),intent(out) :: aa,bb
2029 !arrays
2030  real(dp),intent(in) :: xx(nn),yy(nn)
2031 
2032 !Local variables-------------------------------
2033 !scalars
2034  integer :: ii
2035  real(dp) :: msrt,sx2,sx,sxy,sy,tx,ty
2036 ! *************************************************************************
2037 
2038  sx=zero ; sy=zero ; sxy=zero ; sx2=zero
2039  do ii=1,nn
2040   tx=xx(ii)
2041   ty=yy(ii)
2042   sx=sx+tx
2043   sy=sy+ty
2044   sxy=sxy+tx*ty
2045   sx2=sx2+tx*tx
2046  end do
2047 
2048  aa=(nn*sxy-sx*sy)/(nn*sx2-sx*sx)
2049  bb=sy/nn-sx*aa/nn
2050 
2051  msrt=zero
2052  do ii=1,nn
2053   tx=xx(ii)
2054   ty=yy(ii)
2055   msrt=msrt+(ty-aa*tx-bb)**2
2056  end do
2057  msrt=SQRT(msrt/nn) ; res=msrt
2058 
2059 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

2077 function linfit_spc(nn,xx,zz,aa,bb) result(res)
2078 
2079 !Arguments ------------------------------------
2080 !scalars
2081  integer,intent(in) :: nn
2082  real(dp) :: res
2083  real(dp),intent(in) :: xx(nn)
2084  complex(spc),intent(in) :: zz(nn)
2085  complex(spc),intent(out) :: aa,bb
2086 !arrays
2087 
2088 !Local variables-------------------------------
2089 !scalars
2090  integer :: ii
2091  real(dp) :: sx,sx2,msrt
2092  complex(dpc) :: sz,sxz
2093 ! *************************************************************************
2094 
2095  sx=zero ; sx2=zero ; msrt=zero
2096  sz=czero ; sxz=czero
2097  do ii=1,nn
2098   sx=sx+xx(ii)
2099   sz=sz+zz(ii)
2100   sxz=sxz+xx(ii)*zz(ii)
2101   sx2=sx2+xx(ii)*xx(ii)
2102  end do
2103 
2104  aa=CMPLX((nn*sxz-sx*sz)/(nn*sx2-sx*sx), kind=spc)
2105  bb=CMPLX(sz/nn-sx*aa/nn, kind=spc)
2106 
2107  do ii=1,nn
2108   msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2109  end do
2110  msrt=SQRT(msrt) ; res=msrt
2111 
2112 end function linfit_spc

m_numeric_tools/linspace [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  linspace

FUNCTION

INPUTS

OUTPUT

SOURCE

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

1871 subroutine list2blocks(list,nblocks,blocks)
1872 
1873 !Arguments ------------------------------------
1874 !scalars
1875  integer,intent(out) :: nblocks
1876  integer,intent(in) :: list(:)
1877 !arrays
1878  integer,intent(out),allocatable :: blocks(:,:)
1879 
1880 !Local variables-------------------------------
1881 !scalars
1882  integer :: ii,nitems
1883 !arrays
1884  integer :: work(2,size(list))
1885 
1886 !************************************************************************
1887 
1888  nitems = size(list)
1889 
1890  ! Handle nitems == 1 case
1891  if (nitems == 1) then
1892    ABI_MALLOC(blocks, (2,1))
1893    blocks = 1
1894    return
1895  end if
1896 
1897  nblocks = 1; work(1,1) = 1
1898 
1899  do ii=2,nitems
1900    if (list(ii) /= (list(ii-1) + 1)) then
1901      work(2,nblocks) = ii - 1
1902      nblocks = nblocks + 1
1903      work(1,nblocks) = ii
1904    end if
1905  end do
1906 
1907  work(2,nblocks) = nitems
1908 
1909  ABI_MALLOC(blocks, (2,nblocks))
1910  blocks = work(:,1:nblocks)
1911 
1912 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

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

1938 subroutine mask2blocks(mask,nblocks,blocks)
1939 
1940 !Arguments ------------------------------------
1941 !scalars
1942  integer,intent(out) :: nblocks
1943  logical,intent(in) :: mask(:)
1944 !arrays
1945  integer,allocatable :: blocks(:,:)
1946 
1947 !Local variables-------------------------------
1948 !scalars
1949  integer :: ii,nitems,start
1950  logical :: inblock
1951 !arrays
1952  integer :: work(2,SIZE(mask))
1953 
1954 !************************************************************************
1955 
1956  ! Find first element.
1957  nitems = size(mask); start = 0
1958  do ii=1,nitems
1959    if (mask(ii)) then
1960      start = ii
1961      exit
1962    end if
1963  end do
1964 
1965  ! Handle no true element or just one.
1966  if (start == 0) then
1967    nblocks = 0
1968    ABI_MALLOC(blocks, (0,0))
1969    return
1970  end if
1971  if (start /= 0 .and. nitems == 1) then
1972    nblocks = 1
1973    ABI_MALLOC(blocks, (2,1))
1974    blocks(:,1) = [1,1]
1975  end if
1976 
1977  nblocks = 1; work(1,1) = start; inblock = .True.
1978 
1979  do ii=start+1,nitems
1980    if (.not.mask(ii)) then
1981      if (inblock) then
1982        inblock = .False.
1983        work(2,nblocks) = ii - 1
1984      end if
1985    else
1986      if (.not. inblock) then
1987        inblock = .True.
1988        nblocks = nblocks + 1
1989        work(1,nblocks) = ii
1990      end if
1991    end if
1992  end do
1993 
1994  if (mask(nitems) .and. inblock) work(2,nblocks) = nitems
1995 
1996  ABI_MALLOC(blocks, (2,nblocks))
1997  blocks = work(:,1:nblocks)
1998 
1999 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

2473  recursive subroutine midpoint_(func,nn,xmin,xmax,quad)
2474 
2475 !Arguments ------------------------------------
2476 !scalars
2477  integer,intent(in) :: nn
2478  !real(dp),external :: func
2479  real(dp),intent(in) :: xmin,xmax
2480  real(dp),intent(inout) :: quad
2481 
2482  interface
2483    function func(x)
2484      use defs_basis
2485      real(dp),intent(in) :: x
2486      real(dp) :: func
2487    end function func
2488  end interface
2489 
2490  !interface
2491  ! function func(x)
2492  !  use defs_basis
2493  !  real(dp),intent(in) :: x(:)
2494  !  real(dp) :: func(SIZE(x))
2495  ! end function func
2496  !end interface
2497 
2498 !Local variables-------------------------------
2499 !scalars
2500  integer  :: npt,ix
2501  real(dp) :: space
2502  character(len=500) :: msg
2503 !arrays
2504  real(dp),allocatable :: xx(:)
2505 
2506 !************************************************************************
2507 
2508  select case (nn)
2509 
2510  case (1)
2511    ! === Initial crude estimate done at the middle of the interval
2512    !quad=(xmax-xmin)*SUM(func((/half*(xmin+xmax)/))) !PARALLEL version
2513    quad=(xmax-xmin)*func(half*(xmin+xmax))
2514 
2515  case (2:)
2516    ! === Add npt interior points, they alternate in spacing between space and 2*space ===
2517    ABI_MALLOC(xx,(2*3**(nn-2)))
2518    npt=3**(nn-2) ; space=(xmax-xmin)/(three*npt)
2519    xx(1:2*npt-1:2)=arth(xmin+half*space,three*space,npt)
2520    xx(2:2*npt:2)=xx(1:2*npt-1:2)+two*space
2521    ! === The new sum is combined with the old integral to give a refined integral ===
2522    !quad=quad/three+space*SUM(func(xx))  !PARALLEL version
2523    quad=quad/three
2524    do ix=1,SIZE(xx)
2525      quad=quad+space*func(xx(ix))
2526    end do
2527    ABI_FREE(xx)
2528 
2529  case (:0)
2530    write(msg,'(a,i3)')' wrong value for nn ',nn
2531    ABI_BUG('Wrong value for nn')
2532  end select
2533 
2534 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

4400 integer function mincm(ii,jj)
4401 
4402 !Arguments ------------------------------------
4403 !scalars
4404  integer,intent(in) :: ii,jj
4405 
4406 !************************************************************************
4407 
4408  if (ii==0.or.jj==0) then
4409    ABI_BUG('ii==0 or jj==0')
4410  end if
4411 
4412  mincm=MAX(ii,jj)
4413  do
4414    if ( ((mincm/ii)*ii)==mincm .and. ((mincm/jj)*jj)==mincm ) RETURN
4415    mincm=mincm+1
4416  end do
4417 
4418 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

3378 pure subroutine mkherm(array,ndim)
3379 
3380 !Arguments -------------------------------
3381 !scalars
3382  integer,intent(in) :: ndim
3383 !arrays
3384  real(dp),intent(inout) :: array(2,ndim,ndim)
3385 
3386 !Local variables -------------------------
3387 !scalars
3388  integer :: i1,i2
3389 
3390 ! *********************************************************************
3391 
3392  do i1=1,ndim
3393    do i2=1,i1
3394      array(1,i1,i2)=(array(1,i1,i2)+array(1,i2,i1))*half
3395      array(2,i1,i2)=(array(2,i1,i2)-array(2,i2,i1))*half
3396      array(1,i2,i1)=array(1,i1,i2)
3397      array(2,i2,i1)=-array(2,i1,i2)
3398    end do
3399  end do
3400 
3401 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

5487 subroutine nderiv(hh,yy,zz,ndim,norder)
5488 
5489 !Arguments ---------------------------------------------
5490 !scalars
5491  integer,intent(in) :: ndim,norder
5492  real(dp),intent(in) :: hh
5493 !arrays
5494  real(dp),intent(in) :: yy(ndim)
5495  real(dp),intent(out) :: zz(ndim)
5496 
5497 !Local variables ---------------------------------------
5498 !scalars
5499  integer :: ier,ii
5500  real(dp) :: aa,bb,cc,h1,y1
5501 
5502 ! *********************************************************************
5503 
5504 !Initialization (common to 1st and 2nd derivative)
5505  h1=one/(12.d0*hh)
5506  y1=yy(ndim-4)
5507 
5508 !FIRST DERIVATIVE
5509 !================
5510  if (norder==1) then
5511 
5512 !  Prepare differentiation loop
5513    bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
5514    cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
5515 !  Start differentiation loop
5516    do ii=5,ndim
5517      aa=bb;bb=cc
5518      cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
5519      zz(ii-4)=aa
5520    end do
5521 !  Normal exit
5522    ier=0
5523    aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
5524    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))
5525    zz(ndim-1)=aa
5526    zz(ndim-2)=cc
5527    zz(ndim-3)=bb
5528 
5529 !  SECOND DERIVATIVE
5530 !  =================
5531  else
5532    h1=h1/hh
5533 !  Prepare differentiation loop
5534    bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
5535    cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
5536 !  Start differentiation loop
5537    do ii=5,ndim
5538      aa=bb;bb=cc
5539      cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
5540      zz(ii-4)=aa
5541    end do
5542 !  Normal exit
5543    ier=0
5544    aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
5545    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))
5546    zz(ndim-1)=aa
5547    zz(ndim-2)=cc
5548    zz(ndim-3)=bb
5549 
5550  end if !norder
5551 
5552 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

4154 complex(dp) function newrap_step(z, f, df)
4155 
4156 !Arguments ------------------------------------
4157 !scalars
4158  complex(dpc),intent(in) :: z,f,df
4159 
4160 !Local variables-------------------------------
4161  real(dp) :: dfm2
4162 ! *************************************************************************
4163 
4164  dfm2=ABS(df)*ABS(df)
4165 
4166  newrap_step = z - (f*CONJG(df))/dfm2
4167  !& z-one/(ABS(df)*ABS(df)) * CMPLX( REAL(f)*REAL(df)+AIMAG(f)*AIMAG(df), -REAL(f)*AIMAG(df)+AIMAG(f)*EAL(df) )
4168 
4169 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

3680 subroutine pack_matrix(mat_in, mat_out, N, cplx)
3681 
3682 !Arguments ------------------------------------
3683 !scalars
3684  integer, intent(in) :: N, cplx
3685  real(dp), intent(in) :: mat_in(cplx, N*N)
3686  real(dp), intent(out) :: mat_out(cplx*N*(N+1)/2)
3687 
3688 !Local variables-------------------------------
3689  integer :: isubh, i, j
3690 
3691 ! *************************************************************************
3692 
3693  isubh = 1
3694  do j=1,N
3695    do i=1,j
3696      mat_out(isubh)    = mat_in(1, (j-1)*N+i)
3697      ! bad for vectorization, but it's not performance critical, so ...
3698      if(cplx == 2) then
3699        mat_out(isubh+1)  = mat_in(2, (j-1)*N+i)
3700      end if
3701      isubh=isubh+cplx
3702    end do
3703  end do
3704 
3705 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

4004 function pade(n, z, f, zz)
4005 
4006 !Arguments ------------------------------------
4007 !scalars
4008  integer,intent(in) :: n
4009  complex(dpc),intent(in) :: zz
4010  complex(dpc) :: pade
4011 !arrays
4012  complex(dpc),intent(in) :: z(n), f(n)
4013 
4014 !Local variables-------------------------------
4015 !scalars
4016  complex(dpc) :: a(n)
4017  complex(dpc) :: Az(0:n), Bz(0:n)
4018  integer :: i
4019 ! *************************************************************************
4020 
4021  call calculate_pade_a(a, n, z, f)
4022 
4023  Az(0)=czero
4024  Az(1)=a(1)
4025  Bz(0)=cone
4026  Bz(1)=cone
4027 
4028  do i=1,n-1
4029    Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4030    Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4031  end do
4032  !write(std_out,*) 'Bz(n)',Bz(n)
4033  !if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) ' Bz(n) ',Bz(n)
4034  pade=Az(n)/Bz(n)
4035  !write(std_out,*) 'pade_approx ', pade_approx
4036 
4037 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

4640 subroutine pfactorize(nn,nfactors,pfactors,powers)
4641 
4642 !Arguments ------------------------------------
4643 !scalars
4644  integer,intent(in) :: nn,nfactors
4645  integer,intent(in) :: pfactors(nfactors)
4646  integer,intent(out) :: powers (nfactors+1)
4647 
4648 !Local variables ------------------------------
4649 !scalars
4650  integer :: tnn,ifc,fact,ipow,maxpwr
4651 
4652 ! *************************************************************************
4653 
4654  powers=0; tnn=nn
4655 
4656  fact_loop: do ifc=1,nfactors
4657    fact = pfactors (ifc)
4658    maxpwr = NINT ( LOG(DBLE(tnn))/LOG(DBLE(fact) ) ) + 1
4659    do ipow=1,maxpwr
4660      if (tnn==1) EXIT fact_loop
4661      if ( MOD(tnn,fact)==0 ) then
4662        tnn=tnn/fact
4663        powers(ifc)=powers(ifc) + 1
4664      end if
4665    end do
4666  end do fact_loop
4667 
4668  if ( nn /= tnn * PRODUCT( pfactors**powers(1:nfactors)) ) then
4669    ABI_BUG('nn/=tnn!')
4670  end if
4671 
4672  powers(nfactors+1) = tnn
4673 
4674 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

2298 subroutine polyn_interp(xa,ya,x,y,dy)
2299 
2300 !Arguments ------------------------------------
2301 !scalars
2302  real(dp),intent(in) :: xa(:),ya(:)
2303  real(dp),intent(in) :: x
2304  real(dp),intent(out) :: y,dy
2305 !Local variables-------------------------------
2306 !scalars
2307  integer :: m,n,ns
2308 !arrays
2309  real(dp),dimension(SIZE(xa)) :: c,d,den,ho
2310 ! *************************************************************************
2311 
2312  n = assert_eq(SIZE(xa),SIZE(ya),'Different size in xa and ya',__FILE__,__LINE__)
2313 
2314  ! === Initialize the tables of c and d ===
2315  c(:)=ya(:) ; d(:)=ya(:) ; ho(:)=xa(:)-x
2316  ! === Find closest table entry and initial approximation to y ===
2317  ns=imin_loc(ABS(x-xa)) ; y=ya(ns)
2318  ns=ns-1
2319  !
2320  ! === For each column of the tableau loop over current c and d and up-date them ===
2321  do m=1,n-1
2322   den(1:n-m)=ho(1:n-m)-ho(1+m:n)
2323   if (ANY(den(1:n-m)==zero)) then
2324    ABI_ERROR('Two input xa are identical')
2325   end if
2326 
2327   den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
2328   d(1:n-m)=ho(1+m:n)*den(1:n-m) ! Update c and d
2329   c(1:n-m)=ho(1:n-m)*den(1:n-m)
2330 
2331   if (2*ns<n-m) then  ! Now decide which correction, c or d, we want to add to the
2332    dy=c(ns+1)         ! accumulating value of y, The last dy added is the error indication.
2333   else
2334    dy=d(ns)
2335    ns=ns-1
2336   end if
2337 
2338   y=y+dy
2339  end do
2340 
2341 end subroutine polyn_interp

m_numeric_tools/print_arr1d_dpc [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  print_arr1d_dpc

FUNCTION

INPUTS

OUTPUT

SOURCE

3829 subroutine print_arr1d_dpc(arr, max_r, unit, mode_paral)
3830 
3831 !Arguments ------------------------------------
3832 !scalars
3833  integer,optional,intent(in) :: unit, max_r
3834  character(len=4),optional,intent(in) :: mode_paral
3835 !arrays
3836  complex(dpc),intent(in) :: arr(:)
3837 
3838 !Local variables-------------------------------
3839 !scalars
3840  integer :: unt,ii,nr,mr
3841  character(len=4) :: mode
3842  character(len=500) :: msg
3843  character(len=100) :: fmth,fmt1
3844 ! *************************************************************************
3845 
3846  unt=std_out ; if (PRESENT(unit      )) unt=unit
3847  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3848  mr=15       ; if (PRESENT(max_r     )) mr=max_r
3849 
3850  if (mode/='COLL'.and.mode/='PERS') then
3851   write(msg,'(2a)')' Wrong value of mode_paral ',mode
3852   ABI_BUG(msg)
3853  end if
3854 
3855  ! Print out matrix.
3856  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
3857 
3858  write(fmth,*)'(6x,',mr,'(i2,6x))'
3859  write(fmt1,*)'(3x,',mr,'f8.3)'
3860 
3861  write(msg,fmth)(ii,ii=1,mr)
3862  call wrtout(unt,msg,mode) ! header
3863  write(msg,fmt1)REAL (arr(1:mr))
3864  call wrtout(unt,msg,mode) !real part
3865  write(msg,fmt1)AIMAG(arr(1:mr))
3866  call wrtout(unt,msg,mode) !imag part
3867 
3868 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

3773 subroutine print_arr1d_spc(arr,max_r,unit,mode_paral)
3774 
3775 !Arguments ------------------------------------
3776 !scalars
3777  integer,optional,intent(in) :: unit,max_r
3778  character(len=4),optional,intent(in) :: mode_paral
3779 !arrays
3780  complex(spc),intent(in) :: arr(:)
3781 
3782 !Local variables-------------------------------
3783 !scalars
3784  integer :: unt,ii,nr,mr
3785  character(len=4) :: mode
3786  character(len=500) :: msg
3787  character(len=100) :: fmth,fmt1
3788 ! *************************************************************************
3789 
3790  unt=std_out ; if (PRESENT(unit      )) unt=unit
3791  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3792  mr=15       ; if (PRESENT(max_r     )) mr=max_r
3793 
3794  if (mode/='COLL'.and.mode/='PERS') then
3795   write(msg,'(2a)')' Wrong value of mode_paral ',mode
3796   ABI_BUG(msg)
3797  end if
3798  !
3799  ! === Print out matrix ===
3800  nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
3801 
3802  write(fmth,*)'(6x,',mr,'(i2,6x))'
3803  write(fmt1,*)'(3x,',mr,'f8.3)'
3804 
3805  write(msg,fmth)(ii,ii=1,mr)
3806  call wrtout(unt,msg,mode) !header
3807  write(msg,fmt1)REAL (arr(1:mr))
3808  call wrtout(unt,msg,mode) !real part
3809  write(msg,fmt1)AIMAG(arr(1:mr))
3810  call wrtout(unt,msg,mode) !imag part
3811 
3812 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

3946 subroutine print_arr2d_dpc(arr,max_r,max_c,unit,mode_paral)
3947 
3948 !Arguments ------------------------------------
3949 !scalars
3950  integer,optional,intent(in) :: unit,max_r,max_c
3951  character(len=4),optional,intent(in) :: mode_paral
3952 !arrays
3953  complex(dpc),intent(in) :: arr(:,:)
3954 
3955 !Local variables-------------------------------
3956 !scalars
3957  integer :: unt,ii,jj,nc,nr,mc,mr
3958  character(len=4) :: mode
3959  character(len=500) :: msg
3960  character(len=100) :: fmth,fmt1,fmt2
3961 ! *************************************************************************
3962 
3963  unt =std_out; if (PRESENT(unit      )) unt =unit
3964  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3965  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
3966  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
3967 
3968  if (mode/='COLL'.and.mode/='PERS') then
3969    write(msg,'(2a)')'Wrong value of mode_paral ',mode
3970    ABI_BUG(msg)
3971  end if
3972  !
3973  ! === Print out matrix ===
3974  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3975  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
3976 
3977  write(fmth,*)'(6x,',mc,'(i2,6x))'
3978  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
3979  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
3980 
3981  write(msg,fmth)(jj,jj=1,mc)
3982  call wrtout(unt, msg, mode) ! header
3983  do ii=1,mr
3984    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
3985    call wrtout(unt,msg,mode) ! real part
3986    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
3987    call wrtout(unt,msg,mode) ! imag part
3988  end do
3989 
3990 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

3885 subroutine print_arr2d_spc(arr, max_r, max_c, unit, mode_paral)
3886 
3887 !Arguments ------------------------------------
3888 !scalars
3889  integer,optional,intent(in) :: unit, max_r, max_c
3890  character(len=4),optional,intent(in) :: mode_paral
3891 !arrays
3892  complex(spc),intent(in) :: arr(:,:)
3893 
3894 !Local variables-------------------------------
3895 !scalars
3896  integer :: unt,ii,jj,nc,nr,mc,mr
3897  character(len=4) :: mode
3898  character(len=500) :: msg
3899  character(len=100) :: fmth,fmt1,fmt2
3900 ! *************************************************************************
3901 
3902  unt =std_out; if (PRESENT(unit      )) unt =unit
3903  mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
3904  mc  =9      ; if (PRESENT(max_c     )) mc  =max_c
3905  mr  =9      ; if (PRESENT(max_r     )) mr  =max_r
3906 
3907  if (mode/='COLL'.and.mode/='PERS') then
3908    write(msg,'(2a)')'Wrong value of mode_paral ',mode
3909    ABI_BUG(msg)
3910  end if
3911  !
3912  ! === Print out matrix ===
3913  nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3914  nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
3915 
3916  write(fmth,*)'(6x,',mc,'(i2,6x))'
3917  write(fmt1,*)'(3x,i2,',mc,'f8.3)'
3918  write(fmt2,*)'(5x   ,',mc,'f8.3,a)'
3919 
3920  write(msg,fmth)(jj,jj=1,mc)
3921  call wrtout(unt,msg,mode) !header
3922  do ii=1,mr
3923    write(msg,fmt1)ii,REAL(arr(ii,1:mc))
3924    call wrtout(unt,msg,mode) !real part
3925    write(msg,fmt2)  AIMAG(arr(ii,1:mc)),ch10
3926    call wrtout(unt,msg,mode) !imag part
3927  end do
3928 
3929 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

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

m_numeric_tools/rdp2cdp_2D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_2D

FUNCTION

INPUTS

OUTPUT

SOURCE

1050 pure function rdp2cdp_2D(rr) result(cc)
1051 
1052 !Arguments ------------------------------------
1053 !scalars
1054  real(dp),intent(in) :: rr(:,:,:)
1055  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3))
1056 
1057 ! *********************************************************************
1058 
1059  cc(:,:)=CMPLX(rr(1,:,:),rr(2,:,:), kind=dpc)
1060 
1061 end function rdp2cdp_2D

m_numeric_tools/rdp2cdp_3D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_3D

FUNCTION

INPUTS

OUTPUT

SOURCE

1078 pure function rdp2cdp_3D(rr) result(cc)
1079 
1080 !Arguments ------------------------------------
1081 !scalars
1082  real(dp),intent(in) :: rr(:,:,:,:)
1083  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4))
1084 
1085 ! *********************************************************************
1086 
1087  cc(:,:,:)=CMPLX(rr(1,:,:,:),rr(2,:,:,:), kind=dpc)
1088 
1089 end function rdp2cdp_3D

m_numeric_tools/rdp2cdp_4D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_4D

FUNCTION

INPUTS

OUTPUT

SOURCE

1106 pure function rdp2cdp_4D(rr) result(cc)
1107 
1108 !Arguments ------------------------------------
1109 !scalars
1110  real(dp),intent(in) :: rr(:,:,:,:,:)
1111  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5))
1112 
1113 ! *********************************************************************
1114 
1115  cc(:,:,:,:)=CMPLX(rr(1,:,:,:,:),rr(2,:,:,:,:), kind=dpc)
1116 
1117 end function rdp2cdp_4D

m_numeric_tools/rdp2cdp_5D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_5D

FUNCTION

INPUTS

OUTPUT

SOURCE

1134 pure function rdp2cdp_5D(rr) result(cc)
1135 
1136 !Arguments ------------------------------------
1137 !scalars
1138  real(dp),intent(in) :: rr(:,:,:,:,:,:)
1139  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6))
1140 
1141 ! *********************************************************************
1142 
1143  cc(:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:),rr(2,:,:,:,:,:), kind=dpc)
1144 
1145 end function rdp2cdp_5D

m_numeric_tools/rdp2cdp_6D [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  rdp2cdp_6D

FUNCTION

INPUTS

OUTPUT

SOURCE

1162 pure function rdp2cdp_6D(rr) result(cc)
1163 
1164 !Arguments ------------------------------------
1165 !scalars
1166  real(dp),intent(in) :: rr(:,:,:,:,:,:,:)
1167  complex(dpc) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6),SIZE(rr,7))
1168 
1169 ! *********************************************************************
1170 
1171  cc(:,:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:,:),rr(2,:,:,:,:,:,:), kind=dpc)
1172 
1173 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

4269 subroutine remove_copies(n_in, set_in, n_out, is_equal)
4270 
4271 !Arguments ------------------------------------
4272 !scalars
4273  integer,intent(in) :: n_in
4274  integer,intent(out) :: n_out
4275 !arrays
4276  real(dp),target,intent(inout) :: set_in(3,n_in)
4277 
4278  interface
4279    function is_equal(k1,k2)
4280      use defs_basis
4281      real(dp),intent(in) :: k1(3),k2(3)
4282      logical :: is_equal
4283    end function is_equal
4284  end interface
4285 
4286 !Local variables-------------------------------
4287 !scalars
4288  integer :: ii,jj
4289  logical :: isnew
4290 !arrays
4291  type rdp1d_pt
4292   integer :: idx
4293   real(dp),pointer :: rpt(:)
4294  end type rdp1d_pt
4295  type(rdp1d_pt),allocatable :: Ap(:)
4296 
4297 ! *************************************************************************
4298 
4299  ABI_MALLOC(Ap,(n_in))
4300  Ap(1)%idx = 1
4301  Ap(1)%rpt => set_in(:,1)
4302 
4303  n_out=1
4304  do ii=2,n_in
4305 
4306    isnew=.TRUE.
4307    do jj=1,n_out
4308      if (is_equal(set_in(:,ii),Ap(jj)%rpt(:))) then
4309        isnew=.FALSE.
4310        exit
4311      end if
4312    end do
4313 
4314    if (isnew) then
4315      n_out=n_out+1
4316      Ap(n_out)%rpt => set_in(:,ii)
4317      Ap(n_out)%idx = ii
4318    end if
4319  end do
4320 
4321  ! The n_out inequivalent items are packed first.
4322  if (n_out/=n_in) then
4323    do ii=1,n_out
4324      jj=Ap(ii)%idx
4325      set_in(:,ii) = set_in(:,jj)
4326      !write(std_out,*) Ap(ii)%idx,Ap(ii)%rpt(:)
4327    end do
4328  end if
4329 
4330  ABI_FREE(Ap)
4331 
4332 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

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

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

5257 pure subroutine rhophi(cx, phi, rho)
5258 
5259 !Arguments ------------------------------------
5260 !scalars
5261  real(dp),intent(out) :: phi,rho
5262 !arrays
5263  real(dp),intent(in) :: cx(2)
5264 
5265 ! ***********************************************************************
5266 
5267  rho = sqrt(cx(1)*cx(1) + cx(2)*cx(2))
5268 
5269  if (abs(cx(1)) > tol8) then
5270 
5271    phi = atan(cx(2)/cx(1))
5272 
5273    ! phi is an element of [-pi,pi]
5274    if (cx(1) < zero) then
5275      if (phi < zero) then
5276        phi = phi + pi
5277      else
5278        phi = phi - pi
5279      end if
5280    end if
5281 
5282  else
5283 
5284    if (cx(2) > tol8) then
5285      phi = pi*half
5286    else if (cx(2) < tol8) then
5287      phi = -pi*half
5288    else
5289      phi = zero
5290    end if
5291 
5292  end if
5293 
5294 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

5217 function simpson(step, values) result(res)
5218 
5219 !Arguments ------------------------------------
5220 !scalars
5221  real(dp),intent(in) :: step
5222  real(dp) :: res
5223 !arrays
5224  real(dp),intent(in) :: values(:)
5225 
5226 !Local variables -------------------------
5227 !scalars
5228  real(dp) :: int_values(size(values))
5229 
5230 ! *********************************************************************
5231 
5232  call simpson_int(size(values),step,values,int_values)
5233  res = int_values(size(values))
5234 
5235 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

3150 function simpson_cplx(npts,step,ff)
3151 
3152 !Arguments ------------------------------------
3153 !scalars
3154  integer,intent(in) :: npts
3155  real(dp),intent(in)  :: step
3156  complex(dpc),intent(in) :: ff(npts)
3157  complex(dpc) :: simpson_cplx
3158 
3159 !Local variables ------------------------------
3160 !scalars
3161  integer :: ii,my_n
3162  complex(dpc) :: sum_even, sum_odd
3163 
3164 !************************************************************************
3165 
3166  my_n=npts; if ((npts/2)*2 == npts) my_n=npts-3
3167 
3168  if (my_n<2) then
3169    ABI_ERROR("Too few points")
3170  end if
3171 
3172  sum_odd=czero
3173  do ii=2,my_n-1,2
3174    sum_odd = sum_odd + ff(ii)
3175  end do
3176 
3177  sum_even=zero
3178  do ii=3,my_n-2,2
3179    sum_even = sum_even + ff(ii)
3180  end do
3181 
3182  ! Eq 25.4.6 Abramowitz. Error is O(step^4)
3183  simpson_cplx = step/three * (ff(1) + four*sum_odd + two*sum_even + ff(my_n))
3184 
3185  if (my_n/=npts) then ! Simpson's 3/8 rule. Eq 25.4.13 Abramowitz. Error is O(step^5)
3186   simpson_cplx = simpson_cplx + three*step/eight * (ff(npts-3) + 3*ff(npts-2) + 3*ff(npts-1) + ff(npts))
3187  end if
3188 
3189 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

5151 subroutine simpson_int(npts, step, values, int_values)
5152 
5153 !Arguments ------------------------------------
5154 !scalars
5155  integer,intent(in) :: npts
5156  real(dp),intent(in) :: step
5157 !arrays
5158  real(dp),intent(in) :: values(npts)
5159  real(dp),intent(out) :: int_values(npts)
5160 
5161 !Local variables -------------------------
5162 !scalars
5163  integer :: ii
5164  real(dp),parameter :: coef1 = 0.375_dp                          !9.0_dp  / 24.0_dp
5165  real(dp),parameter :: coef2 = 1.166666666666666666666666667_dp  !28.0_dp / 24.0_dp
5166  real(dp),parameter :: coef3 = 0.958333333333333333333333333_dp  !23.0_dp / 24.0_dp
5167  character(len=500) :: msg
5168 
5169 ! *********************************************************************
5170 
5171  if (npts < 6) then
5172    write(msg,"(a,i0)")"Number of points in integrand function must be >=6 while it is: ",npts
5173    ABI_ERROR(msg)
5174  end if
5175 
5176 !-----------------------------------------------------------------
5177 !Simpson integral of input function
5178 !-----------------------------------------------------------------
5179 
5180 !first point is 0: don t store it
5181 !do integration equivalent to Simpson O(1/N^4) from NumRec in C p 134  NumRec in Fortran p 128
5182  int_values(1) =               coef1*values(1)
5183  int_values(2) = int_values(1) + coef2*values(2)
5184  int_values(3) = int_values(2) + coef3*values(3)
5185 
5186  do ii=4,npts-3
5187    int_values(ii) = int_values(ii-1) + values(ii)
5188  end do
5189 
5190  int_values(npts-2) = int_values(npts-3) + coef3*values(npts-2)
5191  int_values(npts-1) = int_values(npts-2) + coef2*values(npts-1)
5192  int_values(npts  ) = int_values(npts-1) + coef1*values(npts  )
5193 
5194  int_values(:) = int_values(:) * step
5195 
5196 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

5426 subroutine smooth(a, mesh, it)
5427 
5428 !Arguments ------------------------------------
5429 !scalars
5430  integer, intent(in) :: it,mesh
5431  real(dp), intent(inout) :: a(mesh)
5432 !Local variables-------------------------------
5433  integer :: i,k
5434  real(dp) :: asm(mesh)
5435 ! *********************************************************************
5436 
5437  do k=1,it
5438    asm(1)=1.0d0/3.0d0*(a(1)+a(2)+a(3))
5439    asm(2)=0.25d0*(a(1)+a(2)+a(3)+a(4))
5440    asm(3)=0.2d0*(a(1)+a(2)+a(3)+a(4)+a(5))
5441    asm(4)=0.2d0*(a(2)+a(3)+a(4)+a(5)+a(6))
5442    asm(5)=0.2d0*(a(3)+a(4)+a(5)+a(6)+a(7))
5443    asm(mesh-4)=0.2d0*(a(mesh-2)+a(mesh-3)+a(mesh-4)+&
5444 &                   a(mesh-5)+a(mesh-6))
5445    asm(mesh-3)=0.2d0*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
5446 &                   a(mesh-4)+a(mesh-5))
5447    asm(mesh-2)=0.2d0*(a(mesh)+a(mesh-1)+a(mesh-2)+&
5448 &                   a(mesh-3)+a(mesh-4))
5449    asm(mesh-1)=0.25d0*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
5450    asm(mesh)=1.0d0/3.0d0*(a(mesh)+a(mesh-1)+a(mesh-2))
5451 
5452    do i=6,mesh-5
5453      asm(i)=0.1d0*a(i)+0.1d0*(a(i+1)+a(i-1))+&
5454 &             0.1d0*(a(i+2)+a(i-2))+&
5455 &             0.1d0*(a(i+3)+a(i-3))+&
5456 &             0.1d0*(a(i+4)+a(i-4))+&
5457 &             0.05d0*(a(i+5)+a(i-5))
5458    end do
5459 
5460    do i=1,mesh
5461      a(i)=asm(i)
5462    end do
5463  end do
5464 
5465 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

4758 pure function stats_eval(arr) result(stats)
4759 
4760 !Arguments ------------------------------------
4761 !scalars
4762  type(stats_t) :: stats
4763 !arrays
4764  real(dp),intent(in) :: arr(:)
4765 
4766 !Local variables ------------------------------
4767 !scalars
4768  integer :: ii,nn
4769  real(dp) :: xx,x2_sum
4770 
4771 ! *************************************************************************
4772 
4773  stats%min   = +HUGE(one)
4774  stats%max   = -HUGE(one)
4775  stats%mean  = zero
4776 
4777  nn = SIZE(arr)
4778  do ii=1,nn
4779    xx = arr(ii)
4780    stats%max  = MAX(stats%max, xx)
4781    stats%min  = MIN(stats%min, xx)
4782    stats%mean = stats%mean + xx
4783  end do
4784 
4785  stats%mean  = stats%mean/nn
4786 
4787  ! Two-pass algorithm for the variance (more stable than the single-pass one).
4788  x2_sum = zero
4789  do ii=1,nn
4790    xx     = arr(ii)
4791    x2_sum = x2_sum + (xx - stats%mean)*(xx - stats%mean)
4792  end do
4793 
4794  if (nn>1) then
4795    stats%stdev  = x2_sum/(nn-1)
4796    stats%stdev = SQRT(ABS(stats%stdev))
4797  else
4798    stats%stdev = zero
4799  end if
4800 
4801 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

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

3613 subroutine symmetrize_dpc(mat, uplo)
3614 
3615 !Arguments ------------------------------------
3616 !scalars
3617  character(len=*),intent(in) :: uplo
3618 !arrays
3619  complex(dpc),intent(inout) :: mat(:,:)
3620 
3621 !Local variables-------------------------------
3622 !scalars
3623  integer :: nn,ii,jj
3624 !arrays
3625  complex(dpc),allocatable :: tmp(:)
3626 ! *************************************************************************
3627 
3628  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3629 
3630  select case (uplo(1:1))
3631  case ("A","a") ! Full matrix has been calculated.
3632    ABI_MALLOC(tmp,(nn))
3633    do ii=1,nn
3634      do jj=ii,nn
3635        tmp(jj)=half*(mat(ii,jj)+mat(jj,ii))
3636      end do
3637      mat(ii,ii:nn)=tmp(ii:nn)
3638      mat(ii:nn,ii)=tmp(ii:nn)
3639    end do
3640    ABI_FREE(tmp)
3641 
3642  case ("U","u") ! Only the upper triangle is used.
3643    do jj=1,nn
3644      do ii=1,jj-1
3645        mat(jj,ii) = mat(ii,jj)
3646      end do
3647    end do
3648 
3649  case ("L","l") ! Only the lower triangle is used.
3650   do jj=1,nn
3651     do ii=1,jj-1
3652       mat(ii,jj) = mat(jj,ii)
3653     end do
3654   end do
3655 
3656  case default
3657    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3658  end select
3659 
3660 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

3538 subroutine symmetrize_spc(mat,uplo)
3539 
3540 !Arguments ------------------------------------
3541 !scalars
3542  character(len=*),intent(in) :: uplo
3543 !arrays
3544  complex(spc),intent(inout) :: mat(:,:)
3545 
3546 !Local variables-------------------------------
3547 !scalars
3548  integer :: nn,ii,jj
3549 !arrays
3550  complex(spc),allocatable :: tmp(:)
3551 ! *************************************************************************
3552 
3553  nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3554 
3555  select case (uplo(1:1))
3556 
3557  case ("A","a") ! Full matrix has been calculated.
3558    ABI_MALLOC(tmp,(nn))
3559    do ii=1,nn
3560      do jj=ii,nn
3561        tmp(jj)=REAL(half)*(mat(ii,jj)+mat(jj,ii))
3562      end do
3563      mat(ii,ii:nn)=tmp(ii:nn)
3564      mat(ii:nn,ii)=tmp(ii:nn)
3565    end do
3566    ABI_FREE(tmp)
3567 
3568  case ("U","u") ! Only the upper triangle is used.
3569    do jj=1,nn
3570      do ii=1,jj-1
3571        mat(jj,ii) = mat(ii,jj)
3572      end do
3573    end do
3574 
3575  case ("L","l") ! Only the lower triangle is used.
3576   do jj=1,nn
3577     do ii=1,jj-1
3578       mat(ii,jj) = mat(jj,ii)
3579     end do
3580   end do
3581 
3582  case default
3583    ABI_ERROR("Wrong uplo"//TRIM(uplo))
3584  end select
3585 
3586 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

2373 recursive subroutine trapezoidal_(func,nn,xmin,xmax,quad)
2374 
2375 !Arguments ------------------------------------
2376 !scalars
2377  integer,intent(in) :: nn
2378  !real(dp),external :: func
2379  real(dp),intent(in) :: xmin,xmax
2380  real(dp),intent(inout) :: quad
2381 
2382  interface
2383    function func(x)
2384      use defs_basis
2385      real(dp),intent(in) :: x
2386      real(dp) :: func
2387    end function func
2388  end interface
2389 
2390  !interface
2391  ! function func(x)
2392  !  use defs_basis
2393  !  real(dp),intent(in) :: x(:)
2394  !  real(dp) :: func(SIZE(x))
2395  ! end function func
2396  !end interface
2397 
2398 !Local variables-------------------------------
2399 !scalars
2400  integer :: npt,ix
2401  real(dp) :: space,new,yy
2402  character(len=500) :: msg
2403 !arrays
2404  !real(dp),allocatable :: xx(:)
2405 !************************************************************************
2406 
2407  select case (nn)
2408 
2409  case (1)
2410    ! === Initial crude estimate (xmax-xmin)(f1+f2)/2 ===
2411    !quad=half*(xmax-xmin)*SUM(func((/xmin,xmax/)))
2412    quad=half*(xmax-xmin)*(func(xmin)+func(xmax))
2413 
2414  case (2:)
2415    ! === Add npt interior points of spacing space ===
2416    npt=2**(nn-2) ; space=(xmax-xmin)/npt
2417    ! === The new sum is combined with the old integral to give a refined integral ===
2418    !new=SUM(func(arth(xmin+half*space,space,npt))) !PARALLEL version
2419    !allocate(xx(npt))
2420    !xx(:)=arth(xmin+half*space,space,npt)
2421    !xx(1)=xmin+half*space
2422    !do ii=2,nn
2423    ! xx(ii)=xx(ii-1)+space
2424    !end do
2425    new=zero
2426    yy=xmin+half*space
2427    do ix=1,npt
2428     !new=new+func(xx(ix))
2429     new=new+func(yy)
2430     yy=yy+space
2431    end do
2432    !deallocate(xx)
2433    quad=half*(quad+space*new)
2434    !write(std_out,*) 'trapezoidal',quad
2435 
2436  case (:0)
2437    write(msg,'(a,i3)')'Wrong value for nn ',nn
2438    ABI_BUG(msg)
2439  end select
2440 
2441 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

5664 function uniformrandom(seed)
5665 
5666 !Arguments ------------------------------------
5667 !scalars
5668  real(dp) :: uniformrandom
5669  integer,intent(inout) :: seed
5670 
5671 !Local variables ---------------------------------------
5672  integer, parameter :: im1=11979,ia1= 430,ic1=2531
5673  integer, parameter :: im2= 6655,ia2= 936,ic2=1399
5674  integer, parameter :: im3= 6075,ia3=1366,ic3=1283
5675  integer, save :: init=0
5676  integer, save :: ii1,ii2,ii3
5677  integer :: kk
5678  real(dp) :: im1inv,im2inv
5679  real(dp), save :: table(97)
5680  character(len=500) :: msg
5681 
5682 ! *********************************************************************
5683 
5684  im1inv=1.0d0/im1 ; im2inv=1.0d0/im2
5685 
5686 !Initialize on first call or when seed<0:
5687  if (seed<0.or.init==0) then
5688    seed=-abs(seed)
5689 
5690 !  First generator
5691    ii1=mod(ic1-seed,im1)
5692    ii1=mod(ia1*ii1+ic1,im1)
5693 !  Second generator
5694    ii2=mod(ii1,im2)
5695    ii1=mod(ia1*ii1+ic1,im1)
5696 !  Third generator
5697    ii3=mod(ii1,im3)
5698 
5699 !  Fill table
5700    do kk=1,97
5701      ii1=mod(ia1*ii1+ic1,im1)
5702      ii2=mod(ia2*ii2+ic2,im2)
5703      table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
5704    enddo
5705 
5706    init=1 ; seed=1
5707  end if
5708 
5709 !Third generator gives index
5710  ii3=mod(ia3*ii3+ic3,im3)
5711  kk=1+(97*ii3)/im3
5712  if (kk<1.or.kk>97) then
5713    write(msg,'(a,2i0,a)' ) ' trouble in uniformrandom; ii3,kk=',ii3,kk,' =>stop'
5714    ABI_ERROR(msg)
5715  end if
5716  uniformrandom=table(kk)
5717 
5718 !Replace old value, based on generators 1 and 2
5719  ii1=mod(ia1*ii1+ic1,im1)
5720  ii2=mod(ia2*ii2+ic2,im2)
5721  table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
5722 
5723 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

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

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

573 pure subroutine unit_matrix_rdp(matrix)
574 
575 !Arguments ------------------------------------
576  real(dp),intent(inout) :: matrix(:,:)
577 
578 !Local variables-------------------------------
579 !scalars
580  integer :: ii,nn
581 ! *********************************************************************
582 
583  nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
584  matrix(:,:)=zero
585  do ii=1,nn
586    matrix(ii,ii)=one
587  end do
588 
589 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

5318 type(vdiff_t) function vdiff_eval(cplex, nr, f1, f2, volume, vd_max) result(vd)
5319 
5320 !Arguments ------------------------------------
5321 !scalars
5322  integer,intent(in) :: cplex,nr
5323  real(dp),intent(in) :: volume
5324  type(vdiff_t),optional,intent(inout) :: vd_max
5325 !arrays
5326  real(dp),intent(in) :: f1(cplex,nr),f2(cplex,nr)
5327 
5328 !Local variables-------------------------------
5329 !scalars
5330  integer :: ir
5331  real(dp) :: num,den,dr
5332  type(stats_t) :: stats
5333 !arrays
5334  real(dp) :: abs_diff(nr)
5335 ! *********************************************************************
5336 
5337  dr = volume / nr
5338 
5339  if (cplex == 1) then
5340    abs_diff = abs(f1(1,:) - f2(1,:))
5341    num = sum(abs_diff)
5342    den = sum(abs(f2(1,:)))
5343 
5344  else if (cplex == 2) then
5345    do ir=1,nr
5346      abs_diff(ir) = sqrt((f1(1,ir) - f2(1,ir))**2 + (f1(2,ir) - f2(2,ir))**2)
5347    end do
5348    num = sum(abs_diff)
5349    den = zero
5350    do ir=1,nr
5351      den = den + sqrt(f2(1,ir)**2 + f2(2,ir)**2)
5352    end do
5353  end if
5354 
5355  vd%int_adiff = num * dr
5356  call safe_div(num,den,zero,vd%l1_rerr)
5357 
5358  stats = stats_eval(abs_diff)
5359  vd%mean_adiff = stats%mean
5360  vd%stdev_adiff = stats%stdev
5361  vd%min_adiff = stats%min
5362  vd%max_adiff = stats%max
5363 
5364  if (present(vd_max)) then
5365    vd_max%int_adiff =   max(vd_max%int_adiff, vd%int_adiff)
5366    vd_max%mean_adiff =  max(vd_max%mean_adiff, vd%mean_adiff)
5367    vd_max%stdev_adiff = max(vd_max%stdev_adiff, vd%stdev_adiff)
5368    vd_max%min_adiff =   max(vd_max%min_adiff, vd%min_adiff)
5369    vd_max%max_adiff =   max(vd_max%max_adiff, vd%max_adiff)
5370    vd_max%l1_rerr =     max(vd_max%l1_rerr, vd%L1_rerr)
5371  end if
5372 
5373 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

5387 subroutine vdiff_print(vd, unit)
5388 
5389 !Arguments ------------------------------------
5390 !scalars
5391  integer,optional,intent(in) :: unit
5392  type(vdiff_t),intent(in) :: vd
5393 
5394 !Local variables-------------------------------
5395 !scalars
5396  integer :: unt
5397 ! *********************************************************************
5398 
5399  unt = std_out; if (present(unit)) unt = unit
5400  write(unt,"(a,es10.3,a)")"  L1_rerr: ", vd%l1_rerr, ","
5401  write(unt,"(a,es10.3,a)")"  'Integral |f1-f2|dr': ", vd%int_adiff, ","
5402  write(unt,"(a,es10.3,a)")"  'min {|f1-f2|}': ", vd%min_adiff, ","
5403  write(unt,"(a,es10.3,a)")"  'Max {|f1-f2|}': ", vd%max_adiff, ","
5404  write(unt,"(a,es10.3,a)")"  'mean {|f1-f2|}': ", vd%mean_adiff, ","
5405  write(unt,"(a,es10.3,a)")"  'stdev {|f1-f2|}': ", vd%stdev_adiff, ","
5406 
5407 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

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

4864 elemental subroutine wrap2_pmhalf(num,red,shift)
4865 
4866 !Arguments -------------------------------
4867 !scalars
4868  real(dp),intent(in) :: num
4869  real(dp),intent(out) :: red,shift
4870 
4871 ! *********************************************************************
4872 
4873  if (num>zero) then
4874    red=mod((num+half-tol12),one)-half+tol12
4875  else
4876    red=-mod(-(num-half-tol12),one)+half+tol12
4877  end if
4878  if(abs(red)<tol12)red=0.0d0
4879  shift=num-red
4880 
4881 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

4824 elemental subroutine wrap2_zero_one(num, red, shift)
4825 
4826 !Arguments ------------------------------------
4827 !scalars
4828  real(dp),intent(in) :: num
4829  real(dp),intent(out) :: red,shift
4830 
4831 ! *************************************************************************
4832 
4833  if (num>zero) then
4834    red=mod((num+tol12),one)-tol12
4835  else
4836    red=-mod(-(num-one+tol12),one)+one-tol12
4837  end if
4838  if(abs(red)<tol12)red=0.0_dp
4839  shift=num-red
4840 
4841 end subroutine wrap2_zero_one