TABLE OF CONTENTS


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

 Copyright (C) 1992-2018 ABINIT group (MG, GMR, XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

PARENTS

CHILDREN

TODO

  1) Use a function to define the size of the Scalapack block according to some heuristic method.
  2) Define a threshold below which Scalapack is not used although the MPI communicator is passed.
  3) Split MPI communicator for Scalapack (.i.e. use a max size for the Scalapack comm; see abi_linalg_init).
  4) On certain networks, xmpi_sum might crash due to the size of the MPI packet.
     This problem should be solved in hide_mpi (Module containing a private global variable
     defining a threshold above which the input array is split into smaller chunks.

SOURCE

 45 #if defined HAVE_CONFIG_H
 46 #include "config.h"
 47 #endif
 48 
 49 #include "abi_common.h"
 50 
 51 MODULE m_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.

  See also SIDE EFFECTS

SIDE EFFECTS

   A       (input/output) COMPLEX(SPC) array, dimension (LDA,N)
           On entry, the N-by-N matrix A.
           On exit, A has been overwritten.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) COMPLEX(SPC) array, dimension (N, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
          orthonormal eigenvectors of the matrix A.
          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
          or the upper triangle (if UPLO='U') of A, including the
          diagonal, is destroyed.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

217 subroutine wrap_CHEEV(jobz,uplo,n,a,w)
218 
219 
220 !This section has been created automatically by the script Abilint (TD).
221 !Do not modify the following lines by hand.
222 #undef ABI_FUNC
223 #define ABI_FUNC 'wrap_CHEEV'
224 !End of the abilint section
225 
226  implicit none
227 
228 !Arguments ------------------------------------
229 !scalars
230  integer,intent(in) :: n
231  character(len=*),intent(in) :: jobz,uplo
232 !scalars
233  real(sp),intent(out) :: w(n)
234  complex(spc),intent(inout) :: a(n,n)
235 
236 !Local variables ------------------------------
237 !scalars
238  integer :: lwork,info
239  character(len=500) :: msg
240 !arrays
241  real(sp),allocatable :: rwork(:)
242  complex(spc),allocatable :: work(:)
243 
244 !************************************************************************
245 
246  lwork = MAX(1,2*n-1)
247 
248  ABI_MALLOC(work, (lwork))
249  ABI_MALLOC(rwork, (MAX(1,3*n-2)))
250 
251  call CHEEV(jobz,uplo,n,a,n,w,work,lwork,rwork,info)
252 
253  if (info < 0) then
254   write(msg,'(a,i0,a)')"The ",-info,"-th argument of CHEEV had an illegal value."
255   MSG_ERROR(msg)
256  end if
257 
258  if (info > 0) then
259   write(msg,'(2a,i0,a)')&
260 &  "CHEEV: the algorithm failed to converge; ",ch10,&
261 &  info," off-diagonal elements of an intermediate tridiagonal form did not converge to zero. "
262   MSG_ERROR(msg)
263  end if
264 
265  ABI_FREE(rwork)
266  ABI_FREE(work)
267 
268  !TODO scaLAPACK version (complex single precision buffer is needed in matrix_scalapack)
269 
270 end subroutine wrap_CHEEV

m_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.

 See also SIDE EFFECTS

SIDE EFFECTS

  AP      (input/output) COMPLEX(SPC) array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the Hermitian matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

          On exit, AP is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
          and first superdiagonal of the tridiagonal matrix T overwrite
          the corresponding elements of A, and if UPLO = 'L', the
          diagonal and first subdiagonal of T overwrite the
          corresponding elements of A.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) REAL(DP) array, dimension (CPLEX, N, N)
          On entry, the (complex Hermitian|Real symmetric) matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
          orthonormal eigenvectors of the matrix A.
          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
          or the upper triangle (if UPLO='U') of A, including the
          diagonal, is destroyed.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) REAL(DP) array, dimension (CPLEX, N, N)
          On entry, the (real symmetric|complex Hermitian) matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) REAL(DP) array, dimension (CPLEX,N, N)
          On  entry, the (real symmetric|Hermitian) matrix A.  If UPLO = "U", the leading N-by-N upper triangular part of A
          <S-F1>contains the upper triangular part of the matrix A.
          If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.

          On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
          The eigenvectors are normalized as follows:
          if ITYPE =1 or 2:
            Z**T*B*Z = I if CPLEX=1
            Z**H*B*Z = I if CPLEX=2.
          if ITYPE = 3,
             Z**T*inv(B)*Z = I if CPLEX=1
             Z**H*inv(B)*Z = I if CPLEX=2

          If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
          (if UPLO="L") of A, including the diagonal, is destroyed.


  B       (input/output) REAL(DP) array, dimension (CPLEX,N, N)
          On entry, the (real symmetric|Hermitian) positive definite matrix B.
          If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
          If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

          On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
          factor U or L from the Cholesky factorization
          B = U**T*U or B = L*L**T if CPLEX=1
          B = U**H*U or B = L*L**H if CPLEX=2

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) REAL(DP) array, dimension (CPLEX, N, N)
          On entry, the (real symmetric| complex Hermitian) matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = "L",
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.

          On exit, the lower triangle (if UPLO="L") or the upper
          triangle (if UPLO="U") of A, including the diagonal, is
          destroyed.

   B      (input/output) REAL(DP) array, dimension (CPLEX, LDB, N)
          On entry, the (real symmetric| complex Hermitian) matrix B.  If UPLO = "U", the leading N-by-N upper triangular part
          of B contains the upper triangular part  of the matrix B.
          If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

          On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
          U or L from the Cholesky factorization B = U**H*U or B = L*L**H.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

  See also SIDE EFFECTS

