eigenvalues and eigenvectors using mkl lapack libraries in Fortran

[ricoamor]

I am trying to calculate the eigenvalues and eigenvectors of matrices of different sizes. I am using a piece of very simple Fortran90 code and I am compiling it linking to the appropriate Lapack libraries included in the Intel MKL package, available in my machine, which runs in Ubuntu. The code "matrix_diag_01.f90" is attached at the end of the message. The "random" module just includes the "ran" random number generator from Numerical Recipes. The code compiles well using

ifort -I $(MKLPATH) -o matrix_diag_01 matrix_diag_01.f90 
      random.f90 $(MKLPATH)libmkl_lapack95.a -Wl,--start-group  
      $(MKLPATH)libmkl_intel_lp64.a $(MKLPATH)libmkl_lapack.a 
      $(MKLPATH)libmkl_intel_thread.a $(MKLPATH)libmkl_core.a
       -Wl,--end-group -lguide -lpthread

The executable works nicely when smallish matrices are given. However, for matrices of size 3000x3000 it produces some strange behaviour. First it gives this error

MKL ERROR : Parameter 8 was incorrect on entry to SSYEVD

However, there are only 3 parameters in the call to SSYEVD. Second, it returns the eigenvectors but not the eigenvalues. I have checked by compiling in another machine with larger memory but the outcome was the same.

Can anyone please help?

Thanks!

PROGRAM matrix_diag_01

USE random

IMPLICIT NONE

INTERFACE 
   SUBROUTINE diag(mat,n)
      INTEGER n
      REAL,DIMENSION(n,n) :: mat
   END SUBROUTINE
END INTERFACE 

INTEGER n,i,j,iseed
REAL, DIMENSION(:), ALLOCATABLE  :: w
REAL, DIMENSION(:,:), ALLOCATABLE :: mat

write (*,*) ' Please enter size of matrix'
read (*,*) n
write (*,*) ' Please type seed'
read (*,*) iseed

allocate (mat(n,n))

do i = 1,n
   do j = 1,n
      mat(i,j) = ran(iseed)
   end do
end do

call diag(mat,n)

stop
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE diag(mat,n)

  USE mkl95_lapack
  USE mkl95_precision

  IMPLICIT NONE

  CHARACTER(len=1) :: jobz = 'V'

  INTEGER n,i

  REAL,DIMENSION(n,n) :: mat
  REAL,DIMENSION(:,:),ALLOCATABLE :: matt,a
  REAL,DIMENSION(:),ALLOCATABLE :: w

  allocate (matt(n,n),a(n,n),w(n))

  matt = mat*transpose(mat)
  a = sqrt(matt)
  open (unit=7,file="matrix.dat",status="unknown")
  do i = 1,n
     write (7,100) a(i,:)
  end do
  close (unit=7)

  call syevd(a,w,jobz)

  open (unit=8,file="eig_val.dat",status="unknown")
  do i = 1,n
     write (8,100) w(i)
  end do
  close (unit=8)

  open (unit=9,file="eig_vec.dat",status="unknown")
  do i = 1,n
     write (9,100) a(i,:)
  end do
  close (unit=9)

  return
 100 format(5000f16.5)
end

[mecej4]

>First it gives this error

>MKL ERROR : Parameter 8 was incorrect on entry to SSYEVD

>However, there are only 3 parameters in the call to SSYEVD. Second, it >returns the eigenvectors but not the eigenvalues.

You call the Fortran 95 module procedure SYEVD, which is disambiguated and after suitable preparation the Fortran 77 subroutine SSYEVD gets called with its 11 arguments. The error message is from the latter. Seeing such a message is usually a symptom of something having gone wrong earlier.

When an MKL routine aborts prematurely, the output arrays may contain whatever junk values were in them, or some intermediate results. Thus, I don't see how you can say that it returned eigenvectors but not eigenvalues.

An arbitrary symmetric matrix need not have all real distinct eigenvalues.

I removed the file writes (I have no interest in the eigenvectors of random matrices), used a call to RANDOM_NUMBER instead of RAN() and the program ran with no error messages with n=3000 even with the 32-bit Intel compiler.

You should consider checking for status with your array allocations and library routine calls, rather than assuming success.

[ricoamor]

> When an MKL routine aborts prematurely, the output arrays may contain whatever

> junk values were in them, or some intermediate results. Thus, I don't see how you

> can say that it returned eigenvectors but not eigenvalues.

Well, when I wrote the eigenvalues for matrices smaller as small as 2*2 and as large as 2500x2500 eigenvectors and eigenvalues were all written and real. For matrices like 3000x3000 it suddenly starts writing all zeros for the eigenvalues while it outputs real (and positive?) eigenvectors.

>An arbitrary symmetric matrix need not have all real eigenvalues.
Oh, that is surprising! I thought that any symmetric matrix was hermitian and therefore had real eigenvalues. Would you propose using a method that allows for complex eigenvalues?

>I removed the file writes (I have no interest in the eigenvectors of random matrices), used a

> call to RANDOM_NUMBER instead of RAN() and the program ran with no error messages with

> n=3000 even with the 32-bit Intel compiler.

For me the program also works. The behaviour however remains unchanged by using one or other random number generator. Also, I am interested in the result since I will of course use the working piece of code for the analysis of non-random matrices. And in order to check the correctness of the result I write them and compare them with another program that can calculate the same quantity albeit slowly (i.e. MATLAB). The comparison is perfect but for these large matrices I am trying to use, where I get the zero eigenvalues.

>You should consider checking for status with your array allocations and library routine calls, rather

> than assuming success.

I am actually assuming failure, not success :-). Would you mind giving pointers on how to do the things you suggest?

Thanks

[mecej4]

> >An arbitrary symmetric matrix need not have all real eigenvalues.
> Oh, that is surprising!

Sorry, my slip; I have made a correction in #1.

Most of the Lapack routines return an optional 'info' argument. If it is not zero, something went wrong.

Here is a modified version of your program that lets you check the eigenvalues and eigenvectors in a more compact way, without having to write out files and reading them into Matlab. After the call to SYEVD returns, for each eigenvalue-eigenvector pair (, x) it computes ||A x - x|| / ||x||, and prints the largest such value. Remember that you should not expect much precision when calculating eigenvalues in single-precision for large matrices. For n = 4000 I get max(||A x - x|| / ||x||) = 1.83E-03. You could try with the other RNG in place of the built-in one.

[fortran]PROGRAM matrix_diag_01

IMPLICIT NONE

INTERFACE 
   SUBROUTINE diag(mat,n)
      INTEGER n
      REAL,DIMENSION(n,n) :: mat
   END SUBROUTINE
END INTERFACE 

INTEGER n,i,j
REAL, DIMENSION(:), ALLOCATABLE  :: w
REAL, DIMENSION(:,:), ALLOCATABLE :: mat

write (*,'(A,$)') ' Please enter size of matrix : '
read (*,*) n

allocate (mat(n,n))

do j = 1,n
   do i = 1,n
      call random_number(mat(i,j))
   end do
end do

call diag(mat,n)

stop
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE diag(mat,n)

  USE mkl95_lapack
  USE mkl95_precision
  USE mkl95_blas
  IMPLICIT NONE

  CHARACTER(len=1) :: jobz = 'V'

  INTEGER n,i

  REAL,DIMENSION(n,n) :: mat
  REAL,DIMENSION(:,:),ALLOCATABLE :: matt,a
  REAL,DIMENSION(:),ALLOCATABLE :: w,y
  integer :: info
  real :: x

  allocate (matt(n,n),a(n,n),w(n),y(n))

  matt = sqrt(mat*transpose(mat))
  a = matt       !save copy for checking later

  call syevd(a,w,jobz,info=info)
  if(info.ne.0)then
     write(*,*)' Info = ',info,' from SYEVD'
     stop
  endif
  
  x=0
  do i=1,n
     y=a(:,i)
     call gemv(matt,a(:,i),y,1.0,-w(i))
     x=max(x,nrm2(y)/nrm2(a(:,i)))
  end do
  write(*,*)' max(Ax-lambda x) = ',x

  return
end
[/fortran]

[ricoamor]

Thanks again, mecej4.

Apparently there is a problem with a wrapper that is used in my system version of mkl so one needs to use the fortran 77 version. With this the program seems to work better.