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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592
  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.
  219. *>
  220. *> If LWORK = -1, then a workspace query is assumed; the routine
  221. *> only calculates the optimal size of the WORK array, returns
  222. *> this value as the first entry of the WORK array, and no error
  223. *> message related to LWORK is issued by XERBLA.
  224. *> \endverbatim
  225. *>
  226. *> \param[out] RWORK
  227. *> \verbatim
  228. *> RWORK is DOUBLE PRECISION array, dimension (8*N)
  229. *> \endverbatim
  230. *>
  231. *> \param[out] BWORK
  232. *> \verbatim
  233. *> BWORK is LOGICAL array, dimension (N)
  234. *> Not referenced if SORT = 'N'.
  235. *> \endverbatim
  236. *>
  237. *> \param[out] INFO
  238. *> \verbatim
  239. *> INFO is INTEGER
  240. *> = 0: successful exit
  241. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  242. *> =1,...,N:
  243. *> The QZ iteration failed. (A,B) are not in Schur
  244. *> form, but ALPHA(j) and BETA(j) should be correct for
  245. *> j=INFO+1,...,N.
  246. *> > N: =N+1: other than QZ iteration failed in ZLAQZ0
  247. *> =N+2: after reordering, roundoff changed values of
  248. *> some complex eigenvalues so that leading
  249. *> eigenvalues in the Generalized Schur form no
  250. *> longer satisfy SELCTG=.TRUE. This could also
  251. *> be caused due to scaling.
  252. *> =N+3: reordering failed in ZTGSEN.
  253. *> \endverbatim
  254. *
  255. * Authors:
  256. * ========
  257. *
  258. *> \author Univ. of Tennessee
  259. *> \author Univ. of California Berkeley
  260. *> \author Univ. of Colorado Denver
  261. *> \author NAG Ltd.
  262. *
  263. *> \ingroup complex16GEeigen
  264. *
  265. * =====================================================================
  266. SUBROUTINE ZGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
  267. $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
  268. $ WORK, LWORK, RWORK, BWORK, INFO )
  269. *
  270. * -- LAPACK driver routine --
  271. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  272. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  273. *
  274. * .. Scalar Arguments ..
  275. CHARACTER JOBVSL, JOBVSR, SORT
  276. INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
  277. * ..
  278. * .. Array Arguments ..
  279. LOGICAL BWORK( * )
  280. DOUBLE PRECISION RWORK( * )
  281. COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
  282. $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
  283. $ WORK( * )
  284. * ..
  285. * .. Function Arguments ..
  286. LOGICAL SELCTG
  287. EXTERNAL SELCTG
  288. * ..
  289. *
  290. * =====================================================================
  291. *
  292. * .. Parameters ..
  293. DOUBLE PRECISION ZERO, ONE
  294. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
  295. COMPLEX*16 CZERO, CONE
  296. PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
  297. $ CONE = ( 1.0D0, 0.0D0 ) )
  298. * ..
  299. * .. Local Scalars ..
  300. LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
  301. $ LQUERY, WANTST
  302. INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
  303. $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT
  304. DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
  305. $ PVSR, SMLNUM
  306. * ..
  307. * .. Local Arrays ..
  308. INTEGER IDUM( 1 )
  309. DOUBLE PRECISION DIF( 2 )
  310. * ..
  311. * .. External Subroutines ..
  312. EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHD3,
  313. $ ZLAQZ0, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
  314. $ ZUNMQR
  315. * ..
  316. * .. External Functions ..
  317. LOGICAL LSAME
  318. DOUBLE PRECISION DLAMCH, ZLANGE
  319. EXTERNAL LSAME, DLAMCH, ZLANGE
  320. * ..
  321. * .. Intrinsic Functions ..
  322. INTRINSIC MAX, SQRT
  323. * ..
  324. * .. Executable Statements ..
  325. *
  326. * Decode the input arguments
  327. *
  328. IF( LSAME( JOBVSL, 'N' ) ) THEN
  329. IJOBVL = 1
  330. ILVSL = .FALSE.
  331. ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
  332. IJOBVL = 2
  333. ILVSL = .TRUE.
  334. ELSE
  335. IJOBVL = -1
  336. ILVSL = .FALSE.
  337. END IF
  338. *
  339. IF( LSAME( JOBVSR, 'N' ) ) THEN
  340. IJOBVR = 1
  341. ILVSR = .FALSE.
  342. ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
  343. IJOBVR = 2
  344. ILVSR = .TRUE.
  345. ELSE
  346. IJOBVR = -1
  347. ILVSR = .FALSE.
  348. END IF
  349. *
  350. WANTST = LSAME( SORT, 'S' )
  351. *
  352. * Test the input arguments
  353. *
  354. INFO = 0
  355. LQUERY = ( LWORK.EQ.-1 )
  356. IF( IJOBVL.LE.0 ) THEN
  357. INFO = -1
  358. ELSE IF( IJOBVR.LE.0 ) THEN
  359. INFO = -2
  360. ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
  361. INFO = -3
  362. ELSE IF( N.LT.0 ) THEN
  363. INFO = -5
  364. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  365. INFO = -7
  366. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  367. INFO = -9
  368. ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
  369. INFO = -14
  370. ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
  371. INFO = -16
  372. ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
  373. INFO = -18
  374. END IF
  375. *
  376. * Compute workspace
  377. *
  378. IF( INFO.EQ.0 ) THEN
  379. CALL ZGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
  380. LWKOPT = MAX( 1, N + INT ( WORK( 1 ) ) )
  381. CALL ZUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
  382. $ -1, IERR )
  383. LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
  384. IF( ILVSL ) THEN
  385. CALL ZUNGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1, IERR )
  386. LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
  387. END IF
  388. CALL ZGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
  389. $ LDVSL, VSR, LDVSR, WORK, -1, IERR )
  390. LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
  391. CALL ZLAQZ0( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
  392. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, -1,
  393. $ RWORK, 0, IERR )
  394. LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
  395. IF( WANTST ) THEN
  396. CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
  397. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM,
  398. $ PVSL, PVSR, DIF, WORK, -1, IDUM, 1, IERR )
  399. LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
  400. END IF
  401. WORK( 1 ) = DCMPLX( LWKOPT )
  402. END IF
  403. *
  404. IF( INFO.NE.0 ) THEN
  405. CALL XERBLA( 'ZGGES3 ', -INFO )
  406. RETURN
  407. ELSE IF( LQUERY ) THEN
  408. RETURN
  409. END IF
  410. *
  411. * Quick return if possible
  412. *
  413. IF( N.EQ.0 ) THEN
  414. SDIM = 0
  415. RETURN
  416. END IF
  417. *
  418. * Get machine constants
  419. *
  420. EPS = DLAMCH( 'P' )
  421. SMLNUM = DLAMCH( 'S' )
  422. BIGNUM = ONE / SMLNUM
  423. CALL DLABAD( SMLNUM, BIGNUM )
  424. SMLNUM = SQRT( SMLNUM ) / EPS
  425. BIGNUM = ONE / SMLNUM
  426. *
  427. * Scale A if max element outside range [SMLNUM,BIGNUM]
  428. *
  429. ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
  430. ILASCL = .FALSE.
  431. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  432. ANRMTO = SMLNUM
  433. ILASCL = .TRUE.
  434. ELSE IF( ANRM.GT.BIGNUM ) THEN
  435. ANRMTO = BIGNUM
  436. ILASCL = .TRUE.
  437. END IF
  438. *
  439. IF( ILASCL )
  440. $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  441. *
  442. * Scale B if max element outside range [SMLNUM,BIGNUM]
  443. *
  444. BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
  445. ILBSCL = .FALSE.
  446. IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  447. BNRMTO = SMLNUM
  448. ILBSCL = .TRUE.
  449. ELSE IF( BNRM.GT.BIGNUM ) THEN
  450. BNRMTO = BIGNUM
  451. ILBSCL = .TRUE.
  452. END IF
  453. *
  454. IF( ILBSCL )
  455. $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  456. *
  457. * Permute the matrix to make it more nearly triangular
  458. *
  459. ILEFT = 1
  460. IRIGHT = N + 1
  461. IRWRK = IRIGHT + N
  462. CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
  463. $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
  464. *
  465. * Reduce B to triangular form (QR decomposition of B)
  466. *
  467. IROWS = IHI + 1 - ILO
  468. ICOLS = N + 1 - ILO
  469. ITAU = 1
  470. IWRK = ITAU + IROWS
  471. CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  472. $ WORK( IWRK ), LWORK+1-IWRK, IERR )
  473. *
  474. * Apply the orthogonal transformation to matrix A
  475. *
  476. CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
  477. $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
  478. $ LWORK+1-IWRK, IERR )
  479. *
  480. * Initialize VSL
  481. *
  482. IF( ILVSL ) THEN
  483. CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
  484. IF( IROWS.GT.1 ) THEN
  485. CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  486. $ VSL( ILO+1, ILO ), LDVSL )
  487. END IF
  488. CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
  489. $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  490. END IF
  491. *
  492. * Initialize VSR
  493. *
  494. IF( ILVSR )
  495. $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
  496. *
  497. * Reduce to generalized Hessenberg form
  498. *
  499. CALL ZGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
  500. $ LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK, IERR )
  501. *
  502. SDIM = 0
  503. *
  504. * Perform QZ algorithm, computing Schur vectors if desired
  505. *
  506. IWRK = ITAU
  507. CALL ZLAQZ0( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
  508. $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
  509. $ LWORK+1-IWRK, RWORK( IRWRK ), 0, IERR )
  510. IF( IERR.NE.0 ) THEN
  511. IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
  512. INFO = IERR
  513. ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
  514. INFO = IERR - N
  515. ELSE
  516. INFO = N + 1
  517. END IF
  518. GO TO 30
  519. END IF
  520. *
  521. * Sort eigenvalues ALPHA/BETA if desired
  522. *
  523. IF( WANTST ) THEN
  524. *
  525. * Undo scaling on eigenvalues before selecting
  526. *
  527. IF( ILASCL )
  528. $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
  529. IF( ILBSCL )
  530. $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
  531. *
  532. * Select eigenvalues
  533. *
  534. DO 10 I = 1, N
  535. BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
  536. 10 CONTINUE
  537. *
  538. CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
  539. $ BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
  540. $ DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
  541. IF( IERR.EQ.1 )
  542. $ INFO = N + 3
  543. *
  544. END IF
  545. *
  546. * Apply back-permutation to VSL and VSR
  547. *
  548. IF( ILVSL )
  549. $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
  550. $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
  551. IF( ILVSR )
  552. $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
  553. $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
  554. *
  555. * Undo scaling
  556. *
  557. IF( ILASCL ) THEN
  558. CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
  559. CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
  560. END IF
  561. *
  562. IF( ILBSCL ) THEN
  563. CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
  564. CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  565. END IF
  566. *
  567. IF( WANTST ) THEN
  568. *
  569. * Check if reordering is correct
  570. *
  571. LASTSL = .TRUE.
  572. SDIM = 0
  573. DO 20 I = 1, N
  574. CURSL = SELCTG( ALPHA( I ), BETA( I ) )
  575. IF( CURSL )
  576. $ SDIM = SDIM + 1
  577. IF( CURSL .AND. .NOT.LASTSL )
  578. $ INFO = N + 2
  579. LASTSL = CURSL
  580. 20 CONTINUE
  581. *
  582. END IF
  583. *
  584. 30 CONTINUE
  585. *
  586. WORK( 1 ) = DCMPLX( LWKOPT )
  587. *
  588. RETURN
  589. *
  590. * End of ZGGES3
  591. *
  592. END