SIDE EFFECTS

   A       (input/output) COMPLEX(DPC) array, dimension (LDA,N)
           On entry, the N-by-N matrix A.
           On exit, A has been overwritten.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) COMPLEX(DPC) array, dimension (N, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
          orthonormal eigenvectors of the matrix A.
          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
          or the upper triangle (if UPLO='U') of A, including the
          diagonal, is destroyed.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

324 subroutine wrap_ZHEEV(jobz,uplo,n,a,w,comm)
325 
326 
327 !This section has been created automatically by the script Abilint (TD).
328 !Do not modify the following lines by hand.
329 #undef ABI_FUNC
330 #define ABI_FUNC 'wrap_ZHEEV'
331 !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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) COMPLEX(DPC) array, dimension (N, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) COMPLEX(DPC) array, dimension (N, N)
          On  entry, the Hermitian matrix A.  If UPLO = "U", the leading N-by-N upper triangular part of A
          <S-F1>contains the upper triangular part of the matrix A.
          If UPLO = "L", the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A.

          On exit, if JOBZ = "V", then A contains the matrix Z of eigenvectors.
          The eigenvectors are normalized as follows: if ITYPE = 1  or  2,
          Z**H*B*Z  = I; if ITYPE = 3, Z**H*inv(B)*Z = I.
          If JOBZ = "N", then on exit the upper triangle (if UPLO="U") or the lower triangle
          (if UPLO="L") of A, including the diagonal, is destroyed.


  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
          On entry, the Hermitian positive definite matrix B.
          If UPLO = "U", the leading N-by-N upper triangular part of B contains the upper triangular part of the matrix B.
          If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

          On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular
          factor U or L from the Cholesky factorization B = U**H*U or B = L*L**H.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  A       (input/output) COMPLEX(DPC) array, dimension (N, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = "L",
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.

          On exit, the lower triangle (if UPLO="L") or the upper
          triangle (if UPLO="U") of A, including the diagonal, is
          destroyed.

   B      (input/output) COMPLEX(DPC) array, dimension (LDB, N)
          On entry, the Hermitian matrix B.  If UPLO = "U", the leading N-by-N upper triangular part
          of B contains the upper triangular part  of the matrix B.
          If UPLO = "L", the leading N-by-N lower triangular part of B contains the lower triangular part of the matrix B.

          On exit, if INFO <= N, the part of B containing the matrix is overwritten by the triangular factor
          U or L from the Cholesky factorization B = U**H*U or B = L*L**H.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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.

 See also SIDE EFFECTS

SIDE EFFECTS

  AP      (input/output) COMPLEX(DPC) array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the Hermitian matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

          On exit, AP is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
          and first superdiagonal of the tridiagonal matrix T overwrite
          the corresponding elements of A, and if UPLO = 'L', the
          diagonal and first subdiagonal of T overwrite the
          corresponding elements of A. Unchanged if ScaLAPACK is used.

PARENTS

CHILDREN

      cwtime,xginv

SOURCE

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