LAPACK  3.6.1 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 June 2016
185 *
186 * @generated from dgeev.f, fortran d -> s, Tue Apr 19 01:47:44 2016
187 *
188 *> \ingroup realGEeigen
189 *
190 * =====================================================================
191  SUBROUTINE sgeev( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
192  \$ ldvr, work, lwork, info )
193  implicit none
194 *
195 * -- LAPACK driver routine (version 3.6.1) --
196 * -- LAPACK is a software package provided by Univ. of Tennessee, --
197 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198 * June 2016
199 *
200 * .. Scalar Arguments ..
201  CHARACTER JOBVL, JOBVR
202  INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
203 * ..
204 * .. Array Arguments ..
205  REAL A( lda, * ), VL( ldvl, * ), VR( ldvr, * ),
206  \$ wi( * ), work( * ), wr( * )
207 * ..
208 *
209 * =====================================================================
210 *
211 * .. Parameters ..
212  REAL ZERO, ONE
213  parameter ( zero = 0.0e0, one = 1.0e0 )
214 * ..
215 * .. Local Scalars ..
216  LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
217  CHARACTER SIDE
218  INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
219  \$ lwork_trevc, maxwrk, minwrk, nout
220  REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
221  \$ sn
222 * ..
223 * .. Local Arrays ..
224  LOGICAL SELECT( 1 )
225  REAL DUM( 1 )
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL sgebak, sgebal, sgehrd, shseqr, slabad, slacpy,
230  \$ xerbla
231 * ..
232 * .. External Functions ..
233  LOGICAL LSAME
234  INTEGER ISAMAX, ILAENV
235  REAL SLAMCH, SLANGE, SLAPY2, SNRM2
236  EXTERNAL lsame, isamax, ilaenv, slamch, slange, slapy2,
237  \$ snrm2
238 * ..
239 * .. Intrinsic Functions ..
240  INTRINSIC max, sqrt
241 * ..
242 * .. Executable Statements ..
243 *
244 * Test the input arguments
245 *
246  info = 0
247  lquery = ( lwork.EQ.-1 )
248  wantvl = lsame( jobvl, 'V' )
249  wantvr = lsame( jobvr, 'V' )
250  IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
251  info = -1
252  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
253  info = -2
254  ELSE IF( n.LT.0 ) THEN
255  info = -3
256  ELSE IF( lda.LT.max( 1, n ) ) THEN
257  info = -5
258  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
259  info = -9
260  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
261  info = -11
262  END IF
263 *
264 * Compute workspace
265 * (Note: Comments in the code beginning "Workspace:" describe the
266 * minimal amount of workspace needed at that point in the code,
267 * as well as the preferred amount for good performance.
268 * NB refers to the optimal block size for the immediately
269 * following subroutine, as returned by ILAENV.
270 * HSWORK refers to the workspace preferred by SHSEQR, as
271 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
272 * the worst case.)
273 *
274  IF( info.EQ.0 ) THEN
275  IF( n.EQ.0 ) THEN
276  minwrk = 1
277  maxwrk = 1
278  ELSE
279  maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
280  IF( wantvl ) THEN
281  minwrk = 4*n
282  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
283  \$ 'SORGHR', ' ', n, 1, n, -1 ) )
284  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
285  \$ work, -1, info )
286  hswork = int( work(1) )
287  maxwrk = max( maxwrk, n + 1, n + hswork )
288  CALL strevc3( 'L', 'B', SELECT, n, a, lda,
289  \$ vl, ldvl, vr, ldvr, n, nout,
290  \$ work, -1, ierr )
291  lwork_trevc = int( work(1) )
292  maxwrk = max( maxwrk, n + lwork_trevc )
293  maxwrk = max( maxwrk, 4*n )
294  ELSE IF( wantvr ) THEN
295  minwrk = 4*n
296  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
297  \$ 'SORGHR', ' ', n, 1, n, -1 ) )
298  CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
299  \$ work, -1, info )
300  hswork = int( work(1) )
301  maxwrk = max( maxwrk, n + 1, n + hswork )
302  CALL strevc3( 'R', 'B', SELECT, n, a, lda,
303  \$ vl, ldvl, vr, ldvr, n, nout,
304  \$ work, -1, ierr )
305  lwork_trevc = int( work(1) )
306  maxwrk = max( maxwrk, n + lwork_trevc )
307  maxwrk = max( maxwrk, 4*n )
308  ELSE
309  minwrk = 3*n
310  CALL shseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
311  \$ work, -1, info )
312  hswork = int( work(1) )
313  maxwrk = max( maxwrk, n + 1, n + hswork )
314  END IF
315  maxwrk = max( maxwrk, minwrk )
316  END IF
317  work( 1 ) = maxwrk
318 *
319  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
320  info = -13
321  END IF
322  END IF
323 *
324  IF( info.NE.0 ) THEN
325  CALL xerbla( 'SGEEV ', -info )
326  RETURN
327  ELSE IF( lquery ) THEN
328  RETURN
329  END IF
330 *
331 * Quick return if possible
332 *
333  IF( n.EQ.0 )
334  \$ RETURN
335 *
336 * Get machine constants
337 *
338  eps = slamch( 'P' )
339  smlnum = slamch( 'S' )
340  bignum = one / smlnum
341  CALL slabad( smlnum, bignum )
342  smlnum = sqrt( smlnum ) / eps
343  bignum = one / smlnum
344 *
345 * Scale A if max element outside range [SMLNUM,BIGNUM]
346 *
347  anrm = slange( 'M', n, n, a, lda, dum )
348  scalea = .false.
349  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
350  scalea = .true.
351  cscale = smlnum
352  ELSE IF( anrm.GT.bignum ) THEN
353  scalea = .true.
354  cscale = bignum
355  END IF
356  IF( scalea )
357  \$ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
358 *
359 * Balance the matrix
360 * (Workspace: need N)
361 *
362  ibal = 1
363  CALL sgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
364 *
365 * Reduce to upper Hessenberg form
366 * (Workspace: need 3*N, prefer 2*N+N*NB)
367 *
368  itau = ibal + n
369  iwrk = itau + n
370  CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
371  \$ lwork-iwrk+1, ierr )
372 *
373  IF( wantvl ) THEN
374 *
375 * Want left eigenvectors
376 * Copy Householder vectors to VL
377 *
378  side = 'L'
379  CALL slacpy( 'L', n, n, a, lda, vl, ldvl )
380 *
381 * Generate orthogonal matrix in VL
382 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
383 *
384  CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
385  \$ lwork-iwrk+1, ierr )
386 *
387 * Perform QR iteration, accumulating Schur vectors in VL
388 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
389 *
390  iwrk = itau
391  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
392  \$ work( iwrk ), lwork-iwrk+1, info )
393 *
394  IF( wantvr ) THEN
395 *
396 * Want left and right eigenvectors
397 * Copy Schur vectors to VR
398 *
399  side = 'B'
400  CALL slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
401  END IF
402 *
403  ELSE IF( wantvr ) THEN
404 *
405 * Want right eigenvectors
406 * Copy Householder vectors to VR
407 *
408  side = 'R'
409  CALL slacpy( 'L', n, n, a, lda, vr, ldvr )
410 *
411 * Generate orthogonal matrix in VR
412 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
413 *
414  CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
415  \$ lwork-iwrk+1, ierr )
416 *
417 * Perform QR iteration, accumulating Schur vectors in VR
418 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
419 *
420  iwrk = itau
421  CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
422  \$ work( iwrk ), lwork-iwrk+1, info )
423 *
424  ELSE
425 *
426 * Compute eigenvalues only
427 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
428 *
429  iwrk = itau
430  CALL shseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
431  \$ work( iwrk ), lwork-iwrk+1, info )
432  END IF
433 *
434 * If INFO .NE. 0 from SHSEQR, then quit
435 *
436  IF( info.NE.0 )
437  \$ GO TO 50
438 *
439  IF( wantvl .OR. wantvr ) THEN
440 *
441 * Compute left and/or right eigenvectors
442 * (Workspace: need 4*N, prefer N + N + 2*N*NB)
443 *
444  CALL strevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
445  \$ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
446  END IF
447 *
448  IF( wantvl ) THEN
449 *
450 * Undo balancing of left eigenvectors
451 * (Workspace: need N)
452 *
453  CALL sgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
454  \$ ierr )
455 *
456 * Normalize left eigenvectors and make largest component real
457 *
458  DO 20 i = 1, n
459  IF( wi( i ).EQ.zero ) THEN
460  scl = one / snrm2( n, vl( 1, i ), 1 )
461  CALL sscal( n, scl, vl( 1, i ), 1 )
462  ELSE IF( wi( i ).GT.zero ) THEN
463  scl = one / slapy2( snrm2( n, vl( 1, i ), 1 ),
464  \$ snrm2( n, vl( 1, i+1 ), 1 ) )
465  CALL sscal( n, scl, vl( 1, i ), 1 )
466  CALL sscal( n, scl, vl( 1, i+1 ), 1 )
467  DO 10 k = 1, n
468  work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
469  10 CONTINUE
470  k = isamax( n, work( iwrk ), 1 )
471  CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
472  CALL srot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
473  vl( k, i+1 ) = zero
474  END IF
475  20 CONTINUE
476  END IF
477 *
478  IF( wantvr ) THEN
479 *
480 * Undo balancing of right eigenvectors
481 * (Workspace: need N)
482 *
483  CALL sgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
484  \$ ierr )
485 *
486 * Normalize right eigenvectors and make largest component real
487 *
488  DO 40 i = 1, n
489  IF( wi( i ).EQ.zero ) THEN
490  scl = one / snrm2( n, vr( 1, i ), 1 )
491  CALL sscal( n, scl, vr( 1, i ), 1 )
492  ELSE IF( wi( i ).GT.zero ) THEN
493  scl = one / slapy2( snrm2( n, vr( 1, i ), 1 ),
494  \$ snrm2( n, vr( 1, i+1 ), 1 ) )
495  CALL sscal( n, scl, vr( 1, i ), 1 )
496  CALL sscal( n, scl, vr( 1, i+1 ), 1 )
497  DO 30 k = 1, n
498  work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
499  30 CONTINUE
500  k = isamax( n, work( iwrk ), 1 )
501  CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
502  CALL srot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
503  vr( k, i+1 ) = zero
504  END IF
505  40 CONTINUE
506  END IF
507 *
508 * Undo scaling if necessary
509 *
510  50 CONTINUE
511  IF( scalea ) THEN
512  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
513  \$ max( n-info, 1 ), ierr )
514  CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
515  \$ max( n-info, 1 ), ierr )
516  IF( info.GT.0 ) THEN
517  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
518  \$ ierr )
519  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
520  \$ ierr )
521  END IF
522  END IF
523 *
524  work( 1 ) = maxwrk
525  RETURN
526 *
527 * End of SGEEV
528 *
529  END
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine strevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
STREVC3
Definition: strevc3.f:241
subroutine sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:162
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine sgeev(JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: sgeev.f:193