heevr and linux mkl10.0.3.020

[Alin_M_Elena]

Hi

I try to move away from using heev to compute eigenvalues and eigenvectors of a hermitian matrix towards heevr.

I call heev like this

call heev(c%a,lambda,'V','U',info)
info integer
lambda(1:n) real(k_pr)
c%a(1:n,1:n) complex(k_pr)
k_pr=kind(1.0d0)

and I link like this

CMakeFiles/tdtbuj.dir/source/main.f90.o -o bin/tdtbuj -i_dynamic -L/opt/intel/mkl/10.0.3.020/lib/em64t lib/libparser.a lib/liblinalg.a /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_lapack95.a /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_solver_lp64.a /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_blas95.a -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core /opt/intel/mkl/10.0.3.020/lib/em64t/libguide.so -lpthread -Wl,-rpath,/opt/intel/mkl/10.0.3.020/lib/em64t

now I replace heev with heevr

call heevr(c%a,lambda,'U',info)

and I get this nice error
fortcom: Error: /home/alin/tdtbuj/source/LinearAlgebra.F90, line 233: There is no matching specific subroutine for this generic subroutine call. [HEEVR]
call heevr(c%a,lambda,'U',info)

I use the lapack95 module to access this functions
use mkl95_LAPACK, only : heevr
or
use mkl95_LAPACK, only : heev

This makes me to think that I miss a parameter. To be even funnier it seems that the documentation of heevr is a little bit inconsistent.

Fortran 95:

call heevr(a, w [,uplo] [,z] [,vl] [,vu] [,il] [,iu] [,m] [,isuppz] [,abstol] [,info])

in the interface note info is missing and two new parameters are added jobz, range

http://www.intel.com/software/products/mkl/docs/WebHelp/lse/functn_heevr.html

I get rid of the info parameter and it "works".
I get the same nice eigenvalues but different eigenvectors.

Here is the real part of the input matrix, he complex one is identical zero.


-1.423200 0.000000 0.000000 0.000000 -0.500239 -0.500239
&nb sp; 0.000000 -0.813300 0.000000 0.000000 0.026084 0.026084
0.000000 0.000000 -0.813300 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 -0.813300 0.521674 -0.521674
-0.500239 0.026084 0.000000 0.521674 -1.000000 0.000000
-0.500239 0.026084 0.000000 -0.521674 0.000000 -1.000000

the heev call eigenvectors are
0.801665 0.000000 -0.066944 0.000000 0.594013 0.000000
-0.019381 0.000000 -0.996098 0.000000 -0.086101 0.000000
0.000000 0.000000 0.000000 1.000000 0.000000 0.000000
0.000000 -0.661237 0.000000 0.000000 0.000000 -0.750177
0.422467 0.530455 0.040667 0.000000 -0.565568 -0.467565
0.422467 -0.530455 0.040667 0.000000 -0.565568 0.467565
the heevr call eigenvectors are
-0.813300 0.000000 0.005807 -0.999933 -0.371456 -0.691673
0.000000 -0.814035 -0.054115 0.011613 -0.021539 0.036065
0.000000 0.000000 -1.138045 -0.720062 0.000000 0.000000
0.000000 0.000000 0.000000 -0.992637 -0.306076 -0.721310
  ; -0.500239 0.026084 0.000000 0.521674 -1.105083 -0.723231
-0.500239 0.026084 0.000000 -0.521674 0.000000 -1.000000

both of them give the same eigenvalues
-1.95043932
-1.65029017
-0.81542979
-0.81330000
-0.47063089
-0.16300983


Any reasons for this behaviour?

All these happen on a linux box
Linux green 2.6.25.5-1.1-default #1 SMP 2008-06-07 01:55:22 +0200 x86_64 x86_64x86_64 GNU/Linux
with opensuse 11
the fortran compiler is
ifort --version
ifort (IFORT) 10.1 20080602
Copyright (C) 1985-2008 Intel Corporation. All rights reserved.
and the mkl 10.0.3.020
on a Intel Core2 CPU T5300@1.73GHZ with 2GB of RAM


Alin

[Vladimir_Koldakov__I]

Hi, Alin,

Square brackets in Fortran 95 call syntax represent optional arguments. It is better to use argument keywords for optional arguments. Below is the fragment of the Intel Fortran Compiler documentation:

