LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
sgeev.f File Reference

Go to the source code of this file.

Functions/Subroutines

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 More...
 

Function/Subroutine Documentation

subroutine sgeev ( character  JOBVL,
character  JOBVR,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  WR,
real, dimension( * )  WI,
real, dimension( ldvl, * )  VL,
integer  LDVL,
real, dimension( ldvr, * )  VR,
integer  LDVR,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Download SGEEV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SGEEV computes for an N-by-N real nonsymmetric matrix A, the
 eigenvalues and, optionally, the left and/or right eigenvectors.

 The right eigenvector v(j) of A satisfies
                  A * v(j) = lambda(j) * v(j)
 where lambda(j) is its eigenvalue.
 The left eigenvector u(j) of A satisfies
               u(j)**H * A = lambda(j) * u(j)**H
 where u(j)**H denotes the conjugate-transpose of u(j).

 The computed eigenvectors are normalized to have Euclidean norm
 equal to 1 and largest component real.
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N': left eigenvectors of A are not computed;
          = 'V': left eigenvectors of A are computed.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N': right eigenvectors of A are not computed;
          = 'V': right eigenvectors of A are computed.
[in]N
          N is INTEGER
          The order of the matrix A. N >= 0.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, the N-by-N matrix A.
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]WR
          WR is REAL array, dimension (N)
[out]WI
          WI is REAL array, dimension (N)
          WR and WI contain the real and imaginary parts,
          respectively, of the computed eigenvalues.  Complex
          conjugate pairs of eigenvalues appear consecutively
          with the eigenvalue having the positive imaginary part
          first.
[out]VL
          VL is REAL array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored one
          after another in the columns of VL, in the same order
          as their eigenvalues.
          If JOBVL = 'N', VL is not referenced.
          If the j-th eigenvalue is real, then u(j) = VL(:,j),
          the j-th column of VL.
          If the j-th and (j+1)-st eigenvalues form a complex
          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
          u(j+1) = VL(:,j) - i*VL(:,j+1).
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.  LDVL >= 1; if
          JOBVL = 'V', LDVL >= N.
[out]VR
          VR is REAL array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors v(j) are stored one
          after another in the columns of VR, in the same order
          as their eigenvalues.
          If JOBVR = 'N', VR is not referenced.
          If the j-th eigenvalue is real, then v(j) = VR(:,j),
          the j-th column of VR.
          If the j-th and (j+1)-st eigenvalues form a complex
          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
          v(j+1) = VR(:,j) - i*VR(:,j+1).
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1; if
          JOBVR = 'V', LDVR >= N.
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,3*N), and
          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
          performance, LWORK must generally be larger.

          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.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = i, the QR algorithm failed to compute all the
                eigenvalues, and no eigenvectors have been computed;
                elements i+1:N of WR and WI contain eigenvalues which
                have converged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 191 of file sgeev.f.

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 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:132
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:162
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:141
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
Definition: strevc.f:224
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:170
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53

Here is the call graph for this function:

Here is the caller graph for this function: