OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
LAPACK » General Matrices » Eigenvalue

Functions

subroutine sgegs (jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, info)
  SGEGS computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine sgegv (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine sgees (jobvs, sort, select, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info)
  SGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine sgeesx (jobvs, sort, select, sense, n, a, lda, sdim, wr, wi, vs, ldvs, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
  SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
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
subroutine sgeevx (balanc, jobvl, jobvr, sense, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, iwork, info)
  SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine sgges (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
  SGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine sgges3 (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
  SGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)
subroutine sggesx (jobvsl, jobvsr, sort, selctg, sense, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
  SGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine sggev (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  SGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine sggev3 (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  SGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)
subroutine sggevx (balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
  SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Detailed Description

This is the group of real eigenvalue driver functions for GE matrices

Function Documentation

◆ sgees()

subroutine sgees ( character jobvs,
character sort,
external select,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer sdim,
real, dimension( * ) wr,
real, dimension( * ) wi,
real, dimension( ldvs, * ) vs,
integer ldvs,
real, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

SGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
!>
!> SGEES computes for an N-by-N real nonsymmetric matrix A, the
!> eigenvalues, the real Schur form T, and, optionally, the matrix of
!> Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
!>
!> Optionally, it also orders the eigenvalues on the diagonal of the
!> real Schur form so that selected eigenvalues are at the top left.
!> The leading columns of Z then form an orthonormal basis for the
!> invariant subspace corresponding to the selected eigenvalues.
!>
!> A matrix is in real Schur form if it is upper quasi-triangular with
!> 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
!> form
!>         [  a  b  ]
!>         [  c  a  ]
!>
!> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
!>
Parameters
[in]JOBVS
!>          JOBVS is CHARACTER*1
!>          = 'N': Schur vectors are not computed;
!>          = 'V': Schur vectors are computed.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the Schur form.
!>          = 'N': Eigenvalues are not ordered;
!>          = 'S': Eigenvalues are ordered (see SELECT).
!>
[in]SELECT
!>          SELECT is a LOGICAL FUNCTION of two REAL arguments
!>          SELECT must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'S', SELECT is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          If SORT = 'N', SELECT is not referenced.
!>          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
!>          SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
!>          conjugate pair of eigenvalues is selected, then both complex
!>          eigenvalues are selected.
!>          Note that a selected complex eigenvalue may no longer
!>          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
!>          ordering may change the value of complex eigenvalues
!>          (especially if the eigenvalue is ill-conditioned); in this
!>          case INFO is set to N+2 (see INFO below).
!>
[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 by its real Schur form T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>                         for which SELECT is true. (Complex conjugate
!>                         pairs for which SELECT is true for either
!>                         eigenvalue count as 2.)
!>
[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 in the same order
!>          that they appear on the diagonal of the output Schur form T.
!>          Complex conjugate pairs of eigenvalues will appear
!>          consecutively with the eigenvalue having the positive
!>          imaginary part first.
!>
[out]VS
!>          VS is REAL array, dimension (LDVS,N)
!>          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
!>          vectors.
!>          If JOBVS = 'N', VS is not referenced.
!>
[in]LDVS
!>          LDVS is INTEGER
!>          The leading dimension of the array VS.  LDVS >= 1; if
!>          JOBVS = 'V', LDVS >= N.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,3*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]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0: if INFO = i, and i is
!>             <= N: the QR algorithm failed to compute all the
!>                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
!>                   contain those eigenvalues which have converged; if
!>                   JOBVS = 'V', VS contains the matrix which reduces A
!>                   to its partially converged Schur form.
!>             = N+1: the eigenvalues could not be reordered because some
!>                   eigenvalues were too close to separate (the problem
!>                   is very ill-conditioned);
!>             = N+2: after reordering, roundoff changed values of some
!>                   complex eigenvalues so that leading eigenvalues in
!>                   the Schur form no longer satisfy SELECT=.TRUE.  This
!>                   could also be caused by underflow due to scaling.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 214 of file sgees.f.

216*
217* -- LAPACK driver routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 CHARACTER JOBVS, SORT
223 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
224* ..
225* .. Array Arguments ..
226 LOGICAL BWORK( * )
227 REAL A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
228 $ WR( * )
229* ..
230* .. Function Arguments ..
231 LOGICAL SELECT
232 EXTERNAL SELECT
233* ..
234*
235* =====================================================================
236*
237* .. Parameters ..
238 REAL ZERO, ONE
239 parameter( zero = 0.0e0, one = 1.0e0 )
240* ..
241* .. Local Scalars ..
242 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
243 $ WANTVS
244 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
245 $ IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
246 REAL ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
247* ..
248* .. Local Arrays ..
249 INTEGER IDUM( 1 )
250 REAL DUM( 1 )
251* ..
252* .. External Subroutines ..
253 EXTERNAL scopy, sgebak, sgebal, sgehrd, shseqr, slabad,
254 $ slacpy, slascl, sorghr, sswap, strsen, xerbla
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 INTEGER ILAENV
259 REAL SLAMCH, SLANGE
260 EXTERNAL lsame, ilaenv, slamch, slange
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC max, sqrt
264* ..
265* .. Executable Statements ..
266*
267* Test the input arguments
268*
269 info = 0
270 lquery = ( lwork.EQ.-1 )
271 wantvs = lsame( jobvs, 'V' )
272 wantst = lsame( sort, 'S' )
273 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
274 info = -1
275 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
276 info = -2
277 ELSE IF( n.LT.0 ) THEN
278 info = -4
279 ELSE IF( lda.LT.max( 1, n ) ) THEN
280 info = -6
281 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
282 info = -11
283 END IF
284*
285* Compute workspace
286* (Note: Comments in the code beginning "Workspace:" describe the
287* minimal amount of workspace needed at that point in the code,
288* as well as the preferred amount for good performance.
289* NB refers to the optimal block size for the immediately
290* following subroutine, as returned by ILAENV.
291* HSWORK refers to the workspace preferred by SHSEQR, as
292* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
293* the worst case.)
294*
295 IF( info.EQ.0 ) THEN
296 IF( n.EQ.0 ) THEN
297 minwrk = 1
298 maxwrk = 1
299 ELSE
300 maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
301 minwrk = 3*n
302*
303 CALL shseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
304 $ work, -1, ieval )
305 hswork = work( 1 )
306*
307 IF( .NOT.wantvs ) THEN
308 maxwrk = max( maxwrk, n + hswork )
309 ELSE
310 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
311 $ 'SORGHR', ' ', n, 1, n, -1 ) )
312 maxwrk = max( maxwrk, n + hswork )
313 END IF
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( 'SGEES ', -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 ) THEN
332 sdim = 0
333 RETURN
334 END IF
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* Permute the matrix to make it more nearly triangular
360* (Workspace: need N)
361*
362 ibal = 1
363 CALL sgebal( 'P', 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 = n + ibal
369 iwrk = n + itau
370 CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
371 $ lwork-iwrk+1, ierr )
372*
373 IF( wantvs ) THEN
374*
375* Copy Householder vectors to VS
376*
377 CALL slacpy( 'L', n, n, a, lda, vs, ldvs )
378*
379* Generate orthogonal matrix in VS
380* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381*
382 CALL sorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
383 $ lwork-iwrk+1, ierr )
384 END IF
385*
386 sdim = 0
387*
388* Perform QR iteration, accumulating Schur vectors in VS if desired
389* (Workspace: need N+1, prefer N+HSWORK (see comments) )
390*
391 iwrk = itau
392 CALL shseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
393 $ work( iwrk ), lwork-iwrk+1, ieval )
394 IF( ieval.GT.0 )
395 $ info = ieval
396*
397* Sort eigenvalues if desired
398*
399 IF( wantst .AND. info.EQ.0 ) THEN
400 IF( scalea ) THEN
401 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
402 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
403 END IF
404 DO 10 i = 1, n
405 bwork( i ) = SELECT( wr( i ), wi( i ) )
406 10 CONTINUE
407*
408* Reorder eigenvalues and transform Schur vectors
409* (Workspace: none needed)
410*
411 CALL strsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
412 $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
413 $ icond )
414 IF( icond.GT.0 )
415 $ info = n + icond
416 END IF
417*
418 IF( wantvs ) THEN
419*
420* Undo balancing
421* (Workspace: need N)
422*
423 CALL sgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
424 $ ierr )
425 END IF
426*
427 IF( scalea ) THEN
428*
429* Undo scaling for the Schur form of A
430*
431 CALL slascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
432 CALL scopy( n, a, lda+1, wr, 1 )
433 IF( cscale.EQ.smlnum ) THEN
434*
435* If scaling back towards underflow, adjust WI if an
436* offdiagonal element of a 2-by-2 block in the Schur form
437* underflows.
438*
439 IF( ieval.GT.0 ) THEN
440 i1 = ieval + 1
441 i2 = ihi - 1
442 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
443 $ max( ilo-1, 1 ), ierr )
444 ELSE IF( wantst ) THEN
445 i1 = 1
446 i2 = n - 1
447 ELSE
448 i1 = ilo
449 i2 = ihi - 1
450 END IF
451 inxt = i1 - 1
452 DO 20 i = i1, i2
453 IF( i.LT.inxt )
454 $ GO TO 20
455 IF( wi( i ).EQ.zero ) THEN
456 inxt = i + 1
457 ELSE
458 IF( a( i+1, i ).EQ.zero ) THEN
459 wi( i ) = zero
460 wi( i+1 ) = zero
461 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
462 $ zero ) THEN
463 wi( i ) = zero
464 wi( i+1 ) = zero
465 IF( i.GT.1 )
466 $ CALL sswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
467 IF( n.GT.i+1 )
468 $ CALL sswap( n-i-1, a( i, i+2 ), lda,
469 $ a( i+1, i+2 ), lda )
470 IF( wantvs ) THEN
471 CALL sswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
472 END IF
473 a( i, i+1 ) = a( i+1, i )
474 a( i+1, i ) = zero
475 END IF
476 inxt = i + 2
477 END IF
478 20 CONTINUE
479 END IF
480*
481* Undo scaling for the imaginary part of the eigenvalues
482*
483 CALL slascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
484 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
485 END IF
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488*
489* Check if reordering successful
490*
491 lastsl = .true.
492 lst2sl = .true.
493 sdim = 0
494 ip = 0
495 DO 30 i = 1, n
496 cursl = SELECT( wr( i ), wi( i ) )
497 IF( wi( i ).EQ.zero ) THEN
498 IF( cursl )
499 $ sdim = sdim + 1
500 ip = 0
501 IF( cursl .AND. .NOT.lastsl )
502 $ info = n + 2
503 ELSE
504 IF( ip.EQ.1 ) THEN
505*
506* Last eigenvalue of conjugate pair
507*
508 cursl = cursl .OR. lastsl
509 lastsl = cursl
510 IF( cursl )
511 $ sdim = sdim + 2
512 ip = -1
513 IF( cursl .AND. .NOT.lst2sl )
514 $ info = n + 2
515 ELSE
516*
517* First eigenvalue of conjugate pair
518*
519 ip = 1
520 END IF
521 END IF
522 lst2sl = lastsl
523 lastsl = cursl
524 30 CONTINUE
525 END IF
526*
527 work( 1 ) = maxwrk
528 RETURN
529*
530* End of SGEES
531*
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
slacpy
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
slascl
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slange
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
sgehrd
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
Definition sgehrd.f:167
sgebal
subroutine sgebal(job, n, a, lda, ilo, ihi, scale, info)
SGEBAL
Definition sgebal.f:160
sgebak
subroutine sgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
SGEBAK
Definition sgebak.f:130
shseqr
subroutine shseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
SHSEQR
Definition shseqr.f:316
strsen
subroutine strsen(job, compq, select, n, t, ldt, q, ldq, wr, wi, m, s, sep, work, lwork, iwork, liwork, info)
STRSEN
Definition strsen.f:314
sorghr
subroutine sorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
SORGHR
Definition sorghr.f:126
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
max
#define max(a, b)
Definition macros.h:21

◆ sgeesx()

subroutine sgeesx ( character jobvs,
character sort,
external select,
character sense,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer sdim,
real, dimension( * ) wr,
real, dimension( * ) wi,
real, dimension( ldvs, * ) vs,
integer ldvs,
real rconde,
real rcondv,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
logical, dimension( * ) bwork,
integer info )

SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
!>
!> SGEESX computes for an N-by-N real nonsymmetric matrix A, the
!> eigenvalues, the real Schur form T, and, optionally, the matrix of
!> Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
!>
!> Optionally, it also orders the eigenvalues on the diagonal of the
!> real Schur form so that selected eigenvalues are at the top left;
!> computes a reciprocal condition number for the average of the
!> selected eigenvalues (RCONDE); and computes a reciprocal condition
!> number for the right invariant subspace corresponding to the
!> selected eigenvalues (RCONDV).  The leading columns of Z form an
!> orthonormal basis for this invariant subspace.
!>
!> For further explanation of the reciprocal condition numbers RCONDE
!> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
!> these quantities are called s and sep respectively).
!>
!> A real matrix is in real Schur form if it is upper quasi-triangular
!> with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
!> the form
!>           [  a  b  ]
!>           [  c  a  ]
!>
!> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
!>
Parameters
[in]JOBVS
!>          JOBVS is CHARACTER*1
!>          = 'N': Schur vectors are not computed;
!>          = 'V': Schur vectors are computed.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the Schur form.
!>          = 'N': Eigenvalues are not ordered;
!>          = 'S': Eigenvalues are ordered (see SELECT).
!>
[in]SELECT
!>          SELECT is a LOGICAL FUNCTION of two REAL arguments
!>          SELECT must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'S', SELECT is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          If SORT = 'N', SELECT is not referenced.
!>          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
!>          SELECT(WR(j),WI(j)) is true; i.e., if either one of a
!>          complex conjugate pair of eigenvalues is selected, then both
!>          are.  Note that a selected complex eigenvalue may no longer
!>          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
!>          ordering may change the value of complex eigenvalues
!>          (especially if the eigenvalue is ill-conditioned); in this
!>          case INFO may be set to N+3 (see INFO below).
!>
[in]SENSE
!>          SENSE is CHARACTER*1
!>          Determines which reciprocal condition numbers are computed.
!>          = 'N': None are computed;
!>          = 'E': Computed for average of selected eigenvalues only;
!>          = 'V': Computed for selected right invariant subspace only;
!>          = 'B': Computed for both.
!>          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
!>
[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 is overwritten by its real Schur form T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>                         for which SELECT is true. (Complex conjugate
!>                         pairs for which SELECT is true for either
!>                         eigenvalue count as 2.)
!>
[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, in the same order that they
!>          appear on the diagonal of the output Schur form T.  Complex
!>          conjugate pairs of eigenvalues appear consecutively with the
!>          eigenvalue having the positive imaginary part first.
!>
[out]VS
!>          VS is REAL array, dimension (LDVS,N)
!>          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
!>          vectors.
!>          If JOBVS = 'N', VS is not referenced.
!>
[in]LDVS
!>          LDVS is INTEGER
!>          The leading dimension of the array VS.  LDVS >= 1, and if
!>          JOBVS = 'V', LDVS >= N.
!>
[out]RCONDE
!>          RCONDE is REAL
!>          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
!>          condition number for the average of the selected eigenvalues.
!>          Not referenced if SENSE = 'N' or 'V'.
!>
[out]RCONDV
!>          RCONDV is REAL
!>          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
!>          condition number for the selected right invariant subspace.
!>          Not referenced if SENSE = 'N' or 'E'.
!>
[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).
!>          Also, if SENSE = 'E' or 'V' or 'B',
!>          LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
!>          selected eigenvalues computed by this routine.  Note that
!>          N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
!>          returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
!>          'B' this may not be large enough.
!>          For good performance, LWORK must generally be larger.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates upper bounds on the optimal sizes of the
!>          arrays WORK and IWORK, returns these values as the first
!>          entries of the WORK and IWORK arrays, and no error messages
!>          related to LWORK or LIWORK are issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
!>          Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
!>          only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
!>          may not be large enough.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates upper bounds on the optimal sizes of
!>          the arrays WORK and IWORK, returns these values as the first
!>          entries of the WORK and IWORK arrays, and no error messages
!>          related to LWORK or LIWORK are issued by XERBLA.
!>
[out]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0: if INFO = i, and i is
!>             <= N: the QR algorithm failed to compute all the
!>                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
!>                   contain those eigenvalues which have converged; if
!>                   JOBVS = 'V', VS contains the transformation which
!>                   reduces A to its partially converged Schur form.
!>             = N+1: the eigenvalues could not be reordered because some
!>                   eigenvalues were too close to separate (the problem
!>                   is very ill-conditioned);
!>             = N+2: after reordering, roundoff changed values of some
!>                   complex eigenvalues so that leading eigenvalues in
!>                   the Schur form no longer satisfy SELECT=.TRUE.  This
!>                   could also be caused by underflow due to scaling.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 278 of file sgeesx.f.

281*
282* -- LAPACK driver routine --
283* -- LAPACK is a software package provided by Univ. of Tennessee, --
284* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
285*
286* .. Scalar Arguments ..
287 CHARACTER JOBVS, SENSE, SORT
288 INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
289 REAL RCONDE, RCONDV
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 INTEGER IWORK( * )
294 REAL A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
295 $ WR( * )
296* ..
297* .. Function Arguments ..
298 LOGICAL SELECT
299 EXTERNAL SELECT
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 REAL ZERO, ONE
306 parameter( zero = 0.0e0, one = 1.0e0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
310 $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
311 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
312 $ IHI, ILO, INXT, IP, ITAU, IWRK, LWRK, LIWRK,
313 $ MAXWRK, MINWRK
314 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315* ..
316* .. Local Arrays ..
317 REAL DUM( 1 )
318* ..
319* .. External Subroutines ..
320 EXTERNAL scopy, sgebak, sgebal, sgehrd, shseqr, slabad,
321 $ slacpy, slascl, sorghr, sswap, strsen, xerbla
322* ..
323* .. External Functions ..
324 LOGICAL LSAME
325 INTEGER ILAENV
326 REAL SLAMCH, SLANGE
327 EXTERNAL lsame, ilaenv, slamch, slange
328* ..
329* .. Intrinsic Functions ..
330 INTRINSIC max, sqrt
331* ..
332* .. Executable Statements ..
333*
334* Test the input arguments
335*
336 info = 0
337 wantvs = lsame( jobvs, 'V' )
338 wantst = lsame( sort, 'S' )
339 wantsn = lsame( sense, 'N' )
340 wantse = lsame( sense, 'E' )
341 wantsv = lsame( sense, 'V' )
342 wantsb = lsame( sense, 'B' )
343 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
344*
345 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
346 info = -1
347 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
348 info = -2
349 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
350 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
351 info = -4
352 ELSE IF( n.LT.0 ) THEN
353 info = -5
354 ELSE IF( lda.LT.max( 1, n ) ) THEN
355 info = -7
356 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
357 info = -12
358 END IF
359*
360* Compute workspace
361* (Note: Comments in the code beginning "RWorkspace:" describe the
362* minimal amount of real workspace needed at that point in the
363* code, as well as the preferred amount for good performance.
364* IWorkspace refers to integer workspace.
365* NB refers to the optimal block size for the immediately
366* following subroutine, as returned by ILAENV.
367* HSWORK refers to the workspace preferred by SHSEQR, as
368* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
369* the worst case.
370* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
371* depends on SDIM, which is computed by the routine STRSEN later
372* in the code.)
373*
374 IF( info.EQ.0 ) THEN
375 liwrk = 1
376 IF( n.EQ.0 ) THEN
377 minwrk = 1
378 lwrk = 1
379 ELSE
380 maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
381 minwrk = 3*n
382*
383 CALL shseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
384 $ work, -1, ieval )
385 hswork = work( 1 )
386*
387 IF( .NOT.wantvs ) THEN
388 maxwrk = max( maxwrk, n + hswork )
389 ELSE
390 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
391 $ 'SORGHR', ' ', n, 1, n, -1 ) )
392 maxwrk = max( maxwrk, n + hswork )
393 END IF
394 lwrk = maxwrk
395 IF( .NOT.wantsn )
396 $ lwrk = max( lwrk, n + ( n*n )/2 )
397 IF( wantsv .OR. wantsb )
398 $ liwrk = ( n*n )/4
399 END IF
400 iwork( 1 ) = liwrk
401 work( 1 ) = lwrk
402*
403 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
404 info = -16
405 ELSE IF( liwork.LT.1 .AND. .NOT.lquery ) THEN
406 info = -18
407 END IF
408 END IF
409*
410 IF( info.NE.0 ) THEN
411 CALL xerbla( 'SGEESX', -info )
412 RETURN
413 ELSE IF( lquery ) THEN
414 RETURN
415 END IF
416*
417* Quick return if possible
418*
419 IF( n.EQ.0 ) THEN
420 sdim = 0
421 RETURN
422 END IF
423*
424* Get machine constants
425*
426 eps = slamch( 'P' )
427 smlnum = slamch( 'S' )
428 bignum = one / smlnum
429 CALL slabad( smlnum, bignum )
430 smlnum = sqrt( smlnum ) / eps
431 bignum = one / smlnum
432*
433* Scale A if max element outside range [SMLNUM,BIGNUM]
434*
435 anrm = slange( 'M', n, n, a, lda, dum )
436 scalea = .false.
437 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
438 scalea = .true.
439 cscale = smlnum
440 ELSE IF( anrm.GT.bignum ) THEN
441 scalea = .true.
442 cscale = bignum
443 END IF
444 IF( scalea )
445 $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
446*
447* Permute the matrix to make it more nearly triangular
448* (RWorkspace: need N)
449*
450 ibal = 1
451 CALL sgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
452*
453* Reduce to upper Hessenberg form
454* (RWorkspace: need 3*N, prefer 2*N+N*NB)
455*
456 itau = n + ibal
457 iwrk = n + itau
458 CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
459 $ lwork-iwrk+1, ierr )
460*
461 IF( wantvs ) THEN
462*
463* Copy Householder vectors to VS
464*
465 CALL slacpy( 'L', n, n, a, lda, vs, ldvs )
466*
467* Generate orthogonal matrix in VS
468* (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
469*
470 CALL sorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
471 $ lwork-iwrk+1, ierr )
472 END IF
473*
474 sdim = 0
475*
476* Perform QR iteration, accumulating Schur vectors in VS if desired
477* (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
478*
479 iwrk = itau
480 CALL shseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
481 $ work( iwrk ), lwork-iwrk+1, ieval )
482 IF( ieval.GT.0 )
483 $ info = ieval
484*
485* Sort eigenvalues if desired
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488 IF( scalea ) THEN
489 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
490 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
491 END IF
492 DO 10 i = 1, n
493 bwork( i ) = SELECT( wr( i ), wi( i ) )
494 10 CONTINUE
495*
496* Reorder eigenvalues, transform Schur vectors, and compute
497* reciprocal condition numbers
498* (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
499* otherwise, need N )
500* (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
501* otherwise, need 0 )
502*
503 CALL strsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
504 $ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
505 $ iwork, liwork, icond )
506 IF( .NOT.wantsn )
507 $ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
508 IF( icond.EQ.-15 ) THEN
509*
510* Not enough real workspace
511*
512 info = -16
513 ELSE IF( icond.EQ.-17 ) THEN
514*
515* Not enough integer workspace
516*
517 info = -18
518 ELSE IF( icond.GT.0 ) THEN
519*
520* STRSEN failed to reorder or to restore standard Schur form
521*
522 info = icond + n
523 END IF
524 END IF
525*
526 IF( wantvs ) THEN
527*
528* Undo balancing
529* (RWorkspace: need N)
530*
531 CALL sgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
532 $ ierr )
533 END IF
534*
535 IF( scalea ) THEN
536*
537* Undo scaling for the Schur form of A
538*
539 CALL slascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
540 CALL scopy( n, a, lda+1, wr, 1 )
541 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
542 dum( 1 ) = rcondv
543 CALL slascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
544 rcondv = dum( 1 )
545 END IF
546 IF( cscale.EQ.smlnum ) THEN
547*
548* If scaling back towards underflow, adjust WI if an
549* offdiagonal element of a 2-by-2 block in the Schur form
550* underflows.
551*
552 IF( ieval.GT.0 ) THEN
553 i1 = ieval + 1
554 i2 = ihi - 1
555 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
556 $ ierr )
557 ELSE IF( wantst ) THEN
558 i1 = 1
559 i2 = n - 1
560 ELSE
561 i1 = ilo
562 i2 = ihi - 1
563 END IF
564 inxt = i1 - 1
565 DO 20 i = i1, i2
566 IF( i.LT.inxt )
567 $ GO TO 20
568 IF( wi( i ).EQ.zero ) THEN
569 inxt = i + 1
570 ELSE
571 IF( a( i+1, i ).EQ.zero ) THEN
572 wi( i ) = zero
573 wi( i+1 ) = zero
574 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
575 $ zero ) THEN
576 wi( i ) = zero
577 wi( i+1 ) = zero
578 IF( i.GT.1 )
579 $ CALL sswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
580 IF( n.GT.i+1 )
581 $ CALL sswap( n-i-1, a( i, i+2 ), lda,
582 $ a( i+1, i+2 ), lda )
583 IF( wantvs ) THEN
584 CALL sswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
585 END IF
586 a( i, i+1 ) = a( i+1, i )
587 a( i+1, i ) = zero
588 END IF
589 inxt = i + 2
590 END IF
591 20 CONTINUE
592 END IF
593 CALL slascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
594 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
595 END IF
596*
597 IF( wantst .AND. info.EQ.0 ) THEN
598*
599* Check if reordering successful
600*
601 lastsl = .true.
602 lst2sl = .true.
603 sdim = 0
604 ip = 0
605 DO 30 i = 1, n
606 cursl = SELECT( wr( i ), wi( i ) )
607 IF( wi( i ).EQ.zero ) THEN
608 IF( cursl )
609 $ sdim = sdim + 1
610 ip = 0
611 IF( cursl .AND. .NOT.lastsl )
612 $ info = n + 2
613 ELSE
614 IF( ip.EQ.1 ) THEN
615*
616* Last eigenvalue of conjugate pair
617*
618 cursl = cursl .OR. lastsl
619 lastsl = cursl
620 IF( cursl )
621 $ sdim = sdim + 2
622 ip = -1
623 IF( cursl .AND. .NOT.lst2sl )
624 $ info = n + 2
625 ELSE
626*
627* First eigenvalue of conjugate pair
628*
629 ip = 1
630 END IF
631 END IF
632 lst2sl = lastsl
633 lastsl = cursl
634 30 CONTINUE
635 END IF
636*
637 work( 1 ) = maxwrk
638 IF( wantsv .OR. wantsb ) THEN
639 iwork( 1 ) = sdim*(n-sdim)
640 ELSE
641 iwork( 1 ) = 1
642 END IF
643*
644 RETURN
645*
646* End of SGEESX
647*

