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_hide_lapack, only : xginv
29 : :
30 : : implicit none
31 : :
32 : : private
33 : :
34 : : public :: invsqrt_matrix ! inv of Sqrt of Matrix
35 : : public :: blockdiago_fordsyev ! inv of Sqrt of Matrix
36 : : public :: blockdiago_forzheev ! inv of Sqrt of Matrix
37 : : public :: mat33det ! Determinant of a 3x3 matrix
38 : : public :: mati3inv ! Invert and transpose orthogonal 3x3 matrix of INTEGER elements.
39 : : public :: mati3det ! Compute the determinant of a 3x3 matrix of INTEGER elements.
40 : : public :: matr3inv ! Invert and TRANSPOSE general 3x3 matrix of real*8 elements.
41 : : public :: is_unitary ! Returns .TRUE. if the input matrix mat is unitary, i.e. mat^\dagger mat = I.
42 : : public :: is_identity ! Returns .TRUE. if the input matrix mat is the identity matrix.
43 : :
44 : : ! the determinant of a 3*3 matrix
45 : : interface mat33det
46 : : procedure real_mat33det
47 : : procedure int_mat33det
48 : : end interface mat33det
49 : :
50 : :
51 : : CONTAINS !===========================================================
52 : :
53 : : !! FUNCTION
54 : : !! Initialize matrix
55 : : !!
56 : : !! INPUTS
57 : : !! ndim = dimension of matrix
58 : : !! matrix= matrix
59 : : !!
60 : : !! OUTPUT
61 : : !! matrix= square root of the matrix
62 : : !! force_diag = 0 if it no 0 on diagonal
63 : : !! = nb of zeros found otherwise
64 : : !!
65 : : !! SOURCE
66 : :
67 : 298 : subroutine invsqrt_matrix(matrix,tndim,force_diag)
68 : :
69 : : !Arguments ------------------------------------
70 : : !scalars
71 : : integer,intent(in) :: tndim
72 : : complex(dp),intent(inout) :: matrix(tndim,tndim)
73 : : integer, intent(out) :: force_diag
74 : :
75 : : !Local variables-------------------------------
76 : : !scalars
77 : : integer :: im,im1,im2,info,lwork,nb_of_zero
78 : : character(len=500) :: message
79 : : real(dp) :: pawprtvol
80 : : !arrays
81 : 298 : real(dp),allocatable :: eig(:),rwork(:)
82 : 298 : complex(dp),allocatable :: zwork(:),diag(:,:)
83 : 298 : complex(dp),allocatable :: sqrtmat(:,:),zhdp2(:,:),sqrtmatinv(:,:)
84 : 298 : complex(dp),allocatable :: initialmatrix(:,:)
85 : : ! *************************************************************************
86 : :
87 : : !Do not remove this silly print instruction. Seems needed to avoid floating
88 : : !point exception on vm1_gcc51 ...
89 : : #if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
90 : : write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
91 : : #endif
92 : :
93 : : DBG_ENTER("COLL")
94 : 298 : pawprtvol=2
95 : :
96 [ - + + - : 1192 : ABI_MALLOC(initialmatrix,(tndim,tndim))
- + + - +
- + - - +
- + ]
97 [ - + + + : 9212 : initialmatrix=matrix
+ + ]
98 : : ! == First diagonalize matrix and keep the matrix for the change of basis
99 : 298 : lwork=2*tndim-1
100 [ + - - + : 894 : ABI_MALLOC(rwork,(3*tndim-2))
+ - + - +
- - + -
+ ]
101 [ - + + - : 894 : ABI_MALLOC(zwork,(lwork))
+ - + - -
+ - + ]
102 [ - + + - : 894 : ABI_MALLOC(eig,(tndim))
+ - + - -
+ - + ]
103 : :
104 : 298 : call zheev('v','u',tndim,matrix,tndim,eig,zwork,lwork,rwork,info)
105 : : if(pawprtvol>3) then
106 : : write(message,'(2a)') ch10,' - rotation matrix - '
107 : : call wrtout(std_out,message,'COLL')
108 : : do im1=1,tndim
109 : : write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
110 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
111 : : & (matrix(im1,im2),im2=1,tndim)
112 : : call wrtout(std_out,message,'COLL')
113 : : end do
114 : : endif
115 : :
116 : :
117 [ - + ]: 298 : ABI_FREE(zwork)
118 [ - + ]: 298 : ABI_FREE(rwork)
119 [ - + ]: 298 : if(info/=0) then
120 : 0 : message = 'Error in diagonalization of zmat (zheev) ! - '
121 : 0 : ABI_ERROR(message)
122 : : end if
123 : :
124 : : ! == Secondly Compute sqrt(diagonalized matrix)
125 [ - + + - : 894 : ABI_MALLOC(diag,(tndim,tndim))
- + + - -
+ - + ]
126 [ + + + + ]: 8914 : diag=czero
127 : : nb_of_zero=0
128 [ + + ]: 1725 : do im=1,tndim
129 : :
130 [ - + ]: 1725 : if(eig(im)< -tol8) then
131 : 0 : message = " - Eigenvalues from zheev are negative or zero ! - "
132 : 0 : write(std_out,*)
133 : 0 : write(std_out,*) " Eigenvalue=",eig(im)
134 : 0 : write(std_out,*) " Matrix is"
135 [ # # ]: 0 : do im1=1,tndim
136 : 0 : write(std_out,'(100f7.3)') (initialmatrix(im1,im2),im2=1,tndim)
137 : : enddo
138 : 0 : ABI_ERROR(message)
139 [ - + ]: 1427 : else if(abs(eig(im))<tol8) then
140 : 0 : nb_of_zero=nb_of_zero+1
141 : : else
142 : 1427 : diag(im,im)=cmplx(one/sqrt(eig(im)),zero,kind=dp)
143 : : endif
144 : : enddo
145 : 298 : force_diag=nb_of_zero
146 [ - + ]: 298 : ABI_FREE(eig)
147 : : ! write(std_out,*) "sqrt(eig) , diag(1,1)",sqrt(eig(1)),diag(1,1)
148 : : ! write(std_out,*) "cmplx(sqrt(eig(1)),zero,dp) , diag(1,1)",cmplx(sqrt(eig(1)),zero,dp),diag(1,1)
149 : : ! write(std_out,*) "sqrt(cmplx(eig(1),zero,dp)) , diag(1,1)",sqrt(cmplx(eig(1),zero,dp)),diag(1,1)
150 : :
151 : : ! == Thirdly Multiply by matrix for the change of basis
152 [ - + + - : 894 : ABI_MALLOC(sqrtmat,(tndim,tndim))
- + + - -
+ - + ]
153 [ - + + - : 894 : ABI_MALLOC(zhdp2,(tndim,tndim))
- + + - -
+ - + ]
154 : : if(pawprtvol>3) then
155 : : write(message,'(2a)') ch10,' - 1.0/sqrt(Eigenmatrix) - '
156 : : call wrtout(std_out,message,'COLL')
157 : : do im1=1,tndim
158 : : write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
159 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
160 : : & (diag(im1,im2),im2=1,tndim)
161 : : call wrtout(std_out,message,'COLL')
162 : : end do
163 : : endif
164 : : !zgemm(A,B,C) : C = op(A) op(B)
165 [ + - + + : 8914 : call zgemm('n','t',tndim,tndim,tndim,cone,diag,tndim,conjg(matrix),tndim,czero,zhdp2,tndim)
+ + - + ]
166 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,matrix,tndim,zhdp2,tndim,czero,sqrtmat,tndim)
167 : : ! if(abs(pawprtvol)>=3) then
168 : : if(pawprtvol>3) then
169 : : write(message,'(3a)') ch10," - inverse Sqrt root of matrix is - "
170 : : call wrtout(std_out,message,'COLL')
171 : : do im1=1,tndim
172 : : write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
173 : : & (sqrtmat(im1,im2),im2=1,tndim)
174 : : call wrtout(std_out,message,'COLL')
175 : : end do
176 : : endif
177 : : ! endif
178 [ - + ]: 298 : ABI_FREE(diag)
179 : :
180 : : ! == Forthly Compute the inverse of the square root
181 : : ! call matcginv_dpc(sqrtmat,tndim,tndim)
182 : : !call xginv(sqrtmat,tndim)
183 [ - + + - : 894 : ABI_MALLOC(sqrtmatinv,(tndim,tndim))
- + + - -
+ - + ]
184 [ - + + + : 9212 : sqrtmatinv=sqrtmat
+ + ]
185 : : if(pawprtvol>3) then
186 : : write(message,'(2a)') ch10," - inverse Sqrt root of matrix is - "
187 : : call wrtout(std_out,message,'COLL')
188 : : do im1=1,tndim
189 : : write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
190 : : & (sqrtmatinv(im1,im2),im2=1,tndim)
191 : : call wrtout(std_out,message,'COLL')
192 : : end do
193 : : endif
194 [ - + ]: 298 : ABI_FREE(sqrtmat)
195 : :
196 : : ! == Fifthly Check that O^{-0/5} O O{-0/5}=I
197 : : ! zgemm(A,B,C) : C = op(A) op(B)
198 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,initialmatrix,tndim,sqrtmatinv,tndim,czero,zhdp2,tndim)
199 : 298 : call zgemm('n','n',tndim,tndim,tndim,cone,sqrtmatinv,tndim,zhdp2,tndim,czero,initialmatrix,tndim)
200 : : if(pawprtvol>3) then
201 : : write(message,'(3a)') ch10," - O^{-0/5} O O^{-0/5}=I - "
202 : : call wrtout(std_out,message,'COLL')
203 : : do im1=1,tndim
204 : : write(message,'(12(1x,18(1x,"(",f10.6,",",f4.1,")")))')&
205 : : ! write(message,'(12(1x,18(1x,"(",f20.16,",",f20.16,")")))')&
206 : : & (initialmatrix(im1,im2),im2=1,tndim)
207 : : call wrtout(std_out,message,'COLL')
208 : : end do
209 : : endif
210 [ - + ]: 298 : ABI_FREE(zhdp2)
211 [ + + + + ]: 8914 : matrix=sqrtmatinv
212 [ - + ]: 298 : ABI_FREE(sqrtmatinv)
213 [ - + ]: 298 : ABI_FREE(initialmatrix)
214 : :
215 : : DBG_EXIT("COLL")
216 : :
217 : 298 : end subroutine invsqrt_matrix
218 : : !!***
219 : :
220 : : !! FUNCTION
221 : : !! Transform matrix into block diagonal form before diagonalisation
222 : : !!
223 : : !! INPUTS
224 : : !! ndim = dimension of matrix
225 : : !! matrix= matrix
226 : : !!
227 : : !! OUTPUT
228 : : !! matrix= square root of the matrix
229 : : !!
230 : : !! SOURCE
231 : :
232 : 0 : subroutine blockdiago_fordsyev(matrix,tndim,eig)
233 : :
234 : : !Arguments ------------------------------------
235 : : !scalars
236 : : integer,intent(in) :: tndim
237 : : real(dp),intent(inout) :: matrix(tndim,tndim)
238 : : real(dp),intent(inout) :: eig(tndim)
239 : : !arrays
240 : :
241 : : !Local variables-------------------------------
242 : : !scalars
243 : : integer :: im1,im2,im3,info,lwork,im4,indice_formax,shift !im5,
244 : : character(len=500) :: message
245 : : real(dp):: tmpx,maxvalue
246 : : integer(dp):: tmpi,newstarting,current_dege,prtopt
247 : : !arrays
248 : 0 : real(dp),allocatable :: work(:)
249 : 0 : real(dp),allocatable :: Permutcol(:,:)
250 : 0 : real(dp),allocatable :: Apermutcol(:,:)
251 : 0 : real(dp),allocatable :: Apermutline(:,:)
252 : 0 : real(dp),allocatable :: Apermutlineback(:,:)
253 : 0 : real(dp),allocatable :: Permutline(:,:)
254 : 0 : real(dp),allocatable :: matrix_save(:,:) !,W(:)
255 : 0 : integer,allocatable :: nonnul(:)
256 : 0 : integer,allocatable :: nonnuldege(:)
257 : : logical :: testdege,swap
258 : : ! *************************************************************************
259 : :
260 : : !!!Do not remove this silly print instruction. Seems needed to avoid floating
261 : : !!!point exception on vm1_gcc51 ...
262 : : !!#if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
263 : : !! write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
264 : : !!#endif
265 : : DBG_ENTER("COLL")
266 : :
267 : 0 : lwork=10*tndim
268 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
269 [ # # ]: 0 : work = zero
270 : :
271 [ # # # # : 0 : ABI_MALLOC(matrix_save,(tndim,tndim))
# # # # #
# # # # #
# # ]
272 [ # # # # : 0 : matrix_save=matrix
# # ]
273 : :
274 [ # # # # : 0 : ABI_MALLOC(Permutcol,(tndim,tndim))
# # # # #
# # # ]
275 : :
276 [ # # # # ]: 0 : Permutcol=zero
277 [ # # ]: 0 : do im1=1,tndim
278 : 0 : Permutcol(im1,im1)=1.d0
279 : : end do
280 : :
281 : 0 : prtopt=0
282 : :
283 [ # # # # : 0 : ABI_MALLOC(nonnul,(tndim))
# # # # #
# # # ]
284 [ # # ]: 0 : do im1=1,tndim
285 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
286 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
287 [ # # ]: 0 : do im2=1,tndim
288 [ # # # # ]: 0 : if (abs(matrix(im1,im2))>0.000000000001.and.im2>nonnul(im1)) then
289 : 0 : nonnul(im1)=nonnul(im1)+1
290 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
291 : : ! permute
292 [ # # ]: 0 : do im3=1,tndim
293 : 0 : tmpx=matrix(im3,im2)
294 : 0 : matrix(im3,im2)=matrix(im3,nonnul(im1))
295 : 0 : matrix(im3,nonnul(im1))=tmpx
296 : 0 : tmpi=Permutcol(im3,im2)
297 : 0 : Permutcol(im3,im2)=Permutcol(im3,nonnul(im1))
298 : 0 : Permutcol(im3,nonnul(im1))=tmpi
299 : : enddo
300 [ # # ]: 0 : elseif (abs(matrix(im1,im2))<0.000000000001) then
301 : 0 : matrix(im1,im2)=zero
302 : : endif
303 : : enddo
304 : : enddo
305 : : if(prtopt==1) then
306 : : write(std_out,*) "MATRIX AFTER COLUMN PERMUT"
307 : : do im1=1,tndim
308 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
309 : : end do
310 : : write(std_out,*) "Permutcol MATRIX AFTER"
311 : : do im1=1,tndim
312 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Permutcol(im1,im2),im2=1,tndim)
313 : : end do
314 : : endif
315 : :
316 [ # # # # : 0 : ABI_MALLOC(Apermutcol,(tndim,tndim))
# # # # #
# # # ]
317 : : if(prtopt==1) then
318 : : write(std_out,*) "Check product of original matrix by permutation matrix "
319 : : endif
320 [ # # # # ]: 0 : Apermutcol=zero
321 [ # # ]: 0 : do im1=1,tndim
322 [ # # ]: 0 : do im2=1,tndim
323 : 0 : Apermutcol(im1,im2)=zero
324 [ # # ]: 0 : do im3=1,tndim
325 : 0 : Apermutcol(im1,im2)=matrix_save(im1,im3)*Permutcol(im3,im2)+Apermutcol(im1,im2)
326 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
327 : : end do
328 : : end do
329 : : end do
330 : : if(prtopt==1) then
331 : : write(std_out,*) "Asave*Permutcol"
332 : : do im1=1,tndim
333 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcol(im1,im2),im2=1,tndim)
334 : : end do
335 : : endif
336 : :
337 [ # # # # : 0 : ABI_MALLOC(Permutline,(tndim,tndim))
# # # # #
# # # ]
338 [ # # # # ]: 0 : Permutline=zero
339 [ # # ]: 0 : do im1=1,tndim
340 : 0 : Permutline(im1,im1)=1.d0
341 : : end do
342 : :
343 [ # # ]: 0 : do im1=1,tndim
344 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
345 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
346 [ # # ]: 0 : do im2=1,tndim
347 : : ! write(std_out,*) "im1,im2, abs matrix(im2,im1),nonnul(im1)",im1,im2,abs(B(im1,im2)),nonnul(im1)
348 [ # # # # ]: 0 : if (abs(matrix(im2,im1))>0.000000000001.and.im2>nonnul(im1)) then
349 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
350 : 0 : nonnul(im1)=nonnul(im1)+1
351 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
352 : : ! permute
353 [ # # ]: 0 : do im3=1,tndim
354 : 0 : tmpx=matrix(im2,im3)
355 : 0 : matrix(im2,im3)=matrix(nonnul(im1),im3)
356 : 0 : matrix(nonnul(im1),im3)=tmpx
357 : 0 : tmpi=Permutline(im2,im3)
358 : 0 : Permutline(im2,im3)=Permutline(nonnul(im1),im3)
359 : 0 : Permutline(nonnul(im1),im3)=tmpi
360 : : enddo
361 [ # # ]: 0 : elseif (abs(matrix(im2,im1))<0.000000000001) then
362 : 0 : matrix(im2,im1)=zero
363 : : endif
364 : : enddo
365 : : enddo
366 : : if(prtopt==1) then
367 : : write(std_out,*) "matrix AFTER"
368 : : do im1=1,tndim
369 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
370 : : end do
371 : : write(std_out,*) "Permutline MATRIX AFTER"
372 : : do im1=1,tndim
373 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Permutline(im1,im2),im2=1,tndim)
374 : : end do
375 : : endif
376 : :
377 : : if(prtopt==1) then
378 : : write(std_out,*) "Check product of Apermutcol matrix by permutation matrix of the line "
379 : : endif
380 [ # # # # : 0 : ABI_MALLOC(Apermutline,(tndim,tndim))
# # # # #
# # # ]
381 [ # # # # ]: 0 : Apermutline=zero
382 [ # # ]: 0 : do im1=1,tndim
383 [ # # ]: 0 : do im2=1,tndim
384 : 0 : Apermutline(im1,im2)=zero
385 [ # # ]: 0 : do im3=1,tndim
386 : 0 : Apermutline(im1,im2)=Apermutcol(im3,im2)*Permutline(im1,im3)+Apermutline(im1,im2)
387 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
388 : : end do
389 : : end do
390 : : end do
391 : : if(prtopt==1) then
392 : : write(std_out,*) "Permutline*Apermutcol"
393 : : do im1=1,tndim
394 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutline(im1,im2),im2=1,tndim)
395 : : end do
396 : : endif
397 [ # # ]: 0 : work=0.d0
398 : 0 : call dsyev('v','u',tndim,matrix_save,tndim,eig,work,lwork,info)
399 [ # # ]: 0 : if(info/=0) then
400 : 0 : message = 'Error in diagonalization of matrix (dsyev) ! - '
401 : 0 : ABI_ERROR(message)
402 : : end if
403 : : if(prtopt==1) then
404 : : write(std_out,*) 'output',INFO
405 : : write(std_out,*) "Eigenvalues"
406 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
407 : : write(std_out,*) "Eigenvectors"
408 : : do im1=1,tndim
409 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix_save(im1,im2),im2=1,tndim)
410 : : end do
411 : : endif
412 : :
413 : :
414 : : ! call dsyev('v','u',tndim,A,LDA,W,WORKTMP,LWORK,INFO)
415 : : ! write(std_out,*) "optimal lwork",worktmp(1)
416 : :
417 [ # # ]: 0 : work=0.d0
418 : 0 : call dsyev('v','u',tndim,matrix,tndim,eig,work,lwork,info)
419 [ # # ]: 0 : if(info/=0) then
420 : 0 : message = 'Error in diagonalization of matrix (dsyev) ! - '
421 : 0 : ABI_ERROR(message)
422 : : end if
423 : : if(prtopt==1) then
424 : : write(std_out,*) 'output',INFO
425 : : write(std_out,*) "Eigenvalues"
426 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
427 : : write(std_out,*) "Eigenvectors"
428 : : do im1=1,tndim
429 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
430 : : end do
431 : : endif
432 : :
433 : :
434 : : !! REORDER EIGENVECTORS
435 [ # # # # : 0 : ABI_MALLOC(nonnuldege,(tndim))
# # # # ]
436 : 0 : newstarting=1
437 : 0 : current_dege=1
438 [ # # ]: 0 : do im4=2,tndim
439 [ # # ]: 0 : if(im4<tndim) testdege=((eig(im4)-eig(im4-1))<tol12)
440 [ # # ]: 0 : if(im4==tndim) then
441 : 0 : testdege=.false.
442 : 0 : current_dege=current_dege+1
443 : : endif
444 [ # # ]: 0 : if(testdege) then
445 : 0 : current_dege=current_dege+1
446 : : else
447 : : !new set of degenerate state: reorder it: put it into block diagonal
448 : : !form for column
449 : : if(prtopt==1) write(std_out,*) "newstarting, current_dege",newstarting, current_dege
450 : 0 : shift=0
451 [ # # ]: 0 : do im1=1,tndim ! balaye les premiers coefficients puis les autres
452 : :
453 : : ! if(im1==1) nonnuldege(im1)=0
454 : : ! if(im1>1) nonnuldege(im1)=nonnuldege(im1-1)
455 : 0 : maxvalue=0.00000001
456 : 0 : swap=.false.
457 [ # # ]: 0 : do im2=newstarting+shift,newstarting+current_dege-1
458 [ # # ]: 0 : if(abs(matrix(im1,im2))>maxvalue) then
459 : 0 : maxvalue=abs(matrix(im1,im2))
460 : 0 : indice_formax=im2
461 : 0 : swap=.true.
462 : : endif
463 : : enddo
464 : : ! found max value: permute
465 [ # # ]: 0 : if(swap) then
466 [ # # ]: 0 : do im3=1,tndim
467 : 0 : tmpx=matrix(im3,indice_formax)
468 : 0 : matrix(im3,indice_formax)=matrix(im3,newstarting+shift)
469 : 0 : matrix(im3,newstarting+shift)=tmpx
470 : : enddo
471 : 0 : shift=shift+1
472 : : endif
473 : : !write(std_out,*) "Eigenvectors after m1"
474 : : !do im3=1,tndim
475 : : ! write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im3,im5),im5=1,tndim)
476 : : !end do
477 : :
478 : : enddo
479 : : if(prtopt==1) then
480 : : write(std_out,*) "Eigenvectors after set of dege"
481 : : do im2=1,tndim
482 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im2,im3),im3=1,tndim)
483 : : end do
484 : : endif
485 : 0 : newstarting=im4
486 : 0 : current_dege=1
487 : : endif
488 : : enddo
489 [ # # ]: 0 : ABI_FREE(nonnuldege)
490 : : if(prtopt==1) then
491 : : write(std_out,*) "Ordered Eigenvectors"
492 : : do im1=1,tndim
493 : : write(std_out,'(2(1x,18(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
494 : : end do
495 : : endif
496 : :
497 : : if(prtopt==1) then
498 : : write(std_out,*) "inverse operation: reconstitute original matrix: only the line here"
499 : : endif
500 [ # # # # : 0 : ABI_MALLOC(Apermutlineback,(tndim,tndim))
# # # # #
# # # ]
501 [ # # # # ]: 0 : Apermutlineback=zero
502 [ # # ]: 0 : do im1=1,tndim
503 [ # # ]: 0 : do im2=1,tndim
504 : 0 : Apermutlineback(im1,im2)=zero
505 [ # # ]: 0 : do im3=1,tndim
506 : 0 : Apermutlineback(im1,im2)=matrix(im3,im2)*Permutline(im3,im1)+Apermutlineback(im1,im2)
507 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
508 : : end do
509 : : end do
510 : : end do
511 [ # # # # ]: 0 : matrix=Apermutlineback
512 : : if(prtopt==1) then
513 : : write(std_out,*) "t(Permutline)*Apermutcol"
514 : : do im1=1,tndim
515 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
516 : : end do
517 : : endif
518 : :
519 : : ! Now, set the first coefficient of eigenvectors positive.
520 : :
521 [ # # ]: 0 : do im2=1,tndim ! loop over eigenvectors
522 [ # # ]: 0 : do im1=1,tndim ! loop over components
523 [ # # ]: 0 : if(abs(matrix(im1,im2))>tol8) then
524 [ # # ]: 0 : if(matrix(im1,im2)<0) then
525 [ # # ]: 0 : do im3=1,tndim
526 [ # # ]: 0 : if(abs(matrix(im3,im2))>tol8) then
527 : 0 : matrix(im3,im2)=-matrix(im3,im2)
528 : : endif
529 : : enddo
530 : : endif
531 : : exit
532 : : endif
533 : : enddo
534 : : enddo
535 : : if(prtopt==1) then
536 : : write(std_out,*) "Impose first component of eigenvectors is positive"
537 : : do im1=1,tndim
538 : : write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
539 : : end do
540 : : endif
541 : : ! write(std_out,*) "inverse operation: reconstitute original matrix: then the column"
542 : : ! Apermutcolback=zero
543 : : ! do im1=1,tndim
544 : : ! do im2=1,tndim
545 : : ! Apermutcolback(im1,im2)=zero
546 : : ! do im3=1,tndim
547 : : ! Apermutcolback(im1,im2)=Apermutlineback(im1,im3)*Permutcol(im2,im3)+Apermutcolback(im1,im2)
548 : : ! ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
549 : : ! end do
550 : : ! end do
551 : : ! end do
552 : : ! write(std_out,*) "Apermutlineback*t(Permutcol)"
553 : : ! do im1=1,10
554 : : ! write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcolback(im1,im2),im2=1,10)
555 : : ! end do
556 : :
557 : :
558 [ # # ]: 0 : ABI_FREE(Apermutlineback)
559 [ # # ]: 0 : ABI_FREE(Apermutline)
560 [ # # ]: 0 : ABI_FREE(matrix_save)
561 [ # # ]: 0 : ABI_FREE(Apermutcol)
562 [ # # ]: 0 : ABI_FREE(work)
563 [ # # ]: 0 : ABI_FREE(Permutcol)
564 [ # # ]: 0 : ABI_FREE(nonnul)
565 [ # # ]: 0 : ABI_FREE(Permutline)
566 : :
567 : : DBG_EXIT("COLL")
568 : :
569 : 0 : end subroutine blockdiago_fordsyev
570 : : !!***
571 : :
572 : : !! FUNCTION
573 : : !! Transform matrix into block diagonal form before diagonalisation
574 : : !!
575 : : !! INPUTS
576 : : !! ndim = dimension of matrix
577 : : !! matrix= matrix
578 : : !!
579 : : !! OUTPUT
580 : : !! matrix= square root of the matrix
581 : : !!
582 : : !! SOURCE
583 : :
584 : 0 : subroutine blockdiago_forzheev(matrix,tndim,eig)
585 : :
586 : : !Arguments ------------------------------------
587 : : !scalars
588 : : integer,intent(in) :: tndim
589 : : complex(dp),intent(inout) :: matrix(tndim,tndim)
590 : : real(dp),intent(inout) :: eig(tndim)
591 : : !arrays
592 : :
593 : : !Local variables-------------------------------
594 : : !scalars
595 : : integer :: im1,im2,im3,info,lwork
596 : : character(len=500) :: message
597 : : complex(dp):: tmpx
598 : : integer(dp):: tmpi,prtopt
599 : : !arrays
600 : 0 : real(dp),allocatable :: rwork(:)
601 : 0 : complex(dp),allocatable :: work(:)
602 : 0 : real(dp),allocatable :: Permutcol(:,:)
603 : 0 : complex(dp),allocatable :: Apermutcol(:,:)
604 : 0 : complex(dp),allocatable :: Apermutline(:,:)
605 : 0 : complex(dp),allocatable :: Apermutlineback(:,:)
606 : 0 : real(dp),allocatable :: Permutline(:,:)
607 : 0 : complex(dp),allocatable :: matrix_save(:,:) !,W(:)
608 : 0 : integer,allocatable :: nonnul(:)
609 : : ! *************************************************************************
610 : :
611 : : !!!Do not remove this silly print instruction. Seems needed to avoid floating
612 : : !!!point exception on vm1_gcc51 ...
613 : : !!#if __GFORTRAN__ == 1 && __GNUC__ == 5 && (__GNUC_MINOR__ == 1 || __GNUC_MINOR__ == 2)
614 : : !! write(std_out,'(a)')' invsqrt_matrix at m_matrix.F90 : enter ( needed to avoid FPE with GCC5[1,2] )'
615 : : !!#endif
616 : : DBG_ENTER("COLL")
617 : :
618 : : !printing option for debug
619 : 0 : prtopt = 0
620 : :
621 : 0 : lwork=10*tndim
622 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
623 [ # # # # : 0 : ABI_MALLOC(rwork,(3*tndim-2))
# # # # #
# # # #
# ]
624 : :
625 [ # # # # : 0 : ABI_MALLOC(matrix_save,(tndim,tndim))
# # # # #
# # # # #
# # ]
626 [ # # # # : 0 : matrix_save=matrix
# # ]
627 : :
628 [ # # # # : 0 : ABI_MALLOC(Permutcol,(tndim,tndim))
# # # # #
# # # # #
# # ]
629 : :
630 [ # # # # ]: 0 : Permutcol=zero
631 [ # # ]: 0 : do im1=1,tndim
632 : 0 : Permutcol(im1,im1)=1.d0
633 : : end do
634 : :
635 : : if(prtopt == 1) then
636 : : write(std_out,*) "MATRIX"
637 : : do im1=1,tndim
638 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix_save(im1,im2),im2=1,tndim)
639 : : end do
640 : : endif
641 : :
642 [ # # # # : 0 : ABI_MALLOC(nonnul,(tndim))
# # # # #
# # # ]
643 [ # # ]: 0 : do im1=1,tndim
644 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
645 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
646 [ # # ]: 0 : do im2=1,tndim
647 [ # # # # ]: 0 : if (abs(matrix(im1,im2))>0.000000000001.and.im2>nonnul(im1)) then
648 : 0 : nonnul(im1)=nonnul(im1)+1
649 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
650 : : ! permute
651 [ # # ]: 0 : do im3=1,tndim
652 : 0 : tmpx=matrix(im3,im2)
653 : 0 : matrix(im3,im2)=matrix(im3,nonnul(im1))
654 : 0 : matrix(im3,nonnul(im1))=tmpx
655 : 0 : tmpi=Permutcol(im3,im2)
656 : 0 : Permutcol(im3,im2)=Permutcol(im3,nonnul(im1))
657 : 0 : Permutcol(im3,nonnul(im1))=tmpi
658 : : enddo
659 [ # # ]: 0 : elseif (abs(matrix(im1,im2))<0.000000000001) then
660 : 0 : matrix(im1,im2)=czero
661 : : endif
662 : : enddo
663 : : enddo
664 : :
665 : : if (prtopt == 1) then
666 : : write(std_out,*) "MATRIX AFTER COLUMN PERMUT"
667 : : do im1=1,tndim
668 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
669 : : end do
670 : : write(std_out,*) "Permutcol MATRIX AFTER"
671 : : do im1=1,tndim
672 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Permutcol(im1,im2),im2=1,tndim)
673 : : end do
674 : : endif
675 : :
676 [ # # # # : 0 : ABI_MALLOC(Apermutcol,(tndim,tndim))
# # # # #
# # # ]
677 [ # # # # ]: 0 : Apermutcol=czero
678 [ # # ]: 0 : do im1=1,tndim
679 [ # # ]: 0 : do im2=1,tndim
680 : 0 : Apermutcol(im1,im2)=czero
681 [ # # ]: 0 : do im3=1,tndim
682 : 0 : Apermutcol(im1,im2)=matrix_save(im1,im3)*Permutcol(im3,im2)+Apermutcol(im1,im2)
683 : : ! write(std_out,*) im1,im2,im3,Apermutcol(im1,im2)
684 : : end do
685 : : end do
686 : : end do
687 : :
688 : : if(prtopt == 1) then
689 : : write(std_out,*) "Check product of original matrix by permutation matrix "
690 : : write(std_out,*) "Asave*Permutcol"
691 : : do im1=1,tndim
692 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Apermutcol(im1,im2),im2=1,tndim)
693 : : end do
694 : : endif
695 : :
696 : :
697 [ # # # # : 0 : ABI_MALLOC(Permutline,(tndim,tndim))
# # # # #
# # # ]
698 [ # # # # ]: 0 : Permutline=zero
699 [ # # ]: 0 : do im1=1,tndim
700 : 0 : Permutline(im1,im1)=1.d0
701 : : end do
702 : :
703 [ # # ]: 0 : do im1=1,tndim
704 [ # # ]: 0 : if(im1==1) nonnul(im1)=0
705 [ # # ]: 0 : if(im1>1) nonnul(im1)=nonnul(im1-1)
706 [ # # ]: 0 : do im2=1,tndim
707 : : ! write(std_out,*) "im1,im2, abs matrix(im2,im1),nonnul(im1)",im1,im2,abs(B(im1,im2)),nonnul(im1)
708 [ # # # # ]: 0 : if (abs(matrix(im2,im1))>0.000000000001.and.im2>nonnul(im1)) then
709 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
710 : 0 : nonnul(im1)=nonnul(im1)+1
711 : : ! write(std_out,*) "im2,nonnul(im1)",im2,nonnul(im1)
712 : : ! permute
713 [ # # ]: 0 : do im3=1,tndim
714 : 0 : tmpx=matrix(im2,im3)
715 : 0 : matrix(im2,im3)=matrix(nonnul(im1),im3)
716 : 0 : matrix(nonnul(im1),im3)=tmpx
717 : 0 : tmpi=Permutline(im2,im3)
718 : 0 : Permutline(im2,im3)=Permutline(nonnul(im1),im3)
719 : 0 : Permutline(nonnul(im1),im3)=tmpi
720 : : enddo
721 [ # # ]: 0 : elseif (abs(matrix(im2,im1))<0.000000000001) then
722 : 0 : matrix(im2,im1)=czero
723 : : endif
724 : : enddo
725 : : enddo
726 : :
727 : : if(prtopt == 1) then
728 : : write(std_out,*) "matrix AFTER"
729 : : do im1=1,tndim
730 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
731 : : end do
732 : : write(std_out,*) "Permutline MATRIX AFTER"
733 : : do im1=1,tndim
734 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Permutline(im1,im2),im2=1,tndim)
735 : : end do
736 : : endif
737 : :
738 [ # # # # : 0 : ABI_MALLOC(Apermutline,(tndim,tndim))
# # # # #
# # # ]
739 [ # # # # ]: 0 : Apermutline=czero
740 [ # # ]: 0 : do im1=1,tndim
741 [ # # ]: 0 : do im2=1,tndim
742 : 0 : Apermutline(im1,im2)=czero
743 [ # # ]: 0 : do im3=1,tndim
744 : 0 : Apermutline(im1,im2)=Apermutcol(im3,im2)*Permutline(im1,im3)+Apermutline(im1,im2)
745 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
746 : : end do
747 : : end do
748 : : end do
749 : :
750 : : if(prtopt == 1) then
751 : : write(std_out,*) "Check product of Apermutcol matrix by permutation matrix of the line "
752 : : write(std_out,*) "Permutline*Apermutcol"
753 : : do im1=1,tndim
754 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (Apermutline(im1,im2),im2=1,tndim)
755 : : end do
756 : : endif
757 : :
758 [ # # ]: 0 : work=czero
759 : 0 : call zheev('v','u',tndim,matrix_save,tndim,eig,work,lwork,rwork,info)
760 [ # # ]: 0 : if(info/=0) then
761 : 0 : message = 'Error in diagonalization of matrix (zheev) ! - '
762 : 0 : ABI_ERROR(message)
763 : : end if
764 : :
765 : : if(prtopt == 1) then
766 : : write(std_out,*) 'output',INFO
767 : : write(std_out,*) "Eigenvalues"
768 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
769 : : write(std_out,*) "Eigenvectors"
770 : : do im1=1,tndim
771 : : write(std_out,'(2(1x,30(1x,f20.15,f20.15)))') (matrix_save(im1,im2),im2=1,tndim)
772 : : end do
773 : : endif
774 : :
775 : : ! call dsyev('v','u',tndim,A,LDA,W,WORKTMP,LWORK,INFO)
776 : : ! write(std_out,*) "optimal lwork",worktmp(1)
777 : :
778 [ # # ]: 0 : work=czero
779 : 0 : call zheev('v','u',tndim,matrix,tndim,eig,work,lwork,rwork,info)
780 [ # # ]: 0 : if(info/=0) then
781 : 0 : message = 'Error in diagonalization of matrix (zheev) ! - '
782 : 0 : ABI_ERROR(message)
783 : : end if
784 : :
785 : : if(prtopt == 1) then
786 : : write(std_out,*) 'output',INFO
787 : : write(std_out,*) "Eigenvalues"
788 : : write(std_out,'(2x,20f20.15) ') (eig(im1),im1=1,tndim)
789 : : write(std_out,*) "Eigenvectors"
790 : : do im1=1,tndim
791 : : write(std_out,'(2(1x,30(1x,f20.15,f20.15)))') (matrix(im1,im2),im2=1,tndim)
792 : : end do
793 : : endif
794 : :
795 [ # # # # : 0 : ABI_MALLOC(Apermutlineback,(tndim,tndim))
# # # # #
# # # ]
796 [ # # # # ]: 0 : Apermutlineback=czero
797 [ # # ]: 0 : do im1=1,tndim
798 [ # # ]: 0 : do im2=1,tndim
799 : 0 : Apermutlineback(im1,im2)=czero
800 [ # # ]: 0 : do im3=1,tndim
801 : 0 : Apermutlineback(im1,im2)=matrix(im3,im2)*Permutline(im3,im1)+Apermutlineback(im1,im2)
802 : : ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
803 : : end do
804 : : end do
805 : : end do
806 [ # # # # ]: 0 : matrix=Apermutlineback
807 : :
808 : : if(prtopt == 1) then
809 : : write(std_out,*) "inverse operation: reconstitute original matrix: first the line"
810 : : write(std_out,*) "t(Permutline)*Apermutcol"
811 : : do im1=1,tndim
812 : : write(std_out,'(2(1x,30(1x,f22.18,f22.18)))') (matrix(im1,im2),im2=1,tndim)
813 : : end do
814 : : endif
815 : :
816 : : ! write(std_out,*) "inverse operation: reconstitute original matrix: then the column"
817 : : ! Apermutcolback=zero
818 : : ! do im1=1,tndim
819 : : ! do im2=1,tndim
820 : : ! Apermutcolback(im1,im2)=zero
821 : : ! do im3=1,tndim
822 : : ! Apermutcolback(im1,im2)=Apermutlineback(im1,im3)*Permutcol(im2,im3)+Apermutcolback(im1,im2)
823 : : ! ! write(std_out,*) PROD(im1,im2),A(im3,im1),A(im3,im2)
824 : : ! end do
825 : : ! end do
826 : : ! end do
827 : : ! write(std_out,*) "Apermutlineback*t(Permutcol)"
828 : : ! do im1=1,10
829 : : ! write(std_out,'(2(1x,18(1x,f22.18,f22.18)))') (Apermutcolback(im1,im2),im2=1,10)
830 : : ! end do
831 : :
832 : :
833 [ # # ]: 0 : ABI_FREE(Apermutlineback)
834 [ # # ]: 0 : ABI_FREE(Apermutline)
835 [ # # ]: 0 : ABI_FREE(matrix_save)
836 [ # # ]: 0 : ABI_FREE(Apermutcol)
837 [ # # ]: 0 : ABI_FREE(work)
838 [ # # ]: 0 : ABI_FREE(Permutcol)
839 [ # # ]: 0 : ABI_FREE(nonnul)
840 [ # # ]: 0 : ABI_FREE(Permutline)
841 [ # # ]: 0 : ABI_FREE(rwork)
842 : :
843 : : DBG_EXIT("COLL")
844 : :
845 : 0 : end subroutine blockdiago_forzheev
846 : : !!***
847 : :
848 : : !! FUNCTION
849 : : !! Compute the determinant of a 3x3 real matrix
850 : : !!
851 : : !! INPUTS
852 : : !! A = 3x3 matrix
853 : : !!
854 : : !! OUTPUT
855 : : !! det = The determinant
856 : : !!
857 : : !! SOURCE
858 : :
859 : 44 : function real_mat33det(A) result(det)
860 : : real(dp), intent(in) :: A(3,3)
861 : : real(dp) :: det
862 : : DET = A(1,1)*A(2,2)*A(3,3) &
863 : : - A(1,1)*A(2,3)*A(3,2) &
864 : : - A(1,2)*A(2,1)*A(3,3) &
865 : : + A(1,2)*A(2,3)*A(3,1) &
866 : : + A(1,3)*A(2,1)*A(3,2) &
867 : 44 : - A(1,3)*A(2,2)*A(3,1)
868 : 44 : end function real_mat33det
869 : : !!***
870 : :
871 : : !! FUNCTION
872 : : !! Compute the determinant of a 3x3 integer matrix
873 : : !!
874 : : !! INPUTS
875 : : !! A = 3x3 matrix
876 : : !!
877 : : !! OUTPUT
878 : : !! det = The determinant
879 : : !!
880 : : !! SOURCE
881 : :
882 : 6 : function int_mat33det(A) result(det)
883 : : integer, intent(in) :: A(3,3)
884 : : integer :: det
885 : : DET = A(1,1)*A(2,2)*A(3,3) &
886 : : - A(1,1)*A(2,3)*A(3,2) &
887 : : - A(1,2)*A(2,1)*A(3,3) &
888 : : + A(1,2)*A(2,3)*A(3,1) &
889 : : + A(1,3)*A(2,1)*A(3,2) &
890 : 6 : - A(1,3)*A(2,2)*A(3,1)
891 : 6 : end function int_mat33det
892 : : !!***
893 : :
894 : : !!****f* m_matrix/mati3inv
895 : : !! NAME
896 : : !! mati3inv
897 : : !!
898 : : !! FUNCTION
899 : : !! Invert and transpose orthogonal 3x3 matrix of INTEGER elements.
900 : : !!
901 : : !! INPUTS
902 : : !! mm = integer matrix to be inverted
903 : : !!
904 : : !! OUTPUT
905 : : !! mit = inverse of mm input matrix
906 : : !!
907 : : !! NOTES
908 : : !! Used for symmetry operations.
909 : : !! This routine applies to ORTHOGONAL matrices only.
910 : : !! Since these form a group, inverses are also integer arrays.
911 : : !! Returned array is TRANSPOSE of inverse, as needed.
912 : : !! Note use of integer arithmetic.
913 : : !!
914 : : !! SOURCE
915 : :
916 : 10486793 : subroutine mati3inv(mm, mit)
917 : :
918 : : !Arguments ------------------------------------
919 : : !arrays
920 : : integer,intent(in) :: mm(3,3)
921 : : integer,intent(out) :: mit(3,3)
922 : :
923 : : !Local variables-------------------------------
924 : : !scalars
925 : : integer :: dd
926 : : character(len=500) :: msg
927 : : !arrays
928 : : integer :: tt(3,3)
929 : : ! *************************************************************************
930 : :
931 : 10486793 : tt(1,1) = mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)
932 : 10486793 : tt(2,1) = mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)
933 : 10486793 : tt(3,1) = mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3)
934 : 10486793 : tt(1,2) = mm(3,1) * mm(2,3) - mm(2,1) * mm(3,3)
935 : 10486793 : tt(2,2) = mm(1,1) * mm(3,3) - mm(3,1) * mm(1,3)
936 : 10486793 : tt(3,2) = mm(2,1) * mm(1,3) - mm(1,1) * mm(2,3)
937 : 10486793 : tt(1,3) = mm(2,1) * mm(3,2) - mm(3,1) * mm(2,2)
938 : 10486793 : tt(2,3) = mm(3,1) * mm(1,2) - mm(1,1) * mm(3,2)
939 : 10486793 : tt(3,3) = mm(1,1) * mm(2,2) - mm(2,1) * mm(1,2)
940 : 10486793 : dd = mm(1,1) * tt(1,1) + mm(2,1) * tt(2,1) + mm(3,1) * tt(3,1)
941 : :
942 : : ! Make sure matrix is not singular
943 [ + - ]: 10486793 : if (dd /= 0) then
944 [ + + + + ]: 136328309 : mit(:,:)=tt(:,:)/dd
945 : : else
946 : 0 : write(msg, '(2a,2x,9(i0,1x),a)' )'Attempting to invert integer array',ch10,mm,' ==> determinant is zero.'
947 : 0 : ABI_ERROR(msg)
948 : : end if
949 : :
950 : : ! If matrix is orthogonal, determinant must be 1 or -1
951 [ - + ]: 10486793 : if (abs(dd) /= 1) then
952 : 0 : write(msg, '(3a,i0)' )'Absolute value of determinant should be one',ch10,'but determinant= ',dd
953 : 0 : ABI_ERROR(msg)
954 : : end if
955 : :
956 : 10486793 : end subroutine mati3inv
957 : : !!***
958 : :
959 : : !!****f* m_matrix/mati3det
960 : : !! NAME
961 : : !! mati3det
962 : : !!
963 : : !! FUNCTION
964 : : !! Compute the determinant of a 3x3 matrix of INTEGER elements.
965 : : !!
966 : : !! INPUTS
967 : : !! mm = integer matrix
968 : : !!
969 : : !! OUTPUT
970 : : !! det = determinant of the matrix
971 : : !!
972 : : !! SOURCE
973 : :
974 : 2755784 : subroutine mati3det(mm, det)
975 : :
976 : : !Arguments ------------------------------------
977 : : integer,intent(in) :: mm(3,3)
978 : : integer,intent(out) :: det
979 : : ! *************************************************************************
980 : :
981 : : det=mm(1,1)*(mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)) &
982 : : + mm(2,1)*(mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)) &
983 : 2755784 : + mm(3,1)*(mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3))
984 : :
985 : 2755784 : end subroutine mati3det
986 : : !!***
987 : :
988 : : !!****f* m_matrix/matr3inv
989 : : !! NAME
990 : : !! matr3inv
991 : : !!
992 : : !! FUNCTION
993 : : !! Invert and transpose general 3x3 matrix of real*8 elements.
994 : : !!
995 : : !! INPUTS
996 : : !! aa = 3x3 matrix to be inverted
997 : : !!
998 : : !! OUTPUT
999 : : !! ait = inverse of aa input matrix
1000 : : !!
1001 : : !! NOTES
1002 : : !! Returned array is TRANSPOSE of inverse, as needed to get g from r.
1003 : : !!
1004 : : !! SOURCE
1005 : :
1006 : 4147325 : subroutine matr3inv(aa, ait)
1007 : :
1008 : : !Arguments ------------------------------------
1009 : : !arrays
1010 : : real(dp),intent(in) :: aa(3,3)
1011 : : real(dp),intent(out) :: ait(3,3)
1012 : :
1013 : : !Local variables-------------------------------
1014 : : !scalars
1015 : : real(dp) :: dd,det,t1,t2,t3
1016 : : character(len=500) :: msg
1017 : : ! *************************************************************************
1018 : :
1019 : 4147325 : t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
1020 : 4147325 : t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
1021 : 4147325 : t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
1022 : 4147325 : det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
1023 : :
1024 : : !Make sure matrix is not singular
1025 [ + - ]: 4147325 : if (abs(det)>tol16) then
1026 : 4147325 : dd=one/det
1027 : : else
1028 : : write(msg, '(2a,2x,9es16.8,a,a,es16.8,a)' )&
1029 : 0 : 'Attempting to invert real(8) 3x3 array',ch10,aa(:,:),ch10,' ==> determinant=',det,' is zero.'
1030 : 0 : ABI_BUG(msg)
1031 : : end if
1032 : :
1033 : 4147325 : ait(1,1) = t1 * dd
1034 : 4147325 : ait(2,1) = t2 * dd
1035 : 4147325 : ait(3,1) = t3 * dd
1036 : 4147325 : ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
1037 : 4147325 : ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
1038 : 4147325 : ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
1039 : 4147325 : ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
1040 : 4147325 : ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
1041 : 4147325 : ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
1042 : :
1043 : 4147325 : end subroutine matr3inv
1044 : : !!***
1045 : :
1046 : : !!****f* m_matrix/is_unitary
1047 : : !! NAME
1048 : : !! is_unitary
1049 : : !!
1050 : : !! FUNCTION
1051 : : !! Returns .TRUE. if the input complex matrix mat is unitary, i.e. mat^\dagger mat = I.
1052 : : !! Also returns the maximum absolute deviation from the identity matrix.
1053 : : !!
1054 : : !! INPUTS
1055 : : !! nn = Dimension of the matrix.
1056 : : !! mat = Complex matrix of size (nn,nn).
1057 : : !! tol = Tolerance for the maximum deviation.
1058 : : !!
1059 : : !! OUTPUT
1060 : : !! .TRUE. if err < tol, .FALSE. otherwise.
1061 : : !! err = Maximum absolute deviation from the identity matrix.
1062 : : !!
1063 : : !! SOURCE
1064 : :
1065 : 0 : logical function is_unitary(nn, mat, tol, err)
1066 : :
1067 : : !Arguments ------------------------------------
1068 : : integer, intent(in) :: nn
1069 : : complex(dp), intent(in) :: mat(nn,nn)
1070 : : real(dp), intent(in) :: tol
1071 : : real(dp), intent(out) :: err
1072 : :
1073 : : !Local variables-------------------------------
1074 : : integer :: ii
1075 : 0 : complex(dp) :: prod(nn,nn), identity(nn,nn)
1076 : : !----------------------------------------------------------------------
1077 : :
1078 : : ! Compute mat^\dagger mat
1079 [ # # # # : 0 : prod = matmul(conjg(transpose(mat)), mat)
# # # # #
# # # # #
# # # # ]
1080 : :
1081 : : ! Build identity matrix
1082 [ # # # # ]: 0 : identity = czero
1083 [ # # ]: 0 : do ii=1,nn
1084 : 0 : identity(ii,ii) = one
1085 : : end do
1086 : :
1087 : : ! Maximum deviation from identity
1088 [ # # # # : 0 : err = maxval(abs(prod - identity))
# # # # #
# # # ]
1089 : 0 : is_unitary = (err < tol)
1090 : :
1091 : 0 : end function is_unitary
1092 : : !!***
1093 : :
1094 : : !!****f* m_matrix/is_identity
1095 : : !! NAME
1096 : : !! is_identity
1097 : : !!
1098 : : !! FUNCTION
1099 : : !! Returns .TRUE. if the complex input matrix mat is the identity matrix.
1100 : : !! Also returns the maximum absolute deviation from the identity matrix.
1101 : : !!
1102 : : !! INPUTS
1103 : : !! nn = Dimension of the matrix.
1104 : : !! mat = Complex matrix of size (nn,nn).
1105 : : !! tol = Tolerance for the maximum deviation.
1106 : : !!
1107 : : !! OUTPUT
1108 : : !! .TRUE. if err < tol, .FALSE. otherwise.
1109 : : !! err = Maximum absolute deviation from the identity matrix.
1110 : : !!
1111 : : !! SOURCE
1112 : :
1113 : 0 : logical function is_identity(nn, mat, tol, err)
1114 : :
1115 : : !Arguments ------------------------------------
1116 : : integer, intent(in) :: nn
1117 : : complex(dp), intent(in) :: mat(nn,nn)
1118 : : real(dp), intent(in) :: tol
1119 : : real(dp), intent(out) :: err
1120 : :
1121 : : !Local variables-------------------------------
1122 : : integer :: ii
1123 : 0 : complex(dp) :: identity(nn,nn)
1124 : : !----------------------------------------------------------------------
1125 : :
1126 : : ! Build identity matrix
1127 [ # # # # ]: 0 : identity = czero
1128 [ # # ]: 0 : do ii=1,nn
1129 : 0 : identity(ii,ii) = one
1130 : : end do
1131 : :
1132 : : ! Maximum deviation from identity
1133 [ # # # # : 0 : err = maxval(abs(mat - identity))
# # # # #
# # # ]
1134 : 0 : is_identity = (err < tol)
1135 : :
1136 : 0 : end function is_identity
1137 : : !!***
1138 : :
1139 : : end module m_matrix
|
