Help with dsyev lapack function in cuda!!

[maurorodas]

Hi, I'm new to this forum, I am a Ph.D. student in Colombia, in my thesis I'm working on an application in quantum chemistry, and I need to pass the LAPACK functions of my application to Cuda, and found with MAGMA, that as I read has routines I need.

One of the routines I use is dsyev LAPACK for eigenvalues ​​and eigenvectors problems, the application is written in Fortran 90, I think the equivalent in MAGMA routines are "magmaf_ssyevd" and "magmaf_dsyevd".

The problem I have is that the input parameters in LAPACK routine are:

jobz, uplo, n, a, lda, w, work, lwork, info.

But on MAGMA calls me two more parameters:

iwork and liwork

I do not know its function, I read the documentation, try to adapt my code but I get no results, someone could help with the interface to the routine work?, I appreciate your help.

The program will run on a NVIDIA Gforce GTX 480.

I apologize for my English, greetings from Colombia

[mgates3]

The LAPACK routine ssyevd should include the same iwork and liwork parameters. See
http://www.netlib.org/lapack/single/ssyevd.f

The MAGMA function is documented in the function itself, e.g., in magma/src/ssyevd.cpp
which is also available in the online documentation at
http://icl.cs.utk.edu/magma/docs/ssyevd_8cpp.html

Let us know if you have further questions.
-mark

[maurorodas]

I've already read, thank you very much, but what I want to do is spend dsyev lapack routine to respective magma routine.

dsyev of lapack routine has no parameters: iwork and liwork, my question is with these parameters, I keep reading anyway, thank you very much

[maurorodas]

Since the two variables thought to use dsyevd, iwork and liwork, but now I get the following error:

forrtl: severe (174): SIGSEGV, segmentation fault occurred

this code works fine with lapack, but I can not work with magma

this is my interface

Code: Select all

module MagmaInterface_

        use Cuda_Manager_, only: sizeof_real
	implicit none
        
        interface MagmaInterface

subroutine magmaf_dsyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
        character          :: jobz
        character          :: uplo
        integer       :: n
        double precision    :: a(*)
        integer       :: lda
        double precision:: w(*)
        double precision    :: work(*)
        integer       :: lwork
        integer       :: iwork(*)
        integer       :: liwork
        integer       :: info
end subroutine magmaf_dsyevd

      end interface MagmaInterface

end module MagmaInterface_

And this is my code:

Code: Select all

subroutine Matrix_eigen( this, eigenValues, eigenVectors, flags, m, dm )
		implicit none
		type(Matrix), intent(in) :: this
		type(Vector), intent(inout) :: eigenValues
		type(Matrix), intent(inout), optional :: eigenVectors
		integer, intent(in), optional :: flags
		integer, intent(in), optional :: dm
		real(8), intent(in), optional :: m(:,:)

		integer :: lengthWorkSpace
		integer :: matrixSize
		integer :: infoProcess
		real(8), allocatable :: workSpace(:)
		type(Matrix) :: eigenVectorsTmp
		integer :: i

  !! Magma variables

                integer, allocatable :: iwork(:)
                integer :: liwork
                
               
		matrixSize = size( this%values, DIM=1 )
                
		if( flags == SYMMETRIC ) then

			lengthWorkSpace=3*matrixSize-1
                        liwork = -1

			allocate( workSpace( lengthWorkSpace ) )
                

			if( present( eigenVectors ) ) then

				if (present ( dm ) ) then
					eigenvectors%values = m
                                    
				else	

				eigenVectors%values=this%values
                             

				end if

				call dsyevd( &
					COMPUTE_EIGENVALUES_AND_EIGENVECTORS, &
					UPPER_TRIANGLE_IS_STORED, &
					matrixSize, &
					eigenVectors%values, &
                                        matrixSize, &
					eigenValues%values, &
                                   	workSpace, &
					lengthWorkSpace, &
                                        iwork, &
                                        liwork, &
					infoProcess )

    

			else
				call Matrix_copyConstructor( eigenVectorsTmp, this )

				call dsyevd( &
					COMPUTE_EIGENVALUES, &
					UPPER_TRIANGLE_IS_STORED, &
					matrixSize, &
					eigenVectorsTmp%values, &
					matrixSize, &
					eigenValues%values, &
                                       	workSpace, &
					lengthWorkSpace, &
                                        iwork, &
                                        liwork, &
					infoProcess )

				call Matrix_destructor( eigenVectorsTmp )

			end if

			
   WRITE(*,*) 'infoProcess = ', infoProcess
			if ( infoProcess /= 0 )  then

				call Matrix_exception(WARNING, "Diagonalization failed", "Class object Matrix in the getEigen() function")

			end if

			do i=1,size(eigenValues%values)
				if( eigenValues%values(i) == Math_NaN ) then

					call Matrix_exception(WARNING, "Diagonalization failed", "Class object Matrix in the getEigen() function")
				end if
			end do


			deallocate(workSpace)

		end if

	end subroutine Matrix_eigen

appreciate all your help

[ffdd]

Could you solve this issue? I have a quite similar problem now and would be very happy to see a working example of dsyevd magma call in fortran.
Thank you in advance

[mgates3]

MAGMA provides a Fortran interface to dsyevd. Here's an example.

Code: Select all

!! dsyevd.f90
!! example compilation and run:
!!
!! gfortran -I $MAGMA_DIR/include -o dsyevd dsyevd.f90 -L $MAGMA_DIR/lib -lmagma -L $OPENBLASDIR/lib -lopenblas
!! setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${MAGMA_DIR}/lib
!! ./dsyevd

program test
    use magma
    implicit none
    
    integer :: i, n, lda, lwork, liwork, info
    
    integer,          allocatable :: iwork(:)
    double precision, allocatable :: w(:)
    double precision, allocatable :: A(:,:), work(:)
    
    integer          :: itmp(1)
    double precision :: rtmp(1)
    double precision :: wtmp(1)
    
    integer :: seed(4)
    
    n   = 500
    lda = n
    
    call magmaf_init()
    
    !! query for workspace sizes
    lwork  = -1
    liwork = -1
    call magmaf_dsyevd( 'v', 'l', n, A, lda, w, wtmp, lwork, itmp, liwork, info )
    lwork  = wtmp(1)
    liwork = itmp(1)
    print '(a,i10,a,i10,a,i10,a,i10)', 'n', n, ', lwork', lwork, ', liwork', liwork
    
    allocate( A( lda, n ))
    allocate( w( n ))
    allocate(  work(  lwork ))
    allocate( iwork( liwork ))
    
    !! initialize A to random values -- replace with your application code
    seed = (/ 0, 1, 2, 3 /)
    call dlarnv( 1, seed, lda*n, A )
    
    !! compute eigenvalues
    call magmaf_dsyevd( 'v', 'l', n, A, lda, w, work, lwork, iwork, liwork, info )
    print '(a,i10)', 'info=', info
    do i=1,n
        print '(a,i4,a,f10.4)', 'w(', i, ')=', w(i)
    end do
    
    deallocate( A )
    deallocate(  work )
    deallocate( iwork )
    
    call magmaf_finalize()
    
end program

[ffdd]

Thank you. But my code is almost exactly the same, if you compare. At least I am sure now, I am on a right way this code.
But still does not find why it is not working for me.

[haidar]

magma_ssyevd is the equivalent of the LAPACK SSYEVD routine which use Divide and Conquer to solver the eigenvalue problem and is faster than SSYEV.
Now you can use the magma example or what Mark mentioned above and try to integrate your stuff.

By the way SSYEV and SSYEVD and all the SSYXXX are for symmetric matrices I saw that you have a flag in your code that switch on SYMETRY flag " if( flags == SYMMETRIC ) then"
Also when you put liwork=-1 the first run of the code is going to give you on output the required size of your workspace which you already did allocate before the first call.
soemthing os wrong here in the order.
either you read the documentation and allocate the required workspace or you put liwork lwork=-1 and you call the routine to return the correct size required you allocate and then you re run to compute.
look at the magma testing directory.

By the way the SSYEVDX use the GPU for the eigensolver and might be faster than SSYEVD
The SSYEVDX_2stage might also be of interest if your matrix is large.
Azzam