You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

zgges3.f 19 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598
  1. *> \brief <b> ZGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)</b>
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download ZGGES3 + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgges3.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgges3.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgges3.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
  22. * $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
  23. * $ WORK, LWORK, RWORK, BWORK, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER JOBVSL, JOBVSR, SORT
  27. * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
  28. * ..
  29. * .. Array Arguments ..
  30. * LOGICAL BWORK( * )
  31. * DOUBLE PRECISION RWORK( * )
  32. * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
  33. * $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
  34. * $ WORK( * )
  35. * ..
  36. * .. Function Arguments ..
  37. * LOGICAL SELCTG
  38. * EXTERNAL SELCTG
  39. * ..
  40. *
  41. *
  42. *> \par Purpose:
  43. * =============
  44. *>
  45. *> \verbatim
  46. *>
  47. *> ZGGES3 computes for a pair of N-by-N complex nonsymmetric matrices
  48. *> (A,B), the generalized eigenvalues, the generalized complex Schur
  49. *> form (S, T), and optionally left and/or right Schur vectors (VSL
  50. *> and VSR). This gives the generalized Schur factorization
  51. *>
  52. *> (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
  53. *>
  54. *> where (VSR)**H is the conjugate-transpose of VSR.
  55. *>
  56. *> Optionally, it also orders the eigenvalues so that a selected cluster
  57. *> of eigenvalues appears in the leading diagonal blocks of the upper
  58. *> triangular matrix S and the upper triangular matrix T. The leading
  59. *> columns of VSL and VSR then form an unitary basis for the
  60. *> corresponding left and right eigenspaces (deflating subspaces).
  61. *>
  62. *> (If only the generalized eigenvalues are needed, use the driver
  63. *> ZGGEV instead, which is faster.)
  64. *>
  65. *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
  66. *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
  67. *> usually represented as the pair (alpha,beta), as there is a
  68. *> reasonable interpretation for beta=0, and even for both being zero.
  69. *>
  70. *> A pair of matrices (S,T) is in generalized complex Schur form if S
  71. *> and T are upper triangular and, in addition, the diagonal elements
  72. *> of T are non-negative real numbers.
  73. *> \endverbatim
  74. *
  75. * Arguments:
  76. * ==========
  77. *
  78. *> \param[in] JOBVSL
  79. *> \verbatim
  80. *> JOBVSL is CHARACTER*1
  81. *> = 'N': do not compute the left Schur vectors;
  82. *> = 'V': compute the left Schur vectors.
  83. *> \endverbatim
  84. *>
  85. *> \param[in] JOBVSR
  86. *> \verbatim
  87. *> JOBVSR is CHARACTER*1
  88. *> = 'N': do not compute the right Schur vectors;
  89. *> = 'V': compute the right Schur vectors.
  90. *> \endverbatim
  91. *>
  92. *> \param[in] SORT
  93. *> \verbatim
  94. *> SORT is CHARACTER*1
  95. *> Specifies whether or not to order the eigenvalues on the
  96. *> diagonal of the generalized Schur form.
  97. *> = 'N': Eigenvalues are not ordered;
  98. *> = 'S': Eigenvalues are ordered (see SELCTG).
  99. *> \endverbatim
  100. *>
  101. *> \param[in] SELCTG
  102. *> \verbatim
  103. *> SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments
  104. *> SELCTG must be declared EXTERNAL in the calling subroutine.
  105. *> If SORT = 'N', SELCTG is not referenced.
  106. *> If SORT = 'S', SELCTG is used to select eigenvalues to sort
  107. *> to the top left of the Schur form.
  108. *> An eigenvalue ALPHA(j)/BETA(j) is selected if
  109. *> SELCTG(ALPHA(j),BETA(j)) is true.
  110. *>
  111. *> Note that a selected complex eigenvalue may no longer satisfy
  112. *> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
  113. *> ordering may change the value of complex eigenvalues
  114. *> (especially if the eigenvalue is ill-conditioned), in this
  115. *> case INFO is set to N+2 (See INFO below).
  116. *> \endverbatim
  117. *>
  118. *> \param[in] N
  119. *> \verbatim
  120. *> N is INTEGER
  121. *> The order of the matrices A, B, VSL, and VSR. N >= 0.
  122. *> \endverbatim
  123. *>
  124. *> \param[in,out] A
  125. *> \verbatim
  126. *> A is COMPLEX*16 array, dimension (LDA, N)
  127. *> On entry, the first of the pair of matrices.
  128. *> On exit, A has been overwritten by its generalized Schur
  129. *> form S.
  130. *> \endverbatim
  131. *>
  132. *> \param[in] LDA
  133. *> \verbatim
  134. *> LDA is INTEGER
  135. *> The leading dimension of A. LDA >= max(1,N).
  136. *> \endverbatim
  137. *>
  138. *> \param[in,out] B
  139. *> \verbatim
  140. *> B is COMPLEX*16 array, dimension (LDB, N)
  141. *> On entry, the second of the pair of matrices.
  142. *> On exit, B has been overwritten by its generalized Schur
  143. *> form T.
  144. *> \endverbatim
  145. *>
  146. *> \param[in] LDB
  147. *> \verbatim
  148. *> LDB is INTEGER
  149. *> The leading dimension of B. LDB >= max(1,N).
  150. *> \endverbatim
  151. *>
  152. *> \param[out] SDIM
  153. *> \verbatim
  154. *> SDIM is INTEGER
  155. *> If SORT = 'N', SDIM = 0.
  156. *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
  157. *> for which SELCTG is true.
  158. *> \endverbatim
  159. *>
  160. *> \param[out] ALPHA
  161. *> \verbatim
  162. *> ALPHA is COMPLEX*16 array, dimension (N)
  163. *> \endverbatim
  164. *>
  165. *> \param[out] BETA
  166. *> \verbatim
  167. *> BETA is COMPLEX*16 array, dimension (N)
  168. *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
  169. *> generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
  170. *> j=1,...,N are the diagonals of the complex Schur form (A,B)
  171. *> output by ZGGES3. The BETA(j) will be non-negative real.
  172. *>
  173. *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
  174. *> underflow, and BETA(j) may even be zero. Thus, the user
  175. *> should avoid naively computing the ratio alpha/beta.
  176. *> However, ALPHA will be always less than and usually
  177. *> comparable with norm(A) in magnitude, and BETA always less
  178. *> than and usually comparable with norm(B).
  179. *> \endverbatim
  180. *>
  181. *> \param[out] VSL
  182. *> \verbatim
  183. *> VSL is COMPLEX*16 array, dimension (LDVSL,N)
  184. *> If JOBVSL = 'V', VSL will contain the left Schur vectors.
  185. *> Not referenced if JOBVSL = 'N'.
  186. *> \endverbatim
  187. *>
  188. *> \param[in] LDVSL
  189. *> \verbatim
  190. *> LDVSL is INTEGER
  191. *> The leading dimension of the matrix VSL. LDVSL >= 1, and
  192. *> if JOBVSL = 'V', LDVSL >= N.
  193. *> \endverbatim
  194. *>
  195. *> \param[out] VSR
  196. *> \verbatim
  197. *> VSR is COMPLEX*16 array, dimension (LDVSR,N)
  198. *> If JOBVSR = 'V', VSR will contain the right Schur vectors.
  199. *> Not referenced if JOBVSR = 'N'.
  200. *> \endverbatim
  201. *>
  202. *> \param[in] LDVSR
  203. *> \verbatim
  204. *> LDVSR is INTEGER
  205. *> The leading dimension of the matrix VSR. LDVSR >= 1, and
  206. *> if JOBVSR = 'V', LDVSR >= N.
  207. *> \endverbatim
  208. *>
  209. *> \param[out] WORK
  210. *> \verbatim
  211. *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
  212. *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  213. *> \endverbatim
  214. *>
  215. *> \param[in] LWORK
  216. *> \verbatim
  217. *> LWORK is INTEGER
  218. *> The dimension of the array WORK. LWORK >= MAX(1,2*N)
  219. *> For good performance, LWORK must generally be larger.
  220. *>
  221. *> If LWORK = -1, then a workspace query is assumed; the routine
  222. *> only calculates the optimal size of the WORK array, returns
  223. *> this value as the first entry of the WORK array, and no error
  224. *> message related to LWORK is issued by XERBLA.
  225. *> \endverbatim
  226. *>
  227. *> \param[out] RWORK
  228. *> \verbatim
  229. *> RWORK is DOUBLE PRECISION array, dimension (8*N)
  230. *> \endverbatim
  231. *>
  232. *> \param[out] BWORK
  233. *> \verbatim
  234. *> BWORK is LOGICAL array, dimension (N)
  235. *> Not referenced if SORT = 'N'.
  236. *> \endverbatim
  237. *>
  238. *> \param[out] INFO
  239. *> \verbatim
  240. *> INFO is INTEGER
  241. *> = 0: successful exit
  242. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  243. *> =1,...,N:
  244. *> The QZ iteration failed. (A,B) are not in Schur
  245. *> form, but ALPHA(j) and BETA(j) should be correct for
  246. *> j=INFO+1,...,N.
  247. *> > N: =N+1: other than QZ iteration failed in ZLAQZ0
  248. *> =N+2: after reordering, roundoff changed values of
  249. *> some complex eigenvalues so that leading
  250. *> eigenvalues in the Generalized Schur form no
  251. *> longer satisfy SELCTG=.TRUE. This could also
  252. *> be caused due to scaling.
  253. *> =N+3: reordering failed in ZTGSEN.
  254. *> \endverbatim
  255. *
  256. * Authors:
  257. * ========
  258. *
  259. *> \author Univ. of Tennessee
  260. *> \author Univ. of California Berkeley
  261. *> \author Univ. of Colorado Denver
  262. *> \author NAG Ltd.
  263. *
  264. *> \ingroup gges3
  265. *
  266. * =====================================================================
  267. SUBROUTINE ZGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
  268. $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
  269. $ WORK, LWORK, RWORK, BWORK, INFO )
  270. *
  271. * -- LAPACK driver routine --
  272. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  273. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  274. *
  275. * .. Scalar Arguments ..
  276. CHARACTER JOBVSL, JOBVSR, SORT
  277. INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
  278. * ..
  279. * .. Array Arguments ..
  280. LOGICAL BWORK( * )
  281. DOUBLE PRECISION RWORK( * )
  282. COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
  283. $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
  284. $ WORK( * )
  285. * ..
  286. * .. Function Arguments ..
  287. LOGICAL SELCTG
  288. EXTERNAL SELCTG
  289. * ..
  290. *
  291. * =====================================================================
  292. *
  293. * .. Parameters ..
  294. DOUBLE PRECISION ZERO, ONE
  295. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
  296. COMPLEX*16 CZERO, CONE
  297. PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
  298. $ CONE = ( 1.0D0, 0.0D0 ) )
  299. * ..
  300. * .. Local Scalars ..
  301. LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
  302. $ LQUERY, WANTST
  303. INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
  304. $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT,
  305. $ LWKMIN
  306. DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
  307. $ PVSR, SMLNUM
  308. * ..
  309. * .. Local Arrays ..
  310. INTEGER IDUM( 1 )
  311. DOUBLE PRECISION DIF( 2 )
  312. * ..
  313. * .. External Subroutines ..
  314. EXTERNAL XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHD3, ZLAQZ0,
  315. $ ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR, ZUNMQR
  316. * ..
  317. * .. External Functions ..
  318. LOGICAL LSAME
  319. DOUBLE PRECISION DLAMCH, ZLANGE
  320. EXTERNAL LSAME, DLAMCH, ZLANGE
  321. * ..
  322. * .. Intrinsic Functions ..
  323. INTRINSIC MAX, SQRT
  324. * ..
  325. * .. Executable Statements ..
  326. *
  327. * Decode the input arguments
  328. *
  329. IF( LSAME( JOBVSL, 'N' ) ) THEN
  330. IJOBVL = 1
  331. ILVSL = .FALSE.
  332. ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
  333. IJOBVL = 2
  334. ILVSL = .TRUE.
  335. ELSE
  336. IJOBVL = -1
  337. ILVSL = .FALSE.
  338. END IF
  339. *
  340. IF( LSAME( JOBVSR, 'N' ) ) THEN
  341. IJOBVR = 1
  342. ILVSR = .FALSE.
  343. ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
  344. IJOBVR = 2
  345. ILVSR = .TRUE.
  346. ELSE
  347. IJOBVR = -1
  348. ILVSR = .FALSE.
  349. END IF
  350. *
  351. WANTST = LSAME( SORT, 'S' )
  352. *
  353. * Test the input arguments
  354. *
  355. INFO = 0
  356. LQUERY = ( LWORK.EQ.-1 )
  357. LWKMIN = MAX( 1, 2*N )
  358. *
  359. IF( IJOBVL.LE.0 ) THEN
  360. INFO = -1
  361. ELSE IF( IJOBVR.LE.0 ) THEN
  362. INFO = -2
  363. ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
  364. INFO = -3
  365. ELSE IF( N.LT.0 ) THEN
  366. INFO = -5
  367. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  368. INFO = -7
  369. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  370. INFO = -9
  371. ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
  372. INFO = -14
  373. ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
  374. INFO = -16
  375. ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
  376. INFO = -18
  377. END IF
  378. *
  379. * Compute workspace
  380. *
  381. IF( INFO.EQ.0 ) THEN
  382. CALL ZGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
  383. LWKOPT = MAX( LWKMIN, N + INT( WORK( 1 ) ) )
  384. CALL ZUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
  385. $ -1, IERR )
  386. LWKOPT = MAX( LWKOPT, N + INT( WORK( 1 ) ) )
  387. IF( ILVSL ) THEN
  388. CALL ZUNGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1, IERR )
  389. LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
  390. END IF
  391. CALL ZGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
  392. $ LDVSL, VSR, LDVSR, WORK, -1, IERR )
  393. LWKOPT = MAX( LWKOPT, N + INT( WORK( 1 ) ) )
  394. CALL ZLAQZ0( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
  395. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, -1,
  396. $ RWORK, 0, IERR )
  397. LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
  398. IF( WANTST ) THEN
  399. CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
  400. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM,
  401. $ PVSL, PVSR, DIF, WORK, -1, IDUM, 1, IERR )
  402. LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
  403. END IF
  404. IF( N.EQ.0 ) THEN
  405. WORK( 1 ) = 1
  406. ELSE
  407. WORK( 1 ) = DCMPLX( LWKOPT )
  408. END IF
  409. END IF
  410. *
  411. IF( INFO.NE.0 ) THEN
  412. CALL XERBLA( 'ZGGES3 ', -INFO )
  413. RETURN
  414. ELSE IF( LQUERY ) THEN
  415. RETURN
  416. END IF
  417. *
  418. * Quick return if possible
  419. *
  420. IF( N.EQ.0 ) THEN
  421. SDIM = 0
  422. RETURN
  423. END IF
  424. *
  425. * Get machine constants
  426. *
  427. EPS = DLAMCH( 'P' )
  428. SMLNUM = DLAMCH( 'S' )
  429. BIGNUM = ONE / SMLNUM
  430. SMLNUM = SQRT( SMLNUM ) / EPS
  431. BIGNUM = ONE / SMLNUM
  432. *
  433. * Scale A if max element outside range [SMLNUM,BIGNUM]
  434. *
  435. ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
  436. ILASCL = .FALSE.
  437. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  438. ANRMTO = SMLNUM
  439. ILASCL = .TRUE.
  440. ELSE IF( ANRM.GT.BIGNUM ) THEN
  441. ANRMTO = BIGNUM
  442. ILASCL = .TRUE.
  443. END IF
  444. *
  445. IF( ILASCL )
  446. $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  447. *
  448. * Scale B if max element outside range [SMLNUM,BIGNUM]
  449. *
  450. BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
  451. ILBSCL = .FALSE.
  452. IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  453. BNRMTO = SMLNUM
  454. ILBSCL = .TRUE.
  455. ELSE IF( BNRM.GT.BIGNUM ) THEN
  456. BNRMTO = BIGNUM
  457. ILBSCL = .TRUE.
  458. END IF
  459. *
  460. IF( ILBSCL )
  461. $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  462. *
  463. * Permute the matrix to make it more nearly triangular
  464. *
  465. ILEFT = 1
  466. IRIGHT = N + 1
  467. IRWRK = IRIGHT + N
  468. CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
  469. $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
  470. *
  471. * Reduce B to triangular form (QR decomposition of B)
  472. *
  473. IROWS = IHI + 1 - ILO
  474. ICOLS = N + 1 - ILO
  475. ITAU = 1
  476. IWRK = ITAU + IROWS
  477. CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  478. $ WORK( IWRK ), LWORK+1-IWRK, IERR )
  479. *
  480. * Apply the orthogonal transformation to matrix A
  481. *
  482. CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
  483. $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
  484. $ LWORK+1-IWRK, IERR )
  485. *
  486. * Initialize VSL
  487. *
  488. IF( ILVSL ) THEN
  489. CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
  490. IF( IROWS.GT.1 ) THEN
  491. CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  492. $ VSL( ILO+1, ILO ), LDVSL )
  493. END IF
  494. CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
  495. $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  496. END IF
  497. *
  498. * Initialize VSR
  499. *
  500. IF( ILVSR )
  501. $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
  502. *
  503. * Reduce to generalized Hessenberg form
  504. *
  505. CALL ZGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
  506. $ LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK, IERR )
  507. *
  508. SDIM = 0
  509. *
  510. * Perform QZ algorithm, computing Schur vectors if desired
  511. *
  512. IWRK = ITAU
  513. CALL ZLAQZ0( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
  514. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
  515. $ LWORK+1-IWRK, RWORK( IRWRK ), 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 30
  525. END IF
  526. *
  527. * Sort eigenvalues ALPHA/BETA if desired
  528. *
  529. IF( WANTST ) THEN
  530. *
  531. * Undo scaling on eigenvalues before selecting
  532. *
  533. IF( ILASCL )
  534. $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
  535. IF( ILBSCL )
  536. $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
  537. *
  538. * Select eigenvalues
  539. *
  540. DO 10 I = 1, N
  541. BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
  542. 10 CONTINUE
  543. *
  544. CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
  545. $ BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
  546. $ DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
  547. IF( IERR.EQ.1 )
  548. $ INFO = N + 3
  549. *
  550. END IF
  551. *
  552. * Apply back-permutation to VSL and VSR
  553. *
  554. IF( ILVSL )
  555. $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
  556. $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
  557. IF( ILVSR )
  558. $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
  559. $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
  560. *
  561. * Undo scaling
  562. *
  563. IF( ILASCL ) THEN
  564. CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
  565. CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
  566. END IF
  567. *
  568. IF( ILBSCL ) THEN
  569. CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
  570. CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  571. END IF
  572. *
  573. IF( WANTST ) THEN
  574. *
  575. * Check if reordering is correct
  576. *
  577. LASTSL = .TRUE.
  578. SDIM = 0
  579. DO 20 I = 1, N
  580. CURSL = SELCTG( ALPHA( I ), BETA( I ) )
  581. IF( CURSL )
  582. $ SDIM = SDIM + 1
  583. IF( CURSL .AND. .NOT.LASTSL )
  584. $ INFO = N + 2
  585. LASTSL = CURSL
  586. 20 CONTINUE
  587. *
  588. END IF
  589. *
  590. 30 CONTINUE
  591. *
  592. WORK( 1 ) = DCMPLX( LWKOPT )
  593. *
  594. RETURN
  595. *
  596. * End of ZGGES3
  597. *
  598. END