LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

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

Definition at line 190 of file sgeev.f.

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