Jump to content

Talk:Jacobi method for complex Hermitian matrices

Page contents not supported in other languages.
From Wikipedia, the free encyclopedia
This is an old revision of this page, as edited by 141.2.134.104 (talk) at 17:42, 17 August 2015. The present address (URL) is a permanent link to this revision, which may differ significantly from the current revision.

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

Retrieved from "https://en.wikipedia.org/w/index.php?title=Talk:Jacobi_method_for_complex_Hermitian_matrices&oldid=676544071"