TABLE OF CONTENTS


ABINIT/m_matrix [ Modules ]

[ Top ] [ Modules ]

NAME

 m_matrix

FUNCTION

 Module containing some function acting on a matrix 
  (sqrt root)

COPYRIGHT

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

NOTES

PARENTS

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 MODULE m_matrix
31 
32  use defs_basis
33  use m_errors
34  use m_abicore
35 
36  use m_hide_lapack,  only : xginv
37 
38  implicit none
39 
40  private
41 
42 ! public :: init_matrix         ! Main creation method
43  public :: invsqrt_matrix         ! inv of Sqrt of Matrix
44  public :: blockdiago_fordsyev         ! inv of Sqrt of Matrix
45  public :: blockdiago_forzheev         ! inv of Sqrt of Matrix
46 ! public :: inverse_matrix      ! Inverse matrix
47 ! public :: nullify_matrix      ! Nullify the object
48 ! public :: destroy_matrix      ! Frees the allocated memory
49 ! public :: print_matrix        ! Printout of the basic info
50 
51 
52 CONTAINS  !===========================================================

FUNCTION

  Initialize matrix

INPUTS

  ndim = dimension of matrix
  matrix= matrix

OUTPUT

  matrix= square root of the matrix

PARENTS

CHILDREN

SOURCE

 70 subroutine invsqrt_matrix(matrix,tndim)
 71 
 72 
 73 !This section has been created automatically by the script Abilint (TD).
 74 !Do not modify the following lines by hand.
 75 #undef ABI_FUNC
 76 #define ABI_FUNC 'invsqrt_matrix'
 77 !End of the abilint section
 78 
 79  implicit none
 80 
 81 !Arguments ------------------------------------
 82 !scalars
 83  integer,intent(in) :: tndim 
 84  complex(dpc),intent(inout) :: matrix(tndim,tndim)
 85 !arrays
 86 
 87 !Local variables-------------------------------
 88 !scalars
 89  integer :: im,im1,im2,info,lwork
 90  character(len=500) :: message
 91  real(dp) :: pawprtvol
 92 !arrays
 93  real(dp),allocatable :: eig(:),rwork(:)
 94  complex(dpc),allocatable :: zwork(:),diag(:,:)
 95  complex(dpc),allocatable :: sqrtmat(:,:),zhdp2(:,:),sqrtmatinv(:,:)
 96  complex(dpc),allocatable :: initialmatrix(:,:)
 97  
 98 ! *************************************************************************
 99 
