LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
sgeev.f
Go to the documentation of this file.
1 *> \brief <b> SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeev.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeev.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeev.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
22 * LDVR, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBVL, JOBVR
26 * INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
30 * \$ WI( * ), WORK( * ), WR( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SGEEV computes for an N-by-N real nonsymmetric matrix A, the
40 *> eigenvalues and, optionally, the left and/or right eigenvectors.
41 *>
42 *> The right eigenvector v(j) of A satisfies
43 *> A * v(j) = lambda(j) * v(j)
44 *> where lambda(j) is its eigenvalue.
45 *> The left eigenvector u(j) of A satisfies
46 *> u(j)**H * A = lambda(j) * u(j)**H
47 *> where u(j)**H denotes the conjugate-transpose of u(j).
48 *>
49 *> The computed eigenvectors are normalized to have Euclidean norm
50 *> equal to 1 and largest component real.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] JOBVL
57 *> \verbatim
58 *> JOBVL is CHARACTER*1
59 *> = 'N': left eigenvectors of A are not computed;
60 *> = 'V': left eigenvectors of A are computed.
61 *> \endverbatim
62 *>
63 *> \param[in] JOBVR
64 *> \verbatim
65 *> JOBVR is CHARACTER*1
66 *> = 'N': right eigenvectors of A are not computed;
67 *> = 'V': right eigenvectors of A are computed.
68 *> \endverbatim
69 *>
70 *> \param[in] N
71 *> \verbatim
72 *> N is INTEGER
73 *> The order of the matrix A. N >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in,out] A
77 *> \verbatim
78 *> A is REAL array, dimension (LDA,N)
79 *> On entry, the N-by-N matrix A.
80 *> On exit, A has been overwritten.
81 *> \endverbatim
82 *>
83 *> \param[in] LDA
84 *> \verbatim
85 *> LDA is INTEGER
86 *> The leading dimension of the array A. LDA >= max(1,N).
87 *> \endverbatim
88 *>
89 *> \param[out] WR
90 *> \verbatim
91 *> WR is REAL array, dimension (N)
92 *> \endverbatim
93 *>
94 *> \param[out] WI
95 *> \verbatim
96 *> WI is REAL array, dimension (N)
97 *> WR and WI contain the real and imaginary parts,
98 *> respectively, of the computed eigenvalues. Complex
99 *> conjugate pairs of eigenvalues appear consecutively
100 *> with the eigenvalue having the positive imaginary part
101 *> first.
102 *> \endverbatim
103 *>
104 *> \param[out] VL
105 *> \verbatim
106 *> VL is REAL array, dimension (LDVL,N)
107 *> If JOBVL = 'V', the left eigenvectors u(j) are stored one
108 *> after another in the columns of VL, in the same order
109 *> as their eigenvalues.
110 *> If JOBVL = 'N', VL is not referenced.
111 *> If the j-th eigenvalue is real, then u(j) = VL(:,j),
112 *> the j-th column of VL.
113 *> If the j-th and (j+1)-st eigenvalues form a complex
114 *> conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
115 *> u(j+1) = VL(:,j) - i*VL(:,j+1).
116 *> \endverbatim
117 *>
118 *> \param[in] LDVL
119 *> \verbatim
120 *> LDVL is INTEGER
121 *> The leading dimension of the array VL. LDVL >= 1; if
122 *> JOBVL = 'V', LDVL >= N.
123 *> \endverbatim
124 *>
125 *> \param[out] VR
126 *> \verbatim
127 *> VR is REAL array, dimension (LDVR,N)
128 *> If JOBVR = 'V', the right eigenvectors v(j) are stored one
129 *> after another in the columns of VR, in the same order
130 *> as their eigenvalues.
131 *> If JOBVR = 'N', VR is not referenced.
132 *> If the j-th eigenvalue is real, then v(j) = VR(:,j),
133 *> the j-th column of VR.
134 *> If the j-th and (j+1)-st eigenvalues form a complex
135 *> conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
136 *> v(j+1) = VR(:,j) - i*VR(:,j+1).
137 *> \endverbatim
138 *>
139 *> \param[in] LDVR
140 *> \verbatim
141 *> LDVR is INTEGER
142 *> The leading dimension of the array VR. LDVR >= 1; if
143 *> JOBVR = 'V', LDVR >= N.
144 *> \endverbatim
145 *>
146 *> \param[out] WORK
147 *> \verbatim
148 *> WORK is REAL array, dimension (MAX(1,LWORK))
149 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
150 *> \endverbatim
151 *>
152 *> \param[in] LWORK
153 *> \verbatim
154 *> LWORK is INTEGER
155 *> The dimension of the array WORK. LWORK >= max(1,3*N), and
156 *> if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
157 *> performance, LWORK must generally be larger.
158 *>
159 *> If LWORK = -1, then a workspace query is assumed; the routine
160 *> only calculates the optimal size of the WORK array, returns
161 *> this value as the first entry of the WORK array, and no error
162 *> message related to LWORK is issued by XERBLA.
163 *> \endverbatim
164 *>
165 *> \param[out] INFO
166 *> \verbatim
167 *> INFO is INTEGER
168 *> = 0: successful exit
169 *> < 0: if INFO = -i, the i-th argument had an illegal value.
170 *> > 0: if INFO = i, the QR algorithm failed to compute all the
171 *> eigenvalues, and no eigenvectors have been computed;
172 *> elements i+1:N of WR and WI contain eigenvalues which
173 *> have converged.
174 *> \endverbatim
175 *
176 * Authors:
177 * ========
178 *
179 *> \author Univ. of Tennessee
180 *> \author Univ. of California Berkeley
181 *> \author Univ. of Colorado Denver
182 *> \author NAG Ltd.
183 *
184 *> \date September 2012
185 *
186 *> \ingroup realGEeigen
187 *
188 * =====================================================================
189  SUBROUTINE sgeev( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
190  \$ ldvr, work, lwork, info )
191 *
192 * -- LAPACK driver routine (version 3.4.2) --
193 * -- LAPACK is a software package provided by Univ. of Tennessee, --
194 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195 * September 2012
196 *
197 * .. Scalar Arguments ..
198  CHARACTER jobvl, jobvr
199  INTEGER info, lda, ldvl, ldvr, lwork, n
200 * ..
201 * .. Array Arguments ..
202  REAL a( lda, * ), vl( ldvl, * ), vr( ldvr, * ),
203  \$ wi( * ), work( * ), wr( * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  REAL zero, one
210  parameter( zero = 0.0e0, one = 1.0e0 )
211 * ..
212 * .. Local Scalars ..
213  LOGICAL lquery, scalea, wantvl, wantvr
214  CHARACTER side
215  INTEGER hswork, i, ibal, ierr, ihi, ilo, itau, iwrk, k,
216  \$ maxwrk, minwrk, nout
217  REAL anrm, bignum, cs, cscale, eps, r, scl, smlnum,
218  \$ sn
219 * ..
220 * .. Local Arrays ..
221  LOGICAL select( 1 )
222  REAL dum( 1 )
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL sgebak, sgebal, sgehrd, shseqr, slabad, slacpy,
227  \$ xerbla
228 * ..
229 * .. External Functions ..
230  LOGICAL lsame
231  INTEGER ilaenv, isamax
232  REAL slamch, slange, slapy2, snrm2
233  EXTERNAL lsame, ilaenv, isamax, slamch, slange, slapy2,
234  \$ snrm2
235 * ..
236 * .. Intrinsic Functions ..
237  INTRINSIC max, sqrt
238 * ..
239 * .. Executable Statements ..
240 *
241 * Test the input arguments
242 *
243  info = 0
244  lquery = ( lwork.EQ.-1 )
245  wantvl = lsame( jobvl, 'V' )
246  wantvr = lsame( jobvr, 'V' )
247  IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
248  info = -1
249  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
250  info = -2
251  ELSE IF( n.LT.0 ) THEN
252  info = -3
253  ELSE IF( lda.LT.max( 1, n ) ) THEN
254  info = -5
255  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
256  info = -9
257  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
258  info = -11
259  END IF
260 *
261 * Compute workspace
262 * (Note: Comments in the code beginning "Workspace:" describe the
263 * minimal amount of workspace needed at that point in the code,
264 * as well as the preferred amount for good performance.
265 * NB refers to the optimal block size for the immediately
266 * following subroutine, as returned by ILAENV.
267 * HSWORK refers to the workspace preferred by SHSEQR, as
268 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
269 * the worst case.)
270 *
271  IF( info.EQ.0 ) THEN
272  IF( n.EQ.0 ) THEN
273  minwrk = 1
274  maxwrk = 1
275  ELSE
276  maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
277  IF( wantvl ) THEN
278  minwrk = 4*n
279  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
280  \$ 'SORGHR', ' ', n, 1, n, -1 ) )
281  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
282  \$ work, -1, info )
283  hswork = work( 1 )
284  maxwrk = max( maxwrk, n + 1, n + hswork )
285  maxwrk = max( maxwrk, 4*n )
286  ELSE IF( wantvr ) THEN
287  minwrk = 4*n
288  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
289  \$ 'SORGHR', ' ', n, 1, n, -1 ) )
290  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
291  \$ work, -1, info )
292  hswork = work( 1 )
293  maxwrk = max( maxwrk, n + 1, n + hswork )
294  maxwrk = max( maxwrk, 4*n )
295  ELSE
296  minwrk = 3*n
297  CALL shseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
298  \$ work, -1, info )
299  hswork = work( 1 )
300  maxwrk = max( maxwrk, n + 1, n + hswork )
301  END IF
302  maxwrk = max( maxwrk, minwrk )
303  END IF
304  work( 1 ) = maxwrk
305 *
306  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
307  info = -13
308  END IF
309  END IF
310 *
311  IF( info.NE.0 ) THEN
312  CALL xerbla( 'SGEEV ', -info )
313  return
314  ELSE IF( lquery ) THEN
315  return
316  END IF
317 *
318 * Quick return if possible
319 *
320  IF( n.EQ.0 )
321  \$ return
322 *
323 * Get machine constants
324 *
325  eps = slamch( 'P' )
326  smlnum = slamch( 'S' )
327  bignum = one / smlnum
328  CALL slabad( smlnum, bignum )
329  smlnum = sqrt( smlnum ) / eps
330  bignum = one / smlnum
331 *
332 * Scale A if max element outside range [SMLNUM,BIGNUM]
333 *
334  anrm = slange( 'M', n, n, a, lda, dum )
335  scalea = .false.
336  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
337  scalea = .true.
338  cscale = smlnum
339  ELSE IF( anrm.GT.bignum ) THEN
340  scalea = .true.
341  cscale = bignum
342  END IF
343  IF( scalea )
344  \$ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
345 *
346 * Balance the matrix
347 * (Workspace: need N)
348 *
349  ibal = 1
350  CALL sgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
351 *
352 * Reduce to upper Hessenberg form
353 * (Workspace: need 3*N, prefer 2*N+N*NB)
354 *
355  itau = ibal + n
356  iwrk = itau + n
357  CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
358  \$ lwork-iwrk+1, ierr )
359 *
360  IF( wantvl ) THEN
361 *
362 * Want left eigenvectors
363 * Copy Householder vectors to VL
364 *
365  side = 'L'
366  CALL slacpy( 'L', n, n, a, lda, vl, ldvl )
367 *
368 * Generate orthogonal matrix in VL
369 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
370 *
371  CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
372  \$ lwork-iwrk+1, ierr )
373 *
374 * Perform QR iteration, accumulating Schur vectors in VL
375 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
376 *
377  iwrk = itau
378  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
379  \$ work( iwrk ), lwork-iwrk+1, info )
380 *
381  IF( wantvr ) THEN
382 *
383 * Want left and right eigenvectors
384 * Copy Schur vectors to VR
385 *
386  side = 'B'
387  CALL slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
388  END IF
389 *
390  ELSE IF( wantvr ) THEN
391 *
392 * Want right eigenvectors
393 * Copy Householder vectors to VR
394 *
395  side = 'R'
396  CALL slacpy( 'L', n, n, a, lda, vr, ldvr )
397 *
398 * Generate orthogonal matrix in VR
399 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
400 *
401  CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
402  \$ lwork-iwrk+1, ierr )
403 *
404 * Perform QR iteration, accumulating Schur vectors in VR
405 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
406 *
407  iwrk = itau
408  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
409  \$ work( iwrk ), lwork-iwrk+1, info )
410 *
411  ELSE
412 *
413 * Compute eigenvalues only
414 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
415 *
416  iwrk = itau
417  CALL shseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
418  \$ work( iwrk ), lwork-iwrk+1, info )
419  END IF
420 *
421 * If INFO > 0 from SHSEQR, then quit
422 *
423  IF( info.GT.0 )
424  \$ go to 50
425 *
426  IF( wantvl .OR. wantvr ) THEN
427 *
428 * Compute left and/or right eigenvectors
429 * (Workspace: need 4*N)
430 *
431  CALL strevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
432  \$ n, nout, work( iwrk ), ierr )
433  END IF
434 *
435  IF( wantvl ) THEN
436 *
437 * Undo balancing of left eigenvectors
438 * (Workspace: need N)
439 *
440  CALL sgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
441  \$ ierr )
442 *
443 * Normalize left eigenvectors and make largest component real
444 *
445  DO 20 i = 1, n
446  IF( wi( i ).EQ.zero ) THEN
447  scl = one / snrm2( n, vl( 1, i ), 1 )
448  CALL sscal( n, scl, vl( 1, i ), 1 )
449  ELSE IF( wi( i ).GT.zero ) THEN
450  scl = one / slapy2( snrm2( n, vl( 1, i ), 1 ),
451  \$ snrm2( n, vl( 1, i+1 ), 1 ) )
452  CALL sscal( n, scl, vl( 1, i ), 1 )
453  CALL sscal( n, scl, vl( 1, i+1 ), 1 )
454  DO 10 k = 1, n
455  work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
456  10 continue
457  k = isamax( n, work( iwrk ), 1 )
458  CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
459  CALL srot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
460  vl( k, i+1 ) = zero
461  END IF
462  20 continue
463  END IF
464 *
465  IF( wantvr ) THEN
466 *
467 * Undo balancing of right eigenvectors
468 * (Workspace: need N)
469 *
470  CALL sgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
471  \$ ierr )
472 *
473 * Normalize right eigenvectors and make largest component real
474 *
475  DO 40 i = 1, n
476  IF( wi( i ).EQ.zero ) THEN
477  scl = one / snrm2( n, vr( 1, i ), 1 )
478  CALL sscal( n, scl, vr( 1, i ), 1 )
479  ELSE IF( wi( i ).GT.zero ) THEN
480  scl = one / slapy2( snrm2( n, vr( 1, i ), 1 ),
481  \$ snrm2( n, vr( 1, i+1 ), 1 ) )
482  CALL sscal( n, scl, vr( 1, i ), 1 )
483  CALL sscal( n, scl, vr( 1, i+1 ), 1 )
484  DO 30 k = 1, n
485  work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
486  30 continue
487  k = isamax( n, work( iwrk ), 1 )
488  CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
489  CALL srot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
490  vr( k, i+1 ) = zero
491  END IF
492  40 continue
493  END IF
494 *
495 * Undo scaling if necessary
496 *
497  50 continue
498  IF( scalea ) THEN
499  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
500  \$ max( n-info, 1 ), ierr )
501  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
502  \$ max( n-info, 1 ), ierr )
503  IF( info.GT.0 ) THEN
504  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
505  \$ ierr )
506  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
507  \$ ierr )
508  END IF
509  END IF
510 *
511  work( 1 ) = maxwrk
512  return
513 *
514 * End of SGEEV
515 *
516  END