./32_util/hermit.F90

ABINIT/hermit

ABINIT/hermit [ Functions ]

NAME

 hermit

FUNCTION

 Take a matrix in hermitian storage mode (lower triangle stored)
 and redefine diagonal elements to impose Hermiticity
 (diagonal terms have to be real).
 If abs(Im(H(i,i)))>4096*machine precision, print error warning.
 (Typical 64 bit machine precision is 2^-52 or 2.22e-16)

 Copyright (C) 1998-2017 ABINIT group (DCA, XG, GMR).
 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 .

INPUTS

  chmin(n*n+n)=complex hermitian matrix with numerical noise possibly
   rendering Im(diagonal elements) approximately 1e-15 or so
  ndim=size of complex hermitian matrix

OUTPUT

  chmout(n*n+n)=redefined matrix with strictly real diagonal elements.
   May be same storage location as chmin.
  ierr=0 if no problem, 1 if the imaginary part of some element
   too large (at present, stop in this case).

TODO

  Name is misleading, perhaps hermit_force_diago?
  Interface allows aliasing

PARENTS

      extrapwf,subdiago

CHILDREN

      wrtout

SOURCE

 43 #if defined HAVE_CONFIG_H
 44 #include "config.h"
 45 #endif
 46 
 47 #include "abi_common.h"
 48 
 49 
 50 subroutine hermit(chmin,chmout,ierr,ndim)
 51 
 52  use defs_basis
 53  use m_errors
 54 
 55 !This section has been created automatically by the script Abilint (TD).
 56 !Do not modify the following lines by hand.
 57 #undef ABI_FUNC
 58 #define ABI_FUNC 'hermit'
 59  use interfaces_14_hidewrite
 60 !End of the abilint section
 61 
 62  implicit none
 63 
 64 !Arguments ------------------------------------
 65 !scalars
 66  integer,intent(in) :: ndim
 67  integer,intent(out) :: ierr
 68 !arrays
 69  real(dp),intent(inout) :: chmin(ndim*ndim+ndim)
 70  real(dp),intent(inout) :: chmout(ndim*ndim+ndim)
 71 
 72 !Local variables-------------------------------
 73 !scalars
 74  integer,save :: mmesgs=20,nmesgs=0
 75  integer :: idim,merrors,nerrors
 76  real(dp),parameter :: eps=epsilon(0.0d0)
 77  real(dp) :: ch_im,ch_re,moduls,tol
 78  character(len=500) :: message
 79 
 80 ! *************************************************************************
 81 
 82  tol=4096.0d0*eps
 83 
 84  ierr=0
 85  merrors=0
 86 
 87 !Copy matrix into possibly new location
 88  chmout(:)=chmin(:)
 89 
 90 !Loop over diagonal elements of matrix (off-diag not altered)
 91  do idim=1,ndim
 92 
 93    ch_im=chmout(idim*idim+idim  )
 94    ch_re=chmout(idim*idim+idim-1)
 95 
 96 !  check for large absolute Im part and print warning when
 97 !  larger than (some factor)*(machine precision)
 98    nerrors=0
 99    if( abs(ch_im) > tol .and. abs(ch_im) > tol8*abs(ch_re)) nerrors=2
100    if( abs(ch_im) > tol .or. abs(ch_im) > tol8*abs(ch_re)) nerrors=1
101 
102    if( (abs(ch_im) > tol .and. nmesgs<mmesgs) .or. nerrors==2)then
103      write(message, '(3a,i0,a,es20.12,a,es20.12,a)' )&
104 &     'Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
105 &     'for component',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
106      call wrtout(std_out,message,'PERS')
107      nmesgs=nmesgs+1
108    end if
109 
110    if( ( abs(ch_im) > tol8*abs(ch_re) .and. nmesgs<mmesgs) .or. nerrors==2)then
111      write(message, '(3a,i0,a,es20.12,a,es20.12,a)' )&
112 &     'Input Hermitian matrix has nonzero relative Im part on diagonal:',ch10,&
113 &     'for component',idim,' Im part is',ch_im,', Re part is',ch_re,'.'
114      call wrtout(std_out,message,'PERS')
115      nmesgs=nmesgs+1
116    end if
117 
118 !  compute modulus $= (\Re^2+\Im^2)^{1/2}$
119    moduls=sqrt(ch_re**2+ch_im**2)
120 
121 !  set Re part to modulus with sign of original Re part
122    chmout(idim*idim+idim-1)=sign(moduls,ch_re)
123 
124 !  set Im part to 0
125    chmout(idim*idim+idim)=zero
126 
127    merrors=max(merrors,nerrors)
128 
129  end do
130 
131  if(merrors==2)then
132    ierr=1
133    write(message, '(3a)' )&
134 &   'Imaginary part(s) of diagonal Hermitian matrix element(s) is too large.',ch10,&
135 &   'See previous messages.'
136    MSG_BUG(message)
137  end if
138 
139 end subroutine hermit

TABLE OF CONTENTS

ABINIT/hermit [ Functions ]