Branch data Line data Source code
1 : : !!****m* ABINIT/m_matrix
2 : : !! NAME
3 : : !! m_matrix
4 : : !!
5 : : !! FUNCTION
6 : : !! Module containing some function acting on a matrix (sqrt root)
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2009-2026 ABINIT group (BA, XG, MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : MODULE m_matrix
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : :
28 : : use m_numeric_tools, only : euclid_div
29 : : use m_hide_lapack, only : xginv
30 : :
31 : : implicit none
32 : :
33 : : private
34 : :
35 : : public :: invsqrt_matrix ! inv of Sqrt of Matrix
36 : : public :: blockdiago_fordsyev ! inv of Sqrt of Matrix
37 : : public :: blockdiago_forzheev ! inv of Sqrt of Matrix
38 : : public :: mat33det ! Determinant of a 3x3 matrix
39 : : public :: mati3inv ! Invert and transpose orthogonal 3x3 matrix of INTEGER elements.
40 : : public :: mati3det ! Compute the determinant of a 3x3 matrix of INTEGER elements.
41 : : public :: matr3inv ! Invert and TRANSPOSE general 3x3 matrix of real*8 elements.
42 : : public :: is_unitary ! Returns .TRUE. if the input matrix mat is unitary, i.e. mat^\dagger mat = I.
43 : : public :: is_identity ! Returns .TRUE. if the input matrix mat is the identity matrix.
44 : : public :: mati3HNF ! Compute Hermite normal form of a 3x3 matrix INTEGER elements.
45 : :
46 : : ! the determinant of a 3*3 matrix
47 : : interface mat33det
48 : : procedure real_mat33det
49 : : procedure int_mat33det
50 : : end interface mat33det
51 : :
52 : :
53 : : CONTAINS !===========================================================
54 : :
55 : : !! FUNCTION
56 : : !! Initialize matrix
57 : : !!
58 : : !! INPUTS
59 : : !! ndim = dimension of matrix
60 : : !! matrix= matrix
61 : : !!
62 : : !! OUTPUT
63 : : !! matrix= square root of the matrix
64 : : !! force_diag = 0 if it no 0 on diagonal
65 : : !! = nb of zeros found otherwise
66 : : !!
67 : : !! SOURCE
68 : :
69 : 298 : subroutine invsqrt_matrix(matrix,tndim,force_diag)
70 : :
71 : : !Arguments ------------------------------------
72 : : !scalars
73 : : integer,intent(in) :: tndim
74 : : complex(dp),intent(inout) :: matrix(tndim,tndim)
75 : : integer, intent(out) :: force_diag
76 : :
77 : : !Local variables-------------------------------
78 : : !scalars
79 : : integer :: im,im1,im2,info,lwork,nb_of_zero
80 : : character(len=500) :: message
81 : : real(dp) :: pawprtvol
82 : : !arrays
83 : 298 : real(dp),allocatable :: eig(:),rwork(:)
84 : 298 : complex(dp),allocatable :: zwork(:),diag(:,:)
85 : 298 : complex(dp),allocatable :: sqrtmat(:,:),zhdp2(:,:),sqrtmatinv(:,:)
86 : 298 : complex(dp),allocatable :: initialmatrix(:,:)
87 : : ! *************************************************************************
88 : :
89 : : !Do not remove this silly print instruction. Seems needed to avoid floating
90 : : !point exception on vm1_gcc51 ...
91 : : #if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
92 : : write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
93 : : #endif
94 : :
95 : : DBG_ENTER("COLL")
96 : 298 : pawprtvol=2
97 : :
98 [ - + + - : 1192 : ABI_MALLOC(initialmatrix,(tndim,tndim))
- + + - +
- + - - +
- + ]
99 [ - + + + : 9212 : initialmatrix=matrix
+ + ]
100 : : ! == First diagonalize matrix and keep the matrix for the change of basis
101 : 298 : lwork=2*tndim-1
102 [ + - - + : 894 : ABI_MALLOC(rwork,(3*tndim-2))
+ - + - +
- - + -
+ ]
103 [ - + + - : 894 : ABI_MALLOC(zwork,(lwork))
+ - + - -
+ - + ]
104 [ - + + - : 894 : ABI_MALLOC(eig,(tndim))
+ - + - -
+ - + ]
105 : :
106 : 298 : call zheev('v','u',tndim,matrix,tndim,eig,zwork,lwork,rwork,info)
107 : : if(pawprtvol>3) then
108 : : write(message,'(2a)') ch10,' - rotation matrix - '
109 : : call wrtout(std_out,message,'COLL')
110 : : do im1=1,tndim
111 : : write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
112 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
113 : : & (matrix(im1,im2),im2=1,tndim)
114 : : call wrtout(std_out,message,'COLL')
115 : : end do
116 : : endif
117 : :
118 : :
119 [ - + ]: 298 : ABI_FREE(zwork)
120 [ - + ]: 298 : ABI_FREE(rwork)
121 [ - + ]: 298 : if(info/=0) then
122 : 0 : message = 'Error in diagonalization of zmat (zheev) ! - '
123 : 0 : ABI_ERROR(message)
124 : : end if
125 : :
126 : : ! == Secondly Compute sqrt(diagonalized matrix)
127 [ - + + - : 894 : ABI_MALLOC(diag,(tndim,tndim))
- + + - -
+ - + ]
128 [ + + + + ]: 8914 : diag=czero
129 : : nb_of_zero=0
130 [ + + ]: 1725 : do im=1,tndim
131 : :
132 [ - + ]: 1725 : if(eig(im)< -tol8) then
133 : 0 : message = " - Eigenvalues from zheev are negative or zero ! - "
134 : 0 : write(std_out,*)
135 : 0 : write(std_out,*) " Eigenvalue=",eig(im)
136 : 0 : write(std_out,*) " Matrix is"
137 [ # # ]: 0 : do im1=1,tndim
138 : 0 : write(std_out,'(100f7.3)') (initialmatrix(im1,im2),im2=1,tndim)
139 : : enddo
140 : 0 : ABI_ERROR(message)
141 [ - + ]: 1427 : else if(abs(eig(im))<tol8) then
142 : 0 : nb_of_zero=nb_of_zero+1
143 : : else
144 : 1427 : diag(im,im)=cmplx(one/sqrt(eig(im)),zero,kind=dp)
145 : : endif
146 : : enddo
147 : 298 : force_diag=nb_of_zero
148 [ - + ]: 298 : ABI_FREE(eig)
149 : : ! write(std_out,*) "sqrt(eig) , diag(1,1)",sqrt(eig(1)),diag(1,1)
150 : : ! write(std_out,*) "cmplx(sqrt(eig(1)),zero,dp) , diag(1,1)",cmplx(sqrt(eig(1)),zero,dp),diag(1,1)
151 : : ! write(std_out,*) "sqrt(cmplx(eig(1),zero,dp)) , diag(1,1)",sqrt(cmplx(eig(1),zero,dp)),diag(1,1)
152 : :
153 : : ! == Thirdly Multiply by matrix for the change of basis
154 [ - + + - : 894 : ABI_MALLOC(sqrtmat,(tndim,tndim))
- + + - -
+ - + ]
155 [ - + + - : 894 : ABI_MALLOC(zhdp2,(tndim,tndim))
- + + - -
+ - + ]
156 : : if(pawprtvol>3) then
157 : : write(message,'(2a)') ch10,' - 1.0/sqrt(Eigenmatrix) - '
158 : : call wrtout(std_out,message,'COLL')
159 : : do im1=1,tndim
160 : : write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
161 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
162 : : & (diag(im1,im2),im2=1,tndim)
163 : : call wrtout(std_out,message,'COLL')
164 : : end do
165 : : endif
166 : : !zgemm(A,B,C) : C = op(A) op(B)
167 [ + - + + : 8914 : call zgemm('n','t',tndim,tndim,tndim,cone,diag,tndim,conjg(matrix),tndim,czero,zhdp2,tndim)
+ + - + ]
168 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,matrix,tndim,zhdp2,tndim,czero,sqrtmat,tndim)
169 : : ! if(abs(pawprtvol)>=3) then
170 : : if(pawprtvol>3) then
171 : : write(message,'(3a)') ch10," - inverse Sqrt root of matrix is - "
172 : : call wrtout(std_out,message,'COLL')
173 : : do im1=1,tndim
174 : : write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
175 : : & (sqrtmat(im1,im2),im2=1,tndim)
176 : : call wrtout(std_out,message,'COLL')
177 : : end do
178 : : endif
179 : : ! endif
180 [ - + ]: 298 : ABI_FREE(diag)
181 : :
182 : : ! == Forthly Compute the inverse of the square root
183 : : ! call matcginv_dpc(sqrtmat,tndim,tndim)
184 : : !call xginv(sqrtmat,tndim)
185 [ - + + - : 894 : ABI_MALLOC(sqrtmatinv,(tndim,tndim))
- + + - -
+ - + ]
186 [ - + + + : 9212 : sqrtmatinv=sqrtmat
+ + ]
187 : : if(pawprtvol>3) then
188 : : write(message,'(2a)') ch10," - inverse Sqrt root of matrix is - "
189 : : call wrtout(std_out,message,'COLL')
190 : : do im1=1,tndim
191 : : write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
192 : : & (sqrtmatinv(im1,im2),im2=1,tndim)
193 : : call wrtout(std_out,message,'COLL')
194 : : end do
195 : : endif
196 [ - + ]: 298 : ABI_FREE(sqrtmat)
197 : :
198 : : ! == Fifthly Check that O^{-0/5} O O{-0/5}=I
199 : : ! zgemm(A,B,C) : C = op(A) op(B)
200 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,initialmatrix,tndim,sqrtmatinv,tndim,czero,zhdp2,tndim)
201 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,sqrtmatinv,tndim,zhdp2,tndim,czero,initialmatrix,tndim)
202 : : if(pawprtvol>3) then
203 : : write(message,'(3a)') ch10," - O^{-0/5} O O^{-0/5}=I - "
204 : : call wrtout(std_out,message,'COLL')
205 : : do im1=1,tndim
206 : : write(message,'(12(1x,18(1x,"(",f10.6,",",f4.1,")")))')&
207 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
208 : : & (initialmatrix(im1,im2),im2=1,tndim)
209 : : call wrtout(std_out,message,'COLL')
210 : : end do
211 : : endif
212 [ - + ]: 298 : ABI_FREE(zhdp2)
213 [ + + + + ]: 8914 : matrix=sqrtmatinv
214 [ - + ]: 298 : ABI_FREE(sqrtmatinv)
215 [ - + ]: 298 : ABI_FREE(initialmatrix)
216 : :
217 : : DBG_EXIT("COLL")
218 : :
219 : 298 : end subroutine invsqrt_matrix
220 : : !!***
221 : :
222 : : !! FUNCTION
223 : : !! Transform matrix into block diagonal form before diagonalisation
224 : : !!
225 : : !! INPUTS
226 : : !! ndim = dimension of matrix
227 : : !! matrix= matrix
228 : : !!
229 : : !! OUTPUT
230 : : !! matrix= square root of the matrix
231 : : !!
232 : : !! SOURCE
233 : :
234 : 0 : subroutine blockdiago_fordsyev(matrix,tndim,eig)
235 : :
236 : : !Arguments ------------------------------------
237 : : !scalars
238 : : integer,intent(in) :: tndim
239 : : real(dp),intent(inout) :: matrix(tndim,tndim)
240 : : real(dp),intent(inout) :: eig(tndim)
241 : : !arrays
242 : :
243 : : !Local variables-------------------------------
244 : : !scalars
245 : : integer :: im1,im2,im3,info,lwork,im4,indice_formax,shift !im5,
246 : : character(len=500) :: message
247 : : real(dp):: tmpx,maxvalue
248 : : integer(dp):: tmpi,newstarting,current_dege,prtopt
249 : : !arrays
250 : 0 : real(dp),allocatable :: work(:)
251 : 0 : real(dp),allocatable :: Permutcol(:,:)
252 : 0 : real(dp),allocatable :: Apermutcol(:,:)
253 : 0 : real(dp),allocatable :: Apermutline(:,:)
254 : 0 : real(dp),allocatable :: Apermutlineback(:,:)
255 : 0 : real(dp),allocatable :: Permutline(:,:)
256 : 0 : real(dp),allocatable :: matrix_save(:,:) !,W(:)
257 : 0 : integer,allocatable :: nonnul(:)
258 : 0 : integer,allocatable :: nonnuldege(:)
259 : : logical :: testdege,swap
260 : : ! *************************************************************************
261 : :
262 : : !!!Do not remove this silly print instruction. Seems needed to avoid floating
263 : : !!!point exception on vm1_gcc51 ...
264 : : !!#if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
265 : : !! write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
266 : : !!#endif
267 : : DBG_ENTER("COLL")
268 : :
269 : 0 : lwork=10*tndim
270 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
271 [ # # ]: 0 : work = zero
272 : :
273 [ # # # # : 0 : ABI_MALLOC(matrix_save,(tndim,tndim))
# # # # #
# # # # #
# # ]
274 [ # # # # : 0 : matrix_save=matrix
# # ]
275 : :
276 [ # # # # : 0 : ABI_MALLOC(Permutcol,(tndim,tndim))
# # # # #
# # # ]
277 : :
278 [ # # # # ]: 0 : Permutcol=zero
279 [ # # ]: 0 : do im1=1,tndim
280 : 0 : Permutcol(im1,im1)=1.d0
281 : : end do
282 : :
283 : 0 : prtopt=0
284 : :
285 [ # # # # : 0 : ABI_MALLOC(nonnul,(tndim))
# # # # #
# # # ]
286 [ # # ]: 0 : do im1=1,tndim
287 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
288 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
289 [ # # ]: 0 : do im2=1,tndim
290 [ # # # # ]: 0 : if (abs(matrix(im1,im2))>0.000000000001.and.im2>nonnul(im1)) then
291 : 0 : nonnul(im1)=nonnul(im1)+1
292 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
293 : : ! permute
294 [ # # ]: 0 : do im3=1,tndim
295 : 0 : tmpx=matrix(im3,im2)
296 : 0 : matrix(im3,im2)=matrix(im3,nonnul(im1))
297 : 0 : matrix(im3,nonnul(im1))=tmpx
298 : 0 : tmpi=Permutcol(im3,im2)
299 : 0 : Permutcol(im3,im2)=Permutcol(im3,nonnul(im1))
300 : 0 : Permutcol(im3,nonnul(im1))=tmpi
301 : : enddo
302 [ # # ]: 0 : elseif (abs(matrix(im1,im2))<0.000000000001) then
303 : 0 : matrix(im1,im2)=zero
304 : : endif
305 : : enddo
306 : : enddo
307 : : if(prtopt==1) then
308 : : write(std_out,*) "MATRIX AFTER COLUMN PERMUT"
309 : : do im1=1,tndim
310 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
311 : : end do
312 : : write(std_out,*) "Permutcol MATRIX AFTER"
313 : : do im1=1,tndim
314 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Permutcol(im1,im2),im2=1,tndim)
315 : : end do
316 : : endif
317 : :
318 [ # # # # : 0 : ABI_MALLOC(Apermutcol,(tndim,tndim))
# # # # #
# # # ]
319 : : if(prtopt==1) then
320 : : write(std_out,*) "Check product of original matrix by permutation matrix "
321 : : endif
322 [ # # # # ]: 0 : Apermutcol=zero
323 [ # # ]: 0 : do im1=1,tndim
324 [ # # ]: 0 : do im2=1,tndim
325 : 0 : Apermutcol(im1,im2)=zero
326 [ # # ]: 0 : do im3=1,tndim
327 : 0 : Apermutcol(im1,im2)=matrix_save(im1,im3)*Permutcol(im3,im2)+Apermutcol(im1,im2)
328 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
329 : : end do
330 : : end do
331 : : end do
332 : : if(prtopt==1) then
333 : : write(std_out,*) "Asave*Permutcol"
334 : : do im1=1,tndim
335 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcol(im1,im2),im2=1,tndim)
336 : : end do
337 : : endif
338 : :
339 [ # # # # : 0 : ABI_MALLOC(Permutline,(tndim,tndim))
# # # # #
# # # ]
340 [ # # # # ]: 0 : Permutline=zero
341 [ # # ]: 0 : do im1=1,tndim
342 : 0 : Permutline(im1,im1)=1.d0
343 : : end do
344 : :
345 [ # # ]: 0 : do im1=1,tndim
346 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
347 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
348 [ # # ]: 0 : do im2=1,tndim
349 : : ! write(std_out,*) "im1,im2, abs matrix(im2,im1),nonnul(im1)",im1,im2,abs(B(im1,im2)),nonnul(im1)
350 [ # # # # ]: 0 : if (abs(matrix(im2,im1))>0.000000000001.and.im2>nonnul(im1)) then
351 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
352 : 0 : nonnul(im1)=nonnul(im1)+1
353 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
354 : : ! permute
355 [ # # ]: 0 : do im3=1,tndim
356 : 0 : tmpx=matrix(im2,im3)
357 : 0 : matrix(im2,im3)=matrix(nonnul(im1),im3)
358 : 0 : matrix(nonnul(im1),im3)=tmpx
359 : 0 : tmpi=Permutline(im2,im3)
360 : 0 : Permutline(im2,im3)=Permutline(nonnul(im1),im3)
361 : 0 : Permutline(nonnul(im1),im3)=tmpi
362 : : enddo
363 [ # # ]: 0 : elseif (abs(matrix(im2,im1))<0.000000000001) then
364 : 0 : matrix(im2,im1)=zero
365 : : endif
366 : : enddo
367 : : enddo
368 : : if(prtopt==1) then
369 : : write(std_out,*) "matrix AFTER"
370 : : do im1=1,tndim
371 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
372 : : end do
373 : : write(std_out,*) "Permutline MATRIX AFTER"
374 : : do im1=1,tndim
375 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Permutline(im1,im2),im2=1,tndim)
376 : : end do
377 : : endif
378 : :
379 : : if(prtopt==1) then
380 : : write(std_out,*) "Check product of Apermutcol matrix by permutation matrix of the line "
381 : : endif
382 [ # # # # : 0 : ABI_MALLOC(Apermutline,(tndim,tndim))
# # # # #
# # # ]
383 [ # # # # ]: 0 : Apermutline=zero
384 [ # # ]: 0 : do im1=1,tndim
385 [ # # ]: 0 : do im2=1,tndim
386 : 0 : Apermutline(im1,im2)=zero
387 [ # # ]: 0 : do im3=1,tndim
388 : 0 : Apermutline(im1,im2)=Apermutcol(im3,im2)*Permutline(im1,im3)+Apermutline(im1,im2)
389 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
390 : : end do
391 : : end do
392 : : end do
393 : : if(prtopt==1) then
394 : : write(std_out,*) "Permutline*Apermutcol"
395 : : do im1=1,tndim
396 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutline(im1,im2),im2=1,tndim)
397 : : end do
398 : : endif
399 [ # # ]: 0 : work=0.d0
400 : 0 : call dsyev('v','u',tndim,matrix_save,tndim,eig,work,lwork,info)
401 [ # # ]: 0 : if(info/=0) then
402 : 0 : message = 'Error in diagonalization of matrix (dsyev) ! - '
403 : 0 : ABI_ERROR(message)
404 : : end if
405 : : if(prtopt==1) then
406 : : write(std_out,*) 'output',INFO
407 : : write(std_out,*) "Eigenvalues"
408 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
409 : : write(std_out,*) "Eigenvectors"
410 : : do im1=1,tndim
411 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix_save(im1,im2),im2=1,tndim)
412 : : end do
413 : : endif
414 : :
415 : :
416 : : ! call dsyev('v','u',tndim,A,LDA,W,WORKTMP,LWORK,INFO)
417 : : ! write(std_out,*) "optimal lwork",worktmp(1)
418 : :
419 [ # # ]: 0 : work=0.d0
420 : 0 : call dsyev('v','u',tndim,matrix,tndim,eig,work,lwork,info)
421 [ # # ]: 0 : if(info/=0) then
422 : 0 : message = 'Error in diagonalization of matrix (dsyev) ! - '
423 : 0 : ABI_ERROR(message)
424 : : end if
425 : : if(prtopt==1) then
426 : : write(std_out,*) 'output',INFO
427 : : write(std_out,*) "Eigenvalues"
428 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
429 : : write(std_out,*) "Eigenvectors"
430 : : do im1=1,tndim
431 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
432 : : end do
433 : : endif
434 : :
435 : :
436 : : !! REORDER EIGENVECTORS
437 [ # # # # : 0 : ABI_MALLOC(nonnuldege,(tndim))
# # # # ]
438 : 0 : newstarting=1
439 : 0 : current_dege=1
440 [ # # ]: 0 : do im4=2,tndim
441 [ # # ]: 0 : if(im4<tndim) testdege=((eig(im4)-eig(im4-1))<tol12)
442 [ # # ]: 0 : if(im4==tndim) then
443 : 0 : testdege=.false.
444 : 0 : current_dege=current_dege+1
445 : : endif
446 [ # # ]: 0 : if(testdege) then
447 : 0 : current_dege=current_dege+1
448 : : else
449 : : !new set of degenerate state: reorder it: put it into block diagonal
450 : : !form for column
451 : : if(prtopt==1) write(std_out,*) "newstarting, current_dege",newstarting, current_dege
452 : 0 : shift=0
453 [ # # ]: 0 : do im1=1,tndim ! balaye les premiers coefficients puis les autres
454 : :
455 : : ! if(im1==1) nonnuldege(im1)=0
456 : : ! if(im1>1) nonnuldege(im1)=nonnuldege(im1-1)
457 : 0 : maxvalue=0.00000001
458 : 0 : swap=.false.
459 [ # # ]: 0 : do im2=newstarting+shift,newstarting+current_dege-1
460 [ # # ]: 0 : if(abs(matrix(im1,im2))>maxvalue) then
461 : 0 : maxvalue=abs(matrix(im1,im2))
462 : 0 : indice_formax=im2
463 : 0 : swap=.true.
464 : : endif
465 : : enddo
466 : : ! found max value: permute
467 [ # # ]: 0 : if(swap) then
468 [ # # ]: 0 : do im3=1,tndim
469 : 0 : tmpx=matrix(im3,indice_formax)
470 : 0 : matrix(im3,indice_formax)=matrix(im3,newstarting+shift)
471 : 0 : matrix(im3,newstarting+shift)=tmpx
472 : : enddo
473 : 0 : shift=shift+1
474 : : endif
475 : : !write(std_out,*) "Eigenvectors after m1"
476 : : !do im3=1,tndim
477 : : ! write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im3,im5),im5=1,tndim)
478 : : !end do
479 : :
480 : : enddo
481 : : if(prtopt==1) then
482 : : write(std_out,*) "Eigenvectors after set of dege"
483 : : do im2=1,tndim
484 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im2,im3),im3=1,tndim)
485 : : end do
486 : : endif
487 : 0 : newstarting=im4
488 : 0 : current_dege=1
489 : : endif
490 : : enddo
491 [ # # ]: 0 : ABI_FREE(nonnuldege)
492 : : if(prtopt==1) then
493 : : write(std_out,*) "Ordered Eigenvectors"
494 : : do im1=1,tndim
495 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
496 : : end do
497 : : endif
498 : :
499 : : if(prtopt==1) then
500 : : write(std_out,*) "inverse operation: reconstitute original matrix: only the line here"
501 : : endif
502 [ # # # # : 0 : ABI_MALLOC(Apermutlineback,(tndim,tndim))
# # # # #
# # # ]
503 [ # # # # ]: 0 : Apermutlineback=zero
504 [ # # ]: 0 : do im1=1,tndim
505 [ # # ]: 0 : do im2=1,tndim
506 : 0 : Apermutlineback(im1,im2)=zero
507 [ # # ]: 0 : do im3=1,tndim
508 : 0 : Apermutlineback(im1,im2)=matrix(im3,im2)*Permutline(im3,im1)+Apermutlineback(im1,im2)
509 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
510 : : end do
511 : : end do
512 : : end do
513 [ # # # # ]: 0 : matrix=Apermutlineback
514 : : if(prtopt==1) then
515 : : write(std_out,*) "t(Permutline)*Apermutcol"
516 : : do im1=1,tndim
517 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
518 : : end do
519 : : endif
520 : :
521 : : ! Now, set the first coefficient of eigenvectors positive.
522 : :
523 [ # # ]: 0 : do im2=1,tndim ! loop over eigenvectors
524 [ # # ]: 0 : do im1=1,tndim ! loop over components
525 [ # # ]: 0 : if(abs(matrix(im1,im2))>tol8) then
526 [ # # ]: 0 : if(matrix(im1,im2)<0) then
527 [ # # ]: 0 : do im3=1,tndim
528 [ # # ]: 0 : if(abs(matrix(im3,im2))>tol8) then
529 : 0 : matrix(im3,im2)=-matrix(im3,im2)
530 : : endif
531 : : enddo
532 : : endif
533 : : exit
534 : : endif
535 : : enddo
536 : : enddo
537 : : if(prtopt==1) then
538 : : write(std_out,*) "Impose first component of eigenvectors is positive"
539 : : do im1=1,tndim
540 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
541 : : end do
542 : : endif
543 : : ! write(std_out,*) "inverse operation: reconstitute original matrix: then the column"
544 : : ! Apermutcolback=zero
545 : : ! do im1=1,tndim
546 : : ! do im2=1,tndim
547 : : ! Apermutcolback(im1,im2)=zero
548 : : ! do im3=1,tndim
549 : : ! Apermutcolback(im1,im2)=Apermutlineback(im1,im3)*Permutcol(im2,im3)+Apermutcolback(im1,im2)
550 : : ! ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
551 : : ! end do
552 : : ! end do
553 : : ! end do
554 : : ! write(std_out,*) "Apermutlineback*t(Permutcol)"
555 : : ! do im1=1,10
556 : : ! write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcolback(im1,im2),im2=1,10)
557 : : ! end do
558 : :
559 : :
560 [ # # ]: 0 : ABI_FREE(Apermutlineback)
561 [ # # ]: 0 : ABI_FREE(Apermutline)
562 [ # # ]: 0 : ABI_FREE(matrix_save)
563 [ # # ]: 0 : ABI_FREE(Apermutcol)
564 [ # # ]: 0 : ABI_FREE(work)
565 [ # # ]: 0 : ABI_FREE(Permutcol)
566 [ # # ]: 0 : ABI_FREE(nonnul)
567 [ # # ]: 0 : ABI_FREE(Permutline)
568 : :
569 : : DBG_EXIT("COLL")
570 : :
571 : 0 : end subroutine blockdiago_fordsyev
572 : : !!***
573 : :
574 : : !! FUNCTION
575 : : !! Transform matrix into block diagonal form before diagonalisation
576 : : !!
577 : : !! INPUTS
578 : : !! ndim = dimension of matrix
579 : : !! matrix= matrix
580 : : !!
581 : : !! OUTPUT
582 : : !! matrix= square root of the matrix
583 : : !!
584 : : !! SOURCE
585 : :
586 : 0 : subroutine blockdiago_forzheev(matrix,tndim,eig)
587 : :
588 : : !Arguments ------------------------------------
589 : : !scalars
590 : : integer,intent(in) :: tndim
591 : : complex(dp),intent(inout) :: matrix(tndim,tndim)
592 : : real(dp),intent(inout) :: eig(tndim)
593 : : !arrays
594 : :
595 : : !Local variables-------------------------------
596 : : !scalars
597 : : integer :: im1,im2,im3,info,lwork
598 : : character(len=500) :: message
599 : : complex(dp):: tmpx
600 : : integer(dp):: tmpi,prtopt
601 : : !arrays
602 : 0 : real(dp),allocatable :: rwork(:)
603 : 0 : complex(dp),allocatable :: work(:)
604 : 0 : real(dp),allocatable :: Permutcol(:,:)
605 : 0 : complex(dp),allocatable :: Apermutcol(:,:)
606 : 0 : complex(dp),allocatable :: Apermutline(:,:)
607 : 0 : complex(dp),allocatable :: Apermutlineback(:,:)
608 : 0 : real(dp),allocatable :: Permutline(:,:)
609 : 0 : complex(dp),allocatable :: matrix_save(:,:) !,W(:)
610 : 0 : integer,allocatable :: nonnul(:)
611 : : ! *************************************************************************
612 : :
613 : : !!!Do not remove this silly print instruction. Seems needed to avoid floating
614 : : !!!point exception on vm1_gcc51 ...
615 : : !!#if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
616 : : !! write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
617 : : !!#endif
618 : : DBG_ENTER("COLL")
619 : :
620 : : !printing option for debug
621 : 0 : prtopt = 0
622 : :
623 : 0 : lwork=10*tndim
624 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
625 [ # # # # : 0 : ABI_MALLOC(rwork,(3*tndim-2))
# # # # #
# # # #
# ]
626 : :
627 [ # # # # : 0 : ABI_MALLOC(matrix_save,(tndim,tndim))
# # # # #
# # # # #
# # ]
628 [ # # # # : 0 : matrix_save=matrix
# # ]
629 : :
630 [ # # # # : 0 : ABI_MALLOC(Permutcol,(tndim,tndim))
# # # # #
# # # # #
# # ]
631 : :
632 [ # # # # ]: 0 : Permutcol=zero
633 [ # # ]: 0 : do im1=1,tndim
634 : 0 : Permutcol(im1,im1)=1.d0
635 : : end do
636 : :
637 : : if(prtopt == 1) then
638 : : write(std_out,*) "MATRIX"
639 : : do im1=1,tndim
640 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix_save(im1,im2),im2=1,tndim)
641 : : end do
642 : : endif
643 : :
644 [ # # # # : 0 : ABI_MALLOC(nonnul,(tndim))
# # # # #
# # # ]
645 [ # # ]: 0 : do im1=1,tndim
646 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
647 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
648 [ # # ]: 0 : do im2=1,tndim
649 [ # # # # ]: 0 : if (abs(matrix(im1,im2))>0.000000000001.and.im2>nonnul(im1)) then
650 : 0 : nonnul(im1)=nonnul(im1)+1
651 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
652 : : ! permute
653 [ # # ]: 0 : do im3=1,tndim
654 : 0 : tmpx=matrix(im3,im2)
655 : 0 : matrix(im3,im2)=matrix(im3,nonnul(im1))
656 : 0 : matrix(im3,nonnul(im1))=tmpx
657 : 0 : tmpi=Permutcol(im3,im2)
658 : 0 : Permutcol(im3,im2)=Permutcol(im3,nonnul(im1))
659 : 0 : Permutcol(im3,nonnul(im1))=tmpi
660 : : enddo
661 [ # # ]: 0 : elseif (abs(matrix(im1,im2))<0.000000000001) then
662 : 0 : matrix(im1,im2)=czero
663 : : endif
664 : : enddo
665 : : enddo
666 : :
667 : : if (prtopt == 1) then
668 : : write(std_out,*) "MATRIX AFTER COLUMN PERMUT"
669 : : do im1=1,tndim
670 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
671 : : end do
672 : : write(std_out,*) "Permutcol MATRIX AFTER"
673 : : do im1=1,tndim
674 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Permutcol(im1,im2),im2=1,tndim)
675 : : end do
676 : : endif
677 : :
678 [ # # # # : 0 : ABI_MALLOC(Apermutcol,(tndim,tndim))
# # # # #
# # # ]
679 [ # # # # ]: 0 : Apermutcol=czero
680 [ # # ]: 0 : do im1=1,tndim
681 [ # # ]: 0 : do im2=1,tndim
682 : 0 : Apermutcol(im1,im2)=czero
683 [ # # ]: 0 : do im3=1,tndim
684 : 0 : Apermutcol(im1,im2)=matrix_save(im1,im3)*Permutcol(im3,im2)+Apermutcol(im1,im2)
685 : : ! write(std_out,*) im1,im2,im3,Apermutcol(im1,im2)
686 : : end do
687 : : end do
688 : : end do
689 : :
690 : : if(prtopt == 1) then
691 : : write(std_out,*) "Check product of original matrix by permutation matrix "
692 : : write(std_out,*) "Asave*Permutcol"
693 : : do im1=1,tndim
694 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Apermutcol(im1,im2),im2=1,tndim)
695 : : end do
696 : : endif
697 : :
698 : :
699 [ # # # # : 0 : ABI_MALLOC(Permutline,(tndim,tndim))
# # # # #
# # # ]
700 [ # # # # ]: 0 : Permutline=zero
701 [ # # ]: 0 : do im1=1,tndim
702 : 0 : Permutline(im1,im1)=1.d0
703 : : end do
704 : :
705 [ # # ]: 0 : do im1=1,tndim
706 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
707 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
708 [ # # ]: 0 : do im2=1,tndim
709 : : ! write(std_out,*) "im1,im2, abs matrix(im2,im1),nonnul(im1)",im1,im2,abs(B(im1,im2)),nonnul(im1)
710 [ # # # # ]: 0 : if (abs(matrix(im2,im1))>0.000000000001.and.im2>nonnul(im1)) then
711 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
712 : 0 : nonnul(im1)=nonnul(im1)+1
713 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
714 : : ! permute
715 [ # # ]: 0 : do im3=1,tndim
716 : 0 : tmpx=matrix(im2,im3)
717 : 0 : matrix(im2,im3)=matrix(nonnul(im1),im3)
718 : 0 : matrix(nonnul(im1),im3)=tmpx
719 : 0 : tmpi=Permutline(im2,im3)
720 : 0 : Permutline(im2,im3)=Permutline(nonnul(im1),im3)
721 : 0 : Permutline(nonnul(im1),im3)=tmpi
722 : : enddo
723 [ # # ]: 0 : elseif (abs(matrix(im2,im1))<0.000000000001) then
724 : 0 : matrix(im2,im1)=czero
725 : : endif
726 : : enddo
727 : : enddo
728 : :
729 : : if(prtopt == 1) then
730 : : write(std_out,*) "matrix AFTER"
731 : : do im1=1,tndim
732 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
733 : : end do
734 : : write(std_out,*) "Permutline MATRIX AFTER"
735 : : do im1=1,tndim
736 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Permutline(im1,im2),im2=1,tndim)
737 : : end do
738 : : endif
739 : :
740 [ # # # # : 0 : ABI_MALLOC(Apermutline,(tndim,tndim))
# # # # #
# # # ]
741 [ # # # # ]: 0 : Apermutline=czero
742 [ # # ]: 0 : do im1=1,tndim
743 [ # # ]: 0 : do im2=1,tndim
744 : 0 : Apermutline(im1,im2)=czero
745 [ # # ]: 0 : do im3=1,tndim
746 : 0 : Apermutline(im1,im2)=Apermutcol(im3,im2)*Permutline(im1,im3)+Apermutline(im1,im2)
747 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
748 : : end do
749 : : end do
750 : : end do
751 : :
752 : : if(prtopt == 1) then
753 : : write(std_out,*) "Check product of Apermutcol matrix by permutation matrix of the line "
754 : : write(std_out,*) "Permutline*Apermutcol"
755 : : do im1=1,tndim
756 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Apermutline(im1,im2),im2=1,tndim)
757 : : end do
758 : : endif
759 : :
760 [ # # ]: 0 : work=czero
761 : 0 : call zheev('v','u',tndim,matrix_save,tndim,eig,work,lwork,rwork,info)
762 [ # # ]: 0 : if(info/=0) then
763 : 0 : message = 'Error in diagonalization of matrix (zheev) ! - '
764 : 0 : ABI_ERROR(message)
765 : : end if
766 : :
767 : : if(prtopt == 1) then
768 : : write(std_out,*) 'output',INFO
769 : : write(std_out,*) "Eigenvalues"
770 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
771 : : write(std_out,*) "Eigenvectors"
772 : : do im1=1,tndim
773 : : write(std_out,'(2(1x,30(1x,f20.15,f20.15)))') (matrix_save(im1,im2),im2=1,tndim)
774 : : end do
775 : : endif
776 : :
777 : : ! call dsyev('v','u',tndim,A,LDA,W,WORKTMP,LWORK,INFO)
778 : : ! write(std_out,*) "optimal lwork",worktmp(1)
779 : :
780 [ # # ]: 0 : work=czero
781 : 0 : call zheev('v','u',tndim,matrix,tndim,eig,work,lwork,rwork,info)
782 [ # # ]: 0 : if(info/=0) then
783 : 0 : message = 'Error in diagonalization of matrix (zheev) ! - '
784 : 0 : ABI_ERROR(message)
785 : : end if
786 : :
787 : : if(prtopt == 1) then
788 : : write(std_out,*) 'output',INFO
789 : : write(std_out,*) "Eigenvalues"
790 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
791 : : write(std_out,*) "Eigenvectors"
792 : : do im1=1,tndim
793 : : write(std_out,'(2(1x,30(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
794 : : end do
795 : : endif
796 : :
797 [ # # # # : 0 : ABI_MALLOC(Apermutlineback,(tndim,tndim))
# # # # #
# # # ]
798 [ # # # # ]: 0 : Apermutlineback=czero
799 [ # # ]: 0 : do im1=1,tndim
800 [ # # ]: 0 : do im2=1,tndim
801 : 0 : Apermutlineback(im1,im2)=czero
802 [ # # ]: 0 : do im3=1,tndim
803 : 0 : Apermutlineback(im1,im2)=matrix(im3,im2)*Permutline(im3,im1)+Apermutlineback(im1,im2)
804 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
805 : : end do
806 : : end do
807 : : end do
808 [ # # # # ]: 0 : matrix=Apermutlineback
809 : :
810 : : if(prtopt == 1) then
811 : : write(std_out,*) "inverse operation: reconstitute original matrix: first the line"
812 : : write(std_out,*) "t(Permutline)*Apermutcol"
813 : : do im1=1,tndim
814 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
815 : : end do
816 : : endif
817 : :
818 : : ! write(std_out,*) "inverse operation: reconstitute original matrix: then the column"
819 : : ! Apermutcolback=zero
820 : : ! do im1=1,tndim
821 : : ! do im2=1,tndim
822 : : ! Apermutcolback(im1,im2)=zero
823 : : ! do im3=1,tndim
824 : : ! Apermutcolback(im1,im2)=Apermutlineback(im1,im3)*Permutcol(im2,im3)+Apermutcolback(im1,im2)
825 : : ! ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
826 : : ! end do
827 : : ! end do
828 : : ! end do
829 : : ! write(std_out,*) "Apermutlineback*t(Permutcol)"
830 : : ! do im1=1,10
831 : : ! write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcolback(im1,im2),im2=1,10)
832 : : ! end do
833 : :
834 : :
835 [ # # ]: 0 : ABI_FREE(Apermutlineback)
836 [ # # ]: 0 : ABI_FREE(Apermutline)
837 [ # # ]: 0 : ABI_FREE(matrix_save)
838 [ # # ]: 0 : ABI_FREE(Apermutcol)
839 [ # # ]: 0 : ABI_FREE(work)
840 [ # # ]: 0 : ABI_FREE(Permutcol)
841 [ # # ]: 0 : ABI_FREE(nonnul)
842 [ # # ]: 0 : ABI_FREE(Permutline)
843 [ # # ]: 0 : ABI_FREE(rwork)
844 : :
845 : : DBG_EXIT("COLL")
846 : :
847 : 0 : end subroutine blockdiago_forzheev
848 : : !!***
849 : :
850 : : !! FUNCTION
851 : : !! Compute the determinant of a 3x3 real matrix
852 : : !!
853 : : !! INPUTS
854 : : !! A = 3x3 matrix
855 : : !!
856 : : !! OUTPUT
857 : : !! det = The determinant
858 : : !!
859 : : !! SOURCE
860 : :
861 : 92 : function real_mat33det(A) result(det)
862 : : real(dp), intent(in) :: A(3,3)
863 : : real(dp) :: det
864 : : DET = A(1,1)*A(2,2)*A(3,3) &
865 : : - A(1,1)*A(2,3)*A(3,2) &
866 : : - A(1,2)*A(2,1)*A(3,3) &
867 : : + A(1,2)*A(2,3)*A(3,1) &
868 : : + A(1,3)*A(2,1)*A(3,2) &
869 : 92 : - A(1,3)*A(2,2)*A(3,1)
870 : 92 : end function real_mat33det
871 : : !!***
872 : :
873 : : !! FUNCTION
874 : : !! Compute the determinant of a 3x3 integer matrix
875 : : !!
876 : : !! INPUTS
877 : : !! A = 3x3 matrix
878 : : !!
879 : : !! OUTPUT
880 : : !! det = The determinant
881 : : !!
882 : : !! SOURCE
883 : :
884 : 6 : function int_mat33det(A) result(det)
885 : : integer, intent(in) :: A(3,3)
886 : : integer :: det
887 : : DET = A(1,1)*A(2,2)*A(3,3) &
888 : : - A(1,1)*A(2,3)*A(3,2) &
889 : : - A(1,2)*A(2,1)*A(3,3) &
890 : : + A(1,2)*A(2,3)*A(3,1) &
891 : : + A(1,3)*A(2,1)*A(3,2) &
892 : 6 : - A(1,3)*A(2,2)*A(3,1)
893 : 6 : end function int_mat33det
894 : : !!***
895 : :
896 : : !!****f* m_matrix/mati3inv
897 : : !! NAME
898 : : !! mati3inv
899 : : !!
900 : : !! FUNCTION
901 : : !! Invert and transpose orthogonal 3x3 matrix of INTEGER elements.
902 : : !!
903 : : !! INPUTS
904 : : !! mm = integer matrix to be inverted
905 : : !!
906 : : !! OUTPUT
907 : : !! mit = inverse of mm input matrix
908 : : !!
909 : : !! NOTES
910 : : !! Used for symmetry operations.
911 : : !! This routine applies to ORTHOGONAL matrices only.
912 : : !! Since these form a group, inverses are also integer arrays.
913 : : !! Returned array is TRANSPOSE of inverse, as needed.
914 : : !! Note use of integer arithmetic.
915 : : !!
916 : : !! SOURCE
917 : :
918 : 8454861 : subroutine mati3inv(mm, mit)
919 : :
920 : : !Arguments ------------------------------------
921 : : !arrays
922 : : integer,intent(in) :: mm(3,3)
923 : : integer,intent(out) :: mit(3,3)
924 : :
925 : : !Local variables-------------------------------
926 : : !scalars
927 : : integer :: dd
928 : : character(len=500) :: msg
929 : : !arrays
930 : : integer :: tt(3,3)
931 : : ! *************************************************************************
932 : :
933 : 8454861 : tt(1,1) = mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)
934 : 8454861 : tt(2,1) = mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)
935 : 8454861 : tt(3,1) = mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3)
936 : 8454861 : tt(1,2) = mm(3,1) * mm(2,3) - mm(2,1) * mm(3,3)
937 : 8454861 : tt(2,2) = mm(1,1) * mm(3,3) - mm(3,1) * mm(1,3)
938 : 8454861 : tt(3,2) = mm(2,1) * mm(1,3) - mm(1,1) * mm(2,3)
939 : 8454861 : tt(1,3) = mm(2,1) * mm(3,2) - mm(3,1) * mm(2,2)
940 : 8454861 : tt(2,3) = mm(3,1) * mm(1,2) - mm(1,1) * mm(3,2)
941 : 8454861 : tt(3,3) = mm(1,1) * mm(2,2) - mm(2,1) * mm(1,2)
942 : 8454861 : dd = mm(1,1) * tt(1,1) + mm(2,1) * tt(2,1) + mm(3,1) * tt(3,1)
943 : :
944 : : ! Make sure matrix is not singular
945 [ + - ]: 8454861 : if (dd /= 0) then
946 [ + + + + ]: 109913193 : mit(:,:)=tt(:,:)/dd
947 : : else
948 : 0 : write(msg, '(2a,2x,9(i0,1x),a)' )'Attempting to invert integer array',ch10,mm,' ==> determinant is zero.'
949 : 0 : ABI_ERROR(msg)
950 : : end if
951 : :
952 : : ! If matrix is orthogonal, determinant must be 1 or -1
953 [ - + ]: 8454861 : if (abs(dd) /= 1) then
954 : 0 : write(msg, '(3a,i0)' )'Absolute value of determinant should be one',ch10,'but determinant= ',dd
955 : 0 : ABI_ERROR(msg)
956 : : end if
957 : :
958 : 8454861 : end subroutine mati3inv
959 : : !!***
960 : :
961 : : !!****f* m_matrix/mati3det
962 : : !! NAME
963 : : !! mati3det
964 : : !!
965 : : !! FUNCTION
966 : : !! Compute the determinant of a 3x3 matrix of INTEGER elements.
967 : : !!
968 : : !! INPUTS
969 : : !! mm = integer matrix
970 : : !!
971 : : !! OUTPUT
972 : : !! det = determinant of the matrix
973 : : !!
974 : : !! SOURCE
975 : :
976 : 2770649 : subroutine mati3det(mm, det)
977 : :
978 : : !Arguments ------------------------------------
979 : : integer,intent(in) :: mm(3,3)
980 : : integer,intent(out) :: det
981 : : ! *************************************************************************
982 : :
983 : : det=mm(1,1)*(mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)) &
984 : : + mm(2,1)*(mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)) &
985 : 2770649 : + mm(3,1)*(mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3))
986 : :
987 : 2770649 : end subroutine mati3det
988 : : !!***
989 : :
990 : : !!****f* m_matrix/matr3inv
991 : : !! NAME
992 : : !! matr3inv
993 : : !!
994 : : !! FUNCTION
995 : : !! Invert and transpose general 3x3 matrix of real*8 elements.
996 : : !!
997 : : !! INPUTS
998 : : !! aa = 3x3 matrix to be inverted
999 : : !!
1000 : : !! OUTPUT
1001 : : !! ait = inverse of aa input matrix
1002 : : !!
1003 : : !! NOTES
1004 : : !! Returned array is TRANSPOSE of inverse, as needed to get g from r.
1005 : : !!
1006 : : !! SOURCE
1007 : :
1008 : 4068277 : subroutine matr3inv(aa, ait)
1009 : :
1010 : : !Arguments ------------------------------------
1011 : : !arrays
1012 : : real(dp),intent(in) :: aa(3,3)
1013 : : real(dp),intent(out) :: ait(3,3)
1014 : :
1015 : : !Local variables-------------------------------
1016 : : !scalars
1017 : : real(dp) :: dd,det,t1,t2,t3
1018 : : character(len=500) :: msg
1019 : : ! *************************************************************************
1020 : :
1021 : 4068277 : t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
1022 : 4068277 : t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
1023 : 4068277 : t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
1024 : 4068277 : det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
1025 : :
1026 : : !Make sure matrix is not singular
1027 [ + - ]: 4068277 : if (abs(det)>tol16) then
1028 : 4068277 : dd=one/det
1029 : : else
1030 : : write(msg, '(2a,2x,9es16.8,a,a,es16.8,a)' )&
1031 : 0 : 'Attempting to invert real(8) 3x3 array',ch10,aa(:,:),ch10,' ==> determinant=',det,' is zero.'
1032 : 0 : ABI_BUG(msg)
1033 : : end if
1034 : :
1035 : 4068277 : ait(1,1) = t1 * dd
1036 : 4068277 : ait(2,1) = t2 * dd
1037 : 4068277 : ait(3,1) = t3 * dd
1038 : 4068277 : ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
1039 : 4068277 : ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
1040 : 4068277 : ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
1041 : 4068277 : ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
1042 : 4068277 : ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
1043 : 4068277 : ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
1044 : :
1045 : 4068277 : end subroutine matr3inv
1046 : : !!***
1047 : :
1048 : : !!****f* m_matrix/mati3_add_columns
1049 : : !! NAME
1050 : : !! mati3_add_columns
1051 : : !!
1052 : : !! FUNCTION
1053 : : !! Replace columns MM[:, ii] by aa*MM[:, ii] + bb*MM[:, jj]
1054 : : !! and MM[:, jj] by cc*MM[:, ii] + dd*M[:, jj]
1055 : : !!
1056 : : !! INPUTS
1057 : : !!
1058 : : !! OUTPUT
1059 : : !!
1060 : : !! NOTES
1061 : : !!
1062 : : !! SOURCE
1063 : :
1064 : : subroutine mati3_add_columns(MM, ii, jj, aa, bb, cc, dd)
1065 : :
1066 : : !Arguments ------------------------------------
1067 : : !arrays
1068 : : integer, intent(inout) :: MM(3,3)
1069 : : integer, intent(in) :: ii, jj, aa, bb, cc, dd
1070 : :
1071 : : !Local variables-------------------------------
1072 : : !scalars
1073 : : integer :: kk, Ui, Uj
1074 : :
1075 : : ! *************************************************************************
1076 : : do kk=1,3
1077 : : Ui = MM(kk,ii)
1078 : : Uj = MM(kk,jj)
1079 : : MM(kk,ii) = aa*Ui + bb*Uj
1080 : : MM(kk,jj) = cc*Ui + dd*Uj
1081 : : end do
1082 : :
1083 : : end subroutine mati3_add_columns
1084 : : !!***
1085 : :
1086 : : !!****f* m_matrix/mati3_add_rows
1087 : : !! NAME
1088 : : !! mati3_add_rows
1089 : : !!
1090 : : !! FUNCTION
1091 : : !! Replace rows
1092 : : !! MM[ii,:] by aa*MM[ii,:] + bb*MM[jj,:]
1093 : : !! and
1094 : : !! MM[jj,:] by cc*MM[ii,:] + dd*M[jj,:]
1095 : : !!
1096 : : !! INPUTS
1097 : : !!
1098 : : !! OUTPUT
1099 : : !!
1100 : : !! NOTES
1101 : : !!
1102 : : !! SOURCE
1103 : :
1104 : : subroutine mati3_add_rows(MM, ii, jj, aa, bb, cc, dd)
1105 : :
1106 : : !Arguments ------------------------------------
1107 : : !arrays
1108 : : integer, intent(inout) :: MM(3,3)
1109 : : integer, intent(in) :: ii, jj, aa, bb, cc, dd
1110 : :
1111 : : !Local variables-------------------------------
1112 : : !scalars
1113 : : integer :: kk, Ui, Uj
1114 : :
1115 : : ! *************************************************************************
1116 : : do kk=1,3
1117 : : Ui = MM(ii,kk)
1118 : : Uj = MM(jj,kk)
1119 : : MM(ii,kk) = aa*Ui + bb*Uj
1120 : : MM(jj,kk) = cc*Ui + dd*Uj
1121 : : end do
1122 : :
1123 : : end subroutine mati3_add_rows
1124 : : !!***
1125 : :
1126 : :
1127 : : !!****f* m_matrix/xgcd
1128 : : !! NAME
1129 : : !! xgcd
1130 : : !!
1131 : : !! FUNCTION
1132 : : !! Extended Euclidian algorithm for greater common divisor (gcd).
1133 : : !! g = gcd(a,b)
1134 : : !! a*x + b*y = g
1135 : : !!
1136 : : !! INPUTS
1137 : : !!
1138 : : !! OUTPUT
1139 : : !!
1140 : : !! NOTES
1141 : : !!
1142 : : !! SOURCE
1143 : :
1144 : : subroutine xgcd(aa,bb,xx,yy,gg)
1145 : :
1146 : : !Arguments ------------------------------------
1147 : : !arrays
1148 : : integer, intent(in) :: aa,bb
1149 : : integer, intent(out) :: xx,yy,gg
1150 : :
1151 : : !Local variables-------------------------------
1152 : : !scalars
1153 : : integer :: cc, qq
1154 : : integer :: old_xx, old_yy, old_gg
1155 : :
1156 : : ! *************************************************************************
1157 : :
1158 : : old_gg = aa; gg = bb
1159 : : old_xx = 1; xx = 0
1160 : : old_yy = 0; yy = 1
1161 : :
1162 : : do while (gg.ne.0)
1163 : : qq = old_gg / gg
1164 : : cc = gg; gg = old_gg - qq*cc; old_gg = cc
1165 : : cc = xx; xx = old_xx - qq*cc; old_xx = cc
1166 : : cc = yy; yy = old_yy - qq*cc; old_yy = cc
1167 : : end do
1168 : :
1169 : : gg = old_gg
1170 : : xx = old_xx
1171 : : yy = old_yy
1172 : :
1173 : : end subroutine xgcd
1174 : : !!***
1175 : :
1176 : : !!-------------------------------------------------------------
1177 : :
1178 : : !! FUNCTION
1179 : : !! Swap rows i and j of matrix A.
1180 : : !!
1181 : : !! INPUTS
1182 : : !!
1183 : : !! OUTPUT
1184 : : !!
1185 : : !! SOURCE
1186 : :
1187 : 34 : subroutine swap_rows(AA, ii, jj)
1188 : : integer, intent(inout) :: AA(3,3)
1189 : : integer, intent(in) :: ii, jj
1190 : : integer :: row(3)
1191 [ + + ]: 136 : row(:) = AA(ii,:)
1192 [ + + ]: 136 : AA(ii,:) = AA(jj,:)
1193 [ + + ]: 136 : AA(jj,:) = row
1194 : 34 : end subroutine swap_rows
1195 : :
1196 : : !!****f* m_matrix/mati3HNF
1197 : : !! NAME
1198 : : !! mati3HNF
1199 : : !!
1200 : : !! FUNCTION
1201 : : !! Compute the row-style Hermite normal form (HNF) of a 3x3 integer matrix A to express
1202 : : !! H = U * A
1203 : : !! where U is unimodular (det = +1 or -1),
1204 : : !! and H is upper triangular as
1205 : : !! H = H11 H12 H13
1206 : : !! 0 H22 H23
1207 : : !! 0 0 H33
1208 : : !! with
1209 : : !! Hii > 0
1210 : : !! 0 <= Hij < Hii for i < j
1211 : : !!
1212 : : !!
1213 : : !! INPUTS
1214 : : !! A = The input matrix of integers to be decomposed.
1215 : : !!
1216 : : !! OUTPUT
1217 : : !! H = The HNF upper triangular matrix.
1218 : : !!
1219 : : !! NOTES
1220 : : !!
1221 : : !! The matrix U is discarded.
1222 : : !! This subroutine is used to enumerate the primitive lattice vectors
1223 : : !! inside a supercell defined by the supercell matrix SS.
1224 : : !!
1225 : : !! SOURCE
1226 : :
1227 : 46 : subroutine mati3HNF(A,H)
1228 : : integer, intent(in) :: A(3,3)
1229 : : integer, intent(out) :: H(3,3)
1230 : :
1231 : : logical :: done
1232 : : integer :: ii, kk
1233 : : integer :: pivot_row
1234 : : integer :: qq
1235 : :
1236 : 46 : H = A
1237 : :
1238 [ + + ]: 184 : do ii=1,3
1239 : :
1240 : : !------------------------------------------
1241 : : ! Find a nonzero pivot in column ii
1242 : : !------------------------------------------
1243 : 138 : pivot_row = 0
1244 [ + - ]: 155 : do kk = ii, 3
1245 [ + + ]: 155 : if (H(kk,ii) /= 0) then
1246 : 138 : pivot_row = kk
1247 : 138 : exit
1248 : : end if
1249 : : end do
1250 : :
1251 : : ! The remaining column is already zero
1252 [ + - ]: 138 : if (pivot_row == 0) cycle
1253 : :
1254 : : !------------------------------------------
1255 : : ! Swap pivot row into place
1256 : : !------------------------------------------
1257 [ + + ]: 138 : if (pivot_row /= ii) then
1258 : 17 : call swap_rows(H,ii,pivot_row)
1259 : : end if
1260 : :
1261 : : !------------------------------------------
1262 : : ! Euclidean elimination below pivot
1263 : : !------------------------------------------
1264 : : do
1265 : 189 : done = .true.
1266 : :
1267 [ + + ]: 422 : do kk=ii+1,3
1268 : :
1269 [ + + ]: 422 : if (H(kk,ii) /= 0) then
1270 : :
1271 : 64 : qq = H(kk,ii) / H(ii,ii)
1272 : :
1273 [ + + ]: 256 : H(kk,:) = H(kk,:) - qq * H(ii,:)
1274 : :
1275 : : ! If remainder got smaller pivot candidate,
1276 : : ! swap rows and continue Euclidean reduction
1277 [ + - + + ]: 64 : if (abs(H(kk,ii)) < abs(H(ii,ii)) .and. H(kk,ii) /= 0) then
1278 : 17 : call swap_rows(H,ii,kk)
1279 : : end if
1280 : :
1281 : : done = .false.
1282 : : end if
1283 : :
1284 : : end do
1285 : :
1286 [ + + ]: 189 : if (done) exit
1287 : : end do
1288 : :
1289 : : !------------------------------------------
1290 : : ! Make pivot positive
1291 : : !------------------------------------------
1292 [ + + ]: 138 : if (H(ii,ii) < 0) then
1293 [ + + ]: 112 : H(ii,:) = -H(ii,:)
1294 : : end if
1295 : :
1296 : : !------------------------------------------
1297 : : ! Reduce entries above/below pivot
1298 : : !------------------------------------------
1299 [ + + ]: 598 : do kk = 1, 3
1300 : :
1301 [ + + + + ]: 552 : if (kk /= ii .and. H(kk,ii) /= 0) then
1302 : :
1303 : 47 : qq = euclid_div(H(kk,ii), H(ii,ii))
1304 : :
1305 [ + + ]: 188 : H(kk,:) = H(kk,:) - qq * H(ii,:)
1306 : :
1307 : : ! Ensure canonical remainder
1308 [ - + ]: 47 : do while (H(kk,ii) < 0)
1309 [ - - ]: 47 : H(kk,:) = H(kk,:) + H(ii,:)
1310 : : end do
1311 : :
1312 [ - + ]: 47 : do while (H(kk,ii) >= H(ii,ii))
1313 [ - - ]: 47 : H(kk,:) = H(kk,:) - H(ii,:)
1314 : : end do
1315 : :
1316 : : end if
1317 : :
1318 : : end do
1319 : :
1320 : : end do
1321 : :
1322 : 46 : end subroutine mati3HNF
1323 : : !!***
1324 : :
1325 : : !!****f* m_matrix/is_unitary
1326 : : !! NAME
1327 : : !! is_unitary
1328 : : !!
1329 : : !! FUNCTION
1330 : : !! Returns .TRUE. if the input complex matrix mat is unitary, i.e. mat^\dagger mat = I.
1331 : : !! Also returns the maximum absolute deviation from the identity matrix.
1332 : : !!
1333 : : !! INPUTS
1334 : : !! nn = Dimension of the matrix.
1335 : : !! mat = Complex matrix of size (nn,nn).
1336 : : !! tol = Tolerance for the maximum deviation.
1337 : : !!
1338 : : !! OUTPUT
1339 : : !! .TRUE. if err < tol, .FALSE. otherwise.
1340 : : !! err = Maximum absolute deviation from the identity matrix.
1341 : : !!
1342 : : !! SOURCE
1343 : :
1344 : 0 : logical function is_unitary(nn, mat, tol, err)
1345 : :
1346 : : !Arguments ------------------------------------
1347 : : integer, intent(in) :: nn
1348 : : complex(dp), intent(in) :: mat(nn,nn)
1349 : : real(dp), intent(in) :: tol
1350 : : real(dp), intent(out) :: err
1351 : :
1352 : : !Local variables-------------------------------
1353 : : integer :: ii
1354 : 0 : complex(dp) :: prod(nn,nn), identity(nn,nn)
1355 : : !----------------------------------------------------------------------
1356 : :
1357 : : ! Compute mat^\dagger mat
1358 [ # # # # : 0 : prod = matmul(conjg(transpose(mat)), mat)
# # # # #
# # # # #
# # # # ]
1359 : :
1360 : : ! Build identity matrix
1361 [ # # # # ]: 0 : identity = czero
1362 [ # # ]: 0 : do ii=1,nn
1363 : 0 : identity(ii,ii) = one
1364 : : end do
1365 : :
1366 : : ! Maximum deviation from identity
1367 [ # # # # : 0 : err = maxval(abs(prod - identity))
# # # # #
# # # ]
1368 : 0 : is_unitary = (err < tol)
1369 : :
1370 : 0 : end function is_unitary
1371 : : !!***
1372 : :
1373 : : !!****f* m_matrix/is_identity
1374 : : !! NAME
1375 : : !! is_identity
1376 : : !!
1377 : : !! FUNCTION
1378 : : !! Returns .TRUE. if the complex input matrix mat is the identity matrix.
1379 : : !! Also returns the maximum absolute deviation from the identity matrix.
1380 : : !!
1381 : : !! INPUTS
1382 : : !! nn = Dimension of the matrix.
1383 : : !! mat = Complex matrix of size (nn,nn).
1384 : : !! tol = Tolerance for the maximum deviation.
1385 : : !!
1386 : : !! OUTPUT
1387 : : !! .TRUE. if err < tol, .FALSE. otherwise.
1388 : : !! err = Maximum absolute deviation from the identity matrix.
1389 : : !!
1390 : : !! SOURCE
1391 : :
1392 : 0 : logical function is_identity(nn, mat, tol, err)
1393 : :
1394 : : !Arguments ------------------------------------
1395 : : integer, intent(in) :: nn
1396 : : complex(dp), intent(in) :: mat(nn,nn)
1397 : : real(dp), intent(in) :: tol
1398 : : real(dp), intent(out) :: err
1399 : :
1400 : : !Local variables-------------------------------
1401 : : integer :: ii
1402 : 0 : complex(dp) :: identity(nn,nn)
1403 : : !----------------------------------------------------------------------
1404 : :
1405 : : ! Build identity matrix
1406 [ # # # # ]: 0 : identity = czero
1407 [ # # ]: 0 : do ii=1,nn
1408 : 0 : identity(ii,ii) = one
1409 : : end do
1410 : :
1411 : : ! Maximum deviation from identity
1412 [ # # # # : 0 : err = maxval(abs(mat - identity))
# # # # #
# # # ]
1413 : 0 : is_identity = (err < tol)
1414 : :
1415 : 0 : end function is_identity
1416 : : !!***
1417 : :
1418 : : end module m_matrix
|