◆ 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,
227 $ slartg, slascl, sorghr, srot, sscal, strevc3,
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*
slapy2
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:113
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
strevc3
subroutine strevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
STREVC3
Definition strevc3.f:237
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89

◆ sgeevx()

subroutine sgeevx ( character balanc,
character jobvl,
character jobvr,
character sense,
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,
integer ilo,
integer ihi,
real, dimension( * ) scale,
real abnrm,
real, dimension( * ) rconde,
real, dimension( * ) rcondv,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer info )

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

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

Purpose:
!>
!> SGEEVX computes for an N-by-N real nonsymmetric matrix A, the
!> eigenvalues and, optionally, the left and/or right eigenvectors.
!>
!> Optionally also, it computes a balancing transformation to improve
!> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
!> SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
!> (RCONDE), and reciprocal condition numbers for the right
!> eigenvectors (RCONDV).
!>
!> 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.
!>
!> Balancing a matrix means permuting the rows and columns to make it
!> more nearly upper triangular, and applying a diagonal similarity
!> transformation D * A * D**(-1), where D is a diagonal matrix, to
!> make its rows and columns closer in norm and the condition numbers
!> of its eigenvalues and eigenvectors smaller.  The computed
!> reciprocal condition numbers correspond to the balanced matrix.
!> Permuting rows and columns will not change the condition numbers
!> (in exact arithmetic) but diagonal scaling will.  For further
!> explanation of balancing, see section 4.10.2 of the LAPACK
!> Users' Guide.
!>
Parameters
[in]BALANC
!>          BALANC is CHARACTER*1
!>          Indicates how the input matrix should be diagonally scaled
!>          and/or permuted to improve the conditioning of its
!>          eigenvalues.
!>          = 'N': Do not diagonally scale or permute;
!>          = 'P': Perform permutations to make the matrix more nearly
!>                 upper triangular. Do not diagonally scale;
!>          = 'S': Diagonally scale the matrix, i.e. replace A by
!>                 D*A*D**(-1), where D is a diagonal matrix chosen
!>                 to make the rows and columns of A more equal in
!>                 norm. Do not permute;
!>          = 'B': Both diagonally scale and permute A.
!>
!>          Computed reciprocal condition numbers will be for the matrix
!>          after balancing and/or permuting. Permuting does not change
!>          condition numbers (in exact arithmetic), but balancing does.
!>
[in]JOBVL
!>          JOBVL is CHARACTER*1
!>          = 'N': left eigenvectors of A are not computed;
!>          = 'V': left eigenvectors of A are computed.
!>          If SENSE = 'E' or 'B', JOBVL must = 'V'.
!>
[in]JOBVR
!>          JOBVR is CHARACTER*1
!>          = 'N': right eigenvectors of A are not computed;
!>          = 'V': right eigenvectors of A are computed.
!>          If SENSE = 'E' or 'B', JOBVR must = 'V'.
!>
[in]SENSE
!>          SENSE is CHARACTER*1
!>          Determines which reciprocal condition numbers are computed.
!>          = 'N': None are computed;
!>          = 'E': Computed for eigenvalues only;
!>          = 'V': Computed for right eigenvectors only;
!>          = 'B': Computed for eigenvalues and right eigenvectors.
!>
!>          If SENSE = 'E' or 'B', both left and right eigenvectors
!>          must also be computed (JOBVL = 'V' and JOBVR = 'V').
!>
[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.  If JOBVL = 'V' or
!>          JOBVR = 'V', A contains the real Schur form of the balanced
!>          version of the input matrix A.
!>
[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 will 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, and if
!>          JOBVR = 'V', LDVR >= N.
!>
[out]ILO
!>          ILO is INTEGER
!>
[out]IHI
!>          IHI is INTEGER
!>          ILO and IHI are integer values determined when A was
!>          balanced.  The balanced A(i,j) = 0 if I > J and
!>          J = 1,...,ILO-1 or I = IHI+1,...,N.
!>
[out]SCALE
!>          SCALE is REAL array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          when balancing A.  If P(j) is the index of the row and column
!>          interchanged with row and column j, and D(j) is the scaling
!>          factor applied to row and column j, then
!>          SCALE(J) = P(J),    for J = 1,...,ILO-1
!>                   = D(J),    for J = ILO,...,IHI
!>                   = P(J)     for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]ABNRM
!>          ABNRM is REAL
!>          The one-norm of the balanced matrix (the maximum
!>          of the sum of absolute values of elements of any column).
!>
[out]RCONDE
!>          RCONDE is REAL array, dimension (N)
!>          RCONDE(j) is the reciprocal condition number of the j-th
!>          eigenvalue.
!>
[out]RCONDV
!>          RCONDV is REAL array, dimension (N)
!>          RCONDV(j) is the reciprocal condition number of the j-th
!>          right eigenvector.
!>
[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.   If SENSE = 'N' or 'E',
!>          LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
!>          LWORK >= 3*N.  If SENSE = 'V' or 'B', LWORK >= N*(N+6).
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (2*N-2)
!>          If SENSE = 'N' or 'E', not referenced.
!>
[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 or condition numbers
!>                have been computed; elements 1:ILO-1 and 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 303 of file sgeevx.f.

306 implicit none
307*
308* -- LAPACK driver routine --
309* -- LAPACK is a software package provided by Univ. of Tennessee, --
310* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
311*
312* .. Scalar Arguments ..
313 CHARACTER BALANC, JOBVL, JOBVR, SENSE
314 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
315 REAL ABNRM
316* ..
317* .. Array Arguments ..
318 INTEGER IWORK( * )
319 REAL A( LDA, * ), RCONDE( * ), RCONDV( * ),
320 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
321 $ WI( * ), WORK( * ), WR( * )
322* ..
323*
324* =====================================================================
325*
326* .. Parameters ..
327 REAL ZERO, ONE
328 parameter( zero = 0.0e0, one = 1.0e0 )
329* ..
330* .. Local Scalars ..
331 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
332 $ WNTSNN, WNTSNV
333 CHARACTER JOB, SIDE
334 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
335 $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
336 REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
337 $ SN
338* ..
339* .. Local Arrays ..
340 LOGICAL SELECT( 1 )
341 REAL DUM( 1 )
342* ..
343* .. External Subroutines ..
344 EXTERNAL sgebak, sgebal, sgehrd, shseqr, slabad, slacpy,
345 $ slartg, slascl, sorghr, srot, sscal, strevc3,
346 $ strsna, xerbla
347* ..
348* .. External Functions ..
349 LOGICAL LSAME
350 INTEGER ISAMAX, ILAENV
351 REAL SLAMCH, SLANGE, SLAPY2, SNRM2
352 EXTERNAL lsame, isamax, ilaenv, slamch, slange, slapy2,
353 $ snrm2
354* ..
355* .. Intrinsic Functions ..
356 INTRINSIC max, sqrt
357* ..
358* .. Executable Statements ..
359*
360* Test the input arguments
361*
362 info = 0
363 lquery = ( lwork.EQ.-1 )
364 wantvl = lsame( jobvl, 'V' )
365 wantvr = lsame( jobvr, 'V' )
366 wntsnn = lsame( sense, 'N' )
367 wntsne = lsame( sense, 'E' )
368 wntsnv = lsame( sense, 'V' )
369 wntsnb = lsame( sense, 'B' )
370 IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc, 'S' )
371 $ .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
372 $ THEN
373 info = -1
374 ELSE IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
375 info = -2
376 ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
377 info = -3
378 ELSE IF( .NOT.( wntsnn .OR. wntsne .OR. wntsnb .OR. wntsnv ) .OR.
379 $ ( ( wntsne .OR. wntsnb ) .AND. .NOT.( wantvl .AND.
380 $ wantvr ) ) ) THEN
381 info = -4
382 ELSE IF( n.LT.0 ) THEN
383 info = -5
384 ELSE IF( lda.LT.max( 1, n ) ) THEN
385 info = -7
386 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
387 info = -11
388 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
389 info = -13
390 END IF
391*
392* Compute workspace
393* (Note: Comments in the code beginning "Workspace:" describe the
394* minimal amount of workspace needed at that point in the code,
395* as well as the preferred amount for good performance.
396* NB refers to the optimal block size for the immediately
397* following subroutine, as returned by ILAENV.
398* HSWORK refers to the workspace preferred by SHSEQR, as
399* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
400* the worst case.)
401*
402 IF( info.EQ.0 ) THEN
403 IF( n.EQ.0 ) THEN
404 minwrk = 1
405 maxwrk = 1
406 ELSE
407 maxwrk = n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
408*
409 IF( wantvl ) THEN
410 CALL strevc3( 'L', 'B', SELECT, n, a, lda,
411 $ vl, ldvl, vr, ldvr,
412 $ n, nout, work, -1, ierr )
413 lwork_trevc = int( work(1) )
414 maxwrk = max( maxwrk, n + lwork_trevc )
415 CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
416 $ work, -1, info )
417 ELSE IF( wantvr ) THEN
418 CALL strevc3( 'R', 'B', SELECT, n, a, lda,
419 $ vl, ldvl, vr, ldvr,
420 $ n, nout, work, -1, ierr )
421 lwork_trevc = int( work(1) )
422 maxwrk = max( maxwrk, n + lwork_trevc )
423 CALL shseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
424 $ work, -1, info )
425 ELSE
426 IF( wntsnn ) THEN
427 CALL shseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr,
428 $ ldvr, work, -1, info )
429 ELSE
430 CALL shseqr( 'S', 'N', n, 1, n, a, lda, wr, wi, vr,
431 $ ldvr, work, -1, info )
432 END IF
433 END IF
434 hswork = int( work(1) )
435*
436 IF( ( .NOT.wantvl ) .AND. ( .NOT.wantvr ) ) THEN
437 minwrk = 2*n
438 IF( .NOT.wntsnn )
439 $ minwrk = max( minwrk, n*n+6*n )
440 maxwrk = max( maxwrk, hswork )
441 IF( .NOT.wntsnn )
442 $ maxwrk = max( maxwrk, n*n + 6*n )
443 ELSE
444 minwrk = 3*n
445 IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
446 $ minwrk = max( minwrk, n*n + 6*n )
447 maxwrk = max( maxwrk, hswork )
448 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'SORGHR',
449 $ ' ', n, 1, n, -1 ) )
450 IF( ( .NOT.wntsnn ) .AND. ( .NOT.wntsne ) )
451 $ maxwrk = max( maxwrk, n*n + 6*n )
452 maxwrk = max( maxwrk, 3*n )
453 END IF
454 maxwrk = max( maxwrk, minwrk )
455 END IF
456 work( 1 ) = maxwrk
457*
458 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
459 info = -21
460 END IF
461 END IF
462*
463 IF( info.NE.0 ) THEN
464 CALL xerbla( 'SGEEVX', -info )
465 RETURN
466 ELSE IF( lquery ) THEN
467 RETURN
468 END IF
469*
470* Quick return if possible
471*
472 IF( n.EQ.0 )
473 $ RETURN
474*
475* Get machine constants
476*
477 eps = slamch( 'P' )
478 smlnum = slamch( 'S' )
479 bignum = one / smlnum
480 CALL slabad( smlnum, bignum )
481 smlnum = sqrt( smlnum ) / eps
482 bignum = one / smlnum
483*
484* Scale A if max element outside range [SMLNUM,BIGNUM]
485*
486 icond = 0
487 anrm = slange( 'M', n, n, a, lda, dum )
488 scalea = .false.
489 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
490 scalea = .true.
491 cscale = smlnum
492 ELSE IF( anrm.GT.bignum ) THEN
493 scalea = .true.
494 cscale = bignum
495 END IF
496 IF( scalea )
497 $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
498*
499* Balance the matrix and compute ABNRM
500*
501 CALL sgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
502 abnrm = slange( '1', n, n, a, lda, dum )
503 IF( scalea ) THEN
504 dum( 1 ) = abnrm
505 CALL slascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
506 abnrm = dum( 1 )
507 END IF
508*
509* Reduce to upper Hessenberg form
510* (Workspace: need 2*N, prefer N+N*NB)
511*
512 itau = 1
513 iwrk = itau + n
514 CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
515 $ lwork-iwrk+1, ierr )
516*
517 IF( wantvl ) THEN
518*
519* Want left eigenvectors
520* Copy Householder vectors to VL
521*
522 side = 'L'
523 CALL slacpy( 'L', n, n, a, lda, vl, ldvl )
524*
525* Generate orthogonal matrix in VL
526* (Workspace: need 2*N-1, prefer N+(N-1)*NB)
527*
528 CALL sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
529 $ lwork-iwrk+1, ierr )
530*
531* Perform QR iteration, accumulating Schur vectors in VL
532* (Workspace: need 1, prefer HSWORK (see comments) )
533*
534 iwrk = itau
535 CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
536 $ work( iwrk ), lwork-iwrk+1, info )
537*
538 IF( wantvr ) THEN
539*
540* Want left and right eigenvectors
541* Copy Schur vectors to VR
542*
543 side = 'B'
544 CALL slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
545 END IF
546*
547 ELSE IF( wantvr ) THEN
548*
549* Want right eigenvectors
550* Copy Householder vectors to VR
551*
552 side = 'R'
553 CALL slacpy( 'L', n, n, a, lda, vr, ldvr )
554*
555* Generate orthogonal matrix in VR
556* (Workspace: need 2*N-1, prefer N+(N-1)*NB)
557*
558 CALL sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
559 $ lwork-iwrk+1, ierr )
560*
561* Perform QR iteration, accumulating Schur vectors in VR
562* (Workspace: need 1, prefer HSWORK (see comments) )
563*
564 iwrk = itau
565 CALL shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
566 $ work( iwrk ), lwork-iwrk+1, info )
567*
568 ELSE
569*
570* Compute eigenvalues only
571* If condition numbers desired, compute Schur form
572*
573 IF( wntsnn ) THEN
574 job = 'E'
575 ELSE
576 job = 'S'
577 END IF
578*
579* (Workspace: need 1, prefer HSWORK (see comments) )
580*
581 iwrk = itau
582 CALL shseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
583 $ work( iwrk ), lwork-iwrk+1, info )
584 END IF
585*
586* If INFO .NE. 0 from SHSEQR, then quit
587*
588 IF( info.NE.0 )
589 $ GO TO 50
590*
591 IF( wantvl .OR. wantvr ) THEN
592*
593* Compute left and/or right eigenvectors
594* (Workspace: need 3*N, prefer N + 2*N*NB)
595*
596 CALL strevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
597 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
598 END IF
599*
600* Compute condition numbers if desired
601* (Workspace: need N*N+6*N unless SENSE = 'E')
602*
603 IF( .NOT.wntsnn ) THEN
604 CALL strsna( sense, 'A', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
605 $ rconde, rcondv, n, nout, work( iwrk ), n, iwork,
606 $ icond )
607 END IF
608*
609 IF( wantvl ) THEN
610*
611* Undo balancing of left eigenvectors
612*
613 CALL sgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,
614 $ ierr )
615*
616* Normalize left eigenvectors and make largest component real
617*
618 DO 20 i = 1, n
619 IF( wi( i ).EQ.zero ) THEN
620 scl = one / snrm2( n, vl( 1, i ), 1 )
621 CALL sscal( n, scl, vl( 1, i ), 1 )
622 ELSE IF( wi( i ).GT.zero ) THEN
623 scl = one / slapy2( snrm2( n, vl( 1, i ), 1 ),
624 $ snrm2( n, vl( 1, i+1 ), 1 ) )
625 CALL sscal( n, scl, vl( 1, i ), 1 )
626 CALL sscal( n, scl, vl( 1, i+1 ), 1 )
627 DO 10 k = 1, n
628 work( k ) = vl( k, i )**2 + vl( k, i+1 )**2
629 10 CONTINUE
630 k = isamax( n, work, 1 )
631 CALL slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
632 CALL srot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
633 vl( k, i+1 ) = zero
634 END IF
635 20 CONTINUE
636 END IF
637*
638 IF( wantvr ) THEN
639*
640* Undo balancing of right eigenvectors
641*
642 CALL sgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,
643 $ ierr )
644*
645* Normalize right eigenvectors and make largest component real
646*
647 DO 40 i = 1, n
648 IF( wi( i ).EQ.zero ) THEN
649 scl = one / snrm2( n, vr( 1, i ), 1 )
650 CALL sscal( n, scl, vr( 1, i ), 1 )
651 ELSE IF( wi( i ).GT.zero ) THEN
652 scl = one / slapy2( snrm2( n, vr( 1, i ), 1 ),
653 $ snrm2( n, vr( 1, i+1 ), 1 ) )
654 CALL sscal( n, scl, vr( 1, i ), 1 )
655 CALL sscal( n, scl, vr( 1, i+1 ), 1 )
656 DO 30 k = 1, n
657 work( k ) = vr( k, i )**2 + vr( k, i+1 )**2
658 30 CONTINUE
659 k = isamax( n, work, 1 )
660 CALL slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
661 CALL srot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
662 vr( k, i+1 ) = zero
663 END IF
664 40 CONTINUE
665 END IF
666*
667* Undo scaling if necessary
668*
669 50 CONTINUE
670 IF( scalea ) THEN
671 CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
672 $ max( n-info, 1 ), ierr )
673 CALL slascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
674 $ max( n-info, 1 ), ierr )
675 IF( info.EQ.0 ) THEN
676 IF( ( wntsnv .OR. wntsnb ) .AND. icond.EQ.0 )
677 $ CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
678 $ ierr )
679 ELSE
680 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
681 $ ierr )
682 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
683 $ ierr )
684 END IF
685 END IF
686*
687 work( 1 ) = maxwrk
688 RETURN
689*
690* End of SGEEVX
691*
strsna
subroutine strsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, iwork, info)
STRSNA
Definition strsna.f:265

