- #1
- 3,971
- 328
Hello guys, I am in need of a subroutine that diagonalizes a nxn Hermitian matrix (really I'm just looking for the eigenvalues and eigenvectors of course). Looking online I found that LAPACK has a zheev subroutine that presumably does just this. I am a little confused on how to use this subroutine though, I'd appreciate it if you guys could take a look and help me with some of my questions:
So, the above are the comments on what all the inputs and outputs of this subroutine are. There's a few things I'm not sure of.
1) It looks to me like this subroutine will take a nxn matrix A and then replace that matrix with a matrix of its eigenvectors, is that correct?
In this case, I should first create a dummy array for A so that my matrix is not destroyed in the process of calling this subroutine correct?
2) What is this LDA? I can't make heads or tails of it other than it should be "N" since Hermitian matrices are square...I don't understand why someone would write this redundant parameter if that's the case though.
3) I have no idea what's the use of WORK, LWORK, and RWORK. Can anyone make heads or tails of this?
Next, take a look at how this subroutine is defined:
I am completely confused on how the arrays are defined. I am used to definitions of arrays like:
complex,dimension(N,N)::A
How come his subroutine's array definitions don't seem to require a "dimension" declaration? How does the compiler know how much memory to allocate to these arrays?
Not to mention the syntax lacks a "::" now. Zheev should be written in fortran, it's extension is .f. I am not familiar with this syntax at all...
Code:
* Definition:
* ===========
*
* SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
* INFO )
*
* .. Scalar Arguments ..
* CHARACTER JOBZ, UPLO
* INTEGER INFO, LDA, LWORK, N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION RWORK( * ), W( * )
* COMPLEX*16 A( LDA, * ), WORK( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
*> complex Hermitian matrix A.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] JOBZ
*> \verbatim
*> JOBZ is CHARACTER*1
*> = 'N': Compute eigenvalues only;
*> = 'V': Compute eigenvalues and eigenvectors.
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> = 'U': Upper triangle of A is stored;
*> = 'L': Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA, N)
*> On entry, the Hermitian matrix A. If UPLO = 'U', the
*> leading N-by-N upper triangular part of A contains the
*> upper triangular part of the matrix A. If UPLO = 'L',
*> the leading N-by-N lower triangular part of A contains
*> the lower triangular part of the matrix A.
*> On exit, if JOBZ = 'V', then if INFO = 0, 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 destroyed.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] W
*> \verbatim
*> W is DOUBLE PRECISION array, dimension (N)
*> If INFO = 0, the eigenvalues in ascending order.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The length of the array WORK. LWORK >= max(1,2*N-1).
*> For optimal efficiency, LWORK >= (NB+1)*N,
*> where NB is the blocksize for ZHETRD returned by ILAENV.
*>
*> If LWORK = -1, then a workspace query is assumed; the routine
*> only calculates the optimal size of the WORK array, returns
*> this value as the first entry of the WORK array, and no error
*> message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> > 0: if INFO = i, the algorithm failed to converge; i
*> off-diagonal elements of an intermediate tridiagonal
*> form did not converge to zero.
*> \endverbatimSo, the above are the comments on what all the inputs and outputs of this subroutine are. There's a few things I'm not sure of.
1) It looks to me like this subroutine will take a nxn matrix A and then replace that matrix with a matrix of its eigenvectors, is that correct?
In this case, I should first create a dummy array for A so that my matrix is not destroyed in the process of calling this subroutine correct?
2) What is this LDA? I can't make heads or tails of it other than it should be "N" since Hermitian matrices are square...I don't understand why someone would write this redundant parameter if that's the case though.
3) I have no idea what's the use of WORK, LWORK, and RWORK. Can anyone make heads or tails of this?
Next, take a look at how this subroutine is defined:
Code:
SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
$ INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
*
* .. Scalar Arguments ..
CHARACTER JOBZ, UPLO
INTEGER INFO, LDA, LWORK, N
* ..
* .. Array Arguments ..
DOUBLE PRECISION RWORK( * ), W( * )
COMPLEX*16 A( LDA, * ), WORK( * )I am completely confused on how the arrays are defined. I am used to definitions of arrays like:
complex,dimension(N,N)::A
How come his subroutine's array definitions don't seem to require a "dimension" declaration? How does the compiler know how much memory to allocate to these arrays?
Not to mention the syntax lacks a "::" now. Zheev should be written in fortran, it's extension is .f. I am not familiar with this syntax at all...
Last edited: