MKL LAPACK SBGVX example

[Ever_B_]

SBGVX computes a few eigenvalues/eigenvectors of a real symmetric matrix so it should be very useful for FEA where we want only a few eigenvalues/vectors of a very large matrix. So, I implemented a test example following the documentations but I cannot get it to work.  I get an access violation. 

Does anyone has a working example using this routine?

Here is mine:

program Console1
implicit none
    ! Variables
    ! http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm#GUID-659297EE-D35B-428C-ABBE-1A7EBE2B0F6E.htm#GUID-659297EE-D35B-428C-ABBE-1A7EBE2B0F6E 
    !
    ! N: matrix size
    ! NJ: number of upper diagonals + 1 
    integer, parameter :: N=3, NJ=2, ku=NJ-1, kl=0, ldab=ku+1, iwork=7*N
    !compute eigenvalues with index il to iu included
    integer, parameter :: il=1, iu=2, ldz=N, ldq=N
    integer info, m
    integer, dimension (N) :: ifail
    doubleprecision, dimension (N,N) :: KG, KM
    doubleprecision, dimension (ldab,N) :: KGB, KMB
    doubleprecision, dimension (N) :: W !eigenvalues
    doubleprecision, dimension (N,N) :: q, z
    doubleprecision vl, vu  !not referenced if range="I"
    doubleprecision, parameter :: abstol=0.0D-3  !absolute error tolerance for eigenvalues
    doubleprecision, dimension (iwork) :: work
    character*1, parameter :: jobz="V"  !compute eigenvalues
    character*1, parameter :: range="I" !compute in interval il--iu
    character*1, parameter :: uplo="U"  !upper triangular symmetric

    !local
    integer i,j

    m = iu-il+1
    
    !test matrix
    KG = 0
    KM = 0
    KG(1,1) = 100
    KG(2,2) = 100
    KG(1,2) =  25
    KG(2,1) =  25
    do i=1,N
        KM(i,i) = 1
    enddo
    print *, "KG"
    write(*,"(3G9.3)") ((KG(i,j),j=1,N),i=1,N)

    !copy to band form using 
    !http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm#GUID-659297EE-D35B-428C-ABBE-1A7EBE2B0F6E.htm#GUID-659297EE-D35B-428C-ABBE-1A7EBE2B0F6E
    KGB = 0
    KMB = 0
    do i=1,N
        do j=1,N
            if (i.ge.max(1,j-ku).and.i.le.min(N,j+kl)) KGB(ku+1+i-j,j) = KG(i,j)
            if (i.ge.max(1,j-ku).and.i.le.min(N,j+kl)) KMB(ku+1+i-j,j) = KM(i,j)
        enddo
    enddo
    print *, "KGB"
    write(*,"(3G9.3)") ((KGB(i,j),j=1,N),i=1,ku+1)
    print *, "KGM"
    write(*,"(3G9.3)") ((KMB(i,j),j=1,N),i=1,ku+1)
    
    !call sbgvx(KGB, KMB, w) 
    
