1

I'm tryng to learn to use a lapack subroutine but I got stuck. I hope this is the right forum... In this fortran program I'd like as a test to find the Shur form of the matrix ((0,1)(1,0)) using cgees, but the matrix "m" after the call in which the Shur form should be written, is made only of zeros (other subrotunes like dgemm work and I don't get errors so I don't think is a problem with the linking). Do you know where I do wrong?

.........................................................

program example    
implicit none     
integer :: d,sdim,info    
real(kind=8),allocatable ::rwork(:)    
logical, external :: sel    
complex(kind=8), allocatable :: m(:,:),w(:),vs(:,:),work(:)
logical,allocatable ::bwork(:)
d=2    
!    
allocate(m(d,d))    
allocate(w(d))    
allocate(vs(d,d))    
allocate(work(2*d))    
allocate(rwork(d))    
allocate(bwork(d))    
!    
m(1,1)=(0.d0,0.d0)    
m(1,2)=(1.d0,0.d0)    
m(2,1)=(1.d0,0.d0)    
m(2,2)=(0.d0,0.d0)    
!    
print*,'m_before',m    
!    
print*,'w_before',w        
call cgees('V','S',sel,d,m,d,sdim,w,vs,d,work,-1,rwork,bwork,info)    
print*,'m',m    
print*,'w',w    
print*,'vs',vs    
print*,info        
end program example


logical function sel(x)    
  complex(kind=8) ::x    
  if (dimag(x)>1.d0) then    
     sel = .true.       
  else    
     sel = .false.    
  end if    
  return    
end

.........................................

Paul
  • 12,045
  • 7
  • 56
  • 129
Thomas
  • 41
  • 3
  • I can't really find the error. Maybe my lapack is not well installed ?? – Thomas Dec 23 '13 at 23:35
  • Instead of (1.d0,0.d0), try cmplx(1.d0,0.d0) and so forth for your other matrix assignments. Does this work? – Paul Dec 24 '13 at 15:25
  • I put what you said but the result is the same. – Thomas Dec 25 '13 at 10:54
  • Just for information, I was able to run the example here, http://www.amath.unc.edu/sysadmin/DOC4.0/perflib/perflib_ref/chapter1.doc13.html , that uses the subroutine dgees ... But actually I need to understand what's the problem with cgees since I want to deal with complex matrices... – Thomas Dec 26 '13 at 14:02
  • 1
    Thanks to http://scicomp.stackexchange.com/questions/10402/real-eigenvalues-finding I understood that the problem is that cgees works with single precision arithmetic and not double as I thought!!! This solves the problem :) – Thomas Dec 26 '13 at 14:54
  • 1
    I'm glad to hear that you found the answer to your question. It would be a good idea if you could post your comment as an answer and accept it as an answer. – Paul Dec 26 '13 at 15:02

1 Answers1

3

Reading this answer Real eigenvalues finding I understood that the problem is that cgees works with single precision arithmetic and not double as I thought!!! This solves the problem :) ...

In general what I understood is that lapack subroutines are named so that:

  • subroutines beggining with "d" (like DGEMM, DGETRI, DGEES,...) deal with real double precision numbers;

  • subroutines beggining with "c" (like CGEMM, CGETRI, CGEES,...) deal with complex single precision numbers;

  • subroutines beggining with "z" (like ZGEMM, ZGETRI, ZGEES...) deal with complex double precision numbers;

Thx a lot,

Community
  • 1
Thomas
  • 41
  • 3