To call a procedure with an optional argument, you must use an explicit interface. If argument keywords are not used, argument association is positional. The first dummy argument becomes associated with the first actual argument, and so on. If argument keywords are used, arguments are associated by the keyword name, so actual arguments can be in a different order than dummy arguments. A keyword is required for an argument only if a preceding optional argument is omitted or if the argument sequence is changed.

So, if you call heev as

call heev(c%a,lambda,'V','U',info)

its ok, because all arguments are on place.

Fortran 95 syntax:

call heev(a, w [,jobz] [,uplo] [,info])

But if you call heevr as:

call heevr(c%a,lambda,'U',info)

it is wrong, because in this case info argument comes as z argument which is matrix.

Fortran 95 syntax:

call heevr(a, w [,uplo] [,z] [,vl] [,vu] [,il] [,iu] [,m] [,isuppz] [,abstol] [,info])

So, replace your call statement with the next one:

call heevr(c%a,lambda,'U',info=info)

or better:

call heevr(c%a,lambda,uplo='U',info=info)

Thanks,

Vladimir

[Alin_M_Elena]

Thank you Vladimir,

Of course the arguments were wrong. Saturday is not a good day to work you can do silly mistakes.

The problem that puzzles me is the computational aspect and why I get different results.

Alin

[Vladimir_Koldakov__I]

Alin,

Try to add one more argument for heevr call: call heevr(c%a,lambda,z=z,uplo='U',info=info) where z is the matrix.

According to documentation for ?heev, matrix a on exit contains eigenvectors (if jobz=v):

On exit, if jobz = 'V', then if info = 0, array a contains the orthonormal eigenvectors of the matrix A.

If jobz = 'N', then on exit the lower triangle (if uplo = 'L') or the upper triangle (if uplo = 'U') of A, including the diagonal, is overwritten.

On the other hand, for ?heevr matrix a is always broken:

On exit, the lower triangle (if uplo = 'L') or the upper triangle (if uplo = 'U') of A, including the diagonal, is overwritten.

To get eigenvectors with heevryou should use z argument:

COMPLEX for cheevr

DOUBLE COMPLEX for zheevr.

Array z(ldz, *); the second dimension of z must be at least max(1, m).

If jobz = 'V', then if info = 0, the first m columns of z contain the orthonormal eigenvectors of the matrix T corresponding to the selected eigenvalues, with the i-th column of z holding the eigenvector associated with w(i).

If jobz = 'N', then z is not referenced.

Note: you must ensure that at least max(1,m) columns are supplied in the array z; if range = 'V', t he exact value of m is not known in advance and an upper bound must be used.

Thanks,

Vladimir

[Alin_M_Elena]

Thank you, Vladimir!

Still no luck now I get a segmentation fault.

Here is a quick example

program testHeevr
use mkl95_LAPACK, only : heevr
implicit none
integer, parameter :: pr=kind(1.0d0)
complex(pr) :: a(6,6),b(6,6),c(6,6)
real(pr) :: l(6)
integer :: i,j,info
c=cmplx(0.0,0.0,pr)
a(1,:)=(/-1.423200,0.000000,0.000000,0.000000,-0.500239,-0.500239/)
a(2,:)=(/0.000000, -0.813300, 0.000000, 0.000000, 0.026084, 0.026084/)
a(3,:)=(/0.000000, 0.000000, -0.813300, 0.000000, 0.000000, 0.000000/)
a(4,:)=(/0.000000, 0.000000, 0.000000, -0.813300, 0.521674, -0.521674/)
a(5,:)=(/-0.500239, 0.026084, 0.000000, 0.521674, -1.000000, 0.000000/)
a(6,:)=(/-0.500239, 0.026084, 0.000000, -0.521674, 0.000000, -1.000000/)
do i=1,6
write(*,'(6f10.4)') (real(a(i,j)),j=1,6)
enddo
do i=1,6
write(*,'(6f10.4)') (aimag(a(i,j)),j=1,6)
enddo
b=a
c=a
l=0.0_pr
info=-1
call heevr(a=b,w=l,uplo="U",z=c,info=info)
do i=1,6
write(*,'(7f10.4)') l(i),(real(c(i,j)),j=1,6)
enddo
end program testHeevr

this is what I used to compile

ifort -o a test.f90 -i_dynamic-L/opt/intel/mkl/10.0.3.020/lib/em64t -lmkl_lapack95 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core /opt/intel/mkl/10.0.3.020/lib/em64t/libguide.so -lpthread -Wl,-rpath,/opt/intel/mkl/10.0.3.020/lib/em64t