!                  1      2     3  4   5   6    7     8    9    10 11   12  13  14  15  16      17 18 19 20   21    22     23     24    25
    call dsbgvx(jobz, range, uplo, N, ku, ku, KGB, ldab, KMB, ldab, q, ldq, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
    print *, "Eigenvalues found=",m
    print *, "Eigenvalues. Info=",info
    write(*,*) W
    pause "Hit ENTER to continue"

end program Console1

[mecej4]

As of today, I find that the IDZ forums are blanking out the text of the first post in many threads. This problem affects even the display of the online MKL manual; I looked up sbgvx, and was given a page with nothing but the name of the routine.

I am, therefore, replying based only on the title of the thread, so my response may be off the mark. Assuming that you want to see an example of using ?sbgvx():

[bash]c:\\Program Files (x86)\\Intel\Composer XE 2013 SP1\\mkl\\examples>unzip -v examples_f95.zip | grep sbgvx
      97  Defl:N       49  50% 04-15-2008 05:28 f0b2a533  lapack95/data/sbgvx.d

    3255  Defl:N     1249  62% 12-31-2013 01:01 d2fa34b4  lapack95/source/sbgvx.f90[/bash]

Previous versions of MKL may have provided a fixed format Fortran example for this routine.

[mecej4]

The blank display problem has gone behind the scenes, so I can see your post now. It appears to me that you started out trying to use the F95 interface, ran into difficulties and switched to the F77 interface. 

Your program works fine with a small fix-up. Replace[fortran]program Console1[/fortran] with [fortran]include 'lapack.f90'

program Console1

use mkl95_lapack[/fortran]In fact, the "include.." line can be left out if the module file mkl95_lapack.mod is already present in the .../mkl/include/<arch> directory, where <arch> is "ia32" or "intel64".

Then, after commenting out the F77 call and activating the F95 call, I compiled using the command [bash]ifort /Qmkl /MD /traceback xsbgvx.f90 mkl_lapack95.lib[/bash] Running the resulting program gave me[bash]...

 Eigenvalues found=           2

 Eigenvalues. Info=           0

  0.000000000000000E+000   75.0000000000000        125.000000000000[/bash]

[Ever_B_]

Thanks I go it to work. In addition to the use statement, I set up VS as follows:

! Project> Properties> Fortran> General> Additional include directories> C:\Program Files (x86)\Intel\Composer XE 2011 SP1\mkl\include\ia32
! Project> Properties> Fortran> Libraries> Use Intel Math Kernel Library> Sequential
! Project> Properties> Fortran> External Procedures> Calling Conventions> Default
! Project> Properties> Linker> General> Enable Incremental Linking> No
! Project> Properties> Linker> Input> Additional dependencies> mkl_lapack95.lib

Thanks a lot.

[mecej4]

Ever B. wrote:

How do I tell VS2010 where to find "mkl95_lapack.mod" ?

Add the path to the directory containing the module file to the list (possibly initially empty) of locations in the "include files" box.

Is the mod already there? or do I have to generate it? 

Look in the .../mkl/include/ia32 (or .../mkl/include/intel64/lp64, etc.) directory. If the file is not already there, after you once compile a source file containing the "include 'lapack.f90'" statement the module will be generated in the working directory. You can then copy or move the module file to the appropriate MKL include directory.

[Ever_B_]

Now I have problems with pbsvx. I can make the F95 interface work but not the F77 interface. What am I doing wrong? 

! Project> Properties> Fortran> General> Additional include directories> C:\Program Files (x86)\Intel\Composer XE 2011 SP1\mkl\include\ia32
! Project> Properties> Fortran> Libraries> Use Intel Math Kernel Library> Sequential
! Project> Properties> Fortran> External Procedures> Calling Conventions> Default
! Project> Properties> Linker> General> Enable Incremental Linking> No
! Project> Properties> Linker> Input> Additional dependencies> mkl_lapack95.lib

program Main1
    ! test MKL symmetric linear algebraic system solver

    use mkl95_lapack        ! F95 explicit interface names are here
    implicit none

    integer, parameter :: neq=3, nhbw=2
    integer i, j, info, kd, nrhs, n
    integer, dimension(neq) :: iwork
    real(8) rcond, ferr(1), berr(1), work(3*neq)
    real(8), dimension(neq,neq) :: A
    real(8), dimension(neq) :: b, x, s
    real(8), dimension(nhbw,neq) :: ab, afb

    n = neq
    kd = nhbw-1
    nrhs = 1
    
    A = 0.0 
    do i = 1, neq 
        b (i) = 1.0
        a (i, i) = 1.d0/i 
    enddo 

    do i = 1, neq
        do j = 1, neq
            !from full to band storage
            if (i.ge.max(1,j-kd).and.i.le.j) ab(kd+1+i-j,j) = a(i,j)    
        enddo
    enddo

    
    ! F95 interface works fine
    !call pbsvx( ab, b, x, info=info )

    ! F77 interface -> Access violation. Why?
    call dpbsvx( "N", "U", neq, kd, nrhs, ab, kd+1, afb, kd+1, "N", s, b, neq, x, neq, rcond, ferr, berr, work, iwork, info )
    
    write(*,*) x
    pause "Hit Enter to continue"
endprogram

[mecej4]

The F77 interfaces are more complicated and error-prone, involving many more arguments than the F95 interfaces. In particular, no interface declarations are provided for the F77 routines, so the compiler is unable to check if the arguments are of the correct type. If you still want to use the former, please read the documentation carefully.

The tenth argument, EQUED, is an output variable when the first argument, FACT, is not equal to 'F'. You are passing a character constant, which leads to the access violation. The solution is to declare a character(len=1) variable, say, "equed", assign 'N' to it if you want to, and pass the variable as the tenth argument.

[Ever_B_]

THANK YOU mecej4 !!! 

I was trying to use F77 interface (and thanks to you, I can use it now) because using Interoperative systems to call a Fortran DLL from C#

[DLLImport(Path)]
public static extern void foo()

I could not get the linker to find the F95 MKL's. But (again) thanks to your previous suggestion, reproduced below:

! Project> Properties> Fortran> General> Additional include directories> C:\Program Files (x86)\Intel\Composer XE 2011 SP1\mkl\include\ia32

now, the Linker finds the F95 and my problems are OVER !!!

The DLLImport requires "foo" to have the following !DEC$ declaration:

!DEC$ATTRIBUTES REFERENCE, STDCALL, DLLEXPORT:: foo

My remaining question is this: "I do not understand why...???" but you cannot set everything as STDCALL, REFERENCE in Project> Properties> Fortran> External Procedures> Calling Conventions> STDCALL, REFERENCE,    because then it will not find the MKL's.

So what I did was to let Calling Conventions>DEFAULT for the entire solution and force STDCALL only for "foo" using the !DEC$ card. 

A mystery but it works!

[mecej4 · Accepted solution]

Ever B. wrote:
My remaining question is this: "I do not understand why...???" but you cannot set everything as STDCALL, REFERENCE in Project> Properties> Fortran> External Procedures> Calling Conventions> STDCALL, REFERENCE,    because then it will not find the MKL's.

There is not much of a mystery here. With current releases of IFort/IntelC++/MKL, the STDCALL-compatible runtime libraries are not installed unless you specify at installation time that such libraries should be included. If STDCALL library routines have not been installed and you still specify STDCALL in the VS IDE, you will have link-time failures. I think that name decoration being done differently for STDCALL is a good thing because, without this distinction, we would suffer mysterious run-time failures instead of easily observed link-time errors.

[Ever_B_]

MKL Lapack SBGVX, a disappointment

I selected SBGVX because it can solve for "selected" eigenvalues/modes (both positive and negative), thinking it would be efficient, but it is extremely slow when compared with Bathe's subspace iteration eigensolver, which unfortunately can solve only for positive eigenvalues. I need only the smallest positive eigenvalue and its mode.

I am trying to solve the problem (K - e * M) u = 0 with M not +definite, and A that is +definite. Since SBGVX needs the second matrix to be +definite, I multiply by -1/e and rearrange it as (M - e' * K) u = 0, then look for the highest e', which is positive, to then compute the lowest positive eigenvalue e=1/e' . It works, but it is very slow. 

Does anyone know how to modify Bathe's subspace iteration eigensolver to compute the highest eigenvalue instead of the smaller?

Or, how can I transform (K - e * M) u = 0 to solve for eigenvalues e^2 instead of e? so I can use Bathe's solver to find the smallest e. Bathe's solver crashes trying to solve for the smallest e when the problem has negative eigenvalues. 

[martenjan]

Maybe I spotted another issue:

07	    ! NJ: number of upper diagonals + 1
08	    integer, parameter :: N=3, NJ=2, ku=NJ-1, kl=0, ldab=ku+1, iwork=7*N
09	    !compute eigenvalues with index il to iu included

iwork should be an array with dimension 5*n, according to the documentation on dsbgvx.

Good luck!

[mecej4]

martenjan wrote:

Maybe I spotted another issue:

08	    integer, parameter :: N=3, NJ=2, ku=NJ-1, kl=0, ldab=ku+1, iwork=7*Niwork should be an array with dimension 5*n, according to the documentation on dsbgvx.

It is always a good idea to check the documentation for the required array size of MKL routine arguments. In this case, however, there need be no concern, because the array argument is an assumed size argument, for which the documentation specifies the minimum size of the array. Passing an array larger than the minimum size is acceptable, if wasteful. In some cases, the same work array may be passed to more than one routine, and in those cases it is natural to declare the work array to be the largest of those specified in the documentation for each routine.