Talk:Jacobi method for complex Hermitian matrices
![]() | This article is rated Start-class on Wikipedia's content assessment scale. It is of interest to the following WikiProjects: | ||||||||||
|
[Untitled]
editThis article seems to be original research which violated Wikipedia rules (See Wikipedia:No_original_research). The cited reference by Strang does not explain the Jacobi Method for Hermitian matrices (only for real symmetric ones).
A few more notes:
- for Hermitian matrices it would be more appropriate to use unitary transformations rather than two rotations.
- the derivation of the angles is missing
MSoegtropI (talk) 12:16, 1 March 2017 (UTC)
A master student:
I tried to derive the equations for the angles, I came up with the results tan(phi2)=2*abs(Hpq)/(Hqq-Hpp). This means there ought to be a missing minus sign. The second correction should be theta1=(2*phi1+pi)/4, I replaced the minus sign with a positive one. I confirmed these results by programming the Jacobi-Sweep algorithm using these conditions (Fortran 2003 code using gfortran compiler), I receive in every case proper eigenvalues and in many cases (maybe there are a few bugs left) proper eigenvectors:
program get_complex_diagonal implicit none !Constants: integer, parameter:: kr15=selected_real_kind(15) real (kr15), parameter:: pi=4.0_kr15*datan(1.0_kr15), err_tol=0.000000000001_kr15 !Varibales complex (kr15),dimension(:), allocatable :: Dtr complex (kr15),dimension(:,:), allocatable::A,Id,U,Utot,D,Doff,Dtemp,Utemp
real (kr15) :: Ddiff,theta1,theta2,phi1,phi2, Serr real (kr15) :: abDpq,arg1,arg2,rDpq,imDpq
integer :: idx, N, ct=1 integer, dimension(2) :: pq integer, dimension(:), allocatable :: lbl write(*,*) 'Determine the dimension of the N by N hermitian matrix A:'
write(*,*) 'N=' read(*,*) N write(*,'(A)') 'Now give the entries rowwise for A:'
allocate(A(1:N,1:N))
do idx=1, N write(*,'(i5,A)') idx , '. row:' read(*,*) A(idx,:) enddo allocate(Id(1:N,1:N))
!Construction of identity matrix:
Id=0.0_kr15 do idx=1,N
Id(idx,idx)=(1.0_kr15, 0.0_kr15)
enddo
allocate(U(1:N,1:N),Utot(1:N,1:N),D(1:N,1:N),Doff(1:N,1:N))
U=Id Utot=U !initiating D and Doff: D=A Doff=D-(Id*D) !initial phi2s set to zero: phi1=0.0_kr15 phi2=phi1 theta1=phi1 theta2=phi1
Serr=sum(dconjg(Doff)*Doff)
write(*,*) Serr
!!Main loop:
!write(*,*) 'Doff' !write(*,*) dabs(Doff)
do while (Serr .ne. 0.0_kr15)
pq=maxloc(dreal(dconjg(Doff)*Doff))
pq=(/minval(pq),maxval(pq)/)
write(*,*) pq
!real, imaginary and absolute of Dpq: rDpq =real(D(pq(1),pq(2)),kr15) imDpq=dimag(D(pq(1),pq(2))) abDpq=dsqrt(rDpq*rDpq+imDpq*imDpq) arg1=imDpq/rDpq
!phi and theta 1: if (imDpq .eq. 0.0_kr15 .and. rDpq .gt. 0.0_kr15) then phi1= 0.0_kr15 !write(*,*) 'phi1=0', angle elseif (imDpq .gt. 0.0_kr15 .and. rDpq.gt. 0.0_kr15 ) then phi1= datan(arg1) !write(*,*) '0<phi1<90', angle elseif (imDpq .gt. 0.0_kr15 .and. rDpq .eq. 0.0_kr15) then phi1= 0.5_kr15*pi !write(*,*) 'phi1=90', angle elseif (imDpq .gt. 0.0_kr15 .and. rDpq .lt. 0.0_kr15) then phi1= datan(arg1)+pi !write(*,*) '90<phi1<180', angle elseif (imDpq .eq. 0.0_kr15 .and. rDpq.lt. 0.0_kr15 ) then phi1= pi !write(*,*) 'phi1=180' , angle elseif (imDpq .lt. 0.0_kr15 .and. rDpq .lt. 0.0_kr15) then phi1= datan(arg1)+pi !write(*,*) '180<phi1<270', angle elseif (imDpq .lt. 0.0_kr15 .and. rDpq .eq. 0.0_kr15 ) then phi1= 1.5_kr15*pi !write(*,*) 'phi1=270', angle elseif (imDpq .lt. 0.0_kr15 .and. rDpq .gt. 0.0_kr15 ) then phi1= datan(arg1)+2.0_kr15*pi !write(*,*) '270<phi1<360', angle
else write(*,*) 'Error in the definition of phi1!' endif ! theta 1: theta1=0.25_kr15*(2.0_kr15*phi1+pi)
! tan(phi2) Ddiff=real(D(pq(2),pq(2))-D(pq(1),pq(1)),kr15) ! Dqq-Dpp arg2=(2.0_kr15*abDpq)/(Ddiff)
!phi and theta 2: if (abDpq .eq. 0.0_kr15 .and. Ddiff .gt. 0.0_kr15 ) then phi2= 0.0_kr15 !write(*,*) 'phi2=0', phi2 elseif (abDpq .gt. 0.0_kr15 .and. Ddiff.gt. 0.0_kr15 ) then phi2= datan(arg2) !write(*,*) '0<phi2<90', phi2 elseif (abDpq .gt. 0.0_kr15 .and. Ddiff .eq. 0.0_kr15 ) then phi2= 0.5_kr15*pi !write(*,*) 'phi2=90', phi2 elseif (abDpq .gt. 0.0_kr15 .and. Ddiff.lt. 0.0_kr15 ) then phi2= datan(arg2)+pi !write(*,*) '90<phi2<180', phi2 elseif (abDpq .eq. 0.0_kr15 .and. Ddiff.lt. 0.0_kr15 ) then phi2= pi !write(*,*) 'phi2=180' , phi2 elseif (abDpq .lt. 0.0_kr15 .and. Ddiff.lt. 0.0_kr15 ) then phi2= datan(arg2)+pi !write(*,*) '180<phi2<270', phi2 elseif (abDpq .lt. 0.0_kr15 .and. Ddiff.eq. 0.0_kr15 ) then phi2= 1.5_kr15*pi !write(*,*) 'phi2=270', phi2 elseif (abDpq .lt. 0.0_kr15 .and. Ddiff.gt. 0.0_kr15 ) then phi2= datan(arg2)+2.0_kr15*pi !write(*,*) '270<phi2<360', phi2 else write(*,*) 'Error in the definition of phi2!' endif !theta2: theta2=0.5_kr15*phi2
!Definition of Unitary matrix:
U=Id
U(pq(1),pq(1))= zexp(dcmplx(0.0_kr15,theta1))*dcmplx(0.0_kr15,dsin(theta2))
U(pq(2),pq(2))= dconjg(U(pq(1),pq(1)))
U(pq(1),pq(2))= zexp(dcmplx(0.0_kr15,theta1))*dcmplx(dcos(theta2),0.0_kr15)
U(pq(2),pq(1))= zexp(dcmplx(0.0_kr15,-theta1))*dcmplx(-dcos(theta2),0.0_kr15)
!Actual transformation: D=matmul(transpose(dconjg(U)),matmul(D,U))
!Calculation of new total rotational matrix:
Utot=matmul(Utot,U)
!Calculation of new offdiagonalentries of D:
Doff=D-(Id*D)
!Calculation of the error in each step:
Serr=sum(dconjg(Doff)*Doff) write(*,*) Serr
!Increase of counter cycle by 1: ct=ct+1
if (Serr .lt. err_tol) then write(*,'(A,i3,A)') 'Normal termination after ', ct, ' iterations.' exit endif if (ct.gt.100) then write(*,*) 'Non-convergence of the Jacobi-Sweep algorithm.' exit endif
enddo
deallocate(U,Doff,Id)
if (ct.gt.100) then stop endif
!Resorting of D and U correspondint to increasing eigenvalues:
allocate(lbl(1:N), Dtr(1:N), Utemp(1:N,1:N), Dtemp(1:N,1:N)) lbl=0 Utemp=Utot
do idx=1, N
Dtr(idx)=D(idx,idx)
enddo
do idx=1, N
lbl(idx)=int(minloc(dreal(Dtr),1)) Dtr(lbl(idx))=maxval(dreal(Dtr))+1 Utemp(:,idx)=Utot(:,lbl(idx))
enddo
!Relabelling of eigenvectors: Utot=Utemp deallocate(lbl,Dtr,Dtemp,Utemp)
!(Test) diagonalization:
D=matmul(transpose(dconjg(Utot)),matmul(A,Utot)) deallocate(A)
!Resulting diagonal matrix written:
write(*,*) 'Diagonalized matrix:'
do idx=1, N
write(*,*) D(idx,:)
enddo
!Matrix of eigenvector written:
write(*,*) 'Matrix of eigenvectors:'
do idx=1, N
write(*,*) Utot(idx,:)
enddo
!Test of eigenvectors:
write(*,*) 'adj(U)*U:'
Id=matmul(dconjg(transpose(Utot)),Utot)
do idx=1, N
write(*,*) Id(idx,:)
enddo
!deallocate(Id)
deallocate(D,Utot)
stop end program get_complex_diagonal