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 10467278 : 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 10467278 : tt(1,1) = mm(2,2) * mm(3,3) - mm(3,2) * mm(2,3)
932 10467278 : tt(2,1) = mm(3,2) * mm(1,3) - mm(1,2) * mm(3,3)
933 10467278 : tt(3,1) = mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3)
934 10467278 : tt(1,2) = mm(3,1) * mm(2,3) - mm(2,1) * mm(3,3)
935 10467278 : tt(2,2) = mm(1,1) * mm(3,3) - mm(3,1) * mm(1,3)
936 10467278 : tt(3,2) = mm(2,1) * mm(1,3) - mm(1,1) * mm(2,3)
937 10467278 : tt(1,3) = mm(2,1) * mm(3,2) - mm(3,1) * mm(2,2)
938 10467278 : tt(2,3) = mm(3,1) * mm(1,2) - mm(1,1) * mm(3,2)
939 10467278 : tt(3,3) = mm(1,1) * mm(2,2) - mm(2,1) * mm(1,2)
940 10467278 : 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 10467278 : if (dd /= 0) then
944 136074614 : 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 10467278 : 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 10467278 : 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 2735041 : 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 2735041 : + mm(3,1)*(mm(1,2) * mm(2,3) - mm(2,2) * mm(1,3))
984 :
985 2735041 : 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 4107641 : 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 4107641 : t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
1020 4107641 : t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
1021 4107641 : t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
1022 4107641 : det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
1023 :
1024 : !Make sure matrix is not singular
1025 4107641 : if (abs(det)>tol16) then
1026 4107641 : 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 4107641 : ait(1,1) = t1 * dd
1034 4107641 : ait(2,1) = t2 * dd
1035 4107641 : ait(3,1) = t3 * dd
1036 4107641 : ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
1037 4107641 : ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
1038 4107641 : ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
1039 4107641 : ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
1040 4107641 : ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
1041 4107641 : ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
1042 :
1043 4107641 : 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
|