◆ sgegs()

subroutine sgegs ( character jobvsl,
character jobvsr,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvsl, * ) vsl,
integer ldvsl,
real, dimension( ldvsr, * ) vsr,
integer ldvsr,
real, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine SGGES.
!>
!> SGEGS computes the eigenvalues, real Schur form, and, optionally,
!> left and or/right Schur vectors of a real matrix pair (A,B).
!> Given two square matrices A and B, the generalized real Schur
!> factorization has the form
!>
!>   A = Q*S*Z**T,  B = Q*T*Z**T
!>
!> where Q and Z are orthogonal matrices, T is upper triangular, and S
!> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
!> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
!> of eigenvalues of (A,B).  The columns of Q are the left Schur vectors
!> and the columns of Z are the right Schur vectors.
!>
!> If only the eigenvalues of (A,B) are needed, the driver routine
!> SGEGV should be used instead.  See SGEGV for a description of the
!> eigenvalues of the generalized nonsymmetric eigenvalue problem
!> (GNEP).
!>
Parameters
[in]JOBVSL
!>          JOBVSL is CHARACTER*1
!>          = 'N':  do not compute the left Schur vectors;
!>          = 'V':  compute the left Schur vectors (returned in VSL).
!>
[in]JOBVSR
!>          JOBVSR is CHARACTER*1
!>          = 'N':  do not compute the right Schur vectors;
!>          = 'V':  compute the right Schur vectors (returned in VSR).
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the matrix A.
!>          On exit, the upper quasi-triangular matrix S from the
!>          generalized real Schur factorization.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the matrix B.
!>          On exit, the upper triangular matrix T from the generalized
!>          real Schur factorization.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>          The real parts of each scalar alpha defining an eigenvalue
!>          of GNEP.
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>          The imaginary parts of each scalar alpha defining an
!>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
!>          eigenvalue is real; if positive, then the j-th and (j+1)-st
!>          eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) = -ALPHAI(j).
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
!>          beta = BETA(j) represent the j-th eigenvalue of the matrix
!>          pair (A,B), in one of the forms lambda = alpha/beta or
!>          mu = beta/alpha.  Since either lambda or mu may overflow,
!>          they should not, in general, be computed.
!>
[out]VSL
!>          VSL is REAL array, dimension (LDVSL,N)
!>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
!>          Not referenced if JOBVSL = 'N'.
!>
[in]LDVSL
!>          LDVSL is INTEGER
!>          The leading dimension of the matrix VSL. LDVSL >=1, and
!>          if JOBVSL = 'V', LDVSL >= N.
!>
[out]VSR
!>          VSR is REAL array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
!>          Not referenced if JOBVSR = 'N'.
!>
[in]LDVSR
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= 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,4*N).
!>          For good performance, LWORK must generally be larger.
!>          To compute the optimal value of LWORK, call ILAENV to get
!>          blocksizes (for SGEQRF, SORMQR, and SORGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR
!>          The optimal LWORK is  2*N + N*(NB+1).
!>
!>          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.
!>          = 1,...,N:
!>                The QZ iteration failed.  (A,B) are not in Schur
!>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  errors that usually indicate LAPACK problems:
!>                =N+1: error return from SGGBAL
!>                =N+2: error return from SGEQRF
!>                =N+3: error return from SORMQR
!>                =N+4: error return from SORGQR
!>                =N+5: error return from SGGHRD
!>                =N+6: error return from SHGEQZ (other than failed
!>                                                iteration)
!>                =N+7: error return from SGGBAK (computing VSL)
!>                =N+8: error return from SGGBAK (computing VSR)
!>                =N+9: error return from SLASCL (various places)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file sgegs.f.

227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBVSL, JOBVSR
234 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
235* ..
236* .. Array Arguments ..
237 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
238 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
239 $ VSR( LDVSR, * ), WORK( * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 REAL ZERO, ONE
246 parameter( zero = 0.0e0, one = 1.0e0 )
247* ..
248* .. Local Scalars ..
249 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
250 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
251 $ ILO, IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
252 $ LWKOPT, NB, NB1, NB2, NB3
253 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254 $ SAFMIN, SMLNUM
255* ..
256* .. External Subroutines ..
257 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slacpy,
258 $ slascl, slaset, sorgqr, sormqr, xerbla
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER ILAENV
263 REAL SLAMCH, SLANGE
264 EXTERNAL ilaenv, lsame, slamch, slange
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC int, max
268* ..
269* .. Executable Statements ..
270*
271* Decode the input arguments
272*
273 IF( lsame( jobvsl, 'N' ) ) THEN
274 ijobvl = 1
275 ilvsl = .false.
276 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
277 ijobvl = 2
278 ilvsl = .true.
279 ELSE
280 ijobvl = -1
281 ilvsl = .false.
282 END IF
283*
284 IF( lsame( jobvsr, 'N' ) ) THEN
285 ijobvr = 1
286 ilvsr = .false.
287 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
288 ijobvr = 2
289 ilvsr = .true.
290 ELSE
291 ijobvr = -1
292 ilvsr = .false.
293 END IF
294*
295* Test the input arguments
296*
297 lwkmin = max( 4*n, 1 )
298 lwkopt = lwkmin
299 work( 1 ) = lwkopt
300 lquery = ( lwork.EQ.-1 )
301 info = 0
302 IF( ijobvl.LE.0 ) THEN
303 info = -1
304 ELSE IF( ijobvr.LE.0 ) THEN
305 info = -2
306 ELSE IF( n.LT.0 ) THEN
307 info = -3
308 ELSE IF( lda.LT.max( 1, n ) ) THEN
309 info = -5
310 ELSE IF( ldb.LT.max( 1, n ) ) THEN
311 info = -7
312 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
313 info = -12
314 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
315 info = -14
316 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
317 info = -16
318 END IF
319*
320 IF( info.EQ.0 ) THEN
321 nb1 = ilaenv( 1, 'SGEQRF', ' ', n, n, -1, -1 )
322 nb2 = ilaenv( 1, 'SORMQR', ' ', n, n, n, -1 )
323 nb3 = ilaenv( 1, 'SORGQR', ' ', n, n, n, -1 )
324 nb = max( nb1, nb2, nb3 )
325 lopt = 2*n+n*(nb+1)
326 work( 1 ) = lopt
327 END IF
328*
329 IF( info.NE.0 ) THEN
330 CALL xerbla( 'SGEGS ', -info )
331 RETURN
332 ELSE IF( lquery ) THEN
333 RETURN
334 END IF
335*
336* Quick return if possible
337*
338 IF( n.EQ.0 )
339 $ RETURN
340*
341* Get machine constants
342*
343 eps = slamch( 'E' )*slamch( 'B' )
344 safmin = slamch( 'S' )
345 smlnum = n*safmin / eps
346 bignum = one / smlnum
347*
348* Scale A if max element outside range [SMLNUM,BIGNUM]
349*
350 anrm = slange( 'M', n, n, a, lda, work )
351 ilascl = .false.
352 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
353 anrmto = smlnum
354 ilascl = .true.
355 ELSE IF( anrm.GT.bignum ) THEN
356 anrmto = bignum
357 ilascl = .true.
358 END IF
359*
360 IF( ilascl ) THEN
361 CALL slascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda, iinfo )
362 IF( iinfo.NE.0 ) THEN
363 info = n + 9
364 RETURN
365 END IF
366 END IF
367*
368* Scale B if max element outside range [SMLNUM,BIGNUM]
369*
370 bnrm = slange( 'M', n, n, b, ldb, work )
371 ilbscl = .false.
372 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
373 bnrmto = smlnum
374 ilbscl = .true.
375 ELSE IF( bnrm.GT.bignum ) THEN
376 bnrmto = bignum
377 ilbscl = .true.
378 END IF
379*
380 IF( ilbscl ) THEN
381 CALL slascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb, iinfo )
382 IF( iinfo.NE.0 ) THEN
383 info = n + 9
384 RETURN
385 END IF
386 END IF
387*
388* Permute the matrix to make it more nearly triangular
389* Workspace layout: (2*N words -- "work..." not actually used)
390* left_permutation, right_permutation, work...
391*
392 ileft = 1
393 iright = n + 1
394 iwork = iright + n
395 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
396 $ work( iright ), work( iwork ), iinfo )
397 IF( iinfo.NE.0 ) THEN
398 info = n + 1
399 GO TO 10
400 END IF
401*
402* Reduce B to triangular form, and initialize VSL and/or VSR
403* Workspace layout: ("work..." must have at least N words)
404* left_permutation, right_permutation, tau, work...
405*
406 irows = ihi + 1 - ilo
407 icols = n + 1 - ilo
408 itau = iwork
409 iwork = itau + irows
410 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
411 $ work( iwork ), lwork+1-iwork, iinfo )
412 IF( iinfo.GE.0 )
413 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
414 IF( iinfo.NE.0 ) THEN
415 info = n + 2
416 GO TO 10
417 END IF
418*
419 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
420 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
421 $ lwork+1-iwork, iinfo )
422 IF( iinfo.GE.0 )
423 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
424 IF( iinfo.NE.0 ) THEN
425 info = n + 3
426 GO TO 10
427 END IF
428*
429 IF( ilvsl ) THEN
430 CALL slaset( 'Full', n, n, zero, one, vsl, ldvsl )
431 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
432 $ vsl( ilo+1, ilo ), ldvsl )
433 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
434 $ work( itau ), work( iwork ), lwork+1-iwork,
435 $ iinfo )
436 IF( iinfo.GE.0 )
437 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
438 IF( iinfo.NE.0 ) THEN
439 info = n + 4
440 GO TO 10
441 END IF
442 END IF
443*
444 IF( ilvsr )
445 $ CALL slaset( 'Full', n, n, zero, one, vsr, ldvsr )
446*
447* Reduce to generalized Hessenberg form
448*
449 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
450 $ ldvsl, vsr, ldvsr, iinfo )
451 IF( iinfo.NE.0 ) THEN
452 info = n + 5
453 GO TO 10
454 END IF
455*
456* Perform QZ algorithm, computing Schur vectors if desired
457* Workspace layout: ("work..." must have at least 1 word)
458* left_permutation, right_permutation, work...
459*
460 iwork = itau
461 CALL shgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
462 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
463 $ work( iwork ), lwork+1-iwork, iinfo )
464 IF( iinfo.GE.0 )
465 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
466 IF( iinfo.NE.0 ) THEN
467 IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
468 info = iinfo
469 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
470 info = iinfo - n
471 ELSE
472 info = n + 6
473 END IF
474 GO TO 10
475 END IF
476*
477* Apply permutation to VSL and VSR
478*
479 IF( ilvsl ) THEN
480 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
481 $ work( iright ), n, vsl, ldvsl, iinfo )
482 IF( iinfo.NE.0 ) THEN
483 info = n + 7
484 GO TO 10
485 END IF
486 END IF
487 IF( ilvsr ) THEN
488 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
489 $ work( iright ), n, vsr, ldvsr, iinfo )
490 IF( iinfo.NE.0 ) THEN
491 info = n + 8
492 GO TO 10
493 END IF
494 END IF
495*
496* Undo scaling
497*
498 IF( ilascl ) THEN
499 CALL slascl( 'H', -1, -1, anrmto, anrm, n, n, a, lda, iinfo )
500 IF( iinfo.NE.0 ) THEN
501 info = n + 9
502 RETURN
503 END IF
504 CALL slascl( 'G', -1, -1, anrmto, anrm, n, 1, alphar, n,
505 $ iinfo )
506 IF( iinfo.NE.0 ) THEN
507 info = n + 9
508 RETURN
509 END IF
510 CALL slascl( 'G', -1, -1, anrmto, anrm, n, 1, alphai, n,
511 $ iinfo )
512 IF( iinfo.NE.0 ) THEN
513 info = n + 9
514 RETURN
515 END IF
516 END IF
517*
518 IF( ilbscl ) THEN
519 CALL slascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb, iinfo )
520 IF( iinfo.NE.0 ) THEN
521 info = n + 9
522 RETURN
523 END IF
524 CALL slascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n, iinfo )
525 IF( iinfo.NE.0 ) THEN
526 info = n + 9
527 RETURN
528 END IF
529 END IF
530*
531 10 CONTINUE
532 work( 1 ) = lwkopt
533*
534 RETURN
535*
536* End of SGEGS
537*
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
sggbak
subroutine sggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
SGGBAK
Definition sggbak.f:147
sggbal
subroutine sggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
SGGBAL
Definition sggbal.f:177
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
shgeqz
subroutine shgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
SHGEQZ
Definition shgeqz.f:304
sgghrd
subroutine sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
Definition sgghrd.f:207
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128

◆ sgegv()

subroutine sgegv ( character jobvl,
character jobvr,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
real, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine SGGEV.
!>
!> SGEGV computes the eigenvalues and, optionally, the left and/or right
!> eigenvectors of a real matrix pair (A,B).
!> Given two square matrices A and B,
!> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
!> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
!> that
!>
!>    A*x = lambda*B*x.
!>
!> An alternate form is to find the eigenvalues mu and corresponding
!> eigenvectors y such that
!>
!>    mu*A*y = B*y.
!>
!> These two forms are equivalent with mu = 1/lambda and x = y if
!> neither lambda nor mu is zero.  In order to deal with the case that
!> lambda or mu is zero or small, two values alpha and beta are returned
!> for each eigenvalue, such that lambda = alpha/beta and
!> mu = beta/alpha.
!>
!> The vectors x and y in the above equations are right eigenvectors of
!> the matrix pair (A,B).  Vectors u and v satisfying
!>
!>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
!>
!> are left eigenvectors of (A,B).
!>
!> Note: this routine performs  on A and B
!>
Parameters
[in]JOBVL
!>          JOBVL is CHARACTER*1
!>          = 'N':  do not compute the left generalized eigenvectors;
!>          = 'V':  compute the left generalized eigenvectors (returned
!>                  in VL).
!>
[in]JOBVR
!>          JOBVR is CHARACTER*1
!>          = 'N':  do not compute the right generalized eigenvectors;
!>          = 'V':  compute the right generalized eigenvectors (returned
!>                  in VR).
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VL, and VR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the matrix A.
!>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
!>          contains the real Schur form of A from the generalized Schur
!>          factorization of the pair (A,B) after balancing.
!>          If no eigenvectors were computed, then only the diagonal
!>          blocks from the Schur form will be correct.  See SGGHRD and
!>          SHGEQZ for details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the matrix B.
!>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
!>          upper triangular matrix obtained from B in the generalized
!>          Schur factorization of the pair (A,B) after balancing.
!>          If no eigenvectors were computed, then only those elements of
!>          B corresponding to the diagonal blocks from the Schur form of
!>          A will be correct.  See SGGHRD and SHGEQZ for details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>          The real parts of each scalar alpha defining an eigenvalue of
!>          GNEP.
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>          The imaginary parts of each scalar alpha defining an
!>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
!>          eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) = -ALPHAI(j).
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>
!>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
!>          beta = BETA(j) represent the j-th eigenvalue of the matrix
!>          pair (A,B), in one of the forms lambda = alpha/beta or
!>          mu = beta/alpha.  Since either lambda or mu may overflow,
!>          they should not, in general, be computed.
!>
[out]VL
!>          VL is REAL array, dimension (LDVL,N)
!>          If JOBVL = 'V', the left eigenvectors u(j) are stored
!>          in the columns of VL, in the same order as their eigenvalues.
!>          If the j-th eigenvalue is real, then u(j) = VL(:,j).
!>          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).
!>
!>          Each eigenvector is scaled so that its largest component has
!>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
!>          corresponding to an eigenvalue with alpha = beta = 0, which
!>          are set to zero.
!>          Not referenced if JOBVL = 'N'.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the matrix VL. LDVL >= 1, and
!>          if JOBVL = 'V', LDVL >= N.
!>
[out]VR
!>          VR is REAL array, dimension (LDVR,N)
!>          If JOBVR = 'V', the right eigenvectors x(j) are stored
!>          in the columns of VR, in the same order as their eigenvalues.
!>          If the j-th eigenvalue is real, then x(j) = VR(:,j).
!>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
!>          pair, then
!>            x(j) = VR(:,j) + i*VR(:,j+1)
!>          and
!>            x(j+1) = VR(:,j) - i*VR(:,j+1).
!>
!>          Each eigenvector is scaled so that its largest component has
!>          abs(real part) + abs(imag. part) = 1, except for eigenvalues
!>          corresponding to an eigenvalue with alpha = beta = 0, which
!>          are set to zero.
!>          Not referenced if JOBVR = 'N'.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          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,8*N).
!>          For good performance, LWORK must generally be larger.
!>          To compute the optimal value of LWORK, call ILAENV to get
!>          blocksizes (for SGEQRF, SORMQR, and SORGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR;
!>          The optimal LWORK is:
!>              2*N + MAX( 6*N, N*(NB+1) ).
!>
!>          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.
!>          = 1,...,N:
!>                The QZ iteration failed.  No eigenvectors have been
!>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  errors that usually indicate LAPACK problems:
!>                =N+1: error return from SGGBAL
!>                =N+2: error return from SGEQRF
!>                =N+3: error return from SORMQR
!>                =N+4: error return from SORGQR
!>                =N+5: error return from SGGHRD
!>                =N+6: error return from SHGEQZ (other than failed
!>                                                iteration)
!>                =N+7: error return from STGEVC
!>                =N+8: error return from SGGBAK (computing VL)
!>                =N+9: error return from SGGBAK (computing VR)
!>                =N+10: error return from SLASCL (various calls)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Balancing
!>  ---------
!>
!>  This driver calls SGGBAL to both permute and scale rows and columns
!>  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
!>  and PL*B*R will be upper triangular except for the diagonal blocks
!>  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
!>  possible.  The diagonal scaling matrices DL and DR are chosen so
!>  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
!>  one (except for the elements that start out zero.)
!>
!>  After the eigenvalues and eigenvectors of the balanced matrices
!>  have been computed, SGGBAK transforms the eigenvectors back to what
!>  they would have been (in perfect arithmetic) if they had not been
!>  balanced.
!>
!>  Contents of A and B on Exit
!>  -------- -- - --- - -- ----
!>
!>  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
!>  both), then on exit the arrays A and B will contain the real Schur
!>  form[*] of the  versions of A and B.  If no eigenvectors
!>  are computed, then only the diagonal blocks will be correct.
!>
!>  [*] See SHGEQZ, SGEGS, or read the book ,
!>      by Golub & van Loan, pub. by Johns Hopkins U. Press.
!>

