Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- eigenvalues and eigenvectors using mkl lapack libraries in Fortran

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

ricoamor

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-06-2011
01:12 AM

354 Views

eigenvalues and eigenvectors using mkl lapack libraries in Fortran

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
```

Link Copied

4 Replies

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-06-2011
02:44 AM

354 Views

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

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

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-06-2011
04:13 AM

354 Views

> 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?

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-06-2011
06:03 AM

354 Views

> 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

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-06-2011
07:14 AM

354 Views

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.