Line data Source code
1 : !!****m* ABINIT/m_numeric_tools
2 : !! NAME
3 : !! m_numeric_tools
4 : !!
5 : !! FUNCTION
6 : !! This module contains basic tools for numeric computations.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG, GMR, MJV, XG, MVeithen, NH, FJ, MT, DCS, FrD, Olevano, Reining, Sottile, AL)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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 contiguous 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 :: polcoe ! Extract coefficients of polynomial interpolation (Numerical Recipes)
59 : public :: quadrature ! Driver routine for performing quadratures in finite domains using different algorithms
60 : public :: cspint ! Estimates the integral of a tabulated function.
61 : public :: ctrap ! Corrected trapezoidal integral on uniform grid of spacing hh.
62 : public :: coeffs_gausslegint ! Compute the coefficients (supports and weights) for Gauss-Legendre integration.
63 : public :: simpson_cplx ! Integrate a complex function via extended Simpson's rule.
64 : public :: hermitianize ! Force a square matrix to be hermitian
65 : public :: mkherm ! Make the complex array(2,ndim,ndim) hermitian, by adding half of it 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[ ! where 1 is not included (tol12)
83 : public :: wrap2_pmhalf ! Transforms a real number in areduced number in the interval ]-1/2,1/2] ! where -1/2 is not included (tol12)
84 : public :: interpol1d ! Linear interpolation in 1D
85 : public :: interpol1d_c ! Linear interpolation in 1D for complex data
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 to be interpolated in interpol3d
89 : public :: interpolate_denpot ! Liner interpolation of scalar field e.g. density of potential
90 : public :: interpolate_ur ! Liner interpolation of complex wavefunctions (single and double precision version)
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 : public :: geteuler ! Compute the Euler angles corresponding to the spin quantization axis
106 :
107 : !MG FIXME: deprecated: just to avoid updating refs while refactoring.
108 : public :: dotproduct
109 :
110 : interface arth
111 : module procedure arth_int
112 : module procedure arth_rdp
113 : end interface arth
114 :
115 : interface reverse
116 : module procedure reverse_int
117 : module procedure reverse_rdp
118 : end interface reverse
119 :
120 : interface set2unit
121 : module procedure unit_matrix_int
122 : module procedure unit_matrix_rdp
123 : module procedure unit_matrix_cdp
124 : end interface set2unit
125 :
126 : interface get_trace
127 : module procedure get_trace_int
128 : module procedure get_trace_rdp
129 : module procedure get_trace_cdp
130 : end interface get_trace
131 :
132 : interface interpolate_ur
133 : module procedure interpolate_ur_spc
134 : module procedure interpolate_ur_dpc
135 : end interface interpolate_ur
136 :
137 : !interface cart_prod33
138 : ! module procedure cart_prod33_int
139 : ! module procedure cart_prod33_rdp
140 : ! module procedure cart_prod33_cdp
141 : !end interface cart_prod33
142 :
143 : interface get_diag
144 : module procedure get_diag_int
145 : module procedure get_diag_rdp
146 : module procedure get_diag_cdp
147 : end interface get_diag
148 :
149 : interface isdiagmat
150 : module procedure isdiagmat_int
151 : module procedure isdiagmat_rdp
152 : !module procedure isdiagmat_cdp
153 : end interface isdiagmat
154 :
155 : interface inrange
156 : module procedure inrange_int
157 : module procedure inrange_dp
158 : end interface inrange
159 :
160 : interface l2int
161 : module procedure l2int_1D
162 : module procedure l2int_2D
163 : module procedure l2int_3D
164 : end interface l2int
165 :
166 : interface r2c
167 : module procedure rdp2cdp_0D
168 : module procedure rdp2cdp_1D
169 : module procedure rdp2cdp_2D
170 : module procedure rdp2cdp_3D
171 : module procedure rdp2cdp_4D
172 : module procedure rdp2cdp_5D
173 : module procedure rdp2cdp_6D
174 : end interface r2c
175 :
176 : interface c2r
177 : module procedure cdp2rdp_0D
178 : module procedure cdp2rdp_1D
179 : module procedure cdp2rdp_2D
180 : module procedure cdp2rdp_3D
181 : module procedure cdp2rdp_4D
182 : module procedure cdp2rdp_5D
183 : module procedure cdp2rdp_6D
184 : end interface c2r
185 :
186 : interface isinteger
187 : module procedure is_integer_0d
188 : module procedure is_integer_1d
189 : end interface isinteger
190 :
191 : interface is_zero
192 : module procedure is_zero_rdp_0d
193 : module procedure is_zero_rdp_1d
194 : end interface is_zero
195 :
196 : interface bisect
197 : module procedure bisect_rdp
198 : module procedure bisect_int
199 : end interface bisect
200 :
201 : interface imax_loc
202 : module procedure imax_loc_int
203 : module procedure imax_loc_rdp
204 : end interface imax_loc
205 :
206 : interface imin_loc
207 : module procedure imin_loc_int
208 : module procedure imin_loc_rdp
209 : end interface imin_loc
210 :
211 : interface linfit
212 : module procedure linfit_rdp
213 : module procedure linfit_spc
214 : module procedure linfit_dpc
215 : end interface linfit
216 :
217 : interface hermitianize
218 : module procedure hermitianize_spc
219 : module procedure hermitianize_dpc
220 : end interface hermitianize
221 :
222 : interface symmetrize
223 : module procedure symmetrize_spc
224 : module procedure symmetrize_dpc
225 : end interface symmetrize
226 :
227 : interface print_arr !TODO add prtm
228 : module procedure print_arr1d_spc
229 : module procedure print_arr1d_dpc
230 : module procedure print_arr2d_spc
231 : module procedure print_arr2d_dpc
232 : end interface print_arr
233 :
234 : interface operator (.x.)
235 : module procedure cross_product_int
236 : module procedure cross_product_rdp
237 : end interface
238 :
239 : interface l2norm
240 : module procedure l2norm_rdp
241 : end interface l2norm
242 :
243 : interface isordered
244 : module procedure isordered_rdp
245 : end interface isordered
246 : !!***
247 :
248 : !----------------------------------------------------------------------
249 :
250 : !!****t* m_numeric_tools/stats_t
251 : !! NAME
252 : !! stats_t
253 : !!
254 : !! FUNCTION
255 : !! Statistical parameters of a data distribution.
256 : !!
257 : !! SOURCE
258 :
259 : type, public :: stats_t
260 : real(dp) :: mean
261 : real(dp) :: stdev
262 : real(dp) :: min
263 : real(dp) :: max
264 : end type stats_t
265 :
266 : public :: stats_eval ! Calculate statistical parameters of a data distribution.
267 : !!***
268 :
269 : !----------------------------------------------------------------------
270 :
271 : !!****t* m_numeric_tools/vdiff_t
272 : !! NAME
273 : !! vdiff_t
274 : !!
275 : !! FUNCTION
276 : !! Estimate the "distance" between two functions tabulated on a homogeneous grid.
277 : !! Use `vidff` function to construct the object.
278 : !!
279 : !! SOURCE
280 :
281 : type,public :: vdiff_t
282 :
283 : real(dp) :: int_adiff = zero ! \int |f1-f2| dr
284 : real(dp) :: mean_adiff = zero ! Mean {|f1-f2|}
285 : real(dp) :: stdev_adiff = zero ! Standard deviation of {|f1-f2|}
286 : real(dp) :: min_adiff = zero ! Min {|f1-f2|}
287 : real(dp) :: max_adiff = zero ! Max {|f1-f2|}
288 : real(dp) :: l1_rerr = zero ! (\int |f1-f2| dr) / (\int |f2| dr)
289 :
290 : contains
291 :
292 : procedure :: eval => vdiff_eval ! Estimate the "distance" between two functions tabulated on a homogeneous grid.
293 : procedure :: print => vdiff_print ! Print vdiff_t to formatted file.
294 :
295 : end type vdiff_t
296 :
297 : !!***
298 :
299 : !===========================================================
300 : ! A single non-empty bin
301 : !===========================================================
302 : type :: bin_t
303 : integer :: npts = 0
304 : real(dp) :: xmin = 0.0_dp
305 : real(dp) :: xmax = 0.0_dp
306 : integer, allocatable :: idx(:) ! indices into x(:)
307 : !contains
308 : ! procedure :: average => bin_average
309 : ! procedure :: size => bin_size
310 : end type bin_t
311 :
312 : !===========================================================
313 : ! Bins container: only store non-empty bins
314 : !===========================================================
315 : type, public :: bins_t
316 : real(dp), pointer :: xvals(:) => null()
317 :
318 : integer :: total_points = 0
319 : integer :: nbins = 0 ! number of non-empty bins
320 :
321 : type(bin_t), allocatable :: bin(:)
322 :
323 : ! Binning metadata
324 : real(dp) :: xmin = zero
325 : real(dp) :: xmax = zero
326 : real(dp) :: dx = zero
327 : contains
328 : procedure :: init => bins_init
329 : !procedure :: free => bins_free
330 : !procedure :: loop => bins_loop
331 : end type bins_t
332 : !!***
333 :
334 : contains !===========================================================
335 : !!***
336 :
337 : !!****f* m_numeric_tools/arth_int
338 : !! NAME
339 : !! arth_int
340 : !!
341 : !! FUNCTION
342 : !! Returns an array of length nn containing an arithmetic progression whose
343 : !! starting value is start and whose step is step.
344 : !!
345 : !! INPUTS
346 : !! start=initial point
347 : !! step=the increment
348 : !! nn=the number of points
349 : !!
350 : !! OUTPUT
351 : !! arth(nn)=the progression
352 : !!
353 : !! SOURCE
354 :
355 0 : pure function arth_int(start, step, nn)
356 :
357 : !Arguments ------------------------------------
358 : !scalars
359 : integer,intent(in) :: nn
360 : integer,intent(in) :: start,step
361 : integer :: arth_int(nn)
362 :
363 : !Local variables-------------------------------
364 : integer :: ii
365 : ! *********************************************************************
366 :
367 0 : select case (nn)
368 :
369 : case (1:)
370 0 : arth_int(1)=start
371 0 : do ii=2,nn
372 0 : arth_int(ii)=arth_int(ii-1)+step
373 : end do
374 :
375 : case (0)
376 0 : return
377 : end select
378 :
379 : end function arth_int
380 : !!***
381 :
382 : !----------------------------------------------------------------------
383 :
384 : !!****f* m_numeric_tools/arth_rdp
385 : !! NAME
386 : !! arth_rdp
387 : !!
388 : !! FUNCTION
389 : !!
390 : !! INPUTS
391 : !!
392 : !! OUTPUT
393 : !!
394 : !! SOURCE
395 :
396 482590 : pure function arth_rdp(start, step, nn)
397 :
398 : !Arguments ------------------------------------
399 : !scalars
400 : integer,intent(in) :: nn
401 : real(dp),intent(in) :: start,step
402 : real(dp) :: arth_rdp(nn)
403 :
404 : !Local variables-------------------------------
405 : integer :: ii
406 : ! *********************************************************************
407 :
408 482590 : select case (nn)
409 : case (1:)
410 482590 : arth_rdp(1)=start
411 5031754 : do ii=2,nn
412 5031754 : arth_rdp(ii)=arth_rdp(ii-1)+step
413 : end do
414 :
415 : case (0)
416 482590 : return
417 : end select
418 :
419 : end function arth_rdp
420 : !!***
421 :
422 : !----------------------------------------------------------------------
423 :
424 : !!****f* m_numeric_tools/linspace
425 : !! NAME
426 : !! linspace
427 : !!
428 : !! FUNCTION
429 : !!
430 : !! INPUTS
431 : !!
432 : !! OUTPUT
433 : !!
434 : !! SOURCE
435 :
436 12520 : pure function linspace(start, stop, nn)
437 :
438 : !Arguments ------------------------------------
439 : !scalars
440 : integer,intent(in) :: nn
441 : real(dp),intent(in) :: start, stop
442 : real(dp) :: linspace(nn)
443 :
444 : !Local variables-------------------------------
445 : real(dp) :: length
446 : integer :: ii
447 : ! *********************************************************************
448 :
449 12520 : select case (nn)
450 : case (1:)
451 12520 : length = stop - start
452 31953847 : do ii=1,nn
453 31953847 : linspace(ii) = start + length * (ii-1) / (nn-1)
454 : end do
455 :
456 : case (0)
457 12520 : return
458 : end select
459 :
460 : end function linspace
461 : !!***
462 :
463 :
464 : !----------------------------------------------------------------------
465 :
466 : !!****f* m_numeric_tools/geop
467 : !! NAME
468 : !! geop
469 : !!
470 : !! FUNCTION
471 : !! Returns an array of length nn containing a geometric progression whose
472 : !! starting value is start and whose factor is factor!
473 : !!
474 : !! INPUTS
475 : !! start=initial point
476 : !! factor=the factor of the geometric progression
477 : !! nn=the number of points
478 : !!
479 : !! OUTPUT
480 : !! geop(nn)=the progression
481 : !!
482 : !! SOURCE
483 :
484 :
485 0 : pure function geop(start,factor,nn) result(res)
486 :
487 : !Arguments ------------------------------------
488 : !scalars
489 : real(dp),intent(in) :: start,factor
490 : integer,intent(in) :: nn
491 : real(dp) :: res(nn)
492 :
493 : !Local variables-------------------------------
494 : integer :: ii
495 : ! *********************************************************************
496 :
497 0 : if (nn>0) res(1)=start
498 0 : do ii=2,nn
499 0 : res(ii)=res(ii-1)*factor
500 : end do
501 :
502 0 : end function geop
503 : !!***
504 :
505 : !----------------------------------------------------------------------
506 :
507 : !!****f* m_numeric_tools/reverse_int
508 : !! NAME
509 : !! reverse_int
510 : !!
511 : !! FUNCTION
512 : !! Reverse a 1D array *IN PLACE*. Target: INT arrays
513 : !!
514 : !! SOURCE
515 :
516 0 : subroutine reverse_int(arr)
517 :
518 : !Arguments ------------------------------------
519 : !scalars
520 : integer,intent(inout) :: arr(:)
521 : !arrays
522 : integer :: ii,nn,swap
523 : ! *************************************************************************
524 :
525 0 : nn = SIZE(arr)
526 0 : if (nn <= 1) return
527 :
528 0 : do ii=1,nn/2
529 0 : swap = arr(ii)
530 0 : arr(ii) = arr(nn-ii+1)
531 0 : arr(nn-ii+1) = swap
532 : end do
533 :
534 : end subroutine reverse_int
535 : !!***
536 :
537 : !----------------------------------------------------------------------
538 :
539 : !!****f* m_numeric_tools/reverse_rdp
540 : !! NAME
541 : !! reverse_rdp
542 : !!
543 : !! FUNCTION
544 : !! Reverse a 1D array *IN PLACE*. Target: DP arrays
545 : !!
546 : !! SOURCE
547 :
548 0 : subroutine reverse_rdp(arr)
549 :
550 : !Arguments ------------------------------------
551 : !scalars
552 : real(dp),intent(inout) :: arr(:)
553 : !arrays
554 : integer :: ii,nn
555 : real(dp) :: swap
556 : ! *************************************************************************
557 :
558 0 : nn = SIZE(arr)
559 0 : if (nn <= 1) return
560 :
561 0 : do ii=1,nn/2
562 0 : swap = arr(ii)
563 0 : arr(ii) = arr(nn-ii+1)
564 0 : arr(nn-ii+1) = swap
565 : end do
566 :
567 : end subroutine reverse_rdp
568 : !!***
569 :
570 : !----------------------------------------------------------------------
571 :
572 : !!****f* m_numeric_tools/unit_matrix_int
573 : !! NAME
574 : !! unit_matrix_int
575 : !!
576 : !! FUNCTION
577 : !! Set the matrix matrix to be a unit matrix (if it is square).
578 : !!
579 : !! SIDE EFFECTS
580 : !! matrix(:,:)=set to unit on exit
581 : !!
582 : !! SOURCE
583 :
584 17 : pure subroutine unit_matrix_int(matrix)
585 :
586 : !Arguments ------------------------------------
587 : integer,intent(inout) :: matrix(:,:)
588 :
589 : !Local variables-------------------------------
590 : !scalars
591 : integer :: ii,nn
592 : ! *********************************************************************
593 :
594 17 : nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
595 221 : matrix(:,:)=0
596 68 : do ii=1,nn
597 68 : matrix(ii,ii)=1
598 : end do
599 :
600 17 : end subroutine unit_matrix_int
601 : !!***
602 :
603 : !----------------------------------------------------------------------
604 :
605 : !!****f* m_numeric_tools/unit_matrix_rdp
606 : !! NAME
607 : !! unit_matrix_rdp
608 : !!
609 : !! FUNCTION
610 : !! Set the matrix matrix to be a unit matrix (if it is square).
611 : !!
612 : !! SIDE EFFECTS
613 : !! matrix(:,:)=set to unit on exit
614 : !!
615 : !! SOURCE
616 :
617 0 : pure subroutine unit_matrix_rdp(matrix)
618 :
619 : !Arguments ------------------------------------
620 : real(dp),intent(inout) :: matrix(:,:)
621 :
622 : !Local variables-------------------------------
623 : integer :: ii,nn
624 : ! *********************************************************************
625 :
626 0 : nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
627 0 : matrix(:,:)=zero
628 0 : do ii=1,nn
629 0 : matrix(ii,ii)=one
630 : end do
631 :
632 0 : end subroutine unit_matrix_rdp
633 : !!***
634 :
635 : !----------------------------------------------------------------------
636 :
637 : !!****f* m_numeric_tools/unit_matrix_cdp
638 : !! NAME
639 : !! unit_matrix_cdp
640 : !!
641 : !! FUNCTION
642 : !! Set the matrix matrix to be a unit matrix (if it is square).
643 : !!
644 : !! SIDE EFFECTS
645 : !! matrix(:,:)=set to unit on exit
646 : !!
647 : !! SOURCE
648 :
649 0 : pure subroutine unit_matrix_cdp(matrix)
650 :
651 : !Arguments ------------------------------------
652 : complex(dp),intent(inout) :: matrix(:,:)
653 :
654 : !Local variables-------------------------------
655 : integer :: ii,nn
656 : ! *********************************************************************
657 :
658 0 : nn=MIN(SIZE(matrix,DIM=1),SIZE(matrix,DIM=2))
659 0 : matrix=czero
660 0 : do ii=1,nn
661 0 : matrix(ii,ii)=cone
662 : end do
663 :
664 0 : end subroutine unit_matrix_cdp
665 : !!***
666 :
667 : !----------------------------------------------------------------------
668 :
669 : !!****f* m_numeric_tools/get_trace_int
670 : !! NAME
671 : !! get_trace_int
672 : !!
673 : !! FUNCTION
674 : !! Calculate the trace of a square matrix
675 : !!
676 : !! INPUTS
677 : !! matrix(:,:)
678 : !!
679 : !! OUTPUT
680 : !! trace=the trace
681 : !!
682 : !! SOURCE
683 :
684 0 : integer pure function get_trace_int(matrix) result(trace)
685 :
686 : !Arguments ------------------------------------
687 : integer,intent(in) :: matrix(:,:)
688 :
689 : !Local variables-------------------------------
690 : integer :: ii
691 : ! *********************************************************************
692 :
693 0 : trace=0
694 0 : do ii=1,size(matrix,dim=1)
695 0 : trace=trace+matrix(ii,ii)
696 : end do
697 :
698 0 : end function get_trace_int
699 : !!***
700 :
701 : !----------------------------------------------------------------------
702 :
703 : !!****f* m_numeric_tools/get_trace_rdp
704 : !! NAME
705 : !! get_trace_int
706 : !!
707 : !! FUNCTION
708 : !! Calculate the trace of a square matrix (real(dp) version)
709 : !!
710 : !! INPUTS
711 : !! matrix(:,:)
712 : !!
713 : !! OUTPUT
714 : !! trace=the trace
715 : !!
716 : !! SOURCE
717 :
718 4947 : real(dp) pure function get_trace_rdp(matrix) result(trace)
719 :
720 : !Arguments ------------------------------------
721 : real(dp),intent(in) :: matrix(:,:)
722 :
723 : !Local variables-------------------------------
724 : integer :: ii
725 : ! *********************************************************************
726 :
727 4947 : trace=zero
728 19788 : do ii=1,size(matrix,dim=1)
729 19788 : trace=trace+matrix(ii,ii)
730 : end do
731 :
732 4947 : end function get_trace_rdp
733 : !!***
734 :
735 : !----------------------------------------------------------------------
736 :
737 : !!****f* m_numeric_tools/get_trace_cdp
738 : !! NAME
739 : !! get_trace_cdp
740 : !!
741 : !! FUNCTION
742 : !! Calculate the trace of a square matrix (complex(dp) version)
743 : !!
744 : !! INPUTS
745 : !!
746 : !! OUTPUT
747 : !!
748 : !! SOURCE
749 :
750 0 : complex(dp) pure function get_trace_cdp(matrix) result(trace)
751 :
752 : !Arguments ------------------------------------
753 : complex(dp),intent(in) :: matrix(:,:)
754 :
755 : !Local variables-------------------------------
756 : integer :: ii
757 : ! *********************************************************************
758 :
759 0 : trace=czero
760 0 : do ii=1,size(matrix,dim=1)
761 0 : trace=trace+matrix(ii,ii)
762 : end do
763 :
764 0 : end function get_trace_cdp
765 : !!***
766 :
767 : !!****f* m_numeric_tools/get_diag_int
768 : !! NAME
769 : !! get_diag_int
770 : !!
771 : !! FUNCTION
772 : !! Return the diagonal of a square matrix as a vector
773 : !!
774 : !! INPUTS
775 : !! matrix(:,:)
776 : !!
777 : !! OUTPUT
778 : !! diag(:)=the diagonal
779 : !!
780 : !! SOURCE
781 :
782 43 : function get_diag_int(mat) result(diag)
783 :
784 : !Arguments ------------------------------------
785 : !scalars
786 : integer,intent(in) :: mat(:,:)
787 : integer :: diag(SIZE(mat,1))
788 :
789 : !Local variables-------------------------------
790 : integer :: ii
791 : ! *************************************************************************
792 :
793 43 : ii = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
794 :
795 172 : do ii=1,SIZE(mat,1)
796 172 : diag(ii)=mat(ii,ii)
797 : end do
798 :
799 43 : end function get_diag_int
800 : !!***
801 :
802 : !----------------------------------------------------------------------
803 :
804 : !!****f* m_numeric_tools/get_diag_rdp
805 : !! NAME
806 : !! get_diag_rdp
807 : !!
808 : !! FUNCTION
809 : !! Return the diagonal of a square matrix as a vector
810 : !!
811 : !! INPUTS
812 : !! matrix(:,:)
813 : !!
814 : !! OUTPUT
815 : !! diag(:)=the diagonal
816 : !!
817 : !! SOURCE
818 :
819 0 : function get_diag_rdp(mat) result(diag)
820 :
821 : !Arguments ------------------------------------
822 : !scalars
823 : real(dp),intent(in) :: mat(:,:)
824 : real(dp) :: diag(SIZE(mat,1))
825 :
826 : !Local variables-------------------------------
827 : integer :: ii
828 : ! *************************************************************************
829 :
830 0 : ABI_CHECK(SIZE(mat,1) == SIZE(mat,2), 'Matrix not square')
831 :
832 0 : do ii=1,SIZE(mat,1)
833 0 : diag(ii) = mat(ii,ii)
834 : end do
835 :
836 0 : end function get_diag_rdp
837 : !!***
838 :
839 : !----------------------------------------------------------------------
840 :
841 : !!****f* m_numeric_tools/get_diag_cdp
842 : !! NAME
843 : !! get_diag_cdp
844 : !!
845 : !! FUNCTION
846 : !!
847 : !! INPUTS
848 : !!
849 : !! OUTPUT
850 : !!
851 : !! SOURCE
852 :
853 0 : function get_diag_cdp(cmat) result(cdiag)
854 :
855 : !Arguments ------------------------------------
856 : !scalars
857 : complex(dp),intent(in) :: cmat(:,:)
858 : complex(dp) :: cdiag(SIZE(cmat,1))
859 :
860 : !Local variables-------------------------------
861 : integer :: ii
862 : ! *************************************************************************
863 :
864 0 : ABI_CHECK(SIZE(cmat,1) == SIZE(cmat,2), 'Matrix not square')
865 :
866 0 : do ii=1,SIZE(cmat,1)
867 0 : cdiag(ii)=cmat(ii,ii)
868 : end do
869 :
870 0 : end function get_diag_cdp
871 : !!***
872 :
873 : !----------------------------------------------------------------------
874 :
875 : !!****f* m_numeric_tools/isdiagmat_int
876 : !! NAME
877 : !! isdiagmat_int
878 : !!
879 : !! FUNCTION
880 : !! True if matrix mat is diagonal
881 : !!
882 : !! SOURCE
883 :
884 11460 : pure logical function isdiagmat_int(mat) result(ans)
885 :
886 : !Arguments ------------------------------------
887 : !scalars
888 : integer,intent(in) :: mat(:,:)
889 :
890 : !Local variables-------------------------------
891 : integer :: ii,jj
892 : ! *************************************************************************
893 :
894 11460 : ans = .True.
895 37935 : do jj=1,size(mat,dim=2)
896 120132 : do ii=1,size(mat,dim=1)
897 84843 : if (ii == jj) cycle
898 82203 : if (mat(ii,jj) /= 0) then
899 2646 : ans = .False.; return
900 : end if
901 : end do
902 : end do
903 :
904 8814 : end function isdiagmat_int
905 : !!***
906 :
907 : !----------------------------------------------------------------------
908 :
909 : !!****f* m_numeric_tools/isdiagmat_rdp
910 : !! NAME
911 : !! isdiagmat_rdp
912 : !!
913 : !! FUNCTION
914 : !! True if matrix mat is diagonal within the given absolute tolerance (default: tol12)
915 : !!
916 : !! SOURCE
917 :
918 0 : pure logical function isdiagmat_rdp(mat, atol) result(ans)
919 :
920 : !Arguments ------------------------------------
921 : !scalars
922 : real(dp),intent(in) :: mat(:,:)
923 : real(dp),optional,intent(in) :: atol
924 :
925 : !Local variables-------------------------------
926 : integer :: ii,jj
927 : real(dp) :: my_atol
928 : ! *************************************************************************
929 :
930 0 : my_atol = tol12; if (present(atol)) my_atol = atol
931 :
932 0 : ans = .True.
933 0 : do jj=1,size(mat,dim=2)
934 0 : do ii=1,size(mat,dim=1)
935 0 : if (ii == jj) cycle
936 0 : if (abs(mat(ii,jj)) > my_atol) then
937 0 : ans = .False.; return
938 : end if
939 : end do
940 : end do
941 :
942 0 : end function isdiagmat_rdp
943 : !!***
944 :
945 : !----------------------------------------------------------------------
946 :
947 : !!****f* m_numeric_tools/l2int_1D
948 : !! NAME
949 : !! l2int_1D
950 : !!
951 : !! FUNCTION
952 : !! Convert a logical array into an int array (True --> 1, False --> 0)
953 : !!
954 : !! INPUTS
955 : !! larr(:)=the input logical array
956 : !!
957 : !! SOURCE
958 :
959 2 : pure function l2int_1D(larr) result(int_arr)
960 :
961 : !Arguments ------------------------------------
962 : !scalars
963 : logical,intent(in) :: larr(:)
964 : integer :: int_arr(size(larr))
965 :
966 : ! *********************************************************************
967 :
968 10 : where (larr)
969 : int_arr = 1
970 : elsewhere
971 : int_arr = 0
972 : end where
973 :
974 2 : end function l2int_1D
975 : !!***
976 :
977 : !----------------------------------------------------------------------
978 :
979 : !!****f* m_numeric_tools/l2int_2D
980 : !! NAME
981 : !! l2int_2D
982 : !!
983 : !! FUNCTION
984 : !! Convert a logical array into an int array (True --> 1, False --> 0)
985 : !!
986 : !! INPUTS
987 : !! larr(:)=the input logical array
988 : !!
989 : !! SOURCE
990 :
991 15 : pure function l2int_2D(larr) result(int_arr)
992 :
993 : !Arguments ------------------------------------
994 : !scalars
995 : logical,intent(in) :: larr(:,:)
996 : integer :: int_arr(size(larr,1), size(larr,2))
997 : ! *********************************************************************
998 :
999 175 : where (larr)
1000 : int_arr = 1
1001 : elsewhere
1002 : int_arr = 0
1003 : end where
1004 :
1005 15 : end function l2int_2D
1006 : !!***
1007 :
1008 : !----------------------------------------------------------------------
1009 :
1010 : !!****f* m_numeric_tools/l2int_3D
1011 : !! NAME
1012 : !! l2int_3D
1013 : !!
1014 : !! FUNCTION
1015 : !! Convert a logical array into an int array (True --> 1, False --> 0)
1016 : !!
1017 : !! INPUTS
1018 : !! larr(:)=the input logical array
1019 : !!
1020 : !! SOURCE
1021 :
1022 15 : pure function l2int_3D(larr) result(int_arr)
1023 :
1024 : !Arguments ------------------------------------
1025 : !scalars
1026 : logical,intent(in) :: larr(:,:,:)
1027 : integer :: int_arr(size(larr,1), size(larr,2), size(larr,3))
1028 : ! *********************************************************************
1029 :
1030 7231 : where (larr)
1031 : int_arr = 1
1032 : elsewhere
1033 : int_arr = 0
1034 : end where
1035 :
1036 15 : end function l2int_3D
1037 : !!***
1038 :
1039 : !----------------------------------------------------------------------
1040 :
1041 : !!***!!****f* m_numeric_tools/rdp2cdp_0D
1042 : !! NAME
1043 : !! rdp2cdp_0D
1044 : !!
1045 : !! FUNCTION
1046 : !! Create a complex number starting from a real array containing real and imaginary part
1047 : !!
1048 : !! INPUTS
1049 : !! rr(:)=the real array
1050 : !!
1051 : !! OUTPUT
1052 : !! cc(:)=the complex array
1053 : !!
1054 : !! SOURCE
1055 :
1056 0 : pure function rdp2cdp_0D(rr) result(cc)
1057 :
1058 : !Arguments ------------------------------------
1059 : !scalars
1060 : real(dp),intent(in) :: rr(2)
1061 : complex(dp) :: cc
1062 : ! *********************************************************************
1063 :
1064 0 : cc = CMPLX(rr(1), rr(2), kind=dp)
1065 :
1066 0 : end function rdp2cdp_0D
1067 : !!***
1068 :
1069 : !----------------------------------------------------------------------
1070 :
1071 : !!***!!****f* m_numeric_tools/rdp2cdp_1D
1072 : !! NAME
1073 : !! rdp2cdp_1D
1074 : !!
1075 : !! FUNCTION
1076 : !! Create a complex array starting from a real array containing real and imaginary part
1077 : !!
1078 : !! INPUTS
1079 : !! rr(:)=the real array
1080 : !!
1081 : !! OUTPUT
1082 : !! cc(:)=the complex array
1083 : !!
1084 : !! SOURCE
1085 :
1086 0 : pure function rdp2cdp_1D(rr) result(cc)
1087 :
1088 : !Arguments ------------------------------------
1089 : !scalars
1090 : real(dp),intent(in) :: rr(:,:)
1091 : complex(dp) :: cc(SIZE(rr,2))
1092 : ! *********************************************************************
1093 :
1094 0 : cc(:)=CMPLX(rr(1,:),rr(2,:),kind=dp)
1095 :
1096 0 : end function rdp2cdp_1D
1097 : !!***
1098 :
1099 : !----------------------------------------------------------------------
1100 :
1101 : !!****f* m_numeric_tools/rdp2cdp_2D
1102 : !! NAME
1103 : !! rdp2cdp_2D
1104 : !!
1105 : !! FUNCTION
1106 : !!
1107 : !! INPUTS
1108 : !!
1109 : !! OUTPUT
1110 : !!
1111 : !! SOURCE
1112 :
1113 0 : pure function rdp2cdp_2D(rr) result(cc)
1114 :
1115 : !Arguments ------------------------------------
1116 : !scalars
1117 : real(dp),intent(in) :: rr(:,:,:)
1118 : complex(dp) :: cc(SIZE(rr,2),SIZE(rr,3))
1119 : ! *********************************************************************
1120 :
1121 0 : cc(:,:)=CMPLX(rr(1,:,:),rr(2,:,:), kind=dp)
1122 :
1123 0 : end function rdp2cdp_2D
1124 : !!***
1125 :
1126 : !----------------------------------------------------------------------
1127 :
1128 : !!****f* m_numeric_tools/rdp2cdp_3D
1129 : !! NAME
1130 : !! rdp2cdp_3D
1131 : !!
1132 : !! FUNCTION
1133 : !!
1134 : !! INPUTS
1135 : !!
1136 : !! OUTPUT
1137 : !!
1138 : !! SOURCE
1139 :
1140 0 : pure function rdp2cdp_3D(rr) result(cc)
1141 :
1142 : !Arguments ------------------------------------
1143 : !scalars
1144 : real(dp),intent(in) :: rr(:,:,:,:)
1145 : complex(dp) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4))
1146 : ! *********************************************************************
1147 :
1148 0 : cc(:,:,:)=CMPLX(rr(1,:,:,:),rr(2,:,:,:), kind=dp)
1149 :
1150 0 : end function rdp2cdp_3D
1151 : !!***
1152 :
1153 : !----------------------------------------------------------------------
1154 :
1155 : !!****f* m_numeric_tools/rdp2cdp_4D
1156 : !! NAME
1157 : !! rdp2cdp_4D
1158 : !!
1159 : !! FUNCTION
1160 : !!
1161 : !! INPUTS
1162 : !!
1163 : !! OUTPUT
1164 : !!
1165 : !! SOURCE
1166 :
1167 0 : pure function rdp2cdp_4D(rr) result(cc)
1168 :
1169 : !Arguments ------------------------------------
1170 : !scalars
1171 : real(dp),intent(in) :: rr(:,:,:,:,:)
1172 : complex(dp) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5))
1173 : ! *********************************************************************
1174 :
1175 0 : cc(:,:,:,:)=CMPLX(rr(1,:,:,:,:),rr(2,:,:,:,:), kind=dp)
1176 :
1177 0 : end function rdp2cdp_4D
1178 : !!***
1179 :
1180 : !----------------------------------------------------------------------
1181 :
1182 : !!****f* m_numeric_tools/rdp2cdp_5D
1183 : !! NAME
1184 : !! rdp2cdp_5D
1185 : !!
1186 : !! FUNCTION
1187 : !!
1188 : !! INPUTS
1189 : !!
1190 : !! OUTPUT
1191 : !!
1192 : !! SOURCE
1193 :
1194 0 : pure function rdp2cdp_5D(rr) result(cc)
1195 :
1196 : !Arguments ------------------------------------
1197 : !scalars
1198 : real(dp),intent(in) :: rr(:,:,:,:,:,:)
1199 : complex(dp) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6))
1200 : ! *********************************************************************
1201 :
1202 0 : cc(:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:),rr(2,:,:,:,:,:), kind=dp)
1203 :
1204 0 : end function rdp2cdp_5D
1205 : !!***
1206 :
1207 : !----------------------------------------------------------------------
1208 :
1209 : !!****f* m_numeric_tools/rdp2cdp_6D
1210 : !! NAME
1211 : !! rdp2cdp_6D
1212 : !!
1213 : !! FUNCTION
1214 : !!
1215 : !! INPUTS
1216 : !!
1217 : !! OUTPUT
1218 : !!
1219 : !! SOURCE
1220 :
1221 0 : pure function rdp2cdp_6D(rr) result(cc)
1222 :
1223 : !Arguments ------------------------------------
1224 : !scalars
1225 : real(dp),intent(in) :: rr(:,:,:,:,:,:,:)
1226 : complex(dp) :: cc(SIZE(rr,2),SIZE(rr,3),SIZE(rr,4),SIZE(rr,5),SIZE(rr,6),SIZE(rr,7))
1227 : ! *********************************************************************
1228 :
1229 0 : cc(:,:,:,:,:,:)=CMPLX(rr(1,:,:,:,:,:,:),rr(2,:,:,:,:,:,:), kind=dp)
1230 :
1231 0 : end function rdp2cdp_6D
1232 : !!***
1233 :
1234 :
1235 : !----------------------------------------------------------------------
1236 :
1237 : !!****f* m_numeric_tools/cdp2rdp_0D
1238 : !! NAME
1239 : !! cdp2rdp_0D
1240 : !!
1241 : !! FUNCTION
1242 : !! Create a real variable containing real and imaginary part starting from a complex array
1243 : !!
1244 : !! INPUTS
1245 : !! cc=the input complex number
1246 : !!
1247 : !! OUTPUT
1248 : !! rr(2=the real array
1249 : !!
1250 : !! SOURCE
1251 :
1252 0 : pure function cdp2rdp_0D(cc) result(rr)
1253 :
1254 : !Arguments ------------------------------------
1255 : !scalars
1256 : complex(dp),intent(in) :: cc
1257 : real(dp) :: rr(2)
1258 : ! *********************************************************************
1259 :
1260 0 : rr(1)=REAL (cc)
1261 0 : rr(2)=AIMAG(cc)
1262 :
1263 0 : end function cdp2rdp_0D
1264 : !!***
1265 :
1266 : !----------------------------------------------------------------------
1267 :
1268 : !!****f* m_numeric_tools/cdp2rdp_1D
1269 : !! NAME
1270 : !! cdp2rdp_1D
1271 : !!
1272 : !! FUNCTION
1273 : !! Create a real array containing real and imaginary part starting from a complex array
1274 : !!
1275 : !! INPUTS
1276 : !! cc(:)=the input complex array
1277 : !!
1278 : !! OUTPUT
1279 : !! rr(2,:)=the real array
1280 : !!
1281 : !! SOURCE
1282 :
1283 291 : pure function cdp2rdp_1D(cc) result(rr)
1284 :
1285 : !Arguments ------------------------------------
1286 : complex(dp),intent(in) :: cc(:)
1287 : real(dp) :: rr(2,SIZE(cc))
1288 : ! *********************************************************************
1289 :
1290 47958 : rr(1,:)=REAL (cc(:))
1291 47958 : rr(2,:)=AIMAG(cc(:))
1292 :
1293 291 : end function cdp2rdp_1D
1294 : !!***
1295 :
1296 : !----------------------------------------------------------------------
1297 :
1298 : !!****f* m_numeric_tools/cdp2rdp_2D
1299 : !! NAME
1300 : !! cdp2rdp_2D
1301 : !!
1302 : !! FUNCTION
1303 : !! Create a real array containing real and imaginary part starting from a complex array
1304 : !!
1305 : !! SOURCE
1306 :
1307 1480 : pure function cdp2rdp_2D(cc) result(rr)
1308 :
1309 : !Arguments ------------------------------------
1310 : complex(dp),intent(in) :: cc(:,:)
1311 : real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2))
1312 : ! *********************************************************************
1313 :
1314 29464 : rr(1,:,:)=REAL (cc(:,:))
1315 29464 : rr(2,:,:)=AIMAG(cc(:,:))
1316 :
1317 1480 : end function cdp2rdp_2D
1318 : !!***
1319 :
1320 : !----------------------------------------------------------------------
1321 :
1322 : !!****f* m_numeric_tools/cdp2rdp_3D
1323 : !! NAME
1324 : !! cdp2rdp_3D
1325 : !!
1326 : !! FUNCTION
1327 : !! Create a real array containing real and imaginary part starting from a complex array
1328 : !!
1329 : !! SOURCE
1330 :
1331 1097 : pure function cdp2rdp_3D(cc) result(rr)
1332 :
1333 : !Arguments ------------------------------------
1334 : complex(dp),intent(in) :: cc(:,:,:)
1335 : real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3))
1336 : ! *********************************************************************
1337 :
1338 196762 : rr(1,:,:,:)=REAL (cc(:,:,:))
1339 196762 : rr(2,:,:,:)=AIMAG(cc(:,:,:))
1340 :
1341 1097 : end function cdp2rdp_3D
1342 : !!***
1343 :
1344 : !----------------------------------------------------------------------
1345 :
1346 : !!****f* m_numeric_tools/cdp2rdp_4D
1347 : !! NAME
1348 : !! cdp2rdp_4D
1349 : !!
1350 : !! FUNCTION
1351 : !! Create a real array containing real and imaginary part starting from a complex array
1352 : !!
1353 : !! SOURCE
1354 :
1355 832 : pure function cdp2rdp_4D(cc) result(rr)
1356 :
1357 : !Arguments ------------------------------------
1358 : complex(dp),intent(in) :: cc(:,:,:,:)
1359 : real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4))
1360 : ! *********************************************************************
1361 :
1362 708274 : rr(1,:,:,:,:)=REAL (cc(:,:,:,:))
1363 708274 : rr(2,:,:,:,:)=AIMAG(cc(:,:,:,:))
1364 :
1365 832 : end function cdp2rdp_4D
1366 : !!***
1367 :
1368 : !----------------------------------------------------------------------
1369 :
1370 : !!****f* m_numeric_tools/cdp2rdp_5D
1371 : !! NAME
1372 : !! cdp2rdp_5D
1373 : !!
1374 : !! FUNCTION
1375 : !! Create a real array containing real and imaginary part starting from a complex array
1376 : !!
1377 : !! SOURCE
1378 :
1379 2 : pure function cdp2rdp_5D(cc) result(rr)
1380 :
1381 : !Arguments ------------------------------------
1382 : complex(dp),intent(in) :: cc(:,:,:,:,:)
1383 : real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5))
1384 : ! *********************************************************************
1385 :
1386 5496 : rr(1,:,:,:,:,:)=REAL (cc(:,:,:,:,:))
1387 5496 : rr(2,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:))
1388 :
1389 2 : end function cdp2rdp_5D
1390 : !!***
1391 :
1392 : !!****f* m_numeric_tools/cdp2rdp_6D
1393 : !! NAME
1394 : !! cdp2rdp_6D
1395 : !!
1396 : !! FUNCTION
1397 : !! Create a real array containing real and imaginary part starting from a complex array
1398 : !!
1399 : !! SOURCE
1400 :
1401 0 : pure function cdp2rdp_6D(cc) result(rr)
1402 :
1403 : !Arguments ------------------------------------
1404 : complex(dp),intent(in) :: cc(:,:,:,:,:,:)
1405 : real(dp) :: rr(2,SIZE(cc,1),SIZE(cc,2),SIZE(cc,3),SIZE(cc,4),SIZE(cc,5),SIZE(cc,6))
1406 : ! *********************************************************************
1407 :
1408 0 : rr(1,:,:,:,:,:,:)=REAL (cc(:,:,:,:,:,:))
1409 0 : rr(2,:,:,:,:,:,:)=AIMAG(cc(:,:,:,:,:,:))
1410 :
1411 0 : end function cdp2rdp_6D
1412 : !!***
1413 :
1414 : !----------------------------------------------------------------------
1415 :
1416 : !!****f* m_numeric_tools/iseven
1417 : !! NAME
1418 : !! iseven
1419 : !!
1420 : !! FUNCTION
1421 : !! Return .TRUE. if the given integer is even
1422 : !!
1423 : !! SOURCE
1424 :
1425 27 : elemental function iseven(nn)
1426 :
1427 : !Arguments ------------------------------------
1428 : !scalars
1429 : integer,intent(in) :: nn
1430 : logical :: iseven
1431 : ! *********************************************************************
1432 :
1433 27 : iseven = ((nn / 2) * 2 == nn)
1434 :
1435 27 : end function iseven
1436 : !!***
1437 :
1438 : !----------------------------------------------------------------------
1439 :
1440 : !!****f* m_numeric_tools/is_integer_0D
1441 : !! NAME
1442 : !! is_integer_0D
1443 : !!
1444 : !! FUNCTION
1445 : !! Return .TRUE. if all elements differ from an integer by less that tol
1446 : !!
1447 : !! INPUTS
1448 : !! rr=the set of real values to be checked
1449 : !! tol=tolerance on the difference between real and integer
1450 : !!
1451 : !! SOURCE
1452 :
1453 223503976 : pure function is_integer_0d(rr,tol) result(ans)
1454 :
1455 : !Arguments ------------------------------------
1456 : !scalars
1457 : real(dp),intent(in) :: tol
1458 : logical :: ans
1459 : !arrays
1460 : real(dp),intent(in) :: rr
1461 : ! *************************************************************************
1462 :
1463 223503976 : ans=(ABS(rr-NINT(rr))<tol)
1464 :
1465 223503976 : end function is_integer_0d
1466 : !!***
1467 :
1468 : !----------------------------------------------------------------------
1469 :
1470 : !!****f* m_numeric_tools/is_integer_1D
1471 : !! NAME
1472 : !! is_integer_1D
1473 : !!
1474 : !! FUNCTION
1475 : !!
1476 : !! INPUTS
1477 : !!
1478 : !! OUTPUT
1479 : !!
1480 : !! SOURCE
1481 :
1482 467114045 : pure function is_integer_1d(rr,tol) result(ans)
1483 :
1484 : !Arguments ------------------------------------
1485 : !scalars
1486 : real(dp),intent(in) :: tol
1487 : logical :: ans
1488 : !arrays
1489 : real(dp),intent(in) :: rr(:)
1490 : ! *************************************************************************
1491 :
1492 541859825 : ans=ALL((ABS(rr-NINT(rr))<tol))
1493 :
1494 467114045 : end function is_integer_1d
1495 : !!***
1496 :
1497 : !----------------------------------------------------------------------
1498 :
1499 : !!****f* m_numeric_tools/is_zero_rdp_0D
1500 : !! NAME
1501 : !! is_zero_rdp_0D
1502 : !!
1503 : !! FUNCTION
1504 : !! Return .TRUE. if all elements differ from zero by less that tol
1505 : !!
1506 : !! INPUTS
1507 : !! rr=the set of real values to be checked
1508 : !! tol=tolerance
1509 : !!
1510 : !! OUTPUT
1511 : !!
1512 : !! SOURCE
1513 :
1514 0 : function is_zero_rdp_0d(rr,tol) result(ans)
1515 :
1516 : !Arguments ------------------------------------
1517 : !scalars
1518 : real(dp),intent(in) :: tol
1519 : logical :: ans
1520 : !arrays
1521 : real(dp),intent(in) :: rr
1522 : ! *************************************************************************
1523 :
1524 0 : ans=(ABS(rr)<tol)
1525 :
1526 0 : end function is_zero_rdp_0d
1527 : !!***
1528 :
1529 : !----------------------------------------------------------------------
1530 :
1531 : !!****f* m_numeric_tools/is_zero_rdp_1d
1532 : !! NAME
1533 : !! is_zero_rdp_1d
1534 : !!
1535 : !! FUNCTION
1536 : !!
1537 : !! INPUTS
1538 : !!
1539 : !! OUTPUT
1540 : !!
1541 : !! SOURCE
1542 :
1543 12200 : function is_zero_rdp_1d(rr,tol) result(ans)
1544 :
1545 : !Arguments ------------------------------------
1546 : !scalars
1547 : real(dp),intent(in) :: tol
1548 : logical :: ans
1549 : !arrays
1550 : real(dp),intent(in) :: rr(:)
1551 : ! *************************************************************************
1552 :
1553 12780 : ans=ALL(ABS(rr(:))<tol)
1554 :
1555 12200 : end function is_zero_rdp_1d
1556 : !!***
1557 :
1558 : !----------------------------------------------------------------------
1559 :
1560 : !!****f* m_numeric_tools/inrange_int
1561 : !! NAME
1562 : !! inrange_int
1563 : !!
1564 : !! FUNCTION
1565 : !! True if int `xval` is inside the interval [win(1), win(2)]
1566 : !!
1567 : !! SOURCE
1568 :
1569 44 : pure logical function inrange_int(xval, win)
1570 :
1571 : !Arguments ------------------------------------
1572 : !scalars
1573 : integer,intent(in) :: xval,win(2)
1574 : ! *************************************************************************
1575 :
1576 44 : inrange_int = (xval >= win(1) .and. xval <= win(2))
1577 :
1578 44 : end function inrange_int
1579 : !!***
1580 :
1581 : !----------------------------------------------------------------------
1582 :
1583 : !!****f* m_numeric_tools/inrange_dp
1584 : !! NAME
1585 : !! inrange_dp
1586 : !!
1587 : !! FUNCTION
1588 : !! True if float `xval` is inside the interval [win(1), win(2)]
1589 : !!
1590 : !! SOURCE
1591 :
1592 0 : pure logical function inrange_dp(xval, win)
1593 :
1594 : !Arguments ------------------------------------
1595 : !scalars
1596 : real(dp),intent(in) :: xval, win(2)
1597 : ! *************************************************************************
1598 :
1599 0 : inrange_dp = (xval >= win(1) .and. xval <= win(2))
1600 :
1601 0 : end function inrange_dp
1602 : !!***
1603 :
1604 : !!****f* m_numeric_tools/bisect_rdp
1605 : !! NAME
1606 : !! bisect_rdp
1607 : !!
1608 : !! FUNCTION
1609 : !! Given an array AA(1:N), and a value x, returns the index j such that AA(j) <= x <= AA(j + 1).
1610 : !! AA must be monotonic, either increasing or decreasing. j=0 or
1611 : !! j=N is returned to indicate that x is out of range.
1612 : !!
1613 : !! SOURCE
1614 :
1615 6574347 : pure function bisect_rdp(AA, xx) result(loc)
1616 :
1617 : !Arguments ------------------------------------
1618 : !scalars
1619 : real(dp),intent(in) :: AA(:)
1620 : real(dp),intent(in) :: xx
1621 : integer :: loc
1622 :
1623 : !Local variables-------------------------------
1624 : integer :: nn,jl,jm,ju
1625 : logical :: ascnd
1626 : ! *********************************************************************
1627 :
1628 6574347 : nn=SIZE(AA); ascnd=(AA(nn)>=AA(1))
1629 : !
1630 : ! Initialize lower and upper limits
1631 6574347 : jl=0; ju=nn+1
1632 : do
1633 42652741 : if (ju-jl<=1) EXIT
1634 36078394 : jm=(ju+jl)/2 ! Compute a midpoint,
1635 42652741 : if (ascnd.EQV.(xx>=AA(jm))) then
1636 : jl=jm ! Replace lower limit
1637 : else
1638 28413311 : ju=jm ! Replace upper limit
1639 : end if
1640 : end do
1641 : !
1642 : ! Set the output, being careful with the endpoints
1643 6574347 : if (xx==AA(1)) then
1644 : loc=1
1645 6473034 : else if (xx==AA(nn)) then
1646 989 : loc=nn-1
1647 : else
1648 : loc=jl
1649 : end if
1650 :
1651 6574347 : end function bisect_rdp
1652 : !!***
1653 :
1654 : !----------------------------------------------------------------------
1655 :
1656 : !!****f* m_numeric_tools/bisect_int
1657 : !! NAME
1658 : !! bisect_int
1659 : !!
1660 : !! FUNCTION
1661 : !! Given an array AA(1:N), and a value x, returns the index j such that AA(j) <= x <= AA(j + 1).
1662 : !! AA must be monotonic, either increasing or decreasing. j=0 or
1663 : !! j=N is returned to indicate that x is out of range.
1664 : !!
1665 : !! INPUTS
1666 : !!
1667 : !! OUTPUT
1668 : !!
1669 : !! SOURCE
1670 0 : pure function bisect_int(AA,xx) result(loc)
1671 :
1672 : !Arguments ------------------------------------
1673 : !scalars
1674 : integer,intent(in) :: AA(:)
1675 : integer,intent(in) :: xx
1676 : integer :: loc
1677 :
1678 : !Local variables-------------------------------
1679 : integer :: nn,jl,jm,ju
1680 : logical :: ascnd
1681 : ! *********************************************************************
1682 :
1683 0 : nn=SIZE(AA) ; ascnd=(AA(nn)>=AA(1))
1684 :
1685 : ! Initialize lower and upper limits
1686 0 : jl=0 ; ju=nn+1
1687 : do
1688 0 : if (ju-jl<=1) EXIT
1689 0 : jm=(ju+jl)/2 ! Compute a midpoint
1690 0 : if (ascnd.EQV.(xx>=AA(jm))) then
1691 : jl=jm ! Replace lower limit
1692 : else
1693 0 : ju=jm ! Replace upper limit
1694 : end if
1695 : end do
1696 : !
1697 : ! Set the output, being careful with the endpoints
1698 0 : if (xx==AA(1)) then
1699 : loc=1
1700 0 : else if (xx==AA(nn)) then
1701 0 : loc=nn-1
1702 : else
1703 : loc=jl
1704 : end if
1705 :
1706 0 : end function bisect_int
1707 : !!***
1708 :
1709 : !----------------------------------------------------------------------
1710 :
1711 : !!****f* m_numeric_tools/imax_loc_int
1712 : !! NAME
1713 : !! imax_loc_int
1714 : !!
1715 : !! FUNCTION
1716 : !! Index of maxloc on an array returned as scalar instead of array-valued
1717 : !!
1718 : !! SOURCE
1719 :
1720 3 : pure function imax_loc_int(iarr,mask)
1721 :
1722 : !Arguments ------------------------------------
1723 : !scalars
1724 : integer :: imax_loc_int
1725 : !arrays
1726 : integer,intent(in) :: iarr(:)
1727 : logical,optional,intent(in) :: mask(:)
1728 :
1729 : !Local variables-------------------------------
1730 : integer :: imax(1)
1731 : ! *************************************************************************
1732 :
1733 3 : if (PRESENT(mask)) then
1734 0 : imax=MAXLOC(iarr,MASK=mask)
1735 : else
1736 17 : imax=MAXLOC(iarr)
1737 : end if
1738 3 : imax_loc_int=imax(1)
1739 :
1740 3 : end function imax_loc_int
1741 : !!***
1742 :
1743 : !----------------------------------------------------------------------
1744 :
1745 : !!****f* m_numeric_tools/imax_loc_rdp
1746 : !! NAME
1747 : !! imax_loc_rdp
1748 : !!
1749 : !! FUNCTION
1750 : !!
1751 : !! INPUTS
1752 : !!
1753 : !! OUTPUT
1754 : !!
1755 : !! SOURCE
1756 3616 : pure function imax_loc_rdp(arr,mask)
1757 :
1758 : !Arguments ------------------------------------
1759 : !scalars
1760 : integer :: imax_loc_rdp
1761 : !arrays
1762 : real(dp),intent(in) :: arr(:)
1763 : logical,optional,intent(in) :: mask(:)
1764 :
1765 : !Local variables-------------------------------
1766 : integer :: imax(1)
1767 : ! *************************************************************************
1768 :
1769 3616 : if (PRESENT(mask)) then
1770 0 : imax=MAXLOC(arr,MASK=mask)
1771 : else
1772 45926 : imax=MAXLOC(arr)
1773 : end if
1774 3616 : imax_loc_rdp=imax(1)
1775 :
1776 3616 : end function imax_loc_rdp
1777 : !!***
1778 :
1779 : !----------------------------------------------------------------------
1780 :
1781 : !!****f* m_numeric_tools/imin_loc_int
1782 : !! NAME
1783 : !! imin_loc_int
1784 : !!
1785 : !! FUNCTION
1786 : !! Index of minloc on an array returned as scalar instead of array-valued
1787 : !!
1788 : !! SOURCE
1789 :
1790 131598 : pure function imin_loc_int(arr, mask)
1791 :
1792 : !Arguments ------------------------------------
1793 : !scalars
1794 : integer :: imin_loc_int
1795 : !arrays
1796 : integer,intent(in) :: arr(:)
1797 : logical,optional,intent(in) :: mask(:)
1798 :
1799 : !Local variables-------------------------------
1800 : integer :: imin(1)
1801 : ! *************************************************************************
1802 :
1803 131598 : if (PRESENT(mask)) then
1804 699890 : imin=MINLOC(arr,MASK=mask)
1805 : else
1806 12 : imin=MINLOC(arr)
1807 : end if
1808 131598 : imin_loc_int=imin(1)
1809 :
1810 131598 : end function imin_loc_int
1811 : !!***
1812 :
1813 : !----------------------------------------------------------------------
1814 :
1815 : !!****f* m_numeric_tools/imin_loc_rdp
1816 : !! NAME
1817 : !! imin_loc_rdp
1818 : !!
1819 : !! FUNCTION
1820 : !!
1821 : !! INPUTS
1822 : !!
1823 : !! OUTPUT
1824 : !!
1825 : !! SOURCE
1826 :
1827 126890 : pure function imin_loc_rdp(arr,mask)
1828 :
1829 : !Arguments ------------------------------------
1830 : !scalars
1831 : integer :: imin_loc_rdp
1832 : !arrays
1833 : real(dp),intent(in) :: arr(:)
1834 : logical,optional,intent(in) :: mask(:)
1835 :
1836 : !Local variables-------------------------------
1837 : integer :: imin(1)
1838 : ! *************************************************************************
1839 :
1840 126890 : if (PRESENT(mask)) then
1841 0 : imin=MINLOC(arr,MASK=mask)
1842 : else
1843 1044474 : imin=MINLOC(arr)
1844 : end if
1845 :
1846 126890 : imin_loc_rdp=imin(1)
1847 :
1848 126890 : end function imin_loc_rdp
1849 : !!***
1850 :
1851 : !----------------------------------------------------------------------
1852 :
1853 : !!****f* m_numeric_tools/lfind
1854 : !! NAME
1855 : !! lfind
1856 : !!
1857 : !! FUNCTION
1858 : !! Find the index of the first occurrence of .True. in a logical array.
1859 : !! Return -1 if not found. If back is True, the search starts from the
1860 : !! last element of the array (default: False).
1861 : !!
1862 : !! INPUTS
1863 : !! mask(:)=Input logical mask
1864 : !!
1865 : !! SOURCE
1866 :
1867 0 : integer pure function lfind(mask, back)
1868 :
1869 : !Arguments ------------------------------------
1870 : !scalars
1871 : logical,intent(in) :: mask(:)
1872 : logical,optional,intent(in) :: back
1873 : !arrays
1874 :
1875 : !Local variables-------------------------------
1876 : !scalars
1877 : integer :: ii,nitems
1878 : logical :: do_back
1879 : !************************************************************************
1880 :
1881 0 : do_back = .False.; if (present(back)) do_back = back
1882 0 : lfind = -1; nitems = size(mask); if (nitems == 0) return
1883 :
1884 0 : if (do_back) then
1885 : ! Backward search
1886 0 : do ii=nitems,1,-1
1887 0 : if (mask(ii)) then
1888 0 : lfind = ii; return
1889 : end if
1890 : end do
1891 : else
1892 : ! Forward search.
1893 0 : do ii=1,nitems
1894 0 : if (mask(ii)) then
1895 0 : lfind = ii; return
1896 : end if
1897 : end do
1898 : end if
1899 :
1900 : end function lfind
1901 : !!***
1902 :
1903 : !----------------------------------------------------------------------
1904 :
1905 : !!****f* m_numeric_tools/list2blocks
1906 : !! NAME
1907 : !! list2blocks
1908 : !!
1909 : !! FUNCTION
1910 : !! Given a list of integers, find the number of contiguous groups of values.
1911 : !! and returns the set of indices that can be used to loop over these groups
1912 : !! Example list = [1,2,3,5,6] --> blocks = [[1,3], [4,5]]
1913 : !!
1914 : !! INPUTS
1915 : !! list(:)=List of integers
1916 : !!
1917 : !! OUTPUTS
1918 : !! nblocks=Number of blocks
1919 : !! blocks(2,nblocks)=
1920 : !! allocatable array in input
1921 : !! in output:
1922 : !! blocks(1,i) gives the start of the i-th block
1923 : !! blocks(2,i) gives the end of the i-th block
1924 : !!
1925 : !! SOURCE
1926 :
1927 30 : subroutine list2blocks(list,nblocks,blocks)
1928 :
1929 : !Arguments ------------------------------------
1930 : !scalars
1931 : integer,intent(out) :: nblocks
1932 : integer,intent(in) :: list(:)
1933 : !arrays
1934 : integer,intent(out),allocatable :: blocks(:,:)
1935 :
1936 : !Local variables-------------------------------
1937 : !scalars
1938 : integer :: ii,nitems
1939 : !arrays
1940 30 : integer :: work(2,size(list))
1941 : !************************************************************************
1942 :
1943 30 : nitems = size(list)
1944 :
1945 : ! Handle nitems == 1 case
1946 30 : if (nitems == 1) then
1947 0 : ABI_MALLOC(blocks, (2,1))
1948 0 : blocks = 1
1949 : return
1950 : end if
1951 :
1952 30 : nblocks = 1; work(1,1) = 1
1953 :
1954 240 : do ii=2,nitems
1955 240 : if (list(ii) /= (list(ii-1) + 1)) then
1956 0 : work(2,nblocks) = ii - 1
1957 0 : nblocks = nblocks + 1
1958 0 : work(1,nblocks) = ii
1959 : end if
1960 : end do
1961 :
1962 30 : work(2,nblocks) = nitems
1963 :
1964 90 : ABI_MALLOC(blocks, (2,nblocks))
1965 150 : blocks = work(:,1:nblocks)
1966 :
1967 : end subroutine list2blocks
1968 : !!***
1969 :
1970 : !----------------------------------------------------------------------
1971 :
1972 : !!****f* m_numeric_tools/mask2blocks
1973 : !! NAME
1974 : !! mask2blocks
1975 : !!
1976 : !! FUNCTION
1977 : !! Give a logical mask, find the number of contiguous groups of .TRUE. values.
1978 : !! and return the set of indices that can be used to loop over these groups
1979 : !!
1980 : !! INPUTS
1981 : !! mask(:)=Input logical mask
1982 : !!
1983 : !! OUTPUTS
1984 : !! nblocks=Number of blocks
1985 : !!
1986 : !! SIDE EFFECTS
1987 : !! blocks(:,:)= Null pointer in input. blocks(2,nblocks) in output where
1988 : !! blocks(1,i) gives the start of the i-th block
1989 : !! blocks(2,i) gives the end of the i-th block
1990 : !!
1991 : !! SOURCE
1992 :
1993 0 : subroutine mask2blocks(mask,nblocks,blocks)
1994 :
1995 : !Arguments ------------------------------------
1996 : !scalars
1997 : integer,intent(out) :: nblocks
1998 : logical,intent(in) :: mask(:)
1999 : !arrays
2000 : integer,allocatable :: blocks(:,:)
2001 :
2002 : !Local variables-------------------------------
2003 : !scalars
2004 : integer :: ii,nitems,start
2005 : logical :: inblock
2006 : !arrays
2007 0 : integer :: work(2,SIZE(mask))
2008 : !************************************************************************
2009 :
2010 : ! Find first element.
2011 0 : nitems = size(mask); start = 0
2012 0 : do ii=1,nitems
2013 0 : if (mask(ii)) then
2014 : start = ii
2015 : exit
2016 : end if
2017 : end do
2018 :
2019 : ! Handle no true element or just one.
2020 0 : if (start == 0) then
2021 0 : nblocks = 0
2022 0 : ABI_MALLOC(blocks, (0,0))
2023 : return
2024 : end if
2025 0 : if (start /= 0 .and. nitems == 1) then
2026 0 : nblocks = 1
2027 0 : ABI_MALLOC(blocks, (2,1))
2028 0 : blocks(:,1) = [1,1]
2029 : end if
2030 :
2031 0 : nblocks = 1; work(1,1) = start; inblock = .True.
2032 :
2033 0 : do ii=start+1,nitems
2034 0 : if (.not.mask(ii)) then
2035 0 : if (inblock) then
2036 0 : inblock = .False.
2037 0 : work(2,nblocks) = ii - 1
2038 : end if
2039 : else
2040 0 : if (.not. inblock) then
2041 0 : inblock = .True.
2042 0 : nblocks = nblocks + 1
2043 0 : work(1,nblocks) = ii
2044 : end if
2045 : end if
2046 : end do
2047 :
2048 0 : if (mask(nitems) .and. inblock) work(2,nblocks) = nitems
2049 :
2050 0 : ABI_MALLOC(blocks, (2,nblocks))
2051 0 : blocks = work(:,1:nblocks)
2052 :
2053 : end subroutine mask2blocks
2054 : !!***
2055 :
2056 : !----------------------------------------------------------------------
2057 :
2058 : !!****f* m_numeric_tools/linfit_rdp
2059 : !! NAME
2060 : !! linfit_rdp
2061 : !!
2062 : !! FUNCTION
2063 : !! Perform a linear fit, y=ax+b, of data
2064 : !!
2065 : !! INPUTS
2066 : !! xx(nn)=xx coordinates
2067 : !! yy(nn)=yy coordinates
2068 : !!
2069 : !! OUTPUT
2070 : !! aa=coefficient of linear term of fit
2071 : !! bb=coefficient of constant term of fit
2072 : !! res=root mean square of differences between data and fit
2073 : !!
2074 : !! SOURCE
2075 :
2076 1804 : function linfit_rdp(nn,xx,yy,aa,bb) result(res)
2077 :
2078 : !Arguments ------------------------------------
2079 : !scalars
2080 : integer,intent(in) :: nn
2081 : real(dp) :: res
2082 : real(dp),intent(out) :: aa,bb
2083 : !arrays
2084 : real(dp),intent(in) :: xx(nn),yy(nn)
2085 :
2086 : !Local variables-------------------------------
2087 : !scalars
2088 : integer :: ii
2089 : real(dp) :: msrt,sx2,sx,sxy,sy,tx,ty
2090 : ! *************************************************************************
2091 :
2092 1804 : sx=zero ; sy=zero ; sxy=zero ; sx2=zero
2093 23048 : do ii=1,nn
2094 21244 : tx=xx(ii)
2095 21244 : ty=yy(ii)
2096 21244 : sx=sx+tx
2097 21244 : sy=sy+ty
2098 21244 : sxy=sxy+tx*ty
2099 23048 : sx2=sx2+tx*tx
2100 : end do
2101 :
2102 1804 : aa=(nn*sxy-sx*sy)/(nn*sx2-sx*sx)
2103 1804 : bb=sy/nn-sx*aa/nn
2104 :
2105 1804 : msrt=zero
2106 23048 : do ii=1,nn
2107 21244 : tx=xx(ii)
2108 21244 : ty=yy(ii)
2109 23048 : msrt=msrt+(ty-aa*tx-bb)**2
2110 : end do
2111 1804 : msrt=SQRT(msrt/nn) ; res=msrt
2112 :
2113 1804 : end function linfit_rdp
2114 : !!***
2115 :
2116 : !----------------------------------------------------------------------
2117 :
2118 : !!****f* m_numeric_tools/linfit_spc
2119 : !! NAME
2120 : !! linfit_spc
2121 : !!
2122 : !! FUNCTION
2123 : !! Perform a linear fit, y=ax+b, of data
2124 : !!
2125 : !! INPUTS
2126 : !!
2127 : !! OUTPUT
2128 : !!
2129 : !! SOURCE
2130 :
2131 0 : function linfit_spc(nn,xx,zz,aa,bb) result(res)
2132 :
2133 : !Arguments ------------------------------------
2134 : !scalars
2135 : integer,intent(in) :: nn
2136 : real(dp) :: res
2137 : real(dp),intent(in) :: xx(nn)
2138 : complex(sp),intent(in) :: zz(nn)
2139 : complex(sp),intent(out) :: aa,bb
2140 :
2141 : !Local variables-------------------------------
2142 : !scalars
2143 : integer :: ii
2144 : real(dp) :: sx,sx2,msrt
2145 : complex(dp) :: sz,sxz
2146 : ! *************************************************************************
2147 :
2148 0 : sx=zero ; sx2=zero ; msrt=zero
2149 0 : sz=czero ; sxz=czero
2150 0 : do ii=1,nn
2151 0 : sx=sx+xx(ii)
2152 0 : sz=sz+zz(ii)
2153 0 : sxz=sxz+xx(ii)*zz(ii)
2154 0 : sx2=sx2+xx(ii)*xx(ii)
2155 : end do
2156 :
2157 0 : aa=CMPLX((nn*sxz-sx*sz)/(nn*sx2-sx*sx), kind=sp)
2158 0 : bb=CMPLX(sz/nn-sx*aa/nn, kind=sp)
2159 :
2160 0 : do ii=1,nn
2161 0 : msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2162 : end do
2163 0 : msrt=SQRT(msrt) ; res=msrt
2164 :
2165 0 : end function linfit_spc
2166 : !!***
2167 :
2168 : !----------------------------------------------------------------------
2169 :
2170 : !!****f* m_numeric_tools/linfit_dpc
2171 : !! NAME
2172 : !! linfit_dpc
2173 : !!
2174 : !! FUNCTION
2175 : !! Perform a linear fit, y=ax+b, of data
2176 : !!
2177 : !! INPUTS
2178 : !!
2179 : !! OUTPUT
2180 : !!
2181 : !! SOURCE
2182 :
2183 0 : function linfit_dpc(nn,xx,zz,aa,bb) result(res)
2184 :
2185 : !Arguments ------------------------------------
2186 : !scalars
2187 : integer,intent(in) :: nn
2188 : real(dp) :: res
2189 : real(dp),intent(in) :: xx(nn)
2190 : complex(dp),intent(in) :: zz(nn)
2191 : complex(dp),intent(out) :: aa,bb
2192 : !arrays
2193 :
2194 : !Local variables-------------------------------
2195 : !scalars
2196 : integer :: ii
2197 : real(dp) :: sx,sx2,msrt
2198 : complex(dp) :: sz,sxz
2199 : ! *************************************************************************
2200 :
2201 0 : sx=zero ; sx2=zero ; msrt=zero
2202 0 : sz=czero ; sxz=czero
2203 0 : do ii=1,nn
2204 0 : sx=sx+xx(ii)
2205 0 : sz=sz+zz(ii)
2206 0 : sxz=sxz+xx(ii)*zz(ii)
2207 0 : sx2=sx2+xx(ii)*xx(ii)
2208 : end do
2209 :
2210 0 : aa=(nn*sxz-sx*sz)/(nn*sx2-sx*sx)
2211 0 : bb=sz/nn-sx*aa/nn
2212 :
2213 0 : do ii=1,nn
2214 0 : msrt=msrt+ABS(zz(ii)-aa*xx(ii)-bb)**2
2215 : end do
2216 0 : msrt=SQRT(msrt) ; res=msrt
2217 :
2218 0 : end function linfit_dpc
2219 : !!***
2220 :
2221 : !----------------------------------------------------------------------
2222 :
2223 : !!****f* m_numeric_tools/llsfit_svd
2224 : !! NAME
2225 : !! llsfit_svd
2226 : !!
2227 : !! FUNCTION
2228 : !! Given a set of N data points (x,y) with individual standard deviations sigma_i,
2229 : !! use chi-square minimization to determine the M coefficients, par, of a function that
2230 : !! depends linearly on nfuncs functions, i.e f(x) = \sum_i^{nfuncs} par_i * func_i(x).
2231 : !! Solve the fitting equations using singular value decomposition of the design matrix as in Eq 14.3.17
2232 : !! of Numerical Recipes. The program returns values for the M fit parameters par, and chi-square.
2233 : !! The user supplies a subroutine funcs(x,nfuncs) that returns the M basis functions evaluated at xx.
2234 : !!
2235 : !! INPUTS
2236 : !!
2237 : !! OUTPUT
2238 : !!
2239 : !! SOURCE
2240 :
2241 0 : subroutine llsfit_svd(xx,yy,sigma,nfuncs,funcs,chisq,par,var,cov,info)
2242 :
2243 : !Arguments ------------------------------------
2244 : !scalars
2245 : integer,intent(in) :: nfuncs
2246 : integer,intent(out) :: info
2247 : real(dp),intent(out) :: chisq
2248 : !arrays
2249 : real(dp),intent(in) :: xx(:),yy(:),sigma(:)
2250 : real(dp),intent(out) :: par(:),var(:),cov(:,:)
2251 :
2252 : interface
2253 : function funcs(xx,nf)
2254 : use defs_basis
2255 : real(dp),intent(in) :: xx
2256 : integer,intent(in) :: nf
2257 : real(dp) :: funcs(nf)
2258 : end function funcs
2259 : end interface
2260 :
2261 : !Local variables-------------------------------
2262 : integer,parameter :: PAD_=50
2263 : integer :: ii,npts,lwork
2264 : real(dp),parameter :: TOL_=1.0e-5_dp
2265 : !arrays
2266 0 : real(dp),dimension(SIZE(xx)) :: bb,sigm1
2267 0 : real(dp),dimension(SIZE(xx),nfuncs) :: dmat,dmat_save
2268 0 : real(dp) :: tmp(nfuncs)
2269 0 : real(dp),allocatable :: work(:),Vt(:,:),U(:,:),S(:)
2270 : ! *************************************************************************
2271 :
2272 0 : npts = assert_eq(SIZE(xx),SIZE(yy),SIZE(sigma),'Wrong size in xx,yy,sigma', __FILE__, __LINE__)
2273 0 : call assert((npts>=nfuncs),'No. of functions must greater than no. of points', __FILE__, __LINE__)
2274 0 : ii = assert_eq(nfuncs,SIZE(cov,1),SIZE(cov,2),SIZE(var),'Wrong size in covariance', __FILE__, __LINE__)
2275 :
2276 : !
2277 : ! === Calculate design matrix and b vector ===
2278 : ! * dmat_ij=f_j(x_i)/sigma_i, b_i=y_i/sigma_i
2279 0 : sigm1(:)=one/sigma(:) ; bb(:)=yy(:)*sigm1(:)
2280 0 : do ii=1,npts
2281 0 : dmat_save(ii,:)=funcs(xx(ii),nfuncs)
2282 : end do
2283 0 : dmat=dmat_save*SPREAD(sigm1,DIM=2,ncopies=nfuncs)
2284 0 : dmat_save(:,:)=dmat(:,:)
2285 : !
2286 : ! === Singular value decomposition ===
2287 0 : lwork=MAX(3*MIN(npts,nfuncs)+MAX(npts,nfuncs),5*MIN(npts,nfuncs)-4)+PAD_
2288 0 : ABI_MALLOC(work,(lwork))
2289 0 : ABI_MALLOC(U,(npts,npts))
2290 0 : ABI_MALLOC(S,(nfuncs))
2291 0 : ABI_MALLOC(Vt,(nfuncs,nfuncs))
2292 :
2293 0 : call DGESVD('A','A',npts,nfuncs,dmat,npts,S,U,npts,Vt,nfuncs,work,lwork,info)
2294 0 : ABI_FREE(work)
2295 : GOTO 10
2296 : !
2297 : ! === Set to zero small singular values according to TOL_ and find coefficients ===
2298 : WHERE (S>TOL_*MAXVAL(S))
2299 : tmp=MATMUL(bb,U)/S
2300 : ELSEWHERE
2301 : S =zero
2302 : tmp=zero
2303 : END WHERE
2304 : par(:)=MATMUL(tmp,Vt)
2305 : !
2306 : ! === Evaluate chi-square ===
2307 : chisq=l2norm(MATMUL(dmat_save,par)-bb)**2
2308 : !
2309 : ! === Calculate covariance and variance ===
2310 : ! C_jk = V_ji V_ki / S_i^2
2311 : WHERE (S/=zero) S=one/(S*S)
2312 :
2313 : ! check this but should be correct
2314 : cov(:,:)=Vt*SPREAD(S,DIM=2,ncopies=nfuncs)
2315 : cov(:,:)=MATMUL(TRANSPOSE(Vt),cov)
2316 : var(:)=SQRT(get_diag(cov))
2317 :
2318 : 10 continue
2319 0 : ABI_FREE(U)
2320 0 : ABI_FREE(S)
2321 0 : ABI_FREE(Vt)
2322 :
2323 0 : end subroutine llsfit_svd
2324 : !!***
2325 :
2326 : !----------------------------------------------------------------------
2327 :
2328 : !!****f* m_numeric_tools/polyn_interp
2329 : !! NAME
2330 : !! polyn_interp
2331 : !!
2332 : !! FUNCTION
2333 : !! Given arrays xa and ya of length N, and given a value x, return a value y, and an error estimate dy.
2334 : !! 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).
2335 : !!
2336 : !! INPUTS
2337 : !! xa(:)=abscissas in ascending order
2338 : !! ya(:)=ordinates
2339 : !! x=the point where the set of data has to be interpolated
2340 : !!
2341 : !! OUTPUT
2342 : !! y=the interpolated value
2343 : !! dy=error estimate
2344 : !!
2345 : !! NOTES
2346 : !! Based on the polint routine reported in Numerical Recipes
2347 : !!
2348 : !! SOURCE
2349 :
2350 120544 : subroutine polyn_interp(xa,ya,x,y,dy)
2351 :
2352 : !Arguments ------------------------------------
2353 : !scalars
2354 : real(dp),intent(in) :: xa(:),ya(:)
2355 : real(dp),intent(in) :: x
2356 : real(dp),intent(out) :: y,dy
2357 : !Local variables-------------------------------
2358 : !scalars
2359 : integer :: m,n,ns
2360 : !arrays
2361 241088 : real(dp),dimension(SIZE(xa)) :: c,d,den,ho
2362 : ! *************************************************************************
2363 :
2364 120544 : n = assert_eq(SIZE(xa),SIZE(ya),'Different size in xa and ya',__FILE__,__LINE__)
2365 :
2366 : ! === Initialize the tables of c and d ===
2367 1928704 : c(:)=ya(:) ; d(:)=ya(:) ; ho(:)=xa(:)-x
2368 : ! === Find closest table entry and initial approximation to y ===
2369 723264 : ns=imin_loc(ABS(x-xa)) ; y=ya(ns)
2370 120544 : ns=ns-1
2371 : !
2372 : ! === For each column of the tableau loop over current c and d and up-date them ===
2373 602720 : do m=1,n-1
2374 1687616 : den(1:n-m)=ho(1:n-m)-ho(1+m:n)
2375 1687616 : if (ANY(den(1:n-m)==zero)) then
2376 0 : ABI_ERROR('Two input xa are identical')
2377 : end if
2378 :
2379 1687616 : den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
2380 1687616 : d(1:n-m)=ho(1+m:n)*den(1:n-m) ! Update c and d
2381 1687616 : c(1:n-m)=ho(1:n-m)*den(1:n-m)
2382 :
2383 482176 : if (2*ns<n-m) then ! Now decide which correction, c or d, we want to add to the
2384 0 : dy=c(ns+1) ! accumulating value of y, The last dy added is the error indication.
2385 : else
2386 482176 : dy=d(ns)
2387 482176 : ns=ns-1
2388 : end if
2389 :
2390 602720 : y=y+dy
2391 : end do
2392 :
2393 120544 : end subroutine polyn_interp
2394 : !!***
2395 :
2396 : !----------------------------------------------------------------------
2397 :
2398 : !!****f* m_numeric_tools/polcoe
2399 : !! NAME
2400 : !! polcoe
2401 : !!
2402 : !! FUNCTION
2403 : !! Given arrays x(1:n) and y(1:n) containing a tabulated function yi = f (xi ), this routine
2404 : !! returns an array with the coefficients cof(1:n) of a polynomial interpolation.
2405 : !!
2406 : !! INPUTS
2407 : !! x(n)=abscissas in ascending order
2408 : !! y(n)=ordinates
2409 : !! n=number of points given to start the interpolation
2410 : !!
2411 : !! OUTPUT
2412 : !! cof(n)= coefficients array
2413 : !!
2414 : !! NOTES
2415 : !! Based on the polcoe routine reported in Numerical Recipies
2416 : !!
2417 : !! SOURCE
2418 :
2419 3248 : subroutine polcoe(x,y,n,cof)
2420 :
2421 : !Arguments ------------------------------------
2422 : !scalars
2423 : integer, intent(in) :: n
2424 : !arrays
2425 : real(dp), intent(in) :: x(n),y(n)
2426 : real(dp), intent(out) :: cof(n)
2427 : !Local variables ------------------------------
2428 : !scalars
2429 : integer, parameter :: NMAX=15
2430 : integer :: i,j,k
2431 : real(dp) :: b,ff,phi
2432 : !arrays
2433 : real(dp) :: s(NMAX)
2434 : ! *************************************************************************
2435 :
2436 9744 : do 11 i=1,n
2437 6496 : s(i)=0.
2438 6496 : cof(i)=0.
2439 3248 : 11 continue
2440 3248 : s(n)=-x(1)
2441 6496 : do 13 i=2,n
2442 6496 : do 12 j=n+1-i,n-1
2443 3248 : s(j)=s(j)-x(i)*s(j+1)
2444 3248 : 12 continue
2445 3248 : s(n)=s(n)-x(i)
2446 3248 : 13 continue
2447 9744 : do 16 j=1,n
2448 6496 : phi=n
2449 12992 : do 14 k=n-1,1,-1
2450 6496 : phi=k*s(k+1)+x(j)*phi
2451 6496 : 14 continue
2452 6496 : ff=y(j)/phi
2453 6496 : b=1.
2454 19488 : do 15 k=n,1,-1
2455 12992 : cof(k)=cof(k)+b*ff
2456 12992 : b=s(k)+x(j)*b
2457 6496 : 15 continue
2458 3248 : 16 continue
2459 :
2460 3248 : end subroutine polcoe
2461 : !!***
2462 :
2463 : !----------------------------------------------------------------------
2464 :
2465 : !!****f* m_numeric_tools/trapezoidal_
2466 : !! NAME
2467 : !! trapezoidal_ (PRIVATE)
2468 : !!
2469 : !! FUNCTION
2470 : !! Compute the n-th stage of refinement of an extended trapezoidal rule
2471 : !! adding 2^(n-2) additional interior point in the finite range of integration
2472 : !!
2473 : !! INPUTS
2474 : !! func(external)=the name of the function to be integrated
2475 : !! xmin,xmax=the limits of integration
2476 : !! nn=integer defining the refinement of the mesh, each call adds 2^(n-2) additional interior points
2477 : !!
2478 : !! OUTPUT
2479 : !! See SIDE EFFECTS
2480 : !!
2481 : !! SIDE EFFECTS
2482 : !! quad=the integral at the n-th stage.
2483 : !!
2484 : !! NOTES
2485 : !! When called with nn=1, the routine returns the crudest estimate of the integral
2486 : !! Subsequent calls with nn=2,3,... (in that sequential order) will improve the accuracy
2487 : !! by adding 2^(n-2) additional interior points. Note that quad should not be modified between sequential calls.
2488 : !! Subroutine is defined as recursive to allow multi-dimensional integrations
2489 : !!
2490 : !! SOURCE
2491 :
2492 0 : recursive subroutine trapezoidal_(func,nn,xmin,xmax,quad)
2493 :
2494 : !Arguments ------------------------------------
2495 : !scalars
2496 : integer,intent(in) :: nn
2497 : !real(dp),external :: func
2498 : real(dp),intent(in) :: xmin,xmax
2499 : real(dp),intent(inout) :: quad
2500 :
2501 : interface
2502 : function func(x)
2503 : use defs_basis
2504 : real(dp),intent(in) :: x
2505 : real(dp) :: func
2506 : end function func
2507 : end interface
2508 :
2509 : !interface
2510 : ! function func(x)
2511 : ! use defs_basis
2512 : ! real(dp),intent(in) :: x(:)
2513 : ! real(dp) :: func(SIZE(x))
2514 : ! end function func
2515 : !end interface
2516 :
2517 : !Local variables-------------------------------
2518 : !scalars
2519 : integer :: npt,ix
2520 : real(dp) :: space,new,yy
2521 : character(len=500) :: msg
2522 : !arrays
2523 : !real(dp),allocatable :: xx(:)
2524 : !************************************************************************
2525 :
2526 0 : select case (nn)
2527 :
2528 : case (1)
2529 : ! === Initial crude estimate (xmax-xmin)(f1+f2)/2 ===
2530 : !quad=half*(xmax-xmin)*SUM(func((/xmin,xmax/)))
2531 0 : quad=half*(xmax-xmin)*(func(xmin)+func(xmax))
2532 :
2533 : case (2:)
2534 : ! === Add npt interior points of spacing space ===
2535 0 : npt=2**(nn-2) ; space=(xmax-xmin)/npt
2536 : ! === The new sum is combined with the old integral to give a refined integral ===
2537 : !new=SUM(func(arth(xmin+half*space,space,npt))) !PARALLEL version
2538 : !allocate(xx(npt))
2539 : !xx(:)=arth(xmin+half*space,space,npt)
2540 : !xx(1)=xmin+half*space
2541 : !do ii=2,nn
2542 : ! xx(ii)=xx(ii-1)+space
2543 : !end do
2544 0 : new=zero
2545 0 : yy=xmin+half*space
2546 0 : do ix=1,npt
2547 : !new=new+func(xx(ix))
2548 0 : new=new+func(yy)
2549 0 : yy=yy+space
2550 : end do
2551 : !deallocate(xx)
2552 0 : quad=half*(quad+space*new)
2553 : !write(std_out,*) 'trapezoidal',quad
2554 :
2555 : case (:0)
2556 0 : write(msg,'(a,i3)')'Wrong value for nn ',nn
2557 0 : ABI_BUG(msg)
2558 : end select
2559 :
2560 0 : end subroutine trapezoidal_
2561 : !!***
2562 :
2563 : !----------------------------------------------------------------------
2564 :
2565 : !!****f* m_numeric_tools/midpoint_
2566 : !! NAME
2567 : !! midpoint_ (PRIVATE)
2568 : !!
2569 : !! FUNCTION
2570 : !! This routine computes the n-th stage of refinement of an extended midpoint rule.
2571 : !!
2572 : !! INPUTS
2573 : !! func(external)=the name of the function to be integrated
2574 : !! xmin,xmax=the limits of integration
2575 : !! nn=integer defining the refinement of the mesh, each call adds (2/3)*3n-1 additional
2576 : !! interior points between xmin ans xmax
2577 : !!
2578 : !! OUTPUT
2579 : !! See SIDE EFFECTS
2580 : !!
2581 : !! SIDE EFFECTS
2582 : !! quad=the integral at the n-th stage.
2583 : !!
2584 : !! NOTES
2585 : !! When called with nn=1, the routine returns as quad the crudest estimate of the integral
2586 : !! Subsequent calls with nn=2,3,... (in that sequential order) will improve the accuracy of quad by adding
2587 : !! (2/3)*3n-1 additional interior points. quad should not be modified between sequential calls.
2588 : !! Subroutine is defined as recursive to allow multi-dimensional integrations
2589 : !!
2590 : !! SOURCE
2591 :
2592 602608 : recursive subroutine midpoint_(func,nn,xmin,xmax,quad)
2593 :
2594 : !Arguments ------------------------------------
2595 : !scalars
2596 : integer,intent(in) :: nn
2597 : !real(dp),external :: func
2598 : real(dp),intent(in) :: xmin,xmax
2599 : real(dp),intent(inout) :: quad
2600 :
2601 : interface
2602 : function func(x)
2603 : use defs_basis
2604 : real(dp),intent(in) :: x
2605 : real(dp) :: func
2606 : end function func
2607 : end interface
2608 :
2609 : !interface
2610 : ! function func(x)
2611 : ! use defs_basis
2612 : ! real(dp),intent(in) :: x(:)
2613 : ! real(dp) :: func(SIZE(x))
2614 : ! end function func
2615 : !end interface
2616 :
2617 : !Local variables-------------------------------
2618 : !scalars
2619 : integer :: npt,ix
2620 : real(dp) :: space
2621 : character(len=500) :: msg
2622 : !arrays
2623 602608 : real(dp),allocatable :: xx(:)
2624 : !************************************************************************
2625 :
2626 723124 : select case (nn)
2627 :
2628 : case (1)
2629 : ! === Initial crude estimate done at the middle of the interval
2630 : !quad=(xmax-xmin)*SUM(func((/half*(xmin+xmax)/))) !PARALLEL version
2631 120516 : quad=(xmax-xmin)*func(half*(xmin+xmax))
2632 :
2633 : case (2:)
2634 : ! === Add npt interior points, they alternate in spacing between space and 2*space ===
2635 1446276 : ABI_MALLOC(xx,(2*3**(nn-2)))
2636 482092 : npt=3**(nn-2) ; space=(xmax-xmin)/(three*npt)
2637 5305000 : xx(1:2*npt-1:2)=arth(xmin+half*space,three*space,npt)
2638 5305000 : xx(2:2*npt:2)=xx(1:2*npt-1:2)+two*space
2639 : ! === The new sum is combined with the old integral to give a refined integral ===
2640 : !quad=quad/three+space*SUM(func(xx)) !PARALLEL version
2641 482092 : quad=quad/three
2642 10127908 : do ix=1,SIZE(xx)
2643 10127908 : quad=quad+space*func(xx(ix))
2644 : end do
2645 482092 : ABI_FREE(xx)
2646 :
2647 : case (:0)
2648 0 : write(msg,'(a,i3)')' wrong value for nn ',nn
2649 602608 : ABI_BUG('Wrong value for nn')
2650 : end select
2651 :
2652 602608 : end subroutine midpoint_
2653 : !!***
2654 :
2655 : !----------------------------------------------------------------------
2656 :
2657 : !!****f* m_numeric_tools/quadrature
2658 : !! NAME
2659 : !! quadrature
2660 : !!
2661 : !! FUNCTION
2662 : !! Driver routine to perform quadratures in finite domains using different techniques.
2663 : !! The routine improves the resolution of the grid until a given accuracy is reached
2664 : !!
2665 : !! INPUTS
2666 : !! func(external)=the function to be integrated
2667 : !! xmin,xmax=the limits of integration
2668 : !! npts=Initial number of points, only for Gauss-Legendre. At each step this number is doubled
2669 : !! accuracy=fractional accuracy required
2670 : !! ntrial=Max number of attempts
2671 : !! qopt=integer flag defining the algorithm for the quadrature:
2672 : !! 1 for Trapezoidal rule, closed, O(1/N^2)
2673 : !! 2 for Simpson based on trapezoidal,closed, O(1/N^4)
2674 : !! 3 for Midpoint rule, open, O(1/N^2)
2675 : !! 4 for midpoint rule with cancellation of leading error, open, O(1/N^4)
2676 : !! 5 for Romberg integration (closed form) and extrapolation for h-->0 (order 10 is hard-coded)
2677 : !! 6 for Romberg integration with midpoint rule and extrapolation for h-->0 (order 10 is hard-coded)
2678 : !! 7 for Gauss-Legendre
2679 : !!
2680 : !! OUTPUT
2681 : !! quad=the integral
2682 : !! ierr=0 if quadrature converged.
2683 : !!
2684 : !! SOURCE
2685 :
2686 120516 : recursive subroutine quadrature(func,xmin,xmax,qopt,quad,ierr,ntrial,accuracy,npts)
2687 :
2688 : !Arguments ------------------------------------
2689 : !scalars
2690 : integer,intent(in) :: qopt
2691 : integer,intent(out) :: ierr
2692 : integer,optional,intent(in) :: ntrial,npts
2693 : real(dp),intent(in) :: xmin,xmax
2694 : real(dp),optional,intent(in) :: accuracy
2695 : real(dp),intent(out) :: quad
2696 :
2697 : interface
2698 : function func(x)
2699 : use defs_basis
2700 : real(dp),intent(in) :: x
2701 : real(dp) :: func
2702 : end function func
2703 : end interface
2704 :
2705 : !interface
2706 : ! function func(x)
2707 : ! use defs_basis
2708 : ! real(dp),intent(in) :: x(:)
2709 : ! real(dp) :: func(SIZE(x))
2710 : ! end function func
2711 : !end interface
2712 :
2713 : !Local variables-------------------------------
2714 : !scalars
2715 : integer :: K,KM,NT,NX,NX0,it,ix
2716 : real(dp) :: EPS,old_st,st,old_quad,dqromb
2717 : real(dp) :: TOL
2718 : character(len=500) :: msg
2719 : !arrays
2720 120516 : real(dp),allocatable :: h(:),s(:), wx(:),xx(:)
2721 : ! *************************************************************************
2722 :
2723 120516 : ierr = 0
2724 120516 : TOL =tol12
2725 120516 : EPS =tol6 ; if (PRESENT(accuracy)) EPS=accuracy
2726 120516 : NT =20 ; if (PRESENT(ntrial )) NT=ntrial
2727 120516 : quad =zero
2728 :
2729 120516 : select case (qopt)
2730 :
2731 : case (1)
2732 : ! === Trapezoidal, closed form, O(1/N^2)
2733 0 : do it=1,NT
2734 0 : call trapezoidal_(func,it,xmin,xmax,quad)
2735 0 : if (it>5) then ! Avoid spurious early convergence
2736 0 : if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2737 : end if
2738 0 : old_quad=quad
2739 : end do
2740 :
2741 : case (2)
2742 : ! === Extended Simpson rule based on trapezoidal O(1/N^4) ===
2743 0 : do it=1,NT
2744 0 : call trapezoidal_(func,it,xmin,xmax,st)
2745 0 : if (it==1) then
2746 0 : quad=st
2747 : else
2748 0 : quad=(four*st-old_st)/three
2749 : end if
2750 0 : if (it>5) then ! Avoid spurious early convergence
2751 0 : if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2752 : end if
2753 0 : old_quad=quad
2754 0 : old_st=st
2755 : end do
2756 :
2757 : case (3)
2758 : ! === Midpoint rule, open form, O(1/N^2) ===
2759 0 : do it=1,NT
2760 0 : call midpoint_(func,it,xmin,xmax,quad)
2761 0 : if (it>4) then ! Avoid spurious early convergence
2762 0 : if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2763 : end if
2764 0 : old_quad=quad
2765 : end do
2766 :
2767 : case (4)
2768 : ! === Midpoint rule with cancellation of leading 1/N^2 term, open form, O(1/N^4) ===
2769 0 : do it=1,NT
2770 0 : call midpoint_(func,it,xmin,xmax,st)
2771 0 : if (it==1) then
2772 0 : quad=st
2773 : else
2774 0 : quad=(nine*st-old_st)/eight
2775 : end if
2776 0 : if (it>4) then ! Avoid spurious early convergence
2777 0 : if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2778 : end if
2779 0 : old_quad=quad
2780 0 : old_st=st
2781 : end do
2782 :
2783 : case (5)
2784 : ! === Romberg Integration, closed form ===
2785 0 : K=5 ; KM=K-1 ! Order 10
2786 0 : ABI_MALLOC(h,(NT+1))
2787 0 : ABI_MALLOC(s,(NT+1))
2788 0 : h=zero
2789 0 : s=zero
2790 0 : h(1)=one
2791 0 : do it=1,NT
2792 0 : call trapezoidal_(func,it,xmin,xmax,s(it))
2793 : !write(std_out,*) ' romberg-trap at ',ncall,it,s(it)
2794 0 : if (it>=K) then
2795 0 : call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
2796 0 : if (ABS(dqromb)<EPS*ABS(quad)) then
2797 0 : ABI_FREE(h)
2798 0 : ABI_FREE(s)
2799 0 : RETURN
2800 : end if
2801 : end if
2802 0 : s(it+1)=s(it)
2803 0 : h(it+1)=quarter*h(it) ! Quarter makes the extrapolation a polynomial in h^2,
2804 : end do ! This is required to use the Euler-Maclaurin formula
2805 0 : ABI_FREE(h)
2806 0 : ABI_FREE(s)
2807 :
2808 : case (6)
2809 : ! === Romberg Integration, closed form ===
2810 120516 : K=5 ; KM=K-1 ! Order 10
2811 361548 : ABI_MALLOC(h,(NT+1))
2812 241032 : ABI_MALLOC(s,(NT+1))
2813 3856512 : h=zero
2814 3856512 : s=zero
2815 120516 : h(1)=one
2816 602608 : do it=1,NT
2817 602608 : call midpoint_(func,it,xmin,xmax,s(it))
2818 602608 : if (it>=K) then
2819 120544 : call polyn_interp(h(it-KM:it),s(it-KM:it),zero,quad,dqromb)
2820 : !write(std_out,*) quad,dqromb
2821 120544 : if (ABS(dqromb)<EPS*ABS(quad)) then
2822 120516 : ABI_FREE(h)
2823 120516 : ABI_FREE(s)
2824 120516 : RETURN
2825 : end if
2826 : end if
2827 482092 : s(it+1)=s(it)
2828 482092 : h(it+1)=ninth*h(it) ! factor is due to step tripling in midpoint and even error series
2829 : end do
2830 0 : ABI_FREE(h)
2831 0 : ABI_FREE(s)
2832 :
2833 : case (7)
2834 : ! === Gauss-Legendre ===
2835 0 : NX0=5 ; if (PRESENT(npts)) NX0=npts
2836 0 : NX=NX0
2837 0 : do it=1,NT
2838 0 : ABI_MALLOC(wx,(NX))
2839 0 : ABI_MALLOC(xx,(NX))
2840 0 : call coeffs_gausslegint(xmin,xmax,xx,wx,NX)
2841 0 : quad=zero
2842 0 : do ix=1,NX
2843 0 : quad=quad+wx(ix)*func(xx(ix))
2844 : end do
2845 0 : ABI_FREE(wx)
2846 0 : ABI_FREE(xx)
2847 0 : if (it>1) then
2848 : !write(std_out,*) quad
2849 0 : if (ABS(quad-old_quad)<EPS*ABS(old_quad).or.(ABS(quad)<TOL.and.ABS(old_quad)<TOL)) RETURN
2850 : end if
2851 0 : old_quad=quad
2852 0 : NX=NX+NX0
2853 : !NX=2*NX
2854 : end do
2855 :
2856 : case default
2857 0 : write(msg,'(a,i3)')'Wrong value for qopt',qopt
2858 120516 : ABI_BUG(msg)
2859 : end select
2860 :
2861 : write(msg,'(a,i0,2(a,es14.6))')&
2862 0 : "Results are not converged within the given accuracy. ntrial= ",NT,"; EPS= ",EPS,"; TOL= ",TOL
2863 0 : ABI_WARNING(msg)
2864 0 : ierr = -1
2865 :
2866 : end subroutine quadrature
2867 : !!***
2868 :
2869 : !!****f* m_numeric_tools/ctrap
2870 : !! NAME
2871 : !! ctrap
2872 : !!
2873 : !! FUNCTION
2874 : !! Do corrected trapezoidal integral on uniform grid of spacing hh.
2875 : !!
2876 : !! INPUTS
2877 : !! imax=highest index of grid=grid point number of upper limit
2878 : !! ff(imax)=integrand values
2879 : !! hh=spacing between x points
2880 : !!
2881 : !! OUTPUT
2882 : !! ans=resulting integral by corrected trapezoid
2883 : !!
2884 : !! NOTES
2885 : !!
2886 : !! SOURCE
2887 :
2888 14955751 : subroutine ctrap(imax,ff,hh,ans)
2889 :
2890 : !Arguments ------------------------------------
2891 : !scalars
2892 : integer,intent(in) :: imax
2893 : real(dp),intent(in) :: hh
2894 : real(dp),intent(out) :: ans
2895 : !arrays
2896 : real(dp),intent(in) :: ff(imax)
2897 :
2898 : !Local variables-------------------------------
2899 : !scalars
2900 : integer :: ir,ir2
2901 : real(dp) :: endpt,sum
2902 : ! *************************************************************************
2903 :
2904 14955751 : if (imax>=10)then
2905 :
2906 : ! endpt=end point correction terms (low and high ends)
2907 : endpt = (23.75d0*(ff(1)+ff(imax )) &
2908 : & + 95.10d0*(ff(2)+ff(imax-1)) &
2909 : & + 55.20d0*(ff(3)+ff(imax-2)) &
2910 : & + 79.30d0*(ff(4)+ff(imax-3)) &
2911 14955735 : & + 70.65d0*(ff(5)+ff(imax-4)))/ 72.d0
2912 14955735 : ir2 = imax - 5
2913 14955735 : sum=0.00d0
2914 14955735 : if (ir2 > 5) then
2915 4439214466 : do ir=6,ir2
2916 4439214466 : sum = sum + ff(ir)
2917 : end do
2918 : end if
2919 14955735 : ans = (sum + endpt ) * hh
2920 :
2921 16 : else if (imax>=8)then
2922 : endpt = (17.0d0*(ff(1)+ff(imax )) &
2923 : & + 59.0d0*(ff(2)+ff(imax-1)) &
2924 : & + 43.0d0*(ff(3)+ff(imax-2)) &
2925 4 : & + 49.0d0*(ff(4)+ff(imax-3)) )/ 48.d0
2926 4 : sum=0.0d0
2927 4 : if(imax==9)sum=ff(5)
2928 4 : ans = (sum + endpt ) * hh
2929 :
2930 : else if (imax==7)then
2931 : ans = (17.0d0*(ff(1)+ff(imax )) &
2932 : & + 59.0d0*(ff(2)+ff(imax-1)) &
2933 : & + 43.0d0*(ff(3)+ff(imax-2)) &
2934 2 : & + 50.0d0* ff(4) )/ 48.d0 *hh
2935 :
2936 : else if (imax==6)then
2937 : ans = (17.0d0*(ff(1)+ff(imax )) &
2938 : & + 59.0d0*(ff(2)+ff(imax-1)) &
2939 2 : & + 44.0d0*(ff(3)+ff(imax-2)) )/ 48.d0 *hh
2940 :
2941 : else if (imax==5)then
2942 : ans = ( (ff(1)+ff(5)) &
2943 : & + four*(ff(2)+ff(4)) &
2944 2 : & + two * ff(3) )/ three *hh
2945 :
2946 : else if (imax==4)then
2947 : ans = (three*(ff(1)+ff(4)) &
2948 2 : & + nine *(ff(2)+ff(3)) )/ eight *hh
2949 :
2950 : else if (imax==3)then
2951 : ans = ( (ff(1)+ff(3)) &
2952 2 : & + four* ff(2) )/ three *hh
2953 :
2954 : else if (imax==2)then
2955 2 : ans = (ff(1)+ff(2))/ two *hh
2956 :
2957 : else if (imax==1)then
2958 0 : ans = ff(1)*hh
2959 :
2960 : end if
2961 :
2962 14955751 : end subroutine ctrap
2963 : !!***
2964 :
2965 : !!****f* m_numeric_tools/cspint
2966 : !! NAME
2967 : !! cspint
2968 : !!
2969 : !! FUNCTION
2970 : !! Estimates the integral of a tabulated function.
2971 : !!
2972 : !! INPUTS
2973 : !!
2974 : !! OUTPUT
2975 : !!
2976 : !! NOTES
2977 : !!
2978 : !! The routine is given the value of a function F(X) at a set of
2979 : !! nodes XTAB, and estimates
2980 : !!
2981 : !! Integral ( A <= X <= B ) F(X) DX
2982 : !!
2983 : !! by computing the cubic natural spline S(X) that interpolates
2984 : !! F(X) at the nodes, and then computing
2985 : !!
2986 : !! Integral ( A <= X <= B ) S(X) DX
2987 : !!
2988 : !! exactly.
2989 : !!
2990 : !! Other output from the program includes the definite integral
2991 : !! from X(1) to X(I) of S(X), and the coefficients necessary for
2992 : !! the user to evaluate the spline S(X) at any point.
2993 : !!
2994 : !! Modified:
2995 : !!
2996 : !! 30 October 2000
2997 : !!
2998 : !! Reference:
2999 : !!
3000 : !! Philip Davis and Philip Rabinowitz,
3001 : !! Methods of Numerical Integration,
3002 : !! Blaisdell Publishing, 1967.
3003 : !!
3004 : !! Parameters:
3005 : !!
3006 : !! Input, real (dp) FTAB(NTAB), contains the tabulated values of
3007 : !! the function, FTAB(I) = F(XTAB(I)).
3008 : !!
3009 : !! Input, real (dp) XTAB(NTAB), contains the points at which the
3010 : !! function was evaluated. The XTAB's must be distinct and
3011 : !! in ascending order.
3012 : !!
3013 : !! Input, integer NTAB, the number of entries in FTAB and
3014 : !! XTAB. NTAB must be at least 3.
3015 : !!
3016 : !! Input, real (dp) A, lower limit of integration.
3017 : !!
3018 : !! Input, real (dp) B, upper limit of integration.
3019 : !!
3020 : !! Output, real (dp) Y(3,NTAB), will contain the coefficients
3021 : !! of the interpolating natural spline over each subinterval.
3022 : !!
3023 : !! For XTAB(I) <= X <= XTAB(I+1),
3024 : !!
3025 : !! S(X) = FTAB(I) + Y(1,I)*(X-XTAB(I))
3026 : !! + Y(2,I)*(X-XTAB(I))**2
3027 : !! + Y(3,I)*(X-XTAB(I))**3
3028 : !!
3029 : !! Output, real (dp) E(NTAB), E(I) = the definite integral from
3030 : !! XTAB(1) to XTAB(I) of S(X).
3031 : !!
3032 : !! Workspace, real (dp) WORK(NTAB).
3033 : !!
3034 : !! Output, real (dp) RESULT, the estimated value of the integral.
3035 : !!
3036 : !!
3037 : !! SOURCE
3038 :
3039 0 : subroutine cspint ( ftab, xtab, ntab, a, b, y, e, work, result )
3040 :
3041 : !Arguments ------------------------------------
3042 : !scalars
3043 : integer, intent(in) :: ntab
3044 : real(dp), intent(in) :: a
3045 : real(dp), intent(in) :: b
3046 : real(dp), intent(inout) :: e(ntab)
3047 : real(dp), intent(in) :: ftab(ntab)
3048 : real(dp), intent(inout) :: work(ntab)
3049 : real(dp), intent(in) :: xtab(ntab)
3050 : real(dp), intent(inout) :: y(3,ntab)
3051 : real(dp), intent(out) :: result
3052 :
3053 : !Local variables ------------------------------
3054 : !scalars
3055 : integer :: i
3056 : integer :: j
3057 : real(dp) :: r
3058 : real(dp) :: s
3059 : real(dp) :: term
3060 : real(dp) :: u
3061 : !************************************************************************
3062 :
3063 0 : if ( ntab < 3 ) then
3064 0 : write(std_out,'(a)' ) ' '
3065 0 : write(std_out,'(a)' ) 'CSPINT - Fatal error!'
3066 0 : write(std_out,'(a,i6)' ) ' NTAB must be at least 3, but input NTAB = ',ntab
3067 0 : ABI_ERROR("Aborting now")
3068 : end if
3069 :
3070 0 : do i = 1, ntab-1
3071 :
3072 0 : if ( xtab(i+1) <= xtab(i) ) then
3073 0 : write(std_out,'(a)' ) ' '
3074 0 : write(std_out,'(a)' ) 'CSPINT - Fatal error!'
3075 0 : write(std_out,'(a)' ) ' Nodes not in strict increasing order.'
3076 0 : write(std_out,'(a,i6)' ) ' XTAB(I) <= XTAB(I-1) for I=',i
3077 0 : write(std_out,'(a,g14.6)' ) ' XTAB(I) = ',xtab(i)
3078 0 : write(std_out,'(a,g14.6)' ) ' XTAB(I-1) = ',xtab(i-1)
3079 0 : ABI_ERROR("Aborting now")
3080 : end if
3081 :
3082 : end do
3083 :
3084 0 : s = zero
3085 0 : do i = 1, ntab-1
3086 0 : r = ( ftab(i+1) - ftab(i) ) / ( xtab(i+1) - xtab(i) )
3087 0 : y(2,i) = r - s
3088 0 : s = r
3089 : end do
3090 :
3091 0 : result = zero
3092 0 : s = zero
3093 0 : r = zero
3094 0 : y(2,1) = zero
3095 0 : y(2,ntab) = zero
3096 :
3097 0 : do i = 2, ntab-1
3098 0 : y(2,i) = y(2,i) + r * y(2,i-1)
3099 0 : work(i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s
3100 0 : s = xtab(i+1) - xtab(i)
3101 0 : r = s / work(i)
3102 : end do
3103 :
3104 0 : do j = 2, ntab-1
3105 0 : i = ntab+1-j
3106 0 : y(2,i) = ( ( xtab(i+1) - xtab(i) ) * y(2,i+1) - y(2,i) ) / work(i)
3107 : end do
3108 :
3109 0 : do i = 1, ntab-1
3110 0 : s = xtab(i+1) - xtab(i)
3111 0 : r = y(2,i+1) - y(2,i)
3112 0 : y(3,i) = r / s
3113 0 : y(2,i) = three * y(2,i)
3114 0 : y(1,i) = ( ftab(i+1) - ftab(i) ) / s - ( y(2,i) + r ) * s
3115 : end do
3116 :
3117 0 : e(1) = 0.0D+00
3118 0 : do i = 1, ntab-1
3119 0 : s = xtab(i+1)-xtab(i)
3120 : term = ((( y(3,i) * quarter * s + y(2,i) * third ) * s &
3121 0 : + y(1,i) * half ) * s + ftab(i) ) * s
3122 0 : e(i+1) = e(i) + term
3123 : end do
3124 : !
3125 : ! Determine where the endpoints A and B lie in the mesh of XTAB's.
3126 : !
3127 0 : r = a
3128 0 : u = one
3129 :
3130 0 : do j = 1, 2
3131 : !
3132 : ! The endpoint is less than or equal to XTAB(1).
3133 : !
3134 0 : if ( r <= xtab(1) ) then
3135 0 : result = result-u*((r-xtab(1))*y(1,1)*half +ftab(1))*(r-xtab(1))
3136 : !
3137 : ! The endpoint is greater than or equal to XTAB(NTAB).
3138 : !
3139 0 : else if ( xtab(ntab) <= r ) then
3140 :
3141 : result = result -u * ( e(ntab) + ( r - xtab(ntab) ) &
3142 : * ( ftab(ntab) + half * ( ftab(ntab-1) &
3143 : + ( xtab(ntab) - xtab(ntab-1) ) * y(1,ntab-1) ) &
3144 0 : * ( r - xtab(ntab) )))
3145 : !
3146 : ! The endpoint is strictly between XTAB(1) and XTAB(NTAB).
3147 : !
3148 : else
3149 :
3150 0 : do i = 1, ntab-1
3151 :
3152 0 : if ( r <= xtab(i+1) ) then
3153 0 : r = r-xtab(i)
3154 : result = result-u*(e(i)+(((y(3,i)*quarter*r+y(2,i)*third)*r &
3155 0 : +y(1,i)*half )*r+ftab(i))*r)
3156 0 : go to 120
3157 : end if
3158 :
3159 : end do
3160 :
3161 : end if
3162 :
3163 : 120 continue
3164 :
3165 0 : u = -one
3166 0 : r = b
3167 :
3168 : end do
3169 :
3170 0 : end subroutine cspint
3171 : !!***
3172 :
3173 : !!****f* m_numeric_tools/coeffs_gausslegint
3174 : !! NAME
3175 : !! coeffs_gausslegint
3176 : !!
3177 : !! FUNCTION
3178 : !! Compute the coefficients (supports and weights) for Gauss-Legendre integration.
3179 : !! Inspired by a routine due to G. Rybicki.
3180 : !!
3181 : !! INPUTS
3182 : !! xmin=lower bound of integration
3183 : !! xmax=upper bound of integration
3184 : !! n=order of integration
3185 : !!
3186 : !! OUTPUT
3187 : !! x(n)=array of support points
3188 : !! weights(n)=array of integration weights
3189 : !!
3190 : !! SOURCE
3191 :
3192 98 : subroutine coeffs_gausslegint(xmin,xmax,x,weights,n)
3193 :
3194 : !Arguments ------------------------------------
3195 : !scalars
3196 : integer,intent(in) :: n
3197 : real(dp),intent(in) :: xmin,xmax
3198 : real(dp),intent(out) :: x(n),weights(n)
3199 :
3200 : !Local variables ------------------------------
3201 : !scalars
3202 : integer :: i,j
3203 : real(dp),parameter :: tol=1.d-13
3204 : real(dp),parameter :: pi=4.d0*atan(1.d0)
3205 : real(dp) :: z,z1,xmean,p1,p2,p3,pp,xl
3206 : !************************************************************************
3207 :
3208 98 : xl=(xmax-xmin)*0.5d0
3209 98 : xmean=(xmax+xmin)*0.5d0
3210 :
3211 1598 : do i=1,(n+1)/2
3212 1500 : z=cos(pi*(i-0.25d0)/(n+0.5d0))
3213 :
3214 : do
3215 4589 : p1=1.d0
3216 4589 : p2=0.d0
3217 :
3218 392920 : do j=1,n
3219 388331 : p3=p2
3220 388331 : p2=p1
3221 392920 : p1=((2.d0*j - 1.d0)*z*p2 - (j-1.d0)*p3)/j
3222 : end do
3223 :
3224 4589 : pp=n*(p2-z*p1)/(1.0d0-z**2)
3225 4589 : z1=z
3226 4589 : z=z1-p1/pp
3227 :
3228 4589 : if(abs(z-z1) < tol) exit
3229 : end do
3230 :
3231 1500 : x(i)=xmean-xl*z
3232 1500 : x(n+1-i)=xmean+xl*z
3233 1500 : weights(i)=2.d0*xl/((1.d0-z**2)*pp**2)
3234 1598 : weights(n+1-i)=weights(i)
3235 : end do
3236 :
3237 98 : end subroutine coeffs_gausslegint
3238 : !!***
3239 :
3240 : !----------------------------------------------------------------------
3241 :
3242 : !!****f* m_numeric_tools/simpson_cplx
3243 : !! NAME
3244 : !! simpson_cplx
3245 : !!
3246 : !! FUNCTION
3247 : !! Integrate a complex function using extended Simpson's rule.
3248 : !!
3249 : !! INPUTS
3250 : !! npts=Number of points.
3251 : !! step=Step of the mesh.
3252 : !! ff(npts)=Values of the integrand.
3253 : !!
3254 : !! OUTPUT
3255 : !! simpson_cplx=Integral of ff.
3256 : !!
3257 : !! NOTES
3258 : !! If npts is odd, the integration is done with the extended Simpson's rule (error = O^(step^4))
3259 : !! If npts is even, the last 4 four points are integrated separately via Simpson's 3/8 rule. Error = O(step^5)
3260 : !! while the first npts-3 points are integrared with the extended Simpson's rule.
3261 : !!
3262 : !! SOURCE
3263 :
3264 42720 : function simpson_cplx(npts,step,ff)
3265 :
3266 : !Arguments ------------------------------------
3267 : !scalars
3268 : integer,intent(in) :: npts
3269 : real(dp),intent(in) :: step
3270 : complex(dp),intent(in) :: ff(npts)
3271 : complex(dp) :: simpson_cplx
3272 :
3273 : !Local variables ------------------------------
3274 : !scalars
3275 : integer :: ii,my_n
3276 : complex(dp) :: sum_even, sum_odd
3277 : !************************************************************************
3278 :
3279 42720 : my_n=npts; if ((npts/2)*2 == npts) my_n=npts-3
3280 :
3281 42720 : if (my_n<2) then
3282 0 : ABI_ERROR("Too few points")
3283 : end if
3284 :
3285 42720 : sum_odd=czero
3286 42720 : do ii=2,my_n-1,2
3287 22929972 : sum_odd = sum_odd + ff(ii)
3288 : end do
3289 :
3290 42720 : sum_even=zero
3291 42720 : do ii=3,my_n-2,2
3292 22887252 : sum_even = sum_even + ff(ii)
3293 : end do
3294 :
3295 : ! Eq 25.4.6 Abramowitz. Error is O(step^4)
3296 42720 : simpson_cplx = step/three * (ff(1) + four*sum_odd + two*sum_even + ff(my_n))
3297 :
3298 42720 : if (my_n/=npts) then ! Simpson's 3/8 rule. Eq 25.4.13 Abramowitz. Error is O(step^5)
3299 0 : simpson_cplx = simpson_cplx + three*step/eight * (ff(npts-3) + 3*ff(npts-2) + 3*ff(npts-1) + ff(npts))
3300 : end if
3301 :
3302 42720 : end function simpson_cplx
3303 : !!***
3304 :
3305 : !----------------------------------------------------------------------
3306 :
3307 : !!****f* m_numeric_tools/hermitianize_spc
3308 : !! NAME
3309 : !! hermitianize_spc
3310 : !!
3311 : !! FUNCTION
3312 : !! Force a square matrix to be hermitian
3313 : !!
3314 : !! INPUTS
3315 : !! uplo=String describing which part of the matrix has been calculated.
3316 : !! Only the first character is tested (no case sensitive). Possible values are:
3317 : !! "All"= Full matrix is supplied in input
3318 : !! "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
3319 : !! "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.
3320 : !!
3321 : !! OUTPUT
3322 : !! (see side effects)
3323 : !!
3324 : !! SIDE EFFECTS
3325 : !! mat(:,:)=complex input matrix, hermitianized at output
3326 : !!
3327 : !! SOURCE
3328 :
3329 0 : subroutine hermitianize_spc(mat,uplo)
3330 :
3331 : !Arguments ------------------------------------
3332 : !scalars
3333 : character(len=*),intent(in) :: uplo
3334 : !arrays
3335 : complex(sp),intent(inout) :: mat(:,:)
3336 :
3337 : !Local variables-------------------------------
3338 : !scalars
3339 : integer :: nn,ii,jj
3340 : !arrays
3341 0 : complex(sp),allocatable :: tmp(:)
3342 : ! *************************************************************************
3343 :
3344 0 : nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3345 :
3346 0 : select case (uplo(1:1))
3347 :
3348 : case ("A","a") ! Full matrix has been calculated.
3349 0 : ABI_MALLOC(tmp,(nn))
3350 0 : do ii=1,nn
3351 0 : do jj=ii,nn
3352 : ! reference half constant is dp not sp
3353 0 : tmp(jj)=real(half)*(mat(ii,jj)+CONJG(mat(jj,ii)))
3354 : end do
3355 0 : mat(ii,ii:nn)=tmp(ii:nn)
3356 0 : mat(ii:nn,ii)=CONJG(tmp(ii:nn))
3357 : end do
3358 0 : ABI_FREE(tmp)
3359 :
3360 : case ("U","u") ! Only the upper triangle is used.
3361 0 : do jj=1,nn
3362 0 : do ii=1,jj
3363 0 : if (ii/=jj) then
3364 0 : mat(jj,ii) = CONJG(mat(ii,jj))
3365 : else
3366 0 : mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3367 : end if
3368 : end do
3369 : end do
3370 :
3371 : case ("L","l") ! Only the lower triangle is used.
3372 0 : do jj=1,nn
3373 0 : do ii=1,jj
3374 0 : if (ii/=jj) then
3375 0 : mat(ii,jj) = CONJG(mat(jj,ii))
3376 : else
3377 0 : mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),0.0_sp)
3378 : end if
3379 : end do
3380 : end do
3381 :
3382 : case default
3383 0 : ABI_ERROR("Wrong uplo"//TRIM(uplo))
3384 : end select
3385 :
3386 0 : end subroutine hermitianize_spc
3387 : !!***
3388 :
3389 : !----------------------------------------------------------------------
3390 :
3391 : !!****f* m_numeric_tools/hermitianize_dpc
3392 : !! NAME
3393 : !! hermitianize_dpc
3394 : !!
3395 : !! FUNCTION
3396 : !! Force a square matrix to be hermitian
3397 : !!
3398 : !! INPUTS
3399 : !! uplo=String describing which part of the matrix has been calculated.
3400 : !! Only the first character is tested (no case sensitive). Possible values are:
3401 : !! "All"= Full matrix is supplied in input
3402 : !! "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
3403 : !! "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.
3404 : !!
3405 : !! OUTPUT
3406 : !! (see side effects)
3407 : !!
3408 : !! SIDE EFFECTS
3409 : !! mat(:,:)=complex input matrix, hermitianized in output
3410 : !!
3411 : !! SOURCE
3412 :
3413 435 : subroutine hermitianize_dpc(mat,uplo)
3414 :
3415 : !Arguments ------------------------------------
3416 : !scalars
3417 : character(len=*),intent(in) :: uplo
3418 : !arrays
3419 : complex(dp),intent(inout) :: mat(:,:)
3420 :
3421 : !Local variables-------------------------------
3422 : !scalars
3423 : integer :: nn,ii,jj
3424 : !arrays
3425 435 : complex(dp),allocatable :: tmp(:)
3426 : ! *************************************************************************
3427 :
3428 435 : nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3429 :
3430 32 : select case (uplo(1:1))
3431 :
3432 : case ("A","a") ! Full matrix has been calculated.
3433 96 : ABI_MALLOC(tmp,(nn))
3434 656 : do ii=1,nn
3435 7560 : do jj=ii,nn
3436 7560 : tmp(jj)=half*(mat(ii,jj)+DCONJG(mat(jj,ii)))
3437 : end do
3438 7560 : mat(ii,ii:nn)=tmp(ii:nn)
3439 7592 : mat(ii:nn,ii)=DCONJG(tmp(ii:nn))
3440 : end do
3441 32 : ABI_FREE(tmp)
3442 :
3443 : case ("U","u") ! Only the upper triangle is used.
3444 5219 : do jj=1,nn
3445 324591 : do ii=1,jj
3446 324188 : if (ii/=jj) then
3447 314556 : mat(jj,ii) = DCONJG(mat(ii,jj))
3448 : else
3449 4816 : mat(ii,ii) = CMPLX(DBLE(mat(ii,ii)),zero, kind=dp)
3450 : end if
3451 : end do
3452 : end do
3453 :
3454 : case ("L","l") ! Only the lower triangle is used.
3455 0 : do jj=1,nn
3456 0 : do ii=1,jj
3457 0 : if (ii/=jj) then
3458 0 : mat(ii,jj) = DCONJG(mat(jj,ii))
3459 : else
3460 0 : mat(ii,ii) = CMPLX(REAL(mat(ii,ii)),zero, kind=dp)
3461 : end if
3462 : end do
3463 : end do
3464 :
3465 : case default
3466 435 : ABI_ERROR("Wrong uplo"//TRIM(uplo))
3467 : end select
3468 :
3469 435 : end subroutine hermitianize_dpc
3470 : !!***
3471 :
3472 : !----------------------------------------------------------------------
3473 :
3474 : !!****f* m_numeric_tools/mkherm
3475 : !! NAME
3476 : !! mkherm
3477 : !!
3478 : !! FUNCTION
3479 : !! Make the complex array(ndim,ndim) hermitian,
3480 : !! by adding half of it to its hermitian conjugate.
3481 : !!
3482 : !! INPUTS
3483 : !! ndim=dimension of the matrix
3484 : !! array= complex matrix
3485 : !!
3486 : !! SIDE EFFECTS
3487 : !! array= hermitian matrix made by adding half of array to its hermitian conjugate
3488 : !!
3489 : !! SOURCE
3490 :
3491 1071417 : pure subroutine mkherm(array,ndim)
3492 :
3493 : !Arguments -------------------------------
3494 : !scalars
3495 : integer,intent(in) :: ndim
3496 : !arrays
3497 : real(dp),intent(inout) :: array(2,ndim,ndim)
3498 :
3499 : !Local variables -------------------------
3500 : !scalars
3501 : integer :: i1,i2
3502 : ! *********************************************************************
3503 :
3504 7837947 : do i1=1,ndim
3505 43353192 : do i2=1,i1
3506 35515245 : array(1,i1,i2)=(array(1,i1,i2)+array(1,i2,i1))*half
3507 35515245 : array(2,i1,i2)=(array(2,i1,i2)-array(2,i2,i1))*half
3508 35515245 : array(1,i2,i1)=array(1,i1,i2)
3509 42281775 : array(2,i2,i1)=-array(2,i1,i2)
3510 : end do
3511 : end do
3512 :
3513 1071417 : end subroutine mkherm
3514 : !!***
3515 :
3516 : !----------------------------------------------------------------------
3517 :
3518 : !!****f* m_numeric_tools/hermit
3519 : !! NAME
3520 : !! hermit
3521 : !!
3522 : !! FUNCTION
3523 : !! Take a matrix in hermitian storage mode (lower triangle stored)
3524 : !! and redefine diagonal elements to impose Hermiticity
3525 : !! (diagonal terms have to be real).
3526 : !! If abs(Im(H(i,i)))>4096*machine precision, print error warning.
3527 : !! (Typical 64 bit machine precision is 2^-52 or 2.22e-16)
3528 : !!
3529 : !! INPUTS
3530 : !! chmin(n*n+n)=complex hermitian matrix with numerical noise possibly
3531 : !! rendering Im(diagonal elements) approximately 1e-15 or so
3532 : !! ndim=size of complex hermitian matrix
3533 : !!
3534 : !! OUTPUT
3535 : !! chmout(n*n+n)=redefined matrix with strictly real diagonal elements.
3536 : !! May be same storage location as chmin.
3537 : !! ierr=0 if no problem, 1 if the imaginary part of some element
3538 : !! too large (at present, stop in this case).
3539 : !!
3540 : !! TODO
3541 : !! Name is misleading, perhaps hermit_force_diago?
3542 : !! Interface allows aliasing
3543 : !!
3544 : !! SOURCE
3545 :
3546 445671 : subroutine hermit(chmin, chmout, ierr, ndim)
3547 :
3548 : !Arguments ------------------------------------
3549 : !scalars
3550 : integer,intent(in) :: ndim
3551 : integer,intent(out) :: ierr
3552 : !arrays
3553 : real(dp),intent(inout) :: chmin(ndim*ndim+ndim)
3554 : real(dp),intent(inout) :: chmout(ndim*ndim+ndim)
3555 :
3556 : !Local variables-------------------------------
3557 : !scalars
3558 : integer,save :: mmesgs=20,nmesgs=0
3559 : integer :: idim,max_errors,nerrors
3560 : real(dp),parameter :: eps=epsilon(0.0d0)
3561 : real(dp) :: ch_im,ch_re,modules,tol
3562 : character(len=500) :: msg
3563 : ! *************************************************************************
3564 :
3565 445671 : tol=4096.0d0*eps
3566 :
3567 445671 : ierr=0
3568 445671 : max_errors=0
3569 :
3570 : !Copy matrix into possibly new location
3571 62007609 : chmout(:)=chmin(:)
3572 :
3573 : !Loop over diagonal elements of matrix (off-diag not altered)
3574 3634990 : do idim=1,ndim
3575 :
3576 3189319 : ch_im=chmout(idim*idim+idim )
3577 3189319 : ch_re=chmout(idim*idim+idim-1)
3578 :
3579 : ! check for large absolute Im part and print warning when
3580 : ! larger than (some factor)*(machine precision)
3581 3189319 : nerrors=0
3582 3189319 : if( abs(ch_im) > tol .and. abs(ch_im) > tol8*abs(ch_re)) nerrors=2
3583 3189319 : if( abs(ch_im) > tol .or. abs(ch_im) > tol8*abs(ch_re)) nerrors=1
3584 :
3585 3189319 : if( (abs(ch_im) > tol .and. nmesgs<mmesgs) .or. nerrors==2)then
3586 : write(msg, '(3a,i0,a,es20.12,a,es20.12,a)' )&
3587 0 : ' Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
3588 0 : ' for component: ',idim,' Im part is: ',ch_im,', Re part is: ',ch_re,'.'
3589 0 : call wrtout(std_out,msg)
3590 0 : nmesgs=nmesgs+1
3591 : end if
3592 :
3593 3189319 : if( ( abs(ch_im) > tol8*abs(ch_re) .and. nmesgs<mmesgs) .or. nerrors==2)then
3594 : write(msg, '(3a,i0,a,es20.12,a,es20.12,a)' )&
3595 0 : ' Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
3596 0 : ' for component: ',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
3597 0 : call wrtout(std_out,msg)
3598 0 : nmesgs=nmesgs+1
3599 : end if
3600 :
3601 : ! compute modulus $= (\Re^2+\Im^2)^{1/2}$
3602 3189319 : modules=sqrt(ch_re**2+ch_im**2)
3603 :
3604 : ! set Re part to modulus with sign of original Re part
3605 3189319 : chmout(idim*idim+idim-1)=sign(modules,ch_re)
3606 :
3607 : ! set Im part to 0
3608 3189319 : chmout(idim*idim+idim)=zero
3609 :
3610 3634990 : max_errors=max(max_errors,nerrors)
3611 : end do
3612 :
3613 445671 : if (max_errors==2)then
3614 0 : ierr=1
3615 : write(msg, '(3a)' )&
3616 0 : 'Imaginary part(s) of diagonal Hermitian matrix element(s) is too large.',ch10,&
3617 0 : 'See previous messages.'
3618 0 : ABI_BUG(msg)
3619 : end if
3620 :
3621 445671 : end subroutine hermit
3622 : !!***
3623 :
3624 : !----------------------------------------------------------------------
3625 :
3626 : !!****f* m_numeric_tools/symmetrize_spc
3627 : !! NAME
3628 : !! symmetrize_spc
3629 : !!
3630 : !! FUNCTION
3631 : !! Force a square matrix to be symmetric.
3632 : !!
3633 : !! INPUTS
3634 : !! uplo=String describing which part of the matrix has been calculated.
3635 : !! Only the first character is tested (no case sensitive). Possible values are:
3636 : !! "All"= Full matrix is supplied in input
3637 : !! "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
3638 : !! "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.
3639 : !!
3640 : !! OUTPUT
3641 : !! (see side effects)
3642 : !!
3643 : !! SIDE EFFECTS
3644 : !! mat(:,:)=complex input matrix, symmetrized at output
3645 : !!
3646 : !! SOURCE
3647 :
3648 0 : subroutine symmetrize_spc(mat,uplo)
3649 :
3650 : !Arguments ------------------------------------
3651 : !scalars
3652 : character(len=*),intent(in) :: uplo
3653 : !arrays
3654 : complex(sp),intent(inout) :: mat(:,:)
3655 :
3656 : !Local variables-------------------------------
3657 : !scalars
3658 : integer :: nn,ii,jj
3659 : !arrays
3660 0 : complex(sp),allocatable :: tmp(:)
3661 : ! *************************************************************************
3662 :
3663 0 : nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3664 :
3665 0 : select case (uplo(1:1))
3666 :
3667 : case ("A","a") ! Full matrix has been calculated.
3668 0 : ABI_MALLOC(tmp,(nn))
3669 0 : do ii=1,nn
3670 0 : do jj=ii,nn
3671 0 : tmp(jj)=REAL(half)*(mat(ii,jj)+mat(jj,ii))
3672 : end do
3673 0 : mat(ii,ii:nn)=tmp(ii:nn)
3674 0 : mat(ii:nn,ii)=tmp(ii:nn)
3675 : end do
3676 0 : ABI_FREE(tmp)
3677 :
3678 : case ("U","u") ! Only the upper triangle is used.
3679 0 : do jj=1,nn
3680 0 : do ii=1,jj-1
3681 0 : mat(jj,ii) = mat(ii,jj)
3682 : end do
3683 : end do
3684 :
3685 : case ("L","l") ! Only the lower triangle is used.
3686 0 : do jj=1,nn
3687 0 : do ii=1,jj-1
3688 0 : mat(ii,jj) = mat(jj,ii)
3689 : end do
3690 : end do
3691 :
3692 : case default
3693 0 : ABI_ERROR("Wrong uplo"//TRIM(uplo))
3694 : end select
3695 :
3696 0 : end subroutine symmetrize_spc
3697 : !!***
3698 :
3699 : !----------------------------------------------------------------------
3700 :
3701 : !!****f* m_numeric_tools/symmetrize_dpc
3702 : !! NAME
3703 : !! symmetrize_dpc
3704 : !!
3705 : !! FUNCTION
3706 : !! Force a square matrix to be symmetric.
3707 : !!
3708 : !! INPUTS
3709 : !! uplo=String describing which part of the matrix has been calculated.
3710 : !! Only the first character is tested (no case sensitive). Possible values are:
3711 : !! "All"= Full matrix is supplied in input
3712 : !! "Upper"=Upper triangle is in input. Lower triangle is reconstructed by symmetry.
3713 : !! "Lower"=Lower triangle is in input. Upper triangle is reconstructed by symmetry.
3714 : !!
3715 : !! OUTPUT
3716 : !! (see side effects)
3717 : !!
3718 : !! SIDE EFFECTS
3719 : !! mat(:,:)=complex input matrix, symmetrized in output
3720 : !!
3721 : !! SOURCE
3722 :
3723 0 : subroutine symmetrize_dpc(mat, uplo)
3724 :
3725 : !Arguments ------------------------------------
3726 : !scalars
3727 : character(len=*),intent(in) :: uplo
3728 : !arrays
3729 : complex(dp),intent(inout) :: mat(:,:)
3730 :
3731 : !Local variables-------------------------------
3732 : !scalars
3733 : integer :: nn,ii,jj
3734 : !arrays
3735 0 : complex(dp),allocatable :: tmp(:)
3736 : ! *************************************************************************
3737 :
3738 0 : nn = assert_eq(SIZE(mat,1),SIZE(mat,2),'Matrix not square',__FILE__,__LINE__)
3739 :
3740 0 : select case (uplo(1:1))
3741 : case ("A","a") ! Full matrix has been calculated.
3742 0 : ABI_MALLOC(tmp,(nn))
3743 0 : do ii=1,nn
3744 0 : do jj=ii,nn
3745 0 : tmp(jj)=half*(mat(ii,jj)+mat(jj,ii))
3746 : end do
3747 0 : mat(ii,ii:nn)=tmp(ii:nn)
3748 0 : mat(ii:nn,ii)=tmp(ii:nn)
3749 : end do
3750 0 : ABI_FREE(tmp)
3751 :
3752 : case ("U","u") ! Only the upper triangle is used.
3753 0 : do jj=1,nn
3754 0 : do ii=1,jj-1
3755 0 : mat(jj,ii) = mat(ii,jj)
3756 : end do
3757 : end do
3758 :
3759 : case ("L","l") ! Only the lower triangle is used.
3760 0 : do jj=1,nn
3761 0 : do ii=1,jj-1
3762 0 : mat(ii,jj) = mat(jj,ii)
3763 : end do
3764 : end do
3765 :
3766 : case default
3767 0 : ABI_ERROR("Wrong uplo"//TRIM(uplo))
3768 : end select
3769 :
3770 0 : end subroutine symmetrize_dpc
3771 : !!***
3772 :
3773 : !!****f* m_numeric_tools/pack_matrix
3774 : !! NAME
3775 : !! pack_matrix
3776 : !!
3777 : !! FUNCTION
3778 : !! Packs a matrix into hermitian format
3779 : !!
3780 : !! INPUTS
3781 : !! N: size of matrix
3782 : !! cplx: 2 if matrix is complex, 1 for real matrix.
3783 : !! mat_in(cplx, N*N)= matrix to be packed
3784 : !!
3785 : !! OUTPUT
3786 : !! mat_out(cplx*N*N+1/2)= packed matrix (upper triangle)
3787 : !!
3788 : !! SOURCE
3789 :
3790 593 : subroutine pack_matrix(mat_in, mat_out, N, cplx)
3791 :
3792 : !Arguments ------------------------------------
3793 : !scalars
3794 : integer, intent(in) :: N, cplx
3795 : real(dp), intent(in) :: mat_in(cplx, N*N)
3796 : real(dp), intent(out) :: mat_out(cplx*N*(N+1)/2)
3797 :
3798 : !Local variables-------------------------------
3799 : integer :: isubh, i, j
3800 : ! *************************************************************************
3801 :
3802 593 : isubh = 1
3803 10173 : do j=1,N
3804 120111 : do i=1,j
3805 109938 : mat_out(isubh) = mat_in(1, (j-1)*N+i)
3806 : ! bad for vectorization, but it's not performance critical, so ...
3807 109938 : if(cplx == 2) then
3808 109938 : mat_out(isubh+1) = mat_in(2, (j-1)*N+i)
3809 : end if
3810 119518 : isubh=isubh+cplx
3811 : end do
3812 : end do
3813 :
3814 593 : end subroutine pack_matrix
3815 : !!***
3816 :
3817 :
3818 : !!****f* m_numeric_tools/check_vec_conjg
3819 : !! NAME
3820 : !! check_vec_conjg
3821 : !!
3822 : !! FUNCTION
3823 : !! Test whether two complex vectors `vec1` and `vec2` of size `nn` are conjugate of each other.
3824 : !! within an optional tolerance abs_tol (default: 1e-6).
3825 : !! Return max absolute difference for real and imag part in abs_diff(2) and exit status in ierr.
3826 : !!
3827 : !! SOURCE
3828 :
3829 0 : integer function check_vec_conjg(nn, vec1, vec2, abs_diff, abs_tol) result(ierr)
3830 :
3831 : !Arguments ------------------------------------
3832 : integer, intent(in) :: nn
3833 : complex(dp), intent(in) :: vec1(nn), vec2(nn)
3834 : real(dp),intent(out) :: abs_diff(2)
3835 : real(dp),optional,intent(in) :: abs_tol
3836 :
3837 : !Local variables-------------------------------
3838 : integer :: ii
3839 : real(dp) :: my_abs_tol
3840 : ! *************************************************************************
3841 :
3842 0 : my_abs_tol = tol6; if (present(abs_tol)) my_abs_tol = abs_tol
3843 0 : ierr = 0
3844 0 : abs_diff = zero
3845 :
3846 0 : do ii=1,nn
3847 0 : abs_diff(1) = max(abs_diff(1), abs(real(vec1(ii)) - real(vec2(ii))))
3848 0 : abs_diff(2) = max(abs_diff(2), abs(aimag(vec1(ii)) + aimag(vec2(ii))))
3849 : end do
3850 :
3851 0 : if (any(abs_diff > abs_tol)) ierr = 1
3852 :
3853 0 : end function check_vec_conjg
3854 : !!***
3855 :
3856 : !----------------------------------------------------------------------
3857 :
3858 : !!****f* m_numeric_tools/print_arr1d_spc
3859 : !! NAME
3860 : !! print_arr1d_spc
3861 : !!
3862 : !! FUNCTION
3863 : !! Print an array using a nice (?) format
3864 : !!
3865 : !! INPUTS
3866 : !! arr(:)=vector/matrix to be printed
3867 : !! units(:)=unit numbers
3868 : !! max_r,max_c(optional)=Max number of rows and columns to be printed
3869 : !! (DEFAULT is 9, output format assumes to be less that 99, but there might be
3870 : !! problems with wrtout if message size exceeds 500 thus max number of elements should be ~60)
3871 : !!
3872 : !! OUTPUT
3873 : !! (only printing)
3874 : !!
3875 : !! SOURCE
3876 :
3877 0 : subroutine print_arr1d_spc(units, arr, max_r)
3878 :
3879 : !Arguments ------------------------------------
3880 : !scalars
3881 : integer,intent(in) :: units(:)
3882 : integer,optional,intent(in) :: max_r
3883 : !arrays
3884 : complex(sp),intent(in) :: arr(:)
3885 :
3886 : !Local variables-------------------------------
3887 : !scalars
3888 : integer :: ii, nr, mr
3889 : character(len=500) :: msg
3890 : character(len=100) :: fmth,fmt1
3891 : ! *************************************************************************
3892 :
3893 0 : mr=15 ; if (PRESENT(max_r )) mr=max_r
3894 0 : nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3895 :
3896 0 : write(fmth,*)'(6x,',mr,'(i2,6x))'
3897 0 : write(fmt1,*)'(3x,',mr,'f8.3)'
3898 :
3899 0 : write(msg,fmth)(ii,ii=1,mr)
3900 0 : call wrtout(units, msg) !header
3901 0 : write(msg,fmt1)REAL (arr(1:mr))
3902 0 : call wrtout(units, msg) !real part
3903 0 : write(msg,fmt1)AIMAG(arr(1:mr))
3904 0 : call wrtout(units, msg) !imag part
3905 :
3906 0 : end subroutine print_arr1d_spc
3907 : !!***
3908 :
3909 : !----------------------------------------------------------------------
3910 :
3911 : !!****f* m_numeric_tools/print_arr1d_dpc
3912 : !! NAME
3913 : !! print_arr1d_dpc
3914 : !!
3915 : !! FUNCTION
3916 : !!
3917 : !! INPUTS
3918 : !!
3919 : !! OUTPUT
3920 : !!
3921 : !! SOURCE
3922 :
3923 1820 : subroutine print_arr1d_dpc(units, arr, max_r)
3924 :
3925 : !Arguments ------------------------------------
3926 : !scalars
3927 : integer,intent(in) :: units(:)
3928 : integer,optional,intent(in) :: max_r
3929 : !arrays
3930 : complex(dp),intent(in) :: arr(:)
3931 :
3932 : !Local variables-------------------------------
3933 : !scalars
3934 : integer :: ii,nr,mr
3935 : character(len=500) :: msg
3936 : character(len=100) :: fmth,fmt1
3937 : ! *************************************************************************
3938 :
3939 1820 : mr=15 ; if (PRESENT(max_r )) mr=max_r
3940 :
3941 : ! Print out matrix.
3942 1820 : nr=SIZE(arr,DIM=1) ; if (mr>nr) mr=nr
3943 :
3944 1820 : write(fmth,*)'(6x,',mr,'(i2,6x))'
3945 1820 : write(fmt1,*)'(3x,',mr,'f8.3)'
3946 :
3947 18200 : write(msg,fmth)(ii,ii=1,mr)
3948 1820 : call wrtout(units, msg) ! header
3949 18200 : write(msg,fmt1)REAL (arr(1:mr))
3950 1820 : call wrtout(units, msg) !real part
3951 18200 : write(msg,fmt1)AIMAG(arr(1:mr))
3952 1820 : call wrtout(units, msg) !imag part
3953 :
3954 1820 : end subroutine print_arr1d_dpc
3955 : !!***
3956 :
3957 : !----------------------------------------------------------------------
3958 :
3959 : !!****f* m_numeric_tools/print_arr2d_spc
3960 : !! NAME
3961 : !! print_arr2d_spc
3962 : !!
3963 : !! FUNCTION
3964 : !!
3965 : !! INPUTS
3966 : !!
3967 : !! OUTPUT
3968 : !!
3969 : !! SOURCE
3970 :
3971 0 : subroutine print_arr2d_spc(units, arr, max_r, max_c)
3972 :
3973 : !Arguments ------------------------------------
3974 : !scalars
3975 : integer,intent(in) :: units(:)
3976 : integer,optional,intent(in) :: max_r, max_c
3977 : !arrays
3978 : complex(sp),intent(in) :: arr(:,:)
3979 :
3980 : !Local variables-------------------------------
3981 : !scalars
3982 : integer :: ii,jj,nc,nr,mc,mr
3983 : character(len=500) :: msg
3984 : character(len=100) :: fmth,fmt1,fmt2
3985 : ! *************************************************************************
3986 :
3987 0 : mc =9 ; if (PRESENT(max_c )) mc =max_c
3988 0 : mr =9 ; if (PRESENT(max_r )) mr =max_r
3989 :
3990 : ! === Print out matrix ===
3991 0 : nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
3992 0 : nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
3993 :
3994 0 : write(fmth,*)'(6x,',mc,'(i2,6x))'
3995 0 : write(fmt1,*)'(3x,i2,',mc,'f8.3)'
3996 0 : write(fmt2,*)'(5x ,',mc,'f8.3,a)'
3997 :
3998 0 : write(msg,fmth)(jj,jj=1,mc)
3999 0 : call wrtout(units, msg) !header
4000 0 : do ii=1,mr
4001 0 : write(msg,fmt1)ii,REAL(arr(ii,1:mc))
4002 0 : call wrtout(units, msg) !real part
4003 0 : write(msg,fmt2) AIMAG(arr(ii,1:mc)),ch10
4004 0 : call wrtout(units, msg) !imag part
4005 : end do
4006 :
4007 0 : end subroutine print_arr2d_spc
4008 : !!***
4009 :
4010 : !----------------------------------------------------------------------
4011 :
4012 : !!****f* m_numeric_tools/print_arr2d_dpc
4013 : !! NAME
4014 : !! print_arr2d_dpc
4015 : !!
4016 : !! FUNCTION
4017 : !!
4018 : !! INPUTS
4019 : !!
4020 : !! OUTPUT
4021 : !!
4022 : !! SOURCE
4023 :
4024 9290 : subroutine print_arr2d_dpc(units, arr, max_r, max_c)
4025 :
4026 : !Arguments ------------------------------------
4027 : !scalars
4028 : integer,intent(in) :: units(:)
4029 : integer,optional,intent(in) :: max_r,max_c
4030 : !arrays
4031 : complex(dp),intent(in) :: arr(:,:)
4032 :
4033 : !Local variables-------------------------------
4034 : !scalars
4035 : integer :: ii,jj,nc,nr,mc,mr
4036 : character(len=500) :: msg
4037 : character(len=100) :: fmth,fmt1,fmt2
4038 : ! *************************************************************************
4039 :
4040 9290 : mc =9 ; if (PRESENT(max_c )) mc =max_c
4041 9290 : mr =9 ; if (PRESENT(max_r )) mr =max_r
4042 :
4043 : ! === Print out matrix ===
4044 9290 : nr=SIZE(arr,DIM=1); if (mr>nr) mr=nr
4045 9290 : nc=SIZE(arr,DIM=2); if (mc>nc) mc=nc
4046 :
4047 9290 : write(fmth,*)'(6x,',mc,'(i2,6x))'
4048 9290 : write(fmt1,*)'(3x,i2,',mc,'f8.3)'
4049 9290 : write(fmt2,*)'(5x ,',mc,'f8.3,a)'
4050 :
4051 92692 : write(msg,fmth)(jj,jj=1,mc)
4052 9290 : call wrtout(units, msg) ! header
4053 66168 : do ii=1,mr
4054 568324 : write(msg,fmt1)ii,REAL(arr(ii,1:mc))
4055 56878 : call wrtout(units, msg) ! real part
4056 568324 : write(msg,fmt2) AIMAG(arr(ii,1:mc)),ch10
4057 66168 : call wrtout(units, msg) ! imag part
4058 : end do
4059 :
4060 9290 : end subroutine print_arr2d_dpc
4061 : !!***
4062 :
4063 : !----------------------------------------------------------------------
4064 :
4065 : !!****f* m_numeric_tools/pade
4066 : !! NAME
4067 : !! pade
4068 : !!
4069 : !! FUNCTION
4070 : !! Calculate the pade approximant in zz of the function f calculated at the n points z
4071 : !!
4072 : !! SOURCE
4073 :
4074 965 : function pade(n, z, f, zz)
4075 :
4076 : !Arguments ------------------------------------
4077 : !scalars
4078 : integer,intent(in) :: n
4079 : complex(dp),intent(in) :: zz
4080 : complex(dp) :: pade
4081 : !arrays
4082 : complex(dp),intent(in) :: z(n), f(n)
4083 :
4084 : !Local variables-------------------------------
4085 : !scalars
4086 : integer :: i
4087 1930 : complex(dp) :: a(n), Az(0:n), Bz(0:n)
4088 : ! *************************************************************************
4089 :
4090 965 : call calculate_pade_a(a, n, z, f)
4091 :
4092 965 : Az(0)=czero
4093 965 : Az(1)=a(1)
4094 965 : Bz(0)=cone
4095 965 : Bz(1)=cone
4096 :
4097 9970 : do i=1,n-1
4098 9005 : Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4099 9970 : Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4100 : end do
4101 965 : pade=Az(n)/Bz(n)
4102 :
4103 : !write(std_out,*) 'pade_approx ', pade
4104 : !write(std_out,*) 'Bz(n)',Bz(n)
4105 : !if (real(Bz(n))==zero .and. aimag(Bz(n))==zero) write(std_out,*) ' Bz(n) ',Bz(n)
4106 :
4107 965 : end function pade
4108 : !!***
4109 :
4110 : !----------------------------------------------------------------------
4111 :
4112 : !!****f* m_numeric_tools/dpade
4113 : !! NAME
4114 : !! dpade
4115 : !!
4116 : !! FUNCTION
4117 : !! Calculate the derivative of the pade approximant in zz of the function f calculated at the n points z
4118 : !!
4119 : !! SOURCE
4120 :
4121 215 : function dpade(n, z, f, zz)
4122 :
4123 : !Arguments ------------------------------------
4124 : !scalars
4125 : integer,intent(in) :: n
4126 : complex(dp),intent(in) :: zz
4127 : complex(dp) :: dpade
4128 : !arrays
4129 : complex(dp),intent(in) :: z(n),f(n)
4130 :
4131 : !Local variables-------------------------------
4132 : !scalars
4133 : integer :: i
4134 : !arrays
4135 430 : complex(dp) :: a(n), Az(0:n), Bz(0:n), dAz(0:n), dBz(0:n)
4136 : ! *************************************************************************
4137 :
4138 215 : call calculate_pade_a(a, n, z, f)
4139 :
4140 215 : Az(0)=czero
4141 215 : Az(1)=a(1)
4142 215 : Bz(0)=cone
4143 215 : Bz(1)=cone
4144 215 : dAz(0)=czero
4145 215 : dAz(1)=czero
4146 215 : dBz(0)=czero
4147 215 : dBz(1)=czero
4148 :
4149 2470 : do i=1,n-1
4150 2255 : Az(i+1)=Az(i)+(zz-z(i))*a(i+1)*Az(i-1)
4151 2255 : Bz(i+1)=Bz(i)+(zz-z(i))*a(i+1)*Bz(i-1)
4152 2255 : dAz(i+1)=dAz(i)+a(i+1)*Az(i-1)+(zz-z(i))*a(i+1)*dAz(i-1)
4153 2470 : dBz(i+1)=dBz(i)+a(i+1)*Bz(i-1)+(zz-z(i))*a(i+1)*dBz(i-1)
4154 : end do
4155 : !write(std_out,*) 'Bz(n)', Bz(n)
4156 : !if (REAL(Bz(n))==zero.and.AIMAG(Bz(n))==zero) write(std_out,*) 'Bz(n)',Bz(n)
4157 : !pade_approx = Az(n) / Bz(n)
4158 215 : dpade=dAz(n)/Bz(n) -Az(n)*dBz(n)/(Bz(n)*Bz(n))
4159 : !write(std_out,*) 'pade_approx ', pade_approx
4160 :
4161 215 : end function dpade
4162 : !!***
4163 :
4164 : !----------------------------------------------------------------------
4165 :
4166 : !!****f* m_numeric_tools/calculate_pade_a
4167 : !! NAME
4168 : !! calculate_pade_a
4169 : !!
4170 : !! FUNCTION
4171 : !!
4172 : !! INPUTS
4173 : !!
4174 : !! OUTPUT
4175 : !!
4176 : !! SOURCE
4177 :
4178 1180 : subroutine calculate_pade_a(a, n, z, f)
4179 :
4180 : !Arguments ------------------------------------
4181 : !scalars
4182 : integer,intent(in) :: n
4183 : complex(dp),intent(in) :: z(n),f(n)
4184 : complex(dp),intent(out) :: a(n)
4185 :
4186 : !Local variables-------------------------------
4187 : !scalars
4188 : integer :: i,j
4189 : !arrays
4190 1180 : complex(dp) :: g(n,n)
4191 : ! *************************************************************************
4192 :
4193 13620 : g(1,1:n)=f(1:n)
4194 :
4195 12440 : do i=2,n
4196 74820 : do j=i,n
4197 : ! 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)
4198 73640 : g(i,j)=(g(i-1,i-1)-g(i-1,j)) / ((z(j)-z(i-1))*g(i-1,j))
4199 : !write(std_out,*) 'g_i(z_j)',i,j,g(i,j)
4200 : end do
4201 : end do
4202 13620 : do i=1,n
4203 13620 : a(i)=g(i,i)
4204 : end do
4205 : !write(std_out,*) 'a ',a(:)
4206 :
4207 1180 : end subroutine calculate_pade_a
4208 : !!***
4209 :
4210 : !----------------------------------------------------------------------
4211 :
4212 : !!****f* m_numeric_tools/newrap_step
4213 : !! NAME
4214 : !! newrap_step
4215 : !!
4216 : !! FUNCTION
4217 : !! Apply single step newton-raphson method to find the root of a complex function
4218 : !! z_k+1 = z_k - f(z_k) / (df/dz(z_k))
4219 : !!
4220 : !! SOURCE
4221 :
4222 109 : complex(dp) function newrap_step(z, f, df)
4223 :
4224 : !Arguments ------------------------------------
4225 : !scalars
4226 : complex(dp),intent(in) :: z,f,df
4227 :
4228 : !Local variables-------------------------------
4229 : real(dp) :: dfm2
4230 : ! *************************************************************************
4231 :
4232 109 : dfm2=ABS(df)*ABS(df)
4233 :
4234 109 : newrap_step = z - (f*CONJG(df))/dfm2
4235 : !& z-one/(ABS(df)*ABS(df)) * CMPLX( REAL(f)*REAL(df)+AIMAG(f)*AIMAG(df), -REAL(f)*AIMAG(df)+AIMAG(f)*EAL(df) )
4236 :
4237 109 : end function newrap_step
4238 : !!***
4239 :
4240 : !----------------------------------------------------------------------
4241 :
4242 : !!****f* m_numeric_tools/cross_product_int
4243 : !! NAME
4244 : !! cross_product_int
4245 : !!
4246 : !! FUNCTION
4247 : !! Return the cross product of two vectors with integer components.
4248 : !!
4249 3641 : pure function cross_product_int(vec1,vec2) result(res)
4250 :
4251 : !Arguments ------------------------------------
4252 : integer,intent(in) :: vec1(3),vec2(3)
4253 : integer :: res(3)
4254 : ! *************************************************************************
4255 :
4256 3641 : res(1)=vec1(2)*vec2(3)-vec1(3)*vec2(2)
4257 3641 : res(2)=vec1(3)*vec2(1)-vec1(1)*vec2(3)
4258 3641 : res(3)=vec1(1)*vec2(2)-vec1(2)*vec2(1)
4259 :
4260 3641 : end function cross_product_int
4261 : !!***
4262 :
4263 : !----------------------------------------------------------------------
4264 :
4265 : !!****f* m_numeric_tools/cross_product_rdp
4266 : !! NAME
4267 : !! cross_product_rdp
4268 : !!
4269 : !! FUNCTION
4270 : !! Return the cross product of two vectors with real double precision components.
4271 : !!
4272 1 : pure function cross_product_rdp(vec1,vec2) result(res)
4273 :
4274 : !Arguments ------------------------------------
4275 : real(dp),intent(in) :: vec1(3),vec2(3)
4276 : real(dp) :: res(3)
4277 : ! *************************************************************************
4278 :
4279 1 : res(1)=vec1(2)*vec2(3)-vec1(3)*vec2(2)
4280 1 : res(2)=vec1(3)*vec2(1)-vec1(1)*vec2(3)
4281 1 : res(3)=vec1(1)*vec2(2)-vec1(2)*vec2(1)
4282 :
4283 1 : end function cross_product_rdp
4284 : !!***
4285 :
4286 : !----------------------------------------------------------------------
4287 :
4288 : !!****f* m_numeric_tools/l2norm_rdp
4289 : !! NAME
4290 : !! l2norm_rdp
4291 : !!
4292 : !! FUNCTION
4293 : !! Return the length (ordinary L2 norm) of a vector.
4294 : !!
4295 :
4296 1 : pure function l2norm_rdp(vec) result(res)
4297 :
4298 : !Arguments ------------------------------------
4299 : real(dp),intent(in) :: vec(:)
4300 : real(dp) :: res
4301 : ! *************************************************************************
4302 :
4303 4 : res=SQRT(DOT_PRODUCT(vec,vec))
4304 :
4305 1 : end function l2norm_rdp
4306 : !!***
4307 :
4308 : !----------------------------------------------------------------------
4309 :
4310 : !!****f* m_numeric_tools/remove_copies
4311 : !! NAME
4312 : !! remove_copies
4313 : !!
4314 : !! FUNCTION
4315 : !! Given an initial set of elements, set_in, return the subset of inequivalent items
4316 : !! packed in the first n_out positions of set_in. Use the logical function is_equal
4317 : !! to define whether two items are equivalent.
4318 : !!
4319 : !! INPUTS
4320 : !! n_in=Initial number of elements.
4321 : !! is_equal=logical function used to discern if two items are equal.
4322 : !!
4323 : !! OUTPUT
4324 : !! n_out=Number of inequivalent items found.
4325 : !!
4326 : !! SIDE EFFECTS
4327 : !! set_in(3,n_in)=
4328 : !! In input the initial set of n_in elements
4329 : !! In output set_in(3,1:n_out) contains the inequivalent elements found.
4330 : !!
4331 : !! NOTES
4332 : !! The routines only deals with arrays of 3D-vectors although generalizing the
4333 : !! algorithm to nD-space is straightforward.
4334 : !!
4335 : !! SOURCE
4336 :
4337 3 : subroutine remove_copies(n_in, set_in, n_out, is_equal)
4338 :
4339 : !Arguments ------------------------------------
4340 : !scalars
4341 : integer,intent(in) :: n_in
4342 : integer,intent(out) :: n_out
4343 : !arrays
4344 : real(dp),target,intent(inout) :: set_in(3,n_in)
4345 :
4346 : interface
4347 : function is_equal(k1,k2)
4348 : use defs_basis
4349 : real(dp),intent(in) :: k1(3),k2(3)
4350 : logical :: is_equal
4351 : end function is_equal
4352 : end interface
4353 :
4354 : !Local variables-------------------------------
4355 : !scalars
4356 : integer :: ii,jj
4357 : logical :: isnew
4358 : !arrays
4359 : type rdp1d_pt
4360 : integer :: idx
4361 : real(dp),pointer :: rpt(:)
4362 : end type rdp1d_pt
4363 3 : type(rdp1d_pt),allocatable :: Ap(:)
4364 : ! *************************************************************************
4365 :
4366 9 : ABI_MALLOC(Ap,(n_in))
4367 3 : Ap(1)%idx = 1
4368 3 : Ap(1)%rpt => set_in(:,1)
4369 :
4370 3 : n_out=1
4371 21 : do ii=2,n_in
4372 :
4373 18 : isnew=.TRUE.
4374 75 : do jj=1,n_out
4375 75 : if (is_equal(set_in(:,ii),Ap(jj)%rpt(:))) then
4376 3 : isnew=.FALSE.
4377 : exit
4378 : end if
4379 : end do
4380 :
4381 3 : if (isnew) then
4382 15 : n_out=n_out+1
4383 15 : Ap(n_out)%rpt => set_in(:,ii)
4384 15 : Ap(n_out)%idx = ii
4385 : end if
4386 : end do
4387 :
4388 : ! The n_out inequivalent items are packed first.
4389 3 : if (n_out/=n_in) then
4390 7 : do ii=1,n_out
4391 6 : jj=Ap(ii)%idx
4392 25 : set_in(:,ii) = set_in(:,jj)
4393 : !write(std_out,*) Ap(ii)%idx,Ap(ii)%rpt(:)
4394 : end do
4395 : end if
4396 :
4397 3 : ABI_FREE(Ap)
4398 :
4399 3 : end subroutine remove_copies
4400 : !!***
4401 :
4402 : !----------------------------------------------------------------------
4403 :
4404 : !!****f* m_numeric_tools/denominator
4405 : !! NAME
4406 : !! denominator
4407 : !!
4408 : !! FUNCTION
4409 : !! Return the denominator of the rational number dd, sign is not considered.
4410 : !!
4411 : !! INPUTS
4412 : !! dd=The rational number
4413 : !! tolerance=Absolute tolerance
4414 : !!
4415 : !! OUTPUT
4416 : !! ierr=If /=0 the input number is not rational within the given tolerance.
4417 : !!
4418 : !! SOURCE
4419 :
4420 71652 : integer function denominator(dd,ierr,tolerance)
4421 :
4422 : !Arguments ------------------------------------
4423 : !scalars
4424 : integer,intent(out) :: ierr
4425 : real(dp),intent(in) :: dd
4426 : real(dp),optional,intent(in) :: tolerance
4427 :
4428 : !Local variables ------------------------------
4429 : !scalars
4430 : integer,parameter :: largest_integer = HUGE(1)
4431 : integer :: ii
4432 : real(dp) :: my_tol
4433 : !************************************************************************
4434 :
4435 71652 : ii=1
4436 71652 : my_tol=0.0001 ; if (PRESENT(tolerance)) my_tol=ABS(tolerance)
4437 84220 : do
4438 155872 : if (ABS(dd*ii-NINT(dd*ii))<my_tol) then
4439 71652 : denominator=ii
4440 71652 : ierr=0
4441 71652 : RETURN
4442 : end if
4443 : ! Handle the case in which dd is not rational within my_tol.
4444 84220 : if (ii==largest_integer) then
4445 0 : denominator=ii
4446 0 : ierr=-1
4447 0 : RETURN
4448 : end if
4449 84220 : ii=ii+1
4450 : end do
4451 :
4452 : end function denominator
4453 : !!***
4454 :
4455 : !----------------------------------------------------------------------
4456 :
4457 : !!****f* m_numeric_tools/mincm
4458 : !! NAME
4459 : !! mincm
4460 : !!
4461 : !! FUNCTION
4462 : !! Return the minimum common multiple of ii and jj.
4463 : !!
4464 : !! SOURCE
4465 :
4466 1255 : integer function mincm(ii,jj)
4467 :
4468 : !Arguments ------------------------------------
4469 : !scalars
4470 : integer,intent(in) :: ii,jj
4471 : !************************************************************************
4472 :
4473 1255 : if (ii==0.or.jj==0) then
4474 0 : ABI_BUG('ii==0 or jj==0')
4475 : end if
4476 :
4477 1255 : mincm=MAX(ii,jj)
4478 0 : do
4479 1255 : if ( ((mincm/ii)*ii)==mincm .and. ((mincm/jj)*jj)==mincm ) RETURN
4480 0 : mincm=mincm+1
4481 : end do
4482 :
4483 : end function mincm
4484 : !!***
4485 :
4486 : !----------------------------------------------------------------------
4487 :
4488 : !!****f* m_numeric_tools/continued_fract
4489 : !! NAME
4490 : !! continued_fract
4491 : !!
4492 : !! FUNCTION
4493 : !! This routine calculates the continued fraction:
4494 : !!
4495 : !! 1
4496 : !! f(z) = _______________________________
4497 : !! z - a1 - b1^2
4498 : !! _____________________
4499 : !! z - a2 - b2^2
4500 : !! ___________
4501 : !! z -a3 - ........
4502 : !!
4503 : !! INPUTS
4504 : !! nlev=Number of "levels" in the continued fraction.
4505 : !! term_type=Type of the terminator.
4506 : !! 0 --> No terminator.
4507 : !! -1 --> Assume constant coefficients for a_i and b_i for i>nlev with a_inf = a(nlev) and b_inf = b(nleb)
4508 : !! 1 --> Same as above but a_inf and b_inf are obtained by averaging over the nlev values.
4509 : !! aa(nlev)=Set of a_i coefficients.
4510 : !! bb(nlev)=Set of b_i coefficients.
4511 : !! nz=Number of points on the z-mesh.
4512 : !! zpts(nz)=z-mesh.
4513 : !!
4514 : !! OUTPUT
4515 : !! spectrum(nz)=Contains f(z) on the input mesh.
4516 : !!
4517 : !! SOURCE
4518 :
4519 9661 : subroutine continued_fract(nlev,term_type,aa,bb,nz,zpts,spectrum)
4520 :
4521 : !Arguments ------------------------------------
4522 : !scalars
4523 : integer,intent(in) :: nlev,term_type,nz
4524 : !arrays
4525 : real(dp),intent(in) :: bb(nlev)
4526 : complex(dp),intent(in) :: aa(nlev)
4527 : complex(dp),intent(in) :: zpts(nz)
4528 : complex(dp),intent(out) :: spectrum(nz)
4529 :
4530 : !Local variables ------------------------------
4531 : !scalars
4532 : integer :: it
4533 : real(dp) :: bb_inf,bg,bup,swap
4534 : complex(dp) :: aa_inf
4535 : character(len=500) :: msg
4536 : !arrays
4537 9661 : complex(dp),allocatable :: div(:),den(:)
4538 : !************************************************************************
4539 :
4540 28983 : ABI_MALLOC(div,(nz))
4541 19322 : ABI_MALLOC(den,(nz))
4542 :
4543 9661 : select case (term_type)
4544 : case (0) ! No terminator.
4545 805059 : div=czero
4546 : case (-1,1)
4547 3370 : if (term_type==-1) then
4548 0 : bb_inf=bb(nlev)
4549 0 : aa_inf=aa(nlev)
4550 : else
4551 89850 : bb_inf=SUM(bb)/nlev
4552 89850 : aa_inf=SUM(aa)/nlev
4553 : end if
4554 : ! Be careful with the sign of the SQRT.
4555 2032110 : div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
4556 : case (2)
4557 0 : ABI_ERROR("To be tested")
4558 0 : div = zero
4559 0 : if (nlev>4) then
4560 0 : bg=zero; bup=zero
4561 0 : do it=1,nlev,2
4562 0 : if (it+2<nlev) bg = bg + bb(it+2)
4563 0 : bup = bup + bb(it)
4564 : end do
4565 0 : bg = bg/(nlev/2+MOD(nlev,2))
4566 0 : bup = bg/((nlev+1)/2)
4567 : !if (iseven(nlev)) then
4568 0 : if (.not.iseven(nlev)) then
4569 0 : swap = bg
4570 0 : bg = bup
4571 0 : bup = bg
4572 : end if
4573 : !write(std_out,*)nlev,bg,bup
4574 : !Here be careful with the sign of SQRT
4575 0 : do it=1,nz
4576 : div(it) = half/zpts(it) * (bb(nlev)/bup)**2 * &
4577 0 : ( (zpts(it)**2 +bup**2 -bg**2) - SQRT( (zpts(it)**2+bup**2-bg**2)**2 -four*(zpts(it)*bup)**2) )
4578 : end do
4579 : end if
4580 :
4581 : case default
4582 0 : write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
4583 9661 : ABI_ERROR(msg)
4584 : end select
4585 :
4586 609034 : do it=nlev,2,-1
4587 112913309 : den(:) = zpts(:) - aa(it) - div(:)
4588 112922970 : div(:) = (bb(it-1)**2 )/ den(:)
4589 : end do
4590 :
4591 2846830 : den = zpts(:) - aa(1) - div(:)
4592 2846830 : div = one/den(:)
4593 :
4594 2837169 : spectrum = div
4595 9661 : ABI_FREE(div)
4596 9661 : ABI_FREE(den)
4597 :
4598 9661 : end subroutine continued_fract
4599 : !!***
4600 :
4601 : !----------------------------------------------------------------------
4602 :
4603 : !!****f* m_numeric_tools/cmplx_sphcart
4604 : !! NAME
4605 : !! cmplx_sphcart
4606 : !!
4607 : !! FUNCTION
4608 : !! Convert an array of complex values stored in spherical coordinates
4609 : !! to Cartesian coordinates with real and imaginary part or vice versa.
4610 : !!
4611 : !! INPUTS
4612 : !! from=Option specifying the format used to store the complex values. See below.
4613 : !! [units]=Option to specify if angles are given in "Radians" (default) or "Degrees".
4614 : !!
4615 : !! SIDE EFFECTS
4616 : !! carr(:,:):
4617 : !! input: array with complex values in Cartesian form if from="C" or spherical form if from="S"
4618 : !! output: array with values converted to the new representation.
4619 : !!
4620 : !! SOURCE
4621 :
4622 0 : subroutine cmplx_sphcart(carr, from, units)
4623 :
4624 : !Arguments ------------------------------------
4625 : !scalars
4626 : character(len=*),intent(in) :: from
4627 : character(len=*),optional,intent(in) :: units
4628 : !arrays
4629 : complex(dp),intent(inout) :: carr(:,:)
4630 :
4631 : !Local variables-------------------------------
4632 : !scalars
4633 : integer :: jj,ii
4634 : real(dp) :: rho,theta,fact
4635 : character(len=500) :: msg
4636 : ! *************************************************************************
4637 :
4638 0 : select case (from(1:1))
4639 :
4640 : case ("S","s") ! Spherical --> Cartesian
4641 :
4642 0 : fact = one
4643 0 : if (PRESENT(units)) then
4644 0 : if (units(1:1) == "D" .or. units(1:1) == "d") fact = two_pi/360_dp
4645 : end if
4646 :
4647 0 : do jj=1,SIZE(carr,DIM=2)
4648 0 : do ii=1,SIZE(carr,DIM=1)
4649 0 : rho = DBLE(carr(ii,jj))
4650 0 : theta= AIMAG(carr(ii,jj)) * fact
4651 0 : carr(ii,jj) = CMPLX(rho*DCOS(theta), rho*DSIN(theta), kind=dp)
4652 : end do
4653 : end do
4654 :
4655 : case ("C","c") ! Cartesian --> Spherical \theta = 2 arctan(y/(rho+x))
4656 :
4657 0 : fact = one
4658 0 : if (PRESENT(units)) then
4659 0 : if (units(1:1) == "D" .or. units(1:1) == "d") fact = 360_dp/two_pi
4660 : end if
4661 :
4662 0 : do jj=1,SIZE(carr,DIM=2)
4663 0 : do ii=1,SIZE(carr,DIM=1)
4664 0 : rho = SQRT(ABS(carr(ii,jj)))
4665 0 : if (rho > tol16) then
4666 0 : theta= two * ATAN( AIMAG(carr(ii,jj)) / (DBLE(carr(ii,jj)) + rho) )
4667 : else
4668 : theta= zero
4669 : end if
4670 0 : carr(ii,jj) = CMPLX(rho, theta*fact, kind=dp)
4671 : end do
4672 : end do
4673 :
4674 : case default
4675 0 : msg = " Wrong value for from: "//TRIM(from)
4676 0 : ABI_BUG(msg)
4677 : end select
4678 :
4679 0 : end subroutine cmplx_sphcart
4680 : !!***
4681 :
4682 : !----------------------------------------------------------------------
4683 :
4684 : !!****f* m_numeric_tools/pfactorize
4685 : !! NAME
4686 : !! pfactorize
4687 : !!
4688 : !! FUNCTION
4689 : !! Factorize a number in terms of an user-specified set of prime factors
4690 : !! nn = alpha * Prod_i p^i 1)
4691 : !!
4692 : !! INPUTS
4693 : !! nn=The number to be factorized.
4694 : !! nfactors=The number of factors
4695 : !! pfactors(nfactors)=The list of prime number e.g. (/ 2, 3, 5, 7, 11 /)
4696 : !!
4697 : !! OUTPUT
4698 : !! powers(nfactors+1)=
4699 : !! The first nfactors entries are the powers i in Eq.1. powers(nfactors+1) is alpha.
4700 : !!
4701 : !! SOURCE
4702 :
4703 0 : subroutine pfactorize(nn,nfactors,pfactors,powers)
4704 :
4705 : !Arguments ------------------------------------
4706 : !scalars
4707 : integer,intent(in) :: nn,nfactors
4708 : integer,intent(in) :: pfactors(nfactors)
4709 : integer,intent(out) :: powers (nfactors+1)
4710 :
4711 : !Local variables ------------------------------
4712 : !scalars
4713 : integer :: tnn,ifc,fact,ipow,maxpwr
4714 : ! *************************************************************************
4715 :
4716 0 : powers=0; tnn=nn
4717 :
4718 0 : fact_loop: do ifc=1,nfactors
4719 0 : fact = pfactors (ifc)
4720 0 : maxpwr = NINT ( LOG(DBLE(tnn))/LOG(DBLE(fact) ) ) + 1
4721 0 : do ipow=1,maxpwr
4722 0 : if (tnn==1) EXIT fact_loop
4723 0 : if ( MOD(tnn,fact)==0 ) then
4724 0 : tnn=tnn/fact
4725 0 : powers(ifc)=powers(ifc) + 1
4726 : end if
4727 : end do
4728 : end do fact_loop
4729 :
4730 0 : if ( nn /= tnn * PRODUCT( pfactors**powers(1:nfactors)) ) then
4731 0 : ABI_BUG('nn/=tnn!')
4732 : end if
4733 :
4734 0 : powers(nfactors+1) = tnn
4735 :
4736 0 : end subroutine pfactorize
4737 : !!***
4738 :
4739 : !----------------------------------------------------------------------
4740 :
4741 : !!****f* m_numeric_tools/isordered
4742 : !! NAME
4743 : !! isordered
4744 : !!
4745 : !! FUNCTION
4746 : !! Return .TRUE. if values in array arr are ordered.
4747 : !! Consider that two double precision numbers within tolerance tol are equal.
4748 : !!
4749 : !! INPUTS
4750 : !! nn=Size of arr.
4751 : !! arr(nn)=The array with real values to be tested.
4752 : !! direction= ">" for ascending numerical order.
4753 : !! ">" for decreasing numerical order.
4754 : !!
4755 : !! SOURCE
4756 :
4757 64 : function isordered_rdp(nn,arr,direction,tol) result(isord)
4758 :
4759 : !Arguments ------------------------------------
4760 : !scalars
4761 : integer,intent(in) :: nn
4762 : real(dp),intent(in) :: tol
4763 : logical :: isord
4764 : character(len=*),intent(in) :: direction
4765 : !arrays
4766 : real(dp),intent(in) :: arr(nn)
4767 :
4768 : !Local variables ------------------------------
4769 : !scalars
4770 : integer :: ii
4771 : real(dp) :: prev
4772 : character(len=500) :: msg
4773 : ! *************************************************************************
4774 :
4775 64 : prev = arr(1); isord =.TRUE.
4776 :
4777 : SELECT CASE (direction(1:1))
4778 : CASE(">")
4779 : ii=2;
4780 611 : do while (ii<=nn .and. isord)
4781 547 : if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) >= prev)
4782 547 : prev = arr(ii)
4783 611 : ii = ii +1
4784 : end do
4785 :
4786 : CASE("<")
4787 : ii=2;
4788 0 : do while (ii<=nn .and. isord)
4789 0 : if (ABS(arr(ii)-prev) > tol) isord = (arr(ii) <= prev)
4790 0 : prev = arr(ii)
4791 0 : ii = ii +1
4792 : end do
4793 :
4794 : CASE DEFAULT
4795 0 : msg = "Wrong direction: "//TRIM(direction)
4796 64 : ABI_ERROR(msg)
4797 : END SELECT
4798 :
4799 64 : end function isordered_rdp
4800 : !!***
4801 :
4802 : !----------------------------------------------------------------------
4803 :
4804 : !!****f* m_numeric_tools/stats_eval
4805 : !! NAME
4806 : !! stats_eval
4807 : !!
4808 : !! FUNCTION
4809 : !! Helper function used to calculate the statistical parameters of a dataset.
4810 : !!
4811 : !! INPUT
4812 : !! arr(:)=Array with the values.
4813 : !!
4814 : !! OUTPUT
4815 : !! stats<stats_t>=Data type storing the parameters of the data set.
4816 : !!
4817 : !! SOURCE
4818 :
4819 1507 : pure function stats_eval(arr) result(stats)
4820 :
4821 : !Arguments ------------------------------------
4822 : !scalars
4823 : type(stats_t) :: stats
4824 : !arrays
4825 : real(dp),intent(in) :: arr(:)
4826 :
4827 : !Local variables ------------------------------
4828 : !scalars
4829 : integer :: ii,nn
4830 : real(dp) :: xx,x2_sum
4831 : ! *************************************************************************
4832 :
4833 1507 : stats%min = +HUGE(one)
4834 1507 : stats%max = -HUGE(one)
4835 1507 : stats%mean = zero
4836 :
4837 1507 : nn = SIZE(arr)
4838 21341 : do ii=1,nn
4839 19834 : xx = arr(ii)
4840 19834 : stats%min = MIN(stats%min, xx)
4841 19834 : stats%max = MAX(stats%max, xx)
4842 21341 : stats%mean = stats%mean + xx
4843 : end do
4844 :
4845 1507 : stats%mean = stats%mean/nn
4846 :
4847 : ! Two-pass algorithm for the variance (more stable than the single-pass one).
4848 1507 : x2_sum = zero
4849 21341 : do ii=1,nn
4850 19834 : xx = arr(ii)
4851 21341 : x2_sum = x2_sum + (xx - stats%mean)*(xx - stats%mean)
4852 : end do
4853 :
4854 1507 : if (nn > 1) then
4855 1507 : stats%stdev = x2_sum/(nn-1)
4856 1507 : stats%stdev = SQRT(ABS(stats%stdev))
4857 : else
4858 : stats%stdev = zero
4859 : end if
4860 :
4861 1507 : end function stats_eval
4862 : !!***
4863 :
4864 : !----------------------------------------------------------------------
4865 :
4866 : !!****f* m_numeric_tools/wrap2_zero_one
4867 : !! NAME
4868 : !! wrap2_zero_one
4869 : !!
4870 : !! FUNCTION
4871 : !! Transforms a real number (num) in its corresponding reduced number
4872 : !! (red) in the interval [0,1[ where 1 is not included (tol12)
4873 : !! num=red+shift
4874 : !!
4875 : !! INPUTS
4876 : !! num=real number
4877 : !!
4878 : !! OUTPUT
4879 : !! red=reduced number of num in the interval [0,1[ where 1 is not included
4880 : !! shift=num-red
4881 : !!
4882 : !! SOURCE
4883 :
4884 151515 : elemental subroutine wrap2_zero_one(num, red, shift)
4885 :
4886 : !Arguments ------------------------------------
4887 : !scalars
4888 : real(dp),intent(in) :: num
4889 : real(dp),intent(out) :: red,shift
4890 : ! *************************************************************************
4891 :
4892 151515 : if (num>zero) then
4893 74140 : red=mod((num+tol12),one)-tol12
4894 : else
4895 77375 : red=-mod(-(num-one+tol12),one)+one-tol12
4896 : end if
4897 151515 : if(abs(red)<tol12)red=0.0_dp
4898 151515 : shift=num-red
4899 :
4900 151515 : end subroutine wrap2_zero_one
4901 : !!***
4902 :
4903 : !----------------------------------------------------------------------
4904 :
4905 : !!****f* m_numeric_tools/wrap2_pmhalf
4906 : !! NAME
4907 : !! wrap2_pmhalf
4908 : !!
4909 : !! FUNCTION
4910 : !! Transforms a real number (num) in its corresponding reduced number
4911 : !! (red) in the interval ]-1/2,1/2] where -1/2 is not included (tol12)
4912 : !! num=red+shift
4913 : !!
4914 : !! INPUTS
4915 : !! num=real number
4916 : !!
4917 : !! OUTPUT
4918 : !! red=reduced number of num in the interval ]-1/2,1/2] where -1/2 is not included
4919 : !! shift=num-red
4920 : !!
4921 : !! SOURCE
4922 :
4923 85876311 : elemental subroutine wrap2_pmhalf(num,red,shift)
4924 :
4925 : !Arguments -------------------------------
4926 : !scalars
4927 : real(dp),intent(in) :: num
4928 : real(dp),intent(out) :: red,shift
4929 :
4930 : ! *********************************************************************
4931 :
4932 85876311 : if (num>zero) then
4933 42064124 : red=mod((num+half-tol12),one)-half+tol12
4934 : else
4935 43812187 : red=-mod(-(num-half-tol12),one)+half+tol12
4936 : end if
4937 85876311 : if(abs(red)<tol12)red=0.0d0
4938 85876311 : shift=num-red
4939 :
4940 85876311 : end subroutine wrap2_pmhalf
4941 : !!***
4942 :
4943 : !----------------------------------------------------------------------
4944 :
4945 : !!***
4946 : !!****f* m_numeric_tools/linear_interpolation
4947 : !! NAME
4948 : !! interpol1d
4949 : !!
4950 : !! FUNCTION
4951 : !! Perform linear interpolation of a set of points pts_o with values val_o
4952 : !! The points pts_o and pts_i are assumed to be ascending ordered arrays.
4953 : !!
4954 : !! INPUTS
4955 : !! npts_o=Number of points in the original array.
4956 : !! npts_i=Number of points in the interpolated array.
4957 : !! pts_o(npts_o)=Points in the original array.
4958 : !! pts_i(npts_i)=Points in the interpolated array.
4959 : !! val_o(npts_o)=Values at the points pts_o.
4960 : !!
4961 : !! OUTPUT
4962 : !! res(npts_i)=Interpolated values at the points pts_i.
4963 : !!
4964 : !! SOURCE
4965 0 : pure function interpol1d(npts_o,npts_i,pts_o,pts_i,val_o) result(res)
4966 : !Arguments ------------------------------------
4967 : !arrays
4968 : integer,intent(in) :: npts_o, npts_i
4969 : real(dp),intent(in) :: pts_o(npts_o), pts_i(npts_i)
4970 : real(dp),intent(in) :: val_o(npts_o)
4971 : real(dp) :: res(npts_i)
4972 : !Local variables-------------------------------
4973 : !scalars
4974 : integer :: ii, jj
4975 : real(dp) :: x1, x2, y1, y2, slope
4976 : ! *************************************************************************
4977 :
4978 0 : res = zero
4979 :
4980 0 : do ii = 1, npts_i
4981 0 : do jj = 1, npts_o + 1
4982 0 : if (jj == npts_o + 1 .or. pts_i(ii) < pts_o(jj)) exit
4983 : end do
4984 0 : if (jj == 1) then
4985 0 : res(ii) = val_o(1)
4986 0 : else if (jj == npts_o + 1) then
4987 0 : res(ii) = val_o(npts_o)
4988 : else
4989 0 : x1 = pts_o(jj - 1)
4990 0 : x2 = pts_o(jj)
4991 0 : y1 = val_o(jj - 1)
4992 0 : y2 = val_o(jj)
4993 :
4994 0 : if (x2 == x1) then
4995 0 : res(ii) = (y1 + y2)/2.0_dp
4996 : else
4997 0 : slope = (y2 - y1) / (x2 - x1)
4998 0 : res(ii) = y1 + slope * (pts_i(ii) - x1)
4999 : end if
5000 : end if
5001 : end do
5002 :
5003 0 : end function interpol1d
5004 : !!***
5005 :
5006 0 : pure function interpol1d_c(npts_o,npts_i,pts_o,pts_i,val_o) result(res)
5007 : !Arguments ------------------------------------
5008 : !arrays
5009 : integer,intent(in) :: npts_o, npts_i
5010 : real(dp),intent(in) :: pts_o(npts_o), pts_i(npts_i)
5011 : complex(dp),intent(in) :: val_o(npts_o)
5012 : complex(dp) :: res(npts_i)
5013 : !Local variables-------------------------------
5014 0 : real(dp) :: val_o_r(npts_o), val_o_i(npts_o), res_r(npts_i), res_i(npts_i)
5015 : ! *************************************************************************
5016 :
5017 : ! Split the complex values into real and imaginary parts.
5018 0 : val_o_r = REAL(val_o, kind=dp)
5019 0 : val_o_i = AIMAG(val_o)
5020 :
5021 : ! Interpolate the real part.
5022 0 : res_r = interpol1d(npts_o, npts_i, pts_o, pts_i, val_o_r)
5023 :
5024 : ! Interpolate the imaginary part.
5025 0 : res_i = interpol1d(npts_o, npts_i, pts_o, pts_i, val_o_i)
5026 :
5027 : ! Combine the results back into complex form.
5028 0 : res = CMPLX(res_r, res_i, kind=dp)
5029 0 : end function interpol1d_c
5030 :
5031 : !!----------------------------------------------------------------------
5032 : !!****f* m_numeric_tools/interpol3d_0d
5033 : !! NAME
5034 : !! interpol3d_0d
5035 : !!
5036 : !! FUNCTION
5037 : !! Computes the value at any point r by linear interpolation
5038 : !! inside the eight vertices of the surrounding cube
5039 : !! r is presumed to be normalized, in a unit cube for the full grid
5040 : !!
5041 : !! INPUTS
5042 : !! r(3)=point coordinate
5043 : !! nr1=grid size along x
5044 : !! nr2=grid size along y
5045 : !! nr3=grid size along z
5046 : !! grid(nr1,nr2,nr3)=grid matrix
5047 : !!
5048 : !! OUTPUT
5049 : !! res=Interpolated value
5050 : !!
5051 : !! SOURCE
5052 :
5053 2042 : pure function interpol3d_0d(rr, nr1, nr2, nr3, grid) result(res)
5054 :
5055 : !Arguments-------------------------------------------------------------
5056 : !scalars
5057 : integer,intent(in) :: nr1, nr2, nr3
5058 : real(dp) :: res
5059 : !arrays
5060 : real(dp),intent(in) :: grid(nr1,nr2,nr3),rr(3)
5061 :
5062 : !Local variables--------------------------------------------------------
5063 : integer :: ir1,ir2,ir3,pr1,pr2,pr3
5064 : real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8, x1,x2,x3
5065 : ! *************************************************************************
5066 :
5067 2042 : call interpol3d_indices(rr, nr1, nr2, nr3, ir1, ir2, ir3, pr1, pr2, pr3)
5068 :
5069 : ! weight
5070 2042 : x1=one+rr(1)*nr1-real(ir1)
5071 2042 : x2=one+rr(2)*nr2-real(ir2)
5072 2042 : x3=one+rr(3)*nr3-real(ir3)
5073 :
5074 : !calculation of the density value
5075 2042 : res1=grid(ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
5076 2042 : res2=grid(pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
5077 2042 : res3=grid(ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
5078 2042 : res4=grid(ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
5079 2042 : res5=grid(pr1, pr2, ir3) * x1*x2*(one-x3)
5080 2042 : res6=grid(ir1, pr2, pr3) * (one-x1)*x2*x3
5081 2042 : res7=grid(pr1, ir2, pr3) * x1*(one-x2)*x3
5082 2042 : res8=grid(pr1, pr2, pr3) * x1*x2*x3
5083 2042 : res=res1+res2+res3+res4+res5+res6+res7+res8
5084 :
5085 2042 : end function interpol3d_0d
5086 : !!***
5087 :
5088 : !----------------------------------------------------------------------
5089 :
5090 : !!****f* m_numeric_tools/interpol3d_1d
5091 : !! NAME
5092 : !! interpol3d_1d
5093 : !!
5094 : !! FUNCTION
5095 : !! Computes the value at any point r by linear interpolation
5096 : !! inside the eight vertices of the surrounding cube
5097 : !! r is presumed to be normalized, in a unit cube for the full grid
5098 : !!
5099 : !! INPUTS
5100 : !! r(3)=point coordinate
5101 : !! nr1=grid size along x
5102 : !! nr2=grid size along y
5103 : !! nr3=grid size along z
5104 : !! grid(cplex,nr1,nr2,nr3)=grid matrix
5105 : !!
5106 : !! OUTPUT
5107 : !! res(cplex)=Interpolated value
5108 : !!
5109 : !! SOURCE
5110 :
5111 223896 : pure function interpol3d_1d(rr, nr1, nr2, nr3, grid, cplex) result(res)
5112 :
5113 : !Arguments-------------------------------------------------------------
5114 : !scalars
5115 : integer,intent(in) :: nr1, nr2, nr3, cplex
5116 : real(dp) :: res(cplex)
5117 : !arrays
5118 : real(dp),intent(in) :: grid(cplex, nr1, nr2, nr3), rr(3)
5119 :
5120 : !Local variables--------------------------------------------------------
5121 : !scalars
5122 : integer :: id,ir1,ir2,ir3,pr1,pr2,pr3
5123 : real(dp) :: res1,res2,res3,res4,res5,res6,res7,res8,x1,x2,x3
5124 : ! *************************************************************************
5125 :
5126 223896 : call interpol3d_indices(rr, nr1, nr2, nr3, ir1, ir2, ir3, pr1, pr2, pr3)
5127 :
5128 : ! weight
5129 223896 : x1 = one + rr(1)*nr1 -real(ir1)
5130 223896 : x2 = one + rr(2)*nr2 -real(ir2)
5131 223896 : x3 = one + rr(3)*nr3 -real(ir3)
5132 :
5133 : ! calculation of the density value
5134 447792 : do id=1,cplex
5135 223896 : res1 = grid(id,ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
5136 223896 : res2 = grid(id,pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
5137 223896 : res3 = grid(id,ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
5138 223896 : res4 = grid(id,ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
5139 223896 : res5 = grid(id,pr1, pr2, ir3) * x1*x2*(one-x3)
5140 223896 : res6 = grid(id,ir1, pr2, pr3) * (one-x1)*x2*x3
5141 223896 : res7 = grid(id,pr1, ir2, pr3) * x1*(one-x2)*x3
5142 223896 : res8 = grid(id,pr1, pr2, pr3) * x1*x2*x3
5143 447792 : res(id) = res1+res2+res3+res4+res5+res6+res7+res8
5144 : end do
5145 :
5146 223896 : end function interpol3d_1d
5147 : !!***
5148 :
5149 : !----------------------------------------------------------------------
5150 :
5151 : !!****f* m_numeric_tools/interpol3d_indices
5152 : !! NAME
5153 : !! interpol3d_indices
5154 : !!
5155 : !! FUNCTION
5156 : !! Computes the indices in a cube which are neighbors to the point to be
5157 : !! interpolated in interpol3d
5158 : !!
5159 : !! INPUTS
5160 : !! rr(3)=point coordinate
5161 : !! nr1=grid size along x
5162 : !! nr2=grid size along y
5163 : !! nr3=grid size along z
5164 : !!
5165 : !! OUTPUT
5166 : !! ir1,ir2,ir3 = bottom left neighbor
5167 : !! pr1,pr2,pr3 = top right neighbor
5168 : !!
5169 : !! SOURCE
5170 :
5171 518537 : pure subroutine interpol3d_indices(rr,nr1,nr2,nr3,ir1,ir2,ir3,pr1,pr2,pr3)
5172 :
5173 : !Arguments-------------------------------------------------------------
5174 : !scalars
5175 : integer,intent(in) :: nr1, nr2, nr3
5176 : integer,intent(out) :: ir1, ir2, ir3, pr1, pr2, pr3
5177 : !arrays
5178 : real(dp),intent(in) :: rr(3)
5179 :
5180 : !Local variables-------------------------------
5181 : real(dp) :: d1,d2,d3
5182 : ! *************************************************************************
5183 :
5184 : !grid density
5185 518537 : d1=one/nr1
5186 518537 : d2=one/nr2
5187 518537 : d3=one/nr3
5188 :
5189 : !lower left
5190 518537 : ir1=int(rr(1)/d1)+1
5191 518537 : ir2=int(rr(2)/d2)+1
5192 518537 : ir3=int(rr(3)/d3)+1
5193 :
5194 : !upper right
5195 518537 : pr1=mod(ir1+1,nr1)
5196 518537 : pr2=mod(ir2+1,nr2)
5197 518537 : pr3=mod(ir3+1,nr3)
5198 :
5199 518537 : if(ir1==0) ir1=nr1
5200 518537 : if(ir2==0) ir2=nr2
5201 518537 : if(ir3==0) ir3=nr3
5202 :
5203 518537 : if(ir1>nr1) ir1=ir1-nr1
5204 518537 : if(ir2>nr2) ir2=ir2-nr2
5205 518537 : if(ir3>nr3) ir3=ir3-nr3
5206 :
5207 518537 : if(pr1==0) pr1=nr1
5208 518537 : if(pr2==0) pr2=nr2
5209 518537 : if(pr3==0) pr3=nr3
5210 :
5211 518537 : end subroutine interpol3d_indices
5212 : !!***
5213 :
5214 : !----------------------------------------------------------------------
5215 :
5216 : !!****f* m_numeric_tools/interpolate_denpot
5217 : !! NAME
5218 : !! interpolate_denpot
5219 : !!
5220 : !! FUNCTION
5221 : !! Linear interpolation of density/potential given on the real space FFT mesh.
5222 : !! Assumes array on full mesh i.e. no MPI-FFT.
5223 : !!
5224 : !! INPUTS
5225 : !! cplex=1 for real, 2 for complex data.
5226 : !! in_ngfft(3)=Mesh divisions of input array
5227 : !! nspden=Number of density components.
5228 : !! in_rhor(cplex * in_nfftot * nspden)=Input array
5229 : !! out_ngfft(3)=Mesh divisions of output array
5230 : !!
5231 : !! OUTPUT
5232 : !! outrhor(cplex * out_nfftot * nspden)=Output array with interpolated data.
5233 : !!
5234 : !! SOURCE
5235 :
5236 119 : subroutine interpolate_denpot(cplex, in_ngfft, nspden, in_rhor, out_ngfft, out_rhor)
5237 :
5238 : !Arguments-------------------------------------------------------------
5239 : !scalars
5240 : integer,intent(in) :: cplex,nspden
5241 : !arrays
5242 : integer,intent(in) :: in_ngfft(3), out_ngfft(3)
5243 : real(dp),intent(in) :: in_rhor(cplex, product(in_ngfft), nspden)
5244 : real(dp),intent(out) :: out_rhor(cplex, product(out_ngfft), nspden)
5245 :
5246 : !Local variables--------------------------------------------------------
5247 : !scalars
5248 : integer :: ispden, ir1, ir2, ir3, ifft
5249 : real(dp) :: rr(3)
5250 : ! *************************************************************************
5251 :
5252 : ! Linear interpolation.
5253 34 : do ispden=1,nspden
5254 376 : do ir3=0,out_ngfft(3) - 1
5255 342 : rr(3) = dble(ir3) / out_ngfft(3)
5256 8275 : do ir2=0,out_ngfft(2) - 1
5257 7916 : rr(2) = dble(ir2) / out_ngfft(2)
5258 232154 : do ir1=0,out_ngfft(1) - 1
5259 223896 : rr(1) = dble(ir1) / out_ngfft(1)
5260 223896 : ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5261 231812 : out_rhor(1:cplex, ifft, ispden) = interpol3d_1d(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_rhor(:,:,ispden), cplex)
5262 : end do
5263 : end do
5264 : end do
5265 : end do
5266 :
5267 17 : end subroutine interpolate_denpot
5268 : !!***
5269 :
5270 :
5271 : !!****f* m_numeric_tools/interpolate_ur_spc
5272 : !! NAME
5273 : !! interpolate_ur_spc
5274 : !!
5275 : !! FUNCTION
5276 : !! Linear interpolation of complex wavefunctions given on the real space FFT mesh.
5277 : !! Assumes array on full mesh i.e. no MPI-FFT.
5278 : !! Single precision version
5279 : !!
5280 : !! INPUTS
5281 : !! in_ngfft(3)=Mesh divisions of input array
5282 : !! ndat=Number of wavefunctions.
5283 : !! in_ur(in_nfftot * ndat)=Input array
5284 : !! out_ngfft(3)=Mesh divisions of output array
5285 : !!
5286 : !! OUTPUT
5287 : !! out_ur(out_nfftot,ndata)=Output array with interpolated data.
5288 : !!
5289 : !! SOURCE
5290 :
5291 0 : subroutine interpolate_ur_spc(in_ngfft, ndat, in_ur, out_ngfft, out_ur)
5292 :
5293 : !Arguments-------------------------------------------------------------
5294 : !scalars
5295 : integer,intent(in) :: ndat
5296 : !arrays
5297 : integer,intent(in) :: in_ngfft(:), out_ngfft(:)
5298 : complex(sp),intent(in) :: in_ur(product(in_ngfft(1:3)), ndat)
5299 : complex(sp),intent(out) :: out_ur(product(out_ngfft(1:3)), ndat)
5300 :
5301 : !Local variables--------------------------------------------------------
5302 : !scalars
5303 : integer :: idat, ir1, ir2, ir3, ifft
5304 : real(dp) :: rr(3)
5305 : ! *************************************************************************
5306 :
5307 : ! Linear interpolation.
5308 0 : do idat=1,ndat
5309 0 : do ir3=0,out_ngfft(3) - 1
5310 0 : rr(3) = dble(ir3) / out_ngfft(3)
5311 0 : do ir2=0,out_ngfft(2) - 1
5312 0 : rr(2) = dble(ir2) / out_ngfft(2)
5313 0 : do ir1=0,out_ngfft(1) - 1
5314 0 : rr(1) = dble(ir1) / out_ngfft(1)
5315 0 : ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5316 0 : out_ur(ifft, idat) = interpol3d_1d_spc(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_ur(:,idat))
5317 : end do
5318 : end do
5319 : end do
5320 : end do
5321 :
5322 0 : end subroutine interpolate_ur_spc
5323 : !!***
5324 :
5325 : !----------------------------------------------------------------------
5326 :
5327 : !!****f* m_numeric_tools/interpol3d_1d_spc
5328 : !! NAME
5329 : !! interpol3d_1d_spc
5330 : !!
5331 : !! FUNCTION
5332 : !! Computes the value at any point r by linear interpolation
5333 : !! inside the eight vertices of the surrounding cube
5334 : !! r is presumed to be normalized, in a unit cube for the full grid
5335 : !!
5336 : !! INPUTS
5337 : !! r(3)=point coordinate
5338 : !! nr1=grid size along x
5339 : !! nr2=grid size along y
5340 : !! nr3=grid size along z
5341 : !! grid(nr1,nr2,nr3)=grid matrix
5342 : !!
5343 : !! OUTPUT
5344 : !! resInterpolated value
5345 : !!
5346 : !! SOURCE
5347 :
5348 0 : pure complex(sp) function interpol3d_1d_spc(rr, nr1, nr2, nr3, grid) result(res)
5349 :
5350 : !Arguments-------------------------------------------------------------
5351 : !scalars
5352 : integer,intent(in) :: nr1, nr2, nr3
5353 : !arrays
5354 : real(dp),intent(in) :: rr(3)
5355 : complex(sp),intent(in) :: grid(nr1, nr2, nr3)
5356 :
5357 : !Local variables--------------------------------------------------------
5358 : !scalars
5359 : integer :: ir1,ir2,ir3,pr1,pr2,pr3
5360 : complex(sp) :: res1,res2,res3,res4,res5,res6,res7,res8,x1,x2,x3
5361 : ! *************************************************************************
5362 :
5363 0 : call interpol3d_indices(rr, nr1, nr2, nr3, ir1, ir2, ir3, pr1, pr2, pr3)
5364 :
5365 : ! weight
5366 0 : x1 = one + rr(1)*nr1 -real(ir1)
5367 0 : x2 = one + rr(2)*nr2 -real(ir2)
5368 0 : x3 = one + rr(3)*nr3 -real(ir3)
5369 :
5370 : ! calculation of the density value
5371 0 : res1 = grid(ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
5372 0 : res2 = grid(pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
5373 0 : res3 = grid(ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
5374 0 : res4 = grid(ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
5375 0 : res5 = grid(pr1, pr2, ir3) * x1*x2*(one-x3)
5376 0 : res6 = grid(ir1, pr2, pr3) * (one-x1)*x2*x3
5377 0 : res7 = grid(pr1, ir2, pr3) * x1*(one-x2)*x3
5378 0 : res8 = grid(pr1, pr2, pr3) * x1*x2*x3
5379 0 : res = res1+res2+res3+res4+res5+res6+res7+res8
5380 :
5381 0 : end function interpol3d_1d_spc
5382 : !!***
5383 :
5384 : !!****f* m_numeric_tools/interpolate_ur_dpc
5385 : !! NAME
5386 : !! interpolate_ur_dpc
5387 : !!
5388 : !! FUNCTION
5389 : !! Linear interpolation of complex wavefunctions given on the real space FFT mesh.
5390 : !! Assumes array on full mesh i.e. no MPI-FFT.
5391 : !! Double precision version
5392 : !!
5393 : !! INPUTS
5394 : !! in_ngfft(3)=Mesh divisions of input array
5395 : !! ndat=Number of wavefunctions.
5396 : !! in_ur(in_nfftot * ndat)=Input array
5397 : !! out_ngfft(3)=Mesh divisions of output array
5398 : !!
5399 : !! OUTPUT
5400 : !! out_ur(out_nfftot,ndata)=Output array with interpolated data.
5401 : !!
5402 : !! SOURCE
5403 :
5404 1190 : subroutine interpolate_ur_dpc(in_ngfft, ndat, in_ur, out_ngfft, out_ur)
5405 :
5406 : !Arguments-------------------------------------------------------------
5407 : !scalars
5408 : integer,intent(in) :: ndat
5409 : !arrays
5410 : integer,intent(in) :: in_ngfft(:), out_ngfft(:)
5411 : complex(dp),intent(in) :: in_ur(product(in_ngfft(1:3)), ndat)
5412 : complex(dp),intent(out) :: out_ur(product(out_ngfft(1:3)), ndat)
5413 :
5414 : !Local variables--------------------------------------------------------
5415 : !scalars
5416 : integer :: idat, ir1, ir2, ir3, ifft
5417 : real(dp) :: rr(3)
5418 : ! *************************************************************************
5419 :
5420 : ! Linear interpolation.
5421 340 : do idat=1,ndat
5422 2089 : do ir3=0,out_ngfft(3) - 1
5423 1749 : rr(3) = dble(ir3) / out_ngfft(3)
5424 23348 : do ir2=0,out_ngfft(2) - 1
5425 21429 : rr(2) = dble(ir2) / out_ngfft(2)
5426 315777 : do ir1=0,out_ngfft(1) - 1
5427 292599 : rr(1) = dble(ir1) / out_ngfft(1)
5428 292599 : ifft = 1 + ir1 + ir2*out_ngfft(1) + ir3*out_ngfft(1)*out_ngfft(2)
5429 314028 : out_ur(ifft, idat) = interpol3d_1d_dpc(rr, in_ngfft(1), in_ngfft(2), in_ngfft(3), in_ur(:,idat))
5430 : end do
5431 : end do
5432 : end do
5433 : end do
5434 :
5435 170 : end subroutine interpolate_ur_dpc
5436 : !!***
5437 :
5438 : !----------------------------------------------------------------------
5439 :
5440 : !!****f* m_numeric_tools/interpol3d_1d_dpc
5441 : !! NAME
5442 : !! interpol3d_1d_dpc
5443 : !!
5444 : !! FUNCTION
5445 : !! Computes the value at any point r by linear interpolation
5446 : !! inside the eight vertices of the surrounding cube
5447 : !! r is presumed to be normalized, in a unit cube for the full grid
5448 : !!
5449 : !! INPUTS
5450 : !! r(3)=point coordinate
5451 : !! nr1=grid size along x
5452 : !! nr2=grid size along y
5453 : !! nr3=grid size along z
5454 : !! grid(nr1,nr2,nr3)=grid matrix
5455 : !!
5456 : !! OUTPUT
5457 : !! resInterpolated value
5458 : !!
5459 : !! SOURCE
5460 :
5461 292599 : pure complex(dp) function interpol3d_1d_dpc(rr, nr1, nr2, nr3, grid) result(res)
5462 :
5463 : !Arguments-------------------------------------------------------------
5464 : !scalars
5465 : integer,intent(in) :: nr1, nr2, nr3
5466 : !arrays
5467 : real(dp),intent(in) :: rr(3)
5468 : complex(dp),intent(in) :: grid(nr1, nr2, nr3)
5469 :
5470 : !Local variables--------------------------------------------------------
5471 : !scalars
5472 : integer :: ir1,ir2,ir3,pr1,pr2,pr3
5473 : complex(dp) :: res1,res2,res3,res4,res5,res6,res7,res8,x1,x2,x3
5474 : ! *************************************************************************
5475 :
5476 292599 : call interpol3d_indices(rr, nr1, nr2, nr3, ir1, ir2, ir3, pr1, pr2, pr3)
5477 :
5478 : ! weight
5479 292599 : x1 = one + rr(1)*nr1 -real(ir1)
5480 292599 : x2 = one + rr(2)*nr2 -real(ir2)
5481 292599 : x3 = one + rr(3)*nr3 -real(ir3)
5482 :
5483 : ! calculation of the density value
5484 292599 : res1 = grid(ir1, ir2, ir3) * (one-x1)*(one-x2)*(one-x3)
5485 292599 : res2 = grid(pr1, ir2, ir3) * x1*(one-x2)*(one-x3)
5486 292599 : res3 = grid(ir1, pr2, ir3) * (one-x1)*x2*(one-x3)
5487 292599 : res4 = grid(ir1, ir2, pr3) * (one-x1)*(one-x2)*x3
5488 292599 : res5 = grid(pr1, pr2, ir3) * x1*x2*(one-x3)
5489 292599 : res6 = grid(ir1, pr2, pr3) * (one-x1)*x2*x3
5490 292599 : res7 = grid(pr1, ir2, pr3) * x1*(one-x2)*x3
5491 292599 : res8 = grid(pr1, pr2, pr3) * x1*x2*x3
5492 292599 : res = res1+res2+res3+res4+res5+res6+res7+res8
5493 :
5494 292599 : end function interpol3d_1d_dpc
5495 : !!***
5496 :
5497 : !----------------------------------------------------------------------
5498 :
5499 : !!****f* m_numeric_tools/simpson_int
5500 : !! NAME
5501 : !! simpson_int
5502 : !!
5503 : !! FUNCTION
5504 : !! Simpson integral of input function
5505 : !!
5506 : !! INPUTS
5507 : !! npts=max number of points on grid for integral
5508 : !! step = space between integral arguments
5509 : !! values(npts)=integrand function.
5510 : !!
5511 : !! OUTPUT
5512 : !! int_values(npts)=integral of values.
5513 : !!
5514 : !! SOURCE
5515 :
5516 3221169 : subroutine simpson_int(npts, step, values, int_values)
5517 :
5518 : !Arguments ------------------------------------
5519 : !scalars
5520 : integer,intent(in) :: npts
5521 : real(dp),intent(in) :: step
5522 : !arrays
5523 : real(dp),intent(in) :: values(npts)
5524 : real(dp),intent(out) :: int_values(npts)
5525 :
5526 : !Local variables -------------------------
5527 : !scalars
5528 : integer :: ii
5529 : real(dp),parameter :: coef1 = 0.375_dp !9.0_dp / 24.0_dp
5530 : real(dp),parameter :: coef2 = 1.166666666666666666666666667_dp !28.0_dp / 24.0_dp
5531 : real(dp),parameter :: coef3 = 0.958333333333333333333333333_dp !23.0_dp / 24.0_dp
5532 : character(len=500) :: msg
5533 : ! *********************************************************************
5534 :
5535 3221169 : if (npts < 6) then
5536 0 : write(msg,"(a,i0)")"Number of points in integrand function must be >=6 while it is: ",npts
5537 0 : ABI_ERROR(msg)
5538 : end if
5539 :
5540 : !-----------------------------------------------------------------
5541 : !Simpson integral of input function
5542 : !-----------------------------------------------------------------
5543 :
5544 : !first point is 0: don t store it
5545 : !do integration equivalent to Simpson O(1/N^4) from NumRec in C p 134 NumRec in Fortran p 128
5546 3221169 : int_values(1) = coef1*values(1)
5547 3221169 : int_values(2) = int_values(1) + coef2*values(2)
5548 3221169 : int_values(3) = int_values(2) + coef3*values(3)
5549 :
5550 3294184923 : do ii=4,npts-3
5551 3294184923 : int_values(ii) = int_values(ii-1) + values(ii)
5552 : end do
5553 :
5554 3221169 : int_values(npts-2) = int_values(npts-3) + coef3*values(npts-2)
5555 3221169 : int_values(npts-1) = int_values(npts-2) + coef2*values(npts-1)
5556 3221169 : int_values(npts ) = int_values(npts-1) + coef1*values(npts )
5557 :
5558 3313511937 : int_values(:) = int_values(:) * step
5559 :
5560 3221169 : end subroutine simpson_int
5561 : !!***
5562 :
5563 : !----------------------------------------------------------------------
5564 :
5565 : !!****f* m_numeric_tools/simpson
5566 : !! NAME
5567 : !! simpson
5568 : !!
5569 : !! FUNCTION
5570 : !! Simpson integral of input function
5571 : !!
5572 : !! INPUTS
5573 : !! step = space between integral arguments
5574 : !! values(npts)=integrand function.
5575 : !!
5576 : !! OUTPUT
5577 : !! integral of values on the full mesh.
5578 : !!
5579 : !! SOURCE
5580 :
5581 796771 : function simpson(step, values) result(res)
5582 :
5583 : !Arguments ------------------------------------
5584 : !scalars
5585 : real(dp),intent(in) :: step
5586 : real(dp) :: res
5587 : !arrays
5588 : real(dp),intent(in) :: values(:)
5589 :
5590 : !Local variables -------------------------
5591 1593542 : real(dp) :: int_values(size(values))
5592 : ! *********************************************************************
5593 :
5594 796771 : call simpson_int(size(values),step,values,int_values)
5595 796771 : res = int_values(size(values))
5596 :
5597 796771 : end function simpson
5598 : !!***
5599 :
5600 : !----------------------------------------------------------------------
5601 :
5602 : !!****f* m_numeric_tools/rhophi
5603 : !! NAME
5604 : !! rhophi
5605 : !!
5606 : !! FUNCTION
5607 : !! Compute the phase and the module of a complex number.
5608 : !! The phase angle is fold into the interval [-pi,pi]
5609 : !!
5610 : !! INPUTS
5611 : !! cx(2) = complex number
5612 : !!
5613 : !! OUTPUT
5614 : !! phi = phase of cx fold into [-pi,pi]
5615 : !! rho = module of cx
5616 : !!
5617 : !! SOURCE
5618 :
5619 290874588 : pure subroutine rhophi(cx, phi, rho)
5620 :
5621 : !Arguments ------------------------------------
5622 : !scalars
5623 : real(dp),intent(out) :: phi,rho
5624 : !arrays
5625 : real(dp),intent(in) :: cx(2)
5626 : ! ***********************************************************************
5627 :
5628 290874588 : rho = sqrt(cx(1)*cx(1) + cx(2)*cx(2))
5629 :
5630 290874588 : if (abs(cx(1)) > tol8) then
5631 290874364 : phi = atan(cx(2)/cx(1))
5632 :
5633 : ! phi is an element of [-pi,pi]
5634 290874364 : if (cx(1) < zero) then
5635 143789994 : if (phi < zero) then
5636 73777284 : phi = phi + pi
5637 : else
5638 70012710 : phi = phi - pi
5639 : end if
5640 : end if
5641 :
5642 : else
5643 :
5644 224 : if (cx(2) > tol8) then
5645 43 : phi = pi*half
5646 181 : else if (cx(2) < tol8) then
5647 181 : phi = -pi*half
5648 : else
5649 0 : phi = zero
5650 : end if
5651 :
5652 : end if
5653 :
5654 290874588 : end subroutine rhophi
5655 : !!***
5656 :
5657 : !----------------------------------------------------------------------
5658 :
5659 : !!****f* m_numeric_tools/vdiff_eval
5660 : !! NAME
5661 : !! vdiff_eval
5662 : !!
5663 : !! FUNCTION
5664 : !! Estimate the "distance" between two functions tabulated on a homogeneous grid.
5665 : !! See vdiff_t
5666 : !!
5667 : !! INPUTS
5668 : !! cplex=1 if f1 and f2 are real, 2 for complex.
5669 : !! nr=Number of points in the mesh.
5670 : !! f1(cplex,nr), f2(cplex,nr)=Vectors with values
5671 : !! [vd_max]= Compute max value of the different entries.
5672 : !!
5673 : !! OUTPUT
5674 : !! vdiff_t object
5675 : !!
5676 : !! SOURCE
5677 :
5678 0 : subroutine vdiff_eval(vd, cplex, nr, f1, f2, volume, vd_max, unit)
5679 :
5680 : !Arguments ------------------------------------
5681 : !scalars
5682 : class(vdiff_t),intent(out) :: vd
5683 : integer,intent(in) :: cplex,nr
5684 : real(dp),intent(in) :: volume
5685 : type(vdiff_t),optional,intent(inout) :: vd_max
5686 : integer,optional,intent(in) :: unit
5687 : !arrays
5688 : real(dp),intent(in) :: f1(cplex,nr),f2(cplex,nr)
5689 :
5690 : !Local variables-------------------------------
5691 : !scalars
5692 : integer :: ir
5693 : real(dp) :: num,den,dr
5694 : type(stats_t) :: stats
5695 : !arrays
5696 0 : real(dp) :: abs_diff(nr)
5697 : ! *********************************************************************
5698 :
5699 0 : dr = volume / nr
5700 :
5701 0 : if (cplex == 1) then
5702 0 : abs_diff = abs(f1(1,:) - f2(1,:))
5703 0 : num = sum(abs_diff)
5704 0 : den = sum(abs(f2(1,:)))
5705 :
5706 0 : else if (cplex == 2) then
5707 0 : do ir=1,nr
5708 0 : abs_diff(ir) = sqrt((f1(1,ir) - f2(1,ir))**2 + (f1(2,ir) - f2(2,ir))**2)
5709 : end do
5710 0 : num = sum(abs_diff)
5711 0 : den = zero
5712 0 : do ir=1,nr
5713 0 : den = den + sqrt(f2(1,ir)**2 + f2(2,ir)**2)
5714 : end do
5715 : end if
5716 :
5717 0 : vd%int_adiff = num * dr
5718 0 : call safe_div(num,den,zero,vd%l1_rerr)
5719 :
5720 0 : stats = stats_eval(abs_diff)
5721 0 : vd%mean_adiff = stats%mean
5722 0 : vd%stdev_adiff = stats%stdev
5723 0 : vd%min_adiff = stats%min
5724 0 : vd%max_adiff = stats%max
5725 :
5726 0 : if (present(vd_max)) then
5727 0 : vd_max%int_adiff = max(vd_max%int_adiff, vd%int_adiff)
5728 0 : vd_max%mean_adiff = max(vd_max%mean_adiff, vd%mean_adiff)
5729 0 : vd_max%stdev_adiff = max(vd_max%stdev_adiff, vd%stdev_adiff)
5730 0 : vd_max%min_adiff = max(vd_max%min_adiff, vd%min_adiff)
5731 0 : vd_max%max_adiff = max(vd_max%max_adiff, vd%max_adiff)
5732 0 : vd_max%l1_rerr = max(vd_max%l1_rerr, vd%L1_rerr)
5733 : end if
5734 :
5735 0 : if (present(unit)) call vd%print(unit=unit)
5736 :
5737 0 : end subroutine vdiff_eval
5738 : !!***
5739 :
5740 : !----------------------------------------------------------------------
5741 :
5742 : !!****f* m_numeric_tools/vdiff_print
5743 : !! NAME
5744 : !! vdiff_print
5745 : !!
5746 : !! FUNCTION
5747 : !! Print vdiff_t to unit
5748 : !!
5749 : !! SOURCE
5750 :
5751 0 : subroutine vdiff_print(vd, unit)
5752 :
5753 : !Arguments ------------------------------------
5754 : !scalars
5755 : class(vdiff_t),intent(in) :: vd
5756 : integer,optional,intent(in) :: unit
5757 :
5758 : !Local variables-------------------------------
5759 : integer :: unt
5760 : ! *********************************************************************
5761 :
5762 0 : unt = std_out; if (present(unit)) unt = unit
5763 0 : write(unt,"(a,es10.3,a)")" L1_rerr: ", vd%l1_rerr, ","
5764 0 : write(unt,"(a,es10.3,a)")" 'Integral |f1-f2|dr': ", vd%int_adiff, ","
5765 0 : write(unt,"(a,es10.3,a)")" 'min {|f1-f2|}': ", vd%min_adiff, ","
5766 0 : write(unt,"(a,es10.3,a)")" 'Max {|f1-f2|}': ", vd%max_adiff, ","
5767 0 : write(unt,"(a,es10.3,a)")" 'mean {|f1-f2|}': ", vd%mean_adiff, ","
5768 0 : write(unt,"(a,es10.3,a)")" 'stdev {|f1-f2|}': ", vd%stdev_adiff, ","
5769 :
5770 0 : end subroutine vdiff_print
5771 : !!***
5772 :
5773 : !!****f* m_numeric_tools/smooth
5774 : !! NAME
5775 : !! smooth data.
5776 : !!
5777 : !! FUNCTION
5778 : !! smooth
5779 : !!
5780 : !! INPUTS
5781 : !! mesh=Number of points.
5782 : !! it=Number of iterations. <= 0 to return a unchanged.
5783 : !!
5784 : !! SIDE EFFECTS
5785 : !! a(mesh)=Input values, smoothed in output
5786 : !!
5787 : !! SOURCE
5788 :
5789 142 : pure subroutine smooth(a, mesh, it)
5790 :
5791 : !Arguments ------------------------------------
5792 : !scalars
5793 : integer, intent(in) :: it,mesh
5794 : real(dp), intent(inout) :: a(mesh)
5795 : !Local variables-------------------------------
5796 : integer :: i,k
5797 284 : real(dp) :: asm(mesh)
5798 : ! *********************************************************************
5799 :
5800 1172 : do k=1,it
5801 1030 : asm(1)=1.0d0/3.0d0*(a(1)+a(2)+a(3))
5802 1030 : asm(2)=0.25d0*(a(1)+a(2)+a(3)+a(4))
5803 1030 : asm(3)=0.2d0*(a(1)+a(2)+a(3)+a(4)+a(5))
5804 1030 : asm(4)=0.2d0*(a(2)+a(3)+a(4)+a(5)+a(6))
5805 1030 : asm(5)=0.2d0*(a(3)+a(4)+a(5)+a(6)+a(7))
5806 : asm(mesh-4)=0.2d0*(a(mesh-2)+a(mesh-3)+a(mesh-4)+&
5807 1030 : & a(mesh-5)+a(mesh-6))
5808 : asm(mesh-3)=0.2d0*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
5809 1030 : & a(mesh-4)+a(mesh-5))
5810 : asm(mesh-2)=0.2d0*(a(mesh)+a(mesh-1)+a(mesh-2)+&
5811 1030 : & a(mesh-3)+a(mesh-4))
5812 1030 : asm(mesh-1)=0.25d0*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
5813 1030 : asm(mesh)=1.0d0/3.0d0*(a(mesh)+a(mesh-1)+a(mesh-2))
5814 :
5815 2509730 : do i=6,mesh-5
5816 : asm(i)=0.1d0*a(i)+0.1d0*(a(i+1)+a(i-1))+&
5817 : & 0.1d0*(a(i+2)+a(i-2))+&
5818 : & 0.1d0*(a(i+3)+a(i-3))+&
5819 : & 0.1d0*(a(i+4)+a(i-4))+&
5820 2509730 : & 0.05d0*(a(i+5)+a(i-5))
5821 : end do
5822 :
5823 2520172 : do i=1,mesh
5824 2520030 : a(i)=asm(i)
5825 : end do
5826 : end do
5827 :
5828 142 : end subroutine smooth
5829 : !!***
5830 :
5831 : !!****f* m_numeric_tools/nderiv
5832 : !! NAME
5833 : !! nderiv
5834 : !!
5835 : !! FUNCTION
5836 : !! Given an input function y(x) on a regular grid,
5837 : !! compute its first or second derivative.
5838 : !!
5839 : !! INPUTS
5840 : !! hh= radial step
5841 : !! ndim= radial mesh size
5842 : !! yy(ndim)= input function
5843 : !! norder= order of derivation (1 or 2)
5844 : !!
5845 : !! OUTPUT
5846 : !! zz(ndim)= first or second derivative of y
5847 : !!
5848 : !! SOURCE
5849 :
5850 42 : pure subroutine nderiv(hh,yy,zz,ndim,norder)
5851 :
5852 : !Arguments ---------------------------------------------
5853 : !scalars
5854 : integer,intent(in) :: ndim,norder
5855 : real(dp),intent(in) :: hh
5856 : !arrays
5857 : real(dp),intent(in) :: yy(ndim)
5858 : real(dp),intent(out) :: zz(ndim)
5859 :
5860 : !Local variables ---------------------------------------
5861 : !scalars
5862 : integer :: ier,ii
5863 : real(dp) :: aa,bb,cc,h1,y1
5864 : ! *********************************************************************
5865 :
5866 : !Initialization (common to 1st and 2nd derivative)
5867 42 : h1=one/(12.d0*hh)
5868 42 : y1=yy(ndim-4)
5869 :
5870 : !FIRST DERIVATIVE
5871 : !================
5872 42 : if (norder==1) then
5873 :
5874 : ! Prepare differentiation loop
5875 42 : bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
5876 42 : cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
5877 : ! Start differentiation loop
5878 58938 : do ii=5,ndim
5879 58896 : aa=bb;bb=cc
5880 58896 : cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
5881 58938 : zz(ii-4)=aa
5882 : end do
5883 : ! Normal exit
5884 42 : ier=0
5885 42 : aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
5886 42 : 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))
5887 42 : zz(ndim-1)=aa
5888 42 : zz(ndim-2)=cc
5889 42 : zz(ndim-3)=bb
5890 :
5891 : ! SECOND DERIVATIVE
5892 : ! =================
5893 : else
5894 0 : h1=h1/hh
5895 : ! Prepare differentiation loop
5896 0 : bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
5897 0 : cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
5898 : ! Start differentiation loop
5899 0 : do ii=5,ndim
5900 0 : aa=bb;bb=cc
5901 0 : cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
5902 0 : zz(ii-4)=aa
5903 : end do
5904 : ! Normal exit
5905 0 : ier=0
5906 0 : aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
5907 0 : 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))
5908 0 : zz(ndim-1)=aa
5909 0 : zz(ndim-2)=cc
5910 0 : zz(ndim-3)=bb
5911 :
5912 : end if !norder
5913 :
5914 42 : end subroutine nderiv
5915 : !!***
5916 :
5917 : !----------------------------------------------------------------------
5918 :
5919 : !!****f* m_numeric_tools/central_finite_diff
5920 : !! NAME
5921 : !! central_finite_diff
5922 : !!
5923 : !! FUNCTION
5924 : !! Coefficients of the central differences, for several orders of accuracy.
5925 : !! See: https://en.wikipedia.org/wiki/Finite_difference_coefficient
5926 : !!
5927 : !! INPUTS
5928 : !! order=Derivative order.
5929 : !! ipos=Index of the point must be in [1,npts]
5930 : !! npts=Number of points used in finite difference, origin at npts/2 + 1
5931 : !!
5932 : !! OUTPUT
5933 : !! coefficients for central finite difference
5934 : !!
5935 : !! SOURCE
5936 :
5937 684 : real(dp) function central_finite_diff(order, ipos, npts) result(fact)
5938 :
5939 : !Arguments ---------------------------------------------
5940 : !scalars
5941 : integer,intent(in) :: ipos,order,npts
5942 :
5943 : !Local variables ---------------------------------------
5944 : !scalars
5945 : real(dp),parameter :: empty=huge(one)
5946 : ! 1st derivative.
5947 : real(dp),parameter :: d1(9,4) = reshape([ &
5948 : [-1/2._dp, 0._dp, 1/2._dp, empty, empty, empty, empty, empty, empty], &
5949 : [ 1/12._dp, -2/3._dp, 0._dp, 2/3._dp, -1/12._dp, empty, empty, empty, empty], &
5950 : [-1/60._dp, 3/20._dp, -3/4._dp, 0._dp, 3/4._dp, -3/20._dp, 1/60._dp, empty, empty], &
5951 : [ 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])
5952 : ! 2nd derivative.
5953 : real(dp),parameter :: d2(9,4) = reshape([ &
5954 : [ 1._dp, -2._dp, 1._dp, empty, empty, empty, empty, empty, empty], &
5955 : [-1/12._dp, 4/3._dp, -5/2._dp, 4/3._dp, -1/12._dp, empty, empty, empty, empty], &
5956 : [ 1/90._dp, -3/20._dp, 3/2._dp, -49/18._dp, 3/2._dp, -3/20._dp, 1/90._dp, empty, empty], &
5957 : [-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])
5958 : ! 3th derivative.
5959 : real(dp),parameter :: d3(9,3) = reshape([ &
5960 : [-1/2._dp, 1._dp, 0._dp, -1._dp, 1/2._dp, empty, empty, empty, empty], &
5961 : [ 1/8._dp, -1._dp, 13/8._dp, 0._dp, -13/8._dp, 1._dp, -1/8._dp, empty, empty], &
5962 : [ -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]], &
5963 : [9,3])
5964 : ! 4th derivative.
5965 : real(dp),parameter :: d4(9,3) = reshape([ &
5966 : [ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp, empty, empty, empty, empty], &
5967 : [ -1/6._dp, 2._dp, -13/2._dp, 28/3._dp, -13/2._dp, 2._dp, -1/6._dp, empty, empty], &
5968 : [ 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])
5969 : ! 5th derivative.
5970 : real(dp),parameter :: d5(7) = [ -1/2._dp, 2._dp, -5/2._dp, 0._dp, 5/2._dp, -2._dp, 1/2._dp]
5971 : ! 6th derivative.
5972 : real(dp),parameter :: d6(7) = [ 1._dp, -6._dp, 15._dp, -20._dp, 15._dp, -6._dp, 1._dp]
5973 : ! *********************************************************************
5974 :
5975 1368 : select case (order)
5976 : case (1)
5977 684 : if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5978 684 : fact = d1(ipos, npts/2)
5979 : case (2)
5980 0 : if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5981 0 : fact = d2(ipos, npts/2)
5982 : case (3)
5983 0 : if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5984 0 : fact = d3(ipos, npts/2)
5985 : case (4)
5986 0 : if (ipos < 1 .or. ipos > 9 .or. npts < 1 .or. npts > 9) goto 10
5987 0 : fact = d4(ipos, npts/2 - 1)
5988 : case (5)
5989 0 : if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5990 0 : fact = d5(ipos)
5991 : case (6)
5992 0 : if (ipos < 1 .or. ipos > 7 .or. npts /= 7) goto 10
5993 0 : fact = d6(ipos)
5994 : case default
5995 684 : ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
5996 : end select
5997 :
5998 684 : if (fact == empty) then
5999 0 : ABI_ERROR(sjoin("Invalid ipos:",itoa(ipos),"for order", itoa(order), "npts", itoa(npts)))
6000 : end if
6001 0 : return
6002 :
6003 0 : 10 ABI_ERROR(sjoin("No entry for ipos:",itoa(ipos),"order", itoa(order), "npts", itoa(npts)))
6004 :
6005 0 : end function central_finite_diff
6006 : !!***
6007 :
6008 : !!****f* m_numeric_tools/uniformrandom
6009 : !! NAME
6010 : !! uniformrandom
6011 : !!
6012 : !! FUNCTION
6013 : !! Returns a uniform random deviate between 0.0 and 1.0.
6014 : !! Set seed to any value < 0 to initialize or reinitialize sequence.
6015 : !! Parameters are chosen from integer overflow=2**23 (conservative).
6016 : !! For some documentation, see Numerical Recipes, 1986, p196.
6017 : !!
6018 : !! INPUTS
6019 : !!
6020 : !! OUTPUT
6021 : !!
6022 : !! NOTES
6023 : !!
6024 : !! SOURCE
6025 :
6026 1390581 : function uniformrandom(seed)
6027 :
6028 : !Arguments ------------------------------------
6029 : !scalars
6030 : real(dp) :: uniformrandom
6031 : integer,intent(inout) :: seed
6032 :
6033 : !Local variables ---------------------------------------
6034 : integer, parameter :: im1=11979,ia1= 430,ic1=2531
6035 : integer, parameter :: im2= 6655,ia2= 936,ic2=1399
6036 : integer, parameter :: im3= 6075,ia3=1366,ic3=1283
6037 : integer, save :: init=0
6038 : integer, save :: ii1,ii2,ii3
6039 : integer :: kk
6040 : real(dp) :: im1inv,im2inv
6041 : real(dp), save :: table(97)
6042 : character(len=500) :: msg
6043 : ! *********************************************************************
6044 :
6045 1390581 : im1inv=1.0d0/im1 ; im2inv=1.0d0/im2
6046 :
6047 : !Initialize on first call or when seed<0:
6048 1390581 : if (seed<0.or.init==0) then
6049 569 : seed=-abs(seed)
6050 :
6051 : ! First generator
6052 569 : ii1=mod(ic1-seed,im1)
6053 569 : ii1=mod(ia1*ii1+ic1,im1)
6054 : ! Second generator
6055 569 : ii2=mod(ii1,im2)
6056 569 : ii1=mod(ia1*ii1+ic1,im1)
6057 : ! Third generator
6058 569 : ii3=mod(ii1,im3)
6059 :
6060 : ! Fill table
6061 55762 : do kk=1,97
6062 55193 : ii1=mod(ia1*ii1+ic1,im1)
6063 55193 : ii2=mod(ia2*ii2+ic2,im2)
6064 55762 : table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
6065 : enddo
6066 :
6067 569 : init=1 ; seed=1
6068 : end if
6069 :
6070 : !Third generator gives index
6071 1390581 : ii3=mod(ia3*ii3+ic3,im3)
6072 1390581 : kk=1+(97*ii3)/im3
6073 1390581 : if (kk<1.or.kk>97) then
6074 0 : write(msg,'(a,2i0,a)' ) ' trouble in uniformrandom; ii3,kk=',ii3,kk,' =>stop'
6075 0 : ABI_ERROR(msg)
6076 : end if
6077 1390581 : uniformrandom=table(kk)
6078 :
6079 : !Replace old value, based on generators 1 and 2
6080 1390581 : ii1=mod(ia1*ii1+ic1,im1)
6081 1390581 : ii2=mod(ia2*ii2+ic2,im2)
6082 1390581 : table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
6083 :
6084 1390581 : end function uniformrandom
6085 : !!***
6086 :
6087 : !!****f* m_numeric_tools/findmin
6088 : !!
6089 : !! NAME
6090 : !! findmin
6091 : !!
6092 : !! FUNCTION
6093 : !! Compute the minimum of a function whose value and derivative are known at two points.
6094 : !! Also deduce different quantities at this predicted point, and at the two other points
6095 : !! It uses a quartic interpolation, with the supplementary
6096 : !! condition that the second derivative vanishes at one and
6097 : !! only one point. See Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].
6098 : !! For this option, lambda_1 must be 1 (new point),
6099 : !! and lambda_2 must be 0 (old point).
6100 : !! Also, if the derivative at the new point is more negative
6101 : !! than the derivative at the old point, the predicted
6102 : !! point cannot correspond to a minimum, but will be lambda=2.5_dp,
6103 : !! if the energy of the second point is lower than the energy
6104 : !! of the first point.
6105 : !!
6106 : !! INPUTS
6107 : !! etotal_1=first value of the function
6108 : !! etotal_2=second value of the function
6109 : !! dedv_1=first value of the derivative
6110 : !! dedv_2=second value of the derivative
6111 : !! lambda_1=first value of the argument
6112 : !! lambda_2=second value of the argument
6113 : !!
6114 : !! OUTPUT
6115 : !! dedv_predict=predicted value of the derivative (usually zero,
6116 : !! except if choice=4, if it happens that a minimum cannot be located,
6117 : !! and a trial step is taken)
6118 : !! d2edv2_predict=predicted value of the second derivative (not if choice=4)
6119 : !! d2edv2_1=first value of the second derivative (not if choice=4)
6120 : !! d2edv2_2=second value of the second derivative (not if choice=4)
6121 : !! etotal_predict=predicted value of the function
6122 : !! lambda_predict=predicted value of the argument
6123 : !! status= 0 if everything went normally ;
6124 : !! 1 if negative second derivative
6125 : !! 2 if some other problem
6126 : !!
6127 : !! SOURCE
6128 :
6129 247 : subroutine findmin(dedv_1,dedv_2,dedv_predict,&
6130 : & d2edv2_1,d2edv2_2,d2edv2_predict,&
6131 : & etotal_1,etotal_2,etotal_predict,&
6132 : & lambda_1,lambda_2,lambda_predict,status)
6133 :
6134 : !Arguments ------------------------------------
6135 : !scalars
6136 : integer,intent(out) :: status
6137 : real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2
6138 : real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict
6139 : real(dp),intent(out) :: etotal_predict,lambda_predict
6140 :
6141 : !Local variables-------------------------------
6142 : !scalars
6143 : real(dp) :: aa,bb,bbp,cc,ccp,d_lambda,dd
6144 : real(dp) :: discr,ee,eep,lambda_shift,sum1,sum2,sum3,uu
6145 : real(dp) :: uu3,vv,vv3
6146 : character(len=500) :: msg
6147 : ! *************************************************************************
6148 :
6149 : !write(std_out,*)' findmin : enter'
6150 : !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2
6151 :
6152 247 : status=0
6153 247 : d_lambda=lambda_1-lambda_2
6154 :
6155 : !DEBUG
6156 : !do choice=3,1,-1
6157 : !ENDDEBUG
6158 :
6159 247 : if(abs(lambda_1-1.0_dp)>tol12 .or. abs(lambda_2)>tol12) then
6160 0 : ABI_BUG('For choice=4, lambda_1 must be 1 and lambda_2 must be 0.')
6161 : end if
6162 :
6163 : !Evaluate quartic interpolation
6164 : !etotal = aa + bb * lambda + cc * lambda**2 + dd * lambda**3 + ee * lambda**4
6165 : !Impose positive second derivative everywhere, with
6166 : !one point where it vanishes : 3*dd**2=8*cc*ee
6167 247 : aa=etotal_2
6168 247 : bb=dedv_2
6169 247 : sum1=etotal_1-aa-bb
6170 247 : sum2=dedv_1-bb
6171 247 : sum3=sum2-2.0_dp*sum1
6172 :
6173 : !Build the discriminant of the associated 2nd degree equation
6174 247 : discr=sum2**2-3.0_dp*sum3**2
6175 247 : if(discr<0.0_dp .or. sum2<0.0_dp)then
6176 :
6177 : ! jmb init
6178 37 : d2edv2_2=0.0
6179 37 : d2edv2_1=0.0
6180 37 : d2edv2_predict=0.0
6181 :
6182 : ! Even if there is a problem, try to keep going ...
6183 37 : ABI_WARNING('The 2nd degree equation has no positive root (choice=4).')
6184 37 : status=2
6185 37 : if(etotal_1<etotal_2)then
6186 : write(msg, '(a,a,a)' )&
6187 37 : 'Will continue, since the new total energy is lower',ch10,&
6188 74 : 'than the old. Take a larger step in the same direction.'
6189 37 : ABI_COMMENT(msg)
6190 37 : lambda_predict=2.5_dp
6191 : else
6192 : write(msg, '(a,a,a,a,a)' )&
6193 0 : 'There is a problem, since the new total energy is larger',ch10,&
6194 0 : 'than the old (choice=4).',ch10,&
6195 0 : 'I take a point between the old and new, close to the old .'
6196 0 : ABI_COMMENT(msg)
6197 0 : lambda_predict=0.25_dp
6198 : end if
6199 : ! Mimic a zero-gradient lambda, in order to avoid spurious
6200 : ! action of the inverse hessian (the next line would be a realistic estimation)
6201 37 : dedv_predict=0.0_dp
6202 : ! dedv_predict=dedv_2+lambda_predict*(dedv_1-dedv_2)
6203 : ! Uses the energies, and the gradient at lambda_2
6204 : etotal_predict=etotal_2+dedv_2*lambda_predict&
6205 37 : & +(etotal_1-etotal_2-dedv_2)*lambda_predict**2
6206 :
6207 : else
6208 :
6209 : ! Here, there is an acceptable solution to the 2nd degree equation
6210 210 : discr=sqrt(discr)
6211 : ! The root that gives the smallest ee corresponds to -discr
6212 : ! This is the one to be used: one aims at modelling the
6213 : ! behaviour of the function as much as possible with the
6214 : ! lowest orders of the polynomial, not the quartic term.
6215 210 : ee=(sum2-discr)*0.5_dp
6216 210 : dd=sum3-2.0_dp*ee
6217 210 : cc=sum1-dd-ee
6218 :
6219 : ! DEBUG
6220 : ! write(std_out,*)'aa,bb,cc,dd,ee',aa,bb,cc,dd,ee
6221 : ! ENDDEBUG
6222 :
6223 : ! Now, must find the unique root of
6224 : ! 0 = bb + 2*cc * lambda + 3*dd * lambda^2 + 4*ee * lambda^3
6225 : ! This root is unique because it was imposed that the second derivative
6226 : ! of the quartic polynomial is everywhere positive.
6227 : ! First, remove the quadratic term, by a shift of lambda
6228 : ! lambdap=lambda-lambda_shift
6229 : ! 0 = bbp + ccp * lambdap + eep * lambdap^3
6230 210 : eep=4.0_dp*ee
6231 210 : lambda_shift=-dd/(4.0_dp*ee)
6232 210 : ccp=2.0_dp*cc-12.0_dp*ee*lambda_shift**2
6233 210 : bbp=bb+ccp*lambda_shift+eep*lambda_shift**3
6234 :
6235 : ! DEBUG
6236 : ! write(std_out,*)'bbp,ccp,eep,lambda_shift',bbp,ccp,eep,lambda_shift
6237 : ! ENDDEBUG
6238 :
6239 : ! The solution of a cubic polynomial equation is as follows :
6240 210 : discr=(bbp/eep)**2+(4.0_dp/27.0_dp)*(ccp/eep)**3
6241 : ! In the present case, discr will always be positive
6242 210 : discr=sqrt(discr)
6243 210 : uu3=0.5_dp*(-bbp/eep+discr) ; uu=sign((abs(uu3))**(1.0_dp/3.0_dp),uu3)
6244 210 : vv3=0.5_dp*(-bbp/eep-discr) ; vv=sign((abs(vv3))**(1.0_dp/3.0_dp),vv3)
6245 210 : lambda_predict=uu+vv
6246 :
6247 : ! Restore the shift
6248 210 : lambda_predict=lambda_predict+lambda_shift
6249 : etotal_predict=aa+bb*lambda_predict+cc*lambda_predict**2+&
6250 210 : & dd*lambda_predict**3+ee*lambda_predict**4
6251 : dedv_predict=bb+2.0_dp*cc*lambda_predict+3.0_dp*dd*lambda_predict**2+&
6252 210 : & 4.0_dp*ee*lambda_predict**3
6253 210 : d2edv2_1=2*cc+6*dd*lambda_1+12*ee*lambda_1**2
6254 210 : d2edv2_2=2*cc+6*dd*lambda_2+12*ee*lambda_2**2
6255 210 : d2edv2_predict=2*cc+6*dd*lambda_predict+12*ee*lambda_predict**2
6256 :
6257 : end if
6258 :
6259 247 : write(msg, '(a,i3)' )' line minimization, algorithm ',4
6260 247 : call wrtout(std_out,msg,'COLL')
6261 247 : write(msg, '(a,a)' )' lambda etotal ',' dedv d2edv2 '
6262 247 : call wrtout(std_out,msg,'COLL')
6263 247 : write(msg, '(a,es12.4,es18.10,2es12.4)' )' old point :',lambda_2,etotal_2,dedv_2,d2edv2_2
6264 247 : call wrtout(std_out,msg,'COLL')
6265 247 : write(msg, '(a,es12.4,es18.10,2es12.4)' )' new point :',lambda_1,etotal_1,dedv_1,d2edv2_1
6266 247 : call wrtout(std_out,msg,'COLL')
6267 247 : write(msg, '(a,es12.4,es18.10,2es12.4)' )' predicted point :',lambda_predict,etotal_predict,dedv_predict,d2edv2_predict
6268 247 : call wrtout(std_out,msg,'COLL')
6269 247 : write(msg, '(a)' ) ' '
6270 247 : call wrtout(std_out,msg,'COLL')
6271 :
6272 247 : end subroutine findmin
6273 : !!***
6274 :
6275 : !!****f* m_numeric_tools/kramerskronig
6276 : !! NAME
6277 : !! kramerskronig
6278 : !!
6279 : !! FUNCTION
6280 : !! check or apply the Kramers Kronig relation:
6281 : !! Re \epsilon(\omega) = 1 + \frac{2}{\pi}
6282 : !! \int_0^\infty d\omega' frac{\omega'}{\omega'^2 - \omega^2} Im \epsilon(\omega')
6283 : !!
6284 : !! INPUTS
6285 : !! nomega=number of real frequencies
6286 : !! omega(nomega)= real frequencies
6287 : !! eps(nomega)= function on the frequency grid (both real and imaginary part)
6288 : !! real part can be used to check whether the K-K relation is satisfied or not
6289 : !! method=method used to perform the integration
6290 : !! 0= naive integration
6291 : !! 1=simpson rule
6292 : !! only_check= if /=0 the real part of eps is checked against the imaginary part,
6293 : !! a final report in written but the array eps is not modified
6294 : !! if ==0 the real part of eps is overwritten using the
6295 : !! results obtained using the Kramers-Kronig relation
6296 : !!
6297 : !! OUTPUT
6298 : !!
6299 : !! NOTES
6300 : !! Inspired to check_kramerskronig of the DP code
6301 : !!
6302 : !! SOURCE
6303 :
6304 0 : subroutine kramerskronig(nomega,omega,eps,method,only_check)
6305 :
6306 : !Arguments ------------------------------------
6307 : !scalars
6308 : integer,intent(in) :: method,nomega,only_check
6309 : !arrays
6310 : real(dp),intent(in) :: omega(nomega)
6311 : complex(dp),intent(inout) :: eps(nomega)
6312 :
6313 : !Local variables-------------------------------
6314 : !scalars
6315 : integer,save :: enough=0
6316 : integer :: ii,ip
6317 : real(dp) :: acc,domega,eav,kkdif,kkrms,ww,wwp
6318 : character(len=500) :: msg
6319 : !arrays
6320 0 : real(dp) :: e1kk(nomega),intkk(nomega),kk(nomega)
6321 : ! *************************************************************************
6322 :
6323 : !Check whether the frequency grid is linear or not
6324 0 : domega = (omega(nomega) - omega(1)) / (nomega-1)
6325 0 : do ii=2,nomega
6326 0 : if (ABS(domega-(omega(ii)-omega(ii-1))) > 0.001) then
6327 0 : if (only_check/=1) then
6328 0 : ABI_WARNING("Check cannot be performed since the frequency step is not constant")
6329 0 : RETURN
6330 : else
6331 0 : ABI_ERROR('Cannot perform integration since frequency step is not constant')
6332 : end if
6333 : end if
6334 : end do
6335 :
6336 : !Check whether omega(1) is small or not
6337 0 : if (omega(1) > 0.1/Ha_eV) then
6338 0 : if (only_check/=1) then
6339 0 : ABI_WARNING('Check cannot be performed since first frequency on the grid > 0.1 eV')
6340 0 : RETURN
6341 : else
6342 0 : ABI_ERROR('Cannot perform integration since first frequency on the grid > 0.1 eV')
6343 : end if
6344 : end if
6345 :
6346 : !If eps(nomega) is not 0 warn
6347 0 : if (AIMAG(eps(nomega)) > 0.1 .and. enough<50) then
6348 0 : enough=enough+1
6349 : write(msg,'(a,f8.4,3a,f8.2,2a)')&
6350 0 : & 'Im epsilon for omega = ',omega(nomega)*Ha_eV,' eV',ch10,&
6351 0 : & 'is not yet zero, epsilon_2 = ',AIMAG(eps(nomega)),ch10,&
6352 0 : & 'Kramers Kronig could give wrong results'
6353 0 : ABI_WARNING(msg)
6354 0 : if (enough==50) then
6355 0 : write(msg,'(3a)')' sufficient number of WARNINGS-',ch10,' stop writing '
6356 0 : call wrtout(std_out,msg,'COLL')
6357 : end if
6358 : end if
6359 :
6360 : !Perform Kramers-Kronig using naive integration
6361 0 : select case (method)
6362 : case (0)
6363 :
6364 0 : do ii=1,nomega
6365 0 : ww = omega(ii)
6366 0 : acc = 0.0_dp
6367 0 : do ip=1,nomega
6368 0 : if (ip == ii) CYCLE
6369 0 : wwp = omega(ip)
6370 0 : acc = acc + wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6371 : end do
6372 0 : e1kk(ii) = one + two/pi*domega* acc
6373 : end do
6374 :
6375 : ! Perform Kramers-Kronig using Simpson integration
6376 : ! Simpson O(1/N^4), from NumRec in C p 134 NumRec in Fortran p 128
6377 : case (1)
6378 :
6379 0 : kk=zero
6380 :
6381 0 : do ii=1,nomega
6382 0 : ww=omega(ii)
6383 0 : do ip=1,nomega
6384 0 : if (ip == ii) CYCLE
6385 0 : wwp = omega(ip)
6386 0 : kk(ip) = wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
6387 : end do
6388 0 : call simpson_int(nomega,domega,kk,intkk)
6389 0 : e1kk(ii) = one + two/pi * intkk(nomega)
6390 : end do
6391 :
6392 : case default
6393 0 : write(msg,'(a,i0)')' Wrong value for method ',method
6394 0 : ABI_BUG(msg)
6395 : end select
6396 :
6397 : !at this point real part is in e1kk, need to put it into eps
6398 0 : do ii=1,nomega
6399 0 : eps(ii)=CMPLX(e1kk(ii),AIMAG(eps(ii)), kind=dp)
6400 : end do
6401 :
6402 : !Verify Kramers-Kronig
6403 : eav = zero
6404 : kkdif = zero
6405 : kkrms = zero
6406 :
6407 : do ii=1,nomega
6408 : kkdif = kkdif + ABS(REAL(eps(ii)) - e1kk(ii))
6409 : kkrms = kkrms + (REAL(eps(ii)) - e1kk(ii))*(REAL(eps(ii)) - e1kk(ii))
6410 : eav = eav + ABS(REAL(eps(ii)))
6411 : end do
6412 :
6413 0 : eav = eav/nomega
6414 0 : kkdif = (kkdif/nomega) / eav
6415 0 : kkrms = (kkrms/nomega) / (eav*eav)
6416 :
6417 0 : kk = ABS(REAL(eps(1)) - e1kk(1)) / REAL(eps(1))
6418 :
6419 : !Write data
6420 0 : write(msg,'(a,f7.2,a)')' Kramers-Kronig transform is verified within ',MAXVAL(kk)*100,"%"
6421 0 : call wrtout(std_out,msg,'COLL')
6422 :
6423 : end subroutine kramerskronig
6424 : !!***
6425 :
6426 : !!****f* ABINIT/dotproduct
6427 : !! NAME
6428 : !! dotproduct
6429 : !!
6430 : !! FUNCTION
6431 : !! scalar product of two vectors
6432 : !!
6433 : !! INPUTS
6434 : !! v1 and v2: two real(dp) vectors
6435 : !!
6436 : !! OUTPUT
6437 : !! scalar product of the two vectors
6438 : !!
6439 : !! WARNINGS
6440 : !! vector size is not checked
6441 : !!
6442 : !! NOTES
6443 : !! I've benchmarked this to be speedier than the intrinsic dot_product even on
6444 : !! big vectors. The point is that less check is performed.
6445 : !!
6446 : !! MG: FIXME: Well, optized blas1 is for sure better than what you wrote!
6447 : !! Now I don't have time to update ref files
6448 : !!
6449 : !! SOURCE
6450 :
6451 8555 : function dotproduct(nv1,nv2,v1,v2)
6452 :
6453 : !Arguments ------------------------------------
6454 : !scalars
6455 : integer,intent(in) :: nv1,nv2
6456 : real(dp) :: dotproduct
6457 : !arrays
6458 : real(dp),intent(in) :: v1(nv1,nv2),v2(nv1,nv2)
6459 :
6460 : !Local variables-------------------------------
6461 : !scalars
6462 : integer :: i,j
6463 : ! *************************************************************************
6464 :
6465 8555 : dotproduct=zero
6466 17110 : do j=1,nv2
6467 85567110 : do i=1,nv1
6468 85558555 : dotproduct=dotproduct+v1(i,j)*v2(i,j)
6469 : end do
6470 : end do
6471 :
6472 8555 : end function dotproduct
6473 : !!***
6474 :
6475 : !!****f* m_numeric_tools/invcb
6476 : !! NAME
6477 : !! invcb
6478 : !!
6479 : !! FUNCTION
6480 : !! Compute a set of inverse cubic roots as fast as possible:
6481 : !! rspts(:)=rhoarr(:)$^\frac{-1}{3}$
6482 : !!
6483 : !! INPUTS
6484 : !! npts=number of real space points on which density is provided
6485 : !! rhoarr(npts)=input data
6486 : !!
6487 : !! OUTPUT
6488 : !! rspts(npts)=inverse cubic root of rhoarr
6489 : !!
6490 : !! SOURCE
6491 :
6492 1268480 : subroutine invcb(rhoarr,rspts,npts)
6493 :
6494 : !Arguments ------------------------------------
6495 : !scalars
6496 : integer,intent(in) :: npts
6497 : !arrays
6498 : real(dp),intent(in) :: rhoarr(npts)
6499 : real(dp),intent(out) :: rspts(npts)
6500 :
6501 : !Local variables-------------------------------
6502 : !scalars
6503 : integer :: ii,ipts
6504 : real(dp),parameter :: c2_27=2.0e0_dp/27.0e0_dp,c5_9=5.0e0_dp/9.0e0_dp
6505 : real(dp),parameter :: c8_9=8.0e0_dp/9.0e0_dp,m1thrd=-third
6506 : real(dp) :: del,prod,rho,rhom1,rhomtrd
6507 : logical :: test
6508 : !character(len=500) :: message
6509 : ! *************************************************************************
6510 :
6511 : !Loop over points : here, brute force algorithm
6512 : !do ipts=1,npts
6513 : !rspts(ipts)=sign( (abs(rhoarr(ipts)))**m1thrd,rhoarr(ipts))
6514 : !end do
6515 :
6516 1268480 : rhomtrd=sign( (abs(rhoarr(1)))**m1thrd, rhoarr(1) )
6517 1268480 : rhom1=one/rhoarr(1)
6518 1268480 : rspts(1)=rhomtrd
6519 2261450933 : do ipts=2,npts
6520 2260182453 : rho=rhoarr(ipts)
6521 2260182453 : prod=rho*rhom1
6522 : ! If the previous point is too far ...
6523 2260182453 : if(prod < 0.01_dp .or. prod > 10._dp )then
6524 9296719 : rhomtrd=sign( (abs(rho))**m1thrd , rho )
6525 9296719 : rhom1=one/rho
6526 : else
6527 2250885734 : del=prod-one
6528 3761148362 : do ii=1,5
6529 : ! Choose one of the two next lines, the last one is more accurate
6530 : ! rhomtrd=((one+third*del)/(one+two_thirds*del))*rhomtrd
6531 3761148362 : rhomtrd=((one+c5_9*del)/(one+del*(c8_9+c2_27*del)))*rhomtrd
6532 3761148362 : rhom1=rhomtrd*rhomtrd*rhomtrd
6533 3761148362 : del=rho*rhom1-one
6534 : ! write(std_out,*)rhomtrd,del
6535 3761148362 : test = del*del < 1.0e-24_dp
6536 3761148362 : if(test) exit
6537 : end do
6538 2250885734 : if( .not. test) then
6539 0 : rhomtrd=sign( (abs(rho))**m1thrd , rho )
6540 : end if
6541 : end if
6542 2261450933 : rspts(ipts)=rhomtrd
6543 : end do
6544 :
6545 1268480 : end subroutine invcb
6546 : !!***
6547 :
6548 : !!****f* ABINIT/safe_div
6549 : !! NAME
6550 : !! safe_div
6551 : !!
6552 : !! FUNCTION
6553 : !! Subroutine safe_div performs "safe division", that is to prevent overflow,
6554 : !! underflow, NaN, or infinity errors. An alternate value is returned if the
6555 : !! division cannot be performed. (bmy, 2/26/08)
6556 : !!
6557 : !! For more information, see the discussion on:
6558 : !! http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/
6559 : !!
6560 : !! Taken by HM from:
6561 : !! http://wiki.seas.harvard.edu/geos-chem/index.php/Floating_point_math_issues#Safe_floating-point_division
6562 : !!
6563 : !! INPUTS
6564 : !! n : Numerator for the division
6565 : !! d : Divisor for the division
6566 : !! altv : Alternate value to be returned if the division can't be done
6567 : !!
6568 : !! OUTPUT
6569 : !!
6570 : !! SOURCE
6571 :
6572 896914 : elemental subroutine safe_div(n, d, altv, q)
6573 :
6574 : !Arguments ----------------------------------------------
6575 : !scalars
6576 : real(dp),intent(in) :: n, d, altv
6577 : real(dp),intent(out) :: q
6578 : ! *********************************************************************
6579 :
6580 896914 : if ( exponent(n) - exponent(d) >= maxexponent(n) .or. d == zero) then
6581 659320 : q = altv
6582 : else
6583 237594 : q = n / d
6584 : endif
6585 :
6586 896914 : end subroutine safe_div
6587 : !!***
6588 :
6589 : !!****f* ABINIT/bool2index
6590 : !! NAME
6591 : !! bool2index
6592 : !!
6593 : !! FUNCTION
6594 : !! Allocate and return array with the indices in the input boolean array `bool_list` that evaluates to .True.
6595 : !!
6596 : !! SOURCE
6597 :
6598 0 : subroutine bool2index(bool_list, out_index)
6599 :
6600 : !Arguments ----------------------------------------------
6601 : !scalars
6602 : logical,intent(in) :: bool_list(:)
6603 : integer,allocatable,intent(inout) :: out_index(:)
6604 :
6605 : !Local variables-------------------------------
6606 : integer :: ii, cnt
6607 : ! *********************************************************************
6608 :
6609 0 : cnt = count(bool_list)
6610 0 : ABI_REMALLOC(out_index, (cnt))
6611 0 : cnt = 0
6612 0 : do ii=1,size(bool_list)
6613 0 : if (bool_list(ii)) then
6614 0 : cnt = cnt + 1
6615 0 : out_index(cnt) = ii
6616 : end if
6617 : end do
6618 :
6619 0 : end subroutine bool2index
6620 : !!***
6621 :
6622 : !----------------------------------------------------------------------
6623 : !!****f* ABINIT/polynomial_regression
6624 : !! NAME
6625 : !! polynomial_regression
6626 : !!
6627 : !! FUNCTION
6628 : !! Perform a polynomial regression on incoming data points, the
6629 : !! x-values of which are stored in array xvals and the y-values
6630 : !! stored in array yvals. Returns a one dimensional array with
6631 : !! fit coefficients (coeffs) and the unbiased RMS error of the
6632 : !! fit as a scalar (RMSerr).
6633 : !!
6634 : !! INPUTS
6635 : !! degree = order of the polynomial
6636 : !! npts = number of data points
6637 : !! xvals(npts) = x-values of those data points
6638 : !! yvals(npts) = y-values of those data points
6639 : !!
6640 : !! OUTPUT
6641 : !! coeffs(degree+1) = coefficients of the polynomial regression
6642 : !! RMSerr = unbiased RMS error on the fit
6643 : !! RMSerr=\sqrt{\frac{1}{npts-1}*
6644 : !! \sum_i^npts{(fitval-yvals(i))**2}}
6645 : !!
6646 : !! SOURCE
6647 : !!
6648 :
6649 12 : subroutine polynomial_regression(degree,npts,xvals,yvals,coeffs,RMSerr)
6650 :
6651 : !Arguments ------------------------------------
6652 : !scalars
6653 : integer :: degree,npts
6654 : real(dp),intent(out) :: RMSerr
6655 : !arrays
6656 : real(dp),intent(in) :: xvals(1:npts),yvals(1:npts)
6657 : real(dp),intent(out) :: coeffs(degree+1)
6658 :
6659 : !Local variables-------------------------------
6660 : !scalars
6661 : integer :: ncoeffs,icoeff,ipoint,info
6662 : real(dp) :: residual,fitval
6663 : !arrays
6664 12 : integer,allocatable :: tmp(:)
6665 12 : real(dp),allocatable :: tmptwo(:), A(:,:),ATA(:,:)
6666 :
6667 : !####################################################################
6668 : !##################### Get Polynomial Fit #########################
6669 :
6670 12 : ncoeffs=degree+1
6671 :
6672 36 : ABI_MALLOC(tmp,(ncoeffs))
6673 36 : ABI_MALLOC(tmptwo,(ncoeffs))
6674 48 : ABI_MALLOC(A,(npts,ncoeffs))
6675 48 : ABI_MALLOC(ATA,(ncoeffs,ncoeffs))
6676 :
6677 : !Construct a polynomial for all input xvalues
6678 48 : do icoeff=1,ncoeffs
6679 246 : do ipoint=1,npts
6680 234 : if (icoeff==1.and.xvals(ipoint)==0.0) then
6681 12 : A(ipoint,icoeff) = 1.0
6682 : else
6683 186 : A(ipoint,icoeff) = xvals(ipoint)**(icoeff-1)
6684 : end if
6685 : end do
6686 : end do
6687 :
6688 : !Get matrix product of transpose of A and A
6689 954 : ATA = matmul(transpose(A),A)
6690 :
6691 : !Compute LU factorization of ATA
6692 12 : call DGETRF(ncoeffs,ncoeffs,ATA,ncoeffs,tmp,info)
6693 12 : ABI_CHECK(info == 0, sjoin('LAPACK DGETRF in polynomial regression returned:', itoa(info)))
6694 :
6695 : !Compute inverse of the LU factorized version of ATA
6696 12 : call DGETRI(ncoeffs,ATA,ncoeffs,tmp,tmptwo,ncoeffs,info)
6697 12 : ABI_CHECK(info == 0, sjoin('LAPACK DGETRI in polynomial regression returned:', itoa(info)))
6698 :
6699 : !Harvest polynomial coefficients
6700 2164 : coeffs = matmul(matmul(ATA,transpose(A)),yvals)
6701 :
6702 : !####################################################################
6703 : !############## RMS error on the polynomial fit ###################
6704 :
6705 12 : residual=0.0d0
6706 78 : do ipoint=1,npts
6707 66 : fitval=0.0d0
6708 264 : do icoeff=1,ncoeffs
6709 264 : if (icoeff==1.and.xvals(ipoint)==0.0) then
6710 12 : fitval=fitval+coeffs(icoeff)
6711 : else
6712 186 : fitval=fitval+coeffs(icoeff)*xvals(ipoint)**(icoeff-1)
6713 : end if
6714 : end do
6715 78 : residual=residual+(fitval-yvals(ipoint))**2
6716 : end do
6717 12 : RMSerr=sqrt(residual/(real(npts-1,8)))
6718 :
6719 12 : ABI_FREE(A)
6720 12 : ABI_FREE(ATA)
6721 12 : ABI_FREE(tmp)
6722 12 : ABI_FREE(tmptwo)
6723 :
6724 12 : end subroutine polynomial_regression
6725 : !!***
6726 :
6727 : !!****f* ABINIT/blocked_loop
6728 : !! NAME
6729 : !! blocked_loop
6730 : !!
6731 : !! FUNCTION
6732 : !! Helper function to implement blocked algorithms inside do loops i.e. algorithms
6733 : !! operating on multiple items up to a maximum `batch_size`.
6734 : !!
6735 : !! Usage:
6736 : !!
6737 : !! batch_size = 4
6738 : !! allocate(work(..., batch_size)
6739 : !!
6740 : !! do loop_index=1, loop_stop, batch_size
6741 : !! ndat = blocked_loop(loop_index, loop_stop, batch_size)
6742 : !! ! operate on ndat items in work
6743 : !! end do
6744 : !!
6745 : !! SOURCE
6746 :
6747 0 : integer pure function blocked_loop(loop_index, loop_stop, batch_size) result(ndat)
6748 :
6749 : !Arguments ----------------------------------------------
6750 : integer,intent(in) :: loop_index, loop_stop, batch_size
6751 : ! *********************************************************************
6752 :
6753 0 : ndat = merge(batch_size, loop_stop - loop_index + 1, loop_index + batch_size - 1 <= loop_stop)
6754 :
6755 0 : end function blocked_loop
6756 : !!***
6757 :
6758 : !!****f* m_euler/geteuler
6759 : !! NAME
6760 : !! geteuler
6761 : !!
6762 : !! FUNCTION
6763 : !! Compute the Euler angles (alpha, beta) corresponding to the spin quantization axis given in Cartesian coordinates.
6764 : !!
6765 : !! INPUTS
6766 : !! spinaxis(3)=spin quantization axis
6767 : !!
6768 : !! OUTPUT
6769 : !! alpha=Euler angle for rotation around z-axis
6770 : !! beta =Euler angle for rotation around y-axis
6771 : !!
6772 : !! SOURCE
6773 :
6774 2342216 : subroutine geteuler(spinaxis, alpha, beta)
6775 :
6776 : !Arguments -------------------------------
6777 : !scalars
6778 : real(dp),intent(out) :: alpha, beta
6779 : !arrays
6780 : real(dp),intent(in) :: spinaxis(3)
6781 :
6782 : !Local variables -------------------------
6783 : !scalars
6784 : real(dp) :: sx, sy, sz, norm, rxy
6785 : !***********************************************************************
6786 :
6787 2342216 : alpha = zero; beta = zero
6788 9368864 : norm = DOT_PRODUCT(spinaxis, spinaxis)
6789 :
6790 2342216 : if (norm <= tol8*tol8) return
6791 :
6792 1239444 : sx = spinaxis(1); sy = spinaxis(2); sz = spinaxis(3)
6793 1239444 : rxy = sqrt(sx*sx + sy*sy)
6794 1239444 : if (rxy > tol8) alpha = atan2(sy, sx)
6795 1239444 : beta = atan2(rxy, sz)
6796 :
6797 : end subroutine geteuler
6798 : !!***
6799 :
6800 : !====================================================================
6801 : ! bin_t: average x-value inside the bin
6802 : !====================================================================
6803 : !function bin_average(self, xvals) result(avg)
6804 : ! class(bin_t), intent(in) :: self
6805 : ! real(dp), intent(in) :: xvals(:)
6806 : ! real(dp) :: avg
6807 : ! integer :: i
6808 : !
6809 : ! if (self%npts == 0) then
6810 : ! avg = 0.0_dp
6811 : ! return
6812 : ! end if
6813 : !
6814 : ! avg = 0.0_dp
6815 : ! do i = 1, self%npts
6816 : ! avg = avg + xvals(self%idx(i))
6817 : ! end do
6818 : ! avg = avg / real(self%npts, dp)
6819 : !end function bin_average
6820 :
6821 : !====================================================================
6822 : ! bins_t initialisation: only non-empty bins are stored
6823 : !====================================================================
6824 0 : subroutine bins_init(self, nn, xvals, dx)
6825 : class(bins_t), intent(out) :: self
6826 : integer,intent(in) :: nn
6827 : real(dp), target, intent(in) :: xvals(nn)
6828 : real(dp), intent(in) :: dx
6829 :
6830 : integer :: i, k, b, idx_bin
6831 : real(dp) :: xmin, xmax, length
6832 0 : integer, allocatable :: count_(:), map(:)
6833 :
6834 0 : self%xvals => xvals
6835 0 : self%dx = dx
6836 0 : self%total_points = nn
6837 :
6838 0 : xmin = minval(xvals)
6839 0 : xmax = maxval(xvals)
6840 0 : self%xmin = xmin
6841 0 : self%xmax = xmax
6842 :
6843 0 : length = xmax - xmin
6844 :
6845 : ! Number of uniform bins
6846 0 : k = int(floor(length/dx)) + 1
6847 0 : ABI_MALLOC(count_, (k))
6848 0 : count_ = 0
6849 :
6850 : ! First pass: count
6851 0 : do i = 1, nn
6852 0 : b = 1 + int( (xvals(i) - xmin) / dx )
6853 0 : if (b < 1) b = 1
6854 0 : if (b > k) b = k
6855 0 : count_(b) = count_(b) + 1
6856 : end do
6857 :
6858 : ! Map from full bin list to non-empty bins
6859 0 : ABI_MALLOC(map, (k))
6860 0 : map = 0
6861 0 : self%nbins = count(count_ > 0)
6862 0 : ABI_MALLOC(self%bin, (self%nbins))
6863 :
6864 : ! Fill map
6865 0 : idx_bin = 0
6866 0 : do b = 1, k
6867 0 : if (count_(b) > 0) then
6868 0 : idx_bin = idx_bin + 1
6869 0 : map(b) = idx_bin
6870 0 : ABI_MALLOC(self%bin(idx_bin)%idx, (count_(b)))
6871 0 : self%bin(idx_bin)%npts = count_(b)
6872 0 : self%bin(idx_bin)%xmin = xmin + (b-1)*dx
6873 0 : self%bin(idx_bin)%xmax = xmin + b *dx
6874 : end if
6875 : end do
6876 :
6877 : ! Temporary counters
6878 0 : count_ = 0
6879 :
6880 : ! Second pass: put indices into non-empty bins
6881 0 : do i = 1, nn
6882 0 : b = 1 + int( (xvals(i) - xmin)/dx )
6883 0 : if (b < 1) b = 1
6884 0 : if (b > k) b = k
6885 0 : if (map(b) > 0) then
6886 0 : idx_bin = map(b)
6887 0 : count_(b) = count_(b) + 1
6888 0 : self%bin(idx_bin)%idx(count_(b)) = i
6889 : end if
6890 : end do
6891 :
6892 : ! Simple loop interface
6893 : !do i = 1, self%nbins
6894 : ! write(*,*) "Bin", i, ": npts=", self%bin(i)%npts
6895 : !end do
6896 :
6897 0 : end subroutine bins_init
6898 :
6899 : !!****t* m_numeric_tools/bins_free
6900 : !! NAME
6901 : !! bins_free
6902 : !!
6903 : !! FUNCTION
6904 : !!
6905 : !! SOURCE
6906 :
6907 : subroutine bins_free(bins)
6908 :
6909 : class(bins_t),intent(inout) :: bins
6910 : !!************************************************************************
6911 : integer :: ii
6912 : do ii = 1, bins%nbins
6913 : ABI_SFREE(bins%bin(ii)%idx)
6914 : end do
6915 : ABI_SFREE(bins%bin)
6916 : !
6917 : end subroutine bins_free
6918 : !!***
6919 :
6920 12 : end module m_numeric_tools
6921 : !!***
|