Definition at line 304 of file sgegv.f.

306*
307* -- LAPACK driver routine --
308* -- LAPACK is a software package provided by Univ. of Tennessee, --
309* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310*
311* .. Scalar Arguments ..
312 CHARACTER JOBVL, JOBVR
313 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314* ..
315* .. Array Arguments ..
316 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
318 $ VR( LDVR, * ), WORK( * )
319* ..
320*
321* =====================================================================
322*
323* .. Parameters ..
324 REAL ZERO, ONE
325 parameter( zero = 0.0e0, one = 1.0e0 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329 CHARACTER CHTEMP
330 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331 $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
332 $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
333 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334 $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
335 $ SALFAI, SALFAR, SBETA, SCALE, TEMP
336* ..
337* .. Local Arrays ..
338 LOGICAL LDUMMA( 1 )
339* ..
340* .. External Subroutines ..
341 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slacpy,
342 $ slascl, slaset, sorgqr, sormqr, stgevc, xerbla
343* ..
344* .. External Functions ..
345 LOGICAL LSAME
346 INTEGER ILAENV
347 REAL SLAMCH, SLANGE
348 EXTERNAL ilaenv, lsame, slamch, slange
349* ..
350* .. Intrinsic Functions ..
351 INTRINSIC abs, int, max
352* ..
353* .. Executable Statements ..
354*
355* Decode the input arguments
356*
357 IF( lsame( jobvl, 'N' ) ) THEN
358 ijobvl = 1
359 ilvl = .false.
360 ELSE IF( lsame( jobvl, 'V' ) ) THEN
361 ijobvl = 2
362 ilvl = .true.
363 ELSE
364 ijobvl = -1
365 ilvl = .false.
366 END IF
367*
368 IF( lsame( jobvr, 'N' ) ) THEN
369 ijobvr = 1
370 ilvr = .false.
371 ELSE IF( lsame( jobvr, 'V' ) ) THEN
372 ijobvr = 2
373 ilvr = .true.
374 ELSE
375 ijobvr = -1
376 ilvr = .false.
377 END IF
378 ilv = ilvl .OR. ilvr
379*
380* Test the input arguments
381*
382 lwkmin = max( 8*n, 1 )
383 lwkopt = lwkmin
384 work( 1 ) = lwkopt
385 lquery = ( lwork.EQ.-1 )
386 info = 0
387 IF( ijobvl.LE.0 ) THEN
388 info = -1
389 ELSE IF( ijobvr.LE.0 ) THEN
390 info = -2
391 ELSE IF( n.LT.0 ) THEN
392 info = -3
393 ELSE IF( lda.LT.max( 1, n ) ) THEN
394 info = -5
395 ELSE IF( ldb.LT.max( 1, n ) ) THEN
396 info = -7
397 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
398 info = -12
399 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
400 info = -14
401 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
402 info = -16
403 END IF
404*
405 IF( info.EQ.0 ) THEN
406 nb1 = ilaenv( 1, 'SGEQRF', ' ', n, n, -1, -1 )
407 nb2 = ilaenv( 1, 'SORMQR', ' ', n, n, n, -1 )
408 nb3 = ilaenv( 1, 'SORGQR', ' ', n, n, n, -1 )
409 nb = max( nb1, nb2, nb3 )
410 lopt = 2*n + max( 6*n, n*(nb+1) )
411 work( 1 ) = lopt
412 END IF
413*
414 IF( info.NE.0 ) THEN
415 CALL xerbla( 'SGEGV ', -info )
416 RETURN
417 ELSE IF( lquery ) THEN
418 RETURN
419 END IF
420*
421* Quick return if possible
422*
423 IF( n.EQ.0 )
424 $ RETURN
425*
426* Get machine constants
427*
428 eps = slamch( 'E' )*slamch( 'B' )
429 safmin = slamch( 'S' )
430 safmin = safmin + safmin
431 safmax = one / safmin
432 onepls = one + ( 4*eps )
433*
434* Scale A
435*
436 anrm = slange( 'M', n, n, a, lda, work )
437 anrm1 = anrm
438 anrm2 = one
439 IF( anrm.LT.one ) THEN
440 IF( safmax*anrm.LT.one ) THEN
441 anrm1 = safmin
442 anrm2 = safmax*anrm
443 END IF
444 END IF
445*
446 IF( anrm.GT.zero ) THEN
447 CALL slascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
448 IF( iinfo.NE.0 ) THEN
449 info = n + 10
450 RETURN
451 END IF
452 END IF
453*
454* Scale B
455*
456 bnrm = slange( 'M', n, n, b, ldb, work )
457 bnrm1 = bnrm
458 bnrm2 = one
459 IF( bnrm.LT.one ) THEN
460 IF( safmax*bnrm.LT.one ) THEN
461 bnrm1 = safmin
462 bnrm2 = safmax*bnrm
463 END IF
464 END IF
465*
466 IF( bnrm.GT.zero ) THEN
467 CALL slascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
468 IF( iinfo.NE.0 ) THEN
469 info = n + 10
470 RETURN
471 END IF
472 END IF
473*
474* Permute the matrix to make it more nearly triangular
475* Workspace layout: (8*N words -- "work" requires 6*N words)
476* left_permutation, right_permutation, work...
477*
478 ileft = 1
479 iright = n + 1
480 iwork = iright + n
481 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
482 $ work( iright ), work( iwork ), iinfo )
483 IF( iinfo.NE.0 ) THEN
484 info = n + 1
485 GO TO 120
486 END IF
487*
488* Reduce B to triangular form, and initialize VL and/or VR
489* Workspace layout: ("work..." must have at least N words)
490* left_permutation, right_permutation, tau, work...
491*
492 irows = ihi + 1 - ilo
493 IF( ilv ) THEN
494 icols = n + 1 - ilo
495 ELSE
496 icols = irows
497 END IF
498 itau = iwork
499 iwork = itau + irows
500 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
501 $ work( iwork ), lwork+1-iwork, iinfo )
502 IF( iinfo.GE.0 )
503 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
504 IF( iinfo.NE.0 ) THEN
505 info = n + 2
506 GO TO 120
507 END IF
508*
509 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
510 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
511 $ lwork+1-iwork, iinfo )
512 IF( iinfo.GE.0 )
513 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514 IF( iinfo.NE.0 ) THEN
515 info = n + 3
516 GO TO 120
517 END IF
518*
519 IF( ilvl ) THEN
520 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
521 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
522 $ vl( ilo+1, ilo ), ldvl )
523 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
524 $ work( itau ), work( iwork ), lwork+1-iwork,
525 $ iinfo )
526 IF( iinfo.GE.0 )
527 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
528 IF( iinfo.NE.0 ) THEN
529 info = n + 4
530 GO TO 120
531 END IF
532 END IF
533*
534 IF( ilvr )
535 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
536*
537* Reduce to generalized Hessenberg form
538*
539 IF( ilv ) THEN
540*
541* Eigenvectors requested -- work on whole matrix.
542*
543 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
544 $ ldvl, vr, ldvr, iinfo )
545 ELSE
546 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
547 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548 END IF
549 IF( iinfo.NE.0 ) THEN
550 info = n + 5
551 GO TO 120
552 END IF
553*
554* Perform QZ algorithm
555* Workspace layout: ("work..." must have at least 1 word)
556* left_permutation, right_permutation, work...
557*
558 iwork = itau
559 IF( ilv ) THEN
560 chtemp = 'S'
561 ELSE
562 chtemp = 'E'
563 END IF
564 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
565 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
566 $ work( iwork ), lwork+1-iwork, iinfo )
567 IF( iinfo.GE.0 )
568 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
569 IF( iinfo.NE.0 ) THEN
570 IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
571 info = iinfo
572 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
573 info = iinfo - n
574 ELSE
575 info = n + 6
576 END IF
577 GO TO 120
578 END IF
579*
580 IF( ilv ) THEN
581*
582* Compute Eigenvectors (STGEVC requires 6*N words of workspace)
583*
584 IF( ilvl ) THEN
585 IF( ilvr ) THEN
586 chtemp = 'B'
587 ELSE
588 chtemp = 'L'
589 END IF
590 ELSE
591 chtemp = 'R'
592 END IF
593*
594 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
595 $ vr, ldvr, n, in, work( iwork ), iinfo )
596 IF( iinfo.NE.0 ) THEN
597 info = n + 7
598 GO TO 120
599 END IF
600*
601* Undo balancing on VL and VR, rescale
602*
603 IF( ilvl ) THEN
604 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
605 $ work( iright ), n, vl, ldvl, iinfo )
606 IF( iinfo.NE.0 ) THEN
607 info = n + 8
608 GO TO 120
609 END IF
610 DO 50 jc = 1, n
611 IF( alphai( jc ).LT.zero )
612 $ GO TO 50
613 temp = zero
614 IF( alphai( jc ).EQ.zero ) THEN
615 DO 10 jr = 1, n
616 temp = max( temp, abs( vl( jr, jc ) ) )
617 10 CONTINUE
618 ELSE
619 DO 20 jr = 1, n
620 temp = max( temp, abs( vl( jr, jc ) )+
621 $ abs( vl( jr, jc+1 ) ) )
622 20 CONTINUE
623 END IF
624 IF( temp.LT.safmin )
625 $ GO TO 50
626 temp = one / temp
627 IF( alphai( jc ).EQ.zero ) THEN
628 DO 30 jr = 1, n
629 vl( jr, jc ) = vl( jr, jc )*temp
630 30 CONTINUE
631 ELSE
632 DO 40 jr = 1, n
633 vl( jr, jc ) = vl( jr, jc )*temp
634 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
635 40 CONTINUE
636 END IF
637 50 CONTINUE
638 END IF
639 IF( ilvr ) THEN
640 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
641 $ work( iright ), n, vr, ldvr, iinfo )
642 IF( iinfo.NE.0 ) THEN
643 info = n + 9
644 GO TO 120
645 END IF
646 DO 100 jc = 1, n
647 IF( alphai( jc ).LT.zero )
648 $ GO TO 100
649 temp = zero
650 IF( alphai( jc ).EQ.zero ) THEN
651 DO 60 jr = 1, n
652 temp = max( temp, abs( vr( jr, jc ) ) )
653 60 CONTINUE
654 ELSE
655 DO 70 jr = 1, n
656 temp = max( temp, abs( vr( jr, jc ) )+
657 $ abs( vr( jr, jc+1 ) ) )
658 70 CONTINUE
659 END IF
660 IF( temp.LT.safmin )
661 $ GO TO 100
662 temp = one / temp
663 IF( alphai( jc ).EQ.zero ) THEN
664 DO 80 jr = 1, n
665 vr( jr, jc ) = vr( jr, jc )*temp
666 80 CONTINUE
667 ELSE
668 DO 90 jr = 1, n
669 vr( jr, jc ) = vr( jr, jc )*temp
670 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
671 90 CONTINUE
672 END IF
673 100 CONTINUE
674 END IF
675*
676* End of eigenvector calculation
677*
678 END IF
679*
680* Undo scaling in alpha, beta
681*
682* Note: this does not give the alpha and beta for the unscaled
683* problem.
684*
685* Un-scaling is limited to avoid underflow in alpha and beta
686* if they are significant.
687*
688 DO 110 jc = 1, n
689 absar = abs( alphar( jc ) )
690 absai = abs( alphai( jc ) )
691 absb = abs( beta( jc ) )
692 salfar = anrm*alphar( jc )
693 salfai = anrm*alphai( jc )
694 sbeta = bnrm*beta( jc )
695 ilimit = .false.
696 scale = one
697*
698* Check for significant underflow in ALPHAI
699*
700 IF( abs( salfai ).LT.safmin .AND. absai.GE.
701 $ max( safmin, eps*absar, eps*absb ) ) THEN
702 ilimit = .true.
703 scale = ( onepls*safmin / anrm1 ) /
704 $ max( onepls*safmin, anrm2*absai )
705*
706 ELSE IF( salfai.EQ.zero ) THEN
707*
708* If insignificant underflow in ALPHAI, then make the
709* conjugate eigenvalue real.
710*
711 IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
712 alphai( jc-1 ) = zero
713 ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
714 alphai( jc+1 ) = zero
715 END IF
716 END IF
717*
718* Check for significant underflow in ALPHAR
719*
720 IF( abs( salfar ).LT.safmin .AND. absar.GE.
721 $ max( safmin, eps*absai, eps*absb ) ) THEN
722 ilimit = .true.
723 scale = max( scale, ( onepls*safmin / anrm1 ) /
724 $ max( onepls*safmin, anrm2*absar ) )
725 END IF
726*
727* Check for significant underflow in BETA
728*
729 IF( abs( sbeta ).LT.safmin .AND. absb.GE.
730 $ max( safmin, eps*absar, eps*absai ) ) THEN
731 ilimit = .true.
732 scale = max( scale, ( onepls*safmin / bnrm1 ) /
733 $ max( onepls*safmin, bnrm2*absb ) )
734 END IF
735*
736* Check for possible overflow when limiting scaling
737*
738 IF( ilimit ) THEN
739 temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
740 $ abs( sbeta ) )
741 IF( temp.GT.one )
742 $ scale = scale / temp
743 IF( scale.LT.one )
744 $ ilimit = .false.
745 END IF
746*
747* Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748*
749 IF( ilimit ) THEN
750 salfar = ( scale*alphar( jc ) )*anrm
751 salfai = ( scale*alphai( jc ) )*anrm
752 sbeta = ( scale*beta( jc ) )*bnrm
753 END IF
754 alphar( jc ) = salfar
755 alphai( jc ) = salfai
756 beta( jc ) = sbeta
757 110 CONTINUE
758*
759 120 CONTINUE
760 work( 1 ) = lwkopt
761*
762 RETURN
763*
764* End of SGEGV
765*
stgevc
subroutine stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC
Definition stgevc.f:295
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339

◆ sgges()

subroutine sgges ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer sdim,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvsl, * ) vsl,
integer ldvsl,
real, dimension( ldvsr, * ) vsr,
integer ldvsr,
real, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

SGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
!>
!> SGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
!> the generalized eigenvalues, the generalized real Schur form (S,T),
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR). This gives the generalized Schur factorization
!>
!>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T.The
!> leading columns of VSL and VSR then form an orthonormal basis for the
!> corresponding left and right eigenspaces (deflating subspaces).
!>
!> (If only the generalized eigenvalues are needed, use the driver
!> SGGEV instead, which is faster.)
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
!> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
!> usually represented as the pair (alpha,beta), as there is a
!> reasonable interpretation for beta=0 or both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
Parameters
[in]JOBVSL
!>          JOBVSL is CHARACTER*1
!>          = 'N':  do not compute the left Schur vectors;
!>          = 'V':  compute the left Schur vectors.
!>
[in]JOBVSR
!>          JOBVSR is CHARACTER*1
!>          = 'N':  do not compute the right Schur vectors;
!>          = 'V':  compute the right Schur vectors.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the generalized Schur form.
!>          = 'N':  Eigenvalues are not ordered;
!>          = 'S':  Eigenvalues are ordered (see SELCTG);
!>
[in]SELCTG
!>          SELCTG is a LOGICAL FUNCTION of three REAL arguments
!>          SELCTG must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'N', SELCTG is not referenced.
!>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>
!>          Note that in the ill-conditioned case, a selected complex
!>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
!>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
!>          in this case.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the first of the pair of matrices.
!>          On exit, A has been overwritten by its generalized Schur
!>          form S.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the second of the pair of matrices.
!>          On exit, B has been overwritten by its generalized Schur
!>          form T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>          for which SELCTG is true.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
!>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI will be always less than and
!>          usually comparable with norm(A) in magnitude, and BETA always
!>          less than and usually comparable with norm(B).
!>
[out]VSL
!>          VSL is REAL array, dimension (LDVSL,N)
!>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
!>          Not referenced if JOBVSL = 'N'.
!>
[in]LDVSL
!>          LDVSL is INTEGER
!>          The leading dimension of the matrix VSL. LDVSL >=1, and
!>          if JOBVSL = 'V', LDVSL >= N.
!>
[out]VSR
!>          VSR is REAL array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
!>          Not referenced if JOBVSR = 'N'.
!>
[in]LDVSR
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= 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.
!>          If N = 0, LWORK >= 1, else LWORK >= max(8*N,6*N+16).
!>          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]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1,...,N:
!>                The QZ iteration failed.  (A,B) are not in Schur
!>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
!>                =N+2: after reordering, roundoff changed values of
!>                      some complex eigenvalues so that leading
!>                      eigenvalues in the Generalized Schur form no
!>                      longer satisfy SELCTG=.TRUE.  This could also
!>                      be caused due to scaling.
!>                =N+3: reordering failed in STGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 281 of file sgges.f.

