zheev

[clement_fevrier]

Hello,

I installed Magma. I would like to use it for zheev. I didn't find any documentation on it, so I guess it's not present in Magma. What alternative do I have?

Maybe zheevd? But I don't know how to "convert" my code to use it. Here the part of my c++ code:

Code: Select all

	unsigned int info = 0;
	const unsigned int lwork = 2 * ncf;
	std::complex< double > * work = new std::complex< double > [lwork];
	double * rwork = new double [3 * ncf - 2];
	double * w  = new double [ncf];
	std::complex< double > * m  = new std::complex< double > [ncf * ncf];
	for(unsigned int i = 0; i < ncf * ncf; ++i)
	{
		m[i] = Hamiltonian[i]; // Where Hamiltonian is a Hermitan Matrix
	}

	zheev_('V', 'U', ncf, m, ncf, w, work, lwork, rwork, info);

Thank you

[mgates3]

Yes, use zheevd. It's faster than zheev. The LAPACK calls are similar:

subroutine zheev (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)

subroutine zheevd (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)

The only changes are to pass the size of rwork as lrwork, and add an integer workspace iwork with size liwork.
Check the MAGMA documentation for the required lwork, lrwork, liwork sizes.
http://icl.cs.utk.edu/projectsfiles/mag ... river.html

For an example, see testing/testing_zheevd.cpp

Note, info should be int (i.e., signed), not unsigned int, as it can return both positive and negative errors. lwork and other variables passed to zheev should probably also be signed, as that is what Fortran uses.

NB: I recommend using a different variable name than "m" for the matrix (e.g., use A or H), as m is usually reserved for the matrix dimension (m x n).

-mark

[clement_fevrier]

Thanks a lot for the answer.

Here the new part of the code

Code: Select all

	magmaDoubleComplex * A  = new magmaDoubleComplex [ncf * ncf];
	for(unsigned int i = 0; i < ncf * ncf; ++i)
	{
		A[i] = make_cuDoubleComplex(Hamiltonian[i].real(), Hamiltonian[i].imag());
	}

	magma_vec_t jobz = MagmaVec;
	magma_uplo_t uplo = MagmaUpper;
	magma_int_t N = (magma_int_t) ncf;
	magma_int_t info = 0;
	magma_int_t lwork = 2 * N;
	magma_int_t lrwork = lwork;
	magma_int_t liwork;
	magma_int_t lda = N;
	magma_int_t * iwork;
	double * rwork = new double [3 * ncf - 2];
	magmaDoubleComplex * work = new magmaDoubleComplex [2 * ncf];
	double * w  = new double [ncf];
	magma_zheevd(jobz, uplo,
		N,
		A,
		lda,
		w,
		work,
		lwork,
		rwork,
		lrwork,
		iwork,
		liwork,
		&info);

I don't know if it's correct.

[clement_fevrier]

I have this error when I launch it

Code: Select all

Error in magma_getdevice_arch: MAGMA not initialized (call magma_init() first) or bad device
Segmentation fault (core dumped)

I guess I miss something.

[mgates3]

As it says, you need to add magma_init(). See the simple example in magma/example/example.c

Always check info after MAGMA and LAPACK calls. It will tell you if arguments are wrong. In this case, you need to adjust your workspace sizes (lwork, lrwork, liwork) to match MAGMA's documentation. Even in LAPACK, zheevd has different workspace sizes than zheev.

I would use lwork, lrwork, liwork as arguments to new/malloc when you allocate arrays, to ensure that they are consistent. In this case, you set lrwork = lwork = 2*N, but allocate rwork as 3*N - 2. Better to do this:

Code: Select all

lwork = max( N + N*NB, 2*N + N*N );
lrwork = 1 + 5*N + 2*N*N;
liwork = 3 + 5*N;
work = new magmaDoubleComplex[ lwork ];
rwork = new double[ lrwork ];
iwork = new magma_int_t[ liwork ];  // was missing

How is Hamiltonian defined? There's probably no need to copy it to a new matrix A. It might just need a cast between std::complex and magmaDoubleComplex.

[clement_fevrier]

Thanks for the help and the advices. It works!

Here the new part of the code

Code: Select all

	magma_init();
	magmaDoubleComplex * A  = new magmaDoubleComplex [ncf*ncf];
	for(unsigned int i = 0; i < ncf*ncf; ++i)
	{
		A[i] = make_cuDoubleComplex(Hamiltonian[i].real(), Hamiltonian[i].imag());
	}
	magma_vec_t jobz = MagmaVec;
	magma_uplo_t uplo = MagmaUpper;
	magma_int_t N = (magma_int_t) ncf;
	magma_int_t info = 0;
	magma_int_t lwork = std::max( N + N*magma_get_zhetrd_nb(N), 2*N + N*N);
	magma_int_t lrwork = 1 + 5*N + 2 * N*N;;
	magma_int_t liwork = 3 + 5*N;
	magma_int_t lda = N;
	magma_int_t * iwork = new magma_int_t [liwork];
	double * rwork = new double [lrwork];
	magmaDoubleComplex * work = new magmaDoubleComplex [lwork];
	double * w  = new double [ncf];
	magma_zheevd(jobz, uplo,
		N,
		A,
		lda,
		w,
		work,
		lwork,
		rwork,
		lrwork,
		iwork,
		liwork,
		&info);

It would be great if I won't have to copy the Hamiltonian. It's defined as
std::complex< double > * Hamiltonian
But I don't know how I can do it.

[mgates3]

See attached file for ways to cast between magmaDoubleComplex, std::complex, and std::vector.
-mark