Alin

[Vladimir_Koldakov__I]

Alin,

Really, I got a Segmentation fault with your compilation/linking options.

I think this is linking problem. Try to replace the

-lmkl_intel_thread -lmkl_core

options in your command line with the next

-Wl,--start-group /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_intel_thread.a /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_core.a -Wl,--end-group

You can see the command line with similar options at the time of examples run.

Thanks,

Vladimir

[Alin_M_Elena]

Thanks Vladimir,

Sorry for the delayed answer. I have been away for few days.

The solution that you suggest works perfectly but it changes the linkage model from dynamic to static. Should I conclude that I can not use dynamic linkage and fortran 95 interfaces?

I prefer dynamic linkage as my problem is not very big.

Alin

[Vladimir_Koldakov__I]

Alin,

Really, the size is too large.

$ ifort -o a -Ilib/em64t t.f90 -i_dynamic -Llib/em64t -lmkl_lapack95 -L../../lib/em64t -lmkl_intel_lp64 -Wl,--start-group ../../lib/em64t/libmkl_intel_thread.a ../../lib/em64t/libmkl_core.a -Wl,--end-group -lguide -lpthread -Wl,-rpath,../../lib/em64t

t.f90(23): (col. 4) remark: LOOP WAS VECTORIZED.

$ wc -c a

49039385 a

Try the next instead of previous:

-Wl,--start-group /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_intel_thread.so /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_lapack.so /opt/intel/mkl/10.0.3.020/lib/em64t/libmkl_core.so -Wl,--end-group

$ ifort -o a -Ilib/em64t t.f90 -i_dynamic -Llib/em64t -lmkl_lapack95 -L../../lib/em64t -lmkl_intel_lp64 -Wl,--start-group ../../lib/em64t/libmkl_intel_thread.so ../../lib/em64t/libmkl_lapack.so ../../lib/em64t/libmkl_core.so -Wl,--end-group -lguide -lpthread -Wl,-rpath,../../lib/em64t

t.f90(23): (col. 4) remark: LOOP WAS VECTORIZED.

$ wc -c a

30681 a

$ ./a

-1.4232 0.0000 0.0000 0.0000 -0.5002 -0.5002

0.0000 -0.8133 0.0000 0.0000 0.0261 0.0261

0.0000 0.0000 -0.8133 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 -0.8133 0.5217 -0.5217

-0.5002 0.0261 0.0000 0.5217 -1.0000 0.0000

-0.5002 0.0261 0.0000 -0.5217 0.0000 -1.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

-1.9504 0.8017 0.0000& nbsp; -0.0669 0.0000 -0.5940 0.0000

-1.6503 -0.0194 0.0000 -0.9961 0.0000 0.0861 0.0000

-0.8154 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000

-0.8133 0.0000 -0.6612 0.0000 0.0000 0.0000 -0.7502

-0.4706 0.4225 0.5305 0.0407 0.0000 0.5656 -0.4676

-0.1630 0.4225 -0.5305 0.0407 0.0000 0.5656 0.4676

-Vladimir

[Alin_M_Elena]

Thanks Vladimir,

This looks a little bit different from what is in the manual. Maybe an update there is required too.

Alin

[Alin_M_Elena]

Having a second reading of your post I am little bit confused.

I do not understand why do you mix the layered default with layered pure?
mkl_lapack is from layered default
mkl_intel_thread from layered pure?

Alin

[Vladimir_Koldakov__I]

Alin,

You can see the same libraries if run dynamic lapack examples:

cd /examples/lapack

make soem64t

-Vladimir

[Alin_M_Elena]

Thanks Vladimir,

That means that the 5-5 -- 5-8 part of the userguide is a little bit confusing stating other things.


LAPACK, static case:
Layered default: libmkl_lapack.a libmkl_em64t.a
Layered pure: libmkl_intel_lp64.a libmkl_intel_thread.a libmkl_core.a

LAPACK, dynamic case:
Layered default:libmkl_lapack.so libmkl.so
Layered pure: libmkl_intel_lp64.so libmkl_intel_thread.so libmkl_core.so

I need a reliable source of information about linking of mkl as I am a maintainer of BLAS/LAPACk modules for cmake. I added some adhoc support for intel mkl libs and know I intend to do it in a proper way.

Alin

[Todd_R_Intel]

I agree; 5-5 -- 5-8 in the User Guide is confusing on this point. I'll submit a change request.
-Todd