284*
285* -- LAPACK driver routine --
286* -- LAPACK is a software package provided by Univ. of Tennessee, --
287* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288*
289* .. Scalar Arguments ..
290 CHARACTER JOBVSL, JOBVSR, SORT
291 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
292* ..
293* .. Array Arguments ..
294 LOGICAL BWORK( * )
295 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
296 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
297 $ VSR( LDVSR, * ), WORK( * )
298* ..
299* .. Function Arguments ..
300 LOGICAL SELCTG
301 EXTERNAL selctg
302* ..
303*
304* =====================================================================
305*
306* .. Parameters ..
307 REAL ZERO, ONE
308 parameter( zero = 0.0e+0, one = 1.0e+0 )
309* ..
310* .. Local Scalars ..
311 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
312 $ LQUERY, LST2SL, WANTST
313 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
314 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
315 $ MINWRK
316 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
317 $ PVSR, SAFMAX, SAFMIN, SMLNUM
318* ..
319* .. Local Arrays ..
320 INTEGER IDUM( 1 )
321 REAL DIF( 2 )
322* ..
323* .. External Subroutines ..
324 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
325 $ slacpy, slascl, slaset, sorgqr, sormqr, stgsen,
326 $ xerbla
327* ..
328* .. External Functions ..
329 LOGICAL LSAME
330 INTEGER ILAENV
331 REAL SLAMCH, SLANGE
332 EXTERNAL lsame, ilaenv, slamch, slange
333* ..
334* .. Intrinsic Functions ..
335 INTRINSIC abs, max, sqrt
336* ..
337* .. Executable Statements ..
338*
339* Decode the input arguments
340*
341 IF( lsame( jobvsl, 'N' ) ) THEN
342 ijobvl = 1
343 ilvsl = .false.
344 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
345 ijobvl = 2
346 ilvsl = .true.
347 ELSE
348 ijobvl = -1
349 ilvsl = .false.
350 END IF
351*
352 IF( lsame( jobvsr, 'N' ) ) THEN
353 ijobvr = 1
354 ilvsr = .false.
355 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
356 ijobvr = 2
357 ilvsr = .true.
358 ELSE
359 ijobvr = -1
360 ilvsr = .false.
361 END IF
362*
363 wantst = lsame( sort, 'S' )
364*
365* Test the input arguments
366*
367 info = 0
368 lquery = ( lwork.EQ.-1 )
369 IF( ijobvl.LE.0 ) THEN
370 info = -1
371 ELSE IF( ijobvr.LE.0 ) THEN
372 info = -2
373 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
374 info = -3
375 ELSE IF( n.LT.0 ) THEN
376 info = -5
377 ELSE IF( lda.LT.max( 1, n ) ) THEN
378 info = -7
379 ELSE IF( ldb.LT.max( 1, n ) ) THEN
380 info = -9
381 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
382 info = -15
383 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
384 info = -17
385 END IF
386*
387* Compute workspace
388* (Note: Comments in the code beginning "Workspace:" describe the
389* minimal amount of workspace needed at that point in the code,
390* as well as the preferred amount for good performance.
391* NB refers to the optimal block size for the immediately
392* following subroutine, as returned by ILAENV.)
393*
394 IF( info.EQ.0 ) THEN
395 IF( n.GT.0 )THEN
396 minwrk = max( 8*n, 6*n + 16 )
397 maxwrk = minwrk - n +
398 $ n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 )
399 maxwrk = max( maxwrk, minwrk - n +
400 $ n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, -1 ) )
401 IF( ilvsl ) THEN
402 maxwrk = max( maxwrk, minwrk - n +
403 $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, -1 ) )
404 END IF
405 ELSE
406 minwrk = 1
407 maxwrk = 1
408 END IF
409 work( 1 ) = maxwrk
410*
411 IF( lwork.LT.minwrk .AND. .NOT.lquery )
412 $ info = -19
413 END IF
414*
415 IF( info.NE.0 ) THEN
416 CALL xerbla( 'SGGES ', -info )
417 RETURN
418 ELSE IF( lquery ) THEN
419 RETURN
420 END IF
421*
422* Quick return if possible
423*
424 IF( n.EQ.0 ) THEN
425 sdim = 0
426 RETURN
427 END IF
428*
429* Get machine constants
430*
431 eps = slamch( 'P' )
432 safmin = slamch( 'S' )
433 safmax = one / safmin
434 CALL slabad( safmin, safmax )
435 smlnum = sqrt( safmin ) / eps
436 bignum = one / smlnum
437*
438* Scale A if max element outside range [SMLNUM,BIGNUM]
439*
440 anrm = slange( 'M', n, n, a, lda, work )
441 ilascl = .false.
442 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
443 anrmto = smlnum
444 ilascl = .true.
445 ELSE IF( anrm.GT.bignum ) THEN
446 anrmto = bignum
447 ilascl = .true.
448 END IF
449 IF( ilascl )
450 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
451*
452* Scale B if max element outside range [SMLNUM,BIGNUM]
453*
454 bnrm = slange( 'M', n, n, b, ldb, work )
455 ilbscl = .false.
456 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
457 bnrmto = smlnum
458 ilbscl = .true.
459 ELSE IF( bnrm.GT.bignum ) THEN
460 bnrmto = bignum
461 ilbscl = .true.
462 END IF
463 IF( ilbscl )
464 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
465*
466* Permute the matrix to make it more nearly triangular
467* (Workspace: need 6*N + 2*N space for storing balancing factors)
468*
469 ileft = 1
470 iright = n + 1
471 iwrk = iright + n
472 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
473 $ work( iright ), work( iwrk ), ierr )
474*
475* Reduce B to triangular form (QR decomposition of B)
476* (Workspace: need N, prefer N*NB)
477*
478 irows = ihi + 1 - ilo
479 icols = n + 1 - ilo
480 itau = iwrk
481 iwrk = itau + irows
482 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
483 $ work( iwrk ), lwork+1-iwrk, ierr )
484*
485* Apply the orthogonal transformation to matrix A
486* (Workspace: need N, prefer N*NB)
487*
488 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
489 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
490 $ lwork+1-iwrk, ierr )
491*
492* Initialize VSL
493* (Workspace: need N, prefer N*NB)
494*
495 IF( ilvsl ) THEN
496 CALL slaset( 'Full', n, n, zero, one, vsl, ldvsl )
497 IF( irows.GT.1 ) THEN
498 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
499 $ vsl( ilo+1, ilo ), ldvsl )
500 END IF
501 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
502 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
503 END IF
504*
505* Initialize VSR
506*
507 IF( ilvsr )
508 $ CALL slaset( 'Full', n, n, zero, one, vsr, ldvsr )
509*
510* Reduce to generalized Hessenberg form
511* (Workspace: none needed)
512*
513 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
514 $ ldvsl, vsr, ldvsr, ierr )
515*
516* Perform QZ algorithm, computing Schur vectors if desired
517* (Workspace: need N)
518*
519 iwrk = itau
520 CALL shgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
521 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
522 $ work( iwrk ), lwork+1-iwrk, ierr )
523 IF( ierr.NE.0 ) THEN
524 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
525 info = ierr
526 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
527 info = ierr - n
528 ELSE
529 info = n + 1
530 END IF
531 GO TO 40
532 END IF
533*
534* Sort eigenvalues ALPHA/BETA if desired
535* (Workspace: need 4*N+16 )
536*
537 sdim = 0
538 IF( wantst ) THEN
539*
540* Undo scaling on eigenvalues before SELCTGing
541*
542 IF( ilascl ) THEN
543 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
544 $ ierr )
545 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
546 $ ierr )
547 END IF
548 IF( ilbscl )
549 $ CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
550*
551* Select eigenvalues
552*
553 DO 10 i = 1, n
554 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
555 10 CONTINUE
556*
557 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
558 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
559 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
560 $ ierr )
561 IF( ierr.EQ.1 )
562 $ info = n + 3
563*
564 END IF
565*
566* Apply back-permutation to VSL and VSR
567* (Workspace: none needed)
568*
569 IF( ilvsl )
570 $ CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
571 $ work( iright ), n, vsl, ldvsl, ierr )
572*
573 IF( ilvsr )
574 $ CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
575 $ work( iright ), n, vsr, ldvsr, ierr )
576*
577* Check if unscaling would cause over/underflow, if so, rescale
578* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
579* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
580*
581 IF( ilascl )THEN
582 DO 50 i = 1, n
583 IF( alphai( i ).NE.zero ) THEN
584 IF( ( alphar( i )/safmax ).GT.( anrmto/anrm ) .OR.
585 $ ( safmin/alphar( i ) ).GT.( anrm/anrmto ) ) THEN
586 work( 1 ) = abs( a( i, i )/alphar( i ) )
587 beta( i ) = beta( i )*work( 1 )
588 alphar( i ) = alphar( i )*work( 1 )
589 alphai( i ) = alphai( i )*work( 1 )
590 ELSE IF( ( alphai( i )/safmax ).GT.( anrmto/anrm ) .OR.
591 $ ( safmin/alphai( i ) ).GT.( anrm/anrmto ) ) THEN
592 work( 1 ) = abs( a( i, i+1 )/alphai( i ) )
593 beta( i ) = beta( i )*work( 1 )
594 alphar( i ) = alphar( i )*work( 1 )
595 alphai( i ) = alphai( i )*work( 1 )
596 END IF
597 END IF
598 50 CONTINUE
599 END IF
600*
601 IF( ilbscl )THEN
602 DO 60 i = 1, n
603 IF( alphai( i ).NE.zero ) THEN
604 IF( ( beta( i )/safmax ).GT.( bnrmto/bnrm ) .OR.
605 $ ( safmin/beta( i ) ).GT.( bnrm/bnrmto ) ) THEN
606 work( 1 ) = abs(b( i, i )/beta( i ))
607 beta( i ) = beta( i )*work( 1 )
608 alphar( i ) = alphar( i )*work( 1 )
609 alphai( i ) = alphai( i )*work( 1 )
610 END IF
611 END IF
612 60 CONTINUE
613 END IF
614*
615* Undo scaling
616*
617 IF( ilascl ) THEN
618 CALL slascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
619 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
620 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
621 END IF
622*
623 IF( ilbscl ) THEN
624 CALL slascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
625 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
626 END IF
627*
628 IF( wantst ) THEN
629*
630* Check if reordering is correct
631*
632 lastsl = .true.
633 lst2sl = .true.
634 sdim = 0
635 ip = 0
636 DO 30 i = 1, n
637 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
638 IF( alphai( i ).EQ.zero ) THEN
639 IF( cursl )
640 $ sdim = sdim + 1
641 ip = 0
642 IF( cursl .AND. .NOT.lastsl )
643 $ info = n + 2
644 ELSE
645 IF( ip.EQ.1 ) THEN
646*
647* Last eigenvalue of conjugate pair
648*
649 cursl = cursl .OR. lastsl
650 lastsl = cursl
651 IF( cursl )
652 $ sdim = sdim + 2
653 ip = -1
654 IF( cursl .AND. .NOT.lst2sl )
655 $ info = n + 2
656 ELSE
657*
658* First eigenvalue of conjugate pair
659*
660 ip = 1
661 END IF
662 END IF
663 lst2sl = lastsl
664 lastsl = cursl
665 30 CONTINUE
666*
667 END IF
668*
669 40 CONTINUE
670*
671 work( 1 ) = maxwrk
672*
673 RETURN
674*
675* End of SGGES
676*
stgsen
subroutine stgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, liwork, info)
STGSEN
Definition stgsen.f:451

◆ sgges3()

subroutine sgges3 ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer sdim,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvsl, * ) vsl,
integer ldvsl,
real, dimension( ldvsr, * ) vsr,
integer ldvsr,
real, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

SGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)

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

Purpose:
!>
!> SGGES3 computes for a pair of N-by-N real nonsymmetric matrices (A,B),
!> the generalized eigenvalues, the generalized real Schur form (S,T),
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR). This gives the generalized Schur factorization
!>
!>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T.The
!> leading columns of VSL and VSR then form an orthonormal basis for the
!> corresponding left and right eigenspaces (deflating subspaces).
!>
!> (If only the generalized eigenvalues are needed, use the driver
!> SGGEV instead, which is faster.)
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
!> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
!> usually represented as the pair (alpha,beta), as there is a
!> reasonable interpretation for beta=0 or both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
Parameters
[in]JOBVSL
!>          JOBVSL is CHARACTER*1
!>          = 'N':  do not compute the left Schur vectors;
!>          = 'V':  compute the left Schur vectors.
!>
[in]JOBVSR
!>          JOBVSR is CHARACTER*1
!>          = 'N':  do not compute the right Schur vectors;
!>          = 'V':  compute the right Schur vectors.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the generalized Schur form.
!>          = 'N':  Eigenvalues are not ordered;
!>          = 'S':  Eigenvalues are ordered (see SELCTG);
!>
[in]SELCTG
!>          SELCTG is a LOGICAL FUNCTION of three REAL arguments
!>          SELCTG must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'N', SELCTG is not referenced.
!>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>
!>          Note that in the ill-conditioned case, a selected complex
!>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
!>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
!>          in this case.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the first of the pair of matrices.
!>          On exit, A has been overwritten by its generalized Schur
!>          form S.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the second of the pair of matrices.
!>          On exit, B has been overwritten by its generalized Schur
!>          form T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>          for which SELCTG is true.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
!>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI will be always less than and
!>          usually comparable with norm(A) in magnitude, and BETA always
!>          less than and usually comparable with norm(B).
!>
[out]VSL
!>          VSL is REAL array, dimension (LDVSL,N)
!>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
!>          Not referenced if JOBVSL = 'N'.
!>
[in]LDVSL
!>          LDVSL is INTEGER
!>          The leading dimension of the matrix VSL. LDVSL >=1, and
!>          if JOBVSL = 'V', LDVSL >= N.
!>
[out]VSR
!>          VSR is REAL array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
!>          Not referenced if JOBVSR = 'N'.
!>
[in]LDVSR
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= 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.
!>
!>          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]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1,...,N:
!>                The QZ iteration failed.  (A,B) are not in Schur
!>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SLAQZ0.
!>                =N+2: after reordering, roundoff changed values of
!>                      some complex eigenvalues so that leading
!>                      eigenvalues in the Generalized Schur form no
!>                      longer satisfy SELCTG=.TRUE.  This could also
!>                      be caused due to scaling.
!>                =N+3: reordering failed in STGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 279 of file sgges3.f.

282*
283* -- LAPACK driver routine --
284* -- LAPACK is a software package provided by Univ. of Tennessee, --
285* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
286*
287* .. Scalar Arguments ..
288 CHARACTER JOBVSL, JOBVSR, SORT
289 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
294 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
295 $ VSR( LDVSR, * ), WORK( * )
296* ..
297* .. Function Arguments ..
298 LOGICAL SELCTG
299 EXTERNAL selctg
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 REAL ZERO, ONE
306 parameter( zero = 0.0e+0, one = 1.0e+0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
310 $ LQUERY, LST2SL, WANTST
311 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
312 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, LWKOPT
313 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
314 $ PVSR, SAFMAX, SAFMIN, SMLNUM
315* ..
316* .. Local Arrays ..
317 INTEGER IDUM( 1 )
318 REAL DIF( 2 )
319* ..
320* .. External Subroutines ..
321 EXTERNAL sgeqrf, sggbak, sggbal, sgghd3, slaqz0, slabad,
322 $ slacpy, slascl, slaset, sorgqr, sormqr, stgsen,
323 $ xerbla
324* ..
325* .. External Functions ..
326 LOGICAL LSAME
327 REAL SLAMCH, SLANGE
328 EXTERNAL lsame, slamch, slange
329* ..
330* .. Intrinsic Functions ..
331 INTRINSIC abs, max, sqrt
332* ..
333* .. Executable Statements ..
334*
335* Decode the input arguments
336*
337 IF( lsame( jobvsl, 'N' ) ) THEN
338 ijobvl = 1
339 ilvsl = .false.
340 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
341 ijobvl = 2
342 ilvsl = .true.
343 ELSE
344 ijobvl = -1
345 ilvsl = .false.
346 END IF
347*
348 IF( lsame( jobvsr, 'N' ) ) THEN
349 ijobvr = 1
350 ilvsr = .false.
351 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
352 ijobvr = 2
353 ilvsr = .true.
354 ELSE
355 ijobvr = -1
356 ilvsr = .false.
357 END IF
358*
359 wantst = lsame( sort, 'S' )
360*
361* Test the input arguments
362*
363 info = 0
364 lquery = ( lwork.EQ.-1 )
365 IF( ijobvl.LE.0 ) THEN
366 info = -1
367 ELSE IF( ijobvr.LE.0 ) THEN
368 info = -2
369 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
370 info = -3
371 ELSE IF( n.LT.0 ) THEN
372 info = -5
373 ELSE IF( lda.LT.max( 1, n ) ) THEN
374 info = -7
375 ELSE IF( ldb.LT.max( 1, n ) ) THEN
376 info = -9
377 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
378 info = -15
379 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
380 info = -17
381 ELSE IF( lwork.LT.6*n+16 .AND. .NOT.lquery ) THEN
382 info = -19
383 END IF
384*
385* Compute workspace
386*
387 IF( info.EQ.0 ) THEN
388 CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
389 lwkopt = max( 6*n+16, 3*n+int( work( 1 ) ) )
390 CALL sormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
391 $ -1, ierr )
392 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
393 IF( ilvsl ) THEN
394 CALL sorgqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
395 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
396 END IF
397 CALL sgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
398 $ ldvsl, vsr, ldvsr, work, -1, ierr )
399 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
400 CALL slaqz0( 'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
401 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
402 $ work, -1, 0, ierr )
403 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
404 IF( wantst ) THEN
405 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
406 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
407 $ sdim, pvsl, pvsr, dif, work, -1, idum, 1,
408 $ ierr )
409 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
410 END IF
411 work( 1 ) = lwkopt
412 END IF
413*
414 IF( info.NE.0 ) THEN
415 CALL xerbla( 'SGGES3 ', -info )
416 RETURN
417 ELSE IF( lquery ) THEN
418 RETURN
419 END IF
420*
421* Quick return if possible
422*
423 IF( n.EQ.0 ) THEN
424 sdim = 0
425 RETURN
426 END IF
427*
428* Get machine constants
429*
430 eps = slamch( 'P' )
431 safmin = slamch( 'S' )
432 safmax = one / safmin
433 CALL slabad( safmin, safmax )
434 smlnum = sqrt( safmin ) / eps
435 bignum = one / smlnum
436*
437* Scale A if max element outside range [SMLNUM,BIGNUM]
438*
439 anrm = slange( 'M', n, n, a, lda, work )
440 ilascl = .false.
441 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
442 anrmto = smlnum
443 ilascl = .true.
444 ELSE IF( anrm.GT.bignum ) THEN
445 anrmto = bignum
446 ilascl = .true.
447 END IF
448 IF( ilascl )
449 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
450*
451* Scale B if max element outside range [SMLNUM,BIGNUM]
452*
453 bnrm = slange( 'M', n, n, b, ldb, work )
454 ilbscl = .false.
455 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
456 bnrmto = smlnum
457 ilbscl = .true.
458 ELSE IF( bnrm.GT.bignum ) THEN
459 bnrmto = bignum
460 ilbscl = .true.
461 END IF
462 IF( ilbscl )
463 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
464*
465* Permute the matrix to make it more nearly triangular
466*
467 ileft = 1
468 iright = n + 1
469 iwrk = iright + n
470 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
471 $ work( iright ), work( iwrk ), ierr )
472*
473* Reduce B to triangular form (QR decomposition of B)
474*
475 irows = ihi + 1 - ilo
476 icols = n + 1 - ilo
477 itau = iwrk
478 iwrk = itau + irows
479 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
480 $ work( iwrk ), lwork+1-iwrk, ierr )
481*
482* Apply the orthogonal transformation to matrix A
483*
484 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
485 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
486 $ lwork+1-iwrk, ierr )
487*
488* Initialize VSL
489*
490 IF( ilvsl ) THEN
491 CALL slaset( 'Full', n, n, zero, one, vsl, ldvsl )
492 IF( irows.GT.1 ) THEN
493 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
494 $ vsl( ilo+1, ilo ), ldvsl )
495 END IF
496 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
497 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
498 END IF
499*
500* Initialize VSR
501*
502 IF( ilvsr )
503 $ CALL slaset( 'Full', n, n, zero, one, vsr, ldvsr )
504*
505* Reduce to generalized Hessenberg form
506*
507 CALL sgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
508 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
509*
510* Perform QZ algorithm, computing Schur vectors if desired
511*
512 iwrk = itau
513 CALL slaqz0( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
514 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
515 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
516 IF( ierr.NE.0 ) THEN
517 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
518 info = ierr
519 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
520 info = ierr - n
521 ELSE
522 info = n + 1
523 END IF
524 GO TO 40
525 END IF
526*
527* Sort eigenvalues ALPHA/BETA if desired
528*
529 sdim = 0
530 IF( wantst ) THEN
531*
532* Undo scaling on eigenvalues before SELCTGing
533*
534 IF( ilascl ) THEN
535 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
536 $ ierr )
537 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
538 $ ierr )
539 END IF
540 IF( ilbscl )
541 $ CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
542*
543* Select eigenvalues
544*
545 DO 10 i = 1, n
546 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
547 10 CONTINUE
548*
549 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
550 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
551 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
552 $ ierr )
553 IF( ierr.EQ.1 )
554 $ info = n + 3
555*
556 END IF
557*
558* Apply back-permutation to VSL and VSR
559*
560 IF( ilvsl )
561 $ CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
562 $ work( iright ), n, vsl, ldvsl, ierr )
563*
564 IF( ilvsr )
565 $ CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
566 $ work( iright ), n, vsr, ldvsr, ierr )
567*
568* Check if unscaling would cause over/underflow, if so, rescale
569* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
570* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
571*
572 IF( ilascl )THEN
573 DO 50 i = 1, n
574 IF( alphai( i ).NE.zero ) THEN
575 IF( ( alphar( i )/safmax ).GT.( anrmto/anrm ) .OR.
576 $ ( safmin/alphar( i ) ).GT.( anrm/anrmto ) ) THEN
577 work( 1 ) = abs( a( i, i )/alphar( i ) )
578 beta( i ) = beta( i )*work( 1 )
579 alphar( i ) = alphar( i )*work( 1 )
580 alphai( i ) = alphai( i )*work( 1 )
581 ELSE IF( ( alphai( i )/safmax ).GT.( anrmto/anrm ) .OR.
582 $ ( safmin/alphai( i ) ).GT.( anrm/anrmto ) ) THEN
583 work( 1 ) = abs( a( i, i+1 )/alphai( i ) )
584 beta( i ) = beta( i )*work( 1 )
585 alphar( i ) = alphar( i )*work( 1 )
586 alphai( i ) = alphai( i )*work( 1 )
587 END IF
588 END IF
589 50 CONTINUE
590 END IF
591*
592 IF( ilbscl )THEN
593 DO 60 i = 1, n
594 IF( alphai( i ).NE.zero ) THEN
595 IF( ( beta( i )/safmax ).GT.( bnrmto/bnrm ) .OR.
596 $ ( safmin/beta( i ) ).GT.( bnrm/bnrmto ) ) THEN
597 work( 1 ) = abs(b( i, i )/beta( i ))
598 beta( i ) = beta( i )*work( 1 )
599 alphar( i ) = alphar( i )*work( 1 )
600 alphai( i ) = alphai( i )*work( 1 )
601 END IF
602 END IF
603 60 CONTINUE
604 END IF
605*
606* Undo scaling
607*
608 IF( ilascl ) THEN
609 CALL slascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
610 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
611 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
612 END IF
613*
614 IF( ilbscl ) THEN
615 CALL slascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
616 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
617 END IF
618*
619 IF( wantst ) THEN
620*
621* Check if reordering is correct
622*
623 lastsl = .true.
624 lst2sl = .true.
625 sdim = 0
626 ip = 0
627 DO 30 i = 1, n
628 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
629 IF( alphai( i ).EQ.zero ) THEN
630 IF( cursl )
631 $ sdim = sdim + 1
632 ip = 0
633 IF( cursl .AND. .NOT.lastsl )
634 $ info = n + 2
635 ELSE
636 IF( ip.EQ.1 ) THEN
637*
638* Last eigenvalue of conjugate pair
639*
640 cursl = cursl .OR. lastsl
641 lastsl = cursl
642 IF( cursl )
643 $ sdim = sdim + 2
644 ip = -1
645 IF( cursl .AND. .NOT.lst2sl )
646 $ info = n + 2
647 ELSE
648*
649* First eigenvalue of conjugate pair
650*
651 ip = 1
652 END IF
653 END IF
654 lst2sl = lastsl
655 lastsl = cursl
656 30 CONTINUE
657*
658 END IF
659*
660 40 CONTINUE
661*
662 work( 1 ) = lwkopt
663*
664 RETURN
665*
666* End of SGGES3
667*
slaqz0
recursive subroutine slaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
SLAQZ0
Definition slaqz0.f:304
sgghd3
subroutine sgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
SGGHD3
Definition sgghd3.f:230

