Ajout de fvn_d_matevx utilisant la routine lapack dgeevx

+ programme de test pour matev séparé en un test réel et un test complexe

Ajout de fvn_d_matevx utilisant la routine lapack dgeevx
+ programme de test pour matev séparé en un test réel et un test complexe
William Daniau
1 parent 41811905dd
Showing 4 changed files with 157 additions and 4 deletions Side-by-side Diff
fvn_linear/fvn_linear.f90
fvn_test/test_matev.f90
fvn_test/test_matev_c.f90
fvn_test/test_matev_r.f90
@@ -1356,6 +1356,104 @@
 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
-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
-integer :: status,i,j
-
-call init_random_seed()
-call random_number(ra)
-call random_number(ia)
-a=ra+fvn_i*ia
-a=a*100
-call fvn_matev(3,a,evala,eveca,status)
-
-write(*,*) "The matrix :"
-write (*,'(3("(",e12.5,",",e12.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
-
-! 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
+program test_matev
+use fvn_linear
+implicit none
+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()
+call random_number(ra)
+call random_number(ia)
+a=ra+fvn_i*ia
+a=a*100
+call fvn_matev(3,a,evala,eveca,status)
+
+write(*,*) "The matrix :"
+write (*,'(3("(",e12.5,",",e12.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
+
+! 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
+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
...	...	@@ -1356,6 +1356,104 @@
1356	1356	end subroutine
1357	1357
1358	1358
	1359	+! This version use dgeevx
	1360	+subroutine fvn_d_matevx(d,a,evala,eveca,status,sortval)
	1361	+ !
	1362	+ ! integer d (in) : matrice rank
	1363	+ ! double precision a(d,d) (in) : The Matrix
	1364	+ ! double complex evala(d) (out) : eigenvalues
	1365	+ ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	1366	+ ! integer (out) : status =0 if something went wrong
	1367	+ !
	1368	+ ! interfacing Lapack routine DGEEVX
	1369	+ implicit none
	1370	+ integer, intent(in) :: d
	1371	+ double precision, intent(in) :: a(d,d)
	1372	+ double complex, intent(out) :: evala(d)
	1373	+ double complex, intent(out) :: eveca(d,d)
	1374	+ integer, intent(out), optional :: status
	1375	+ logical, intent(in), optional :: sortval
	1376	+
	1377	+ double precision, allocatable :: wc_a(:,:) ! a working copy of a
	1378	+ integer :: info
	1379	+ integer :: lwork
	1380	+ double precision, allocatable :: wr(:),wi(:)
	1381	+ double precision, allocatable :: vl(:,:)
	1382	+ double precision, allocatable :: vr(:,:)
	1383	+ double precision, allocatable :: work(:)
	1384	+ double precision :: twork(1)
	1385	+ integer i
	1386	+ integer j
	1387	+
	1388	+ integer :: ilo,ihi
	1389	+ double precision, allocatable, dimension(:) :: scal,rconde,rcondv
	1390	+ double precision :: abnrm
	1391	+ integer, allocatable, dimension(:) :: iwork
	1392	+
	1393	+ if (present(status)) status=1
	1394	+
	1395	+ ! making a working copy of a
	1396	+ allocate(wc_a(d,d))
	1397	+ !call dcopy(d*d,a,1,wc_a,1)
	1398	+ wc_a(:,:)=a(:,:)
	1399	+
	1400	+ allocate(wr(d))
	1401	+ allocate(wi(d))
	1402	+ allocate(vr(d,d),vl(d,d))
	1403	+ allocate(scal(d),rconde(d),rcondv(d),iwork(2*d-2))
	1404	+ ! query optimal work size
	1405	+ 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)
	1406	+ lwork=int(twork(1))
	1407	+ allocate(work(lwork))
	1408	+ 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)
	1409	+
	1410	+ if (info /= 0) then
	1411	+ if (present(status)) status=0
	1412	+ deallocate(work)
	1413	+ deallocate(iwork)
	1414	+ deallocate(rcondv)
	1415	+ deallocate(rconde)
	1416	+ deallocate(scal)
	1417	+ deallocate(vr)
	1418	+ deallocate(vl)
	1419	+ deallocate(wi)
	1420	+ deallocate(wr)
	1421	+ deallocate(wc_a)
	1422	+ return
	1423	+ end if
	1424	+
	1425	+ ! now fill in the results
	1426	+ i=1
	1427	+ do while(i<=d)
	1428	+ evala(i)=dcmplx(wr(i),wi(i))
	1429	+ if (wi(i) == 0.) then ! eigenvalue is real
	1430	+ eveca(:,i)=dcmplx(vr(:,i),0.)
	1431	+ else ! eigenvalue is complex
	1432	+ evala(i+1)=dcmplx(wr(i+1),wi(i+1))
	1433	+ eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
	1434	+ eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
	1435	+ i=i+1
	1436	+ end if
	1437	+ i=i+1
	1438	+ enddo
	1439	+
	1440	+ deallocate(work)
	1441	+ deallocate(iwork)
	1442	+ deallocate(rcondv)
	1443	+ deallocate(rconde)
	1444	+ deallocate(scal)
	1445	+ deallocate(vr)
	1446	+ deallocate(vl)
	1447	+ deallocate(wi)
	1448	+ deallocate(wr)
	1449	+ deallocate(wc_a)
	1450	+
	1451	+ ! sorting
	1452	+ if (present(sortval) .and. sortval) then
	1453	+ call fvn_z_sort_eigen(d,evala,eveca)
	1454	+ end if
	1455	+end subroutine
	1456	+
1359	1457	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1360	1458	!
1361	1459	! Least square problem
1		-program test_matev
2		-use fvn_linear
3		-implicit none
4		-complex(kind=dp_kind),dimension(3,3) :: a
5		-real(kind=dp_kind),dimension(3,3) :: ra,ia
6		-complex(kind=dp_kind),dimension(3) :: evala
7		-complex(kind=dp_kind),dimension(3,3) :: eveca
8		-integer :: status,i,j
9		-
10		-call init_random_seed()
11		-call random_number(ra)
12		-call random_number(ia)
13		-a=ra+fvn_i*ia
14		-a=a*100
15		-call fvn_matev(3,a,evala,eveca,status)
16		-
17		-write(,) "The matrix :"
18		-write (*,'(3("(",e12.5,",",e12.5,")"))') a
19		-write (,)
20		-do i=1,3
21		- write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
22		- write(*,'("Modulus : ",e12.5)') abs(evala(i))
23		- write(,) "Associated Eigenvector :"
24		- do j=1,3
25		- write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
26		- end do
27		- write(,)
28		-end do
29		-
30		-! tri
31		-write(,) "With sort option"
32		-call fvn_matev(3,a,evala,eveca,status,.true.)
33		-do i=1,3
34		- write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
35		- write(*,'("Modulus : ",e12.5)') abs(evala(i))
36		- write(,) "Associated Eigenvector :"
37		- do j=1,3
38		- write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
39		- end do
40		- write(,)
41		-end do
42		-
43		-
44		-end program
	1	+program test_matev
	2	+use fvn_linear
	3	+implicit none
	4	+complex(8),dimension(3,3) :: a
	5	+real(8),dimension(3,3) :: ra,ia
	6	+complex(8),dimension(3) :: evala
	7	+complex(8),dimension(3,3) :: eveca
	8	+integer :: status,i,j
	9	+
	10	+call init_random_seed()
	11	+call random_number(ra)
	12	+call random_number(ia)
	13	+a=ra+fvn_i*ia
	14	+a=a*100
	15	+call fvn_matev(3,a,evala,eveca,status)
	16	+
	17	+write(,) "The matrix :"
	18	+write (*,'(3("(",e12.5,",",e12.5,")"))') a
	19	+write (,)
	20	+do i=1,3
	21	+ write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
	22	+ write(*,'("Modulus : ",e12.5)') abs(evala(i))
	23	+ write(,) "Associated Eigenvector :"
	24	+ do j=1,3
	25	+ write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
	26	+ end do
	27	+ write(,)
	28	+end do
	29	+
	30	+! tri
	31	+write(,) "With sort option"
	32	+call fvn_matev(3,a,evala,eveca,status,.true.)
	33	+do i=1,3
	34	+ write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
	35	+ write(*,'("Modulus : ",e12.5)') abs(evala(i))
	36	+ write(,) "Associated Eigenvector :"
	37	+ do j=1,3
	38	+ write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
	39	+ end do
	40	+ write(,)
	41	+end do
	42	+
	43	+
	44	+end program