LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sgeev()

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
June 2016

Definition at line 193 of file sgeev.f.

193  implicit none
194 *
195 * -- LAPACK driver routine (version 3.7.0) --
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 *
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:132
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:94
subroutine strevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
STREVC3
Definition: strevc3.f:241
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
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:76
subroutine sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
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
Here is the call graph for this function:
Here is the caller graph for this function: