Talk:Jacobi method for complex Hermitian matrices

[Untitled]

This 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:

1. for Hermitian matrices it would be more appropriate to use unitary transformations rather than two rotations.
2. 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