From 6dd35851650ad66ba50c59e7125ef43b8f3cae45 Mon Sep 17 00:00:00 2001 From: William Daniau Date: Sat, 17 Oct 2009 00:18:14 +0200 Subject: [PATCH] =?UTF-8?q?Ajout=20de=20fvn=5Fd=5Fmatevx=20utilisant=20la?= =?UTF-8?q?=20routine=20lapack=20dgeevx=20+=20programme=20de=20test=20pour?= =?UTF-8?q?=20matev=20s=C3=A9par=C3=A9=20en=20un=20test=20r=C3=A9el=20et?= =?UTF-8?q?=20un=20test=20complexe?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- fvn_linear/fvn_linear.f90 | 98 +++++++++++++++++++++++++++ fvn_test/{test_matev.f90 => test_matev_c.f90} | 8 +-- fvn_test/test_matev_r.f90 | 55 +++++++++++++++ 3 files changed, 157 insertions(+), 4 deletions(-) rename fvn_test/{test_matev.f90 => test_matev_c.f90} (82%) create mode 100644 fvn_test/test_matev_r.f90 diff --git a/fvn_linear/fvn_linear.f90 b/fvn_linear/fvn_linear.f90 index 52a8acf..e645c27 100644 --- a/fvn_linear/fvn_linear.f90 +++ b/fvn_linear/fvn_linear.f90 @@ -1356,6 +1356,104 @@ subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) end subroutine +! This version use dgeevx +subroutine fvn_d_matevx(d,a,evala,eveca,status,sortval) + ! + ! integer d (in) : matrice rank + ! double precision a(d,d) (in) : The Matrix + ! double complex evala(d) (out) : eigenvalues + ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector + ! integer (out) : status =0 if something went wrong + ! + ! interfacing Lapack routine DGEEVX + implicit none + integer, intent(in) :: d + double precision, intent(in) :: a(d,d) + double complex, intent(out) :: evala(d) + double complex, intent(out) :: eveca(d,d) + integer, intent(out), optional :: status + logical, intent(in), optional :: sortval + + double precision, allocatable :: wc_a(:,:) ! a working copy of a + integer :: info + integer :: lwork + double precision, allocatable :: wr(:),wi(:) + double precision, allocatable :: vl(:,:) + double precision, allocatable :: vr(:,:) + double precision, allocatable :: work(:) + double precision :: twork(1) + integer i + integer j + + integer :: ilo,ihi + double precision, allocatable, dimension(:) :: scal,rconde,rcondv + double precision :: abnrm + integer, allocatable, dimension(:) :: iwork + + if (present(status)) status=1 + + ! making a working copy of a + allocate(wc_a(d,d)) + !call dcopy(d*d,a,1,wc_a,1) + wc_a(:,:)=a(:,:) + + allocate(wr(d)) + allocate(wi(d)) + allocate(vr(d,d),vl(d,d)) + allocate(scal(d),rconde(d),rcondv(d),iwork(2*d-2)) + ! query optimal work size + call dgeevx('B','V','V','B',d,wc_a,d,wr,wi,vl,d,vr,d,ilo,ihi,scal,abnrm,rconde,rcondv,twork,-1,iwork,info) + lwork=int(twork(1)) + allocate(work(lwork)) + call dgeevx('B','V','V','B',d,wc_a,d,wr,wi,vl,d,vr,d,ilo,ihi,scal,abnrm,rconde,rcondv,work,lwork,iwork,info) + + if (info /= 0) then + if (present(status)) status=0 + deallocate(work) + deallocate(iwork) + deallocate(rcondv) + deallocate(rconde) + deallocate(scal) + deallocate(vr) + deallocate(vl) + deallocate(wi) + deallocate(wr) + deallocate(wc_a) + return + end if + + ! now fill in the results + i=1 + do while(i<=d) + evala(i)=dcmplx(wr(i),wi(i)) + if (wi(i) == 0.) then ! eigenvalue is real + eveca(:,i)=dcmplx(vr(:,i),0.) + else ! eigenvalue is complex + evala(i+1)=dcmplx(wr(i+1),wi(i+1)) + eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) + eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) + i=i+1 + end if + i=i+1 + enddo + + deallocate(work) + deallocate(iwork) + deallocate(rcondv) + deallocate(rconde) + deallocate(scal) + deallocate(vr) + deallocate(vl) + deallocate(wi) + deallocate(wr) + deallocate(wc_a) + + ! sorting + if (present(sortval) .and. sortval) then + call fvn_z_sort_eigen(d,evala,eveca) + end if +end subroutine + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Least square problem diff --git a/fvn_test/test_matev.f90 b/fvn_test/test_matev_c.f90 similarity index 82% rename from fvn_test/test_matev.f90 rename to fvn_test/test_matev_c.f90 index 28cd173..70a6705 100644 --- a/fvn_test/test_matev.f90 +++ b/fvn_test/test_matev_c.f90 @@ -1,10 +1,10 @@ program test_matev use fvn_linear implicit none -complex(kind=dp_kind),dimension(3,3) :: a -real(kind=dp_kind),dimension(3,3) :: ra,ia -complex(kind=dp_kind),dimension(3) :: evala -complex(kind=dp_kind),dimension(3,3) :: eveca +complex(8),dimension(3,3) :: a +real(8),dimension(3,3) :: ra,ia +complex(8),dimension(3) :: evala +complex(8),dimension(3,3) :: eveca integer :: status,i,j call init_random_seed() diff --git a/fvn_test/test_matev_r.f90 b/fvn_test/test_matev_r.f90 new file mode 100644 index 0000000..e626695 --- /dev/null +++ b/fvn_test/test_matev_r.f90 @@ -0,0 +1,55 @@ +program test_matev +use fvn_linear +implicit none +real(8),dimension(3,3) :: a +complex(8),dimension(3) :: evala +complex(8),dimension(3,3) :: eveca +integer :: status,i,j + +call init_random_seed() +call random_number(a) +a=a*100 +call fvn_matev(3,a,evala,eveca,status) + +write(*,*) "The matrix :" +write (*,'(3e12.5)') a +write (*,*) +do i=1,3 + write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i) + write(*,'("Modulus : ",e12.5)') abs(evala(i)) + write(*,*) "Associated Eigenvector :" + do j=1,3 + write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i) + end do + write(*,*) +end do + +write(*,*) "Computation using matevx" +call fvn_d_matevx(3,a,evala,eveca,status) +do i=1,3 + write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i) + write(*,'("Modulus : ",e12.5)') abs(evala(i)) + write(*,*) "Associated Eigenvector :" + do j=1,3 + write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i) + end do + write(*,*) +end do + + + +! tri +write(*,*) "With sort option" +call fvn_matev(3,a,evala,eveca,status,.true.) +do i=1,3 + write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i) + write(*,'("Modulus : ",e12.5)') abs(evala(i)) + write(*,*) "Associated Eigenvector :" + do j=1,3 + write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i) + end do + write(*,*) +end do + + +end program -- 2.16.4