From 68642e5d21b0cb88b82aa6c7a681e3b448197363 Mon Sep 17 00:00:00 2001 From: William Daniau Date: Tue, 3 Nov 2020 21:09:33 +0100 Subject: [PATCH] Using [SDCZ]GEEVX instead of [SDCZ]GEEV to avoid use of [SDCZ]GEBAL which cause NaN in some circumstances --- fvn_linear/fvn_linear.f90 | 58 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 49 insertions(+), 9 deletions(-) diff --git a/fvn_linear/fvn_linear.f90 b/fvn_linear/fvn_linear.f90 index 8bfa9cf..a07f22b 100644 --- a/fvn_linear/fvn_linear.f90 +++ b/fvn_linear/fvn_linear.f90 @@ -1095,6 +1095,7 @@ subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) ! integer(kind=ip_kind) (out) : status =0 if something went wrong ! ! interfacing Lapack routine SGEEV + ! 2020 11 03 : SGEEVX with BALANC='N' implicit none integer(kind=ip_kind), intent(in) :: d real(kind=sp_kind), intent(in) :: a(d,d) @@ -1114,6 +1115,10 @@ subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) integer(kind=ip_kind) i integer(kind=ip_kind) j + integer(kind=ip_kind) :: ilo,ihi, iwork + real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv + real(kind=sp_kind) :: abnrm + if (present(status)) status=1 ! making a working copy of a @@ -1124,12 +1129,16 @@ subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) allocate(wr(d)) allocate(wi(d)) allocate(vr(d,d)) + allocate(scal(d),rconde(d),rcondv(d)) ! query optimal work size - call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) + ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) + call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info) lwork=int(twork(1)) allocate(work(lwork)) - call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) + ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) + call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info) + if (info /= 0) then if (present(status)) status=0 deallocate(work) @@ -1137,6 +1146,7 @@ subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) deallocate(wi) deallocate(wr) deallocate(wc_a) + deallocate(scal,rconde,rcondv) return end if @@ -1159,6 +1169,7 @@ subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) deallocate(wi) deallocate(wr) deallocate(wc_a) + deallocate(scal,rconde,rcondv) ! sorting if (present(sortval) .and. sortval) then @@ -1176,6 +1187,7 @@ subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) ! integer(kind=ip_kind) (out) : status =0 if something went wrong ! ! interfacing Lapack routine DGEEV + ! 2020 11 03 : DGEEVX with BALANC='N' implicit none integer(kind=ip_kind), intent(in) :: d real(kind=dp_kind), intent(in) :: a(d,d) @@ -1195,6 +1207,10 @@ subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) integer(kind=ip_kind) i integer(kind=ip_kind) j + integer(kind=ip_kind) :: ilo,ihi, iwork + real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv + real(kind=dp_kind) :: abnrm + if (present(status)) status=1 ! making a working copy of a @@ -1205,11 +1221,14 @@ subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) allocate(wr(d)) allocate(wi(d)) allocate(vr(d,d)) + allocate(scal(d),rconde(d),rcondv(d)) ! query optimal work size - call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) + ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) + call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info) lwork=int(twork(1)) allocate(work(lwork)) - call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) + ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) + call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info) if (info /= 0) then if (present(status)) status=0 @@ -1218,6 +1237,7 @@ subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) deallocate(wi) deallocate(wr) deallocate(wc_a) + deallocate(scal,rconde,rcondv) return end if @@ -1241,6 +1261,7 @@ subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) deallocate(wi) deallocate(wr) deallocate(wc_a) + deallocate(scal,rconde,rcondv) ! sorting if (present(sortval) .and. sortval) then @@ -1258,6 +1279,7 @@ subroutine fvn_c_matev(d,a,evala,eveca,status,sortval) ! integer(kind=ip_kind) (out) : status =0 if something went wrong ! ! interfacing Lapack routine CGEEV + ! 2020 11 03 : CGEEVX with BALANC='N' implicit none integer(kind=ip_kind), intent(in) :: d complex(kind=sp_kind), intent(in) :: a(d,d) @@ -1274,6 +1296,10 @@ subroutine fvn_c_matev(d,a,evala,eveca,status,sortval) real(kind=sp_kind), allocatable :: rwork(:) complex(kind=sp_kind) :: vl ! unused but necessary for the call + integer(kind=ip_kind) :: ilo,ihi + real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv + real(kind=sp_kind) :: abnrm + if (present(status)) status=1 ! making a working copy of a @@ -1283,11 +1309,14 @@ subroutine fvn_c_matev(d,a,evala,eveca,status,sortval) ! rwork must be allocated before query allocate(rwork(2*d)) + allocate(scal(d),rconde(d),rcondv(d)) ! query optimal work size - call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) + ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) + call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info) lwork=int(twork(1)) allocate(work(lwork)) - call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) + ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) + call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info) if (info /= 0) then if (present(status)) status=0 @@ -1295,6 +1324,7 @@ subroutine fvn_c_matev(d,a,evala,eveca,status,sortval) deallocate(rwork) deallocate(work) deallocate(wc_a) + deallocate(scal,rconde,rcondv) ! sorting if (present(sortval) .and. sortval) then @@ -1312,6 +1342,7 @@ subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) ! integer(kind=ip_kind) (out) : status =0 if something went wrong ! ! interfacing Lapack routine ZGEEV + ! 2020 11 03 : ZGEEVX with BALANC='N' implicit none integer(kind=ip_kind), intent(in) :: d complex(kind=dp_kind), intent(in) :: a(d,d) @@ -1328,6 +1359,11 @@ subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) real(kind=dp_kind), allocatable :: rwork(:) complex(kind=dp_kind) :: vl ! unused but necessary for the call + integer(kind=ip_kind) :: ilo,ihi + real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv + real(kind=dp_kind) :: abnrm + + if (present(status)) status=1 ! making a working copy of a @@ -1337,11 +1373,14 @@ subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) ! rwork must be allocated before query allocate(rwork(2*d)) + allocate(scal(d),rconde(d),rcondv(d)) ! query optimal work size - call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) + ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) + call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info) lwork=int(twork(1)) allocate(work(lwork)) - call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) + ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) + call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info) if (info /= 0) then if (present(status)) status=0 @@ -1349,6 +1388,7 @@ subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) deallocate(rwork) deallocate(work) deallocate(wc_a) + deallocate(scal,rconde,rcondv) ! sorting if (present(sortval) .and. sortval) then @@ -1600,4 +1640,4 @@ deallocate(mat,bmat,singval) end subroutine -end module fvn_linear \ No newline at end of file +end module fvn_linear -- 2.16.4