◆ sggesx()

subroutine sggesx ( character jobvsl,
character jobvsr,
character sort,
external selctg,
character sense,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer sdim,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvsl, * ) vsl,
integer ldvsl,
real, dimension( ldvsr, * ) vsr,
integer ldvsr,
real, dimension( 2 ) rconde,
real, dimension( 2 ) rcondv,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
logical, dimension( * ) bwork,
integer info )

SGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
!>
!> SGGESX computes for a pair of N-by-N real nonsymmetric matrices
!> (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR).  This gives the generalized Schur factorization
!>
!>      (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T; computes
!> a reciprocal condition number for the average of the selected
!> eigenvalues (RCONDE); and computes a reciprocal condition number for
!> the right and left deflating subspaces corresponding to the selected
!> eigenvalues (RCONDV). The leading columns of VSL and VSR then form
!> an orthonormal basis for the corresponding left and right eigenspaces
!> (deflating subspaces).
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
!> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
!> usually represented as the pair (alpha,beta), as there is a
!> reasonable interpretation for beta=0 or for both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
Parameters
[in]JOBVSL
!>          JOBVSL is CHARACTER*1
!>          = 'N':  do not compute the left Schur vectors;
!>          = 'V':  compute the left Schur vectors.
!>
[in]JOBVSR
!>          JOBVSR is CHARACTER*1
!>          = 'N':  do not compute the right Schur vectors;
!>          = 'V':  compute the right Schur vectors.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the generalized Schur form.
!>          = 'N':  Eigenvalues are not ordered;
!>          = 'S':  Eigenvalues are ordered (see SELCTG).
!>
[in]SELCTG
!>          SELCTG is a LOGICAL FUNCTION of three REAL arguments
!>          SELCTG must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'N', SELCTG is not referenced.
!>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>          Note that a selected complex eigenvalue may no longer satisfy
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
!>          since ordering may change the value of complex eigenvalues
!>          (especially if the eigenvalue is ill-conditioned), in this
!>          case INFO is set to N+3.
!>
[in]SENSE
!>          SENSE is CHARACTER*1
!>          Determines which reciprocal condition numbers are computed.
!>          = 'N':  None are computed;
!>          = 'E':  Computed for average of selected eigenvalues only;
!>          = 'V':  Computed for selected deflating subspaces only;
!>          = 'B':  Computed for both.
!>          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the first of the pair of matrices.
!>          On exit, A has been overwritten by its generalized Schur
!>          form S.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the second of the pair of matrices.
!>          On exit, B has been overwritten by its generalized Schur
!>          form T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>          for which SELCTG is true.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i
!>          and BETA(j),j=1,...,N  are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI will be always less than and
!>          usually comparable with norm(A) in magnitude, and BETA always
!>          less than and usually comparable with norm(B).
!>
[out]VSL
!>          VSL is REAL array, dimension (LDVSL,N)
!>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
!>          Not referenced if JOBVSL = 'N'.
!>
[in]LDVSL
!>          LDVSL is INTEGER
!>          The leading dimension of the matrix VSL. LDVSL >=1, and
!>          if JOBVSL = 'V', LDVSL >= N.
!>
[out]VSR
!>          VSR is REAL array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
!>          Not referenced if JOBVSR = 'N'.
!>
[in]LDVSR
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= N.
!>
[out]RCONDE
!>          RCONDE is REAL array, dimension ( 2 )
!>          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
!>          reciprocal condition numbers for the average of the selected
!>          eigenvalues.
!>          Not referenced if SENSE = 'N' or 'V'.
!>
[out]RCONDV
!>          RCONDV is REAL array, dimension ( 2 )
!>          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
!>          reciprocal condition numbers for the selected deflating
!>          subspaces.
!>          Not referenced if SENSE = 'N' or 'E'.
!>
[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.
!>          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
!>          LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
!>          LWORK >= max( 8*N, 6*N+16 ).
!>          Note that 2*SDIM*(N-SDIM) <= N*N/2.
!>          Note also that an error is only returned if
!>          LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
!>          this may not be large enough.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the bound on the optimal size of the WORK
!>          array and the minimum size of the IWORK array, returns these
!>          values as the first entries of the WORK and IWORK arrays, and
!>          no error message related to LWORK or LIWORK is issued by
!>          XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
!>          LIWORK >= N+6.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the bound on the optimal size of the
!>          WORK array and the minimum size of the IWORK array, returns
!>          these values as the first entries of the WORK and IWORK
!>          arrays, and no error message related to LWORK or LIWORK is
!>          issued by XERBLA.
!>
[out]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1,...,N:
!>                The QZ iteration failed.  (A,B) are not in Schur
!>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SHGEQZ
!>                =N+2: after reordering, roundoff changed values of
!>                      some complex eigenvalues so that leading
!>                      eigenvalues in the Generalized Schur form no
!>                      longer satisfy SELCTG=.TRUE.  This could also
!>                      be caused due to scaling.
!>                =N+3: reordering failed in STGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  An approximate (asymptotic) bound on the average absolute error of
!>  the selected eigenvalues is
!>
!>       EPS * norm((A, B)) / RCONDE( 1 ).
!>
!>  An approximate (asymptotic) bound on the maximum angular error in
!>  the computed deflating subspaces is
!>
!>       EPS * norm((A, B)) / RCONDV( 2 ).
!>
!>  See LAPACK User's Guide, section 4.11 for more information.
!>

Definition at line 361 of file sggesx.f.

365*
366* -- LAPACK driver routine --
367* -- LAPACK is a software package provided by Univ. of Tennessee, --
368* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
369*
370* .. Scalar Arguments ..
371 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
372 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
373 $ SDIM
374* ..
375* .. Array Arguments ..
376 LOGICAL BWORK( * )
377 INTEGER IWORK( * )
378 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
379 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
380 $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
381 $ WORK( * )
382* ..
383* .. Function Arguments ..
384 LOGICAL SELCTG
385 EXTERNAL selctg
386* ..
387*
388* =====================================================================
389*
390* .. Parameters ..
391 REAL ZERO, ONE
392 parameter( zero = 0.0e+0, one = 1.0e+0 )
393* ..
394* .. Local Scalars ..
395 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
396 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
397 $ WANTSV
398 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
399 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
400 $ LIWMIN, LWRK, MAXWRK, MINWRK
401 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
402 $ PR, SAFMAX, SAFMIN, SMLNUM
403* ..
404* .. Local Arrays ..
405 REAL DIF( 2 )
406* ..
407* .. External Subroutines ..
408 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
409 $ slacpy, slascl, slaset, sorgqr, sormqr, stgsen,
410 $ xerbla
411* ..
412* .. External Functions ..
413 LOGICAL LSAME
414 INTEGER ILAENV
415 REAL SLAMCH, SLANGE
416 EXTERNAL lsame, ilaenv, slamch, slange
417* ..
418* .. Intrinsic Functions ..
419 INTRINSIC abs, max, sqrt
420* ..
421* .. Executable Statements ..
422*
423* Decode the input arguments
424*
425 IF( lsame( jobvsl, 'N' ) ) THEN
426 ijobvl = 1
427 ilvsl = .false.
428 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
429 ijobvl = 2
430 ilvsl = .true.
431 ELSE
432 ijobvl = -1
433 ilvsl = .false.
434 END IF
435*
436 IF( lsame( jobvsr, 'N' ) ) THEN
437 ijobvr = 1
438 ilvsr = .false.
439 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
440 ijobvr = 2
441 ilvsr = .true.
442 ELSE
443 ijobvr = -1
444 ilvsr = .false.
445 END IF
446*
447 wantst = lsame( sort, 'S' )
448 wantsn = lsame( sense, 'N' )
449 wantse = lsame( sense, 'E' )
450 wantsv = lsame( sense, 'V' )
451 wantsb = lsame( sense, 'B' )
452 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
453 IF( wantsn ) THEN
454 ijob = 0
455 ELSE IF( wantse ) THEN
456 ijob = 1
457 ELSE IF( wantsv ) THEN
458 ijob = 2
459 ELSE IF( wantsb ) THEN
460 ijob = 4
461 END IF
462*
463* Test the input arguments
464*
465 info = 0
466 IF( ijobvl.LE.0 ) THEN
467 info = -1
468 ELSE IF( ijobvr.LE.0 ) THEN
469 info = -2
470 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
471 info = -3
472 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
473 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
474 info = -5
475 ELSE IF( n.LT.0 ) THEN
476 info = -6
477 ELSE IF( lda.LT.max( 1, n ) ) THEN
478 info = -8
479 ELSE IF( ldb.LT.max( 1, n ) ) THEN
480 info = -10
481 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
482 info = -16
483 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
484 info = -18
485 END IF
486*
487* Compute workspace
488* (Note: Comments in the code beginning "Workspace:" describe the
489* minimal amount of workspace needed at that point in the code,
490* as well as the preferred amount for good performance.
491* NB refers to the optimal block size for the immediately
492* following subroutine, as returned by ILAENV.)
493*
494 IF( info.EQ.0 ) THEN
495 IF( n.GT.0) THEN
496 minwrk = max( 8*n, 6*n + 16 )
497 maxwrk = minwrk - n +
498 $ n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 )
499 maxwrk = max( maxwrk, minwrk - n +
500 $ n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, -1 ) )
501 IF( ilvsl ) THEN
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, -1 ) )
504 END IF
505 lwrk = maxwrk
506 IF( ijob.GE.1 )
507 $ lwrk = max( lwrk, n*n/2 )
508 ELSE
509 minwrk = 1
510 maxwrk = 1
511 lwrk = 1
512 END IF
513 work( 1 ) = lwrk
514 IF( wantsn .OR. n.EQ.0 ) THEN
515 liwmin = 1
516 ELSE
517 liwmin = n + 6
518 END IF
519 iwork( 1 ) = liwmin
520*
521 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
522 info = -22
523 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
524 info = -24
525 END IF
526 END IF
527*
528 IF( info.NE.0 ) THEN
529 CALL xerbla( 'SGGESX', -info )
530 RETURN
531 ELSE IF (lquery) THEN
532 RETURN
533 END IF
534*
535* Quick return if possible
536*
537 IF( n.EQ.0 ) THEN
538 sdim = 0
539 RETURN
540 END IF
541*
542* Get machine constants
543*
544 eps = slamch( 'P' )
545 safmin = slamch( 'S' )
546 safmax = one / safmin
547 CALL slabad( safmin, safmax )
548 smlnum = sqrt( safmin ) / eps
549 bignum = one / smlnum
550*
551* Scale A if max element outside range [SMLNUM,BIGNUM]
552*
553 anrm = slange( 'M', n, n, a, lda, work )
554 ilascl = .false.
555 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
556 anrmto = smlnum
557 ilascl = .true.
558 ELSE IF( anrm.GT.bignum ) THEN
559 anrmto = bignum
560 ilascl = .true.
561 END IF
562 IF( ilascl )
563 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
564*
565* Scale B if max element outside range [SMLNUM,BIGNUM]
566*
567 bnrm = slange( 'M', n, n, b, ldb, work )
568 ilbscl = .false.
569 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
570 bnrmto = smlnum
571 ilbscl = .true.
572 ELSE IF( bnrm.GT.bignum ) THEN
573 bnrmto = bignum
574 ilbscl = .true.
575 END IF
576 IF( ilbscl )
577 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
578*
579* Permute the matrix to make it more nearly triangular
580* (Workspace: need 6*N + 2*N for permutation parameters)
581*
582 ileft = 1
583 iright = n + 1
584 iwrk = iright + n
585 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
586 $ work( iright ), work( iwrk ), ierr )
587*
588* Reduce B to triangular form (QR decomposition of B)
589* (Workspace: need N, prefer N*NB)
590*
591 irows = ihi + 1 - ilo
592 icols = n + 1 - ilo
593 itau = iwrk
594 iwrk = itau + irows
595 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
596 $ work( iwrk ), lwork+1-iwrk, ierr )
597*
598* Apply the orthogonal transformation to matrix A
599* (Workspace: need N, prefer N*NB)
600*
601 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
602 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
603 $ lwork+1-iwrk, ierr )
604*
605* Initialize VSL
606* (Workspace: need N, prefer N*NB)
607*
608 IF( ilvsl ) THEN
609 CALL slaset( 'Full', n, n, zero, one, vsl, ldvsl )
610 IF( irows.GT.1 ) THEN
611 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
612 $ vsl( ilo+1, ilo ), ldvsl )
613 END IF
614 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
615 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
616 END IF
617*
618* Initialize VSR
619*
620 IF( ilvsr )
621 $ CALL slaset( 'Full', n, n, zero, one, vsr, ldvsr )
622*
623* Reduce to generalized Hessenberg form
624* (Workspace: none needed)
625*
626 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
627 $ ldvsl, vsr, ldvsr, ierr )
628*
629 sdim = 0
630*
631* Perform QZ algorithm, computing Schur vectors if desired
632* (Workspace: need N)
633*
634 iwrk = itau
635 CALL shgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
636 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
637 $ work( iwrk ), lwork+1-iwrk, ierr )
638 IF( ierr.NE.0 ) THEN
639 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
640 info = ierr
641 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
642 info = ierr - n
643 ELSE
644 info = n + 1
645 END IF
646 GO TO 50
647 END IF
648*
649* Sort eigenvalues ALPHA/BETA and compute the reciprocal of
650* condition number(s)
651* (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
652* otherwise, need 8*(N+1) )
653*
654 IF( wantst ) THEN
655*
656* Undo scaling on eigenvalues before SELCTGing
657*
658 IF( ilascl ) THEN
659 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
660 $ ierr )
661 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
662 $ ierr )
663 END IF
664 IF( ilbscl )
665 $ CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
666*
667* Select eigenvalues
668*
669 DO 10 i = 1, n
670 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
671 10 CONTINUE
672*
673* Reorder eigenvalues, transform Generalized Schur vectors, and
674* compute reciprocal condition numbers
675*
676 CALL stgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
677 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
678 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
679 $ iwork, liwork, ierr )
680*
681 IF( ijob.GE.1 )
682 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
683 IF( ierr.EQ.-22 ) THEN
684*
685* not enough real workspace
686*
687 info = -22
688 ELSE
689 IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
690 rconde( 1 ) = pl
691 rconde( 2 ) = pr
692 END IF
693 IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
694 rcondv( 1 ) = dif( 1 )
695 rcondv( 2 ) = dif( 2 )
696 END IF
697 IF( ierr.EQ.1 )
698 $ info = n + 3
699 END IF
700*
701 END IF
702*
703* Apply permutation to VSL and VSR
704* (Workspace: none needed)
705*
706 IF( ilvsl )
707 $ CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
708 $ work( iright ), n, vsl, ldvsl, ierr )
709*
710 IF( ilvsr )
711 $ CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
712 $ work( iright ), n, vsr, ldvsr, ierr )
713*
714* Check if unscaling would cause over/underflow, if so, rescale
715* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
716* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
717*
718 IF( ilascl ) THEN
719 DO 20 i = 1, n
720 IF( alphai( i ).NE.zero ) THEN
721 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
722 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
723 $ THEN
724 work( 1 ) = abs( a( i, i ) / alphar( i ) )
725 beta( i ) = beta( i )*work( 1 )
726 alphar( i ) = alphar( i )*work( 1 )
727 alphai( i ) = alphai( i )*work( 1 )
728 ELSE IF( ( alphai( i ) / safmax ).GT.( anrmto / anrm )
729 $ .OR. ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
730 $ THEN
731 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
732 beta( i ) = beta( i )*work( 1 )
733 alphar( i ) = alphar( i )*work( 1 )
734 alphai( i ) = alphai( i )*work( 1 )
735 END IF
736 END IF
737 20 CONTINUE
738 END IF
739*
740 IF( ilbscl ) THEN
741 DO 25 i = 1, n
742 IF( alphai( i ).NE.zero ) THEN
743 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
744 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
745 work( 1 ) = abs( b( i, i ) / beta( i ) )
746 beta( i ) = beta( i )*work( 1 )
747 alphar( i ) = alphar( i )*work( 1 )
748 alphai( i ) = alphai( i )*work( 1 )
749 END IF
750 END IF
751 25 CONTINUE
752 END IF
753*
754* Undo scaling
755*
756 IF( ilascl ) THEN
757 CALL slascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
758 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
759 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
760 END IF
761*
762 IF( ilbscl ) THEN
763 CALL slascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
764 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
765 END IF
766*
767 IF( wantst ) THEN
768*
769* Check if reordering is correct
770*
771 lastsl = .true.
772 lst2sl = .true.
773 sdim = 0
774 ip = 0
775 DO 40 i = 1, n
776 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
777 IF( alphai( i ).EQ.zero ) THEN
778 IF( cursl )
779 $ sdim = sdim + 1
780 ip = 0
781 IF( cursl .AND. .NOT.lastsl )
782 $ info = n + 2
783 ELSE
784 IF( ip.EQ.1 ) THEN
785*
786* Last eigenvalue of conjugate pair
787*
788 cursl = cursl .OR. lastsl
789 lastsl = cursl
790 IF( cursl )
791 $ sdim = sdim + 2
792 ip = -1
793 IF( cursl .AND. .NOT.lst2sl )
794 $ info = n + 2
795 ELSE
796*
797* First eigenvalue of conjugate pair
798*
799 ip = 1
800 END IF
801 END IF
802 lst2sl = lastsl
803 lastsl = cursl
804 40 CONTINUE
805*
806 END IF
807*
808 50 CONTINUE
809*
810 work( 1 ) = maxwrk
811 iwork( 1 ) = liwmin
812*
813 RETURN
814*
815* End of SGGESX
816*

◆ sggev()

subroutine sggev ( character jobvl,
character jobvr,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
real, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> SGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)
!> the generalized eigenvalues, and optionally, the left and/or right
!> generalized eigenvectors.
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
!> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
!> singular. It is usually represented as the pair (alpha,beta), as
!> there is a reasonable interpretation for beta=0, and even for both
!> being zero.
!>
!> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  A * v(j) = lambda(j) * B * v(j).
!>
!> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
!>
!> where u(j)**H is the conjugate-transpose of u(j).
!>
!>
Parameters
[in]JOBVL
!>          JOBVL is CHARACTER*1
!>          = 'N':  do not compute the left generalized eigenvectors;
!>          = 'V':  compute the left generalized eigenvectors.
!>
[in]JOBVR
!>          JOBVR is CHARACTER*1
!>          = 'N':  do not compute the right generalized eigenvectors;
!>          = 'V':  compute the right generalized eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VL, and VR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the matrix A in the pair (A,B).
!>          On exit, A has been overwritten.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the matrix B in the pair (A,B).
!>          On exit, B has been overwritten.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio
!>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
!>          than and usually comparable with norm(A) in magnitude, and
!>          BETA always less than and usually comparable with norm(B).
!>
[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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part)+abs(imag. part)=1.
!>          Not referenced if JOBVL = 'N'.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the matrix VL. LDVL >= 1, and
!>          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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part)+abs(imag. part)=1.
!>          Not referenced if JOBVR = 'N'.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          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,8*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.
!>          = 1,...,N:
!>                The QZ iteration failed.  No eigenvectors have been
!>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
!>                =N+2: error return from STGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file sggev.f.

