- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This looks a little bit different from what is in the manual. Maybe an update there is required too.
Alin
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Alin,
You can see the same libraries if run dynamic lapack examples:
cd
make soem64t
-Vladimir
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
-Todd
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page