KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsygv.f
Go to the documentation of this file.
1  SUBROUTINE dsygv( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
2  $ LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * June 30, 1999
8 *
9 * .. Scalar Arguments ..
10  CHARACTER JOBZ, UPLO
11  INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSYGV computes all the eigenvalues, and optionally, the eigenvectors
21 * of a real generalized symmetric-definite eigenproblem, of the form
22 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
23 * Here A and B are assumed to be symmetric and B is also
24 * positive definite.
25 *
26 * Arguments
27 * =========
28 *
29 * ITYPE (input) INTEGER
30 * Specifies the problem type to be solved:
31 * = 1: A*x = (lambda)*B*x
32 * = 2: A*B*x = (lambda)*x
33 * = 3: B*A*x = (lambda)*x
34 *
35 * JOBZ (input) CHARACTER*1
36 * = 'N': Compute eigenvalues only;
37 * = 'V': Compute eigenvalues and eigenvectors.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangles of A and B are stored;
41 * = 'L': Lower triangles of A and B are stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrices A and B. N >= 0.
45 *
46 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
47 * On entry, the symmetric matrix A. If UPLO = 'U', the
48 * leading N-by-N upper triangular part of A contains the
49 * upper triangular part of the matrix A. If UPLO = 'L',
50 * the leading N-by-N lower triangular part of A contains
51 * the lower triangular part of the matrix A.
52 *
53 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
54 * matrix Z of eigenvectors. The eigenvectors are normalized
55 * as follows:
56 * if ITYPE = 1 or 2, Z**T*B*Z = I;
57 * if ITYPE = 3, Z**T*inv(B)*Z = I.
58 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
59 * or the lower triangle (if UPLO='L') of A, including the
60 * diagonal, is destroyed.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of the array A. LDA >= max(1,N).
64 *
65 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
66 * On entry, the symmetric positive definite matrix B.
67 * If UPLO = 'U', the leading N-by-N upper triangular part of B
68 * contains the upper triangular part of the matrix B.
69 * If UPLO = 'L', the leading N-by-N lower triangular part of B
70 * contains the lower triangular part of the matrix B.
71 *
72 * On exit, if INFO <= N, the part of B containing the matrix is
73 * overwritten by the triangular factor U or L from the Cholesky
74 * factorization B = U**T*U or B = L*L**T.
75 *
76 * LDB (input) INTEGER
77 * The leading dimension of the array B. LDB >= max(1,N).
78 *
79 * W (output) DOUBLE PRECISION array, dimension (N)
80 * If INFO = 0, the eigenvalues in ascending order.
81 *
82 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
83 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84 *
85 * LWORK (input) INTEGER
86 * The length of the array WORK. LWORK >= max(1,3*N-1).
87 * For optimal efficiency, LWORK >= (NB+2)*N,
88 * where NB is the blocksize for DSYTRD returned by ILAENV.
89 *
90 * If LWORK = -1, then a workspace query is assumed; the routine
91 * only calculates the optimal size of the WORK array, returns
92 * this value as the first entry of the WORK array, and no error
93 * message related to LWORK is issued by XERBLA.
94 *
95 * INFO (output) INTEGER
96 * = 0: successful exit
97 * < 0: if INFO = -i, the i-th argument had an illegal value
98 * > 0: DPOTRF or DSYEV returned an error code:
99 * <= N: if INFO = i, DSYEV failed to converge;
100 * i off-diagonal elements of an intermediate
101 * tridiagonal form did not converge to zero;
102 * > N: if INFO = N + i, for 1 <= i <= N, then the leading
103 * minor of order i of B is not positive definite.
104 * The factorization of B could not be completed and
105 * no eigenvalues or eigenvectors were computed.
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110  DOUBLE PRECISION ONE
111  parameter( one = 1.0d+0 )
112 * ..
113 * .. Local Scalars ..
114  LOGICAL LQUERY, UPPER, WANTZ
115  CHARACTER TRANS
116  INTEGER LWKOPT, NB, NEIG
117 * ..
118 * .. External Functions ..
119  LOGICAL LSAME
120  INTEGER ILAENV
121  EXTERNAL lsame, ilaenv
122 * ..
123 * .. External Subroutines ..
124  EXTERNAL dpotrf, dsyev, dsygst, dtrmm, dtrsm, xerbla
125 * ..
126 * .. Intrinsic Functions ..
127  INTRINSIC max
128 * ..
129 * .. Executable Statements ..
130 *
131 * Test the input parameters.
132 *
133  wantz = lsame( jobz, 'V' )
134  upper = lsame( uplo, 'U' )
135  lquery = ( lwork.EQ.-1 )
136 *
137  info = 0
138  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
139  info = -1
140  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
141  info = -2
142  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
143  info = -3
144  ELSE IF( n.LT.0 ) THEN
145  info = -4
146  ELSE IF( lda.LT.max( 1, n ) ) THEN
147  info = -6
148  ELSE IF( ldb.LT.max( 1, n ) ) THEN
149  info = -8
150  ELSE IF( lwork.LT.max( 1, 3*n-1 ) .AND. .NOT.lquery ) THEN
151  info = -11
152  END IF
153 *
154  IF( info.EQ.0 ) THEN
155  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
156  lwkopt = ( nb+2 )*n
157  work( 1 ) = lwkopt
158  END IF
159 *
160  IF( info.NE.0 ) THEN
161  CALL xerbla( 'DSYGV ', -info )
162  RETURN
163  ELSE IF( lquery ) THEN
164  RETURN
165  END IF
166 *
167 * Quick return if possible
168 *
169  IF( n.EQ.0 )
170  $ RETURN
171 *
172 * Form a Cholesky factorization of B.
173 *
174  CALL dpotrf( uplo, n, b, ldb, info )
175  IF( info.NE.0 ) THEN
176  info = n + info
177  RETURN
178  END IF
179 *
180 * Transform problem to standard eigenvalue problem and solve.
181 *
182  CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
183  CALL dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
184 *
185  IF( wantz ) THEN
186 *
187 * Backtransform eigenvectors to the original problem.
188 *
189  neig = n
190  IF( info.GT.0 )
191  $ neig = info - 1
192  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
193 *
194 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
195 * backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
196 *
197  IF( upper ) THEN
198  trans = 'N'
199  ELSE
200  trans = 'T'
201  END IF
202 *
203  CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
204  $ b, ldb, a, lda )
205 *
206  ELSE IF( itype.EQ.3 ) THEN
207 *
208 * For B*A*x=(lambda)*x;
209 * backtransform eigenvectors: x = L*y or U'*y
210 *
211  IF( upper ) THEN
212  trans = 'T'
213  ELSE
214  trans = 'N'
215  END IF
216 *
217  CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
218  $ b, ldb, a, lda )
219  END IF
220  END IF
221 *
222  work( 1 ) = lwkopt
223  RETURN
224 *
225 * End of DSYGV
226 *
227  END
subroutine dpotrf(UPLO, N, A, LDA, INFO)
Definition: dpotrf.f:2
subroutine dsyev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)
Definition: dsyev.f:2
subroutine dsygst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
Definition: dsygst.f:2
subroutine dsygv(ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO)
Definition: dsygv.f:3
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrmm.f:3
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrsm.f:3
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2