## ABINIT/jacobi [ Functions ]

[ Top ] [ 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

```3602 subroutine jacobi(a,n,np,d,v,nrot)
3603
3604
3605 !This section has been created automatically by the script Abilint (TD).
3606 !Do not modify the following lines by hand.
3607 #undef ABI_FUNC
3608 #define ABI_FUNC 'jacobi'
3609 !End of the abilint section
3610
3611  implicit none
3612 !Arguments
3613  integer :: n,np,nrot
3614  real*8 :: a(np,np),d(np),v(np,np)
3615 !Local variables
3616  integer, parameter :: NMAX=500
3617  integer i,ip,iq,j
3618  real*8 c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
3619  do ip=1,n
3620    do iq=1,n
3621      v(ip,iq)=0.
3622    enddo
3623    v(ip,ip)=1.
3624  enddo
3625  do ip=1,n
3626    b(ip)=a(ip,ip)
3627    d(ip)=b(ip)
3628    z(ip)=0.
3629  enddo
3630  nrot=0
3631  do i=1,50
3632    sm=0.
3633    do ip=1,n-1
3634      do iq=ip+1,n
3635        sm=sm+abs(a(ip,iq))
3636      enddo
3637    enddo
3638    if(sm.eq.0.)return
3639    if(i.lt.4)then
3640      tresh=0.2*sm/n**2
3641    else
3642      tresh=0.
3643    endif
3644    do ip=1,n-1
3645      do iq=ip+1,n
3646        g=100.*abs(a(ip,iq))
3647        if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) &
3648 &          .and.(abs(d(iq))+g.eq.abs(d(iq))))then
3649             a(ip,iq)=0.
3650        else if(abs(a(ip,iq)).gt.tresh)then
3651          h=d(iq)-d(ip)
3652          if(abs(h)+g.eq.abs(h))then
3653            t=a(ip,iq)/h
3654          else
3655            theta=0.5*h/a(ip,iq)
3656            t=1./(abs(theta)+sqrt(1.+theta**2))
3657            if(theta.lt.0.)t=-t
3658          endif
3659          c=1./sqrt(1+t**2)
3660          s=t*c
3661          tau=s/(1.+c)
3662          h=t*a(ip,iq)
3663          z(ip)=z(ip)-h
3664          z(iq)=z(iq)+h
3665          d(ip)=d(ip)-h
3666          d(iq)=d(iq)+h
3667          a(ip,iq)=0.
3668          do j=1,ip-1
3669            g=a(j,ip)
3670            h=a(j,iq)
3671            a(j,ip)=g-s*(h+g*tau)
3672            a(j,iq)=h+s*(g-h*tau)
3673          enddo
3674          do j=ip+1,iq-1
3675            g=a(ip,j)
3676            h=a(j,iq)
3677            a(ip,j)=g-s*(h+g*tau)
3678            a(j,iq)=h+s*(g-h*tau)
3679          enddo
3680          do j=iq+1,n
3681            g=a(ip,j)
3682            h=a(iq,j)
3683            a(ip,j)=g-s*(h+g*tau)
3684            a(iq,j)=h+s*(g-h*tau)
3685          enddo
3686          do j=1,n
3687            g=v(j,ip)
3688            h=v(j,iq)
3689            v(j,ip)=g-s*(h+g*tau)
3690            v(j,iq)=h+s*(g-h*tau)
3691          enddo
3692          nrot=nrot+1
3693        endif
3694      enddo
3695    enddo
3696    do ip=1,n
3697      b(ip)=b(ip)+z(ip)
3698      d(ip)=b(ip)
3699      z(ip)=0.
3700    enddo
3701  enddo
3702  write(std_out,*) 'too many iterations in jacobi'
3703
3704 end subroutine jacobi
```

## ABINIT/m_hide_lapack [ Modules ]

[ Top ] [ Modules ]

NAME

```  m_hide_lapack
```

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 (C) 1992-2018 ABINIT group (MG, GMR, XG)
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_hide_lapack
52
53  use defs_basis
54  use m_abicore
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_hide_lapack/cginv [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```2823 subroutine cginv(a,n,comm)
2824
2825
2826 !This section has been created automatically by the script Abilint (TD).
2827 !Do not modify the following lines by hand.
2828 #undef ABI_FUNC
2829 #define ABI_FUNC 'cginv'
2830 !End of the abilint section
2831
2832  implicit none
2833
2834 !Arguments ------------------------------------
2835 !scalars
2836  integer,intent(in) :: n
2837  integer,optional,intent(in) :: comm
2838 !arrays
2839  complex(spc),intent(inout) :: a(n,n)
2840
2841 !Local variables-------------------------------
2842 !scalars
2843  integer :: lwork,info,nprocs
2844  logical :: use_scalapack
2845  character(len=500) :: msg
2846 !arrays
2847  integer,allocatable :: ipiv(:)
2848  complex(spc),allocatable :: work(:)
2849 #ifdef HAVE_LINALG_SCALAPACK
2850  integer :: ierr,istwf_k,ipiv_size,liwork,tbloc
2851  integer,allocatable :: iwork(:)
2852  type(matrix_scalapack)    :: Slk_mat
2853  type(processor_scalapack) :: Slk_processor
2854 #endif
2855
2856 ! *************************************************************************
2857
2858  use_scalapack=.FALSE.
2859  if (PRESENT(comm)) then
2860   nprocs = xmpi_comm_size(comm)
2861 #ifdef HAVE_LINALG_SCALAPACK
2862   use_scalapack = (nprocs>1)
2863 #endif
2864  end if
2865
2866  SELECT CASE(use_scalapack)
2867
2868  CASE (.FALSE.)
2869
2870   ABI_MALLOC(ipiv,(n))
2871
2872   call CGETRF(n,n,a,n,ipiv,info) ! P* L* U  Factorization.
2873
2874   if (info < 0) then
2875    write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRF had an illegal value."
2876    MSG_ERROR(msg)
2877   end if
2878
2879   if (info > 0) then
2880    write(msg,'(3a,i0,4a)')&
2881 &   "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,&
2882 &   "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,&
2883 &   "The factorization has been completed but the factor U is exactly singular.",ch10,&
2884 &   "Division by zero will occur if it is used to solve a system of equations."
2885    MSG_ERROR(msg)
2886   end if
2887
2888   lwork=MAX(1,n)
2889
2890   ABI_MALLOC(work,(lwork))
2891
2892   call CGETRI(n,a,n,ipiv,work,lwork,info) ! Inverts U and the computes inv(A)
2893
2894   if (info < 0) then
2895    write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGETRI had an illegal value."
2896    MSG_ERROR(msg)
2897   end if
2898
2899   if (info > 0) then
2900    write(msg,'(3a,i0,a)')&
2901 &   "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,&
2902     "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed."
2903    MSG_ERROR(msg)
2904   end if
2905
2906   ABI_FREE(ipiv)
2907   ABI_FREE(work)
2908
2909   RETURN
2910
2911  CASE (.TRUE.)
2912
2913 #if 0
2914 ! FIXME matrix_scalapack does not have a single precision complex buffer
2915
2916 #ifdef HAVE_LINALG_SCALAPACK
2917   call init_scalapack(Slk_processor,comm)
2918   istwf_k=1
2919
2920   ! Initialize and fill Scalapack matrix from the global one.
2921   tbloc=SLK_BLOCK_SIZE
2922   call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2923
2924   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
2925   call wrtout(std_out,msg,"COLL")
2926
2927   ! IMPORTANT NOTE: PZGETRF requires square block decomposition i.e.,  MB_A = NB_A.
2928   if ( Slk_mat%descript%tab(MB_)/=Slk_mat%descript%tab(NB_) ) then
2929    msg ="PZGETRF requires square block decomposition i.e.,  MB_A = NB_A."
2930    MSG_ERROR(msg)
2931   end if
2932
2933   !!call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
2934
2935   ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_)
2936   ABI_MALLOC(ipiv,(ipiv_size))
2937
2938   call PCGETRF(Slk_mat%sizeb_global(1),Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx_sp,&
2939 &   1,1,Slk_mat%descript%tab,ipiv,info) ! P * L * U  Factorization.
2940
2941   if (info/=0) then
2942    write(msg,'(a,i0)')"PCGETRF returned info= ",info
2943    MSG_ERROR(msg)
2944   end if
2945
2946   ! Get optimal size of workspace for PCGETRI.
2947   lwork=-1; liwork=-1
2948   ABI_MALLOC(work,(1))
2949   ABI_MALLOC(iwork,(1))
2950
2951   call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,&
2952 &  work,lwork,iwork,liwork,info)
2953
2954   ABI_CHECK(info==0,"PZGETRI: Error during compuation of workspace size")
2955
2956   lwork = NINT(DBLE(work(1))); liwork=iwork(1)
2957   ABI_FREE(work)
2958   ABI_FREE(iwork)
2959
2960   ! Solve the problem.
2961   ABI_MALLOC(work,(lwork))
2962   ABI_MALLOC(iwork,(liwork))
2963
2964   call PCGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx_sp,1,1,Slk_mat%descript%tab,ipiv,&
2965 &  work,lwork,iwork,liwork,info)
2966
2967   if (info/=0) then
2968    write(msg,'(a,i0)')"PZGETRI returned info= ",info
2969    MSG_ERROR(msg)
2970   end if
2971
2972   ABI_FREE(work)
2973   ABI_FREE(iwork)
2974   ABI_FREE(ipiv)
2975
2976   ! Reconstruct the global matrix from the distributed one.
2977   a = czero
2978   !! call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a)  ! Fill the entries calculated by this node.
2979   call destruction_matrix_scalapack(Slk_mat)
2980
2981   call xmpi_sum(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
2982   call end_scalapack(Slk_processor)
2983
2984   RETURN
2985 #endif
2986
2987 #endif
2988
2989   MSG_BUG("You should not be here!")
2990
2991  END SELECT
2992
2993 end subroutine cginv
```

## m_hide_lapack/dzegdi [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3918 subroutine dzgedi(a,lda,n,ipvt,det,work,job)
3919
3920
3921 !This section has been created automatically by the script Abilint (TD).
3922 !Do not modify the following lines by hand.
3923 #undef ABI_FUNC
3924 #define ABI_FUNC 'dzgedi'
3925 !End of the abilint section
3926
3927       implicit none
3928
3929       integer :: lda,n,ipvt(n),job
3930       real*8 :: a(2,lda,n),det(2,2),work(2,n)
3931 !
3932 !     zgedi computes the determinant and inverse of a matrix
3933 !     using the factors computed by zgeco or zgefa.
3934 !
3935 !     on entry
3936 !
3937 !        a       complex*16(lda, n)
3938 !                the output from zgeco or zgefa.
3939 !
3940 !        lda     integer
3941 !                the leading dimension of the array  a .
3942 !
3943 !        n       integer
3944 !                the order of the matrix  a .
3945 !
3946 !        ipvt    integer(n)
3947 !                the pivot vector from zgeco or zgefa.
3948 !
3949 !        work    complex*16(n)
3950 !                work vector.  contents destroyed.
3951 !
3952 !        job     integer
3953 !                = 11   both determinant and inverse.
3954 !                = 01   inverse only.
3955 !                = 10   determinant only.
3956 !
3957 !     on return
3958 !
3959 !        a       inverse of original matrix if requested.
3960 !                otherwise unchanged.
3961 !
3962 !        det     complex*16(2)
3963 !                determinant of original matrix if requested.
3964 !                otherwise not referenced.
3965 !                determinant = det(1) * 10.0**det(2)
3966 !                with  1.0 .le. cabs1(det(1)) .lt. 10.0
3967 !                or  det(1) .eq. 0.0 .
3968 !
3969 !     error condition
3970 !
3971 !        a division by zero will occur if the input factor contains
3972 !        a zero on the diagonal and the inverse is requested.
3973 !        it will not occur if the subroutines are called correctly
3974 !        and if zgeco has set rcond .gt. 0.0 or zgefa has set
3975 !        info .eq. 0 .
3976 !
3977 !     linpack. this version dated 08/14/78 .
3978 !     cleve moler, university of new mexico, argonne national lab.
3979 !
3980 !     subroutines and functions
3981 !
3982 !     internal variables
3983 !
3984       double precision :: r(2),rk(2),rkj(2)
3985       double precision :: ten,rinv2,rabs
3986       integer :: i,j,k,kb,kp1,l,nm1
3987 !
3988 !     compute determinant
3989 !
3990       if (job/10 .eq. 0) go to 70
3991          det(1,1) = 1.0d0; det(2,1) = 0.0d0
3992          det(1,2) = 0.0d0; det(2,2) = 0.0d0
3993          ten = 10.0d0
3994          do i = 1, n
3995             if (ipvt(i) .ne. i) then
3996                 det(1,1) = -det(1,1)
3997                 det(2,1) = -det(2,1)
3998             end if
3999             r(1)=det(1,1); r(2)=det(2,1)
4000             det(1,1) = r(1)*a(1,i,i)-r(2)*a(2,i,i)
4001             det(2,1) = r(2)*a(1,i,i)+r(1)*a(2,i,i)
4002 !        ...exit
4003             rabs = abs(det(1,1))+abs(det(2,1))
4004             if (rabs .eq. 0.0d0) go to 60
4005    10       continue
4006             rabs = abs(det(1,1))+abs(det(2,1))
4007             if (rabs .ge. 1.0d0) go to 20
4008                det(1,1) = ten*det(1,1); det(2,1) = ten*det(2,1)
4009                det(1,2) = det(1,2) - 1.0d0
4010             go to 10
4011    20       continue
4012    30       continue
4013             rabs = abs(det(1,1))+abs(det(2,1))
4014             if (rabs .lt. ten) go to 40
4015                det(1,1) = det(1,1)/ten; det(2,1) = det(2,1)/ten
4016                det(1,2) = det(1,2) + 1.0d0
4017             go to 30
4018    40       continue
4019          end do
4020    60    continue
4021    70 continue
4022 !
4023 !     compute inverse(u)
4024 !
4025       if (mod(job,10) .eq. 0) go to 150
4026          do 100 k = 1, n
4027             !a(k,k) = (1.0d0,0.0d0)/a(k,k)
4028             !t = -a(k,k)
4029             !call zscal(k-1,t,a(1,k),1)
4030             rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2)
4031             a(1,k,k) =  rinv2*a(1,k,k)
4032             a(2,k,k) = -rinv2*a(2,k,k)
4033             rk(1) = -a(1,k,k); rk(2) = -a(2,k,k)
4034             do i=1,k-1
4035                r(1)=a(1,i,k)
4036                r(2)=a(2,i,k)
4037                a(1,i,k)=rk(1)*r(1)-rk(2)*r(2)
4038                a(2,i,k)=rk(1)*r(2)+rk(2)*r(1)
4039             end do
4040             kp1 = k + 1
4041             if (n .lt. kp1) go to 90
4042             do 80 j = kp1, n
4043                !t = a(k,j)
4044                !a(k,j) = (0.0d0,0.0d0)
4045                !call zaxpy(k,t,a(1,k),1,a(1,j),1)
4046                rkj(1) = a(1,k,j); rkj(2) = a(2,k,j)
4047                a(1,k,j) = 0.d0; a(2,k,j) = 0.d0
4048                do i=1,k
4049                   a(1,i,j)=rkj(1)*a(1,i,k)-rkj(2)*a(2,i,k)+a(1,i,j)
4050                   a(2,i,j)=rkj(2)*a(1,i,k)+rkj(1)*a(2,i,k)+a(2,i,j)
4051                end do
4052    80       continue
4053    90       continue
4054   100    continue
4055   do i=1,n
4056   end do
4057 !
4058 !        form inverse(u)*inverse(l)
4059 !
4060          nm1 = n - 1
4061          if (nm1 .lt. 1) go to 140
4062          do 130 kb = 1, nm1
4063             k = n - kb
4064             kp1 = k + 1
4065             do 110 i = kp1, n
4066                work(1,i) = a(1,i,k); work(2,i) = a(2,i,k)
4067                a(1,i,k) = 0.0d0; a(2,i,k) = 0.d0
4068   110       continue
4069             do 120 j = kp1, n
4070                r(1) = work(1,j); r(2) = work(2,j)
4071                !call zaxpy(n,t,a(1,j),1,a(1,k),1)
4072                do i=1,n
4073                   a(1,i,k)=r(1)*a(1,i,j)-r(2)*a(2,i,j)+a(1,i,k)
4074                   a(2,i,k)=r(2)*a(1,i,j)+r(1)*a(2,i,j)+a(2,i,k)
4075                end do
4076   120       continue
4077             l = ipvt(k)
4078             if (l .ne. k) then
4079                !call zswap(n,a(1,k),1,a(1,l),1)
4080                do i=1,n
4081                   r(1) = a(1,i,k); r(2) = a(2,i,k)
4082                   a(1,i,k) = a(1,i,l); a(2,i,k) = a(2,i,l)
4083                   a(1,i,l) = r(1); a(2,i,l) = r(2)
4084                end do
4085             end if
4086   130    continue
4087   140    continue
4088   150 continue
4089
4090 end subroutine dzgedi
```

## m_hide_lapack/dzgefa [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```4116 subroutine dzgefa(a,lda,n,ipvt,info)
4117
4118  use m_linalg_interfaces
4119
4120 !This section has been created automatically by the script Abilint (TD).
4121 !Do not modify the following lines by hand.
4122 #undef ABI_FUNC
4123 #define ABI_FUNC 'dzgefa'
4124 !End of the abilint section
4125
4126  implicit none
4127
4128 !Arguments
4129  integer :: lda,n,ipvt(n),info
4130  real*8  :: a(2,lda,n)
4131 !
4132 !     zgefa factors a complex*16 matrix by gaussian elimination.
4133 !
4134 !     dzgefa is usually called by zgeco, but it can be called
4135 !     directly with a saving in time if  rcond  is not needed.
4136 !     (time for zgeco) = (1 + 9/n)*(time for zgefa) .
4137 !
4138 !     on entry
4139 !
4140 !        a       complex*16(lda, n)
4141 !                the matrix to be factored.
4142 !
4143 !        lda     integer
4144 !                the leading dimension of the array  a .
4145 !
4146 !        n       integer
4147 !                the order of the matrix  a .
4148 !
4149 !     on return
4150 !
4151 !        a       an upper triangular matrix and the multipliers
4152 !                which were used to obtain it.
4153 !                the factorization can be written  a = l*u  where
4154 !                l  is a product of permutation and unit lower
4155 !                triangular matrices and  u  is upper triangular.
4156 !
4157 !        ipvt    integer(n)
4158 !                an integer vector of pivot indices.
4159 !
4160 !        info    integer
4161 !                = 0  normal value.
4162 !                = k  if  u(k,k) .eq. 0.0 .  this is not an error
4163 !                     condition for this subroutine, but it does
4164 !                     indicate that zgesl or zgedi will divide by zero
4165 !                     if called.  use  rcond  in zgeco for a reliable
4166 !                     indication of singularity.
4167 !
4168 !     linpack. this version dated 08/14/78 .
4169 !     cleve moler, university of new mexico, argonne national lab.
4170 !
4171 !     subroutines and functions
4172 !
4173 !     internal variables
4174 !
4175 !Local variables
4176  real*8 :: r(2),rk(2),rlj(2)
4177  real*8 :: rinv2,rmax,rabs
4178  integer :: i,j,k,kp1,l,nm1
4179
4180 !
4181 !     gaussian elimination with partial pivoting
4182 !
4183       info = 0
4184       nm1 = n - 1
4185       if (nm1 .lt. 1) go to 70
4186       do 60 k = 1, nm1
4187          kp1 = k + 1
4188 !
4189 !        find l = pivot index
4190 !
4191          !l = izamax(n-k+1,a(k,k),1) + k - 1
4192          rmax=0.d0
4193          l=0
4194          do i=k,n
4195             rabs=abs(a(1,i,k))+abs(a(2,i,k))
4196             if(rmax<=rabs) then
4197               rmax=rabs
4198               l=i
4199             end if
4200          end do
4201          ipvt(k) = l
4202 !
4203 !        zero pivot implies this column already triangularized
4204 !
4205          if (abs(a(1,l,k))+abs(a(2,l,k)) .eq. 0.0d0) go to 40
4206 !
4207 !           interchange if necessary
4208 !
4209             if (l .eq. k) go to 10
4210                r(1) = a(1,l,k); r(2) = a(2,l,k)
4211                a(1,l,k) = a(1,k,k); a(2,l,k) = a(2,k,k)
4212                a(1,k,k) = r(1); a(2,k,k) = r(2)
4213    10       continue
4214 !
4215 !           compute multipliers
4216 !
4217             rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2)
4218             rk(1) = -rinv2*a(1,k,k)
4219             rk(2) =  rinv2*a(2,k,k)
4220             !call zscal(n-k,t,a(k+1,k),1)
4221             do i=k+1,n
4222                r(1)=a(1,i,k)
4223                r(2)=a(2,i,k)
4224                a(1,i,k)=rk(1)*r(1)-rk(2)*r(2)
4225                a(2,i,k)=rk(1)*r(2)+rk(2)*r(1)
4226             end do
4227 !
4228 !           row elimination with column indexing
4229 !
4230             do j = kp1, n
4231                rlj(1) = a(1,l,j); rlj(2) = a(2,l,j)
4232                if (l .eq. k) go to 20
4233                   a(1,l,j) = a(1,k,j); a(2,l,j) = a(2,k,j)
4234                   a(1,k,j) = rlj(1); a(2,k,j) = rlj(2)
4235    20          continue
4236                !call zaxpy(n-k,t,a(1,k+1,k),1,a(1,k+1,j),1)
4237                do i=k+1,n
4238                   a(1,i,j)=rlj(1)*a(1,i,k)-rlj(2)*a(2,i,k)+a(1,i,j)
4239                   a(2,i,j)=rlj(2)*a(1,i,k)+rlj(1)*a(2,i,k)+a(2,i,j)
4240                end do
4241             end do
4242          go to 50
4243    40    continue
4244             info = k
4245    50    continue
4246    60 continue
4247    70 continue
4248       ipvt(n) = n
4249       if (abs(a(1,n,n))+abs(a(2,n,n)) .eq. 0.0d0) info = n
4250
4251 end subroutine dzgefa
```

## m_hide_lapack/lubksb [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3848 SUBROUTINE lubksb(a,n,np,indx,b)
3849
3850
3851 !This section has been created automatically by the script Abilint (TD).
3852 !Do not modify the following lines by hand.
3853 #undef ABI_FUNC
3854 #define ABI_FUNC 'lubksb'
3855 !End of the abilint section
3856
3857       IMPLICIT NONE
3858
3859       INTEGER n,np,indx(n)
3860       REAL*8 a(np,np),b(n)
3861
3862       INTEGER i,ii,j,ll
3863       REAL*8 sum
3864 !      write(std_out,*) 'ENTERING LUBKSB...'
3865 !      write(std_out,201) ((a(i,j),j=1,n),i=1,n)
3866 ! 201  FORMAT('A in lubksb ',/,3F16.8,/,3F16.8,/,3F16.8)
3867
3868       ii=0
3869       do i=1,n
3870         ll=indx(i)
3871         sum=b(ll)
3872         b(ll)=b(i)
3873         if (ii.ne.0)then
3874           do j=ii,i-1
3875             sum=sum-a(i,j)*b(j)
3876           enddo
3877         else if (sum.ne.0.) then
3878           ii=i
3879         endif
3880         b(i)=sum
3881       enddo
3882       do i=n,1,-1
3883         sum=b(i)
3884         do j=i+1,n
3885           sum=sum-a(i,j)*b(j)
3886         enddo
3887         b(i)=sum/a(i,i)
3888       enddo
3889 !      write(std_out,*) 'LEAVING LUBKSB...'
3890       return
3891
3892 END SUBROUTINE LUBKSB
```

