LAPACK 3.12.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, slacpy, slartg,
228* ..
229* .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ISAMAX, ILAENV
232 REAL SLAMCH, SLANGE, SLAPY2, SNRM2, SROUNDUP_LWORK
233 EXTERNAL lsame, isamax, ilaenv, slamch, slange, slapy2,
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 = int( work(1) )
284 maxwrk = max( maxwrk, n + 1, n + hswork )
285 CALL strevc3( 'L', 'B', SELECT, n, a, lda,
286 $ vl, ldvl, vr, ldvr, n, nout,
287 $ work, -1, ierr )
288 lwork_trevc = int( work(1) )
289 maxwrk = max( maxwrk, n + lwork_trevc )
290 maxwrk = max( maxwrk, 4*n )
291 ELSE IF( wantvr ) THEN
292 minwrk = 4*n
293 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
294 $ 'SORGHR', ' ', n, 1, n, -1 ) )
295 CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
296 $ work, -1, info )
297 hswork = int( work(1) )
298 maxwrk = max( maxwrk, n + 1, n + hswork )
299 CALL strevc3( 'R', 'B', SELECT, n, a, lda,
300 $ vl, ldvl, vr, ldvr, n, nout,
301 $ work, -1, ierr )
302 lwork_trevc = int( work(1) )
303 maxwrk = max( maxwrk, n + lwork_trevc )
304 maxwrk = max( maxwrk, 4*n )
305 ELSE
306 minwrk = 3*n
307 CALL shseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
308 $ work, -1, info )
309 hswork = int( work(1) )
310 maxwrk = max( maxwrk, n + 1, n + hswork )
311 END IF
312 maxwrk = max( maxwrk, minwrk )
313 END IF
314 work( 1 ) = sroundup_lwork(maxwrk)
315*
316 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319 END IF
320*
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'SGEEV ', -info )
323 RETURN
324 ELSE IF( lquery ) THEN
325 RETURN
326 END IF
327*
328* Quick return if possible
329*
330 IF( n.EQ.0 )
331 $ RETURN
332*
333* Get machine constants
334*
335 eps = slamch( 'P' )
336 smlnum = slamch( 'S' )
337 bignum = one / smlnum
338 smlnum = sqrt( smlnum ) / eps
339 bignum = one / smlnum
340*
341* Scale A if max element outside range [SMLNUM,BIGNUM]
342*
343 anrm = slange( 'M', n, n, a, lda, dum )
344 scalea = .false.
345 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
346 scalea = .true.
347 cscale = smlnum
348 ELSE IF( anrm.GT.bignum ) THEN
349 scalea = .true.
350 cscale = bignum
351 END IF
352 IF( scalea )
353 $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
354*
355* Balance the matrix
356* (Workspace: need N)
357*
358 ibal = 1
359 CALL sgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
360*
361* Reduce to upper Hessenberg form
362* (Workspace: need 3*N, prefer 2*N+N*NB)
363*
364 itau = ibal + n
365 iwrk = itau + n
366 CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
367 $ lwork-iwrk+1, ierr )
368*
369 IF( wantvl ) THEN
370*
371* Want left eigenvectors
372* Copy Householder vectors to VL
373*
374 side = 'L'
375 CALL slacpy( 'L', n, n, a, lda, vl, ldvl )
376*
377* Generate orthogonal matrix in VL
378* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
379*
380 CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
381 $ lwork-iwrk+1, ierr )
382*
383* Perform QR iteration, accumulating Schur vectors in VL
384* (Workspace: need N+1, prefer N+HSWORK (see comments) )
385*
386 iwrk = itau
387 CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
388 $ work( iwrk ), lwork-iwrk+1, info )
389*
390 IF( wantvr ) THEN
391*
392* Want left and right eigenvectors
393* Copy Schur vectors to VR
394*
395 side = 'B'
396 CALL slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
397 END IF
398*
399 ELSE IF( wantvr ) THEN
400*
401* Want right eigenvectors
402* Copy Householder vectors to VR
403*
404 side = 'R'
405 CALL slacpy( 'L', n, n, a, lda, vr, ldvr )
406*
407* Generate orthogonal matrix in VR
408* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
409*
410 CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
411 $ lwork-iwrk+1, ierr )
412*
413* Perform QR iteration, accumulating Schur vectors in VR
414* (Workspace: need N+1, prefer N+HSWORK (see comments) )
415*
416 iwrk = itau
417 CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
418 $ work( iwrk ), lwork-iwrk+1, info )
419*
420 ELSE
421*
422* Compute eigenvalues only
423* (Workspace: need N+1, prefer N+HSWORK (see comments) )
424*
425 iwrk = itau
426 CALL shseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
427 $ work( iwrk ), lwork-iwrk+1, info )
428 END IF
429*
430* If INFO .NE. 0 from SHSEQR, then quit
431*
432 IF( info.NE.0 )
433 $ GO TO 50
434*
435 IF( wantvl .OR. wantvr ) THEN
436*
437* Compute left and/or right eigenvectors
438* (Workspace: need 4*N, prefer N + N + 2*N*NB)
439*
440 CALL strevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
441 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
442 END IF
443*
444 IF( wantvl ) THEN
445*
446* Undo balancing of left eigenvectors
447* (Workspace: need N)
448*
449 CALL sgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
450 $ ierr )
451*
452* Normalize left eigenvectors and make largest component real
453*
454 DO 20 i = 1, n
455 IF( wi( i ).EQ.zero ) THEN
456 scl = one / snrm2( n, vl( 1, i ), 1 )
457 CALL sscal( n, scl, vl( 1, i ), 1 )
458 ELSE IF( wi( i ).GT.zero ) THEN
459 scl = one / slapy2( snrm2( n, vl( 1, i ), 1 ),
460 $ snrm2( n, vl( 1, i+1 ), 1 ) )
461 CALL sscal( n, scl, vl( 1, i ), 1 )
462 CALL sscal( n, scl, vl( 1, i+1 ), 1 )
463 DO 10 k = 1, n
464 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
465 10 CONTINUE
466 k = isamax( n, work( iwrk ), 1 )
467 CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
468 CALL srot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
469 vl( k, i+1 ) = zero
470 END IF
471 20 CONTINUE
472 END IF
473*
474 IF( wantvr ) THEN
475*
476* Undo balancing of right eigenvectors
477* (Workspace: need N)
478*
479 CALL sgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
480 $ ierr )
481*
482* Normalize right eigenvectors and make largest component real
483*
484 DO 40 i = 1, n
485 IF( wi( i ).EQ.zero ) THEN
486 scl = one / snrm2( n, vr( 1, i ), 1 )
487 CALL sscal( n, scl, vr( 1, i ), 1 )
488 ELSE IF( wi( i ).GT.zero ) THEN
489 scl = one / slapy2( snrm2( n, vr( 1, i ), 1 ),
490 $ snrm2( n, vr( 1, i+1 ), 1 ) )
491 CALL sscal( n, scl, vr( 1, i ), 1 )
492 CALL sscal( n, scl, vr( 1, i+1 ), 1 )
493 DO 30 k = 1, n
494 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
495 30 CONTINUE
496 k = isamax( n, work( iwrk ), 1 )
497 CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
498 CALL srot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
499 vr( k, i+1 ) = zero
500 END IF
501 40 CONTINUE
502 END IF
503*
504* Undo scaling if necessary
505*
506 50 CONTINUE
507 IF( scalea ) THEN
508 CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
509 $ max( n-info, 1 ), ierr )
510 CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
511 $ max( n-info, 1 ), ierr )
512 IF( info.GT.0 ) THEN
513 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
514 $ ierr )
515 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
516 $ ierr )
517 END IF
518 END IF
519*
520 work( 1 ) = sroundup_lwork(maxwrk)
521 RETURN
522*
523* End of SGEEV
524*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
SGEBAK
Definition sgebak.f:130
subroutine sgebal(job, n, a, lda, ilo, ihi, scale, info)
SGEBAL
Definition sgebal.f:163
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
Definition sgehrd.f:167
subroutine shseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
SHSEQR
Definition shseqr.f:316
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine strevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
STREVC3
Definition strevc3.f:237
subroutine sorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
SORGHR
Definition sorghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: