TABLE OF CONTENTS
ABINIT/m_matrix [ Modules ]
NAME
m_matrix
FUNCTION
Module containing some function acting on a matrix (sqrt root)
COPYRIGHT
Copyright (C) 2009-2024 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
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 MODULE m_matrix 26 27 use defs_basis 28 use m_errors 29 use m_abicore 30 31 use m_hide_lapack, only : xginv 32 33 implicit none 34 35 private 36 37 ! public :: init_matrix ! Main creation method 38 public :: invsqrt_matrix ! inv of Sqrt of Matrix 39 public :: blockdiago_fordsyev ! inv of Sqrt of Matrix 40 public :: blockdiago_forzheev ! inv of Sqrt of Matrix 41 ! public :: inverse_matrix ! Inverse matrix 42 ! public :: nullify_matrix ! Nullify the object 43 ! public :: destroy_matrix ! Frees the allocated memory 44 ! public :: print_matrix ! Printout of the basic info 45 46 47 CONTAINS !===========================================================
FUNCTION
Initialize matrix
INPUTS
ndim = dimension of matrix matrix= matrix
OUTPUT
matrix= square root of the matrix
force_diag = 0 if it no 0 on diagonal
= nb of zeros found otherwise
SOURCE
63 subroutine invsqrt_matrix(matrix,tndim,force_diag) 64 65 !Arguments ------------------------------------ 66 !scalars 67 integer,intent(in) :: tndim 68 complex(dpc),intent(inout) :: matrix(tndim,tndim) 69 integer, intent(out) :: force_diag 70 !arrays 71 72 !Local variables------------------------------- 73 !scalars 74 integer :: im,im1,im2,info,lwork,nb_of_zero 75 character(len=500) :: message 76 real(dp) :: pawprtvol 77 !arrays 78 real(dp),allocatable :: eig(:),rwork(:) 79 complex(dpc),allocatable :: zwork(:),diag(:,:) 80 complex(dpc),allocatable :: sqrtmat(:,:),zhdp2(:,:),sqrtmatinv(:,:) 81 complex(dpc),allocatable :: initialmatrix(:,:) 82 83 ! ************************************************************************* 84 85 !Do not remove this silly print instruction. Seems needed to avoid floating 86 !point exception on vm1_gcc51 ... 87 #if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2) 88 write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )' 89 #endif 90 91 DBG_ENTER("COLL") 92 pawprtvol=2 93 94 ABI_MALLOC(initialmatrix,(tndim,tndim)) 95 initialmatrix=matrix 96 ! == First diagonalize matrix and keep the matrix for the change of basis 97 lwork=2*tndim-1 98 ABI_MALLOC(rwork,(3*tndim-2)) 99 ABI_MALLOC(zwork,(lwork)) 100 ABI_MALLOC(eig,(tndim)) 101 102 call zheev('v','u',tndim,matrix,tndim,eig,zwork,lwork,rwork,info) 103 if(pawprtvol>3) then 104 write(message,'(2a)') ch10,' - rotation matrix - ' 105 call wrtout(std_out,message,'COLL') 106 do im1=1,tndim 107 write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')& 108 ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')& 109 & (matrix(im1,im2),im2=1,tndim) 110 call wrtout(std_out,message,'COLL') 111 end do 112 endif 113 114 115 ABI_FREE(zwork) 116 ABI_FREE(rwork) 117 if(info/=0) then 118 message = 'Error in diagonalization of zmat (zheev) ! - ' 119 ABI_ERROR(message) 120 end if 121 122 ! == Secondly Compute sqrt(diagonalized matrix) 123 ABI_MALLOC(diag,(tndim,tndim)) 124 diag=czero 125 nb_of_zero=0 126 do im=1,tndim 127 128 if(eig(im)< -tol8) then 129 message = " - Eigenvalues from zheev are negative or zero ! - " 130 write(std_out,*) 131 write(std_out,*) " Eigenvalue=",eig(im) 132 write(std_out,*) " Matrix is" 133 do im1=1,tndim 134 write(std_out,'(100f7.3)') (initialmatrix(im1,im2),im2=1,tndim) 135 enddo 136 ABI_ERROR(message) 137 else if(abs(eig(im))<tol8) then 138 nb_of_zero=nb_of_zero+1 139 else 140 diag(im,im)=cmplx(one/sqrt(eig(im)),zero,kind=dp) 141 endif 142 enddo 143 force_diag=nb_of_zero 144 ABI_FREE(eig) 145 ! write(std_out,*) "sqrt(eig) , diag(1,1)",sqrt(eig(1)),diag(1,1) 146 ! write(std_out,*) "cmplx(sqrt(eig(1)),zero,dp) , diag(1,1)",cmplx(sqrt(eig(1)),zero,dp),diag(1,1) 147 ! write(std_out,*) "sqrt(cmplx(eig(1),zero,dp)) , diag(1,1)",sqrt(cmplx(eig(1),zero,dp)),diag(1,1) 148 149 ! == Thirdly Multiply by matrix for the change of basis 150 ABI_MALLOC(sqrtmat,(tndim,tndim)) 151 ABI_MALLOC(zhdp2,(tndim,tndim)) 152 if(pawprtvol>3) then 153 write(message,'(2a)') ch10,' - 1.0/sqrt(Eigenmatrix) - ' 154 call wrtout(std_out,message,'COLL') 155 do im1=1,tndim 156 write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')& 157 ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')& 158 & (diag(im1,im2),im2=1,tndim) 159 call wrtout(std_out,message,'COLL') 160 end do 161 endif 162 !zgemm(A,B,C) : C = op(A) op(B) 163 call zgemm('n','t',tndim,tndim,tndim,cone,diag,tndim,conjg(matrix),tndim,czero,zhdp2,tndim) 164 call zgemm('n','n',tndim,tndim,tndim,cone,matrix,tndim,zhdp2,tndim,czero,sqrtmat,tndim) 165 ! if(abs(pawprtvol)>=3) then 166 if(pawprtvol>3) then 167 write(message,'(3a)') ch10," - inverse Sqrt root of matrix is - " 168 call wrtout(std_out,message,'COLL') 169 do im1=1,tndim 170 write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')& 171 & (sqrtmat(im1,im2),im2=1,tndim) 172 call wrtout(std_out,message,'COLL') 173 end do 174 endif 175 ! endif 176 ABI_FREE(diag) 177 178 ! == Forthly Compute the inverse of the square root 179 ! call matcginv_dpc(sqrtmat,tndim,tndim) 180 !call xginv(sqrtmat,tndim) 181 ABI_MALLOC(sqrtmatinv,(tndim,tndim)) 182 sqrtmatinv=sqrtmat 183 if(pawprtvol>3) then 184 write(message,'(2a)') ch10," - inverse Sqrt root of matrix is - " 185 call wrtout(std_out,message,'COLL') 186 do im1=1,tndim 187 write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')& 188 & (sqrtmatinv(im1,im2),im2=1,tndim) 189 call wrtout(std_out,message,'COLL') 190 end do 191 endif 192 ABI_FREE(sqrtmat) 193 194 ! == Fifthly Check that O^{-0/5} O O{-0/5}=I 195 ! zgemm(A,B,C) : C = op(A) op(B) 196 call zgemm('n','n',tndim,tndim,tndim,cone,initialmatrix,tndim,sqrtmatinv,tndim,czero,zhdp2,tndim) 197 call zgemm('n','n',tndim,tndim,tndim,cone,sqrtmatinv,tndim,zhdp2,tndim,czero,initialmatrix,tndim) 198 if(pawprtvol>3) then 199 write(message,'(3a)') ch10," - O^{-0/5} O O^{-0/5}=I - " 200 call wrtout(std_out,message,'COLL') 201 do im1=1,tndim 202 write(message,'(12(1x,18(1x,"(",f10.6,",",f4.1,")")))')& 203 ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')& 204 & (initialmatrix(im1,im2),im2=1,tndim) 205 call wrtout(std_out,message,'COLL') 206 end do 207 endif 208 ABI_FREE(zhdp2) 209 matrix=sqrtmatinv 210 ABI_FREE(sqrtmatinv) 211 ABI_FREE(initialmatrix) 212 213 DBG_EXIT("COLL") 214 215 end subroutine invsqrt_matrix