226*
227* -- LAPACK driver routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 CHARACTER JOBVL, JOBVR
233 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234* ..
235* .. Array Arguments ..
236 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
238 $ VR( LDVR, * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 REAL ZERO, ONE
245 parameter( zero = 0.0e+0, one = 1.0e+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249 CHARACTER CHTEMP
250 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
252 $ MINWRK
253 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254 $ SMLNUM, TEMP
255* ..
256* .. Local Arrays ..
257 LOGICAL LDUMMA( 1 )
258* ..
259* .. External Subroutines ..
260 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
261 $ slacpy, slascl, slaset, sorgqr, sormqr, stgevc,
262 $ xerbla
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 INTEGER ILAENV
267 REAL SLAMCH, SLANGE
268 EXTERNAL lsame, ilaenv, slamch, slange
269* ..
270* .. Intrinsic Functions ..
271 INTRINSIC abs, max, sqrt
272* ..
273* .. Executable Statements ..
274*
275* Decode the input arguments
276*
277 IF( lsame( jobvl, 'N' ) ) THEN
278 ijobvl = 1
279 ilvl = .false.
280 ELSE IF( lsame( jobvl, 'V' ) ) THEN
281 ijobvl = 2
282 ilvl = .true.
283 ELSE
284 ijobvl = -1
285 ilvl = .false.
286 END IF
287*
288 IF( lsame( jobvr, 'N' ) ) THEN
289 ijobvr = 1
290 ilvr = .false.
291 ELSE IF( lsame( jobvr, 'V' ) ) THEN
292 ijobvr = 2
293 ilvr = .true.
294 ELSE
295 ijobvr = -1
296 ilvr = .false.
297 END IF
298 ilv = ilvl .OR. ilvr
299*
300* Test the input arguments
301*
302 info = 0
303 lquery = ( lwork.EQ.-1 )
304 IF( ijobvl.LE.0 ) THEN
305 info = -1
306 ELSE IF( ijobvr.LE.0 ) THEN
307 info = -2
308 ELSE IF( n.LT.0 ) THEN
309 info = -3
310 ELSE IF( lda.LT.max( 1, n ) ) THEN
311 info = -5
312 ELSE IF( ldb.LT.max( 1, n ) ) THEN
313 info = -7
314 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
315 info = -12
316 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317 info = -14
318 END IF
319*
320* Compute workspace
321* (Note: Comments in the code beginning "Workspace:" describe the
322* minimal amount of workspace needed at that point in the code,
323* as well as the preferred amount for good performance.
324* NB refers to the optimal block size for the immediately
325* following subroutine, as returned by ILAENV. The workspace is
326* computed assuming ILO = 1 and IHI = N, the worst case.)
327*
328 IF( info.EQ.0 ) THEN
329 minwrk = max( 1, 8*n )
330 maxwrk = max( 1, n*( 7 +
331 $ ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) ) )
332 maxwrk = max( maxwrk, n*( 7 +
333 $ ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) ) )
334 IF( ilvl ) THEN
335 maxwrk = max( maxwrk, n*( 7 +
336 $ ilaenv( 1, 'SORGQR', ' ', n, 1, n, -1 ) ) )
337 END IF
338 work( 1 ) = maxwrk
339*
340 IF( lwork.LT.minwrk .AND. .NOT.lquery )
341 $ info = -16
342 END IF
343*
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'SGGEV ', -info )
346 RETURN
347 ELSE IF( lquery ) THEN
348 RETURN
349 END IF
350*
351* Quick return if possible
352*
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Get machine constants
357*
358 eps = slamch( 'P' )
359 smlnum = slamch( 'S' )
360 bignum = one / smlnum
361 CALL slabad( smlnum, bignum )
362 smlnum = sqrt( smlnum ) / eps
363 bignum = one / smlnum
364*
365* Scale A if max element outside range [SMLNUM,BIGNUM]
366*
367 anrm = slange( 'M', n, n, a, lda, work )
368 ilascl = .false.
369 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
370 anrmto = smlnum
371 ilascl = .true.
372 ELSE IF( anrm.GT.bignum ) THEN
373 anrmto = bignum
374 ilascl = .true.
375 END IF
376 IF( ilascl )
377 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
378*
379* Scale B if max element outside range [SMLNUM,BIGNUM]
380*
381 bnrm = slange( 'M', n, n, b, ldb, work )
382 ilbscl = .false.
383 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
384 bnrmto = smlnum
385 ilbscl = .true.
386 ELSE IF( bnrm.GT.bignum ) THEN
387 bnrmto = bignum
388 ilbscl = .true.
389 END IF
390 IF( ilbscl )
391 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
392*
393* Permute the matrices A, B to isolate eigenvalues if possible
394* (Workspace: need 6*N)
395*
396 ileft = 1
397 iright = n + 1
398 iwrk = iright + n
399 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
400 $ work( iright ), work( iwrk ), ierr )
401*
402* Reduce B to triangular form (QR decomposition of B)
403* (Workspace: need N, prefer N*NB)
404*
405 irows = ihi + 1 - ilo
406 IF( ilv ) THEN
407 icols = n + 1 - ilo
408 ELSE
409 icols = irows
410 END IF
411 itau = iwrk
412 iwrk = itau + irows
413 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
414 $ work( iwrk ), lwork+1-iwrk, ierr )
415*
416* Apply the orthogonal transformation to matrix A
417* (Workspace: need N, prefer N*NB)
418*
419 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
420 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
421 $ lwork+1-iwrk, ierr )
422*
423* Initialize VL
424* (Workspace: need N, prefer N*NB)
425*
426 IF( ilvl ) THEN
427 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
428 IF( irows.GT.1 ) THEN
429 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
430 $ vl( ilo+1, ilo ), ldvl )
431 END IF
432 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
433 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
434 END IF
435*
436* Initialize VR
437*
438 IF( ilvr )
439 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
440*
441* Reduce to generalized Hessenberg form
442* (Workspace: none needed)
443*
444 IF( ilv ) THEN
445*
446* Eigenvectors requested -- work on whole matrix.
447*
448 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449 $ ldvl, vr, ldvr, ierr )
450 ELSE
451 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
452 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
453 END IF
454*
455* Perform QZ algorithm (Compute eigenvalues, and optionally, the
456* Schur forms and Schur vectors)
457* (Workspace: need N)
458*
459 iwrk = itau
460 IF( ilv ) THEN
461 chtemp = 'S'
462 ELSE
463 chtemp = 'E'
464 END IF
465 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
466 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
467 $ work( iwrk ), lwork+1-iwrk, ierr )
468 IF( ierr.NE.0 ) THEN
469 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
470 info = ierr
471 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
472 info = ierr - n
473 ELSE
474 info = n + 1
475 END IF
476 GO TO 110
477 END IF
478*
479* Compute Eigenvectors
480* (Workspace: need 6*N)
481*
482 IF( ilv ) THEN
483 IF( ilvl ) THEN
484 IF( ilvr ) THEN
485 chtemp = 'B'
486 ELSE
487 chtemp = 'L'
488 END IF
489 ELSE
490 chtemp = 'R'
491 END IF
492 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
493 $ vr, ldvr, n, in, work( iwrk ), ierr )
494 IF( ierr.NE.0 ) THEN
495 info = n + 2
496 GO TO 110
497 END IF
498*
499* Undo balancing on VL and VR and normalization
500* (Workspace: none needed)
501*
502 IF( ilvl ) THEN
503 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
504 $ work( iright ), n, vl, ldvl, ierr )
505 DO 50 jc = 1, n
506 IF( alphai( jc ).LT.zero )
507 $ GO TO 50
508 temp = zero
509 IF( alphai( jc ).EQ.zero ) THEN
510 DO 10 jr = 1, n
511 temp = max( temp, abs( vl( jr, jc ) ) )
512 10 CONTINUE
513 ELSE
514 DO 20 jr = 1, n
515 temp = max( temp, abs( vl( jr, jc ) )+
516 $ abs( vl( jr, jc+1 ) ) )
517 20 CONTINUE
518 END IF
519 IF( temp.LT.smlnum )
520 $ GO TO 50
521 temp = one / temp
522 IF( alphai( jc ).EQ.zero ) THEN
523 DO 30 jr = 1, n
524 vl( jr, jc ) = vl( jr, jc )*temp
525 30 CONTINUE
526 ELSE
527 DO 40 jr = 1, n
528 vl( jr, jc ) = vl( jr, jc )*temp
529 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
530 40 CONTINUE
531 END IF
532 50 CONTINUE
533 END IF
534 IF( ilvr ) THEN
535 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
536 $ work( iright ), n, vr, ldvr, ierr )
537 DO 100 jc = 1, n
538 IF( alphai( jc ).LT.zero )
539 $ GO TO 100
540 temp = zero
541 IF( alphai( jc ).EQ.zero ) THEN
542 DO 60 jr = 1, n
543 temp = max( temp, abs( vr( jr, jc ) ) )
544 60 CONTINUE
545 ELSE
546 DO 70 jr = 1, n
547 temp = max( temp, abs( vr( jr, jc ) )+
548 $ abs( vr( jr, jc+1 ) ) )
549 70 CONTINUE
550 END IF
551 IF( temp.LT.smlnum )
552 $ GO TO 100
553 temp = one / temp
554 IF( alphai( jc ).EQ.zero ) THEN
555 DO 80 jr = 1, n
556 vr( jr, jc ) = vr( jr, jc )*temp
557 80 CONTINUE
558 ELSE
559 DO 90 jr = 1, n
560 vr( jr, jc ) = vr( jr, jc )*temp
561 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
562 90 CONTINUE
563 END IF
564 100 CONTINUE
565 END IF
566*
567* End of eigenvector calculation
568*
569 END IF
570*
571* Undo scaling if necessary
572*
573 110 CONTINUE
574*
575 IF( ilascl ) THEN
576 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
577 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
578 END IF
579*
580 IF( ilbscl ) THEN
581 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
582 END IF
583*
584 work( 1 ) = maxwrk
585 RETURN
586*
587* End of SGGEV
588*

◆ sggev3()

subroutine sggev3 ( character jobvl,
character jobvr,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
real, dimension( * ) work,
integer lwork,
integer info )

SGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)

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

Purpose:
!>
!> SGGEV3 computes for a pair of N-by-N real nonsymmetric matrices (A,B)
!> the generalized eigenvalues, and optionally, the left and/or right
!> generalized eigenvectors.
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
!> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
!> singular. It is usually represented as the pair (alpha,beta), as
!> there is a reasonable interpretation for beta=0, and even for both
!> being zero.
!>
!> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  A * v(j) = lambda(j) * B * v(j).
!>
!> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
!>
!> where u(j)**H is the conjugate-transpose of u(j).
!>
!>
Parameters
[in]JOBVL
!>          JOBVL is CHARACTER*1
!>          = 'N':  do not compute the left generalized eigenvectors;
!>          = 'V':  compute the left generalized eigenvectors.
!>
[in]JOBVR
!>          JOBVR is CHARACTER*1
!>          = 'N':  do not compute the right generalized eigenvectors;
!>          = 'V':  compute the right generalized eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VL, and VR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the matrix A in the pair (A,B).
!>          On exit, A has been overwritten.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the matrix B in the pair (A,B).
!>          On exit, B has been overwritten.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio
!>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
!>          than and usually comparable with norm(A) in magnitude, and
!>          BETA always less than and usually comparable with norm(B).
!>
[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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part)+abs(imag. part)=1.
!>          Not referenced if JOBVL = 'N'.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the matrix VL. LDVL >= 1, and
!>          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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part)+abs(imag. part)=1.
!>          Not referenced if JOBVR = 'N'.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          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
!>
!>          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.
!>          = 1,...,N:
!>                The QZ iteration failed.  No eigenvectors have been
!>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SLAQZ0.
!>                =N+2: error return from STGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file sggev3.f.

226*
227* -- LAPACK driver routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 CHARACTER JOBVL, JOBVR
233 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234* ..
235* .. Array Arguments ..
236 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
238 $ VR( LDVR, * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 REAL ZERO, ONE
245 parameter( zero = 0.0e+0, one = 1.0e+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249 CHARACTER CHTEMP
250 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, LWKOPT
252 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
253 $ SMLNUM, TEMP
254* ..
255* .. Local Arrays ..
256 LOGICAL LDUMMA( 1 )
257* ..
258* .. External Subroutines ..
259 EXTERNAL sgeqrf, sggbak, sggbal, sgghd3, slaqz0, slabad,
260 $ slacpy, slascl, slaset, sorgqr, sormqr, stgevc,
261 $ xerbla
262* ..
263* .. External Functions ..
264 LOGICAL LSAME
265 REAL SLAMCH, SLANGE
266 EXTERNAL lsame, slamch, slange
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC abs, max, sqrt
270* ..
271* .. Executable Statements ..
272*
273* Decode the input arguments
274*
275 IF( lsame( jobvl, 'N' ) ) THEN
276 ijobvl = 1
277 ilvl = .false.
278 ELSE IF( lsame( jobvl, 'V' ) ) THEN
279 ijobvl = 2
280 ilvl = .true.
281 ELSE
282 ijobvl = -1
283 ilvl = .false.
284 END IF
285*
286 IF( lsame( jobvr, 'N' ) ) THEN
287 ijobvr = 1
288 ilvr = .false.
289 ELSE IF( lsame( jobvr, 'V' ) ) THEN
290 ijobvr = 2
291 ilvr = .true.
292 ELSE
293 ijobvr = -1
294 ilvr = .false.
295 END IF
296 ilv = ilvl .OR. ilvr
297*
298* Test the input arguments
299*
300 info = 0
301 lquery = ( lwork.EQ.-1 )
302 IF( ijobvl.LE.0 ) THEN
303 info = -1
304 ELSE IF( ijobvr.LE.0 ) THEN
305 info = -2
306 ELSE IF( n.LT.0 ) THEN
307 info = -3
308 ELSE IF( lda.LT.max( 1, n ) ) THEN
309 info = -5
310 ELSE IF( ldb.LT.max( 1, n ) ) THEN
311 info = -7
312 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
313 info = -12
314 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
315 info = -14
316 ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
317 info = -16
318 END IF
319*
320* Compute workspace
321*
322 IF( info.EQ.0 ) THEN
323 CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
324 lwkopt = max( 1, 8*n, 3*n+int( work( 1 ) ) )
325 CALL sormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
326 $ -1, ierr )
327 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
328 CALL sgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl, ldvl,
329 $ vr, ldvr, work, -1, ierr )
330 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
331 IF( ilvl ) THEN
332 CALL sorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
333 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
334 CALL slaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
335 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
336 $ work, -1, 0, ierr )
337 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
338 ELSE
339 CALL slaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
340 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
341 $ work, -1, 0, ierr )
342 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
343 END IF
344 work( 1 ) = real( lwkopt )
345*
346 END IF
347*
348 IF( info.NE.0 ) THEN
349 CALL xerbla( 'SGGEV3 ', -info )
350 RETURN
351 ELSE IF( lquery ) THEN
352 RETURN
353 END IF
354*
355* Quick return if possible
356*
357 IF( n.EQ.0 )
358 $ RETURN
359*
360* Get machine constants
361*
362 eps = slamch( 'P' )
363 smlnum = slamch( 'S' )
364 bignum = one / smlnum
365 CALL slabad( smlnum, bignum )
366 smlnum = sqrt( smlnum ) / eps
367 bignum = one / smlnum
368*
369* Scale A if max element outside range [SMLNUM,BIGNUM]
370*
371 anrm = slange( 'M', n, n, a, lda, work )
372 ilascl = .false.
373 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
374 anrmto = smlnum
375 ilascl = .true.
376 ELSE IF( anrm.GT.bignum ) THEN
377 anrmto = bignum
378 ilascl = .true.
379 END IF
380 IF( ilascl )
381 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
382*
383* Scale B if max element outside range [SMLNUM,BIGNUM]
384*
385 bnrm = slange( 'M', n, n, b, ldb, work )
386 ilbscl = .false.
387 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
388 bnrmto = smlnum
389 ilbscl = .true.
390 ELSE IF( bnrm.GT.bignum ) THEN
391 bnrmto = bignum
392 ilbscl = .true.
393 END IF
394 IF( ilbscl )
395 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
396*
397* Permute the matrices A, B to isolate eigenvalues if possible
398*
399 ileft = 1
400 iright = n + 1
401 iwrk = iright + n
402 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
403 $ work( iright ), work( iwrk ), ierr )
404*
405* Reduce B to triangular form (QR decomposition of B)
406*
407 irows = ihi + 1 - ilo
408 IF( ilv ) THEN
409 icols = n + 1 - ilo
410 ELSE
411 icols = irows
412 END IF
413 itau = iwrk
414 iwrk = itau + irows
415 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
416 $ work( iwrk ), lwork+1-iwrk, ierr )
417*
418* Apply the orthogonal transformation to matrix A
419*
420 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
421 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
422 $ lwork+1-iwrk, ierr )
423*
424* Initialize VL
425*
426 IF( ilvl ) THEN
427 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
428 IF( irows.GT.1 ) THEN
429 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
430 $ vl( ilo+1, ilo ), ldvl )
431 END IF
432 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
433 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
434 END IF
435*
436* Initialize VR
437*
438 IF( ilvr )
439 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
440*
441* Reduce to generalized Hessenberg form
442*
443 IF( ilv ) THEN
444*
445* Eigenvectors requested -- work on whole matrix.
446*
447 CALL sgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
448 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
449 ELSE
450 CALL sgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
451 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
452 $ work( iwrk ), lwork+1-iwrk, ierr )
453 END IF
454*
455* Perform QZ algorithm (Compute eigenvalues, and optionally, the
456* Schur forms and Schur vectors)
457*
458 iwrk = itau
459 IF( ilv ) THEN
460 chtemp = 'S'
461 ELSE
462 chtemp = 'E'
463 END IF
464 CALL slaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
465 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
466 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
467 IF( ierr.NE.0 ) THEN
468 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
469 info = ierr
470 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
471 info = ierr - n
472 ELSE
473 info = n + 1
474 END IF
475 GO TO 110
476 END IF
477*
478* Compute Eigenvectors
479*
480 IF( ilv ) THEN
481 IF( ilvl ) THEN
482 IF( ilvr ) THEN
483 chtemp = 'B'
484 ELSE
485 chtemp = 'L'
486 END IF
487 ELSE
488 chtemp = 'R'
489 END IF
490 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
491 $ vr, ldvr, n, in, work( iwrk ), ierr )
492 IF( ierr.NE.0 ) THEN
493 info = n + 2
494 GO TO 110
495 END IF
496*
497* Undo balancing on VL and VR and normalization
498*
499 IF( ilvl ) THEN
500 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
501 $ work( iright ), n, vl, ldvl, ierr )
502 DO 50 jc = 1, n
503 IF( alphai( jc ).LT.zero )
504 $ GO TO 50
505 temp = zero
506 IF( alphai( jc ).EQ.zero ) THEN
507 DO 10 jr = 1, n
508 temp = max( temp, abs( vl( jr, jc ) ) )
509 10 CONTINUE
510 ELSE
511 DO 20 jr = 1, n
512 temp = max( temp, abs( vl( jr, jc ) )+
513 $ abs( vl( jr, jc+1 ) ) )
514 20 CONTINUE
515 END IF
516 IF( temp.LT.smlnum )
517 $ GO TO 50
518 temp = one / temp
519 IF( alphai( jc ).EQ.zero ) THEN
520 DO 30 jr = 1, n
521 vl( jr, jc ) = vl( jr, jc )*temp
522 30 CONTINUE
523 ELSE
524 DO 40 jr = 1, n
525 vl( jr, jc ) = vl( jr, jc )*temp
526 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
527 40 CONTINUE
528 END IF
529 50 CONTINUE
530 END IF
531 IF( ilvr ) THEN
532 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
533 $ work( iright ), n, vr, ldvr, ierr )
534 DO 100 jc = 1, n
535 IF( alphai( jc ).LT.zero )
536 $ GO TO 100
537 temp = zero
538 IF( alphai( jc ).EQ.zero ) THEN
539 DO 60 jr = 1, n
540 temp = max( temp, abs( vr( jr, jc ) ) )
541 60 CONTINUE
542 ELSE
543 DO 70 jr = 1, n
544 temp = max( temp, abs( vr( jr, jc ) )+
545 $ abs( vr( jr, jc+1 ) ) )
546 70 CONTINUE
547 END IF
548 IF( temp.LT.smlnum )
549 $ GO TO 100
550 temp = one / temp
551 IF( alphai( jc ).EQ.zero ) THEN
552 DO 80 jr = 1, n
553 vr( jr, jc ) = vr( jr, jc )*temp
554 80 CONTINUE
555 ELSE
556 DO 90 jr = 1, n
557 vr( jr, jc ) = vr( jr, jc )*temp
558 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
559 90 CONTINUE
560 END IF
561 100 CONTINUE
562 END IF
563*
564* End of eigenvector calculation
565*
566 END IF
567*
568* Undo scaling if necessary
569*
570 110 CONTINUE
571*
572 IF( ilascl ) THEN
573 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
574 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
575 END IF
576*
577 IF( ilbscl ) THEN
578 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
579 END IF
580*
581 work( 1 ) = real( lwkopt )
582 RETURN
583*
584* End of SGGEV3
585*