## m_hide_lapack/ludcmp [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3734 SUBROUTINE ludcmp(a,n,np,indx,id,info)
3735
3736
3737 !This section has been created automatically by the script Abilint (TD).
3738 !Do not modify the following lines by hand.
3739 #undef ABI_FUNC
3740 #define ABI_FUNC 'ludcmp'
3741 !End of the abilint section
3742
3743       implicit none
3744
3745       INTEGER n,np,indx(n),NMAX,id,info
3746       REAL*8 a(np,np),TINY
3747       PARAMETER (NMAX=500,TINY=1.0e-20)
3748
3749       INTEGER i,imax,j,k
3750       REAL*8 aamax,dum,sum,vv(NMAX)
3751
3752 !      write(std_out,*) 'ENTERING LUDCMP...'
3753 !      write(std_out,*) 'in ludcmp n=',n,' np=',np
3754 !      write(std_out,201) ((a(i,j),j=1,n),i=1,n)
3755 ! 201  FORMAT('A in ludcmp ',/,3F16.8,/,3F16.8,/,3F16.8)
3756       id=1
3757       info=0
3758       do i=1,n
3759         aamax=0.
3760         do j=1,n
3761           if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
3762         enddo
3763         if (aamax.eq.0.) then
3764           write(std_out,*) 'LUDCMP: singular matrix !!!'
3765           do j=1,3
3766             write(std_out,*) (a(j,k),k=1,3)
3767           enddo
3768           info=1
3769           return
3770 !          stop 'singular matrix in ludcmp'
3771         endif
3772         vv(i)=1./aamax
3773       enddo
3774       do j=1,n
3775         do i=1,j-1
3776           sum=a(i,j)
3777           do k=1,i-1
3778             sum=sum-a(i,k)*a(k,j)
3779           enddo
3780           a(i,j)=sum
3781         enddo
3782         aamax=0.
3783         do i=j,n
3784           sum=a(i,j)
3785           do k=1,j-1
3786             sum=sum-a(i,k)*a(k,j)
3787           enddo
3788           a(i,j)=sum
3789           dum=vv(i)*abs(sum)
3790           if (dum.ge.aamax) then
3791             imax=i
3792             aamax=dum
3793           endif
3794         enddo
3795         if (j.ne.imax)then
3796           do  k=1,n
3797             dum=a(imax,k)
3798             a(imax,k)=a(j,k)
3799             a(j,k)=dum
3800           enddo
3801           id=-id
3802           vv(imax)=vv(j)
3803         endif
3804         indx(j)=imax
3805         if(a(j,j).eq.0.)a(j,j)=TINY
3806         if(j.ne.n)then
3807           dum=1./a(j,j)
3808           do i=j+1,n
3809             a(i,j)=a(i,j)*dum
3810           enddo
3811         endif
3812       enddo
3813 !      write(std_out,*) 'LEAVING LUDCMP...'
3814       return
3815 END SUBROUTINE ludcmp
```

## m_hide_lapack/matr3eigval [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3535 subroutine matr3eigval(eigval,matr)
3536
3537
3538 !This section has been created automatically by the script Abilint (TD).
3539 !Do not modify the following lines by hand.
3540 #undef ABI_FUNC
3541 #define ABI_FUNC 'matr3eigval'
3542 !End of the abilint section
3543
3544  implicit none
3545
3546 !Arguments ------------------------------------
3547 !arrays
3548  real(dp),intent(in) :: matr(3,3)
3549  real(dp),intent(out) :: eigval(3)
3550
3551 !Local variables-------------------------------
3552 !scalars
3553  integer :: ier
3554 !arrays
3555  real(dp) :: eigvec(2,3,3),matrx(2,6),zhpev1(2,2*3-1),zhpev2(3*3-2)
3556
3557 ! *************************************************************************
3558
3559  matrx(1,1)=matr(1,1)
3560  matrx(1,2)=matr(1,2)
3561  matrx(1,3)=matr(2,2)
3562  matrx(1,4)=matr(1,3)
3563  matrx(1,5)=matr(2,3)
3564  matrx(1,6)=matr(3,3)
3565  matrx(2,:)=zero
3566
3567  call ZHPEV ('V','U',3,matrx,eigval,eigvec,3,zhpev1,zhpev2,ier)
3568 !write(std_out,*)' eigval=',eigval
3569
3570 end subroutine matr3eigval
```

## m_hide_lapack/matrginv [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3402 subroutine matrginv(a,lda,n)
3403
3404
3405 !This section has been created automatically by the script Abilint (TD).
3406 !Do not modify the following lines by hand.
3407 #undef ABI_FUNC
3408 #define ABI_FUNC 'matrginv'
3409 !End of the abilint section
3410
3411  implicit none
3412
3413 !Arguments ------------------------------------
3414 !scalars
3415  integer,intent(in) :: lda,n
3416 !arrays
3417  real(dp),intent(inout) :: a(lda,n)
3418
3419 !Local variables-------------------------------
3420 !scalars
3421  integer :: ierr,nwork
3422 #if defined HAVE_LINALG_ESSL
3423  real(dp) :: rcond
3424 #endif
3425  character(len=500) :: message
3426 !arrays
3427  integer,allocatable :: ipvt(:)
3428 #if defined HAVE_LINALG_ESSL
3429  real(dp) :: det(2)
3430 #elif defined HAVE_LINALG_ASL
3431  real(dp) :: det(2)
3432 #endif
3433  real(dp),allocatable :: work(:)
3434
3435 ! *************************************************************************
3436
3437 #if defined HAVE_LINALG_ESSL
3438  nwork=200*n
3439 #else
3440  nwork=n
3441 #endif
3442
3443  ABI_ALLOCATE(work,(nwork))
3444  ABI_ALLOCATE(ipvt,(n))
3445
3446
3447 #if defined HAVE_LINALG_ESSL
3448
3449  call dgeicd(a,lda,n,0,rcond,det,work,nwork)
3450  if(abs(rcond)==zero) then
3451    write(message, '(10a)' ) ch10,&
3452 &   ' matrginv : BUG -',ch10,&
3453 &   '  The matrix that has been passed in argument of this subroutine',ch10,&
3454 &   '  is probably either singular or nearly singular.',ch10,&
3455 &   '  The ESSL routine dgeicd failed.',ch10,&
3456 &   '  Action: Contact ABINIT group '
3457    MSG_ERROR(message)
3458  end if
3459
3460 #elif defined HAVE_LINALG_ASL
3461
3462  call dbgmlu(a,lda,n,ipvt,ierr)
3463  if(ierr /= 0) then
3464    write(message, '(10a)' ) ch10,&
3465 &   ' matrginv : BUG -',ch10,&
3466 &   '  The matrix that has been passed in argument of this subroutine',ch10,&
3467 &   '  is probably either singular or nearly singular.',ch10,&
3468 &   '  The ASL routine dbgmlu failed.',ch10,&
3469 &   '  Action: Contact ABINIT group '
3470    MSG_ERROR(message)
3471  end if
3472  call dbgmdi(a,lda,n,ipvt,det,-1,work,ierr)
3473  if(ierr /= 0) then
3474    write(message, '(10a)' ) ch10,&
3475 &   ' matrginv : BUG -',ch10,&
3476 &   '  The matrix that has been passed in argument of this subroutine',ch10,&
3477 &   '  is probably either singular or nearly singular.',ch10,&
3478 &   '  The ASL routine dbgmdi failed.',ch10,&
3479 &   '  Action: Contact ABINIT group '
3480    MSG_ERROR(message)
3481  end if
3482
3483 #else
3484
3485  call dgetrf(n,n,a,lda,ipvt,ierr)
3486  if(ierr /= 0) then
3487    write(message, '(10a)' ) ch10,&
3488 &   ' matrginv : BUG -',ch10,&
3489 &   '  The matrix that has been passed in argument of this subroutine',ch10,&
3490 &   '  is probably either singular or nearly singular.',ch10,&
3491 &   '  The LAPACK routine dgetrf failed.',ch10,&
3492 &   '  Action: Contact ABINIT group '
3493    MSG_ERROR(message)
3494  end if
3495  call dgetri(n,a,lda,ipvt,work,n,ierr)
3496  if(ierr /= 0) then
3497    write(message, '(10a)' ) ch10,&
3498 &   ' matrginv : BUG -',ch10,&
3499 &   '  The matrix that has been passed in argument of this subroutine',ch10,&
3500 &   '  is probably either singular or nearly singular.',ch10,&
3501 &   '  The LAPACK routine dgetri failed.',ch10,&
3502 &   '  Action: Contact ABINIT group '
3503    MSG_ERROR(message)
3504  end if
3505
3506 #endif
3507
3508  ABI_DEALLOCATE(work)
3509  ABI_DEALLOCATE(ipvt)
3510
3511 end subroutine matrginv
```

## m_hide_lapack/test_xginv [ Functions ]

[ Top ] [ m_hide_lapack ] [ Functions ]

NAME

```  test_xginv
```

FUNCTION

## m_hide_lapack/wrap_CGEEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```2593 subroutine wrap_CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr)
2594
2595
2596 !This section has been created automatically by the script Abilint (TD).
2597 !Do not modify the following lines by hand.
2598 #undef ABI_FUNC
2599 #define ABI_FUNC 'wrap_CGEEV'
2600 !End of the abilint section
2601
2602  implicit none
2603
2604 !Arguments ------------------------------------
2605 !scalars
2606  integer,intent(in) :: n,lda,ldvl,ldvr
2607  character(len=*),intent(in) ::  jobvl,jobvr
2608 !arrays
2609  complex(spc),intent(inout) :: a(lda,n)
2610  complex(spc),intent(out) :: w(n)
2611  complex(spc),intent(out) :: vl(ldvl,n)
2612  complex(spc),intent(out) :: vr(ldvr,n)
2613
2614 !Local variables ------------------------------
2615 !scalars
2616  integer :: info,lwork
2617  character(len=500) :: msg
2618 !arrays
2619  real(sp),allocatable :: rwork(:)
2620  complex(spc),allocatable :: work(:)
2621
2622 !************************************************************************
2623
2624  lwork = MAX(1,2*n)
2625
2626  ABI_MALLOC(work,(lwork))
2627  ABI_MALLOC(rwork,(2*n))
2628
2629  call CGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info)
2630
2631  if (info < 0) then
2632   write(msg,'(a,i0,a)')" The ",-info,"-th argument of CGEEV had an illegal value."
2633   MSG_ERROR(msg)
2634  end if
2635
2636  if (info > 0) then
2637   write(msg,'(3a,i0,a,i0,a)')&
2638 &  "CGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,&
2639 &  "Elements ",info+1,":",n," of W contain eigenvalues which have converged. "
2640   MSG_ERROR(msg)
2641  end if
2642
2643  ABI_FREE(work)
2644  ABI_FREE(rwork)
2645
2646 end subroutine wrap_CGEEV
```

## m_hide_lapack/wrap_CHEEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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_hide_lapack/wrap_CHPEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```709 subroutine wrap_CHPEV(jobz,uplo,n,ap,w,z,ldz)
710
711
712 !This section has been created automatically by the script Abilint (TD).
713 !Do not modify the following lines by hand.
714 #undef ABI_FUNC
715 #define ABI_FUNC 'wrap_CHPEV'
716 !End of the abilint section
717
718  implicit none
719
720 !Arguments ------------------------------------
721 !scalars
722  integer,intent(in) :: n,ldz
723  character(len=*),intent(in) :: jobz,uplo
724 !arrays
725  real(sp),intent(out) :: w(n)
726  complex(spc),intent(inout) :: ap(n*(n+1)/2)
727  complex(spc),intent(out) :: z(ldz,n)
728
729 !Local variables ------------------------------
730 !scalars
731  integer :: info
732  character(len=500) :: msg
733 !arrays
734  real(sp),allocatable :: rwork(:)
735  complex(spc),allocatable :: work(:)
736
737 !************************************************************************
738
739  ABI_MALLOC(work, (MAX(1,2*n-1)))
740  ABI_MALLOC(rwork, (MAX(1,3*n-2)))
741
742  call CHPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, RWORK, INFO )
743
744  if (info < 0) then
745   write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value."
746   MSG_ERROR(msg)
747  end if
748
749  if (info > 0) then
750   write(msg,'(2a,i0,a)')&
751 &  "ZHPEV: the algorithm failed to converge; ",ch10,&
752 &  info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
753   MSG_ERROR(msg)
754  end if
755
756  ABI_FREE(rwork)
757  ABI_FREE(work)
758
759 end subroutine wrap_CHPEV
```

## m_hide_lapack/wrap_DSYEV_ZHEEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```496 subroutine wrap_DSYEV_ZHEEV(jobz,uplo,cplex,n,a,w,comm)
497
498
499 !This section has been created automatically by the script Abilint (TD).
500 !Do not modify the following lines by hand.
501 #undef ABI_FUNC
502 #define ABI_FUNC 'wrap_DSYEV_ZHEEV'
503 !End of the abilint section
504
505  implicit none
506
507 !Arguments ------------------------------------
508 !scalars
509  integer,intent(in) :: n,cplex
510  integer,optional,intent(in) :: comm
511  character(len=*),intent(in) :: jobz,uplo
512 !arrays
513  real(dp),intent(inout) :: a(cplex,n,n)
514  real(dp),intent(out) :: w(n)
515
516 !Local variables ------------------------------
517 !scalars
518  integer :: lwork,info,nprocs
519  logical :: use_scalapack
520  character(len=500) :: msg
521 !arrays
522  real(dp),allocatable :: rwork(:)
523  real(dp),allocatable :: work_real(:)
524  complex(dpc),allocatable :: work_cplx(:)
525 #ifdef HAVE_LINALG_SCALAPACK
526  integer :: ierr,istwf_k,tbloc
527  logical :: want_eigenvectors
528  type(matrix_scalapack)    :: Slk_mat,Slk_vec
529  type(processor_scalapack) :: Slk_processor
530 #endif
531 !************************************************************************
532
533  use_scalapack=.FALSE.
534  if (PRESENT(comm)) then
535   nprocs = xmpi_comm_size(comm)
536 #ifdef HAVE_LINALG_SCALAPACK
537   use_scalapack = (nprocs>1)
538 #endif
539  end if
540
541  if (ALL(cplex/=(/1,2/))) then
542   write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
543   MSG_BUG(msg)
544  end if
545
546  SELECT CASE(use_scalapack)
547
548  CASE (.FALSE.)
549
550   if (cplex==1) then  ! Real symmetric case.
551
552    lwork = MAX(1,3*n-1)
553
554    ABI_MALLOC(work_real,(lwork))
555
556    call DSYEV(jobz,uplo,n,a,n,w,work_real,lwork,info)
557
558    if (info < 0) then
559     write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYEV had an illegal value."
560     MSG_ERROR(msg)
561    end if
562
563    if (info > 0) then
564     write(msg,'(2a,i0,a)')&
565 &    "DSYEV: the algorithm failed to converge; ",ch10,&
566 &    info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
567     MSG_ERROR(msg)
568    end if
569
570    ABI_FREE(work_real)
571
572    RETURN
573
574   else                ! Hermitian case.
575
576    lwork = MAX(1,2*n-1)
577
578    ABI_MALLOC(work_cplx, (lwork))
579    ABI_MALLOC(rwork, (MAX(1,3*n-2)))
580
581    call ZHEEV(jobz,uplo,n,a,n,w,work_cplx,lwork,rwork,info)
582
583    if (info < 0) then
584     write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value."
585     MSG_ERROR(msg)
586    end if
587
588    if (info > 0) then
589     write(msg,'(2a,i0,a)')&
590 &    "ZHEEV: the algorithm failed to converge; ",ch10,&
591 &    info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
592     MSG_ERROR(msg)
593    end if
594
595    ABI_FREE(rwork)
596    ABI_FREE(work_cplx)
597
598    RETURN
599   end if ! cplex
600
601  CASE (.TRUE.)
602
603 #ifdef HAVE_LINALG_SCALAPACK
604
605   MSG_ERROR("Not coded yet")
606
607 !  call init_scalapack(Slk_processor,comm)
608 !  istwf_k=1
609 !
610 !  ! Initialize and fill Scalapack matrix from the global one.
611 !  tbloc=SLK_BLOCK_SIZE
612 !  call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
613 !
614 !  write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
615 !  call wrtout(std_out,msg,"PERS")
616 !
617 !  call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
618 !
619 !  want_eigenvectors = firstchar(jobz,(/"V","v"/))
620 !  if (want_eigenvectors) then ! Initialize the distributed vectors.
621 !   call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
622 !  end if
623 !
624 !  ! Solve the problem with scaLAPACK.
625 !  call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
626 !
627 !  call destruction_matrix_scalapack(Slk_mat)
628 !
629 !  if (want_eigenvectors) then ! A is overwritten with the eigenvectors
630 !   a = czero
631 !   call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
632 !   call destruction_matrix_scalapack(Slk_vec)
633 !   call xmpi_sum(a,comm,ierr)                        ! Fill the remaing entries of the global matrix
634 !  end if
635 !
636 !  call end_scalapack(Slk_processor)
637
638   RETURN
639 #endif
640
641   MSG_BUG("You should not be here!")
642
643  END SELECT
644
645 end subroutine wrap_DSYEV_ZHEEV
```

## m_hide_lapack/wrap_DSYEVX_ZHEEVX [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```1759 subroutine wrap_DSYEVX_ZHEEVX(jobz,range,uplo,cplex,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
1760
1761
1762 !This section has been created automatically by the script Abilint (TD).
1763 !Do not modify the following lines by hand.
1764 #undef ABI_FUNC
1765 #define ABI_FUNC 'wrap_DSYEVX_ZHEEVX'
1766 !End of the abilint section
1767
1768  implicit none
1769
1770 !Arguments ------------------------------------
1771 !scalars
1772  integer,intent(in) :: il,iu,ldz,n,cplex
1773  integer,optional,intent(in) :: comm
1774  integer,intent(inout) :: m
1775  real(dp),intent(in) :: abstol,vl,vu
1776  character(len=*),intent(in) :: jobz,range,uplo
1777 !arrays
1778  real(dp),intent(out) :: w(n)
1779  !real(dp),intent(out) :: z(cplex,ldz,n)
1780  real(dp),intent(out) :: z(cplex,ldz,m)
1781  real(dp),intent(inout) :: a(cplex,n,n)
1782
1783 !Local variables ------------------------------
1784 !scalars
1785  integer :: lwork,info,nprocs
1786  logical :: use_scalapack
1787  character(len=500) :: msg
1788 !arrays
1789  integer,allocatable :: ifail(:),iwork(:)
1790  real(dp),allocatable :: rwork(:)
1791  real(dp),allocatable :: work_real(:)
1792  complex(dpc),allocatable :: work_cplx(:)
1793 #ifdef HAVE_LINALG_SCALAPACK
1794  integer :: ierr,istwf_k,tbloc
1795  logical :: want_eigenvectors
1796  type(matrix_scalapack)    :: Slk_mat,Slk_vec
1797  type(processor_scalapack) :: Slk_processor
1798 #endif
1799
1800 !************************************************************************
1801
1802  use_scalapack=.FALSE.
1803  if (PRESENT(comm)) then
1804   nprocs = xmpi_comm_size(comm)
1805 #ifdef HAVE_LINALG_SCALAPACK
1806   use_scalapack = (nprocs>1)
1807 #endif
1808  end if
1809
1810  if (ALL(cplex/=(/1,2/))) then
1811   write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
1812   MSG_ERROR(msg)
1813  end if
1814
1815  SELECT CASE(use_scalapack)
1816
1817  CASE (.FALSE.) ! Standard LAPACK call.
1818
1819   if (cplex==1) then      ! Real symmetric case
1820
1821    lwork = MAX(1,8*n)
1822
1823    ABI_MALLOC(work_real,(lwork))
1824    ABI_MALLOC(iwork,(5*n))
1825    ABI_MALLOC(ifail,(n))
1826
1827    call DSYEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_real,lwork,iwork,ifail,info)
1828
1829    if (info < 0) then
1830     write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYEVX had an illegal value."
1831     MSG_ERROR(msg)
1832    end if
1833
1834    if (info > 0) then
1835     write(msg,'(2a,i0,a)')&
1836 &    "DSYEVX: the algorithm failed to converge; ",ch10,&
1837 &    info,"eigenvectors failed to converge. "
1838     MSG_ERROR(msg)
1839    end if
1840
1841    ABI_FREE(work_real)
1842    ABI_FREE(iwork)
1843    ABI_FREE(ifail)
1844
1845    RETURN
1846
1847   else                    ! Complex Hermitian case.
1848
1849    lwork = MAX(1,2*n)
1850
1851    ABI_MALLOC(work_cplx,(lwork))
1852    ABI_MALLOC(rwork,(7*n))
1853    ABI_MALLOC(iwork,(5*n))
1854    ABI_MALLOC(ifail,(n))
1855
1856    call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work_cplx,lwork,rwork,iwork,ifail,info)
1857
1858    if (info < 0) then
1859     write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value."
1860     MSG_ERROR(msg)
1861    end if
1862
1863    if (info > 0) then
1864     write(msg,'(2a,i0,a)')&
1865 &    "ZHEEVX: the algorithm failed to converge; ",ch10,&
1866 &    info,"eigenvectors failed to converge. "
1867     MSG_ERROR(msg)
1868    end if
1869
1870    ABI_FREE(iwork)
1871    ABI_FREE(ifail)
1872    ABI_FREE(rwork)
1873    ABI_FREE(work_cplx)
1874
1875    RETURN
1876
1877   end if
1878
1879  CASE (.TRUE.)
1880
1881 #ifdef HAVE_LINALG_SCALAPACK
1882
1883   MSG_ERROR("Not coded yet")
1884   ! call init_scalapack(Slk_processor,comm)
1885   ! istwf_k=1
1886
1887   ! ! Initialize and fill Scalapack matrix from the global one.
1888   ! tbloc=SLK_BLOCK_SIZE
1889   ! call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1890
1891   ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
1892   ! call wrtout(std_out,msg,"COLL")
1893
1894   ! call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
1895
1896   ! want_eigenvectors = firstchar(jobz,(/"V","v"/))
1897   ! if (want_eigenvectors) then ! Initialize the distributed vectors.
1898   !  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1899   ! end if
1900
1901   ! ! Solve the problem.
1902   ! call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w)
1903
1904   ! call destruction_matrix_scalapack(Slk_mat)
1905   !
1906   ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors
1907   !  z = czero
1908   !  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
1909   !  call destruction_matrix_scalapack(Slk_vec)
1910   !  call xmpi_sum(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
1911   ! end if
1912
1913   ! call end_scalapack(Slk_processor)
1914
1915   RETURN
1916 #endif
1917
1918   MSG_BUG("You should not be here!")
1919
1920  END SELECT
1921
1922 end subroutine wrap_DSYEVX_ZHEEVX
```

## m_hide_lapack/wrap_DSYGV_ZHEGV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```1222 subroutine wrap_DSYGV_ZHEGV(itype,jobz,uplo,cplex,n,a,b,w,comm)
1223
1224
1225 !This section has been created automatically by the script Abilint (TD).
1226 !Do not modify the following lines by hand.
1227 #undef ABI_FUNC
1228 #define ABI_FUNC 'wrap_DSYGV_ZHEGV'
1229 !End of the abilint section
1230
1231  implicit none
1232
1233 !Arguments ------------------------------------
1234 !scalars
1235  integer,intent(in) :: n,itype,cplex
1236  integer,optional,intent(in) :: comm
1237  character(len=*),intent(in) :: jobz,uplo
1238 !arrays
1239  real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n)
1240  real(dp),intent(out) :: w(n)
1241
1242 !Local variables ------------------------------
1243 !scalars
1244  integer :: lwork,info,nprocs,ii
1245  logical :: use_scalapack
1246  character(len=500) :: msg
1247 !arrays
1248  real(dp),allocatable :: rwork(:)
1249  real(dp),allocatable :: work_real(:)
1250  complex(dpc),allocatable :: work_cplx(:)
1251 #ifdef HAVE_LINALG_SCALAPACK
1252  integer :: ierr,istwf_k,tbloc
1253  type(matrix_scalapack)    :: Slk_matA,Slk_matB
1254  type(processor_scalapack) :: Slk_processor
1255 #endif
1256 !************************************************************************
1257
1258  use_scalapack=.FALSE.
1259  if (PRESENT(comm)) then
1260   nprocs = xmpi_comm_size(comm)
1261 #ifdef HAVE_LINALG_SCALAPACK
1262   use_scalapack = (nprocs>1)
1263 #endif
1264  end if
1265
1266  if (ALL(cplex/=(/1,2/))) then
1267   write(msg,'(a,i0)')"Wrong value for cplex: ",cplex
1268   MSG_ERROR(msg)
1269  end if
1270
1271  SELECT CASE(use_scalapack)
1272
1273  CASE (.FALSE.)
1274
1275   if (cplex==1) then ! Real symmetric case.
1276
1277    lwork = MAX(1,3*n-1)
1278
1279    ABI_MALLOC(work_real,(lwork))
1280
1281    call DSYGV(itype,jobz,uplo,n,a,n,b,n,w,work_real,lwork,info)
1282
1283    if (info < 0) then
1284     write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYGV had an illegal value."
1285     MSG_ERROR(msg)
1286    end if
1287
1288    if (info > 0) then
1289     if (info<= n) then
1290      write(msg,'(2a,i0,a)')&
1291 &     " DSYGV failed to converge: ",ch10,&
1292 &     info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
1293     else
1294      ii = info -n
1295      write(msg,'(3a,i0,3a)')&
1296 &    "DSYGV failed to converge: ",ch10,&
1297 &    "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1298 &    "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1299     end if
1300     MSG_ERROR(msg)
1301    end if
1302
1303    ABI_FREE(work_real)
1304
1305    RETURN
1306
1307   else               ! complex Hermitian case
1308
1309    lwork = MAX(1,2*n-1)
1310
1311    ABI_MALLOC(work_cplx,(lwork))
1312    ABI_MALLOC(rwork,(MAX(1,3*n-2)))
1313
1314    call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work_cplx,lwork,rwork,info)
1315
1316    if (info < 0) then
1317     write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value."
1318     MSG_ERROR(msg)
1319    end if
1320
1321    if (info > 0) then
1322     if (info<= n) then
1323      write(msg,'(2a,i0,a)')&
1324 &     "ZHEEV failed to converge: ",ch10,&
1325 &     info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
1326     else
1327      ii = info -n
1328      write(msg,'(3a,i0,3a)')&
1329 &    "ZHEEV failed to converge: ",ch10,&
1330 &    "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1331 &    "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1332     end if
1333     MSG_ERROR(msg)
1334    end if
1335
1336    ABI_FREE(rwork)
1337    ABI_FREE(work_cplx)
1338
1339    RETURN
1340
1341   end if ! cplex
1342
1343  CASE (.TRUE.)
1344
1345 #ifdef HAVE_LINALG_SCALAPACK
1346
1347   MSG_ERROR("Not coded yet")
1348   ! call init_scalapack(Slk_processor,comm)
1349   ! istwf_k=1
1350
1351   ! ! Initialize and fill Scalapack matrix from the global one.
1352   ! tbloc=SLK_BLOCK_SIZE
1353
1354   ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
1355   ! call wrtout(std_out,msg,"COLL")
1356
1357   ! call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1358   ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
1359
1360   ! call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1361   ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
1362
1363   ! ! Solve the problem with scaLAPACK.
1364   ! MSG_ERROR("slk_pZHEGV not yet coded")
1365   ! ! TODO
1366   ! call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
1367
1368   ! call destruction_matrix_scalapack(Slk_matB)
1369   !
1370   ! if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors
1371   !  a = czero
1372   !  call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
1373   !  call xmpi_sum(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
1374   ! end if
1375
1376   ! call destruction_matrix_scalapack(Slk_matA)
1377
1378   ! call end_scalapack(Slk_processor)
1379
1380   RETURN
1381 #endif
1382
1383   MSG_BUG("You should not be here!")
1384
1385  END SELECT
1386
1387 end subroutine wrap_DSYGV_ZHEGV
```

## m_hide_lapack/wrap_DSYGVX_ZHEGVX [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```2334 subroutine wrap_DSYGVX_ZHEGVX(itype,jobz,range,uplo,cplex,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
2335
2336
2337 !This section has been created automatically by the script Abilint (TD).
2338 !Do not modify the following lines by hand.
2339 #undef ABI_FUNC
2340 #define ABI_FUNC 'wrap_DSYGVX_ZHEGVX'
2341 !End of the abilint section
2342
2343  implicit none
2344
2345 !Arguments ------------------------------------
2346 !scalars
2347  integer,intent(in) :: il,iu,ldz,n,itype,cplex
2348  integer,optional,intent(in) :: comm
2349  integer,intent(inout) :: m
2350  real(dp),intent(in) :: abstol,vl,vu
2351  character(len=*),intent(in) :: jobz,range,uplo
2352 !arrays
2353  real(dp),intent(out) :: w(n)
2354  !real(dp),intent(out) :: z(cplex,ldz,n)
2355  real(dp),intent(out) :: z(cplex,ldz,m)
2356  real(dp),intent(inout) :: a(cplex,n,n),b(cplex,n,n)
2357
2358 !Local variables ------------------------------
2359 !scalars
2360  integer :: lwork,info,nprocs,ii
2361  logical :: use_scalapack
2362  character(len=500) :: msg
2363 !arrays
2364  integer,allocatable :: ifail(:),iwork(:)
2365  real(dp),allocatable :: rwork(:)
2366  real(dp),allocatable :: work_real(:)
2367  complex(dpc),allocatable :: work_cplx(:)
2368 #ifdef HAVE_LINALG_SCALAPACK
2369  integer :: ierr,istwf_k,tbloc
2370  logical :: want_eigenvectors
2371  type(matrix_scalapack)    :: Slk_matA,Slk_matB,Slk_vec
2372  type(processor_scalapack) :: Slk_processor
2373 #endif
2374
2375 !************************************************************************
2376
2377  use_scalapack=.FALSE.
2378  if (PRESENT(comm)) then
2379   nprocs = xmpi_comm_size(comm)
2380 #ifdef HAVE_LINALG_SCALAPACK
2381   use_scalapack = (nprocs>1)
2382 #endif
2383  end if
2384
2385  if (ALL(cplex/=(/1,2/))) then
2386   write(msg,'(a,i0)')" Wrong value for cplex: ",cplex
2387   MSG_ERROR(msg)
2388  end if
2389
2390  SELECT CASE(use_scalapack)
2391
2392  CASE (.FALSE.) ! Standard LAPACK call.
2393
2394   if (cplex==1) then  ! Real symmetric case
2395
2396    lwork = MAX(1,8*n)
2397
2398    ABI_MALLOC(work_real,(lwork))
2399    ABI_MALLOC(iwork,(5*n))
2400    ABI_MALLOC(ifail,(n))
2401
2402    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)
2403
2404    if (info < 0) then
2405     write(msg,'(a,i0,a)')" The ",-info,"-th argument of DSYGVX had an illegal value."
2406     MSG_ERROR(msg)
2407    end if
2408
2409    if (info > 0) then
2410     if (info<= n) then
2411      write(msg,'(a,i0,a)')&
2412 &     " DSYGVX failed to converge: ",info," eigenvectors failed to converge. "
2413     else
2414      ii = info -n
2415      write(msg,'(3a,i0,3a)')&
2416 &    " DSYGVX failed to converge: ",ch10,&
2417 &    " The leading minor of order ",ii," of B is not positive definite. ",ch10,&
2418 &    " The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
2419     end if
2420     MSG_ERROR(msg)
2421    end if
2422
2423    ABI_FREE(iwork)
2424    ABI_FREE(ifail)
2425    ABI_FREE(work_real)
2426
2427    RETURN
2428
2429   else                ! Complex Hermitian case.
2430
2431    lwork = MAX(1,2*n)
2432
2433    ABI_MALLOC(work_cplx,(lwork))
2434    ABI_MALLOC(rwork,(7*n))
2435    ABI_MALLOC(iwork,(5*n))
2436    ABI_MALLOC(ifail,(n))
2437
2438    !write(std_out,*)"Calling ZHEGVX"
2439
2440    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)
2441
2442    if (info < 0) then
2443     write(msg,'(a,i0,a)')"The ",-info,"-th argument of ZHEGVX had an illegal value."
2444     MSG_ERROR(msg)
2445    end if
2446
2447    if (info > 0) then
2448     if (info<= n) then
2449      write(msg,'(a,i0,a)')&
2450 &     "ZHEGVX failed to converge: ",info," eigenvectors failed to converge. "
2451     else
2452      ii = info -n
2453      write(msg,'(3a,i0,3a)')&
2454 &    "ZHEEVX failed to converge: ",ch10,&
2455 &    "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
2456 &    "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
2457     end if
2458     MSG_ERROR(msg)
2459    end if
2460
2461    ABI_FREE(iwork)
2462    ABI_FREE(ifail)
2463    ABI_FREE(rwork)
2464    ABI_FREE(work_cplx)
2465
2466    RETURN
2467
2468   end if ! cplex
2469
2470  CASE (.TRUE.)
2471
2472 #ifdef HAVE_LINALG_SCALAPACK
2473
2474   MSG_ERROR("not coded yet")
2475   ! call init_scalapack(Slk_processor,comm)
2476   ! istwf_k=1
2477
2478   ! ! Initialize and fill Scalapack matrix from the global one.
2479   ! tbloc=SLK_BLOCK_SIZE
2480
2481   ! write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
2482   ! call wrtout(std_out,msg,"COLL")
2483
2484   ! call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2485   ! call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
2486
2487   ! call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2488   ! call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
2489
2490   ! want_eigenvectors = firstchar(jobz,(/"V","v"/))
2491   ! if (want_eigenvectors) then ! Initialize the distributed vectors.
2492   !  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2493   ! end if
2494
2495   ! ! Solve the problem.
2496   ! MSG_ERROR("slk_pZHEGVX not coded yet")
2497   ! ! TODO write the scaLAPACK wrapper.
2498   ! call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
2499
2500   ! call destruction_matrix_scalapack(Slk_matA)
2501   ! call destruction_matrix_scalapack(Slk_matB)
2502   !
2503   ! if (want_eigenvectors) then ! A is overwritten with the eigenvectors
2504   !  z = czero
2505   !  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
2506   !  call destruction_matrix_scalapack(Slk_vec)
2507   !  call xmpi_sum(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
2508   ! end if
2509
2510   ! call end_scalapack(Slk_processor)
2511
2512   ! RETURN
2513 #endif
2514
2515   MSG_BUG("You should not be here!")
2516
2517  END SELECT
2518
2519 end subroutine wrap_DSYGVX_ZHEGVX
```

## m_hide_lapack/wrap_ZGEEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```2721 subroutine wrap_ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr)
2722
2723
2724 !This section has been created automatically by the script Abilint (TD).
2725 !Do not modify the following lines by hand.
2726 #undef ABI_FUNC
2727 #define ABI_FUNC 'wrap_ZGEEV'
2728 !End of the abilint section
2729
2730  implicit none
2731
2732 !Arguments ------------------------------------
2733 !scalars
2734  integer,intent(in) :: n,lda,ldvl,ldvr
2735  character(len=*),intent(in) ::  jobvl,jobvr
2736 !arrays
2737  complex(dpc),intent(inout) :: a(lda,n)
2738  complex(dpc),intent(out) :: w(n)
2739  complex(dpc),intent(out) :: vl(ldvl,n)
2740  complex(dpc),intent(out) :: vr(ldvr,n)
2741
2742 !Local variables ------------------------------
2743 !scalars
2744  integer :: info,lwork
2745  logical :: use_scalapack
2746  character(len=500) :: msg
2747 !arrays
2748  real(dp),allocatable :: rwork(:)
2749  complex(dpc),allocatable :: work(:)
2750
2751 !************************************************************************
2752
2753  use_scalapack=.FALSE.
2754
2755  SELECT CASE(use_scalapack)
2756
2757  CASE (.FALSE.)
2758
2759   lwork = MAX(1,2*n)
2760
2761   ABI_MALLOC(work,(lwork))
2762   ABI_MALLOC(rwork,(2*n))
2763
2764   call ZGEEV(jobvl,jobvr,n,a,lda,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info)
2765
2766   if (info < 0) then
2767    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGEEV had an illegal value."
2768    MSG_ERROR(msg)
2769   end if
2770
2771   if (info > 0) then
2772    write(msg,'(3a,i0,a,i0,a)')&
2773 &   "ZGEEV: The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed;",ch10,&
2774 &   "Elements ",info+1,":",n," of W contain eigenvalues which have converged. "
2775    MSG_ERROR(msg)
2776   end if
2777
2778   ABI_FREE(work)
2779   ABI_FREE(rwork)
2780
2781   RETURN
2782
2783  CASE (.TRUE.)
2784
2785   MSG_BUG("You should not be here!")
2786
2787  END SELECT
2788
2789 end subroutine wrap_ZGEEV
```

## m_hide_lapack/wrap_ZHEEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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 !End of the abilint section
332
333  implicit none
334
335 !Arguments ------------------------------------
336 !scalars
337  integer,intent(in) :: n
338  integer,optional,intent(in) :: comm
339  character(len=*),intent(in) :: jobz,uplo
340 !arrays
341  complex(dpc),intent(inout) :: a(n,n)
342  real(dp),intent(out) :: w(n)
343
344 !Local variables ------------------------------
345 !scalars
346  integer :: lwork,info,nprocs
347  logical :: use_scalapack
348  character(len=500) :: msg
349 !arrays
350  real(dp),allocatable :: rwork(:)
351  complex(dpc),allocatable :: work(:)
352 #ifdef HAVE_LINALG_SCALAPACK
353  integer :: ierr,istwf_k,tbloc
354  logical :: want_eigenvectors
355  type(matrix_scalapack)    :: Slk_mat,Slk_vec
356  type(processor_scalapack) :: Slk_processor
357 #endif
358 !************************************************************************
359
360  use_scalapack=.FALSE.
361  if (PRESENT(comm)) then
362   nprocs = xmpi_comm_size(comm)
363 #ifdef HAVE_LINALG_SCALAPACK
364   use_scalapack = (nprocs>1)
365 #endif
366  end if
367
368  SELECT CASE(use_scalapack)
369
370  CASE (.FALSE.)
371
372   lwork = MAX(1,2*n-1)
373
374   ABI_MALLOC(work, (lwork))
375   ABI_MALLOC(rwork, (MAX(1,3*n-2)))
376
377   call ZHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info)
378
379   if (info < 0) then
380    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEV had an illegal value."
381    MSG_ERROR(msg)
382   end if
383
384   if (info > 0) then
385    write(msg,'(2a,i0,a)')&
386 &   "ZHEEV: the algorithm failed to converge; ",ch10,&
387 &   info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
388    MSG_ERROR(msg)
389   end if
390
391   ABI_FREE(rwork)
392   ABI_FREE(work)
393
394   RETURN
395
396  CASE (.TRUE.)
397
398 #ifdef HAVE_LINALG_SCALAPACK
399   call init_scalapack(Slk_processor,comm)
400   istwf_k=1
401
402   ! Initialize and fill Scalapack matrix from the global one.
403   tbloc=SLK_BLOCK_SIZE
404   call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
405
406   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
407   call wrtout(std_out,msg,"COLL")
408
409   call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
410
411   want_eigenvectors = firstchar(jobz,(/"V","v"/))
412   if (want_eigenvectors) then ! Initialize the distributed vectors.
413    call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
414   end if
415
416   ! Solve the problem with scaLAPACK.
417   call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
418
419   call destruction_matrix_scalapack(Slk_mat)
420
421   if (want_eigenvectors) then ! A is overwritten with the eigenvectors
422    a = czero
423    call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
424    call destruction_matrix_scalapack(Slk_vec)
425    call xmpi_sum(a,comm,ierr)                        ! Fill the remaing entries of the global matrix
426   end if
427
428   call end_scalapack(Slk_processor)
429
430   RETURN
431 #endif
432
433   MSG_BUG("You should not be here!")
434
435  END SELECT
436
437 end subroutine wrap_ZHEEV
```

## m_hide_lapack/wrap_ZHEEVX [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```1510 subroutine wrap_ZHEEVX(jobz,range,uplo,n,a,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
1511
1512
1513 !This section has been created automatically by the script Abilint (TD).
1514 !Do not modify the following lines by hand.
1515 #undef ABI_FUNC
1516 #define ABI_FUNC 'wrap_ZHEEVX'
1517 !End of the abilint section
1518
1519  implicit none
1520
1521 !Arguments ------------------------------------
1522 !scalars
1523  integer,intent(in) :: il,iu,ldz,n
1524  integer,optional,intent(in) :: comm
1525  integer,intent(inout) :: m
1526  real(dp),intent(in) :: abstol,vl,vu
1527  character(len=*),intent(in) :: jobz,range,uplo
1528 !arrays
1529  real(dp),intent(out) :: w(n)
1530  complex(dpc),intent(out) :: z(ldz,m)
1531  complex(dpc),intent(inout) :: a(n,n)
1532
1533 !Local variables ------------------------------
1534 !scalars
1535  integer :: lwork,info,nprocs
1536  logical :: use_scalapack
1537  character(len=500) :: msg
1538 !arrays
1539  integer,allocatable :: ifail(:),iwork(:)
1540  real(dp),allocatable :: rwork(:)
1541  complex(dpc),allocatable :: work(:)
1542 #ifdef HAVE_LINALG_SCALAPACK
1543  integer :: ierr,istwf_k,tbloc
1544  logical :: want_eigenvectors
1545  type(matrix_scalapack)    :: Slk_mat,Slk_vec
1546  type(processor_scalapack) :: Slk_processor
1547 #endif
1548
1549 !************************************************************************
1550
1551  use_scalapack=.FALSE.
1552  if (PRESENT(comm)) then
1553   nprocs = xmpi_comm_size(comm)
1554 #ifdef HAVE_LINALG_SCALAPACK
1555   use_scalapack = (nprocs>1)
1556 #endif
1557  end if
1558
1559  SELECT CASE(use_scalapack)
1560
1561  CASE (.FALSE.) ! Standard LAPACK call.
1562
1563   lwork = MAX(1,2*n)
1564
1565   ABI_MALLOC(work,(lwork))
1566   ABI_MALLOC(rwork,(7*n))
1567   ABI_MALLOC(iwork,(5*n))
1568   ABI_MALLOC(ifail,(n))
1569
1570   call ZHEEVX(jobz,range,uplo,n,a,n,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info)
1571
1572   if (info < 0) then
1573    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEEVX had an illegal value."
1574    MSG_ERROR(msg)
1575   end if
1576
1577   if (info > 0) then
1578    write(msg,'(2a,i0,a)')&
1579 &   "ZHEEVX: the algorithm failed to converge; ",ch10,&
1580 &   info,"eigenvectors failed to converge. "
1581    MSG_ERROR(msg)
1582   end if
1583
1584   ABI_FREE(iwork)
1585   ABI_FREE(ifail)
1586   ABI_FREE(rwork)
1587   ABI_FREE(work)
1588
1589   RETURN
1590
1591  CASE (.TRUE.)
1592
1593 #ifdef HAVE_LINALG_SCALAPACK
1594   call init_scalapack(Slk_processor,comm)
1595   istwf_k=1
1596
1597   ! Initialize and fill Scalapack matrix from the global one.
1598   tbloc=SLK_BLOCK_SIZE
1599   call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1600
1601   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
1602   call wrtout(std_out,msg,"COLL")
1603
1604   call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
1605
1606   want_eigenvectors = firstchar(jobz,(/"V","v"/))
1607   if (want_eigenvectors) then ! Initialize the distributed vectors.
1608    call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1609   end if
1610
1611   ! Solve the problem.
1612   call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w)
1613
1614   call destruction_matrix_scalapack(Slk_mat)
1615
1616   if (want_eigenvectors) then ! A is overwritten with the eigenvectors
1617    z = czero
1618    call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
1619    call destruction_matrix_scalapack(Slk_vec)
1620    call xmpi_sum(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
1621   end if
1622
1623   call end_scalapack(Slk_processor)
1624
1625   RETURN
1626 #endif
1627
1628   MSG_BUG("You should not be here!")
1629
1630  END SELECT
1631
1632 end subroutine wrap_ZHEEVX
```

## m_hide_lapack/wrap_ZHEGV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```1015 subroutine wrap_ZHEGV(itype,jobz,uplo,n,a,b,w,comm)
1016
1017
1018 !This section has been created automatically by the script Abilint (TD).
1019 !Do not modify the following lines by hand.
1020 #undef ABI_FUNC
1021 #define ABI_FUNC 'wrap_ZHEGV'
1022 !End of the abilint section
1023
1024  implicit none
1025
1026 !Arguments ------------------------------------
1027 !scalars
1028  integer,intent(in) :: n,itype
1029  integer,optional,intent(in) :: comm
1030  character(len=*),intent(in) :: jobz,uplo
1031 !arrays
1032  complex(dpc),intent(inout) :: a(n,n),b(n,n)
1033  real(dp),intent(out) :: w(n)
1034
1035 !Local variables ------------------------------
1036 !scalars
1037  integer :: lwork,info,nprocs,ii
1038  logical :: use_scalapack
1039  character(len=500) :: msg
1040 !arrays
1041  real(dp),allocatable :: rwork(:)
1042  complex(dpc),allocatable :: work(:)
1043 #ifdef HAVE_LINALG_SCALAPACK
1044  integer :: ierr,istwf_k,tbloc
1045  type(matrix_scalapack)    :: Slk_matA,Slk_matB
1046  type(processor_scalapack) :: Slk_processor
1047 #endif
1048 !************************************************************************
1049
1050  use_scalapack=.FALSE.
1051  if (PRESENT(comm)) then
1052   nprocs = xmpi_comm_size(comm)
1053 #ifdef HAVE_LINALG_SCALAPACK
1054   use_scalapack = (nprocs>1)
1055 #endif
1056  end if
1057
1058  SELECT CASE(use_scalapack)
1059
1060  CASE (.FALSE.)
1061
1062   lwork = MAX(1,2*n-1)
1063
1064   ABI_MALLOC(work, (lwork))
1065   ABI_MALLOC(rwork, (MAX(1,3*n-2)))
1066
1067   call ZHEGV(itype,jobz,uplo,n,a,n,b,n,w,work,lwork,rwork,info)
1068
1069   if (info < 0) then
1070    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGV had an illegal value."
1071    MSG_ERROR(msg)
1072   end if
1073
1074   if (info > 0) then
1075    if (info<= n) then
1076     write(msg,'(2a,i0,a)')&
1077 &    "ZHEGV failed to converge: ",ch10,&
1078 &    info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
1079    else
1080     ii = info -n
1081     write(msg,'(3a,i0,3a)')&
1082 &   "ZHEGV failed to converge: ",ch10,&
1083 &   "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
1084 &   "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
1085    end if
1086    MSG_ERROR(msg)
1087   end if
1088
1089   ABI_FREE(rwork)
1090   ABI_FREE(work)
1091
1092   RETURN
1093
1094  CASE (.TRUE.)
1095
1096 #ifdef HAVE_LINALG_SCALAPACK
1097   call init_scalapack(Slk_processor,comm)
1098   istwf_k=1
1099
1100   ! Initialize and fill Scalapack matrix from the global one.
1101   tbloc=SLK_BLOCK_SIZE
1102
1103   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
1104   call wrtout(std_out,msg,"COLL")
1105
1106   call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1107   call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
1108
1109   call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
1110   call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
1111
1112   ! Solve the problem with scaLAPACK.
1113   MSG_ERROR("slk_pZHEGV not yet coded")
1114   ! TODO
1115   !% call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
1116
1117   call destruction_matrix_scalapack(Slk_matB)
1118
1119   if (firstchar(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors
1120    a = czero
1121    call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
1122    call xmpi_sum(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
1123   end if
1124
1125   call destruction_matrix_scalapack(Slk_matA)
1126
1127   call end_scalapack(Slk_processor)
1128
1129   RETURN
1130 #endif
1131
1132   MSG_BUG("You should not be here!")
1133
1134  END SELECT
1135
1136 end subroutine wrap_ZHEGV
```

## m_hide_lapack/wrap_ZHEGVX [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```2056 subroutine wrap_ZHEGVX(itype,jobz,range,uplo,n,a,b,vl,vu,il,iu,abstol,m,w,z,ldz,comm)
2057
2058
2059 !This section has been created automatically by the script Abilint (TD).
2060 !Do not modify the following lines by hand.
2061 #undef ABI_FUNC
2062 #define ABI_FUNC 'wrap_ZHEGVX'
2063 !End of the abilint section
2064
2065  implicit none
2066
2067 !Arguments ------------------------------------
2068 !scalars
2069  integer,intent(in) :: il,iu,ldz,n,itype
2070  integer,optional,intent(in) :: comm
2071  integer,intent(inout) :: m
2072  real(dp),intent(in) :: abstol,vl,vu
2073  character(len=*),intent(in) :: jobz,range,uplo
2074 !arrays
2075  real(dp),intent(out) :: w(n)
2076  !complex(dpc),intent(out) :: z(ldz,n)
2077  complex(dpc),intent(out) :: z(ldz,m)
2078  complex(dpc),intent(inout) :: a(n,n),b(n,n)
2079
2080 !Local variables ------------------------------
2081 !scalars
2082  integer :: lwork,info,nprocs,ii
2083  logical :: use_scalapack
2084  character(len=500) :: msg
2085 !arrays
2086  integer,allocatable :: ifail(:),iwork(:)
2087  real(dp),allocatable :: rwork(:)
2088  complex(dpc),allocatable :: work(:)
2089 #ifdef HAVE_LINALG_SCALAPACK
2090  integer :: ierr,istwf_k,tbloc
2091  logical :: want_eigenvectors
2092  type(matrix_scalapack)    :: Slk_matA,Slk_matB,Slk_vec
2093  type(processor_scalapack) :: Slk_processor
2094 #endif
2095
2096 !************************************************************************
2097
2098  use_scalapack=.FALSE.
2099  if (PRESENT(comm)) then
2100   nprocs = xmpi_comm_size(comm)
2101 #ifdef HAVE_LINALG_SCALAPACK
2102   use_scalapack = (nprocs>1)
2103 #endif
2104  end if
2105
2106  SELECT CASE(use_scalapack)
2107
2108  CASE (.FALSE.) ! Standard LAPACK call.
2109
2110   lwork = MAX(1,2*n)
2111
2112   ABI_MALLOC(work,(lwork))
2113   ABI_MALLOC(rwork,(7*n))
2114   ABI_MALLOC(iwork,(5*n))
2115   ABI_MALLOC(ifail,(n))
2116
2117   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)
2118
2119   if (info < 0) then
2120    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHEGVX had an illegal value."
2121    MSG_ERROR(msg)
2122   end if
2123
2124   if (info > 0) then
2125    if (info<= n) then
2126     write(msg,'(a,i0,a)')&
2127 &    "ZHEGVX failed to converge: ",info," eigenvectors failed to converge. "
2128    else
2129     ii = info -n
2130     write(msg,'(3a,i0,3a)')&
2131 &   "ZHEEVX failed to converge: ",ch10,&
2132 &   "The leading minor of order ",ii," of B is not positive definite. ",ch10,&
2133 &   "The factorization of B could not be completed and no eigenvalues or eigenvectors were computed."
2134    end if
2135    MSG_ERROR(msg)
2136   end if
2137
2138   ABI_FREE(iwork)
2139   ABI_FREE(ifail)
2140   ABI_FREE(rwork)
2141   ABI_FREE(work)
2142
2143   RETURN
2144
2145  CASE (.TRUE.)
2146
2147 #ifdef HAVE_LINALG_SCALAPACK
2148   call init_scalapack(Slk_processor,comm)
2149   istwf_k=1
2150
2151   ! Initialize and fill Scalapack matrix from the global one.
2152   tbloc=SLK_BLOCK_SIZE
2153
2154   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
2155   call wrtout(std_out,msg,"COLL")
2156
2157   call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2158   call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
2159
2160   call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2161   call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
2162
2163   want_eigenvectors = firstchar(jobz,(/"V","v"/))
2164   if (want_eigenvectors) then ! Initialize the distributed vectors.
2165    call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
2166   end if
2167
2168   ! Solve the problem.
2169   MSG_ERROR("slk_pZHEGVX not coded yet")
2170   ! TODO write the scaLAPACK wrapper.
2171   !call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
2172
2173   call destruction_matrix_scalapack(Slk_matA)
2174   call destruction_matrix_scalapack(Slk_matB)
2175
2176   if (want_eigenvectors) then ! A is overwritten with the eigenvectors
2177    z = czero
2178    call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
2179    call destruction_matrix_scalapack(Slk_vec)
2180    call xmpi_sum(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
2181   end if
2182
2183   call end_scalapack(Slk_processor)
2184
2185   RETURN
2186 #endif
2187
2188   MSG_BUG("You should not be here!")
2189
2190  END SELECT
2191
2192 end subroutine wrap_ZHEGVX
```

## m_hide_lapack/wrap_ZHPEV [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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.

```

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

```829 subroutine wrap_ZHPEV(jobz,uplo,n,ap,w,z,ldz,comm)
830
831
832 !This section has been created automatically by the script Abilint (TD).
833 !Do not modify the following lines by hand.
834 #undef ABI_FUNC
835 #define ABI_FUNC 'wrap_ZHPEV'
836 !End of the abilint section
837
838  implicit none
839
840 !Arguments ------------------------------------
841 !scalars
842  integer,intent(in) :: n,ldz
843  integer,optional,intent(in) :: comm
844  character(len=*),intent(in) :: jobz,uplo
845 !arrays
846  real(dp),intent(out) :: w(n)
847  complex(dpc),intent(inout) :: ap(n*(n+1)/2)
848  complex(dpc),intent(out) :: z(ldz,n)
849
850 !Local variables ------------------------------
851 !scalars
852  integer :: info,nprocs
853  logical :: use_scalapack
854  character(len=500) :: msg
855 !arrays
856  real(dp),allocatable :: rwork(:)
857  complex(dpc),allocatable :: work(:)
858 #ifdef HAVE_LINALG_SCALAPACK
859  integer :: ierr,istwf_k,tbloc
860  logical :: want_eigenvectors
861  type(matrix_scalapack)    :: Slk_mat,Slk_vec
862  type(processor_scalapack) :: Slk_processor
863 #endif
864
865 !************************************************************************
866
867  use_scalapack=.FALSE.
868  if (PRESENT(comm)) then
869   nprocs = xmpi_comm_size(comm)
870 #ifdef HAVE_LINALG_SCALAPACK
871   use_scalapack = (nprocs>1)
872 #endif
873  end if
874
875  SELECT CASE(use_scalapack)
876
877  CASE (.FALSE.)
878
879   ABI_MALLOC(work, (MAX(1,2*n-1)))
880   ABI_MALLOC(rwork, (MAX(1,3*n-2)))
881
882   call ZHPEV(jobz,uplo,n,ap,w,z,ldz,work,rwork,info)
883
884   if (info < 0) then
885    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZHPEV had an illegal value."
886    MSG_ERROR(msg)
887   end if
888
889   if (info > 0) then
890    write(msg,'(2a,i0,a)')&
891 &   "ZHPEV: the algorithm failed to converge; ",ch10,&
892 &   info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
893    MSG_ERROR(msg)
894   end if
895
896   ABI_FREE(rwork)
897   ABI_FREE(work)
898   RETURN
899
900  CASE (.TRUE.)
901
902 #ifdef HAVE_LINALG_SCALAPACK
903    call init_scalapack(Slk_processor,comm)
904    istwf_k=1
905
906    ! Initialize and fill Scalapack matrix from the global one.
907    tbloc=SLK_BLOCK_SIZE
908    call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
909
910    write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
911    call wrtout(std_out,msg,"COLL")
912
913    call slk_matrix_from_global_dpc_1Dp(Slk_mat,uplo,ap)
914
915    want_eigenvectors = firstchar(jobz,(/"V","v"/))
916    if (want_eigenvectors) then ! Initialize the distributed vectors.
917     call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc)
918    end if
919
920    ! Solve the problem with scaLAPACK.
921    call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
922
923    call destruction_matrix_scalapack(Slk_mat)
924
925    if (want_eigenvectors) then ! Collect the eigenvectors.
926     z = zero
927     call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
928     call destruction_matrix_scalapack(Slk_vec)
929     call xmpi_sum(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
930    end if
931
932    call end_scalapack(Slk_processor)
933
934   RETURN
935 #endif
936
937   MSG_BUG("You should not be here!")
938
939  END SELECT
940
941 end subroutine wrap_ZHPEV
```

## m_hide_lapack/zginv [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3024 subroutine zginv(a,n,comm)
3025
3026
3027 !This section has been created automatically by the script Abilint (TD).
3028 !Do not modify the following lines by hand.
3029 #undef ABI_FUNC
3030 #define ABI_FUNC 'zginv'
3031 !End of the abilint section
3032
3033  implicit none
3034
3035 !Arguments ------------------------------------
3036 !scalars
3037  integer,intent(in) :: n
3038  integer,optional,intent(in) :: comm
3039 !arrays
3040  complex(dpc),intent(inout) :: a(n,n)
3041
3042 !Local variables-------------------------------
3043 !scalars
3044  integer :: lwork,info,nprocs
3045  logical :: use_scalapack
3046  character(len=500) :: msg
3047 !arrays
3048  integer,allocatable :: ipiv(:)
3049  complex(dpc),allocatable :: work(:)
3050 #ifdef HAVE_LINALG_SCALAPACK
3051  integer :: istwf_k,tbloc,ierr
3052  type(matrix_scalapack)    :: Slk_mat
3053  type(processor_scalapack) :: Slk_processor
3054 #endif
3055
3056 ! *************************************************************************
3057
3058  use_scalapack=.FALSE.
3059  if (PRESENT(comm)) then
3060   nprocs = xmpi_comm_size(comm)
3061 #ifdef HAVE_LINALG_SCALAPACK
3062   use_scalapack = (nprocs>1)
3063 #endif
3064  end if
3065
3066  SELECT CASE(use_scalapack)
3067
3068  CASE (.FALSE.)
3069
3070   ABI_MALLOC(ipiv,(n))
3071   call ZGETRF(n,n,a,n,ipiv,info) ! P* L* U  Factorization.
3072
3073   if (info < 0) then
3074    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRF had an illegal value."
3075    MSG_ERROR(msg)
3076   end if
3077
3078   if (info > 0) then
3079    write(msg,'(3a,i0,4a)')&
3080 &   "The matrix that has been passed in argument is probably either singular or nearly singular.",ch10,&
3081 &   "U(i,i) in the P*L*U factorization is exactly zero for i = ",info,ch10,&
3082 &   "The factorization has been completed but the factor U is exactly singular.",ch10,&
3083 &   "Division by zero will occur if it is used to solve a system of equations."
3084    MSG_ERROR(msg)
3085   end if
3086
3087   lwork=MAX(1,n)
3088
3089   ABI_MALLOC(work,(lwork))
3090
3091   call ZGETRI(n,a,n,ipiv,work,lwork,info) ! Invert U and then compute inv(A)
3092
3093   if (info < 0) then
3094    write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZGETRI had an illegal value."
3095    MSG_ERROR(msg)
3096   end if
3097
3098   if (info > 0) then
3099    write(msg,'(3a,i0,a)')&
3100 &   "The matrix that has been passed to this subroutine is probably either singular or nearly singular.",ch10,&
3101 &   "U(i,i) for i= ",info," is exactly zero; the matrix is singular and its inverse could not be computed."
3102    MSG_ERROR(msg)
3103   end if
3104
3105   ABI_FREE(ipiv)
3106   ABI_FREE(work)
3107
3108   RETURN
3109
3110  CASE (.TRUE.)
3111
3112 #ifdef HAVE_LINALG_SCALAPACK
3113   call init_scalapack(Slk_processor,comm)
3114   istwf_k=1
3115
3116   ! Initialize and fill Scalapack matrix from the global one.
3117   tbloc=SLK_BLOCK_SIZE
3118   call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
3119
3120   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
3121   call wrtout(std_out,msg,"COLL")
3122
3123   call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
3124
3125   ! Perform the calculation with scaLAPACK.
3126   call slk_zinvert(Slk_mat)
3127
3128   ! Reconstruct the global matrix from the distributed one.
3129   a = czero
3130   call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a)  ! Fill the entries calculated by this node.
3131   call destruction_matrix_scalapack(Slk_mat)
3132
3133   call xmpi_sum(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
3134   call end_scalapack(Slk_processor)
3135
3136   RETURN
3137 #endif
3138
3139   MSG_BUG("You should not be here!")
3140
3141  END SELECT
3142
3143 end subroutine zginv
```

## m_hide_lapack/zhpd_invert [ Functions ]

[ Top ] [ m_hide_lapack ] [ 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

```3182 subroutine zhpd_invert(uplo,a,n,comm)
3183
3184
3185 !This section has been created automatically by the script Abilint (TD).
3186 !Do not modify the following lines by hand.
3187 #undef ABI_FUNC
3188 #define ABI_FUNC 'zhpd_invert'
3189 !End of the abilint section
3190
3191  implicit none
3192
3193 !Arguments ------------------------------------
3194 !scalars
3195  character(len=*),intent(in) :: uplo
3196  integer,intent(in) :: n
3197  integer,optional,intent(in) :: comm
3198 !arrays
3199  complex(dpc),intent(inout) :: a(n,n)
3200
3201 !Local variables-------------------------------
3202 !scalars
3203  integer :: info,nprocs
3204  logical :: use_scalapack
3205  character(len=500) :: msg
3206 !arrays
3207 #ifdef HAVE_LINALG_SCALAPACK
3208  integer :: istwf_k,tbloc,ierr
3209  type(matrix_scalapack)    :: Slk_mat
3210  type(processor_scalapack) :: Slk_processor
3211 #endif
3212
3213 ! *************************************************************************
3214
3215  use_scalapack=.FALSE.
3216  if (PRESENT(comm)) then
3217   nprocs = xmpi_comm_size(comm)
3218 #ifdef HAVE_LINALG_SCALAPACK
3219   use_scalapack = (nprocs>1)
3220 #endif
3221  end if
3222
3223  SELECT CASE(use_scalapack)
3224
3225  CASE (.FALSE.)
3226    ! *  ZPOTRF computes the Cholesky factorization of a complex Hermitian positive definite.
3227    ! *     A = U**H * U,  if UPLO = 'U', or
3228    ! *     A = L  * L**H,  if UPLO = 'L',
3229    call ZPOTRF(uplo,n,a,n,info)
3230
3231    if (info < 0) then
3232      write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRF had an illegal value."
3233      MSG_ERROR(msg)
3234    end if
3235
3236    if (info > 0) then
3237     write(msg,'(a,i0,3a)')&
3238 &     "The leading minor of order ",info," is not positive definite, ",ch10,&
3239 &     "and the factorization could not be completed."
3240     MSG_ERROR(msg)
3241    end if
3242    !
3243    ! *  ZPOTRI computes the inverse of a complex Hermitian positive definite
3244    ! *  matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
3245    ! *  computed by ZPOTRF.
3246    ! *  On exit, the upper or lower triangle of the (Hermitian)
3247    ! *  inverse of A, overwriting the input factor U or L.
3248    call ZPOTRI(uplo,n,a,n,info)
3249
3250   if (info < 0) then
3251     write(msg,'(a,i0,a)')" The ",-info,"-th argument of ZPOTRI had an illegal value."
3252     MSG_ERROR(msg)
3253   end if
3254
3255   if (info > 0) then
3256     write(msg,'(a,2(1x,i0),a)')&
3257 &     "The ( ",info,info,")element of the factor U or L is zero, and the inverse could not be computed."
3258     MSG_ERROR(msg)
3259   end if
3260
3261   RETURN
3262
3263  CASE (.TRUE.)
3264
3265 #ifdef HAVE_LINALG_SCALAPACK
3266   call init_scalapack(Slk_processor,comm)
3267   istwf_k=1
3268
3269   ! Initialize and fill Scalapack matrix from the global one.
3270   tbloc=SLK_BLOCK_SIZE
3271   call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
3272
3273   write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs = ",nprocs,"; block size = ",tbloc
3274   call wrtout(std_out,msg,"COLL")
3275
3276   call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
3277
3278   ! Perform the calculation with scaLAPACK.
3279   call slk_zdhp_invert(Slk_mat,uplo)
3280
3281   ! Reconstruct the global matrix from the distributed one.
3282   a = czero
3283   call slk_matrix_to_global_dpc_2D(Slk_mat,uplo,a)  ! Fill the entries calculated by this node.
3284   call destruction_matrix_scalapack(Slk_mat)
3285
3286   call xmpi_sum(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
3287   call end_scalapack(Slk_processor)
3288
3289   RETURN
3290 #endif
3291    MSG_BUG("You should not be here!")
3292  END SELECT
3293
3294 end subroutine zhpd_invert
```