TABLE OF CONTENTS
- ABINIT/blocked_loop
- ABINIT/bool2index
- ABINIT/dotproduct
- ABINIT/m_numeric_tools
- ABINIT/polynomial_regression
- ABINIT/safe_div
- m_numeric_tools/arth_int
- m_numeric_tools/arth_rdp
- m_numeric_tools/bisect_int
- m_numeric_tools/bisect_rdp
- m_numeric_tools/calculate_pade_a
- m_numeric_tools/cdp2rdp_0D
- m_numeric_tools/cdp2rdp_1D
- m_numeric_tools/cdp2rdp_2D
- m_numeric_tools/cdp2rdp_3D
- m_numeric_tools/cdp2rdp_4D
- m_numeric_tools/cdp2rdp_5D
- m_numeric_tools/central_finite_diff
- m_numeric_tools/check_vec_conjg
- m_numeric_tools/cmplx_sphcart
- m_numeric_tools/coeffs_gausslegint
- m_numeric_tools/continued_fract
- m_numeric_tools/cross_product_int
- m_numeric_tools/cross_product_rdp
- m_numeric_tools/cspint
- m_numeric_tools/ctrap
- m_numeric_tools/denominator
- m_numeric_tools/dpade
- m_numeric_tools/findmin
- m_numeric_tools/geop
- m_numeric_tools/get_diag_cdp
- m_numeric_tools/get_diag_int
- m_numeric_tools/get_diag_rdp
- m_numeric_tools/get_trace_cdp
- m_numeric_tools/get_trace_int
- m_numeric_tools/get_trace_rdp
- m_numeric_tools/hermit
- m_numeric_tools/hermitianize_dpc
- m_numeric_tools/hermitianize_spc
- m_numeric_tools/imax_loc_int
- m_numeric_tools/imax_loc_rdp
- m_numeric_tools/imin_loc_int
- m_numeric_tools/imin_loc_rdp
- m_numeric_tools/inrange_dp
- m_numeric_tools/inrange_int
- m_numeric_tools/interpol3d_0d
- m_numeric_tools/interpol3d_1d
- m_numeric_tools/interpol3d_indices
- m_numeric_tools/interpolate_denpot
- m_numeric_tools/invcb
- m_numeric_tools/is_integer_0D
- m_numeric_tools/is_integer_1D
- m_numeric_tools/is_zero_rdp_0D
- m_numeric_tools/is_zero_rdp_1d
- m_numeric_tools/isdiagmat_int
- m_numeric_tools/isdiagmat_rdp
- m_numeric_tools/iseven
- m_numeric_tools/isordered
- m_numeric_tools/kramerskronig
- m_numeric_tools/l2int_1D
- m_numeric_tools/l2int_2D
- m_numeric_tools/l2int_3D
- m_numeric_tools/l2norm_rdp
- m_numeric_tools/lfind
- m_numeric_tools/linfit_dpc
- m_numeric_tools/linfit_rdp
- m_numeric_tools/linfit_spc
- m_numeric_tools/linspace
- m_numeric_tools/list2blocks
- m_numeric_tools/llsfit_svd
- m_numeric_tools/mask2blocks
- m_numeric_tools/midpoint_
- m_numeric_tools/mincm
- m_numeric_tools/mkherm
- m_numeric_tools/nderiv
- m_numeric_tools/newrap_step
- m_numeric_tools/pack_matrix
- m_numeric_tools/pade
- m_numeric_tools/pfactorize
- m_numeric_tools/polyn_interp
- m_numeric_tools/print_arr1d_dpc
- m_numeric_tools/print_arr1d_spc
- m_numeric_tools/print_arr2d_dpc
- m_numeric_tools/print_arr2d_spc
- m_numeric_tools/quadrature
- m_numeric_tools/rdp2cdp_2D
- m_numeric_tools/rdp2cdp_3D
- m_numeric_tools/rdp2cdp_4D
- m_numeric_tools/rdp2cdp_5D
- m_numeric_tools/rdp2cdp_6D
- m_numeric_tools/remove_copies
- m_numeric_tools/reverse_int
- m_numeric_tools/reverse_rdp
- m_numeric_tools/rhophi
- m_numeric_tools/simpson
- m_numeric_tools/simpson_cplx
- m_numeric_tools/simpson_int
- m_numeric_tools/smooth
- m_numeric_tools/stats_eval
- m_numeric_tools/stats_t
- m_numeric_tools/symmetrize_dpc
- m_numeric_tools/symmetrize_spc
- m_numeric_tools/trapezoidal_
- m_numeric_tools/uniformrandom
- m_numeric_tools/unit_matrix_cdp
- m_numeric_tools/unit_matrix_int
- m_numeric_tools/unit_matrix_rdp
- m_numeric_tools/vdiff_eval
- m_numeric_tools/vdiff_print
- m_numeric_tools/vdiff_t
- m_numeric_tools/wrap2_pmhalf
- m_numeric_tools/wrap2_zero_one
ABINIT/blocked_loop [ 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 ]
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 ]
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 ]
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 ]
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 ]
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