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

◆ 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, real Schur form, and, optionally, the left and/or right Schur vectors of a real matrix pair (A,B)

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 222 of file sgegs.f.

225*
226* -- LAPACK driver routine --
227* -- LAPACK is a software package provided by Univ. of Tennessee, --
228* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
229*
230* .. Scalar Arguments ..
231 CHARACTER JOBVSL, JOBVSR
232 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
233* ..
234* .. Array Arguments ..
235 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
236 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
237 $ VSR( LDVSR, * ), WORK( * )
238* ..
239*
240* =====================================================================
241*
242* .. Parameters ..
243 REAL ZERO, ONE
244 parameter( zero = 0.0e0, one = 1.0e0 )
245* ..
246* .. Local Scalars ..
247 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
248 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
249 $ ILO, IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
250 $ LWKOPT, NB, NB1, NB2, NB3
251 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
252 $ SAFMIN, SMLNUM
253* ..
254* .. External Subroutines ..
255 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slacpy,
256 $ slascl, slaset, sorgqr, sormqr, xerbla
257* ..
258* .. External Functions ..
259 LOGICAL LSAME
260 INTEGER ILAENV
261 REAL SLAMCH, SLANGE
262 EXTERNAL ilaenv, lsame, slamch, slange
263* ..
264* .. Intrinsic Functions ..
265 INTRINSIC int, max
266* ..
267* .. Executable Statements ..
268*
269* Decode the input arguments
270*
271 IF( lsame( jobvsl, 'N' ) ) THEN
272 ijobvl = 1
273 ilvsl = .false.
274 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
275 ijobvl = 2
276 ilvsl = .true.
277 ELSE
278 ijobvl = -1
279 ilvsl = .false.
280 END IF
281*
282 IF( lsame( jobvsr, 'N' ) ) THEN
283 ijobvr = 1
284 ilvsr = .false.
285 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
286 ijobvr = 2
287 ilvsr = .true.
288 ELSE
289 ijobvr = -1
290 ilvsr = .false.
291 END IF
292*
293* Test the input arguments
294*
295 lwkmin = max( 4*n, 1 )
296 lwkopt = lwkmin
297 work( 1 ) = lwkopt
298 lquery = ( lwork.EQ.-1 )
299 info = 0
300 IF( ijobvl.LE.0 ) THEN
301 info = -1
302 ELSE IF( ijobvr.LE.0 ) THEN
303 info = -2
304 ELSE IF( n.LT.0 ) THEN
305 info = -3
306 ELSE IF( lda.LT.max( 1, n ) ) THEN
307 info = -5
308 ELSE IF( ldb.LT.max( 1, n ) ) THEN
309 info = -7
310 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
311 info = -12
312 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
313 info = -14
314 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
315 info = -16
316 END IF
317*
318 IF( info.EQ.0 ) THEN
319 nb1 = ilaenv( 1, 'SGEQRF', ' ', n, n, -1, -1 )
320 nb2 = ilaenv( 1, 'SORMQR', ' ', n, n, n, -1 )
321 nb3 = ilaenv( 1, 'SORGQR', ' ', n, n, n, -1 )
322 nb = max( nb1, nb2, nb3 )
323 lopt = 2*n+n*(nb+1)
324 work( 1 ) = lopt
325 END IF
326*
327 IF( info.NE.0 ) THEN
328 CALL xerbla( 'SGEGS ', -info )
329 RETURN
330 ELSE IF( lquery ) THEN
331 RETURN
332 END IF
333*
334* Quick return if possible
335*
336 IF( n.EQ.0 )
337 $ RETURN
338*
339* Get machine constants
340*
341 eps = slamch( 'E' )*slamch( 'B' )
342 safmin = slamch( 'S' )
343 smlnum = n*safmin / eps
344 bignum = one / smlnum
345*
346* Scale A if max element outside range [SMLNUM,BIGNUM]
347*
348 anrm = slange( 'M', n, n, a, lda, work )
349 ilascl = .false.
350 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
351 anrmto = smlnum
352 ilascl = .true.
353 ELSE IF( anrm.GT.bignum ) THEN
354 anrmto = bignum
355 ilascl = .true.
356 END IF
357*
358 IF( ilascl ) THEN
359 CALL slascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda,
360 $ iinfo )
361 IF( iinfo.NE.0 ) THEN
362 info = n + 9
363 RETURN
364 END IF
365 END IF
366*
367* Scale B if max element outside range [SMLNUM,BIGNUM]
368*
369 bnrm = slange( 'M', n, n, b, ldb, work )
370 ilbscl = .false.
371 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
372 bnrmto = smlnum
373 ilbscl = .true.
374 ELSE IF( bnrm.GT.bignum ) THEN
375 bnrmto = bignum
376 ilbscl = .true.
377 END IF
378*
379 IF( ilbscl ) THEN
380 CALL slascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb,
381 $ 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,
500 $ iinfo )
501 IF( iinfo.NE.0 ) THEN
502 info = n + 9
503 RETURN
504 END IF
505 CALL slascl( 'G', -1, -1, anrmto, anrm, n, 1, alphar, n,
506 $ iinfo )
507 IF( iinfo.NE.0 ) THEN
508 info = n + 9
509 RETURN
510 END IF
511 CALL slascl( 'G', -1, -1, anrmto, anrm, n, 1, alphai, n,
512 $ iinfo )
513 IF( iinfo.NE.0 ) THEN
514 info = n + 9
515 RETURN
516 END IF
517 END IF
518*
519 IF( ilbscl ) THEN
520 CALL slascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb,
521 $ iinfo )
522 IF( iinfo.NE.0 ) THEN
523 info = n + 9
524 RETURN
525 END IF
526 CALL slascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n,
527 $ iinfo )
528 IF( iinfo.NE.0 ) THEN
529 info = n + 9
530 RETURN
531 END IF
532 END IF
533*
534 10 CONTINUE
535 work( 1 ) = lwkopt
536*
537 RETURN
538*
539* End of SGEGS
540*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:144
sggbak
subroutine sggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
SGGBAK
Definition sggbak.f:146
sggbal
subroutine sggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
SGGBAL
Definition sggbal.f:175
sgghrd
subroutine sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
Definition sgghrd.f:206
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:303
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
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:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:112
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:142
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:126
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166
Here is the call graph for this function: