diff --git a/fvn_linear/fvn_linear.f90 b/fvn_linear/fvn_linear.f90
index 2e8bc79..31a6272 100644
--- a/fvn_linear/fvn_linear.f90
+++ b/fvn_linear/fvn_linear.f90
@@ -1025,6 +1025,67 @@ end subroutine
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+! August 2009
+! William Daniau
+! Adding sorting of eigenvalues and vectors
+subroutine fvn_z_sort_eigen(d,v,vec)
+! this routine takes in input :
+! v : a vector containing unsorted eigenvalues
+! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
+! 
+! At the end of subroutine the eigenvalues are sorted by decreasing order of 
+! modulus so as the vectors.
+implicit none
+integer :: d,i,j
+complex(8),dimension(d) :: v
+complex(8),dimension(d,d) :: vec
+complex(8) :: cur_v
+complex(8), dimension(d) :: cur_vec 
+
+do i=2,d
+  cur_v=v(i)
+  cur_vec=vec(:,i)
+  j=i-1
+  do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
+    v(j+1)=v(j)
+    vec(:,j+1)=vec(:,j)
+    j=j-1
+  end do
+  v(j+1)=cur_v
+  vec(:,j+1)=cur_vec
+end do
+
+end subroutine
+
+subroutine fvn_c_sort_eigen(d,v,vec)
+! this routine takes in input :
+! v : a vector containing unsorted eigenvalues
+! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
+! 
+! At the end of subroutine the eigenvalues are sorted by decreasing order of 
+! modulus so as the vectors.
+implicit none
+integer :: d,i,j
+complex(4),dimension(d) :: v
+complex(4),dimension(d,d) :: vec
+complex(4) :: cur_v
+complex(4), dimension(d) :: cur_vec 
+
+do i=2,d
+  cur_v=v(i)
+  cur_vec=vec(:,i)
+  j=i-1
+  do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
+    v(j+1)=v(j)
+    vec(:,j+1)=vec(:,j)
+    j=j-1
+  end do
+  v(j+1)=cur_v
+  vec(:,j+1)=cur_vec
+end do
+
+end subroutine
+
 subroutine fvn_s_matev(d,a,evala,eveca,status)
     ! 
     ! integer d (in) : matrice rank
@@ -1098,6 +1159,9 @@ subroutine fvn_s_matev(d,a,evala,eveca,status)
     deallocate(wr)
     deallocate(wc_a)
 
+    ! sorting
+    call fvn_c_sort_eigen(d,evala,eveca)
+    
 end subroutine
 
 subroutine fvn_d_matev(d,a,evala,eveca,status)
@@ -1174,6 +1238,9 @@ subroutine fvn_d_matev(d,a,evala,eveca,status)
     deallocate(wr)
     deallocate(wc_a)
 
+    ! sorting
+    call fvn_z_sort_eigen(d,evala,eveca)
+
 end subroutine
 
 subroutine fvn_c_matev(d,a,evala,eveca,status)
@@ -1222,6 +1289,9 @@ subroutine fvn_c_matev(d,a,evala,eveca,status)
     deallocate(work)
     deallocate(wc_a)
 
+    ! sorting
+    call fvn_c_sort_eigen(d,evala,eveca)
+
 end subroutine
 
 subroutine fvn_z_matev(d,a,evala,eveca,status)
@@ -1269,6 +1339,9 @@ subroutine fvn_z_matev(d,a,evala,eveca,status)
     deallocate(work)
     deallocate(wc_a)
 
+    ! sorting
+    call fvn_z_sort_eigen(d,evala,eveca)
+
 end subroutine