TABLE OF CONTENTS


ABINIT/m_gaussian_quadrature [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gaussian_quadrature

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_gaussian_quadrature
28 !----------------------------------------------------------------------------------------------------
29 !
30 ! This module contains routines to generate the weights and positions at which to evaluate
31 ! an integrand in order to produce a Gaussian quadrature of a given type.
32 !
33 ! The code here was obtained from the internet (see header to subroutine cdgqf for further details).
34 ! The code was slightly modified to use doubles instead of floats in order to reach greater
35 ! accuracy.
36 !
37 !----------------------------------------------------------------------------------------------------
38 
39  use m_abicore
40  use m_errors
41 
42  use defs_basis, only : dp, tol20, std_out
43  use m_io_tools, only : open_file
44 
45  implicit none
46 
47  private

m_hamiltonian/cdgqf [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  cdgqf

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

383 subroutine cdgqf ( nt, gaussian_kind, alpha, beta, t, wts )
384 
385 !*****************************************************************************80
386 !
387 !! CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
388 !
389 !  Discussion:
390 !
391 !    This routine computes all the knots and weights of a Gauss quadrature
392 !    formula with a classical weight function with default values for A and B,
393 !    and only simple knots.
394 !
395 !    There are no moments checks and no printing is done.
396 !
397 !    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
398 !
399 !  Licensing:
400 !
401 !    This code is distributed under the GNU LGPL license.
402 !
403 !  Modified:
404 !
405 !    04 January 2010
406 !
407 !  Author:
408 !
409 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
410 !    FORTRAN90 version by John Burkardt.
411 !
412 !  Reference:
413 !
414 !    Sylvan Elhay, Jaroslav Kautsky,
415 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
416 !    Interpolatory Quadrature,
417 !    ACM Transactions on Mathematical Software,
418 !    Volume 13, Number 4, December 1987, pages 399-415.
419 !
420 !  Parameters:
421 !
422 !    Input, integer ( kind = 4 ) NT, the number of knots.
423 !
424 !    Input, integer ( kind = 4 ) KIND, the rule.
425 !    1, Legendre,             (a,b)       1.0
426 !    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
427 !    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
428 !    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
429 !    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
430 !    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
431 !    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
432 !    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
433 !
434 !    Input, real ( dp ) ALPHA, the value of Alpha, if needed.
435 !
436 !    Input, real ( dp ) BETA, the value of Beta, if needed.
437 !
438 !    Output, real ( dp ) T(NT), the knots.
439 !
440 !    Output, real ( dp ) WTS(NT), the weights.
441 !
442 
443 !This section has been created automatically by the script Abilint (TD).
444 !Do not modify the following lines by hand.
445 #undef ABI_FUNC
446 #define ABI_FUNC 'cdgqf'
447 !End of the abilint section
448 
449   implicit none
450 
451   integer ( kind = 4 ) nt
452 
453   real ( dp ) aj(nt)
454   real ( dp ) alpha
455   real ( dp ) beta
456   real ( dp ) bj(nt)
457   integer ( kind = 4 ) gaussian_kind
458   real ( dp ) t(nt)
459   real ( dp ) wts(nt)
460   real ( dp ) zemu
461 
462 ! *************************************************************************
463 
464   call parchk ( gaussian_kind, 2 * nt, alpha, beta )
465 !
466 !  Get the Jacobi matrix and zero-th moment.
467 !
468   call class_matrix ( gaussian_kind, nt, alpha, beta, aj, bj, zemu )
469 !
470 !  Compute the knots and weights.
471 !
472   call sgqf ( nt, aj, bj, zemu, t, wts )
473 
474   return
475 end subroutine cdgqf

m_hamiltonian/cgqf [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  cgqf

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_efmas,m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

497 subroutine cgqf ( nt, gaussian_kind, alpha, beta, a, b, t, wts )
498 
499 !*****************************************************************************80
500 !
501 !! CGQF computes knots and weights of a Gauss quadrature formula.
502 !
503 !  Discussion:
504 !
505 !    The user may specify the interval (A,B).
506 !
507 !    Only simple knots are produced.
508 !
509 !    Use routine EIQFS to evaluate this quadrature formula.
510 !
511 !  Licensing:
512 !
513 !    This code is distributed under the GNU LGPL license.
514 !
515 !  Modified:
516 !
517 !    16 February 2010
518 !
519 !  Author:
520 !
521 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
522 !    FORTRAN90 version by John Burkardt.
523 !
524 !  Reference:
525 !
526 !    Sylvan Elhay, Jaroslav Kautsky,
527 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
528 !    Interpolatory Quadrature,
529 !    ACM Transactions on Mathematical Software,
530 !    Volume 13, Number 4, December 1987, pages 399-415.
531 !
532 !  Parameters:
533 !
534 !    Input, integer ( kind = 4 ) NT, the number of knots.
535 !
536 !    Input, integer ( kind = 4 ) KIND, the rule.
537 !    1, Legendre,             (a,b)       1.0
538 !    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^-0.5)
539 !    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
540 !    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
541 !    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
542 !    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
543 !    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
544 !    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
545 !    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
546 !
547 !    Input, real ( dp ) ALPHA, the value of Alpha, if needed.
548 !
549 !    Input, real ( dp ) BETA, the value of Beta, if needed.
550 !
551 !    Input, real ( dp ) A, B, the interval endpoints, or
552 !    other parameters.
553 !
554 !    Output, real ( dp ) T(NT), the knots.
555 !
556 !    Output, real ( dp ) WTS(NT), the weights.
557 !
558 
559 !This section has been created automatically by the script Abilint (TD).
560 !Do not modify the following lines by hand.
561 #undef ABI_FUNC
562 #define ABI_FUNC 'cgqf'
563 !End of the abilint section
564 
565   implicit none
566 
567   integer ( kind = 4 ) nt
568 
569   real ( dp ) a
570   real ( dp ) alpha
571   real ( dp ) b
572   real ( dp ) beta
573   integer ( kind = 4 ) i
574   integer ( kind = 4 ) gaussian_kind
575   integer ( kind = 4 ), allocatable :: mlt(:)
576   integer ( kind = 4 ), allocatable :: ndx(:)
577   real ( dp ) t(nt)
578   real ( dp ) wts(nt)
579 
580 ! *************************************************************************
581 
582 !
583 !  Compute the Gauss quadrature formula for default values of A and B.
584 !
585   call cdgqf ( nt, gaussian_kind, alpha, beta, t, wts )
586 !
587 !  Prepare to scale the quadrature formula to other weight function with
588 !  valid A and B.
589 !
590   ABI_ALLOCATE ( mlt, (1:nt) )
591 
592   mlt(1:nt) = 1
593 
594   ABI_ALLOCATE ( ndx, (1:nt) )
595 
596   do i = 1, nt
597     ndx(i) = i
598   end do
599 
600   call scqf ( nt, t, mlt, wts, nt, ndx, wts, t, gaussian_kind, alpha, beta, a, b )
601 
602   ABI_DEALLOCATE ( mlt )
603   ABI_DEALLOCATE ( ndx )
604 
605   return
606 end subroutine cgqf

m_hamiltonian/ch_cap [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  ch_cap

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

628 subroutine ch_cap ( c )
629 
630 !*****************************************************************************80
631 !
632 !! CH_CAP capitalizes a single character.
633 !
634 !  Licensing:
635 !
636 !    This code is distributed under the GNU LGPL license.
637 !
638 !  Modified:
639 !
640 !    19 July 1998
641 !
642 !  Author:
643 !
644 !    John Burkardt
645 !
646 !  Parameters:
647 !
648 !    Input/output, character C, the character to capitalize.
649 !
650 
651 !This section has been created automatically by the script Abilint (TD).
652 !Do not modify the following lines by hand.
653 #undef ABI_FUNC
654 #define ABI_FUNC 'ch_cap'
655 !End of the abilint section
656 
657   implicit none
658 
659   character              c
660   integer ( kind = 4 ) itemp
661 
662 ! *************************************************************************
663 
664   itemp = ichar ( c )
665 
666   if ( 97 <= itemp .and. itemp <= 122 ) then
667     c = char ( itemp - 32 )
668   end if
669 
670   return
671 end subroutine ch_cap

m_hamiltonian/ch_eqi [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  ch_eqi

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

693 function ch_eqi ( c1, c2 )
694 
695 !*****************************************************************************80
696 !
697 !! CH_EQI is a case insensitive comparison of two characters for equality.
698 !
699 !  Example:
700 !
701 !    CH_EQI ( 'A', 'a' ) is .TRUE.
702 !
703 !  Licensing:
704 !
705 !    This code is distributed under the GNU LGPL license.
706 !
707 !  Modified:
708 !
709 !    28 July 2000
710 !
711 !  Author:
712 !
713 !    John Burkardt
714 !
715 !  Parameters:
716 !
717 !    Input, character C1, C2, the characters to compare.
718 !
719 !    Output, logical CH_EQI, the result of the comparison.
720 !
721 
722 !This section has been created automatically by the script Abilint (TD).
723 !Do not modify the following lines by hand.
724 #undef ABI_FUNC
725 #define ABI_FUNC 'ch_eqi'
726 !End of the abilint section
727 
728   implicit none
729 
730   logical ch_eqi
731   character c1
732   character c1_cap
733   character c2
734   character c2_cap
735 
736 ! *************************************************************************
737 
738   c1_cap = c1
739   c2_cap = c2
740 
741   call ch_cap ( c1_cap )
742   call ch_cap ( c2_cap )
743 
744   if ( c1_cap == c2_cap ) then
745     ch_eqi = .true.
746   else
747     ch_eqi = .false.
748   end if
749 
750   return
751 end function ch_eqi

m_hamiltonian/ch_to_digit [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  ch_to_digit

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

773 subroutine ch_to_digit ( c, digit )
774 
775 !*****************************************************************************80
776 !
777 !! CH_TO_DIGIT returns the integer value of a base 10 digit.
778 !
779 !  Example:
780 !
781 !     C   DIGIT
782 !    ---  -----
783 !    '0'    0
784 !    '1'    1
785 !    ...  ...
786 !    '9'    9
787 !    ' '    0
788 !    'X'   -1
789 !
790 !  Licensing:
791 !
792 !    This code is distributed under the GNU LGPL license.
793 !
794 !  Modified:
795 !
796 !    04 August 1999
797 !
798 !  Author:
799 !
800 !    John Burkardt
801 !
802 !  Parameters:
803 !
804 !    Input, character C, the decimal digit, '0' through '9' or blank
805 !    are legal.
806 !
807 !    Output, integer ( kind = 4 ) DIGIT, the corresponding integer value.
808 !    If C was 'illegal', then DIGIT is -1.
809 !
810 
811 !This section has been created automatically by the script Abilint (TD).
812 !Do not modify the following lines by hand.
813 #undef ABI_FUNC
814 #define ABI_FUNC 'ch_to_digit'
815 !End of the abilint section
816 
817   implicit none
818 
819   character c
820   integer ( kind = 4 ) digit
821 
822 ! *************************************************************************
823 
824   if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then
825 
826     digit = ichar ( c ) - 48
827 
828   else if ( c == ' ' ) then
829 
830     digit = 0
831 
832   else
833 
834     digit = -1
835 
836   end if
837 
838   return
839 end subroutine ch_to_digit

m_hamiltonian/class_matrix [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  class_matrix

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

 861 subroutine class_matrix ( gaussian_kind, m, alpha, beta, aj, bj, zemu )
 862 
 863 !*****************************************************************************80
 864 !
 865 !! CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
 866 !
 867 !  Discussion:
 868 !
 869 !    This routine computes the diagonal AJ and sub-diagonal BJ
 870 !    elements of the order M tridiagonal symmetric Jacobi matrix
 871 !    associated with the polynomials orthogonal with respect to
 872 !    the weight function specified by KIND.
 873 !
 874 !    For weight functions 1-7, M elements are defined in BJ even
 875 !    though only M-1 are needed.  For weight function 8, BJ(M) is
 876 !    set to zero.
 877 !
 878 !    The zero-th moment of the weight function is returned in ZEMU.
 879 !
 880 !  Licensing:
 881 !
 882 !    This code is distributed under the GNU LGPL license.
 883 !
 884 !  Modified:
 885 !
 886 !    27 December 2009
 887 !
 888 !  Author:
 889 !
 890 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
 891 !    FORTRAN90 version by John Burkardt.
 892 !
 893 !  Reference:
 894 !
 895 !    Sylvan Elhay, Jaroslav Kautsky,
 896 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
 897 !    Interpolatory Quadrature,
 898 !    ACM Transactions on Mathematical Software,
 899 !    Volume 13, Number 4, December 1987, pages 399-415.
 900 !
 901 !  Parameters:
 902 !
 903 !    Input, integer ( kind = 4 ) KIND, the rule.
 904 !    1, Legendre,             (a,b)       1.0
 905 !    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
 906 !    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
 907 !    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
 908 !    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
 909 !    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
 910 !    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
 911 !    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
 912 !
 913 !    Input, integer ( kind = 4 ) M, the order of the Jacobi matrix.
 914 !
 915 !    Input, real ( dp ) ALPHA, the value of Alpha, if needed.
 916 !
 917 !    Input, real ( dp ) BETA, the value of Beta, if needed.
 918 !
 919 !    Output, real ( dp ) AJ(M), BJ(M), the diagonal and subdiagonal
 920 !    of the Jacobi matrix.
 921 !
 922 !    Output, real ( dp ) ZEMU, the zero-th moment.
 923 !
 924 
 925 !This section has been created automatically by the script Abilint (TD).
 926 !Do not modify the following lines by hand.
 927 #undef ABI_FUNC
 928 #define ABI_FUNC 'class_matrix'
 929 !End of the abilint section
 930 
 931   implicit none
 932 
 933   integer ( kind = 4 ) m
 934 
 935   real ( dp ) a2b2
 936   real ( dp ) ab
 937   real ( dp ) aba
 938   real ( dp ) abi
 939   real ( dp ) abj
 940   real ( dp ) abti
 941   real ( dp ) aj(m)
 942   real ( dp ) alpha
 943   real ( dp ) apone
 944   real ( dp ) beta
 945   real ( dp ) bj(m)
 946   integer ( kind = 4 ) i
 947   integer ( kind = 4 ) gaussian_kind
 948   real ( dp ), parameter :: pi = 3.14159265358979323846264338327950D+00
 949   real ( dp ) temp
 950   real ( dp ) temp2
 951   real ( dp ) zemu
 952 
 953 ! *************************************************************************
 954 
 955   temp = epsilon ( temp )
 956 
 957   call parchk ( gaussian_kind, 2 * m - 1, alpha, beta )
 958 
 959   temp2 = 0.5D+00
 960 
 961   if ( 500.0D+00 * temp < abs ( ( r8_gamma_gq ( temp2 ) )**2 - pi ) ) then
 962     MSG_ERROR('Gamma function does not match machine parameters.')
 963   end if
 964 
 965   if ( gaussian_kind == 1 ) then
 966 
 967     ab = 0.0D+00
 968 
 969     zemu = 2.0D+00 / ( ab + 1.0D+00 )
 970 
 971     aj(1:m) = 0.0D+00
 972 
 973     do i = 1, m
 974       abi = i + ab * mod ( i, 2 )
 975       abj = 2 * i + ab
 976       bj(i) = abi * abi / ( abj * abj - 1.0D+00 )
 977     end do
 978     bj(1:m) =  sqrt ( bj(1:m) )
 979 
 980   else if ( gaussian_kind == 2 ) then
 981 
 982     zemu = pi
 983 
 984     aj(1:m) = 0.0D+00
 985 
 986     bj(1) =  sqrt ( 0.5D+00 )
 987     bj(2:m) = 0.5D+00
 988 
 989   else if ( gaussian_kind == 3 ) then
 990 
 991     ab = alpha * 2.0D+00
 992     zemu = 2.0D+00**( ab + 1.0D+00 ) * r8_gamma_gq ( alpha + 1.0D+00 )**2 &
 993       / r8_gamma_gq ( ab + 2.0D+00 )
 994 
 995     aj(1:m) = 0.0D+00
 996     bj(1) = 1.0D+00 / ( 2.0D+00 * alpha + 3.0D+00 )
 997     do i = 2, m
 998       bj(i) = i * ( i + ab ) / ( 4.0D+00 * ( i + alpha )**2 - 1.0D+00 )
 999     end do
1000     bj(1:m) =  sqrt ( bj(1:m) )
1001 
1002   else if ( gaussian_kind == 4 ) then
1003 
1004     ab = alpha + beta
1005     abi = 2.0D+00 + ab
1006     zemu = 2.0D+00**( ab + 1.0D+00 ) * r8_gamma_gq ( alpha + 1.0D+00 ) &
1007       * r8_gamma_gq ( beta + 1.0D+00 ) / r8_gamma_gq ( abi )
1008     aj(1) = ( beta - alpha ) / abi
1009     bj(1) = 4.0D+00 * ( 1.0 + alpha ) * ( 1.0D+00 + beta ) &
1010       / ( ( abi + 1.0D+00 ) * abi * abi )
1011     a2b2 = beta * beta - alpha * alpha
1012 
1013     do i = 2, m
1014       abi = 2.0D+00 * i + ab
1015       aj(i) = a2b2 / ( ( abi - 2.0D+00 ) * abi )
1016       abi = abi**2
1017       bj(i) = 4.0D+00 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) &
1018         / ( ( abi - 1.0D+00 ) * abi )
1019     end do
1020     bj(1:m) =  sqrt ( bj(1:m) )
1021 
1022   else if ( gaussian_kind == 5 ) then
1023 
1024     zemu = r8_gamma_gq ( alpha + 1.0D+00 )
1025 
1026     do i = 1, m
1027       aj(i) = 2.0D+00 * i - 1.0D+00 + alpha
1028       bj(i) = i * ( i + alpha )
1029     end do
1030     bj(1:m) =  sqrt ( bj(1:m) )
1031 
1032   else if ( gaussian_kind == 6 ) then
1033 
1034     zemu = r8_gamma_gq ( ( alpha + 1.0D+00 ) / 2.0D+00 )
1035 
1036     aj(1:m) = 0.0D+00
1037 
1038     do i = 1, m
1039       bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0D+00
1040     end do
1041     bj(1:m) =  sqrt ( bj(1:m) )
1042 
1043   else if ( gaussian_kind == 7 ) then
1044 
1045     ab = alpha
1046     zemu = 2.0D+00 / ( ab + 1.0D+00 )
1047 
1048     aj(1:m) = 0.0D+00
1049 
1050     do i = 1, m
1051       abi = i + ab * mod ( i, 2 )
1052       abj = 2 * i + ab
1053       bj(i) = abi * abi / ( abj * abj - 1.0D+00 )
1054     end do
1055     bj(1:m) =  sqrt ( bj(1:m) )
1056 
1057   else if ( gaussian_kind == 8 ) then
1058 
1059     ab = alpha + beta
1060     zemu = r8_gamma_gq ( alpha + 1.0D+00 ) * r8_gamma_gq ( - ( ab + 1.0D+00 ) ) &
1061       / r8_gamma_gq ( - beta )
1062     apone = alpha + 1.0D+00
1063     aba = ab * apone
1064     aj(1) = - apone / ( ab + 2.0D+00 )
1065     bj(1) = - aj(1) * ( beta + 1.0D+00 ) / ( ab + 2.0D+00 ) / ( ab + 3.0D+00 )
1066     do i = 2, m
1067       abti = ab + 2.0D+00 * i
1068       aj(i) = aba + 2.0D+00 * ( ab + i ) * ( i - 1 )
1069       aj(i) = - aj(i) / abti / ( abti - 2.0D+00 )
1070     end do
1071 
1072     do i = 2, m - 1
1073       abti = ab + 2.0D+00 * i
1074       bj(i) = i * ( alpha + i ) / ( abti - 1.0D+00 ) * ( beta + i ) &
1075         / ( abti**2 ) * ( ab + i ) / ( abti + 1.0D+00 )
1076     end do
1077 
1078     bj(m) = 0.0D+00
1079     bj(1:m) =  sqrt ( bj(1:m) )
1080 
1081   end if
1082 
1083   return
1084 end subroutine class_matrix

m_hamiltonian/gaussian_quadrature_gegenbauer [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  gaussian_quadrature_gegenbauer

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

      imtqlx

SOURCE

 77  subroutine gaussian_quadrature_gegenbauer(integrand,integral,order)
 78  !----------------------------------------------------------------------------------------------------
 79  !
 80  ! This subroutine performs a gaussian quadrature of given order of the input function integrand
 81  ! and stores the result in the variable integral.
 82  !
 83  ! The integrand function is assumed to be defined on the range from 0 to infinity, and it is
 84  ! assumed to to behave as x^2 near x=0, and 1/x^4 as x-> infinity.
 85  !
 86  ! The Gegenbauer integration method is of the form:
 87  !
 88  ! int_{-1}^{1} dt (1-t)^2(1+t)^2 f(t) =  sum_{i} w_i f(t_i)
 89  !
 90  !
 91  ! In order to perform such an integral, make a change of variable x = (1-t)/(1+t); this implies
 92  !
 93  !              int_0^infty dx I(x) = int_{-1}^{1} dt g(t),     g(t) = 2/(1+t)^2 I((1-t)/(1+t))
 94  !
 95  ! It is easy to show that g(1-delta) ~ delta^2 and g(-1+delta) ~ delta^2 for delta small.
 96  ! Thus, we can define
 97  !
 98  !              int_{-1}^{1} dt g(t) = int_{-1}^{1} dt (1+t)^2(1-t)^2 f(t), f(t) = g(t)/(1-t^2)^2
 99  !
100  !
101  !----------------------------------------------------------------------------------------------------
102 
103 !This section has been created automatically by the script Abilint (TD).
104 !Do not modify the following lines by hand.
105 #undef ABI_FUNC
106 #define ABI_FUNC 'gaussian_quadrature_gegenbauer'
107 !End of the abilint section
108 
109   implicit none
110 
111 
112   ! input and output variables
113   integer,   intent(in)  :: order
114 
115   interface
116         function integrand(x)
117         integer, parameter :: dp=kind(1.0d0)
118         real(dp) integrand
119         real(dp),intent(in) :: x
120         end function integrand
121   end interface
122 
123   real(dp),  intent(out) :: integral
124 
125 
126   real(dp) :: t(order), w(order)
127 
128 
129   integer  :: i
130 
131   integer  :: gaussian_kind
132   real(dp) :: alpha, beta, a, b
133   real(dp) :: x, y, f
134 
135 ! *************************************************************************
136 
137 
138   ! Set the arguments and weights
139   gaussian_kind = 3
140   alpha=  2.0_dp
141   beta =  2.0_dp
142   a    = -1.0_dp
143   b    =  1.0_dp
144   call cgqf ( order, gaussian_kind, alpha, beta, a, b, t, w )
145 
146 
147 
148   ! call the integrant the right number of times
149 
150   integral = 0.0_dp
151 
152   do i = 1, order
153         x = (1.0_dp-t(i))/(1.0_dp+t(i))
154 
155         y = integrand(x)
156 
157         f = 2.0_dp/(1.0_dp+t(i))**2/(1.0_dp-t(i)**2)**2*y
158 
159         integral = integral + w(i)*f
160 
161   end do ! i
162 
163  end subroutine gaussian_quadrature_gegenbauer

m_hamiltonian/gaussian_quadrature_legendre [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  gaussian_quadrature_legendre

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

      imtqlx

SOURCE

184  subroutine gaussian_quadrature_legendre(integrand,integral,order)
185  !----------------------------------------------------------------------------------------------------
186  !
187  ! This subroutine performs a gaussian quadrature of given order of the input function integrand
188  ! and stores the result in the variable integral.
189  !
190  ! The integrand function is assumed to be defined on the range from 0 to infinity, and it is
191  ! assumed to to behave as x^2 near x=0, and 1/x^4 as x-> infinity.
192  !
193  ! The Legendre integration method is of the form:
194  !
195  ! int_{-1}^{1} dt f(t) =  sum_{i} w_i f(t_i)
196  !
197  !
198  ! In order to perform such an integral, make a change of variable x = (1-t)/(1+t); this implies
199  !
200  !              int_0^infty dx I(x) = int_{-1}^{1} dt f(t),     f(t) = 2/(1+t)^2 I((1-t)/(1+t))
201  !
202  !
203  !
204  !----------------------------------------------------------------------------------------------------
205 
206 !This section has been created automatically by the script Abilint (TD).
207 !Do not modify the following lines by hand.
208 #undef ABI_FUNC
209 #define ABI_FUNC 'gaussian_quadrature_legendre'
210 !End of the abilint section
211 
212   implicit none
213 
214 
215   ! input and output variables
216   integer,   intent(in)  :: order
217 
218   interface
219         function integrand(x)
220         integer, parameter :: dp=kind(1.0d0)
221         real(dp) integrand
222         real(dp),intent(in) :: x
223         end function integrand
224   end interface
225 
226   real(dp),  intent(out) :: integral
227 
228 
229   real(dp) :: t(order), w(order)
230 
231 
232   integer  :: i
233 
234   integer  :: gaussian_kind
235   real(dp) :: alpha, beta, a, b
236   real(dp) :: x, y, f
237 
238 ! *************************************************************************
239 
240 
241   ! Set the arguments and weights
242   gaussian_kind = 1
243   alpha=  0.0_dp
244   beta =  0.0_dp
245   a    = -1.0_dp
246   b    =  1.0_dp
247   call cgqf ( order, gaussian_kind, alpha, beta, a, b, t, w )
248 
249 
250 
251   ! call the integrant the right number of times
252 
253   integral = 0.0_dp
254 
255   do i = 1, order
256         x = (1.0_dp-t(i))/(1.0_dp+t(i))
257 
258         y = integrand(x)
259 
260         f = 2.0_dp/(1.0_dp+t(i))**2*y
261 
262         integral = integral + w(i)*f
263 
264   end do ! i
265 
266  end subroutine gaussian_quadrature_legendre

m_hamiltonian/get_frequencies_and_weights_legendre [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  get_frequencies_and_weights_legendre

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_DielectricArray

CHILDREN

      imtqlx

SOURCE

289  subroutine get_frequencies_and_weights_legendre(order,frequencies,weights)
290  !----------------------------------------------------------------------------------------------------
291  !
292  ! This subroutine generates the frequencies and weights to integrate using legendre Gaussian
293  ! quadrature for a given order.
294  !
295  ! The integrand function is assumed to be defined on the range from 0 to infinity, and it is
296  ! assumed to to behave as x^2 near x=0, and 1/x^4 as x-> infinity.
297  !
298  ! The Legendre integration method is of the form:
299  !
300  ! int_{-1}^{1} dt f(t) =  sum_{i} w_i f(t_i)
301  !
302  ! In order to perform such an integral, make a change of variable
303  !
304  !              omega = (1-t)/(1+t)
305  !
306  ! and thus
307  !       int_0^infty domega I(omega) = int_{-1}^{1} dt f(t),    f(t) = 2/(1+t)^2 I((1-t)/(1+t))
308  !
309  ! and so the weights are defined as weights(i) = 2/(1+ti)^2*wi.
310  !
311  !
312  !----------------------------------------------------------------------------------------------------
313 
314 !This section has been created automatically by the script Abilint (TD).
315 !Do not modify the following lines by hand.
316 #undef ABI_FUNC
317 #define ABI_FUNC 'get_frequencies_and_weights_legendre'
318 !End of the abilint section
319 
320   implicit none
321   ! input and output variables
322   integer,   intent(in)   :: order
323   real(dp),  intent(out)  :: frequencies(order)
324   real(dp),  intent(out)  :: weights(order)
325 
326 
327   real(dp) :: t(order), w(order)
328 
329 
330   integer  :: i
331 
332   integer  :: gaussian_kind
333   real(dp) :: alpha, beta, a, b
334   real(dp) :: x,  f
335 
336 ! *************************************************************************
337 
338 
339   ! Set the arguments and weights
340   gaussian_kind = 1
341   alpha =  0.0_dp
342   beta  =  0.0_dp
343   a     = -1.0_dp
344   b     =  1.0_dp
345 
346   call cgqf ( order, gaussian_kind, alpha, beta, a, b, t, w )
347 
348 
349   ! call the integrant the right number of times
350 
351   do i = 1, order
352 
353         x = (1.0_dp-t(i))/(1.0_dp+t(i))
354         f = 2.0_dp/(1.0_dp+t(i))**2
355 
356         frequencies(i) = x
357         weights(i)     = w(i)*f
358 
359   end do ! i
360 
361  end subroutine get_frequencies_and_weights_legendre

m_hamiltonian/imtqlx [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  imtqlx

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

1106 subroutine imtqlx ( n, d, e, z )
1107 
1108 !*****************************************************************************80
1109 !
1110 !! IMTQLX diagonalizes a symmetric tridiagonal matrix.
1111 !
1112 !  Discussion:
1113 !
1114 !    This routine is a slightly modified version of the EISPACK routine to
1115 !    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
1116 !
1117 !    The authors thank the authors of EISPACK for permission to use this
1118 !    routine.
1119 !
1120 !    It has been modified to produce the product Q' * Z, where Z is an input
1121 !    vector and Q is the orthogonal matrix diagonalizing the input matrix.
1122 !    The changes consist (essentially) of applying the orthogonal
1123 !    transformations directly to Z as they are generated.
1124 !
1125 !  Licensing:
1126 !
1127 !    This code is distributed under the GNU LGPL license.
1128 !
1129 !  Modified:
1130 !
1131 !    27 December 2009
1132 !
1133 !  Author:
1134 !
1135 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
1136 !    FORTRAN90 version by John Burkardt.
1137 !
1138 !  Reference:
1139 !
1140 !    Sylvan Elhay, Jaroslav Kautsky,
1141 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
1142 !    Interpolatory Quadrature,
1143 !    ACM Transactions on Mathematical Software,
1144 !    Volume 13, Number 4, December 1987, pages 399-415.
1145 !
1146 !    Roger Martin, James Wilkinson,
1147 !    The Implicit QL Algorithm,
1148 !    Numerische Mathematik,
1149 !    Volume 12, Number 5, December 1968, pages 377-383.
1150 !
1151 !  Parameters:
1152 !
1153 !    Input, integer ( kind = 4 ) N, the order of the matrix.
1154 !
1155 !    Input/output, real ( dp ) D(N), the diagonal entries of the matrix.
1156 !    On output, the information in D has been overwritten.
1157 !
1158 !    Input/output, real ( dp ) E(N), the subdiagonal entries of the
1159 !    matrix, in entries E(1) through E(N-1).  On output, the information in
1160 !    E has been overwritten.
1161 !
1162 !    Input/output, real ( dp ) Z(N).  On input, a vector.  On output,
1163 !    the value of Q' * Z, where Q is the matrix that diagonalizes the
1164 !    input symmetric tridiagonal matrix.
1165 !
1166 
1167 !This section has been created automatically by the script Abilint (TD).
1168 !Do not modify the following lines by hand.
1169 #undef ABI_FUNC
1170 #define ABI_FUNC 'imtqlx'
1171 !End of the abilint section
1172 
1173   implicit none
1174 
1175   integer ( kind = 4 ) n
1176 
1177   real ( dp ) b
1178   real ( dp ) c
1179   real ( dp ) d(n)
1180   real ( dp ) e(n)
1181   real ( dp ) f
1182   real ( dp ) g
1183   integer ( kind = 4 ) i
1184   integer ( kind = 4 ) ii
1185   integer ( kind = 4 ), parameter :: itn = 30
1186   integer ( kind = 4 ) j
1187   integer ( kind = 4 ) k
1188   integer ( kind = 4 ) l
1189   integer ( kind = 4 ) m
1190   integer ( kind = 4 ) mml
1191   real ( dp ) p
1192   real ( dp ) prec
1193   real ( dp ) r
1194   real ( dp ) s
1195   real ( dp ) z(n)
1196 
1197 ! *************************************************************************
1198 
1199   prec = epsilon ( prec )
1200 
1201   if ( n == 1 ) then
1202     return
1203   end if
1204 
1205   e(n) = 0.0D+00
1206 
1207   do l = 1, n
1208 
1209     j = 0
1210 
1211     do
1212 
1213       do m = l, n
1214 
1215         if ( m == n ) then
1216           exit
1217         end if
1218 
1219         if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) ) then
1220           exit
1221         end if
1222 
1223       end do
1224 
1225       p = d(l)
1226 
1227       if ( m == l ) then
1228         exit
1229       end if
1230 
1231       if ( itn <= j ) then
1232         write ( std_out, '(a)' ) ' '
1233         write ( std_out, '(a)' ) 'IMTQLX - Fatal error!'
1234         write ( std_out, '(a)' ) '  Iteration limit exceeded.'
1235         write ( std_out, '(a,i8)' ) '  J = ', j
1236         write ( std_out, '(a,i8)' ) '  L = ', l
1237         write ( std_out, '(a,i8)' ) '  M = ', m
1238         write ( std_out, '(a,i8)' ) '  N = ', n
1239         MSG_ERROR("Aborting now")
1240       end if
1241 
1242       j = j + 1
1243       g = ( d(l+1) - p ) / ( 2.0D+00 * e(l) )
1244       r =  sqrt ( g * g + 1.0D+00 )
1245       g = d(m) - p + e(l) / ( g + sign ( r, g ) )
1246       s = 1.0D+00
1247       c = 1.0D+00
1248       p = 0.0D+00
1249       mml = m - l
1250 
1251       do ii = 1, mml
1252 
1253         i = m - ii
1254         f = s * e(i)
1255         b = c * e(i)
1256 
1257         if ( abs ( g ) <= abs ( f ) ) then
1258           c = g / f
1259           r =  sqrt ( c * c + 1.0D+00 )
1260           e(i+1) = f * r
1261           s = 1.0D+00 / r
1262           c = c * s
1263         else
1264           s = f / g
1265           r =  sqrt ( s * s + 1.0D+00 )
1266           e(i+1) = g * r
1267           c = 1.0D+00 / r
1268           s = s * c
1269         end if
1270 
1271         g = d(i+1) - p
1272         r = ( d(i) - g ) * s + 2.0D+00 * c * b
1273         p = s * r
1274         d(i+1) = g + p
1275         g = c * r - b
1276         f = z(i+1)
1277         z(i+1) = s * z(i) + c * f
1278         z(i) = c * z(i) - s * f
1279 
1280       end do
1281 
1282       d(l) = d(l) - p
1283       e(l) = g
1284       e(m) = 0.0D+00
1285 
1286     end do
1287 
1288   end do
1289 !
1290 !  Sorting.
1291 !
1292   do ii = 2, n
1293 
1294     i = ii - 1
1295     k = i
1296     p = d(i)
1297 
1298     do j = ii, n
1299       if ( d(j) < p ) then
1300         k = j
1301         p = d(j)
1302       end if
1303     end do
1304 
1305     if ( k /= i ) then
1306       d(k) = d(i)
1307       d(i) = p
1308       p = z(i)
1309       z(i) = z(k)
1310       z(k) = p
1311     end if
1312 
1313   end do
1314 
1315   return
1316 end subroutine imtqlx

m_hamiltonian/parchk [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  parchk

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

1338 subroutine parchk ( gaussian_kind, m, alpha, beta )
1339 
1340 !*****************************************************************************80
1341 !
1342 !! PARCHK checks parameters ALPHA and BETA for classical weight functions.
1343 !
1344 !  Licensing:
1345 !
1346 !    This code is distributed under the GNU LGPL license.
1347 !
1348 !  Modified:
1349 !
1350 !    27 December 2009
1351 !
1352 !  Author:
1353 !
1354 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
1355 !    FORTRAN90 version by John Burkardt.
1356 !
1357 !  Reference:
1358 !
1359 !    Sylvan Elhay, Jaroslav Kautsky,
1360 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
1361 !    Interpolatory Quadrature,
1362 !    ACM Transactions on Mathematical Software,
1363 !    Volume 13, Number 4, December 1987, pages 399-415.
1364 !
1365 !  Parameters:
1366 !
1367 !    Input, integer ( kind = 4 ) KIND, the rule.
1368 !    1, Legendre,             (a,b)       1.0
1369 !    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
1370 !    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
1371 !    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
1372 !    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
1373 !    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
1374 !    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
1375 !    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
1376 !
1377 !    Input, integer ( kind = 4 ) M, the order of the highest moment to
1378 !    be calculated.  This value is only needed when KIND = 8.
1379 !
1380 !    Input, real ( dp ) ALPHA, BETA, the parameters, if required
1381 !    by the value of KIND.
1382 !
1383 
1384 !This section has been created automatically by the script Abilint (TD).
1385 !Do not modify the following lines by hand.
1386 #undef ABI_FUNC
1387 #define ABI_FUNC 'parchk'
1388 !End of the abilint section
1389 
1390   implicit none
1391 
1392   real ( dp ) alpha
1393   real ( dp ) beta
1394   integer ( kind = 4 ) gaussian_kind
1395   integer ( kind = 4 ) m
1396   real ( dp ) tmp
1397 
1398 ! *************************************************************************
1399 
1400   if ( gaussian_kind <= 0 ) then
1401     MSG_ERROR('PARCHK - Fatal error: KIND <= 0.')
1402   end if
1403 !
1404 !  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
1405 !
1406   if ( 3 <= gaussian_kind .and. alpha <= -1.0D+00 ) then
1407     MSG_ERROR('PARCHK - Fatal error! 3 <= KIND and ALPHA <= -1.')
1408   end if
1409 !
1410 !  Check BETA for Jacobi.
1411 !
1412   if ( gaussian_kind == 4 .and. beta <= -1.0D+00 ) then
1413     MSG_ERROR('PARCHK - Fatal error: KIND == 4 and BETA <= -1.0.')
1414   end if
1415 !
1416 !  Check ALPHA and BETA for rational.
1417 !
1418   if ( gaussian_kind == 8 ) then
1419     tmp = alpha + beta + m + 1.0D+00
1420     if ( 0.0D+00 <= tmp .or. tmp <= beta ) then
1421       MSG_ERROR('PARCHK - Fatal error!  KIND == 8 but condition on ALPHA and BETA fails.')
1422     end if
1423   end if
1424 
1425   return
1426 end subroutine parchk

m_hamiltonian/r8_gamma_gq [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  r8_gamma_gq

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

1449 function r8_gamma_gq ( x )
1450 
1451 !*****************************************************************************80
1452 !
1453 !! R8_GAMMA evaluates Gamma(X) for a real argument.
1454 !
1455 !  Discussion:
1456 !
1457 !    This routine calculates the gamma function for a real argument X.
1458 !
1459 !    Computation is based on an algorithm outlined in reference 1.
1460 !    The program uses rational functions that approximate the gamma
1461 !    function to at least 20 significant decimal digits.  Coefficients
1462 !    for the approximation over the interval (1,2) are unpublished.
1463 !    Those for the approximation for 12 <= X are from reference 2.
1464 !
1465 !  Licensing:
1466 !
1467 !    This code is distributed under the GNU LGPL license.
1468 !
1469 !  Modified:
1470 !
1471 !    11 February 2008
1472 !
1473 !  Author:
1474 !
1475 !    Original FORTRAN77 version by William Cody, Laura Stoltz.
1476 !    FORTRAN90 version by John Burkardt.
1477 !
1478 !  Reference:
1479 !
1480 !    William Cody,
1481 !    An Overview of Software Development for Special Functions,
1482 !    in Numerical Analysis Dundee, 1975,
1483 !    edited by GA Watson,
1484 !    Lecture Notes in Mathematics 506,
1485 !    Springer, 1976.
1486 !
1487 !    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
1488 !    Charles Mesztenyi, John Rice, Henry Thatcher,
1489 !    Christoph Witzgall,
1490 !    Computer Approximations,
1491 !    Wiley, 1968,
1492 !    LC: QA297.C64.
1493 !
1494 !  Parameters:
1495 !
1496 !    Input, real ( dp ) X, the argument of the function.
1497 !
1498 !    Output, real ( dp ) R8_GAMMA, the value of the function.
1499 !
1500 
1501 !This section has been created automatically by the script Abilint (TD).
1502 !Do not modify the following lines by hand.
1503 #undef ABI_FUNC
1504 #define ABI_FUNC 'r8_gamma_gq'
1505 !End of the abilint section
1506 
1507   implicit none
1508 
1509   real ( dp ), dimension ( 7 ) :: c = (/ &
1510    -1.910444077728D-03, &
1511     8.4171387781295D-04, &
1512    -5.952379913043012D-04, &
1513     7.93650793500350248D-04, &
1514    -2.777777777777681622553D-03, &
1515     8.333333333333333331554247D-02, &
1516     5.7083835261D-03 /)
1517   real ( dp ), parameter :: eps = 2.22D-16
1518   real ( dp ) fact
1519   integer ( kind = 4 ) i
1520   integer ( kind = 4 ) n
1521   real ( dp ), dimension ( 8 ) :: p = (/ &
1522     -1.71618513886549492533811D+00, &
1523      2.47656508055759199108314D+01, &
1524     -3.79804256470945635097577D+02, &
1525      6.29331155312818442661052D+02, &
1526      8.66966202790413211295064D+02, &
1527     -3.14512729688483675254357D+04, &
1528     -3.61444134186911729807069D+04, &
1529      6.64561438202405440627855D+04 /)
1530   logical parity
1531   real ( dp ), parameter :: pi = 3.1415926535897932384626434D+00
1532   real ( dp ), dimension ( 8 ) :: q = (/ &
1533     -3.08402300119738975254353D+01, &
1534      3.15350626979604161529144D+02, &
1535     -1.01515636749021914166146D+03, &
1536     -3.10777167157231109440444D+03, &
1537      2.25381184209801510330112D+04, &
1538      4.75584627752788110767815D+03, &
1539     -1.34659959864969306392456D+05, &
1540     -1.15132259675553483497211D+05 /)
1541   real ( dp ) r8_gamma_gq
1542   real ( dp ) res
1543   real ( dp ), parameter :: sqrtpi = 0.9189385332046727417803297D+00
1544   real ( dp ) sum
1545   real ( dp ) x
1546   real ( dp ), parameter :: xbig = 171.624D+00
1547   real ( dp ) xden
1548   real ( dp ), parameter :: xinf = 1.0D+30
1549   real ( dp ), parameter :: xminin = 2.23D-308
1550   real ( dp ) xnum
1551   real ( dp ) y
1552   real ( dp ) y1
1553   real ( dp ) ysq
1554   real ( dp ) z
1555 
1556 ! *************************************************************************
1557 
1558   parity = .false.
1559   fact = 1.0D+00
1560   n = 0
1561   y = x
1562 !
1563 !  Argument is negative.
1564 !
1565   if ( y <= 0.0D+00 ) then
1566 
1567     y = - x
1568     y1 = aint ( y )
1569     res = y - y1
1570 
1571     if (abs(res) > tol20) then !if ( res /= 0.0D+00 ) then
1572 
1573       if ( abs(y1 - aint ( y1 * 0.5D+00 ) * 2.0D+00) > tol20 ) then ! if ( y1 /= aint ( y1 * 0.5D+00 ) * 2.0D+00 ) then
1574         parity = .true.
1575       end if
1576 
1577       fact = - pi / sin ( pi * res )
1578       y = y + 1.0D+00
1579 
1580     else
1581 
1582       res = xinf
1583       r8_gamma_gq = res
1584       return
1585 
1586     end if
1587 
1588   end if
1589 !
1590 !  Argument is positive.
1591 !
1592   if ( y < eps ) then
1593 !
1594 !  Argument < EPS.
1595 !
1596     if ( xminin <= y ) then
1597       res = 1.0D+00 / y
1598     else
1599       res = xinf
1600       r8_gamma_gq = res
1601       return
1602     end if
1603 
1604   else if ( y < 12.0D+00 ) then
1605 
1606     y1 = y
1607 !
1608 !  0.0 < argument < 1.0.
1609 !
1610     if ( y < 1.0D+00 ) then
1611 
1612       z = y
1613       y = y + 1.0D+00
1614 !
1615 !  1.0 < argument < 12.0.
1616 !  Reduce argument if necessary.
1617 !
1618     else
1619 
1620       n = int ( y ) - 1
1621       y = y - real ( n, dp )
1622       z = y - 1.0D+00
1623 
1624     end if
1625 !
1626 !  Evaluate approximation for 1.0 < argument < 2.0.
1627 !
1628     xnum = 0.0D+00
1629     xden = 1.0D+00
1630     do i = 1, 8
1631       xnum = ( xnum + p(i) ) * z
1632       xden = xden * z + q(i)
1633     end do
1634 
1635     res = xnum / xden + 1.0D+00
1636 !
1637 !  Adjust result for case  0.0 < argument < 1.0.
1638 !
1639     if ( y1 < y ) then
1640 
1641       res = res / y1
1642 !
1643 !  Adjust result for case 2.0 < argument < 12.0.
1644 !
1645     else if ( y < y1 ) then
1646 
1647       do i = 1, n
1648         res = res * y
1649         y = y + 1.0D+00
1650       end do
1651 
1652     end if
1653 
1654   else
1655 !
1656 !  Evaluate for 12.0 <= argument.
1657 !
1658     if ( y <= xbig ) then
1659 
1660       ysq = y * y
1661       sum = c(7)
1662       do i = 1, 6
1663         sum = sum / ysq + c(i)
1664       end do
1665       sum = sum / y - y + sqrtpi
1666       sum = sum + ( y - 0.5D+00 ) * log ( y )
1667       res = exp ( sum )
1668 
1669     else
1670 
1671       res = xinf
1672       r8_gamma_gq = res
1673       return
1674 
1675     end if
1676 
1677   end if
1678 !
1679 !  Final adjustments and return.
1680 !
1681   if ( parity ) then
1682     res = - res
1683   end if
1684 
1685   if ( abs(fact - 1.0D+00) > tol20 ) then !if ( fact /= 1.0D+00 ) then
1686     res = fact / res
1687   end if
1688 
1689   r8_gamma_gq = res
1690 
1691   return
1692 end function r8_gamma_gq

m_hamiltonian/r8mat_write [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  r8mat_write

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

1714 subroutine r8mat_write ( output_filename, m, n, table )
1715 
1716 !*****************************************************************************80
1717 !
1718 !! R8MAT_WRITE writes an R8MAT file.
1719 !
1720 !  Licensing:
1721 !
1722 !    This code is distributed under the GNU LGPL license.
1723 !
1724 !  Modified:
1725 !
1726 !    31 May 2009
1727 !
1728 !  Author:
1729 !
1730 !    John Burkardt
1731 !
1732 !  Parameters:
1733 !
1734 !    Input, character ( len = * ) OUTPUT_FILENAME, the output file name.
1735 !
1736 !    Input, integer ( kind = 4 ) M, the spatial dimension.
1737 !
1738 !    Input, integer ( kind = 4 ) N, the number of points.
1739 !
1740 !    Input, real ( dp ) TABLE(M,N), the table data.
1741 !
1742 
1743 !This section has been created automatically by the script Abilint (TD).
1744 !Do not modify the following lines by hand.
1745 #undef ABI_FUNC
1746 #define ABI_FUNC 'r8mat_write'
1747 !End of the abilint section
1748 
1749   implicit none
1750 
1751   integer ( kind = 4 ) m
1752   integer ( kind = 4 ) n
1753 
1754   integer ( kind = 4 ) j
1755   character ( len = * ) output_filename
1756   character(len=500) :: msg
1757   integer ( kind = 4 ) output_unit
1758   character ( len = 30 ) string
1759   real ( dp ) table(m,n)
1760 
1761 ! *************************************************************************
1762 
1763 !
1764 !  Open the file.
1765 !
1766   if (open_file(output_filename, msg, newunit=output_unit, status='replace') /= 0) then
1767     MSG_ERROR(msg)
1768   end if
1769 
1770 !
1771 !  Create a format string.
1772 !
1773 !  For less precision in the output file, try:
1774 !
1775 !                                            '(', m, 'g', 14, '.', 6, ')'
1776 !
1777   if ( 0 < m .and. 0 < n ) then
1778 
1779     write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')'
1780 !
1781 !  Write the data.
1782 !
1783     do j = 1, n
1784       write ( output_unit, string ) table(1:m,j)
1785     end do
1786 
1787   end if
1788 !
1789 !  Close the file.
1790 !
1791   close ( unit = output_unit )
1792 
1793 end subroutine r8mat_write

m_hamiltonian/rule_write [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  rule_write

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

      imtqlx

SOURCE

1815 subroutine rule_write ( order, x, w, r, filename )
1816 
1817 !*****************************************************************************80
1818 !
1819 !! RULE_WRITE writes a quadrature rule to a file.
1820 !
1821 !  Licensing:
1822 !
1823 !    This code is distributed under the GNU LGPL license.
1824 !
1825 !  Modified:
1826 !
1827 !    18 February 2010
1828 !
1829 !  Author:
1830 !
1831 !    John Burkardt
1832 !
1833 !  Parameters:
1834 !
1835 !    Input, integer ( kind = 4 ) ORDER, the order of the rule.
1836 !
1837 !    Input, real ( dp ) X(ORDER), the abscissas.
1838 !
1839 !    Input, real ( dp ) W(ORDER), the weights.
1840 !
1841 !    Input, real ( dp ) R(2), defines the region.
1842 !
1843 !    Input, character ( len = * ) FILENAME, specifies the output.
1844 !    'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining weights,
1845 !    abscissas, and region.
1846 !
1847 
1848 !This section has been created automatically by the script Abilint (TD).
1849 !Do not modify the following lines by hand.
1850 #undef ABI_FUNC
1851 #define ABI_FUNC 'rule_write'
1852 !End of the abilint section
1853 
1854   implicit none
1855 
1856   integer ( kind = 4 ) order
1857 
1858   character ( len = * ) filename
1859   character ( len = 255 ) filename_r
1860   character ( len = 255 ) filename_w
1861   character ( len = 255 ) filename_x
1862   real ( dp ) r(2)
1863   real ( dp ) w(order)
1864   real ( dp ) x(order)
1865 
1866 ! *************************************************************************
1867 
1868   filename_w = trim ( filename ) // '_w.txt'
1869   filename_x = trim ( filename ) // '_x.txt'
1870   filename_r = trim ( filename ) // '_r.txt'
1871 
1872   write ( std_out, '(a)' ) ' '
1873   write ( std_out, '(a)' ) '  Creating quadrature files.'
1874   write ( std_out, '(a)' ) ' '
1875   write ( std_out, '(a)' ) '  "Root" file name is   "' // trim ( filename ) // '".'
1876   write ( std_out, '(a)' ) ' '
1877   write ( std_out, '(a)' ) '  Weight file will be   "' // trim ( filename_w ) // '".'
1878   write ( std_out, '(a)' ) '  Abscissa file will be "' // trim ( filename_x ) // '".'
1879   write ( std_out, '(a)' ) '  Region file will be   "' // trim ( filename_r ) // '".'
1880 
1881   call r8mat_write ( filename_w, 1, order, w )
1882   call r8mat_write ( filename_x, 1, order, x )
1883   call r8mat_write ( filename_r, 1, 2,     r )
1884 
1885   return
1886 end subroutine rule_write

m_hamiltonian/s_to_i4 [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  s_to_i4

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

      imtqlx

SOURCE

1908 subroutine s_to_i4 ( s, ival, ierror, length )
1909 
1910 !*****************************************************************************80
1911 !
1912 !! S_TO_I4 reads an I4 from a string.
1913 !
1914 !  Licensing:
1915 !
1916 !    This code is distributed under the GNU LGPL license.
1917 !
1918 !  Modified:
1919 !
1920 !    28 June 2000
1921 !
1922 !  Author:
1923 !
1924 !    John Burkardt
1925 !
1926 !  Parameters:
1927 !
1928 !    Input, character ( len = * ) S, a string to be examined.
1929 !
1930 !    Output, integer ( kind = 4 ) IVAL, the integer value read from the string.
1931 !    If the string is blank, then IVAL will be returned 0.
1932 !
1933 !    Output, integer ( kind = 4 ) IERROR, an error flag.
1934 !    0, no error.
1935 !    1, an error occurred.
1936 !
1937 !    Output, integer ( kind = 4 ) LENGTH, the number of characters of S
1938 !    used to make IVAL.
1939 !
1940 
1941 !This section has been created automatically by the script Abilint (TD).
1942 !Do not modify the following lines by hand.
1943 #undef ABI_FUNC
1944 #define ABI_FUNC 's_to_i4'
1945 !End of the abilint section
1946 
1947   implicit none
1948 
1949   character c
1950   integer ( kind = 4 ) i
1951   integer ( kind = 4 ) ierror
1952   integer ( kind = 4 ) isgn
1953   integer ( kind = 4 ) istate
1954   integer ( kind = 4 ) ival
1955   integer ( kind = 4 ) length
1956   character ( len = * ) s
1957 
1958 ! *************************************************************************
1959 
1960   ierror = 0
1961   istate = 0
1962   isgn = 1
1963   ival = 0
1964 
1965   do i = 1, len_trim ( s )
1966 
1967     c = s(i:i)
1968 !
1969 !  Haven't read anything.
1970 !
1971     if ( istate == 0 ) then
1972 
1973       if ( c == ' ' ) then
1974 
1975       else if ( c == '-' ) then
1976         istate = 1
1977         isgn = -1
1978       else if ( c == '+' ) then
1979         istate = 1
1980         isgn = + 1
1981       else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
1982         istate = 2
1983         ival = ichar ( c ) - ichar ( '0' )
1984       else
1985         ierror = 1
1986         return
1987       end if
1988 !
1989 !  Have read the sign, expecting digits.
1990 !
1991     else if ( istate == 1 ) then
1992 
1993       if ( c == ' ' ) then
1994 
1995       else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
1996         istate = 2
1997         ival = ichar ( c ) - ichar ( '0' )
1998       else
1999         ierror = 1
2000         return
2001       end if
2002 !
2003 !  Have read at least one digit, expecting more.
2004 !
2005     else if ( istate == 2 ) then
2006 
2007       if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
2008         ival = 10 * ival + ichar ( c ) - ichar ( '0' )
2009       else
2010         ival = isgn * ival
2011         length = i - 1
2012         return
2013       end if
2014 
2015     end if
2016 
2017   end do
2018 !
2019 !  If we read all the characters in the string, see if we're OK.
2020 !
2021   if ( istate == 2 ) then
2022     ival = isgn * ival
2023     length = len_trim ( s )
2024   else
2025     ierror = 1
2026     length = 0
2027   end if
2028 
2029   return
2030 end subroutine s_to_i4

m_hamiltonian/s_to_r8 [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  s_to_r8

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

      imtqlx

SOURCE

2051 subroutine s_to_r8 ( s, dval, ierror, length )
2052 
2053 !*****************************************************************************80
2054 !
2055 !! S_TO_R8 reads an R8 from a string.
2056 !
2057 !  Discussion:
2058 !
2059 !    The routine will read as many characters as possible until it reaches
2060 !    the end of the string, or encounters a character which cannot be
2061 !    part of the number.
2062 !
2063 !    Legal input is:
2064 !
2065 !       1 blanks,
2066 !       2 '+' or '-' sign,
2067 !       2.5 blanks
2068 !       3 integer part,
2069 !       4 decimal point,
2070 !       5 fraction part,
2071 !       6 'E' or 'e' or 'D' or 'd', exponent marker,
2072 !       7 exponent sign,
2073 !       8 exponent integer part,
2074 !       9 exponent decimal point,
2075 !      10 exponent fraction part,
2076 !      11 blanks,
2077 !      12 final comma or semicolon,
2078 !
2079 !    with most quantities optional.
2080 !
2081 !  Example:
2082 !
2083 !    S                 DVAL
2084 !
2085 !    '1'               1.0
2086 !    '     1   '       1.0
2087 !    '1A'              1.0
2088 !    '12,34,56'        12.0
2089 !    '  34 7'          34.0
2090 !    '-1E2ABCD'        -100.0
2091 !    '-1X2ABCD'        -1.0
2092 !    ' 2E-1'           0.2
2093 !    '23.45'           23.45
2094 !    '-4.2E+2'         -420.0
2095 !    '17d2'            1700.0
2096 !    '-14e-2'         -0.14
2097 !    'e2'              100.0
2098 !    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
2099 !
2100 !  Licensing:
2101 !
2102 !    This code is distributed under the GNU LGPL license.
2103 !
2104 !  Modified:
2105 !
2106 !    07 September 2004
2107 !
2108 !  Author:
2109 !
2110 !    John Burkardt
2111 !
2112 !  Parameters:
2113 !
2114 !    Input, character ( len = * ) S, the string containing the
2115 !    data to be read.  Reading will begin at position 1 and
2116 !    terminate at the end of the string, or when no more
2117 !    characters can be read to form a legal real.  Blanks,
2118 !    commas, or other nonnumeric data will, in particular,
2119 !    cause the conversion to halt.
2120 !
2121 !    Output, real ( dp ) DVAL, the value read from the string.
2122 !
2123 !    Output, integer ( kind = 4 ) IERROR, error flag.
2124 !    0, no errors occurred.
2125 !    1, 2, 6 or 7, the input number was garbled.  The
2126 !    value of IERROR is the last type of input successfully
2127 !    read.  For instance, 1 means initial blanks, 2 means
2128 !    a plus or minus sign, and so on.
2129 !
2130 !    Output, integer ( kind = 4 ) LENGTH, the number of characters read
2131 !    to form the number, including any terminating
2132 !    characters such as a trailing comma or blanks.
2133 !
2134 
2135 !This section has been created automatically by the script Abilint (TD).
2136 !Do not modify the following lines by hand.
2137 #undef ABI_FUNC
2138 #define ABI_FUNC 's_to_r8'
2139 !End of the abilint section
2140 
2141   implicit none
2142 
2143   character c
2144   real ( dp ) dval
2145   integer ( kind = 4 ) ierror
2146   integer ( kind = 4 ) ihave
2147   integer ( kind = 4 ) isgn
2148   integer ( kind = 4 ) iterm
2149   integer ( kind = 4 ) jbot
2150   integer ( kind = 4 ) jsgn
2151   integer ( kind = 4 ) jtop
2152   integer ( kind = 4 ) length
2153   integer ( kind = 4 ) nchar
2154   integer ( kind = 4 ) ndig
2155   real ( dp ) rbot
2156   real ( dp ) rexp
2157   real ( dp ) rtop
2158   character ( len = * ) s
2159 
2160 ! *************************************************************************
2161 
2162   nchar = len_trim ( s )
2163 
2164   ierror = 0
2165   dval = 0.0D+00
2166   length = -1
2167   isgn = 1
2168   rtop = 0
2169   rbot = 1
2170   jsgn = 1
2171   jtop = 0
2172   jbot = 1
2173   ihave = 1
2174   iterm = 0
2175 
2176   do
2177 
2178     length = length + 1
2179 
2180     if ( nchar < length+1 ) then
2181       exit
2182     end if
2183 
2184     c = s(length+1:length+1)
2185 !
2186 !  Blank character.
2187 !
2188     if ( c == ' ' ) then
2189 
2190       if ( ihave == 2 ) then
2191 
2192       else if ( ihave == 6 .or. ihave == 7 ) then
2193         iterm = 1
2194       else if ( 1 < ihave ) then
2195         ihave = 11
2196       end if
2197 !
2198 !  Comma.
2199 !
2200     else if ( c == ',' .or. c == ';' ) then
2201 
2202       if ( ihave /= 1 ) then
2203         iterm = 1
2204         ihave = 12
2205         length = length + 1
2206       end if
2207 !
2208 !  Minus sign.
2209 !
2210     else if ( c == '-' ) then
2211 
2212       if ( ihave == 1 ) then
2213         ihave = 2
2214         isgn = -1
2215       else if ( ihave == 6 ) then
2216         ihave = 7
2217         jsgn = -1
2218       else
2219         iterm = 1
2220       end if
2221 !
2222 !  Plus sign.
2223 !
2224     else if ( c == '+' ) then
2225 
2226       if ( ihave == 1 ) then
2227         ihave = 2
2228       else if ( ihave == 6 ) then
2229         ihave = 7
2230       else
2231         iterm = 1
2232       end if
2233 !
2234 !  Decimal point.
2235 !
2236     else if ( c == '.' ) then
2237 
2238       if ( ihave < 4 ) then
2239         ihave = 4
2240       else if ( 6 <= ihave .and. ihave <= 8 ) then
2241         ihave = 9
2242       else
2243         iterm = 1
2244       end if
2245 !
2246 !  Scientific notation exponent marker.
2247 !
2248     else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then
2249 
2250       if ( ihave < 6 ) then
2251         ihave = 6
2252       else
2253         iterm = 1
2254       end if
2255 !
2256 !  Digit.
2257 !
2258     else if (  ihave < 11 .and. lle ( '0', c ) .and. lle ( c, '9' ) ) then
2259 
2260       if ( ihave <= 2 ) then
2261         ihave = 3
2262       else if ( ihave == 4 ) then
2263         ihave = 5
2264       else if ( ihave == 6 .or. ihave == 7 ) then
2265         ihave = 8
2266       else if ( ihave == 9 ) then
2267         ihave = 10
2268       end if
2269 
2270       call ch_to_digit ( c, ndig )
2271 
2272       if ( ihave == 3 ) then
2273         rtop = 10.0D+00 * rtop + real ( ndig, dp )
2274       else if ( ihave == 5 ) then
2275         rtop = 10.0D+00 * rtop + real ( ndig, dp )
2276         rbot = 10.0D+00 * rbot
2277       else if ( ihave == 8 ) then
2278         jtop = 10 * jtop + ndig
2279       else if ( ihave == 10 ) then
2280         jtop = 10 * jtop + ndig
2281         jbot = 10 * jbot
2282       end if
2283 !
2284 !  Anything else is regarded as a terminator.
2285 !
2286     else
2287       iterm = 1
2288     end if
2289 !
2290 !  If we haven't seen a terminator, and we haven't examined the
2291 !  entire string, go get the next character.
2292 !
2293     if ( iterm == 1 ) then
2294       exit
2295     end if
2296 
2297   end do
2298 !
2299 !  If we haven't seen a terminator, and we have examined the
2300 !  entire string, then we're done, and LENGTH is equal to NCHAR.
2301 !
2302   if ( iterm /= 1 .and. length+1 == nchar ) then
2303     length = nchar
2304   end if
2305 !
2306 !  Number seems to have terminated.  Have we got a legal number?
2307 !  Not if we terminated in states 1, 2, 6 or 7!
2308 !
2309   if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then
2310     ierror = ihave
2311     write ( std_out, '(a)' ) ' '
2312     write ( std_out, '(a)' ) 'S_TO_R8 - Serious error!'
2313     write ( std_out, '(a)' ) '  Illegal or nonnumeric input:'
2314     write ( std_out, '(a)' ) '    ' // trim ( s )
2315     return
2316   end if
2317 !
2318 !  Number seems OK.  Form it.
2319 !
2320   if ( jtop == 0 ) then
2321     rexp = 1.0D+00
2322   else
2323     if ( jbot == 1 ) then
2324       rexp = 10.0D+00 ** ( jsgn * jtop )
2325     else
2326       rexp = 10.0D+00 ** ( real ( jsgn * jtop, dp ) &
2327         / real ( jbot, dp ) )
2328     end if
2329   end if
2330 
2331   dval = real ( isgn, dp ) * rexp * rtop / rbot
2332 
2333   return
2334 end subroutine s_to_r8

m_hamiltonian/scqf [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  scqf

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

2356 subroutine scqf ( nt, t, mlt, wts, nwts, ndx, swts, st, gaussian_kind, alpha, beta, a, &
2357   b )
2358 
2359 !*****************************************************************************80
2360 !
2361 !! SCQF scales a quadrature formula to a nonstandard interval.
2362 !
2363 !  Discussion:
2364 !
2365 !    The arrays WTS and SWTS may coincide.
2366 !
2367 !    The arrays T and ST may coincide.
2368 !
2369 !  Licensing:
2370 !
2371 !    This code is distributed under the GNU LGPL license.
2372 !
2373 !  Modified:
2374 !
2375 !    27 December 2009
2376 !
2377 !  Author:
2378 !
2379 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
2380 !    FORTRAN90 version by John Burkardt.
2381 !
2382 !  Reference:
2383 !
2384 !    Sylvan Elhay, Jaroslav Kautsky,
2385 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
2386 !    Interpolatory Quadrature,
2387 !    ACM Transactions on Mathematical Software,
2388 !    Volume 13, Number 4, December 1987, pages 399-415.
2389 !
2390 !  Parameters:
2391 !
2392 !    Input, integer ( kind = 4 ) NT, the number of knots.
2393 !
2394 !    Input, real ( dp ) T(NT), the original knots.
2395 !
2396 !    Input, integer ( kind = 4 ) MLT(NT), the multiplicity of the knots.
2397 !
2398 !    Input, real ( dp ) WTS(NWTS), the weights.
2399 !
2400 !    Input, integer ( kind = 4 ) NWTS, the number of weights.
2401 !
2402 !    Input, integer ( kind = 4 ) NDX(NT), used to index the array WTS.
2403 !    For more details see the comments in CAWIQ.
2404 !
2405 !    Output, real ( dp ) SWTS(NWTS), the scaled weights.
2406 !
2407 !    Output, real ( dp ) ST(NT), the scaled knots.
2408 !
2409 !    Input, integer ( kind = 4 ) KIND, the rule.
2410 !    1, Legendre,             (a,b)       1.0
2411 !    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
2412 !    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
2413 !    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
2414 !    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
2415 !    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
2416 !    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
2417 !    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
2418 !    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
2419 !
2420 !    Input, real ( dp ) ALPHA, the value of Alpha, if needed.
2421 !
2422 !    Input, real ( dp ) BETA, the value of Beta, if needed.
2423 !
2424 !    Input, real ( dp ) A, B, the interval endpoints.
2425 !
2426 
2427 !This section has been created automatically by the script Abilint (TD).
2428 !Do not modify the following lines by hand.
2429 #undef ABI_FUNC
2430 #define ABI_FUNC 'scqf'
2431 !End of the abilint section
2432 
2433   implicit none
2434 
2435   integer ( kind = 4 ) nt
2436   integer ( kind = 4 ) nwts
2437 
2438   real ( dp ) a
2439   real ( dp ) al
2440   real ( dp ) alpha
2441   real ( dp ) b
2442   real ( dp ) be
2443   real ( dp ) beta
2444   integer ( kind = 4 ) i
2445   integer ( kind = 4 ) k
2446   integer ( kind = 4 ) gaussian_kind
2447   integer ( kind = 4 ) l
2448   integer ( kind = 4 ) mlt(nt)
2449   integer ( kind = 4 ) ndx(nt)
2450   real ( dp ) p
2451   real ( dp ) shft
2452   real ( dp ) slp
2453   real ( dp ) st(nt)
2454   real ( dp ) swts(nwts)
2455   real ( dp ) t(nt)
2456   real ( dp ) temp
2457   real ( dp ) tmp
2458   real ( dp ) wts(nwts)
2459 
2460 ! *************************************************************************
2461 
2462   temp = epsilon ( temp )
2463 
2464   call parchk ( gaussian_kind, 1, alpha, beta )
2465 
2466   if ( gaussian_kind == 1 ) then
2467 
2468     al = 0.0D+00
2469     be = 0.0D+00
2470 
2471     if ( abs ( b - a ) <= temp ) then
2472       MSG_ERROR('SCQF - Fatal error! |B - A| too small.')
2473     end if
2474 
2475     shft = ( a + b ) / 2.0D+00
2476     slp = ( b - a ) / 2.0D+00
2477 
2478   else if ( gaussian_kind == 2 ) then
2479 
2480     al = -0.5D+00
2481     be = -0.5D+00
2482 
2483     if ( abs ( b - a ) <= temp ) then
2484       MSG_ERROR('SCQF - Fatal error! |B - A| too small.')
2485     end if
2486 
2487     shft = ( a + b ) / 2.0D+00
2488     slp = ( b - a ) / 2.0D+00
2489 
2490   else if ( gaussian_kind == 3 ) then
2491 
2492     al = alpha
2493     be = alpha
2494 
2495     if ( abs ( b - a ) <= temp ) then
2496       MSG_ERROR('SCQF - Fatal error!  |B - A| too small.')
2497     end if
2498 
2499     shft = ( a + b ) / 2.0D+00
2500     slp = ( b - a ) / 2.0D+00
2501 
2502   else if ( gaussian_kind == 4 ) then
2503 
2504     al = alpha
2505     be = beta
2506 
2507     if ( abs ( b - a ) <= temp ) then
2508       MSG_ERROR("SCQF - Fatal error! |B - A| too small.")
2509     end if
2510 
2511     shft = ( a + b ) / 2.0D+00
2512     slp = ( b - a ) / 2.0D+00
2513 
2514   else if ( gaussian_kind == 5 ) then
2515 
2516     if ( b <= 0.0D+00 ) then
2517       MSG_ERROR('SCQF - Fatal error!  B <= 0')
2518     end if
2519 
2520     shft = a
2521     slp = 1.0D+00 / b
2522     al = alpha
2523     be = 0.0D+00
2524 
2525   else if ( gaussian_kind == 6 ) then
2526 
2527     if ( b <= 0.0D+00 ) then
2528       MSG_ERROR('SCQF - Fatal error! B <= 0.')
2529     end if
2530 
2531     shft = a
2532     slp = 1.0D+00 / sqrt ( b )
2533     al = alpha
2534     be = 0.0D+00
2535 
2536   else if ( gaussian_kind == 7 ) then
2537 
2538     al = alpha
2539     be = 0.0D+00
2540 
2541     if ( abs ( b - a ) <= temp ) then
2542       MSG_ERROR('|B - A| too small.')
2543     end if
2544 
2545     shft = ( a + b ) / 2.0D+00
2546     slp = ( b - a ) / 2.0D+00
2547 
2548   else if ( gaussian_kind == 8 ) then
2549 
2550     if ( a + b <= 0.0D+00 ) then
2551       MSG_ERROR('  A + B <= 0.')
2552     end if
2553 
2554     shft = a
2555     slp = a + b
2556     al = alpha
2557     be = beta
2558 
2559   else if ( gaussian_kind == 9 ) then
2560 
2561     al = 0.5D+00
2562     be = 0.5D+00
2563 
2564     if ( abs ( b - a ) <= temp ) then
2565       MSG_ERROR('|B - A| too small.')
2566     end if
2567 
2568     shft = ( a + b ) / 2.0D+00
2569     slp = ( b - a ) / 2.0D+00
2570 
2571   end if
2572 
2573   p = slp**( al + be + 1.0D+00 )
2574 
2575   do k = 1, nt
2576 
2577     st(k) = shft + slp * t(k)
2578     l = abs ( ndx(k) )
2579 
2580     if ( l /= 0 ) then
2581       tmp = p
2582       do i = l, l + mlt(k) - 1
2583         swts(i) = wts(i) * tmp
2584         tmp = tmp * slp
2585       end do
2586     end if
2587 
2588   end do
2589 
2590   return
2591 end subroutine scqf

m_hamiltonian/sgqf [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  sgqf

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gaussian_quadrature

CHILDREN

      imtqlx

SOURCE

2614 subroutine sgqf ( nt, aj, bj, zemu, t, wts )
2615 
2616 !*****************************************************************************80
2617 !
2618 !! SGQF computes knots and weights of a Gauss Quadrature formula.
2619 !
2620 !  Discussion:
2621 !
2622 !    This routine computes all the knots and weights of a Gauss quadrature
2623 !    formula with simple knots from the Jacobi matrix and the zero-th
2624 !    moment of the weight function, using the Golub-Welsch technique.
2625 !
2626 !  Licensing:
2627 !
2628 !    This code is distributed under the GNU LGPL license.
2629 !
2630 !  Modified:
2631 !
2632 !    04 January 2010
2633 !
2634 !  Author:
2635 !
2636 !    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
2637 !    FORTRAN90 version by John Burkardt.
2638 !
2639 !  Reference:
2640 !
2641 !    Sylvan Elhay, Jaroslav Kautsky,
2642 !    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
2643 !    Interpolatory Quadrature,
2644 !    ACM Transactions on Mathematical Software,
2645 !    Volume 13, Number 4, December 1987, pages 399-415.
2646 !
2647 !  Parameters:
2648 !
2649 !    Input, integer ( kind = 4 ) NT, the number of knots.
2650 !
2651 !    Input, real ( dp ) AJ(NT), the diagonal of the Jacobi matrix.
2652 !
2653 !    Input/output, real ( dp ) BJ(NT), the subdiagonal of the Jacobi
2654 !    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
2655 !
2656 !    Input, real ( dp ) ZEMU, the zero-th moment of the weight function.
2657 !
2658 !    Output, real ( dp ) T(NT), the knots.
2659 !
2660 !    Output, real ( dp ) WTS(NT), the weights.
2661 !
2662 
2663 !This section has been created automatically by the script Abilint (TD).
2664 !Do not modify the following lines by hand.
2665 #undef ABI_FUNC
2666 #define ABI_FUNC 'sgqf'
2667 !End of the abilint section
2668 
2669   implicit none
2670 
2671   integer ( kind = 4 ) nt
2672 
2673   real ( dp ) aj(nt)
2674   real ( dp ) bj(nt)
2675   real ( dp ) t(nt)
2676   real ( dp ) wts(nt)
2677   real ( dp ) zemu
2678 
2679 ! *************************************************************************
2680 
2681 !
2682 !  Exit if the zero-th moment is not positive.
2683 !
2684   if ( zemu <= 0.0D+00 ) then
2685     MSG_ERROR('SGQF - Fatal error! ZEMU <= 0.')
2686   end if
2687 !
2688 !  Set up vectors for IMTQLX.
2689 !
2690   t(1:nt) = aj(1:nt)
2691 
2692   wts(1) = sqrt ( zemu )
2693   wts(2:nt) = 0.0D+00
2694 !
2695 !  Diagonalize the Jacobi matrix.
2696 !
2697   call imtqlx ( nt, t, bj, wts )
2698 
2699   wts(1:nt) = wts(1:nt)**2
2700 
2701   return
2702 end subroutine sgqf