TABLE OF CONTENTS
- ABINIT/jacobi
- ABINIT/m_abilasi
- m_abilasi/cginv
- m_abilasi/dzegdi
- m_abilasi/dzgefa
- m_abilasi/lubksb
- m_abilasi/ludcmp
- m_abilasi/matr3eigval
- m_abilasi/matrginv
- m_abilasi/test_xginv
- m_abilasi/wrap_CGEEV
- m_abilasi/wrap_CHEEV
- m_abilasi/wrap_CHPEV
- m_abilasi/wrap_DSYEV_ZHEEV
- m_abilasi/wrap_DSYEVX_ZHEEVX
- m_abilasi/wrap_DSYGV_ZHEGV
- m_abilasi/wrap_DSYGVX_ZHEGVX
- m_abilasi/wrap_ZGEEV
- m_abilasi/wrap_ZHEEV
- m_abilasi/wrap_ZHEEVX
- m_abilasi/wrap_ZHEGV
- m_abilasi/wrap_ZHEGVX
- m_abilasi/wrap_ZHPEV
- m_abilasi/zginv
- m_abilasi/zhpd_invert
ABINIT/jacobi [ Functions ]
NAME
jacobi
FUNCTION
Computes all eigenvalues and eigenvectors of a real symmetric matrix a, which is of size n by n, stored in a physical np by np array. On output, elements of a above the diagonal are destroyed. d returns the eigenvalues of a in its first n elements. v is a matrix with the same logical and physical dimensions as a, whose columns contain, on output, the normalized eigenvectors of a. nrot returns the number of Jacobi rotations that were required.
INPUTS
OUTPUT
NOTES
This routine is deprecated, use Lapack API
PARENTS
conducti_nc,critic
CHILDREN
SOURCE
3610 subroutine jacobi(a,n,np,d,v,nrot) 3611 3612 3613 !This section has been created automatically by the script Abilint (TD). 3614 !Do not modify the following lines by hand. 3615 #undef ABI_FUNC 3616 #define ABI_FUNC 'jacobi' 3617 !End of the abilint section 3618 3619 implicit none 3620 !Arguments 3621 integer :: n,np,nrot 3622 real*8 :: a(np,np),d(np),v(np,np) 3623 !Local variables 3624 integer, parameter :: NMAX=500 3625 integer i,ip,iq,j 3626 real*8 c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) 3627 do ip=1,n 3628 do iq=1,n 3629 v(ip,iq)=0. 3630 enddo 3631 v(ip,ip)=1. 3632 enddo 3633 do ip=1,n 3634 b(ip)=a(ip,ip) 3635 d(ip)=b(ip) 3636 z(ip)=0. 3637 enddo 3638 nrot=0 3639 do i=1,50 3640 sm=0. 3641 do ip=1,n-1 3642 do iq=ip+1,n 3643 sm=sm+abs(a(ip,iq)) 3644 enddo 3645 enddo 3646 if(sm.eq.0.)return 3647 if(i.lt.4)then 3648 tresh=0.2*sm/n**2 3649 else 3650 tresh=0. 3651 endif 3652 do ip=1,n-1 3653 do iq=ip+1,n 3654 g=100.*abs(a(ip,iq)) 3655 if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) & 3656 & .and.(abs(d(iq))+g.eq.abs(d(iq))))then 3657 a(ip,iq)=0. 3658 else if(abs(a(ip,iq)).gt.tresh)then 3659 h=d(iq)-d(ip) 3660 if(abs(h)+g.eq.abs(h))then 3661 t=a(ip,iq)/h 3662 else 3663 theta=0.5*h/a(ip,iq) 3664 t=1./(abs(theta)+sqrt(1.+theta**2)) 3665 if(theta.lt.0.)t=-t 3666 endif 3667 c=1./sqrt(1+t**2) 3668 s=t*c 3669 tau=s/(1.+c) 3670 h=t*a(ip,iq) 3671 z(ip)=z(ip)-h 3672 z(iq)=z(iq)+h 3673 d(ip)=d(ip)-h 3674 d(iq)=d(iq)+h 3675 a(ip,iq)=0. 3676 do j=1,ip-1 3677 g=a(j,ip) 3678 h=a(j,iq) 3679 a(j,ip)=g-s*(h+g*tau) 3680 a(j,iq)=h+s*(g-h*tau) 3681 enddo 3682 do j=ip+1,iq-1 3683 g=a(ip,j) 3684 h=a(j,iq) 3685 a(ip,j)=g-s*(h+g*tau) 3686 a(j,iq)=h+s*(g-h*tau) 3687 enddo 3688 do j=iq+1,n 3689 g=a(ip,j) 3690 h=a(iq,j) 3691 a(ip,j)=g-s*(h+g*tau) 3692 a(iq,j)=h+s*(g-h*tau) 3693 enddo 3694 do j=1,n 3695 g=v(j,ip) 3696 h=v(j,iq) 3697 v(j,ip)=g-s*(h+g*tau) 3698 v(j,iq)=h+s*(g-h*tau) 3699 enddo 3700 nrot=nrot+1 3701 endif 3702 enddo 3703 enddo 3704 do ip=1,n 3705 b(ip)=b(ip)+z(ip) 3706 d(ip)=b(ip) 3707 z(ip)=0. 3708 enddo 3709 enddo 3710 write(std_out,*) 'too many iterations in jacobi' 3711 3712 end subroutine jacobi
ABINIT/m_abilasi [ Modules ]
NAME
m_abilasi
FUNCTION
ABInit Linear Algebra Subroutine Interfaces. This modules provides interfaces performing the overloading of commonly used Lapack routines. The main purpose of this module is to create a layer between abinit routines and Lapack procedures. This layer can be used to hide the parallel Scalapack version. In this case, only the MPI commutator has to be provided in input as the wrapper will take care of the initialization of the Scalapack grid as well as of the distribution of the matrix. Note that this allows one to reduce the CPU time per processor but not the memory allocated since the entire matrix has to be provided in input. The interfaces are very similar to the Lapack F77 version (neither F90 constructs nor F90 assumed size arrays are used). The main simplification with respect to the F77 version of Lapack is that the work arrays are allocated inside the wrapper with optimal size thus reducing the number of input argcomm_scalapackuments that has to be passed. Leading dimensions have been removed from the interface whenever possible. In F90 one can pass the array descriptor if the routines should operate on a slice of the local array (seldom done in abinit). Using array descriptor is OK but will it likely slow-down the calculation as some compilers perform a copy of the input-output data. If efficiency is a concern, then the F77 call should be used
COPYRIGHT
Copyright (C) 1992-2018 ABINIT group (MG, GMR, XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
PARENTS
CHILDREN
TODO
1) Use a function to define the size of the Scalapack block according to some heuristic method. 2) Define a threshold below which Scalapack is not used although the MPI communicator is passed. 3) Split MPI communicator for Scalapack (.i.e. use a max size for the Scalapack comm; see abi_linalg_init). 4) On certain networks, xmpi_sum might crash due to the size of the MPI packet. This problem should be solved in hide_mpi (Module containing a private global variable defining a threshold above which the input array is split into smaller chunks.
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 MODULE m_abilasi 52 53 use defs_basis 54 use m_profiling_abi 55 use m_xmpi 56 use m_errors 57 use m_slk 58 use m_linalg_interfaces 59 60 use m_time, only : cwtime 61 use m_fstrings, only : firstchar 62 63 implicit none 64 65 private 66 67 ! Procedures for complex Hermitian matrices 68 69 public :: xheev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix. 70 71 public :: xhpev ! Computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix 72 ! in packed storage (Scalapack version not available) 73 74 public :: xhegv ! Compute all the eigenvalues, and optionally, the eigenvectors of a complex generalized 75 ! Hermitian-definite eigenproblem, of the form: 76 ! A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x 77 78 public :: xheevx ! Computes selected eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. 79 ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of 80 ! indices for the desired eigenvalues. 81 82 public :: xhegvx ! Computes selected eigenvalues, and optionally, eigenvectors of a complex generalized Hermitian-definite eigenproblem, 83 ! of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. 84 ! Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of 85 ! indices for the desired eigenvalues. 86 87 ! Procedures for complex non-symmetric matrices 88 89 public :: xgeev ! Computes for a complex nonsymmetric matrix A, the eigenvalues and, optionally, 90 ! the left and/or right eigenvectors. 91 92 public :: xginv ! Invert a general matrix of complex elements by means of LU factorization. 93 94 95 public :: xhdp_invert ! Invert a Hermitian positive definite matrix. 96 97 interface xheev 98 module procedure wrap_CHEEV 99 module procedure wrap_ZHEEV 100 module procedure wrap_DSYEV_ZHEEV 101 end interface xheev 102 103 interface xhpev 104 module procedure wrap_CHPEV 105 module procedure wrap_ZHPEV 106 end interface xhpev 107 108 interface xhegv 109 module procedure wrap_ZHEGV 110 module procedure wrap_DSYGV_ZHEGV 111 end interface xhegv 112 113 interface xheevx 114 module procedure wrap_ZHEEVX 115 module procedure wrap_DSYEVX_ZHEEVX 116 end interface xheevx 117 118 interface xhegvx 119 module procedure wrap_ZHEGVX 120 module procedure wrap_DSYGVX_ZHEGVX 121 end interface xhegvx 122 123 interface xgeev 124 module procedure wrap_CGEEV 125 module procedure wrap_ZGEEV 126 end interface xgeev 127 128 interface xginv 129 module procedure cginv 130 module procedure zginv 131 end interface xginv 132 133 interface xhdp_invert 134 module procedure zhpd_invert 135 end interface xhdp_invert 136 137 public :: matrginv ! Invert a general matrix of real*8 elements. 138 public :: matr3eigval ! Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode. 139 140 !FIXME This procedures are deprecated, use lapack API 141 public :: jacobi ! Computes all eigenvalues and eigenvectors of a real symmetric matrix a, 142 public :: ludcmp 143 public :: lubksb 144 public :: dzgedi 145 public :: dzgefa 146 147 !---------------------------------------------------------------------- 148 ! support for unitary tests and profiling. 149 150 type,public :: latime_t 151 character(len=500) :: testname 152 integer :: msize 153 real(dp) :: ctime 154 real(dp) :: wtime 155 real(dp) :: max_abserr=-one 156 real(dp) :: gflops 157 end type latime_t 158 159 public :: test_xginv 160 161 !---------------------------------------------------------------------- 162 ! private variables 163 164 integer,private,parameter :: SLK_BLOCK_SIZE = 24 165 ! Default block size for Scalapack distribution. 166 ! As recommended by Intel MKL, a more sensible default than the previous value of 40 167 168 CONTAINS !=========================================================================================================================
m_abilasi/cginv [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
cginv
FUNCTION
Invert a general matrix of complex elements in single precision. CGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges. CGETRI computes the inverse of a matrix using the LU factorization computed by CGETRF.
INPUTS
n=size of complex matrix a
a=matrix of complex elements
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)= array of complex elements, input, inverted at output
TODO
Add Scalapack version, matrix_scalapack has to be modified by adding a single precision complex buffer.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
2828 subroutine cginv(a,n,comm) 2829 2830 2831 !This section has been created automatically by the script Abilint (TD). 2832 !Do not modify the following lines by hand. 2833 #undef ABI_FUNC 2834 #define ABI_FUNC 'cginv' 2835 use interfaces_14_hidewrite 2836 !End of the abilint section 2837 2838 implicit none 2839 2840 !Arguments ------------------------------------ 2841 !scalars 2842 integer,intent(in) :: n 2843 integer,optional,intent(in) :: comm 2844 !arrays 2845 complex(spc),intent(inout) :: a(n,n) 2846 2847 !Local variables------------------------------- 2848 !scalars 2849 integer :: lwork,info,nprocs 2850 logical :: use_scalapack 2851 character(len=500) :: msg 2852 !arrays 2853 integer,allocatable :: ipiv(:) 2854 complex(spc),allocatable :: work(:) 2855 #ifdef HAVE_LINALG_SCALAPACK 2856 integer :: ierr,istwf_k,ipiv_size,liwork,tbloc 2857 integer,allocatable :: iwork(:) 2858 type(matrix_scalapack) :: Slk_mat 2859 type(processor_scalapack) :: Slk_processor 2860 #endif 2861 2862 ! ************************************************************************* 2863 2864 use_scalapack=.FALSE. 2865 if (PRESENT(comm)) then 2866 nprocs = xmpi_comm_size(comm) 2867 #ifdef HAVE_LINALG_SCALAPACK 2868 use_scalapack = (nprocs>1) 2869 #endif 2870 end if 2871 2872 SELECT CASE(use_scalapack) 2873 2874 CASE (.FALSE.) 2875 2876 ABI_MALLOC(ipiv,(n)) 2877 2878 call CGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization. 2879 2880 if (info < 0) then 2881 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRF had an illegal value." 2882 MSG_ERROR(msg) 2883 end if 2884 2885 if (info > 0) then 2886 write(msg,'(3a,i0,4a)')& 2887 & "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,& 2888 & "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,& 2889 & "The factorization has been completed but the factor U is exactly singular.",ch10,& 2890 & "Division by zero will occur if it is used to solve a system of equations." 2891 MSG_ERROR(msg) 2892 end if 2893 2894 lwork=MAX(1,n) 2895 2896 ABI_MALLOC(work,(lwork)) 2897 2898 call CGETRI(n,a,n,ipiv,work,lwork,info) ! Inverts U and the computes inv(A) 2899 2900 if (info < 0) then 2901 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRI had an illegal value." 2902 MSG_ERROR(msg) 2903 end if 2904 2905 if (info > 0) then 2906 write(msg,'(3a,i0,a)')& 2907 & "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,& 2908 "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed." 2909 MSG_ERROR(msg) 2910 end if 2911 2912 ABI_FREE(ipiv) 2913 ABI_FREE(work) 2914 2915 RETURN 2916 2917 CASE (.TRUE.) 2918 2919 #if 0 2920 ! FIXME matrix_scalapack does not have a single precision complex buffer 2921 2922 #ifdef HAVE_LINALG_SCALAPACK 2923 call init_scalapack(Slk_processor,comm) 2924 istwf_k=1 2925 2926 ! Initialize and fill Scalapack matrix from the global one. 2927 tbloc=SLK_BLOCK_SIZE 2928 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2929 2930 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 2931 call wrtout(std_out,msg,"COLL") 2932 2933 ! IMPORTANT NOTE: PZGETRF requires square block decomposition i.e., MB_A = NB_A. 2934 if ( Slk_mat%descript%tab(MB_)/=Slk_mat%descript%tab(NB_) ) then 2935 msg ="PZGETRF requires square block decomposition i.e., MB_A = NB_A." 2936 MSG_ERROR(msg) 2937 end if 2938 2939 !!call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a) 2940 2941 ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_) 2942 ABI_MALLOC(ipiv,(ipiv_size)) 2943 2944 call PCGETRF(Slk_mat%sizeb_global(1),Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx_sp,& 2945 & 1,1,Slk_mat%descript%tab,ipiv,info) ! P * L * U Factorization. 2946 2947 if (info/=0) then 2948 write(msg,'(a,i0)')"PCGETRF returned info= ",info 2949 MSG_ERROR(msg) 2950 end if 2951 2952 ! Get optimal size of workspace for PCGETRI. 2953 lwork=-1; liwork=-1 2954 ABI_MALLOC(work,(1)) 2955 ABI_MALLOC(iwork,(1)) 2956 2957 call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,& 2958 & work,lwork,iwork,liwork,info) 2959 2960 ABI_CHECK(info==0,"PZGETRI: Error during compuation of workspace size") 2961 2962 lwork = NINT(DBLE(work(1))); liwork=iwork(1) 2963 ABI_FREE(work) 2964 ABI_FREE(iwork) 2965 2966 ! Solve the problem. 2967 ABI_MALLOC(work,(lwork)) 2968 ABI_MALLOC(iwork,(liwork)) 2969 2970 call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,& 2971 & work,lwork,iwork,liwork,info) 2972 2973 if (info/=0) then 2974 write(msg,'(a,i0)')"PZGETRI returned info= ",info 2975 MSG_ERROR(msg) 2976 end if 2977 2978 ABI_FREE(work) 2979 ABI_FREE(iwork) 2980 ABI_FREE(ipiv) 2981 2982 ! Reconstruct the global matrix from the distributed one. 2983 a = czero 2984 !! call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node. 2985 call destruction_matrix_scalapack(Slk_mat) 2986 2987 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 2988 call end_scalapack(Slk_processor) 2989 2990 RETURN 2991 #endif 2992 2993 #endif 2994 2995 MSG_BUG("You should not be here!") 2996 2997 END SELECT 2998 2999 end subroutine cginv
m_abilasi/dzegdi [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
dzgedi
FUNCTION
This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16 for the purpose of ABINIT
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
berryphase,dfptnl_mv,qmatrix,relaxpol,smatrix,uderiv
CHILDREN
SOURCE
3926 subroutine dzgedi(a,lda,n,ipvt,det,work,job) 3927 3928 3929 !This section has been created automatically by the script Abilint (TD). 3930 !Do not modify the following lines by hand. 3931 #undef ABI_FUNC 3932 #define ABI_FUNC 'dzgedi' 3933 !End of the abilint section 3934 3935 implicit none 3936 3937 integer :: lda,n,ipvt(n),job 3938 real*8 :: a(2,lda,n),det(2,2),work(2,n) 3939 ! 3940 ! zgedi computes the determinant and inverse of a matrix 3941 ! using the factors computed by zgeco or zgefa. 3942 ! 3943 ! on entry 3944 ! 3945 ! a complex*16(lda, n) 3946 ! the output from zgeco or zgefa. 3947 ! 3948 ! lda integer 3949 ! the leading dimension of the array a . 3950 ! 3951 ! n integer 3952 ! the order of the matrix a . 3953 ! 3954 ! ipvt integer(n) 3955 ! the pivot vector from zgeco or zgefa. 3956 ! 3957 ! work complex*16(n) 3958 ! work vector. contents destroyed. 3959 ! 3960 ! job integer 3961 ! = 11 both determinant and inverse. 3962 ! = 01 inverse only. 3963 ! = 10 determinant only. 3964 ! 3965 ! on return 3966 ! 3967 ! a inverse of original matrix if requested. 3968 ! otherwise unchanged. 3969 ! 3970 ! det complex*16(2) 3971 ! determinant of original matrix if requested. 3972 ! otherwise not referenced. 3973 ! determinant = det(1) * 10.0**det(2) 3974 ! with 1.0 .le. cabs1(det(1)) .lt. 10.0 3975 ! or det(1) .eq. 0.0 . 3976 ! 3977 ! error condition 3978 ! 3979 ! a division by zero will occur if the input factor contains 3980 ! a zero on the diagonal and the inverse is requested. 3981 ! it will not occur if the subroutines are called correctly 3982 ! and if zgeco has set rcond .gt. 0.0 or zgefa has set 3983 ! info .eq. 0 . 3984 ! 3985 ! linpack. this version dated 08/14/78 . 3986 ! cleve moler, university of new mexico, argonne national lab. 3987 ! 3988 ! subroutines and functions 3989 ! 3990 ! internal variables 3991 ! 3992 double precision :: r(2),rk(2),rkj(2) 3993 double precision :: ten,rinv2,rabs 3994 integer :: i,j,k,kb,kp1,l,nm1 3995 ! 3996 ! compute determinant 3997 ! 3998 if (job/10 .eq. 0) go to 70 3999 det(1,1) = 1.0d0; det(2,1) = 0.0d0 4000 det(1,2) = 0.0d0; det(2,2) = 0.0d0 4001 ten = 10.0d0 4002 do i = 1, n 4003 if (ipvt(i) .ne. i) then 4004 det(1,1) = -det(1,1) 4005 det(2,1) = -det(2,1) 4006 end if 4007 r(1)=det(1,1); r(2)=det(2,1) 4008 det(1,1) = r(1)*a(1,i,i)-r(2)*a(2,i,i) 4009 det(2,1) = r(2)*a(1,i,i)+r(1)*a(2,i,i) 4010 ! ...exit 4011 rabs = abs(det(1,1))+abs(det(2,1)) 4012 if (rabs .eq. 0.0d0) go to 60 4013 10 continue 4014 rabs = abs(det(1,1))+abs(det(2,1)) 4015 if (rabs .ge. 1.0d0) go to 20 4016 det(1,1) = ten*det(1,1); det(2,1) = ten*det(2,1) 4017 det(1,2) = det(1,2) - 1.0d0 4018 go to 10 4019 20 continue 4020 30 continue 4021 rabs = abs(det(1,1))+abs(det(2,1)) 4022 if (rabs .lt. ten) go to 40 4023 det(1,1) = det(1,1)/ten; det(2,1) = det(2,1)/ten 4024 det(1,2) = det(1,2) + 1.0d0 4025 go to 30 4026 40 continue 4027 end do 4028 60 continue 4029 70 continue 4030 ! 4031 ! compute inverse(u) 4032 ! 4033 if (mod(job,10) .eq. 0) go to 150 4034 do 100 k = 1, n 4035 !a(k,k) = (1.0d0,0.0d0)/a(k,k) 4036 !t = -a(k,k) 4037 !call zscal(k-1,t,a(1,k),1) 4038 rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2) 4039 a(1,k,k) = rinv2*a(1,k,k) 4040 a(2,k,k) = -rinv2*a(2,k,k) 4041 rk(1) = -a(1,k,k); rk(2) = -a(2,k,k) 4042 do i=1,k-1 4043 r(1)=a(1,i,k) 4044 r(2)=a(2,i,k) 4045 a(1,i,k)=rk(1)*r(1)-rk(2)*r(2) 4046 a(2,i,k)=rk(1)*r(2)+rk(2)*r(1) 4047 end do 4048 kp1 = k + 1 4049 if (n .lt. kp1) go to 90 4050 do 80 j = kp1, n 4051 !t = a(k,j) 4052 !a(k,j) = (0.0d0,0.0d0) 4053 !call zaxpy(k,t,a(1,k),1,a(1,j),1) 4054 rkj(1) = a(1,k,j); rkj(2) = a(2,k,j) 4055 a(1,k,j) = 0.d0; a(2,k,j) = 0.d0 4056 do i=1,k 4057 a(1,i,j)=rkj(1)*a(1,i,k)-rkj(2)*a(2,i,k)+a(1,i,j) 4058 a(2,i,j)=rkj(2)*a(1,i,k)+rkj(1)*a(2,i,k)+a(2,i,j) 4059 end do 4060 80 continue 4061 90 continue 4062 100 continue 4063 do i=1,n 4064 end do 4065 ! 4066 ! form inverse(u)*inverse(l) 4067 ! 4068 nm1 = n - 1 4069 if (nm1 .lt. 1) go to 140 4070 do 130 kb = 1, nm1 4071 k = n - kb 4072 kp1 = k + 1 4073 do 110 i = kp1, n 4074 work(1,i) = a(1,i,k); work(2,i) = a(2,i,k) 4075 a(1,i,k) = 0.0d0; a(2,i,k) = 0.d0 4076 110 continue 4077 do 120 j = kp1, n 4078 r(1) = work(1,j); r(2) = work(2,j) 4079 !call zaxpy(n,t,a(1,j),1,a(1,k),1) 4080 do i=1,n 4081 a(1,i,k)=r(1)*a(1,i,j)-r(2)*a(2,i,j)+a(1,i,k) 4082 a(2,i,k)=r(2)*a(1,i,j)+r(1)*a(2,i,j)+a(2,i,k) 4083 end do 4084 120 continue 4085 l = ipvt(k) 4086 if (l .ne. k) then 4087 !call zswap(n,a(1,k),1,a(1,l),1) 4088 do i=1,n 4089 r(1) = a(1,i,k); r(2) = a(2,i,k) 4090 a(1,i,k) = a(1,i,l); a(2,i,k) = a(2,i,l) 4091 a(1,i,l) = r(1); a(2,i,l) = r(2) 4092 end do 4093 end if 4094 130 continue 4095 140 continue 4096 150 continue 4097 4098 end subroutine dzgedi
m_abilasi/dzgefa [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
dzgefa
FUNCTION
This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16 for the purpose of ABINIT (2008,TD)
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
berryphase,dfptnl_mv,qmatrix,relaxpol,smatrix,uderiv
CHILDREN
SOURCE
4124 subroutine dzgefa(a,lda,n,ipvt,info) 4125 4126 use m_linalg_interfaces 4127 4128 !This section has been created automatically by the script Abilint (TD). 4129 !Do not modify the following lines by hand. 4130 #undef ABI_FUNC 4131 #define ABI_FUNC 'dzgefa' 4132 !End of the abilint section 4133 4134 implicit none 4135 4136 !Arguments 4137 integer :: lda,n,ipvt(n),info 4138 real*8 :: a(2,lda,n) 4139 ! 4140 ! zgefa factors a complex*16 matrix by gaussian elimination. 4141 ! 4142 ! dzgefa is usually called by zgeco, but it can be called 4143 ! directly with a saving in time if rcond is not needed. 4144 ! (time for zgeco) = (1 + 9/n)*(time for zgefa) . 4145 ! 4146 ! on entry 4147 ! 4148 ! a complex*16(lda, n) 4149 ! the matrix to be factored. 4150 ! 4151 ! lda integer 4152 ! the leading dimension of the array a . 4153 ! 4154 ! n integer 4155 ! the order of the matrix a . 4156 ! 4157 ! on return 4158 ! 4159 ! a an upper triangular matrix and the multipliers 4160 ! which were used to obtain it. 4161 ! the factorization can be written a = l*u where 4162 ! l is a product of permutation and unit lower 4163 ! triangular matrices and u is upper triangular. 4164 ! 4165 ! ipvt integer(n) 4166 ! an integer vector of pivot indices. 4167 ! 4168 ! info integer 4169 ! = 0 normal value. 4170 ! = k if u(k,k) .eq. 0.0 . this is not an error 4171 ! condition for this subroutine, but it does 4172 ! indicate that zgesl or zgedi will divide by zero 4173 ! if called. use rcond in zgeco for a reliable 4174 ! indication of singularity. 4175 ! 4176 ! linpack. this version dated 08/14/78 . 4177 ! cleve moler, university of new mexico, argonne national lab. 4178 ! 4179 ! subroutines and functions 4180 ! 4181 ! internal variables 4182 ! 4183 !Local variables 4184 real*8 :: r(2),rk(2),rlj(2) 4185 real*8 :: rinv2,rmax,rabs 4186 integer :: i,j,k,kp1,l,nm1 4187 4188 ! 4189 ! gaussian elimination with partial pivoting 4190 ! 4191 info = 0 4192 nm1 = n - 1 4193 if (nm1 .lt. 1) go to 70 4194 do 60 k = 1, nm1 4195 kp1 = k + 1 4196 ! 4197 ! find l = pivot index 4198 ! 4199 !l = izamax(n-k+1,a(k,k),1) + k - 1 4200 rmax=0.d0 4201 l=0 4202 do i=k,n 4203 rabs=abs(a(1,i,k))+abs(a(2,i,k)) 4204 if(rmax<=rabs) then 4205 rmax=rabs 4206 l=i 4207 end if 4208 end do 4209 ipvt(k) = l 4210 ! 4211 ! zero pivot implies this column already triangularized 4212 ! 4213 if (abs(a(1,l,k))+abs(a(2,l,k)) .eq. 0.0d0) go to 40 4214 ! 4215 ! interchange if necessary 4216 ! 4217 if (l .eq. k) go to 10 4218 r(1) = a(1,l,k); r(2) = a(2,l,k) 4219 a(1,l,k) = a(1,k,k); a(2,l,k) = a(2,k,k) 4220 a(1,k,k) = r(1); a(2,k,k) = r(2) 4221 10 continue 4222 ! 4223 ! compute multipliers 4224 ! 4225 rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2) 4226 rk(1) = -rinv2*a(1,k,k) 4227 rk(2) = rinv2*a(2,k,k) 4228 !call zscal(n-k,t,a(k+1,k),1) 4229 do i=k+1,n 4230 r(1)=a(1,i,k) 4231 r(2)=a(2,i,k) 4232 a(1,i,k)=rk(1)*r(1)-rk(2)*r(2) 4233 a(2,i,k)=rk(1)*r(2)+rk(2)*r(1) 4234 end do 4235 ! 4236 ! row elimination with column indexing 4237 ! 4238 do j = kp1, n 4239 rlj(1) = a(1,l,j); rlj(2) = a(2,l,j) 4240 if (l .eq. k) go to 20 4241 a(1,l,j) = a(1,k,j); a(2,l,j) = a(2,k,j) 4242 a(1,k,j) = rlj(1); a(2,k,j) = rlj(2) 4243 20 continue 4244 !call zaxpy(n-k,t,a(1,k+1,k),1,a(1,k+1,j),1) 4245 do i=k+1,n 4246 a(1,i,j)=rlj(1)*a(1,i,k)-rlj(2)*a(2,i,k)+a(1,i,j) 4247 a(2,i,j)=rlj(2)*a(1,i,k)+rlj(1)*a(2,i,k)+a(2,i,j) 4248 end do 4249 end do 4250 go to 50 4251 40 continue 4252 info = k 4253 50 continue 4254 60 continue 4255 70 continue 4256 ipvt(n) = n 4257 if (abs(a(1,n,n))+abs(a(2,n,n)) .eq. 0.0d0) info = n 4258 4259 end subroutine dzgefa
m_abilasi/lubksb [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
lubksb
FUNCTION
Solves the set of n linear equations A . X = B. Here a is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx is input as the permutation vector returned by ludcmp. b(1:n) is input as the right-hand side vector B, and returns with the solution vector X. a, n, np, and indx are not modified by this routine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that b will begin with many zero elements, so it is efficient for use in matrix inversion.
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
This routine is deprecated, use lapack API
PARENTS
CHILDREN
SOURCE
3856 SUBROUTINE lubksb(a,n,np,indx,b) 3857 3858 3859 !This section has been created automatically by the script Abilint (TD). 3860 !Do not modify the following lines by hand. 3861 #undef ABI_FUNC 3862 #define ABI_FUNC 'lubksb' 3863 !End of the abilint section 3864 3865 IMPLICIT NONE 3866 3867 INTEGER n,np,indx(n) 3868 REAL*8 a(np,np),b(n) 3869 3870 INTEGER i,ii,j,ll 3871 REAL*8 sum 3872 ! write(std_out,*) 'ENTERING LUBKSB...' 3873 ! write(std_out,201) ((a(i,j),j=1,n),i=1,n) 3874 ! 201 FORMAT('A in lubksb ',/,3F16.8,/,3F16.8,/,3F16.8) 3875 3876 ii=0 3877 do i=1,n 3878 ll=indx(i) 3879 sum=b(ll) 3880 b(ll)=b(i) 3881 if (ii.ne.0)then 3882 do j=ii,i-1 3883 sum=sum-a(i,j)*b(j) 3884 enddo 3885 else if (sum.ne.0.) then 3886 ii=i 3887 endif 3888 b(i)=sum 3889 enddo 3890 do i=n,1,-1 3891 sum=b(i) 3892 do j=i+1,n 3893 sum=sum-a(i,j)*b(j) 3894 enddo 3895 b(i)=sum/a(i,i) 3896 enddo 3897 ! write(std_out,*) 'LEAVING LUBKSB...' 3898 return 3899 3900 END SUBROUTINE LUBKSB
m_abilasi/ludcmp [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
ludcmp
FUNCTION
Given a matrix a(1:n,1:n), with physical dimension np by np, this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above; indx(1:n) is an output vector that records the row permutation effected by the partial pivoting; id is output as +- 1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix.
INPUTS
OUTPUT
NOTES
This routine is depreacted, use lapack API
PARENTS
CHILDREN
SOURCE
3742 SUBROUTINE ludcmp(a,n,np,indx,id,info) 3743 3744 3745 !This section has been created automatically by the script Abilint (TD). 3746 !Do not modify the following lines by hand. 3747 #undef ABI_FUNC 3748 #define ABI_FUNC 'ludcmp' 3749 !End of the abilint section 3750 3751 implicit none 3752 3753 INTEGER n,np,indx(n),NMAX,id,info 3754 REAL*8 a(np,np),TINY 3755 PARAMETER (NMAX=500,TINY=1.0e-20) 3756 3757 INTEGER i,imax,j,k 3758 REAL*8 aamax,dum,sum,vv(NMAX) 3759 3760 ! write(std_out,*) 'ENTERING LUDCMP...' 3761 ! write(std_out,*) 'in ludcmp n=',n,' np=',np 3762 ! write(std_out,201) ((a(i,j),j=1,n),i=1,n) 3763 ! 201 FORMAT('A in ludcmp ',/,3F16.8,/,3F16.8,/,3F16.8) 3764 id=1 3765 info=0 3766 do i=1,n 3767 aamax=0. 3768 do j=1,n 3769 if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 3770 enddo 3771 if (aamax.eq.0.) then 3772 write(std_out,*) 'LUDCMP: singular matrix !!!' 3773 do j=1,3 3774 write(std_out,*) (a(j,k),k=1,3) 3775 enddo 3776 info=1 3777 return 3778 ! stop 'singular matrix in ludcmp' 3779 endif 3780 vv(i)=1./aamax 3781 enddo 3782 do j=1,n 3783 do i=1,j-1 3784 sum=a(i,j) 3785 do k=1,i-1 3786 sum=sum-a(i,k)*a(k,j) 3787 enddo 3788 a(i,j)=sum 3789 enddo 3790 aamax=0. 3791 do i=j,n 3792 sum=a(i,j) 3793 do k=1,j-1 3794 sum=sum-a(i,k)*a(k,j) 3795 enddo 3796 a(i,j)=sum 3797 dum=vv(i)*abs(sum) 3798 if (dum.ge.aamax) then 3799 imax=i 3800 aamax=dum 3801 endif 3802 enddo 3803 if (j.ne.imax)then 3804 do k=1,n 3805 dum=a(imax,k) 3806 a(imax,k)=a(j,k) 3807 a(j,k)=dum 3808 enddo 3809 id=-id 3810 vv(imax)=vv(j) 3811 endif 3812 indx(j)=imax 3813 if(a(j,j).eq.0.)a(j,j)=TINY 3814 if(j.ne.n)then 3815 dum=1./a(j,j) 3816 do i=j+1,n 3817 a(i,j)=a(i,j)*dum 3818 enddo 3819 endif 3820 enddo 3821 ! write(std_out,*) 'LEAVING LUDCMP...' 3822 return 3823 END SUBROUTINE ludcmp
m_abilasi/matr3eigval [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
matr3eigval
FUNCTION
Find the eigenvalues of a real symmetric 3x3 matrix, entered in full storage mode.
INPUTS
matr(3,3)=real symmetric 3x3 matrix
OUTPUT
eigval(3)=three eigenvalues
PARENTS
chkdilatmx
CHILDREN
zhpev
SOURCE
3543 subroutine matr3eigval(eigval,matr) 3544 3545 3546 !This section has been created automatically by the script Abilint (TD). 3547 !Do not modify the following lines by hand. 3548 #undef ABI_FUNC 3549 #define ABI_FUNC 'matr3eigval' 3550 !End of the abilint section 3551 3552 implicit none 3553 3554 !Arguments ------------------------------------ 3555 !arrays 3556 real(dp),intent(in) :: matr(3,3) 3557 real(dp),intent(out) :: eigval(3) 3558 3559 !Local variables------------------------------- 3560 !scalars 3561 integer :: ier 3562 !arrays 3563 real(dp) :: eigvec(2,3,3),matrx(2,6),zhpev1(2,2*3-1),zhpev2(3*3-2) 3564 3565 ! ************************************************************************* 3566 3567 matrx(1,1)=matr(1,1) 3568 matrx(1,2)=matr(1,2) 3569 matrx(1,3)=matr(2,2) 3570 matrx(1,4)=matr(1,3) 3571 matrx(1,5)=matr(2,3) 3572 matrx(1,6)=matr(3,3) 3573 matrx(2,:)=zero 3574 3575 call ZHPEV ('V','U',3,matrx,eigval,eigvec,3,zhpev1,zhpev2,ier) 3576 !write(std_out,*)' eigval=',eigval 3577 3578 end subroutine matr3eigval
m_abilasi/matrginv [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
matrginv
FUNCTION
Invert a general matrix of real*8 elements.
INPUTS
lda=leading dimension of complex matrix a n=size of complex matrix a a=matrix of real elements
OUTPUT
a=inverse of a input matrix
SIDE EFFECTS
a(lda,n)= array of real elements, input, inverted at output
PARENTS
calc_optical_mels,ddb_elast,ddb_piezo,get_tau_k,linear_optics_paw m_haydock,m_vcoul,matpointsym,mka2f_tr,mlwfovlp_ylmfar,setup_bse strainsym
CHILDREN
dbgmdi,dbgmlu,dgeicd,dgetrf,dgetri
SOURCE
3410 subroutine matrginv(a,lda,n) 3411 3412 3413 !This section has been created automatically by the script Abilint (TD). 3414 !Do not modify the following lines by hand. 3415 #undef ABI_FUNC 3416 #define ABI_FUNC 'matrginv' 3417 !End of the abilint section 3418 3419 implicit none 3420 3421 !Arguments ------------------------------------ 3422 !scalars 3423 integer,intent(in) :: lda,n 3424 !arrays 3425 real(dp),intent(inout) :: a(lda,n) 3426 3427 !Local variables------------------------------- 3428 !scalars 3429 integer :: ierr,nwork 3430 #if defined HAVE_LINALG_ESSL 3431 real(dp) :: rcond 3432 #endif 3433 character(len=500) :: message 3434 !arrays 3435 integer,allocatable :: ipvt(:) 3436 #if defined HAVE_LINALG_ESSL 3437 real(dp) :: det(2) 3438 #elif defined HAVE_LINALG_ASL 3439 real(dp) :: det(2) 3440 #endif 3441 real(dp),allocatable :: work(:) 3442 3443 ! ************************************************************************* 3444 3445 #if defined HAVE_LINALG_ESSL 3446 nwork=200*n 3447 #else 3448 nwork=n 3449 #endif 3450 3451 ABI_ALLOCATE(work,(nwork)) 3452 ABI_ALLOCATE(ipvt,(n)) 3453 3454 3455 #if defined HAVE_LINALG_ESSL 3456 3457 call dgeicd(a,lda,n,0,rcond,det,work,nwork) 3458 if(abs(rcond)==zero) then 3459 write(message, '(10a)' ) ch10,& 3460 & ' matrginv : BUG -',ch10,& 3461 & ' The matrix that has been passed in argument of this subroutine',ch10,& 3462 & ' is probably either singular or nearly singular.',ch10,& 3463 & ' The ESSL routine dgeicd failed.',ch10,& 3464 & ' Action : Contact ABINIT group ' 3465 MSG_ERROR(message) 3466 end if 3467 3468 #elif defined HAVE_LINALG_ASL 3469 3470 call dbgmlu(a,lda,n,ipvt,ierr) 3471 if(ierr /= 0) then 3472 write(message, '(10a)' ) ch10,& 3473 & ' matrginv : BUG -',ch10,& 3474 & ' The matrix that has been passed in argument of this subroutine',ch10,& 3475 & ' is probably either singular or nearly singular.',ch10,& 3476 & ' The ASL routine dbgmlu failed.',ch10,& 3477 & ' Action : Contact ABINIT group ' 3478 MSG_ERROR(message) 3479 end if 3480 call dbgmdi(a,lda,n,ipvt,det,-1,work,ierr) 3481 if(ierr /= 0) then 3482 write(message, '(10a)' ) ch10,& 3483 & ' matrginv : BUG -',ch10,& 3484 & ' The matrix that has been passed in argument of this subroutine',ch10,& 3485 & ' is probably either singular or nearly singular.',ch10,& 3486 & ' The ASL routine dbgmdi failed.',ch10,& 3487 & ' Action : Contact ABINIT group ' 3488 MSG_ERROR(message) 3489 end if 3490 3491 #else 3492 3493 call dgetrf(n,n,a,lda,ipvt,ierr) 3494 if(ierr /= 0) then 3495 write(message, '(10a)' ) ch10,& 3496 & ' matrginv : BUG -',ch10,& 3497 & ' The matrix that has been passed in argument of this subroutine',ch10,& 3498 & ' is probably either singular or nearly singular.',ch10,& 3499 & ' The LAPACK routine dgetrf failed.',ch10,& 3500 & ' Action : Contact ABINIT group ' 3501 MSG_ERROR(message) 3502 end if 3503 call dgetri(n,a,lda,ipvt,work,n,ierr) 3504 if(ierr /= 0) then 3505 write(message, '(10a)' ) ch10,& 3506 & ' matrginv : BUG -',ch10,& 3507 & ' The matrix that has been passed in argument of this subroutine',ch10,& 3508 & ' is probably either singular or nearly singular.',ch10,& 3509 & ' The LAPACK routine dgetri failed.',ch10,& 3510 & ' Action : Contact ABINIT group ' 3511 MSG_ERROR(message) 3512 end if 3513 3514 #endif 3515 3516 ABI_DEALLOCATE(work) 3517 ABI_DEALLOCATE(ipvt) 3518 3519 end subroutine matrginv
m_abilasi/test_xginv [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
test_xginv
FUNCTION
m_abilasi/wrap_CGEEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_CGEEV
FUNCTION
wrap_CGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors using single precision arithmetic. [PRIVATE] The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.
INPUTS
JOBVL (input) CHARACTER*1
= 'N': left eigenvectors of A are not computed;
= 'V': left eigenvectors of are computed.
JOBVR (input) CHARACTER*1
= 'N': right eigenvectors of A are not computed;
= 'V': right eigenvectors of A are computed.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
LDVL (input) INTEGER
The leading dimension of the array VL. LDVL >= 1; if
JOBVL = 'V', LDVL >= N.
LDVR (input) INTEGER
The leading dimension of the array VR. LDVR >= 1; if
JOBVR = 'V', LDVR >= N.
OUTPUT
W (output) COMPLEX(SPC) array, dimension (N)
W contains the computed eigenvalues.
VL (output) COMPLEX(SCP) array, dimension (LDVL,N)
If JOBVL = 'V', the left eigenvectors u(j) are stored one
after another in the columns of VL, in the same order
as their eigenvalues.
If JOBVL = 'N', VL is not referenced.
u(j) = VL(:,j), the j-th column of VL.
VR (output) COMPLEX(SPC) array, dimension (LDVR,N)
If JOBVR = 'V', the right eigenvectors v(j) are stored one
after another in the columns of VR, in the same order
as their eigenvalues.
If JOBVR = 'N', VR is not referenced.
v(j) = VR(:,j), the j-th column of VR.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(SPC) array, dimension (LDA,N)
On entry, the N-by-N matrix A.
On exit, A has been overwritten.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
2598 subroutine wrap_CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr) 2599 2600 2601 !This section has been created automatically by the script Abilint (TD). 2602 !Do not modify the following lines by hand. 2603 #undef ABI_FUNC 2604 #define ABI_FUNC 'wrap_CGEEV' 2605 !End of the abilint section 2606 2607 implicit none 2608 2609 !Arguments ------------------------------------ 2610 !scalars 2611 integer,intent(in) :: n,lda,ldvl,ldvr 2612 character(len=*),intent(in) :: jobvl,jobvr 2613 !arrays 2614 complex(spc),intent(inout) :: a(lda,n) 2615 complex(spc),intent(out) :: w(n) 2616 complex(spc),intent(out) :: vl(ldvl,n) 2617 complex(spc),intent(out) :: vr(ldvr,n) 2618 2619 !Local variables ------------------------------ 2620 !scalars 2621 integer :: info,lwork 2622 character(len=500) :: msg 2623 !arrays 2624 real(sp),allocatable :: rwork(:) 2625 complex(spc),allocatable :: work(:) 2626 2627 !************************************************************************ 2628 2629 lwork = MAX(1,2*n) 2630 2631 ABI_MALLOC(work,(lwork)) 2632 ABI_MALLOC(rwork,(2*n)) 2633 2634 call CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info) 2635 2636 if (info < 0) then 2637 write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGEEV had an illegal value." 2638 MSG_ERROR(msg) 2639 end if 2640 2641 if (info > 0) then 2642 write(msg,'(3a,i0,a,i0,a)')& 2643 & "CGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,& 2644 & "Elements ",info+1,":",n," of W contain eigenvalues which have converged. " 2645 MSG_ERROR(msg) 2646 end if 2647 2648 ABI_FREE(work) 2649 ABI_FREE(rwork) 2650 2651 end subroutine wrap_CGEEV
m_abilasi/wrap_CHEEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_CHEEV
FUNCTION
wrap_CHEEV computes the eigenvalues and, optionally, the eigenvectors of a complex Hermitian matrix in single precision. [PRIVATE]
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
OUTPUT
W (output) REAL(SP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(SPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
217 subroutine wrap_CHEEV(jobz,uplo,n,a,w) 218 219 220 !This section has been created automatically by the script Abilint (TD). 221 !Do not modify the following lines by hand. 222 #undef ABI_FUNC 223 #define ABI_FUNC 'wrap_CHEEV' 224 !End of the abilint section 225 226 implicit none 227 228 !Arguments ------------------------------------ 229 !scalars 230 integer,intent(in) :: n 231 character(len=*),intent(in) :: jobz,uplo 232 !scalars 233 real(sp),intent(out) :: w(n) 234 complex(spc),intent(inout) :: a(n,n) 235 236 !Local variables ------------------------------ 237 !scalars 238 integer :: lwork,info 239 character(len=500) :: msg 240 !arrays 241 real(sp),allocatable :: rwork(:) 242 complex(spc),allocatable :: work(:) 243 244 !************************************************************************ 245 246 lwork = MAX(1,2*n-1) 247 248 ABI_MALLOC(work, (lwork)) 249 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 250 251 call CHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info) 252 253 if (info < 0) then 254 write(msg,'(a,i0,a)')"The ",-info,"-th argument of CHEEV had an illegal value." 255 MSG_ERROR(msg) 256 end if 257 258 if (info > 0) then 259 write(msg,'(2a,i0,a)')& 260 & "CHEEV: the algorithm failed to converge; ",ch10,& 261 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 262 MSG_ERROR(msg) 263 end if 264 265 ABI_FREE(rwork) 266 ABI_FREE(work) 267 268 !TODO scaLAPACK version (complex single precision buffer is needed in matrix_scalapack) 269 270 end subroutine wrap_CHEEV
m_abilasi/wrap_CHPEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_CHPEV
FUNCTION
wrap_CHPEV computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
OUTPUT
W (output) REAL(SP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
Z (output) COMPLEX(SPC) array, dimension (LDZ, N)
If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
eigenvectors of the matrix A, with the i-th column of Z
holding the eigenvector associated with W(i).
If JOBZ = 'N', then Z is not referenced.
See also SIDE EFFECTS
SIDE EFFECTS
AP (input/output) COMPLEX(SPC) array, dimension (N*(N+1)/2)
On entry, the upper or lower triangle of the Hermitian matrix
A, packed columnwise in a linear array. The j-th column of A
is stored in the array AP as follows:
if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
On exit, AP is overwritten by values generated during the
reduction to tridiagonal form. If UPLO = 'U', the diagonal
and first superdiagonal of the tridiagonal matrix T overwrite
the corresponding elements of A, and if UPLO = 'L', the
diagonal and first subdiagonal of T overwrite the
corresponding elements of A.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
710 subroutine wrap_CHPEV(jobz,uplo,n,ap,w,z,ldz) 711 712 713 !This section has been created automatically by the script Abilint (TD). 714 !Do not modify the following lines by hand. 715 #undef ABI_FUNC 716 #define ABI_FUNC 'wrap_CHPEV' 717 !End of the abilint section 718 719 implicit none 720 721 !Arguments ------------------------------------ 722 !scalars 723 integer,intent(in) :: n,ldz 724 character(len=*),intent(in) :: jobz,uplo 725 !arrays 726 real(sp),intent(out) :: w(n) 727 complex(spc),intent(inout) :: ap(n*(n+1)/2) 728 complex(spc),intent(out) :: z(ldz,n) 729 730 !Local variables ------------------------------ 731 !scalars 732 integer :: info 733 character(len=500) :: msg 734 !arrays 735 real(sp),allocatable :: rwork(:) 736 complex(spc),allocatable :: work(:) 737 738 !************************************************************************ 739 740 ABI_MALLOC(work, (MAX(1,2*n-1))) 741 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 742 743 call CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, INFO ) 744 745 if (info < 0) then 746 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value." 747 MSG_ERROR(msg) 748 end if 749 750 if (info > 0) then 751 write(msg,'(2a,i0,a)')& 752 & "ZHPEV: the algorithm failed to converge; ",ch10,& 753 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 754 MSG_ERROR(msg) 755 end if 756 757 ABI_FREE(rwork) 758 ABI_FREE(work) 759 760 end subroutine wrap_CHPEV
m_abilasi/wrap_DSYEV_ZHEEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_DSYEV_ZHEEV
FUNCTION
wrap_DSYEV_ZHEEV computes the eigenvalues and, optionally, the eigenvectors of a (complex Hermitian| real symmetric) matrix in double precision. [PRIVATE]
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX= Size of the first dimension of the A matrix.
1 for a real symmetric matrix.
2 for complex Hermitian matrix stored in a real array with real and imaginary part.
N (input) INTEGER
The order of the matrix A. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (complex Hermitian|Real symmetric) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
497 subroutine wrap_DSYEV_ZHEEV(jobz,uplo,cplex,n,a,w,comm) 498 499 500 !This section has been created automatically by the script Abilint (TD). 501 !Do not modify the following lines by hand. 502 #undef ABI_FUNC 503 #define ABI_FUNC 'wrap_DSYEV_ZHEEV' 504 !End of the abilint section 505 506 implicit none 507 508 !Arguments ------------------------------------ 509 !scalars 510 integer,intent(in) :: n,cplex 511 integer,optional,intent(in) :: comm 512 character(len=*),intent(in) :: jobz,uplo 513 !arrays 514 real(dp),intent(inout) :: a(cplex,n,n) 515 real(dp),intent(out) :: w(n) 516 517 !Local variables ------------------------------ 518 !scalars 519 integer :: lwork,info,nprocs 520 logical :: use_scalapack 521 character(len=500) :: msg 522 !arrays 523 real(dp),allocatable :: rwork(:) 524 real(dp),allocatable :: work_real(:) 525 complex(dpc),allocatable :: work_cplx(:) 526 #ifdef HAVE_LINALG_SCALAPACK 527 integer :: ierr,istwf_k,tbloc 528 logical :: want_eigenvectors 529 type(matrix_scalapack) :: Slk_mat,Slk_vec 530 type(processor_scalapack) :: Slk_processor 531 #endif 532 !************************************************************************ 533 534 use_scalapack=.FALSE. 535 if (PRESENT(comm)) then 536 nprocs = xmpi_comm_size(comm) 537 #ifdef HAVE_LINALG_SCALAPACK 538 use_scalapack = (nprocs>1) 539 #endif 540 end if 541 542 if (ALL(cplex/=(/1,2/))) then 543 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 544 MSG_BUG(msg) 545 end if 546 547 SELECT CASE(use_scalapack) 548 549 CASE (.FALSE.) 550 551 if (cplex==1) then ! Real symmetric case. 552 553 lwork = MAX(1,3*n-1) 554 555 ABI_MALLOC(work_real,(lwork)) 556 557 call DSYEV(jobz,uplo,n,a,n,w,work_real,lwork,info) 558 559 if (info < 0) then 560 write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYEV had an illegal value." 561 MSG_ERROR(msg) 562 end if 563 564 if (info > 0) then 565 write(msg,'(2a,i0,a)')& 566 & "DSYEV: the algorithm failed to converge; ",ch10,& 567 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 568 MSG_ERROR(msg) 569 end if 570 571 ABI_FREE(work_real) 572 573 RETURN 574 575 else ! Hermitian case. 576 577 lwork = MAX(1,2*n-1) 578 579 ABI_MALLOC(work_cplx, (lwork)) 580 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 581 582 call ZHEEV(jobz,uplo,n,a,n,w,work_cplx,lwork,rwork,info) 583 584 if (info < 0) then 585 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value." 586 MSG_ERROR(msg) 587 end if 588 589 if (info > 0) then 590 write(msg,'(2a,i0,a)')& 591 & "ZHEEV: the algorithm failed to converge; ",ch10,& 592 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 593 MSG_ERROR(msg) 594 end if 595 596 ABI_FREE(rwork) 597 ABI_FREE(work_cplx) 598 599 RETURN 600 end if ! cplex 601 602 CASE (.TRUE.) 603 604 #ifdef HAVE_LINALG_SCALAPACK 605 606 MSG_ERROR("Not coded yet") 607 608 ! call init_scalapack(Slk_processor,comm) 609 ! istwf_k=1 610 ! 611 ! ! Initialize and fill Scalapack matrix from the global one. 612 ! tbloc=SLK_BLOCK_SIZE 613 ! call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 614 ! 615 ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 616 ! call wrtout(std_out,msg,"PERS") 617 ! 618 ! call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 619 ! 620 ! want_eigenvectors = firstchar(jobz,(/"V","v"/)) 621 ! if (want_eigenvectors) then ! Initialize the distributed vectors. 622 ! call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 623 ! end if 624 ! 625 ! ! Solve the problem with scaLAPACK. 626 ! call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w) 627 ! 628 ! call destruction_matrix_scalapack(Slk_mat) 629 ! 630 ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors 631 ! a = czero 632 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node. 633 ! call destruction_matrix_scalapack(Slk_vec) 634 ! call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 635 ! end if 636 ! 637 ! call end_scalapack(Slk_processor) 638 639 RETURN 640 #endif 641 642 MSG_BUG("You should not be here!") 643 644 END SELECT 645 646 end subroutine wrap_DSYEV_ZHEEV
m_abilasi/wrap_DSYEVX_ZHEEVX [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_DSYEVX_ZHEEVX
FUNCTION
wrap_DSYEVX_ZHEEVX computes selected eigenvalues and, optionally, eigenvectors of a (real symmetric|complex Hermitian) matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX Size of the first dimension of the matrix A.
1 for real symmetric matrix
2 for complex Hermitian matrix.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
See "Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," by Demmel and
Kahan, LAPACK Working Note #3.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) REAL(DP) array, dimension (CPLEX, LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (real symmetric|complex Hermitian) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO='L') or the upper
triangle (if UPLO='U') of A, including the diagonal, is
destroyed.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
1763 subroutine wrap_DSYEVX_ZHEEVX(jobz,range,uplo,cplex,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 1764 1765 1766 !This section has been created automatically by the script Abilint (TD). 1767 !Do not modify the following lines by hand. 1768 #undef ABI_FUNC 1769 #define ABI_FUNC 'wrap_DSYEVX_ZHEEVX' 1770 !End of the abilint section 1771 1772 implicit none 1773 1774 !Arguments ------------------------------------ 1775 !scalars 1776 integer,intent(in) :: il,iu,ldz,n,cplex 1777 integer,optional,intent(in) :: comm 1778 integer,intent(inout) :: m 1779 real(dp),intent(in) :: abstol,vl,vu 1780 character(len=*),intent(in) :: jobz,range,uplo 1781 !arrays 1782 real(dp),intent(out) :: w(n) 1783 !real(dp),intent(out) :: z(cplex,ldz,n) 1784 real(dp),intent(out) :: z(cplex,ldz,m) 1785 real(dp),intent(inout) :: a(cplex,n,n) 1786 1787 !Local variables ------------------------------ 1788 !scalars 1789 integer :: lwork,info,nprocs 1790 logical :: use_scalapack 1791 character(len=500) :: msg 1792 !arrays 1793 integer,allocatable :: ifail(:),iwork(:) 1794 real(dp),allocatable :: rwork(:) 1795 real(dp),allocatable :: work_real(:) 1796 complex(dpc),allocatable :: work_cplx(:) 1797 #ifdef HAVE_LINALG_SCALAPACK 1798 integer :: ierr,istwf_k,tbloc 1799 logical :: want_eigenvectors 1800 type(matrix_scalapack) :: Slk_mat,Slk_vec 1801 type(processor_scalapack) :: Slk_processor 1802 #endif 1803 1804 !************************************************************************ 1805 1806 use_scalapack=.FALSE. 1807 if (PRESENT(comm)) then 1808 nprocs = xmpi_comm_size(comm) 1809 #ifdef HAVE_LINALG_SCALAPACK 1810 use_scalapack = (nprocs>1) 1811 #endif 1812 end if 1813 1814 if (ALL(cplex/=(/1,2/))) then 1815 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 1816 MSG_ERROR(msg) 1817 end if 1818 1819 SELECT CASE(use_scalapack) 1820 1821 CASE (.FALSE.) ! Standard LAPACK call. 1822 1823 if (cplex==1) then ! Real symmetric case 1824 1825 lwork = MAX(1,8*n) 1826 1827 ABI_MALLOC(work_real,(lwork)) 1828 ABI_MALLOC(iwork,(5*n)) 1829 ABI_MALLOC(ifail,(n)) 1830 1831 call DSYEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,info) 1832 1833 if (info < 0) then 1834 write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYEVX had an illegal value." 1835 MSG_ERROR(msg) 1836 end if 1837 1838 if (info > 0) then 1839 write(msg,'(2a,i0,a)')& 1840 & "DSYEVX: the algorithm failed to converge; ",ch10,& 1841 & info,"eigenvectors failed to converge. " 1842 MSG_ERROR(msg) 1843 end if 1844 1845 ABI_FREE(work_real) 1846 ABI_FREE(iwork) 1847 ABI_FREE(ifail) 1848 1849 RETURN 1850 1851 else ! Complex Hermitian case. 1852 1853 lwork = MAX(1,2*n) 1854 1855 ABI_MALLOC(work_cplx,(lwork)) 1856 ABI_MALLOC(rwork,(7*n)) 1857 ABI_MALLOC(iwork,(5*n)) 1858 ABI_MALLOC(ifail,(n)) 1859 1860 call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,info) 1861 1862 if (info < 0) then 1863 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value." 1864 MSG_ERROR(msg) 1865 end if 1866 1867 if (info > 0) then 1868 write(msg,'(2a,i0,a)')& 1869 & "ZHEEVX: the algorithm failed to converge; ",ch10,& 1870 & info,"eigenvectors failed to converge. " 1871 MSG_ERROR(msg) 1872 end if 1873 1874 ABI_FREE(iwork) 1875 ABI_FREE(ifail) 1876 ABI_FREE(rwork) 1877 ABI_FREE(work_cplx) 1878 1879 RETURN 1880 1881 end if 1882 1883 CASE (.TRUE.) 1884 1885 #ifdef HAVE_LINALG_SCALAPACK 1886 1887 MSG_ERROR("Not coded yet") 1888 ! call init_scalapack(Slk_processor,comm) 1889 ! istwf_k=1 1890 1891 ! ! Initialize and fill Scalapack matrix from the global one. 1892 ! tbloc=SLK_BLOCK_SIZE 1893 ! call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1894 1895 ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 1896 ! call wrtout(std_out,msg,"COLL") 1897 1898 ! call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 1899 1900 ! want_eigenvectors = firstchar(jobz,(/"V","v"/)) 1901 ! if (want_eigenvectors) then ! Initialize the distributed vectors. 1902 ! call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1903 ! end if 1904 1905 ! ! Solve the problem. 1906 ! call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w) 1907 1908 ! call destruction_matrix_scalapack(Slk_mat) 1909 ! 1910 ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors 1911 ! z = czero 1912 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 1913 ! call destruction_matrix_scalapack(Slk_vec) 1914 ! call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 1915 ! end if 1916 1917 ! call end_scalapack(Slk_processor) 1918 1919 RETURN 1920 #endif 1921 1922 MSG_BUG("You should not be here!") 1923 1924 END SELECT 1925 1926 end subroutine wrap_DSYEVX_ZHEEVX
m_abilasi/wrap_DSYGV_ZHEGV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_DSYGV_ZHEGV
FUNCTION
wrap_DSYGV_ZHEGV computes all the eigenvalues, and optionally, the eigenvectors of a
(real generalized symmetric-definite| complex generalized Hermitian-definite)
eigenproblem, of the form
A*x=(lambda)*B*x (1),
A*Bx=(lambda)*x, (2), or
B*A*x=(lambda)*x (3).
Here A and B are assumed to be (symmetric|Hermitian) and B is also positive definite.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= "N": Compute eigenvalues only;
= "V": Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= "U": Upper triangle of A and B are stored;
= "L": Lower triangle of A and B are stored.
CPLEX Size of the first dimension of the A and B matrices.
1 for a real symmetric matrix.
2 for a complex Hermitian matrix.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX,N, N)
On entry, the (real symmetric|Hermitian) matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
<S-F1>contains the upper triangular part of the matrix A.
If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
The eigenvectors are normalized as follows:
if ITYPE =1 or 2:
Z**T*B*Z = I if CPLEX=1
Z**H*B*Z = I if CPLEX=2.
if ITYPE = 3,
Z**T*inv(B)*Z = I if CPLEX=1
Z**H*inv(B)*Z = I if CPLEX=2
If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
(if UPLO="L") of A, including the diagonal, is destroyed.
B (input/output) REAL(DP) array, dimension (CPLEX,N, N)
On entry, the (real symmetric|Hermitian) positive definite matrix B.
If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
factor U or L from the Cholesky factorization
B = U**T*U or B = L*L**T if CPLEX=1
B = U**H*U or B = L*L**H if CPLEX=2
PARENTS
CHILDREN
cwtime,xginv
SOURCE
1225 subroutine wrap_DSYGV_ZHEGV(itype,jobz,uplo,cplex,n,a,b,w,comm) 1226 1227 1228 !This section has been created automatically by the script Abilint (TD). 1229 !Do not modify the following lines by hand. 1230 #undef ABI_FUNC 1231 #define ABI_FUNC 'wrap_DSYGV_ZHEGV' 1232 !End of the abilint section 1233 1234 implicit none 1235 1236 !Arguments ------------------------------------ 1237 !scalars 1238 integer,intent(in) :: n,itype,cplex 1239 integer,optional,intent(in) :: comm 1240 character(len=*),intent(in) :: jobz,uplo 1241 !arrays 1242 real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n) 1243 real(dp),intent(out) :: w(n) 1244 1245 !Local variables ------------------------------ 1246 !scalars 1247 integer :: lwork,info,nprocs,ii 1248 logical :: use_scalapack 1249 character(len=500) :: msg 1250 !arrays 1251 real(dp),allocatable :: rwork(:) 1252 real(dp),allocatable :: work_real(:) 1253 complex(dpc),allocatable :: work_cplx(:) 1254 #ifdef HAVE_LINALG_SCALAPACK 1255 integer :: ierr,istwf_k,tbloc 1256 type(matrix_scalapack) :: Slk_matA,Slk_matB 1257 type(processor_scalapack) :: Slk_processor 1258 #endif 1259 !************************************************************************ 1260 1261 use_scalapack=.FALSE. 1262 if (PRESENT(comm)) then 1263 nprocs = xmpi_comm_size(comm) 1264 #ifdef HAVE_LINALG_SCALAPACK 1265 use_scalapack = (nprocs>1) 1266 #endif 1267 end if 1268 1269 if (ALL(cplex/=(/1,2/))) then 1270 write(msg,'(a,i0)')"Wrong value for cplex: ",cplex 1271 MSG_ERROR(msg) 1272 end if 1273 1274 SELECT CASE(use_scalapack) 1275 1276 CASE (.FALSE.) 1277 1278 if (cplex==1) then ! Real symmetric case. 1279 1280 lwork = MAX(1,3*n-1) 1281 1282 ABI_MALLOC(work_real,(lwork)) 1283 1284 call DSYGV(itype,jobz,uplo,n,a,n,b,n,w,work_real,lwork,info) 1285 1286 if (info < 0) then 1287 write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYGV had an illegal value." 1288 MSG_ERROR(msg) 1289 end if 1290 1291 if (info > 0) then 1292 if (info<= n) then 1293 write(msg,'(2a,i0,a)')& 1294 & " DSYGV failed to converge: ",ch10,& 1295 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 1296 else 1297 ii = info -n 1298 write(msg,'(3a,i0,3a)')& 1299 & "DSYGV failed to converge: ",ch10,& 1300 & "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1301 & "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1302 end if 1303 MSG_ERROR(msg) 1304 end if 1305 1306 ABI_FREE(work_real) 1307 1308 RETURN 1309 1310 else ! complex Hermitian case 1311 1312 lwork = MAX(1,2*n-1) 1313 1314 ABI_MALLOC(work_cplx,(lwork)) 1315 ABI_MALLOC(rwork,(MAX(1,3*n-2))) 1316 1317 call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work_cplx,lwork,rwork,info) 1318 1319 if (info < 0) then 1320 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value." 1321 MSG_ERROR(msg) 1322 end if 1323 1324 if (info > 0) then 1325 if (info<= n) then 1326 write(msg,'(2a,i0,a)')& 1327 & "ZHEEV failed to converge: ",ch10,& 1328 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 1329 else 1330 ii = info -n 1331 write(msg,'(3a,i0,3a)')& 1332 & "ZHEEV failed to converge: ",ch10,& 1333 & "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1334 & "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1335 end if 1336 MSG_ERROR(msg) 1337 end if 1338 1339 ABI_FREE(rwork) 1340 ABI_FREE(work_cplx) 1341 1342 RETURN 1343 1344 end if ! cplex 1345 1346 CASE (.TRUE.) 1347 1348 #ifdef HAVE_LINALG_SCALAPACK 1349 1350 MSG_ERROR("Not coded yet") 1351 ! call init_scalapack(Slk_processor,comm) 1352 ! istwf_k=1 1353 1354 ! ! Initialize and fill Scalapack matrix from the global one. 1355 ! tbloc=SLK_BLOCK_SIZE 1356 1357 ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 1358 ! call wrtout(std_out,msg,"COLL") 1359 1360 ! call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1361 ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 1362 1363 ! call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1364 ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 1365 1366 ! ! Solve the problem with scaLAPACK. 1367 ! MSG_ERROR("slk_pZHEGV not yet coded") 1368 ! ! TODO 1369 ! call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w) 1370 1371 ! call destruction_matrix_scalapack(Slk_matB) 1372 ! 1373 ! if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 1374 ! a = czero 1375 ! call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node. 1376 ! call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 1377 ! end if 1378 1379 ! call destruction_matrix_scalapack(Slk_matA) 1380 1381 ! call end_scalapack(Slk_processor) 1382 1383 RETURN 1384 #endif 1385 1386 MSG_BUG("You should not be here!") 1387 1388 END SELECT 1389 1390 end subroutine wrap_DSYGV_ZHEGV
m_abilasi/wrap_DSYGVX_ZHEGVX [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_DSYGVX_ZHEGVX
FUNCTION
wrap_DSYGVX_ZHEGVX - compute selected eigenvalues, and optionally, eigenvectors of a (real symmetric-definite|complex generalized Hermitian-definite) eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and B are assumed to be (real symmetric|complex Hermitian) and B is also positive definite. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
CPLEX Size of the first dimension of the matrices A and B
1 for Real symmetric matrices
2 for complex Hermitianmatrices
N (input) INTEGER
The order of the matrices A and B. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) REAL(DP) array, dimension (CPLEX ,LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
The eigenvectors are normalized as follows:
if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) REAL(DP) array, dimension (CPLEX, N, N)
On entry, the (real symmetric| complex Hermitian) matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = "L",
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO="L") or the upper
triangle (if UPLO="U") of A, including the diagonal, is
destroyed.
B (input/output) REAL(DP) array, dimension (CPLEX, LDB, N)
On entry, the (real symmetric| complex Hermitian) matrix B. If UPLO = "U", the leading N-by-N upper triangular part
of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
2339 subroutine wrap_DSYGVX_ZHEGVX(itype,jobz,range,uplo,cplex,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 2340 2341 2342 !This section has been created automatically by the script Abilint (TD). 2343 !Do not modify the following lines by hand. 2344 #undef ABI_FUNC 2345 #define ABI_FUNC 'wrap_DSYGVX_ZHEGVX' 2346 !End of the abilint section 2347 2348 implicit none 2349 2350 !Arguments ------------------------------------ 2351 !scalars 2352 integer,intent(in) :: il,iu,ldz,n,itype,cplex 2353 integer,optional,intent(in) :: comm 2354 integer,intent(inout) :: m 2355 real(dp),intent(in) :: abstol,vl,vu 2356 character(len=*),intent(in) :: jobz,range,uplo 2357 !arrays 2358 real(dp),intent(out) :: w(n) 2359 !real(dp),intent(out) :: z(cplex,ldz,n) 2360 real(dp),intent(out) :: z(cplex,ldz,m) 2361 real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n) 2362 2363 !Local variables ------------------------------ 2364 !scalars 2365 integer :: lwork,info,nprocs,ii 2366 logical :: use_scalapack 2367 character(len=500) :: msg 2368 !arrays 2369 integer,allocatable :: ifail(:),iwork(:) 2370 real(dp),allocatable :: rwork(:) 2371 real(dp),allocatable :: work_real(:) 2372 complex(dpc),allocatable :: work_cplx(:) 2373 #ifdef HAVE_LINALG_SCALAPACK 2374 integer :: ierr,istwf_k,tbloc 2375 logical :: want_eigenvectors 2376 type(matrix_scalapack) :: Slk_matA,Slk_matB,Slk_vec 2377 type(processor_scalapack) :: Slk_processor 2378 #endif 2379 2380 !************************************************************************ 2381 2382 use_scalapack=.FALSE. 2383 if (PRESENT(comm)) then 2384 nprocs = xmpi_comm_size(comm) 2385 #ifdef HAVE_LINALG_SCALAPACK 2386 use_scalapack = (nprocs>1) 2387 #endif 2388 end if 2389 2390 if (ALL(cplex/=(/1,2/))) then 2391 write(msg,'(a,i0)')" Wrong value for cplex: ",cplex 2392 MSG_ERROR(msg) 2393 end if 2394 2395 SELECT CASE(use_scalapack) 2396 2397 CASE (.FALSE.) ! Standard LAPACK call. 2398 2399 if (cplex==1) then ! Real symmetric case 2400 2401 lwork = MAX(1,8*n) 2402 2403 ABI_MALLOC(work_real,(lwork)) 2404 ABI_MALLOC(iwork,(5*n)) 2405 ABI_MALLOC(ifail,(n)) 2406 2407 call DSYGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,info) 2408 2409 if (info < 0) then 2410 write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYGVX had an illegal value." 2411 MSG_ERROR(msg) 2412 end if 2413 2414 if (info > 0) then 2415 if (info<= n) then 2416 write(msg,'(a,i0,a)')& 2417 & " DSYGVX failed to converge: ",info," eigenvectors failed to converge. " 2418 else 2419 ii = info -n 2420 write(msg,'(3a,i0,3a)')& 2421 & " DSYGVX failed to converge: ",ch10,& 2422 & " The leading minor of order ",ii," of B is not positive definite. ",ch10,& 2423 & " The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 2424 end if 2425 MSG_ERROR(msg) 2426 end if 2427 2428 ABI_FREE(iwork) 2429 ABI_FREE(ifail) 2430 ABI_FREE(work_real) 2431 2432 RETURN 2433 2434 else ! Complex Hermitian case. 2435 2436 lwork = MAX(1,2*n) 2437 2438 ABI_MALLOC(work_cplx,(lwork)) 2439 ABI_MALLOC(rwork,(7*n)) 2440 ABI_MALLOC(iwork,(5*n)) 2441 ABI_MALLOC(ifail,(n)) 2442 2443 !write(std_out,*)"Calling ZHEGVX" 2444 2445 call ZHEGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,info) 2446 2447 if (info < 0) then 2448 write(msg,'(a,i0,a)')"The ",-info,"-th argument of ZHEGVX had an illegal value." 2449 MSG_ERROR(msg) 2450 end if 2451 2452 if (info > 0) then 2453 if (info<= n) then 2454 write(msg,'(a,i0,a)')& 2455 & "ZHEGVX failed to converge: ",info," eigenvectors failed to converge. " 2456 else 2457 ii = info -n 2458 write(msg,'(3a,i0,3a)')& 2459 & "ZHEEVX failed to converge: ",ch10,& 2460 & "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 2461 & "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 2462 end if 2463 MSG_ERROR(msg) 2464 end if 2465 2466 ABI_FREE(iwork) 2467 ABI_FREE(ifail) 2468 ABI_FREE(rwork) 2469 ABI_FREE(work_cplx) 2470 2471 RETURN 2472 2473 end if ! cplex 2474 2475 CASE (.TRUE.) 2476 2477 #ifdef HAVE_LINALG_SCALAPACK 2478 2479 MSG_ERROR("not coded yet") 2480 ! call init_scalapack(Slk_processor,comm) 2481 ! istwf_k=1 2482 2483 ! ! Initialize and fill Scalapack matrix from the global one. 2484 ! tbloc=SLK_BLOCK_SIZE 2485 2486 ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 2487 ! call wrtout(std_out,msg,"COLL") 2488 2489 ! call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2490 ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 2491 2492 ! call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2493 ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 2494 2495 ! want_eigenvectors = firstchar(jobz,(/"V","v"/)) 2496 ! if (want_eigenvectors) then ! Initialize the distributed vectors. 2497 ! call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2498 ! end if 2499 2500 ! ! Solve the problem. 2501 ! MSG_ERROR("slk_pZHEGVX not coded yet") 2502 ! ! TODO write the scaLAPACK wrapper. 2503 ! call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w) 2504 2505 ! call destruction_matrix_scalapack(Slk_matA) 2506 ! call destruction_matrix_scalapack(Slk_matB) 2507 ! 2508 ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors 2509 ! z = czero 2510 ! call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 2511 ! call destruction_matrix_scalapack(Slk_vec) 2512 ! call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 2513 ! end if 2514 2515 ! call end_scalapack(Slk_processor) 2516 2517 ! RETURN 2518 #endif 2519 2520 MSG_BUG("You should not be here!") 2521 2522 END SELECT 2523 2524 end subroutine wrap_DSYGVX_ZHEGVX
m_abilasi/wrap_ZGEEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZGEEV
FUNCTION
wrap_ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors using double precision arithmetic. [PRIVATE] The right eigenvector v(j) of A satisfies: A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j). The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. No scalapack version is available (PZGEEV is not provided by the Scalapack team)
INPUTS
JOBVL (input) CHARACTER*1
= 'N': left eigenvectors of A are not computed;
= 'V': left eigenvectors of are computed.
JOBVR (input) CHARACTER*1
= 'N': right eigenvectors of A are not computed;
= 'V': right eigenvectors of A are computed.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
LDVL (input) INTEGER
The leading dimension of the array VL. LDVL >= 1; if
JOBVL = 'V', LDVL >= N.
LDVR (input) INTEGER
The leading dimension of the array VR. LDVR >= 1; if
JOBVR = 'V', LDVR >= N.
OUTPUT
W (output) COMPLEX(DPC) array, dimension (N)
W contains the computed eigenvalues.
VL (output) COMPLEX(DPC) array, dimension (LDVL,N)
If JOBVL = 'V', the left eigenvectors u(j) are stored one
after another in the columns of VL, in the same order
as their eigenvalues.
If JOBVL = 'N', VL is not referenced.
u(j) = VL(:,j), the j-th column of VL.
VR (output) COMPLEX(DPC) array, dimension (LDVR,N)
If JOBVR = 'V', the right eigenvectors v(j) are stored one
after another in the columns of VR, in the same order
as their eigenvalues.
If JOBVR = 'N', VR is not referenced.
v(j) = VR(:,j), the j-th column of VR.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (LDA,N)
On entry, the N-by-N matrix A.
On exit, A has been overwritten.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
2726 subroutine wrap_ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr) 2727 2728 2729 !This section has been created automatically by the script Abilint (TD). 2730 !Do not modify the following lines by hand. 2731 #undef ABI_FUNC 2732 #define ABI_FUNC 'wrap_ZGEEV' 2733 !End of the abilint section 2734 2735 implicit none 2736 2737 !Arguments ------------------------------------ 2738 !scalars 2739 integer,intent(in) :: n,lda,ldvl,ldvr 2740 character(len=*),intent(in) :: jobvl,jobvr 2741 !arrays 2742 complex(dpc),intent(inout) :: a(lda,n) 2743 complex(dpc),intent(out) :: w(n) 2744 complex(dpc),intent(out) :: vl(ldvl,n) 2745 complex(dpc),intent(out) :: vr(ldvr,n) 2746 2747 !Local variables ------------------------------ 2748 !scalars 2749 integer :: info,lwork 2750 logical :: use_scalapack 2751 character(len=500) :: msg 2752 !arrays 2753 real(dp),allocatable :: rwork(:) 2754 complex(dpc),allocatable :: work(:) 2755 2756 !************************************************************************ 2757 2758 use_scalapack=.FALSE. 2759 2760 SELECT CASE(use_scalapack) 2761 2762 CASE (.FALSE.) 2763 2764 lwork = MAX(1,2*n) 2765 2766 ABI_MALLOC(work,(lwork)) 2767 ABI_MALLOC(rwork,(2*n)) 2768 2769 call ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info) 2770 2771 if (info < 0) then 2772 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGEEV had an illegal value." 2773 MSG_ERROR(msg) 2774 end if 2775 2776 if (info > 0) then 2777 write(msg,'(3a,i0,a,i0,a)')& 2778 & "ZGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,& 2779 & "Elements ",info+1,":",n," of W contain eigenvalues which have converged. " 2780 MSG_ERROR(msg) 2781 end if 2782 2783 ABI_FREE(work) 2784 ABI_FREE(rwork) 2785 2786 RETURN 2787 2788 CASE (.TRUE.) 2789 2790 MSG_BUG("You should not be here!") 2791 2792 END SELECT 2793 2794 end subroutine wrap_ZGEEV
m_abilasi/wrap_ZHEEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZHEEV
FUNCTION
wrap_ZHEEV computes the eigenvalues and, optionally, the eigenvectors of a complex Hermitian matrix in double precision. [PRIVATE]
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, if JOBZ = 'V', then if INFO = 0, A contains the
orthonormal eigenvectors of the matrix A.
If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
or the upper triangle (if UPLO='U') of A, including the
diagonal, is destroyed.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
324 subroutine wrap_ZHEEV(jobz,uplo,n,a,w,comm) 325 326 327 !This section has been created automatically by the script Abilint (TD). 328 !Do not modify the following lines by hand. 329 #undef ABI_FUNC 330 #define ABI_FUNC 'wrap_ZHEEV' 331 use interfaces_14_hidewrite 332 !End of the abilint section 333 334 implicit none 335 336 !Arguments ------------------------------------ 337 !scalars 338 integer,intent(in) :: n 339 integer,optional,intent(in) :: comm 340 character(len=*),intent(in) :: jobz,uplo 341 !arrays 342 complex(dpc),intent(inout) :: a(n,n) 343 real(dp),intent(out) :: w(n) 344 345 !Local variables ------------------------------ 346 !scalars 347 integer :: lwork,info,nprocs 348 logical :: use_scalapack 349 character(len=500) :: msg 350 !arrays 351 real(dp),allocatable :: rwork(:) 352 complex(dpc),allocatable :: work(:) 353 #ifdef HAVE_LINALG_SCALAPACK 354 integer :: ierr,istwf_k,tbloc 355 logical :: want_eigenvectors 356 type(matrix_scalapack) :: Slk_mat,Slk_vec 357 type(processor_scalapack) :: Slk_processor 358 #endif 359 !************************************************************************ 360 361 use_scalapack=.FALSE. 362 if (PRESENT(comm)) then 363 nprocs = xmpi_comm_size(comm) 364 #ifdef HAVE_LINALG_SCALAPACK 365 use_scalapack = (nprocs>1) 366 #endif 367 end if 368 369 SELECT CASE(use_scalapack) 370 371 CASE (.FALSE.) 372 373 lwork = MAX(1,2*n-1) 374 375 ABI_MALLOC(work, (lwork)) 376 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 377 378 call ZHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info) 379 380 if (info < 0) then 381 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value." 382 MSG_ERROR(msg) 383 end if 384 385 if (info > 0) then 386 write(msg,'(2a,i0,a)')& 387 & "ZHEEV: the algorithm failed to converge; ",ch10,& 388 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 389 MSG_ERROR(msg) 390 end if 391 392 ABI_FREE(rwork) 393 ABI_FREE(work) 394 395 RETURN 396 397 CASE (.TRUE.) 398 399 #ifdef HAVE_LINALG_SCALAPACK 400 call init_scalapack(Slk_processor,comm) 401 istwf_k=1 402 403 ! Initialize and fill Scalapack matrix from the global one. 404 tbloc=SLK_BLOCK_SIZE 405 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 406 407 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 408 call wrtout(std_out,msg,"COLL") 409 410 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 411 412 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 413 if (want_eigenvectors) then ! Initialize the distributed vectors. 414 call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 415 end if 416 417 ! Solve the problem with scaLAPACK. 418 call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w) 419 420 call destruction_matrix_scalapack(Slk_mat) 421 422 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 423 a = czero 424 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node. 425 call destruction_matrix_scalapack(Slk_vec) 426 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 427 end if 428 429 call end_scalapack(Slk_processor) 430 431 RETURN 432 #endif 433 434 MSG_BUG("You should not be here!") 435 436 END SELECT 437 438 end subroutine wrap_ZHEEV
m_abilasi/wrap_ZHEEVX [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZHEEVX
FUNCTION
wrap_ZHEEVX computes selected eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
See "Computing Small Singular Values of Bidiagonal Matrices
with Guaranteed High Relative Accuracy," by Demmel and
Kahan, LAPACK Working Note #3.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L',
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO='L') or the upper
triangle (if UPLO='U') of A, including the diagonal, is
destroyed.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
1513 subroutine wrap_ZHEEVX(jobz,range,uplo,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 1514 1515 1516 !This section has been created automatically by the script Abilint (TD). 1517 !Do not modify the following lines by hand. 1518 #undef ABI_FUNC 1519 #define ABI_FUNC 'wrap_ZHEEVX' 1520 use interfaces_14_hidewrite 1521 !End of the abilint section 1522 1523 implicit none 1524 1525 !Arguments ------------------------------------ 1526 !scalars 1527 integer,intent(in) :: il,iu,ldz,n 1528 integer,optional,intent(in) :: comm 1529 integer,intent(inout) :: m 1530 real(dp),intent(in) :: abstol,vl,vu 1531 character(len=*),intent(in) :: jobz,range,uplo 1532 !arrays 1533 real(dp),intent(out) :: w(n) 1534 complex(dpc),intent(out) :: z(ldz,m) 1535 complex(dpc),intent(inout) :: a(n,n) 1536 1537 !Local variables ------------------------------ 1538 !scalars 1539 integer :: lwork,info,nprocs 1540 logical :: use_scalapack 1541 character(len=500) :: msg 1542 !arrays 1543 integer,allocatable :: ifail(:),iwork(:) 1544 real(dp),allocatable :: rwork(:) 1545 complex(dpc),allocatable :: work(:) 1546 #ifdef HAVE_LINALG_SCALAPACK 1547 integer :: ierr,istwf_k,tbloc 1548 logical :: want_eigenvectors 1549 type(matrix_scalapack) :: Slk_mat,Slk_vec 1550 type(processor_scalapack) :: Slk_processor 1551 #endif 1552 1553 !************************************************************************ 1554 1555 use_scalapack=.FALSE. 1556 if (PRESENT(comm)) then 1557 nprocs = xmpi_comm_size(comm) 1558 #ifdef HAVE_LINALG_SCALAPACK 1559 use_scalapack = (nprocs>1) 1560 #endif 1561 end if 1562 1563 SELECT CASE(use_scalapack) 1564 1565 CASE (.FALSE.) ! Standard LAPACK call. 1566 1567 lwork = MAX(1,2*n) 1568 1569 ABI_MALLOC(work,(lwork)) 1570 ABI_MALLOC(rwork,(7*n)) 1571 ABI_MALLOC(iwork,(5*n)) 1572 ABI_MALLOC(ifail,(n)) 1573 1574 call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info) 1575 1576 if (info < 0) then 1577 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value." 1578 MSG_ERROR(msg) 1579 end if 1580 1581 if (info > 0) then 1582 write(msg,'(2a,i0,a)')& 1583 & "ZHEEVX: the algorithm failed to converge; ",ch10,& 1584 & info,"eigenvectors failed to converge. " 1585 MSG_ERROR(msg) 1586 end if 1587 1588 ABI_FREE(iwork) 1589 ABI_FREE(ifail) 1590 ABI_FREE(rwork) 1591 ABI_FREE(work) 1592 1593 RETURN 1594 1595 CASE (.TRUE.) 1596 1597 #ifdef HAVE_LINALG_SCALAPACK 1598 call init_scalapack(Slk_processor,comm) 1599 istwf_k=1 1600 1601 ! Initialize and fill Scalapack matrix from the global one. 1602 tbloc=SLK_BLOCK_SIZE 1603 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1604 1605 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 1606 call wrtout(std_out,msg,"COLL") 1607 1608 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 1609 1610 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 1611 if (want_eigenvectors) then ! Initialize the distributed vectors. 1612 call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1613 end if 1614 1615 ! Solve the problem. 1616 call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w) 1617 1618 call destruction_matrix_scalapack(Slk_mat) 1619 1620 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 1621 z = czero 1622 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 1623 call destruction_matrix_scalapack(Slk_vec) 1624 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 1625 end if 1626 1627 call end_scalapack(Slk_processor) 1628 1629 RETURN 1630 #endif 1631 1632 MSG_BUG("You should not be here!") 1633 1634 END SELECT 1635 1636 end subroutine wrap_ZHEEVX
m_abilasi/wrap_ZHEGV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZHEGV
FUNCTION
wrap_ZHEGV computes all the eigenvalues, and optionally, the eigenvectors of a complex generalized
Hermitian-definite eigenproblem, of the form
A*x=(lambda)*B*x (1),
A*Bx=(lambda)*x, (2), or
B*A*x=(lambda)*x (3).
Here A and B are assumed to be Hermitian and B is also positive definite.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= "N": Compute eigenvalues only;
= "V": Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= "U": Upper triangle of A and B are stored;
= "L": Lower triangle of A and B are stored.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = "U", the leading N-by-N upper triangular part of A
<S-F1>contains the upper triangular part of the matrix A.
If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.
On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
The eigenvectors are normalized as follows: if ITYPE = 1 or 2,
Z**H*B*Z = I; if ITYPE = 3, Z**H*inv(B)*Z = I.
If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
(if UPLO="L") of A, including the diagonal, is destroyed.
B (input/output) COMPLEX*16 array, dimension (LDB, N)
On entry, the Hermitian positive definite matrix B.
If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
factor U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
1017 subroutine wrap_ZHEGV(itype,jobz,uplo,n,a,b,w,comm) 1018 1019 1020 !This section has been created automatically by the script Abilint (TD). 1021 !Do not modify the following lines by hand. 1022 #undef ABI_FUNC 1023 #define ABI_FUNC 'wrap_ZHEGV' 1024 use interfaces_14_hidewrite 1025 !End of the abilint section 1026 1027 implicit none 1028 1029 !Arguments ------------------------------------ 1030 !scalars 1031 integer,intent(in) :: n,itype 1032 integer,optional,intent(in) :: comm 1033 character(len=*),intent(in) :: jobz,uplo 1034 !arrays 1035 complex(dpc),intent(inout) :: a(n,n),b(n,n) 1036 real(dp),intent(out) :: w(n) 1037 1038 !Local variables ------------------------------ 1039 !scalars 1040 integer :: lwork,info,nprocs,ii 1041 logical :: use_scalapack 1042 character(len=500) :: msg 1043 !arrays 1044 real(dp),allocatable :: rwork(:) 1045 complex(dpc),allocatable :: work(:) 1046 #ifdef HAVE_LINALG_SCALAPACK 1047 integer :: ierr,istwf_k,tbloc 1048 type(matrix_scalapack) :: Slk_matA,Slk_matB 1049 type(processor_scalapack) :: Slk_processor 1050 #endif 1051 !************************************************************************ 1052 1053 use_scalapack=.FALSE. 1054 if (PRESENT(comm)) then 1055 nprocs = xmpi_comm_size(comm) 1056 #ifdef HAVE_LINALG_SCALAPACK 1057 use_scalapack = (nprocs>1) 1058 #endif 1059 end if 1060 1061 SELECT CASE(use_scalapack) 1062 1063 CASE (.FALSE.) 1064 1065 lwork = MAX(1,2*n-1) 1066 1067 ABI_MALLOC(work, (lwork)) 1068 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 1069 1070 call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work,lwork,rwork,info) 1071 1072 if (info < 0) then 1073 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value." 1074 MSG_ERROR(msg) 1075 end if 1076 1077 if (info > 0) then 1078 if (info<= n) then 1079 write(msg,'(2a,i0,a)')& 1080 & "ZHEGV failed to converge: ",ch10,& 1081 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 1082 else 1083 ii = info -n 1084 write(msg,'(3a,i0,3a)')& 1085 & "ZHEGV failed to converge: ",ch10,& 1086 & "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 1087 & "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 1088 end if 1089 MSG_ERROR(msg) 1090 end if 1091 1092 ABI_FREE(rwork) 1093 ABI_FREE(work) 1094 1095 RETURN 1096 1097 CASE (.TRUE.) 1098 1099 #ifdef HAVE_LINALG_SCALAPACK 1100 call init_scalapack(Slk_processor,comm) 1101 istwf_k=1 1102 1103 ! Initialize and fill Scalapack matrix from the global one. 1104 tbloc=SLK_BLOCK_SIZE 1105 1106 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 1107 call wrtout(std_out,msg,"COLL") 1108 1109 call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1110 call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 1111 1112 call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 1113 call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 1114 1115 ! Solve the problem with scaLAPACK. 1116 MSG_ERROR("slk_pZHEGV not yet coded") 1117 ! TODO 1118 !% call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w) 1119 1120 call destruction_matrix_scalapack(Slk_matB) 1121 1122 if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 1123 a = czero 1124 call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node. 1125 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 1126 end if 1127 1128 call destruction_matrix_scalapack(Slk_matA) 1129 1130 call end_scalapack(Slk_processor) 1131 1132 RETURN 1133 #endif 1134 1135 MSG_BUG("You should not be here!") 1136 1137 END SELECT 1138 1139 end subroutine wrap_ZHEGV
m_abilasi/wrap_ZHEGVX [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZHEGVX
FUNCTION
wrap_ZHEGVX - compute selected eigenvalues, and optionally, eigenvectors of a complex generalized Hermitian-definite eigenproblem, of the form A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.
INPUTS
ITYPE (input) INTEGER Specifies the problem type to be solved:
= 1: A*x = (lambda)*B*x
= 2: A*B*x = (lambda)*x
= 3: B*A*x = (lambda)*x
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval (VL,VU]
will be found.
= 'I': the IL-th through IU-th eigenvalues will be found.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrices A and B. N >= 0.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).
VL (input) REAL(DP)
VU (input) REAL(DP)
If RANGE='V', the lower and upper bounds of the interval to
be searched for eigenvalues. VL < VU.
Not referenced if RANGE = 'A' or 'I'.
IL (input) INTEGER
IU (input) INTEGER
If RANGE='I', the indices (in ascending order) of the
smallest and largest eigenvalues to be returned.
1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
Not referenced if RANGE = 'A' or 'V'.
ABSTOL (input) REAL(DP)
The absolute error tolerance for the eigenvalues.
An approximate eigenvalue is accepted as converged
when it is determined to lie in an interval [a,b]
of width less than or equal to
ABSTOL + EPS * max( |a|,|b| ) ,
where EPS is the machine precision. If ABSTOL is less than
or equal to zero, then EPS*|T| will be used in its place,
where |T| is the 1-norm of the tridiagonal matrix obtained
by reducing A to tridiagonal form.
Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero.
If this routine returns with INFO>0, indicating that some
eigenvectors did not converge, try setting ABSTOL to
2*DLAMCH('S').
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called.
OUTPUT
M (output) INTEGER
The total number of eigenvalues found. 0 <= M <= N.
If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
W (output) REAL(DP) array, dimension (N)
On normal exit, the first M elements contain the selected
eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, max(1,M))
If JOBZ = 'V', then if INFO = 0, the first M columns of Z
contain the orthonormal eigenvectors of the matrix A
corresponding to the selected eigenvalues, with the i-th
column of Z holding the eigenvector associated with W(i).
The eigenvectors are normalized as follows: if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
If an eigenvector fails to converge, then that column of Z
contains the latest approximation to the eigenvector, and the
index of the eigenvector is returned in IFAIL.
If JOBZ = 'N', then Z is not referenced.
Note: the user must ensure that at least max(1,M) columns are
supplied in the array Z; if RANGE = 'V', the exact value of M
is not known in advance and an upper bound must be used.
See also SIDE EFFECTS
SIDE EFFECTS
A (input/output) COMPLEX(DPC) array, dimension (N, N)
On entry, the Hermitian matrix A. If UPLO = 'U', the
leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = "L",
the leading N-by-N lower triangular part of A contains
the lower triangular part of the matrix A.
On exit, the lower triangle (if UPLO="L") or the upper
triangle (if UPLO="U") of A, including the diagonal, is
destroyed.
B (input/output) COMPLEX(DPC) array, dimension (LDB, N)
On entry, the Hermitian matrix B. If UPLO = "U", the leading N-by-N upper triangular part
of B contains the upper triangular part of the matrix B.
If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.
On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
U or L from the Cholesky factorization B = U**H*U or B = L*L**H.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
2060 subroutine wrap_ZHEGVX(itype,jobz,range,uplo,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm) 2061 2062 2063 !This section has been created automatically by the script Abilint (TD). 2064 !Do not modify the following lines by hand. 2065 #undef ABI_FUNC 2066 #define ABI_FUNC 'wrap_ZHEGVX' 2067 use interfaces_14_hidewrite 2068 !End of the abilint section 2069 2070 implicit none 2071 2072 !Arguments ------------------------------------ 2073 !scalars 2074 integer,intent(in) :: il,iu,ldz,n,itype 2075 integer,optional,intent(in) :: comm 2076 integer,intent(inout) :: m 2077 real(dp),intent(in) :: abstol,vl,vu 2078 character(len=*),intent(in) :: jobz,range,uplo 2079 !arrays 2080 real(dp),intent(out) :: w(n) 2081 !complex(dpc),intent(out) :: z(ldz,n) 2082 complex(dpc),intent(out) :: z(ldz,m) 2083 complex(dpc),intent(inout) :: a(n,n),b(n,n) 2084 2085 !Local variables ------------------------------ 2086 !scalars 2087 integer :: lwork,info,nprocs,ii 2088 logical :: use_scalapack 2089 character(len=500) :: msg 2090 !arrays 2091 integer,allocatable :: ifail(:),iwork(:) 2092 real(dp),allocatable :: rwork(:) 2093 complex(dpc),allocatable :: work(:) 2094 #ifdef HAVE_LINALG_SCALAPACK 2095 integer :: ierr,istwf_k,tbloc 2096 logical :: want_eigenvectors 2097 type(matrix_scalapack) :: Slk_matA,Slk_matB,Slk_vec 2098 type(processor_scalapack) :: Slk_processor 2099 #endif 2100 2101 !************************************************************************ 2102 2103 use_scalapack=.FALSE. 2104 if (PRESENT(comm)) then 2105 nprocs = xmpi_comm_size(comm) 2106 #ifdef HAVE_LINALG_SCALAPACK 2107 use_scalapack = (nprocs>1) 2108 #endif 2109 end if 2110 2111 SELECT CASE(use_scalapack) 2112 2113 CASE (.FALSE.) ! Standard LAPACK call. 2114 2115 lwork = MAX(1,2*n) 2116 2117 ABI_MALLOC(work,(lwork)) 2118 ABI_MALLOC(rwork,(7*n)) 2119 ABI_MALLOC(iwork,(5*n)) 2120 ABI_MALLOC(ifail,(n)) 2121 2122 call ZHEGVX(itype,jobz,range,uplo,n,a,n,b,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info) 2123 2124 if (info < 0) then 2125 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGVX had an illegal value." 2126 MSG_ERROR(msg) 2127 end if 2128 2129 if (info > 0) then 2130 if (info<= n) then 2131 write(msg,'(a,i0,a)')& 2132 & "ZHEGVX failed to converge: ",info," eigenvectors failed to converge. " 2133 else 2134 ii = info -n 2135 write(msg,'(3a,i0,3a)')& 2136 & "ZHEEVX failed to converge: ",ch10,& 2137 & "The leading minor of order ",ii," of B is not positive definite. ",ch10,& 2138 & "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed." 2139 end if 2140 MSG_ERROR(msg) 2141 end if 2142 2143 ABI_FREE(iwork) 2144 ABI_FREE(ifail) 2145 ABI_FREE(rwork) 2146 ABI_FREE(work) 2147 2148 RETURN 2149 2150 CASE (.TRUE.) 2151 2152 #ifdef HAVE_LINALG_SCALAPACK 2153 call init_scalapack(Slk_processor,comm) 2154 istwf_k=1 2155 2156 ! Initialize and fill Scalapack matrix from the global one. 2157 tbloc=SLK_BLOCK_SIZE 2158 2159 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 2160 call wrtout(std_out,msg,"COLL") 2161 2162 call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2163 call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a) 2164 2165 call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2166 call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b) 2167 2168 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 2169 if (want_eigenvectors) then ! Initialize the distributed vectors. 2170 call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 2171 end if 2172 2173 ! Solve the problem. 2174 MSG_ERROR("slk_pZHEGVX not coded yet") 2175 ! TODO write the scaLAPACK wrapper. 2176 !call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w) 2177 2178 call destruction_matrix_scalapack(Slk_matA) 2179 call destruction_matrix_scalapack(Slk_matB) 2180 2181 if (want_eigenvectors) then ! A is overwritten with the eigenvectors 2182 z = czero 2183 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 2184 call destruction_matrix_scalapack(Slk_vec) 2185 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 2186 end if 2187 2188 call end_scalapack(Slk_processor) 2189 2190 RETURN 2191 #endif 2192 2193 MSG_BUG("You should not be here!") 2194 2195 END SELECT 2196 2197 end subroutine wrap_ZHEGVX
m_abilasi/wrap_ZHPEV [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
wrap_ZHPEV
FUNCTION
wrap_ZHPEV computes all the eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix in packed storage. Scalapack version is not available. [PRIVATE].
INPUTS
JOBZ (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
UPLO (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
N (input) INTEGER
The order of the matrix A. N >= 0.
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1, and if
JOBZ = 'V', LDZ >= max(1,N).
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1,
in this case the sequential LAPACK routine is called. Note that scalapack does not provide native
support for packed symmetric matrices. Threfore we have to distribute the full matrix among the nodes.
in order to perform the calculation in parallel.
OUTPUT
W (output) REAL(DP) array, dimension (N)
If INFO = 0, the eigenvalues in ascending order.
Z (output) COMPLEX(DPC) array, dimension (LDZ, N)
If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
eigenvectors of the matrix A, with the i-th column of Z
holding the eigenvector associated with W(i).
If JOBZ = 'N', then Z is not referenced.
See also SIDE EFFECTS
SIDE EFFECTS
AP (input/output) COMPLEX(DPC) array, dimension (N*(N+1)/2)
On entry, the upper or lower triangle of the Hermitian matrix
A, packed columnwise in a linear array. The j-th column of A
is stored in the array AP as follows:
if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
On exit, AP is overwritten by values generated during the
reduction to tridiagonal form. If UPLO = 'U', the diagonal
and first superdiagonal of the tridiagonal matrix T overwrite
the corresponding elements of A, and if UPLO = 'L', the
diagonal and first subdiagonal of T overwrite the
corresponding elements of A. Unchanged if ScaLAPACK is used.
PARENTS
CHILDREN
cwtime,xginv
SOURCE
830 subroutine wrap_ZHPEV(jobz,uplo,n,ap,w,z,ldz,comm) 831 832 833 !This section has been created automatically by the script Abilint (TD). 834 !Do not modify the following lines by hand. 835 #undef ABI_FUNC 836 #define ABI_FUNC 'wrap_ZHPEV' 837 use interfaces_14_hidewrite 838 !End of the abilint section 839 840 implicit none 841 842 !Arguments ------------------------------------ 843 !scalars 844 integer,intent(in) :: n,ldz 845 integer,optional,intent(in) :: comm 846 character(len=*),intent(in) :: jobz,uplo 847 !arrays 848 real(dp),intent(out) :: w(n) 849 complex(dpc),intent(inout) :: ap(n*(n+1)/2) 850 complex(dpc),intent(out) :: z(ldz,n) 851 852 !Local variables ------------------------------ 853 !scalars 854 integer :: info,nprocs 855 logical :: use_scalapack 856 character(len=500) :: msg 857 !arrays 858 real(dp),allocatable :: rwork(:) 859 complex(dpc),allocatable :: work(:) 860 #ifdef HAVE_LINALG_SCALAPACK 861 integer :: ierr,istwf_k,tbloc 862 logical :: want_eigenvectors 863 type(matrix_scalapack) :: Slk_mat,Slk_vec 864 type(processor_scalapack) :: Slk_processor 865 #endif 866 867 !************************************************************************ 868 869 use_scalapack=.FALSE. 870 if (PRESENT(comm)) then 871 nprocs = xmpi_comm_size(comm) 872 #ifdef HAVE_LINALG_SCALAPACK 873 use_scalapack = (nprocs>1) 874 #endif 875 end if 876 877 SELECT CASE(use_scalapack) 878 879 CASE (.FALSE.) 880 881 ABI_MALLOC(work, (MAX(1,2*n-1))) 882 ABI_MALLOC(rwork, (MAX(1,3*n-2))) 883 884 call ZHPEV(jobz,uplo,n,ap,w,z,ldz,work,rwork,info) 885 886 if (info < 0) then 887 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHPEV had an illegal value." 888 MSG_ERROR(msg) 889 end if 890 891 if (info > 0) then 892 write(msg,'(2a,i0,a)')& 893 & "ZHPEV: the algorithm failed to converge; ",ch10,& 894 & info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. " 895 MSG_ERROR(msg) 896 end if 897 898 ABI_FREE(rwork) 899 ABI_FREE(work) 900 RETURN 901 902 CASE (.TRUE.) 903 904 #ifdef HAVE_LINALG_SCALAPACK 905 call init_scalapack(Slk_processor,comm) 906 istwf_k=1 907 908 ! Initialize and fill Scalapack matrix from the global one. 909 tbloc=SLK_BLOCK_SIZE 910 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 911 912 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 913 call wrtout(std_out,msg,"COLL") 914 915 call slk_matrix_from_global_dpc_1Dp(Slk_mat,uplo,ap) 916 917 want_eigenvectors = firstchar(jobz,(/"V","v"/)) 918 if (want_eigenvectors) then ! Initialize the distributed vectors. 919 call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 920 end if 921 922 ! Solve the problem with scaLAPACK. 923 call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w) 924 925 call destruction_matrix_scalapack(Slk_mat) 926 927 if (want_eigenvectors) then ! Collect the eigenvectors. 928 z = zero 929 call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node. 930 call destruction_matrix_scalapack(Slk_vec) 931 call xmpi_sum(z,comm,ierr) ! Fill the remaing entries of the global matrix 932 end if 933 934 call end_scalapack(Slk_processor) 935 936 RETURN 937 #endif 938 939 MSG_BUG("You should not be here!") 940 941 END SELECT 942 943 end subroutine wrap_ZHPEV
m_abilasi/zginv [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
zginv
FUNCTION
Invert a general matrix of complex elements in double precision. ZGETRF computes an LU factorization of a general N-by-N matrix A using partial pivoting with row interchanges. ZGETRI computes the inverse of a matrix using the LU factorization computed by ZGETRF.
INPUTS
n=size of complex matrix a
a=matrix of complex elements
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1.
In this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)= array of complex elements, input, inverted at output
PARENTS
CHILDREN
cwtime,xginv
SOURCE
3030 subroutine zginv(a,n,comm) 3031 3032 3033 !This section has been created automatically by the script Abilint (TD). 3034 !Do not modify the following lines by hand. 3035 #undef ABI_FUNC 3036 #define ABI_FUNC 'zginv' 3037 use interfaces_14_hidewrite 3038 !End of the abilint section 3039 3040 implicit none 3041 3042 !Arguments ------------------------------------ 3043 !scalars 3044 integer,intent(in) :: n 3045 integer,optional,intent(in) :: comm 3046 !arrays 3047 complex(dpc),intent(inout) :: a(n,n) 3048 3049 !Local variables------------------------------- 3050 !scalars 3051 integer :: lwork,info,nprocs 3052 logical :: use_scalapack 3053 character(len=500) :: msg 3054 !arrays 3055 integer,allocatable :: ipiv(:) 3056 complex(dpc),allocatable :: work(:) 3057 #ifdef HAVE_LINALG_SCALAPACK 3058 integer :: istwf_k,tbloc,ierr 3059 type(matrix_scalapack) :: Slk_mat 3060 type(processor_scalapack) :: Slk_processor 3061 #endif 3062 3063 ! ************************************************************************* 3064 3065 use_scalapack=.FALSE. 3066 if (PRESENT(comm)) then 3067 nprocs = xmpi_comm_size(comm) 3068 #ifdef HAVE_LINALG_SCALAPACK 3069 use_scalapack = (nprocs>1) 3070 #endif 3071 end if 3072 3073 SELECT CASE(use_scalapack) 3074 3075 CASE (.FALSE.) 3076 3077 ABI_MALLOC(ipiv,(n)) 3078 call ZGETRF(n,n,a,n,ipiv,info) ! P* L* U Factorization. 3079 3080 if (info < 0) then 3081 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRF had an illegal value." 3082 MSG_ERROR(msg) 3083 end if 3084 3085 if (info > 0) then 3086 write(msg,'(3a,i0,4a)')& 3087 & "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,& 3088 & "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,& 3089 & "The factorization has been completed but the factor U is exactly singular.",ch10,& 3090 & "Division by zero will occur if it is used to solve a system of equations." 3091 MSG_ERROR(msg) 3092 end if 3093 3094 lwork=MAX(1,n) 3095 3096 ABI_MALLOC(work,(lwork)) 3097 3098 call ZGETRI(n,a,n,ipiv,work,lwork,info) ! Invert U and then compute inv(A) 3099 3100 if (info < 0) then 3101 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRI had an illegal value." 3102 MSG_ERROR(msg) 3103 end if 3104 3105 if (info > 0) then 3106 write(msg,'(3a,i0,a)')& 3107 & "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,& 3108 & "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed." 3109 MSG_ERROR(msg) 3110 end if 3111 3112 ABI_FREE(ipiv) 3113 ABI_FREE(work) 3114 3115 RETURN 3116 3117 CASE (.TRUE.) 3118 3119 #ifdef HAVE_LINALG_SCALAPACK 3120 call init_scalapack(Slk_processor,comm) 3121 istwf_k=1 3122 3123 ! Initialize and fill Scalapack matrix from the global one. 3124 tbloc=SLK_BLOCK_SIZE 3125 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 3126 3127 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 3128 call wrtout(std_out,msg,"COLL") 3129 3130 call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a) 3131 3132 ! Perform the calculation with scaLAPACK. 3133 call slk_zinvert(Slk_mat) 3134 3135 ! Reconstruct the global matrix from the distributed one. 3136 a = czero 3137 call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a) ! Fill the entries calculated by this node. 3138 call destruction_matrix_scalapack(Slk_mat) 3139 3140 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 3141 call end_scalapack(Slk_processor) 3142 3143 RETURN 3144 #endif 3145 3146 MSG_BUG("You should not be here!") 3147 3148 END SELECT 3149 3150 end subroutine zginv
m_abilasi/zhpd_invert [ Functions ]
[ Top ] [ m_abilasi ] [ Functions ]
NAME
zhpd_invert
FUNCTION
Invert a Hermitian positive definite matrix of complex elements in double precision.
INPUTS
uplo= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
n=size of complex matrix a
a=matrix of complex elements
[comm]=MPI communicator for ScaLAPACK inversion. Only available if the code has been compiled with Scalapack support.
To avoid wasting CPU time the scalapack initialization is avoided if the number of processors in 1.
In this case the sequential LAPACK routine is called.
SIDE EFFECTS
a(n,n)=
On entry, the Hermitian matrix A. If UPLO = 'U', the leading
N-by-N upper triangular part of A contains the upper
triangular part of the matrix A, and the strictly lower
triangular part of A is not referenced. If UPLO = 'L', the
leading N-by-N lower triangular part of A contains the lower
triangular part of the matrix A, and the strictly upper
triangular part of A is not referenced.
On exit, the upper or lower triangle of the (Hermitian) inverse of A
PARENTS
CHILDREN
cwtime,xginv
SOURCE
3189 subroutine zhpd_invert(uplo,a,n,comm) 3190 3191 3192 !This section has been created automatically by the script Abilint (TD). 3193 !Do not modify the following lines by hand. 3194 #undef ABI_FUNC 3195 #define ABI_FUNC 'zhpd_invert' 3196 use interfaces_14_hidewrite 3197 !End of the abilint section 3198 3199 implicit none 3200 3201 !Arguments ------------------------------------ 3202 !scalars 3203 character(len=*),intent(in) :: uplo 3204 integer,intent(in) :: n 3205 integer,optional,intent(in) :: comm 3206 !arrays 3207 complex(dpc),intent(inout) :: a(n,n) 3208 3209 !Local variables------------------------------- 3210 !scalars 3211 integer :: info,nprocs 3212 logical :: use_scalapack 3213 character(len=500) :: msg 3214 !arrays 3215 #ifdef HAVE_LINALG_SCALAPACK 3216 integer :: istwf_k,tbloc,ierr 3217 type(matrix_scalapack) :: Slk_mat 3218 type(processor_scalapack) :: Slk_processor 3219 #endif 3220 3221 ! ************************************************************************* 3222 3223 use_scalapack=.FALSE. 3224 if (PRESENT(comm)) then 3225 nprocs = xmpi_comm_size(comm) 3226 #ifdef HAVE_LINALG_SCALAPACK 3227 use_scalapack = (nprocs>1) 3228 #endif 3229 end if 3230 3231 SELECT CASE(use_scalapack) 3232 3233 CASE (.FALSE.) 3234 ! * ZPOTRF computes the Cholesky factorization of a complex Hermitian positive definite. 3235 ! * A = U**H * U, if UPLO = 'U', or 3236 ! * A = L * L**H, if UPLO = 'L', 3237 call ZPOTRF(uplo,n,a,n,info) 3238 3239 if (info < 0) then 3240 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRF had an illegal value." 3241 MSG_ERROR(msg) 3242 end if 3243 3244 if (info > 0) then 3245 write(msg,'(a,i0,3a)')& 3246 & "The leading minor of order ",info," is not positive definite, ",ch10,& 3247 & "and the factorization could not be completed." 3248 MSG_ERROR(msg) 3249 end if 3250 ! 3251 ! * ZPOTRI computes the inverse of a complex Hermitian positive definite 3252 ! * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 3253 ! * computed by ZPOTRF. 3254 ! * On exit, the upper or lower triangle of the (Hermitian) 3255 ! * inverse of A, overwriting the input factor U or L. 3256 call ZPOTRI(uplo,n,a,n,info) 3257 3258 if (info < 0) then 3259 write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value." 3260 MSG_ERROR(msg) 3261 end if 3262 3263 if (info > 0) then 3264 write(msg,'(a,2(1x,i0),a)')& 3265 & "The ( ",info,info,")element of the factor U or L is zero, and the inverse could not be computed." 3266 MSG_ERROR(msg) 3267 end if 3268 3269 RETURN 3270 3271 CASE (.TRUE.) 3272 3273 #ifdef HAVE_LINALG_SCALAPACK 3274 call init_scalapack(Slk_processor,comm) 3275 istwf_k=1 3276 3277 ! Initialize and fill Scalapack matrix from the global one. 3278 tbloc=SLK_BLOCK_SIZE 3279 call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc) 3280 3281 write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc 3282 call wrtout(std_out,msg,"COLL") 3283 3284 call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a) 3285 3286 ! Perform the calculation with scaLAPACK. 3287 call slk_zdhp_invert(Slk_mat,uplo) 3288 3289 ! Reconstruct the global matrix from the distributed one. 3290 a = czero 3291 call slk_matrix_to_global_dpc_2D(Slk_mat,uplo,a) ! Fill the entries calculated by this node. 3292 call destruction_matrix_scalapack(Slk_mat) 3293 3294 call xmpi_sum(a,comm,ierr) ! Fill the remaing entries of the global matrix 3295 call end_scalapack(Slk_processor) 3296 3297 RETURN 3298 #endif 3299 MSG_BUG("You should not be here!") 3300 END SELECT 3301 3302 end subroutine zhpd_invert