100 !Do not remove this silly print instruction. Seems needed to avoid floating
101 !point exception on vm1_gcc51 ...
102 #if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
103  write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
104 #endif
105 
106  DBG_ENTER("COLL")
107  pawprtvol=2
108 
109  ABI_ALLOCATE(initialmatrix,(tndim,tndim))
110  initialmatrix=matrix
111 !  == First diagonalize matrix and keep the matrix for the change of basis
112  lwork=2*tndim-1
113  ABI_ALLOCATE(rwork,(3*tndim-2))
114  ABI_ALLOCATE(zwork,(lwork))
115  ABI_ALLOCATE(eig,(tndim))
116  
117  call zheev('v','u',tndim,matrix,tndim,eig,zwork,lwork,rwork,info)
118  
119  
120  ABI_DEALLOCATE(zwork)
121  ABI_DEALLOCATE(rwork)
122  if(info/=0) then
123   message = 'Error in diagonalization of zmat (zheev) ! - '
124   MSG_ERROR(message)
125  end if
126 
127 !  == Secondly Compute sqrt(diagonalized matrix)
128  ABI_ALLOCATE(diag,(tndim,tndim))
129  diag=czero
130  do im=1,tndim
131 
132    if(eig(im)<tol16) then
133      message = "  - Eigenvalues from zheev are negative or zero ! - "
134      MSG_ERROR(message)
135    else
136      diag(im,im)=cmplx(sqrt(eig(im)),zero,kind=dp)
137    endif
138  enddo
139  ABI_DEALLOCATE(eig)
140 ! write(std_out,*) "sqrt(eig)                , diag(1,1)",sqrt(eig(1)),diag(1,1)
141 ! write(std_out,*) "cmplx(sqrt(eig(1)),zero,dp) , diag(1,1)",cmplx(sqrt(eig(1)),zero,dp),diag(1,1)
142 ! write(std_out,*) "sqrt(cmplx(eig(1),zero,dp)) , diag(1,1)",sqrt(cmplx(eig(1),zero,dp)),diag(1,1)
143 
144 !  == Thirdly Multiply by  matrix for the change of basis
145  ABI_ALLOCATE(sqrtmat,(tndim,tndim))
146  ABI_ALLOCATE(zhdp2,(tndim,tndim))
147  if(pawprtvol>3) then
148    write(message,'(2a)') ch10,'  - sqrt(Eigenmatrix) - '
149    call wrtout(std_out,message,'COLL')
150    do im1=1,tndim
151      write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
152 !     write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
153 &     (diag(im1,im2),im2=1,tndim)
154     call wrtout(std_out,message,'COLL')
155    end do
156  endif
157 !zgemm(A,B,C) : C = op(A) op(B)
158  call zgemm('n','t',tndim,tndim,tndim,cone,diag,tndim,conjg(matrix),tndim,czero,zhdp2,tndim)
159  call zgemm('n','n',tndim,tndim,tndim,cone,matrix,tndim,zhdp2,tndim,czero,sqrtmat,tndim)
160 ! if(abs(pawprtvol)>=3) then
161  if(pawprtvol>3) then
162    write(message,'(3a)') ch10,"  - Sqrt root of matrix is - "
163    call wrtout(std_out,message,'COLL')
164    do im1=1,tndim
165      write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
166 &     (sqrtmat(im1,im2),im2=1,tndim)
167      call wrtout(std_out,message,'COLL')
168    end do
169  endif
170 ! endif
171  ABI_DEALLOCATE(diag)
172 
173 !  == Forthly Compute the inverse of the square root
174 ! call matcginv_dpc(sqrtmat,tndim,tndim)
175  call xginv(sqrtmat,tndim)
176  ABI_ALLOCATE(sqrtmatinv,(tndim,tndim))
177  sqrtmatinv=sqrtmat
178  if(pawprtvol>3) then
179    write(message,'(2a)') ch10,"  - inverse Sqrt root of matrix is - "
180    call wrtout(std_out,message,'COLL')
181    do im1=1,tndim
182      write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
183 &     (sqrtmatinv(im1,im2),im2=1,tndim)
184      call wrtout(std_out,message,'COLL')
185    end do
186  endif
187  ABI_DEALLOCATE(sqrtmat)
188 
189 !  == Fifthly Check that O^{-0/5} O O{-0/5}=I
190 !  zgemm(A,B,C) : C = op(A) op(B)
191  call zgemm('n','n',tndim,tndim,tndim,cone,initialmatrix,tndim,sqrtmatinv,tndim,czero,zhdp2,tndim)
192  call zgemm('n','n',tndim,tndim,tndim,cone,sqrtmatinv,tndim,zhdp2,tndim,czero,initialmatrix,tndim)
193  if(pawprtvol>3) then
194    write(message,'(3a)') ch10,"  - O^{-0/5} O O^{-0/5}=I - "
195    call wrtout(std_out,message,'COLL')
196    do im1=1,tndim
197      write(message,'(12(1x,18(1x,"(",f10.6,",",f4.1,")")))')&
198 !     write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
199 &     (initialmatrix(im1,im2),im2=1,tndim)
200      call wrtout(std_out,message,'COLL')
201    end do
202  endif
203  ABI_DEALLOCATE(zhdp2)
204  matrix=sqrtmatinv
205  ABI_DEALLOCATE(sqrtmatinv)
206  ABI_DEALLOCATE(initialmatrix)
207 
208  DBG_EXIT("COLL")
209 
210 end subroutine invsqrt_matrix