◆ sggevx()

subroutine sggevx ( character balanc,
character jobvl,
character jobvr,
character sense,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
integer ilo,
integer ihi,
real, dimension( * ) lscale,
real, dimension( * ) rscale,
real abnrm,
real bbnrm,
real, dimension( * ) rconde,
real, dimension( * ) rcondv,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> SGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
!> the generalized eigenvalues, and optionally, the left and/or right
!> generalized eigenvectors.
!>
!> Optionally also, it computes a balancing transformation to improve
!> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
!> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
!> the eigenvalues (RCONDE), and reciprocal condition numbers for the
!> right eigenvectors (RCONDV).
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
!> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
!> singular. It is usually represented as the pair (alpha,beta), as
!> there is a reasonable interpretation for beta=0, and even for both
!> being zero.
!>
!> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  A * v(j) = lambda(j) * B * v(j) .
!>
!> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  u(j)**H * A  = lambda(j) * u(j)**H * B.
!>
!> where u(j)**H is the conjugate-transpose of u(j).
!>
!>
Parameters
[in]BALANC
!>          BALANC is CHARACTER*1
!>          Specifies the balance option to be performed.
!>          = 'N':  do not diagonally scale or permute;
!>          = 'P':  permute only;
!>          = 'S':  scale only;
!>          = 'B':  both permute and scale.
!>          Computed reciprocal condition numbers will be for the
!>          matrices after permuting and/or balancing. Permuting does
!>          not change condition numbers (in exact arithmetic), but
!>          balancing does.
!>
[in]JOBVL
!>          JOBVL is CHARACTER*1
!>          = 'N':  do not compute the left generalized eigenvectors;
!>          = 'V':  compute the left generalized eigenvectors.
!>
[in]JOBVR
!>          JOBVR is CHARACTER*1
!>          = 'N':  do not compute the right generalized eigenvectors;
!>          = 'V':  compute the right generalized eigenvectors.
!>
[in]SENSE
!>          SENSE is CHARACTER*1
!>          Determines which reciprocal condition numbers are computed.
!>          = 'N': none are computed;
!>          = 'E': computed for eigenvalues only;
!>          = 'V': computed for eigenvectors only;
!>          = 'B': computed for eigenvalues and eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VL, and VR.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the matrix A in the pair (A,B).
!>          On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
!>          or both, then A contains the first part of the real Schur
!>          form of the  versions of the input A and B.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the matrix B in the pair (A,B).
!>          On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
!>          or both, then B contains the second part of the real Schur
!>          form of the  versions of the input A and B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio
!>          ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
!>          than and usually comparable with norm(A) in magnitude, and
!>          BETA always less than and usually comparable with norm(B).
!>
[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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector will be scaled so the largest component have
!>          abs(real part) + abs(imag. part) = 1.
!>          Not referenced if JOBVL = 'N'.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the matrix VL. LDVL >= 1, and
!>          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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th 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).
!>          Each eigenvector will be scaled so the largest component have
!>          abs(real part) + abs(imag. part) = 1.
!>          Not referenced if JOBVR = 'N'.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          if JOBVR = 'V', LDVR >= N.
!>
[out]ILO
!>          ILO is INTEGER
!>
[out]IHI
!>          IHI is INTEGER
!>          ILO and IHI are integer values such that on exit
!>          A(i,j) = 0 and B(i,j) = 0 if i > j and
!>          j = 1,...,ILO-1 or i = IHI+1,...,N.
!>          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
!>
[out]LSCALE
!>          LSCALE is REAL array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the left side of A and B.  If PL(j) is the index of the
!>          row interchanged with row j, and DL(j) is the scaling
!>          factor applied to row j, then
!>            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
!>                      = DL(j)  for j = ILO,...,IHI
!>                      = PL(j)  for j = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]RSCALE
!>          RSCALE is REAL array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the right side of A and B.  If PR(j) is the index of the
!>          column interchanged with column j, and DR(j) is the scaling
!>          factor applied to column j, then
!>            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
!>                      = DR(j)  for j = ILO,...,IHI
!>                      = PR(j)  for j = IHI+1,...,N
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]ABNRM
!>          ABNRM is REAL
!>          The one-norm of the balanced matrix A.
!>
[out]BBNRM
!>          BBNRM is REAL
!>          The one-norm of the balanced matrix B.
!>
[out]RCONDE
!>          RCONDE is REAL array, dimension (N)
!>          If SENSE = 'E' or 'B', the reciprocal condition numbers of
!>          the eigenvalues, stored in consecutive elements of the array.
!>          For a complex conjugate pair of eigenvalues two consecutive
!>          elements of RCONDE are set to the same value. Thus RCONDE(j),
!>          RCONDV(j), and the j-th columns of VL and VR all correspond
!>          to the j-th eigenpair.
!>          If SENSE = 'N' or 'V', RCONDE is not referenced.
!>
[out]RCONDV
!>          RCONDV is REAL array, dimension (N)
!>          If SENSE = 'V' or 'B', the estimated reciprocal condition
!>          numbers of the eigenvectors, stored in consecutive elements
!>          of the array. For a complex eigenvector two consecutive
!>          elements of RCONDV are set to the same value. If the
!>          eigenvalues cannot be reordered to compute RCONDV(j),
!>          RCONDV(j) is set to 0; this can only occur when the true
!>          value would be very small anyway.
!>          If SENSE = 'N' or 'E', RCONDV is not referenced.
!>
[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,2*N).
!>          If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
!>          LWORK >= max(1,6*N).
!>          If SENSE = 'E', LWORK >= max(1,10*N).
!>          If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (N+6)
!>          If SENSE = 'E', IWORK is not referenced.
!>
[out]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          If SENSE = 'N', BWORK is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1,...,N:
!>                The QZ iteration failed.  No eigenvectors have been
!>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
!>                =N+2: error return from STGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Balancing a matrix pair (A,B) includes, first, permuting rows and
!>  columns to isolate eigenvalues, second, applying diagonal similarity
!>  transformation to the rows and columns to make the rows and columns
!>  as close in norm as possible. The computed reciprocal condition
!>  numbers correspond to the balanced matrix. Permuting rows and columns
!>  will not change the condition numbers (in exact arithmetic) but
!>  diagonal scaling will.  For further explanation of balancing, see
!>  section 4.11.1.2 of LAPACK Users' Guide.
!>
!>  An approximate error bound on the chordal distance between the i-th
!>  computed generalized eigenvalue w and the corresponding exact
!>  eigenvalue lambda is
!>
!>       chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
!>
!>  An approximate error bound for the angle between the i-th computed
!>  eigenvector VL(i) or VR(i) is given by
!>
!>       EPS * norm(ABNRM, BBNRM) / DIF(i).
!>
!>  For further explanation of the reciprocal condition numbers RCONDE
!>  and RCONDV, see section 4.11 of LAPACK User's Guide.
!>

Definition at line 387 of file sggevx.f.

391*
392* -- LAPACK driver routine --
393* -- LAPACK is a software package provided by Univ. of Tennessee, --
394* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
395*
396* .. Scalar Arguments ..
397 CHARACTER BALANC, JOBVL, JOBVR, SENSE
398 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
399 REAL ABNRM, BBNRM
400* ..
401* .. Array Arguments ..
402 LOGICAL BWORK( * )
403 INTEGER IWORK( * )
404 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
405 $ B( LDB, * ), BETA( * ), LSCALE( * ),
406 $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
407 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
408* ..
409*
410* =====================================================================
411*
412* .. Parameters ..
413 REAL ZERO, ONE
414 parameter( zero = 0.0e+0, one = 1.0e+0 )
415* ..
416* .. Local Scalars ..
417 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
418 $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
419 CHARACTER CHTEMP
420 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
421 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
422 $ MINWRK, MM
423 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
424 $ SMLNUM, TEMP
425* ..
426* .. Local Arrays ..
427 LOGICAL LDUMMA( 1 )
428* ..
429* .. External Subroutines ..
430 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
431 $ slacpy, slascl, slaset, sorgqr, sormqr, stgevc,
432 $ stgsna, xerbla
433* ..
434* .. External Functions ..
435 LOGICAL LSAME
436 INTEGER ILAENV
437 REAL SLAMCH, SLANGE
438 EXTERNAL lsame, ilaenv, slamch, slange
439* ..
440* .. Intrinsic Functions ..
441 INTRINSIC abs, max, sqrt
442* ..
443* .. Executable Statements ..
444*
445* Decode the input arguments
446*
447 IF( lsame( jobvl, 'N' ) ) THEN
448 ijobvl = 1
449 ilvl = .false.
450 ELSE IF( lsame( jobvl, 'V' ) ) THEN
451 ijobvl = 2
452 ilvl = .true.
453 ELSE
454 ijobvl = -1
455 ilvl = .false.
456 END IF
457*
458 IF( lsame( jobvr, 'N' ) ) THEN
459 ijobvr = 1
460 ilvr = .false.
461 ELSE IF( lsame( jobvr, 'V' ) ) THEN
462 ijobvr = 2
463 ilvr = .true.
464 ELSE
465 ijobvr = -1
466 ilvr = .false.
467 END IF
468 ilv = ilvl .OR. ilvr
469*
470 noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
471 wantsn = lsame( sense, 'N' )
472 wantse = lsame( sense, 'E' )
473 wantsv = lsame( sense, 'V' )
474 wantsb = lsame( sense, 'B' )
475*
476* Test the input arguments
477*
478 info = 0
479 lquery = ( lwork.EQ.-1 )
480 IF( .NOT.( noscl .OR. lsame( balanc, 'S' ) .OR.
481 $ lsame( balanc, 'B' ) ) ) THEN
482 info = -1
483 ELSE IF( ijobvl.LE.0 ) THEN
484 info = -2
485 ELSE IF( ijobvr.LE.0 ) THEN
486 info = -3
487 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
488 $ THEN
489 info = -4
490 ELSE IF( n.LT.0 ) THEN
491 info = -5
492 ELSE IF( lda.LT.max( 1, n ) ) THEN
493 info = -7
494 ELSE IF( ldb.LT.max( 1, n ) ) THEN
495 info = -9
496 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
497 info = -14
498 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
499 info = -16
500 END IF
501*
502* Compute workspace
503* (Note: Comments in the code beginning "Workspace:" describe the
504* minimal amount of workspace needed at that point in the code,
505* as well as the preferred amount for good performance.
506* NB refers to the optimal block size for the immediately
507* following subroutine, as returned by ILAENV. The workspace is
508* computed assuming ILO = 1 and IHI = N, the worst case.)
509*
510 IF( info.EQ.0 ) THEN
511 IF( n.EQ.0 ) THEN
512 minwrk = 1
513 maxwrk = 1
514 ELSE
515 IF( noscl .AND. .NOT.ilv ) THEN
516 minwrk = 2*n
517 ELSE
518 minwrk = 6*n
519 END IF
520 IF( wantse ) THEN
521 minwrk = 10*n
522 ELSE IF( wantsv .OR. wantsb ) THEN
523 minwrk = 2*n*( n + 4 ) + 16
524 END IF
525 maxwrk = minwrk
526 maxwrk = max( maxwrk,
527 $ n + n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) )
528 maxwrk = max( maxwrk,
529 $ n + n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) )
530 IF( ilvl ) THEN
531 maxwrk = max( maxwrk, n +
532 $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, 0 ) )
533 END IF
534 END IF
535 work( 1 ) = maxwrk
536*
537 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
538 info = -26
539 END IF
540 END IF
541*
542 IF( info.NE.0 ) THEN
543 CALL xerbla( 'SGGEVX', -info )
544 RETURN
545 ELSE IF( lquery ) THEN
546 RETURN
547 END IF
548*
549* Quick return if possible
550*
551 IF( n.EQ.0 )
552 $ RETURN
553*
554*
555* Get machine constants
556*
557 eps = slamch( 'P' )
558 smlnum = slamch( 'S' )
559 bignum = one / smlnum
560 CALL slabad( smlnum, bignum )
561 smlnum = sqrt( smlnum ) / eps
562 bignum = one / smlnum
563*
564* Scale A if max element outside range [SMLNUM,BIGNUM]
565*
566 anrm = slange( 'M', n, n, a, lda, work )
567 ilascl = .false.
568 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
569 anrmto = smlnum
570 ilascl = .true.
571 ELSE IF( anrm.GT.bignum ) THEN
572 anrmto = bignum
573 ilascl = .true.
574 END IF
575 IF( ilascl )
576 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
577*
578* Scale B if max element outside range [SMLNUM,BIGNUM]
579*
580 bnrm = slange( 'M', n, n, b, ldb, work )
581 ilbscl = .false.
582 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
583 bnrmto = smlnum
584 ilbscl = .true.
585 ELSE IF( bnrm.GT.bignum ) THEN
586 bnrmto = bignum
587 ilbscl = .true.
588 END IF
589 IF( ilbscl )
590 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
591*
592* Permute and/or balance the matrix pair (A,B)
593* (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
594*
595 CALL sggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
596 $ work, ierr )
597*
598* Compute ABNRM and BBNRM
599*
600 abnrm = slange( '1', n, n, a, lda, work( 1 ) )
601 IF( ilascl ) THEN
602 work( 1 ) = abnrm
603 CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
604 $ ierr )
605 abnrm = work( 1 )
606 END IF
607*
608 bbnrm = slange( '1', n, n, b, ldb, work( 1 ) )
609 IF( ilbscl ) THEN
610 work( 1 ) = bbnrm
611 CALL slascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
612 $ ierr )
613 bbnrm = work( 1 )
614 END IF
615*
616* Reduce B to triangular form (QR decomposition of B)
617* (Workspace: need N, prefer N*NB )
618*
619 irows = ihi + 1 - ilo
620 IF( ilv .OR. .NOT.wantsn ) THEN
621 icols = n + 1 - ilo
622 ELSE
623 icols = irows
624 END IF
625 itau = 1
626 iwrk = itau + irows
627 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
628 $ work( iwrk ), lwork+1-iwrk, ierr )
629*
630* Apply the orthogonal transformation to A
631* (Workspace: need N, prefer N*NB)
632*
633 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
634 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
635 $ lwork+1-iwrk, ierr )
636*
637* Initialize VL and/or VR
638* (Workspace: need N, prefer N*NB)
639*
640 IF( ilvl ) THEN
641 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
642 IF( irows.GT.1 ) THEN
643 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
644 $ vl( ilo+1, ilo ), ldvl )
645 END IF
646 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
647 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
648 END IF
649*
650 IF( ilvr )
651 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
652*
653* Reduce to generalized Hessenberg form
654* (Workspace: none needed)
655*
656 IF( ilv .OR. .NOT.wantsn ) THEN
657*
658* Eigenvectors requested -- work on whole matrix.
659*
660 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
661 $ ldvl, vr, ldvr, ierr )
662 ELSE
663 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
664 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
665 END IF
666*
667* Perform QZ algorithm (Compute eigenvalues, and optionally, the
668* Schur forms and Schur vectors)
669* (Workspace: need N)
670*
671 IF( ilv .OR. .NOT.wantsn ) THEN
672 chtemp = 'S'
673 ELSE
674 chtemp = 'E'
675 END IF
676*
677 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
678 $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
679 $ lwork, ierr )
680 IF( ierr.NE.0 ) THEN
681 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
682 info = ierr
683 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
684 info = ierr - n
685 ELSE
686 info = n + 1
687 END IF
688 GO TO 130
689 END IF
690*
691* Compute Eigenvectors and estimate condition numbers if desired
692* (Workspace: STGEVC: need 6*N
693* STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
694* need N otherwise )
695*
696 IF( ilv .OR. .NOT.wantsn ) THEN
697 IF( ilv ) THEN
698 IF( ilvl ) THEN
699 IF( ilvr ) THEN
700 chtemp = 'B'
701 ELSE
702 chtemp = 'L'
703 END IF
704 ELSE
705 chtemp = 'R'
706 END IF
707*
708 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
709 $ ldvl, vr, ldvr, n, in, work, ierr )
710 IF( ierr.NE.0 ) THEN
711 info = n + 2
712 GO TO 130
713 END IF
714 END IF
715*
716 IF( .NOT.wantsn ) THEN
717*
718* compute eigenvectors (STGEVC) and estimate condition
719* numbers (STGSNA). Note that the definition of the condition
720* number is not invariant under transformation (u,v) to
721* (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
722* Schur form (S,T), Q and Z are orthogonal matrices. In order
723* to avoid using extra 2*N*N workspace, we have to recalculate
724* eigenvectors and estimate one condition numbers at a time.
725*
726 pair = .false.
727 DO 20 i = 1, n
728*
729 IF( pair ) THEN
730 pair = .false.
731 GO TO 20
732 END IF
733 mm = 1
734 IF( i.LT.n ) THEN
735 IF( a( i+1, i ).NE.zero ) THEN
736 pair = .true.
737 mm = 2
738 END IF
739 END IF
740*
741 DO 10 j = 1, n
742 bwork( j ) = .false.
743 10 CONTINUE
744 IF( mm.EQ.1 ) THEN
745 bwork( i ) = .true.
746 ELSE IF( mm.EQ.2 ) THEN
747 bwork( i ) = .true.
748 bwork( i+1 ) = .true.
749 END IF
750*
751 iwrk = mm*n + 1
752 iwrk1 = iwrk + mm*n
753*
754* Compute a pair of left and right eigenvectors.
755* (compute workspace: need up to 4*N + 6*N)
756*
757 IF( wantse .OR. wantsb ) THEN
758 CALL stgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
759 $ work( 1 ), n, work( iwrk ), n, mm, m,
760 $ work( iwrk1 ), ierr )
761 IF( ierr.NE.0 ) THEN
762 info = n + 2
763 GO TO 130
764 END IF
765 END IF
766*
767 CALL stgsna( sense, 'S', bwork, n, a, lda, b, ldb,
768 $ work( 1 ), n, work( iwrk ), n, rconde( i ),
769 $ rcondv( i ), mm, m, work( iwrk1 ),
770 $ lwork-iwrk1+1, iwork, ierr )
771*
772 20 CONTINUE
773 END IF
774 END IF
775*
776* Undo balancing on VL and VR and normalization
777* (Workspace: none needed)
778*
779 IF( ilvl ) THEN
780 CALL sggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
781 $ ldvl, ierr )
782*
783 DO 70 jc = 1, n
784 IF( alphai( jc ).LT.zero )
785 $ GO TO 70
786 temp = zero
787 IF( alphai( jc ).EQ.zero ) THEN
788 DO 30 jr = 1, n
789 temp = max( temp, abs( vl( jr, jc ) ) )
790 30 CONTINUE
791 ELSE
792 DO 40 jr = 1, n
793 temp = max( temp, abs( vl( jr, jc ) )+
794 $ abs( vl( jr, jc+1 ) ) )
795 40 CONTINUE
796 END IF
797 IF( temp.LT.smlnum )
798 $ GO TO 70
799 temp = one / temp
800 IF( alphai( jc ).EQ.zero ) THEN
801 DO 50 jr = 1, n
802 vl( jr, jc ) = vl( jr, jc )*temp
803 50 CONTINUE
804 ELSE
805 DO 60 jr = 1, n
806 vl( jr, jc ) = vl( jr, jc )*temp
807 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
808 60 CONTINUE
809 END IF
810 70 CONTINUE
811 END IF
812 IF( ilvr ) THEN
813 CALL sggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
814 $ ldvr, ierr )
815 DO 120 jc = 1, n
816 IF( alphai( jc ).LT.zero )
817 $ GO TO 120
818 temp = zero
819 IF( alphai( jc ).EQ.zero ) THEN
820 DO 80 jr = 1, n
821 temp = max( temp, abs( vr( jr, jc ) ) )
822 80 CONTINUE
823 ELSE
824 DO 90 jr = 1, n
825 temp = max( temp, abs( vr( jr, jc ) )+
826 $ abs( vr( jr, jc+1 ) ) )
827 90 CONTINUE
828 END IF
829 IF( temp.LT.smlnum )
830 $ GO TO 120
831 temp = one / temp
832 IF( alphai( jc ).EQ.zero ) THEN
833 DO 100 jr = 1, n
834 vr( jr, jc ) = vr( jr, jc )*temp
835 100 CONTINUE
836 ELSE
837 DO 110 jr = 1, n
838 vr( jr, jc ) = vr( jr, jc )*temp
839 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
840 110 CONTINUE
841 END IF
842 120 CONTINUE
843 END IF
844*
845* Undo scaling if necessary
846*
847 130 CONTINUE
848*
849 IF( ilascl ) THEN
850 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
851 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
852 END IF
853*
854 IF( ilbscl ) THEN
855 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
856 END IF
857*
858 work( 1 ) = maxwrk
859 RETURN
860*
861* End of SGGEVX
862*
stgsna
subroutine stgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
STGSNA
Definition stgsna.f:381