PageRenderTime 47ms CodeModel.GetById 20ms RepoModel.GetById 0ms app.codeStats 0ms

/lapack/double/dggevx.f

http://github.com/pymc-devs/pymc
FORTRAN Legacy | 718 lines | 320 code | 0 blank | 398 comment | 0 complexity | bc76a1ae782b553ed7cee0f8de7911cb MD5 | raw file
Possible License(s): CC-BY-SA-3.0
  1. SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
  2. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
  3. $ IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
  4. $ RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
  5. *
  6. * -- LAPACK driver routine (version 3.1) --
  7. * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  8. * November 2006
  9. *
  10. * .. Scalar Arguments ..
  11. CHARACTER BALANC, JOBVL, JOBVR, SENSE
  12. INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  13. DOUBLE PRECISION ABNRM, BBNRM
  14. * ..
  15. * .. Array Arguments ..
  16. LOGICAL BWORK( * )
  17. INTEGER IWORK( * )
  18. DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  19. $ B( LDB, * ), BETA( * ), LSCALE( * ),
  20. $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
  21. $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
  22. * ..
  23. *
  24. * Purpose
  25. * =======
  26. *
  27. * DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
  28. * the generalized eigenvalues, and optionally, the left and/or right
  29. * generalized eigenvectors.
  30. *
  31. * Optionally also, it computes a balancing transformation to improve
  32. * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
  33. * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
  34. * the eigenvalues (RCONDE), and reciprocal condition numbers for the
  35. * right eigenvectors (RCONDV).
  36. *
  37. * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
  38. * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
  39. * singular. It is usually represented as the pair (alpha,beta), as
  40. * there is a reasonable interpretation for beta=0, and even for both
  41. * being zero.
  42. *
  43. * The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
  44. * of (A,B) satisfies
  45. *
  46. * A * v(j) = lambda(j) * B * v(j) .
  47. *
  48. * The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
  49. * of (A,B) satisfies
  50. *
  51. * u(j)**H * A = lambda(j) * u(j)**H * B.
  52. *
  53. * where u(j)**H is the conjugate-transpose of u(j).
  54. *
  55. *
  56. * Arguments
  57. * =========
  58. *
  59. * BALANC (input) CHARACTER*1
  60. * Specifies the balance option to be performed.
  61. * = 'N': do not diagonally scale or permute;
  62. * = 'P': permute only;
  63. * = 'S': scale only;
  64. * = 'B': both permute and scale.
  65. * Computed reciprocal condition numbers will be for the
  66. * matrices after permuting and/or balancing. Permuting does
  67. * not change condition numbers (in exact arithmetic), but
  68. * balancing does.
  69. *
  70. * JOBVL (input) CHARACTER*1
  71. * = 'N': do not compute the left generalized eigenvectors;
  72. * = 'V': compute the left generalized eigenvectors.
  73. *
  74. * JOBVR (input) CHARACTER*1
  75. * = 'N': do not compute the right generalized eigenvectors;
  76. * = 'V': compute the right generalized eigenvectors.
  77. *
  78. * SENSE (input) CHARACTER*1
  79. * Determines which reciprocal condition numbers are computed.
  80. * = 'N': none are computed;
  81. * = 'E': computed for eigenvalues only;
  82. * = 'V': computed for eigenvectors only;
  83. * = 'B': computed for eigenvalues and eigenvectors.
  84. *
  85. * N (input) INTEGER
  86. * The order of the matrices A, B, VL, and VR. N >= 0.
  87. *
  88. * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  89. * On entry, the matrix A in the pair (A,B).
  90. * On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
  91. * or both, then A contains the first part of the real Schur
  92. * form of the "balanced" versions of the input A and B.
  93. *
  94. * LDA (input) INTEGER
  95. * The leading dimension of A. LDA >= max(1,N).
  96. *
  97. * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
  98. * On entry, the matrix B in the pair (A,B).
  99. * On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
  100. * or both, then B contains the second part of the real Schur
  101. * form of the "balanced" versions of the input A and B.
  102. *
  103. * LDB (input) INTEGER
  104. * The leading dimension of B. LDB >= max(1,N).
  105. *
  106. * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
  107. * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
  108. * BETA (output) DOUBLE PRECISION array, dimension (N)
  109. * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
  110. * be the generalized eigenvalues. If ALPHAI(j) is zero, then
  111. * the j-th eigenvalue is real; if positive, then the j-th and
  112. * (j+1)-st eigenvalues are a complex conjugate pair, with
  113. * ALPHAI(j+1) negative.
  114. *
  115. * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
  116. * may easily over- or underflow, and BETA(j) may even be zero.
  117. * Thus, the user should avoid naively computing the ratio
  118. * ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
  119. * than and usually comparable with norm(A) in magnitude, and
  120. * BETA always less than and usually comparable with norm(B).
  121. *
  122. * VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
  123. * If JOBVL = 'V', the left eigenvectors u(j) are stored one
  124. * after another in the columns of VL, in the same order as
  125. * their eigenvalues. If the j-th eigenvalue is real, then
  126. * u(j) = VL(:,j), the j-th column of VL. If the j-th and
  127. * (j+1)-th eigenvalues form a complex conjugate pair, then
  128. * u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
  129. * Each eigenvector will be scaled so the largest component have
  130. * abs(real part) + abs(imag. part) = 1.
  131. * Not referenced if JOBVL = 'N'.
  132. *
  133. * LDVL (input) INTEGER
  134. * The leading dimension of the matrix VL. LDVL >= 1, and
  135. * if JOBVL = 'V', LDVL >= N.
  136. *
  137. * VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
  138. * If JOBVR = 'V', the right eigenvectors v(j) are stored one
  139. * after another in the columns of VR, in the same order as
  140. * their eigenvalues. If the j-th eigenvalue is real, then
  141. * v(j) = VR(:,j), the j-th column of VR. If the j-th and
  142. * (j+1)-th eigenvalues form a complex conjugate pair, then
  143. * v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
  144. * Each eigenvector will be scaled so the largest component have
  145. * abs(real part) + abs(imag. part) = 1.
  146. * Not referenced if JOBVR = 'N'.
  147. *
  148. * LDVR (input) INTEGER
  149. * The leading dimension of the matrix VR. LDVR >= 1, and
  150. * if JOBVR = 'V', LDVR >= N.
  151. *
  152. * ILO (output) INTEGER
  153. * IHI (output) INTEGER
  154. * ILO and IHI are integer values such that on exit
  155. * A(i,j) = 0 and B(i,j) = 0 if i > j and
  156. * j = 1,...,ILO-1 or i = IHI+1,...,N.
  157. * If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
  158. *
  159. * LSCALE (output) DOUBLE PRECISION array, dimension (N)
  160. * Details of the permutations and scaling factors applied
  161. * to the left side of A and B. If PL(j) is the index of the
  162. * row interchanged with row j, and DL(j) is the scaling
  163. * factor applied to row j, then
  164. * LSCALE(j) = PL(j) for j = 1,...,ILO-1
  165. * = DL(j) for j = ILO,...,IHI
  166. * = PL(j) for j = IHI+1,...,N.
  167. * The order in which the interchanges are made is N to IHI+1,
  168. * then 1 to ILO-1.
  169. *
  170. * RSCALE (output) DOUBLE PRECISION array, dimension (N)
  171. * Details of the permutations and scaling factors applied
  172. * to the right side of A and B. If PR(j) is the index of the
  173. * column interchanged with column j, and DR(j) is the scaling
  174. * factor applied to column j, then
  175. * RSCALE(j) = PR(j) for j = 1,...,ILO-1
  176. * = DR(j) for j = ILO,...,IHI
  177. * = PR(j) for j = IHI+1,...,N
  178. * The order in which the interchanges are made is N to IHI+1,
  179. * then 1 to ILO-1.
  180. *
  181. * ABNRM (output) DOUBLE PRECISION
  182. * The one-norm of the balanced matrix A.
  183. *
  184. * BBNRM (output) DOUBLE PRECISION
  185. * The one-norm of the balanced matrix B.
  186. *
  187. * RCONDE (output) DOUBLE PRECISION array, dimension (N)
  188. * If SENSE = 'E' or 'B', the reciprocal condition numbers of
  189. * the eigenvalues, stored in consecutive elements of the array.
  190. * For a complex conjugate pair of eigenvalues two consecutive
  191. * elements of RCONDE are set to the same value. Thus RCONDE(j),
  192. * RCONDV(j), and the j-th columns of VL and VR all correspond
  193. * to the j-th eigenpair.
  194. * If SENSE = 'N or 'V', RCONDE is not referenced.
  195. *
  196. * RCONDV (output) DOUBLE PRECISION array, dimension (N)
  197. * If SENSE = 'V' or 'B', the estimated reciprocal condition
  198. * numbers of the eigenvectors, stored in consecutive elements
  199. * of the array. For a complex eigenvector two consecutive
  200. * elements of RCONDV are set to the same value. If the
  201. * eigenvalues cannot be reordered to compute RCONDV(j),
  202. * RCONDV(j) is set to 0; this can only occur when the true
  203. * value would be very small anyway.
  204. * If SENSE = 'N' or 'E', RCONDV is not referenced.
  205. *
  206. * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  207. * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  208. *
  209. * LWORK (input) INTEGER
  210. * The dimension of the array WORK. LWORK >= max(1,2*N).
  211. * If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
  212. * LWORK >= max(1,6*N).
  213. * If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
  214. * If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
  215. *
  216. * If LWORK = -1, then a workspace query is assumed; the routine
  217. * only calculates the optimal size of the WORK array, returns
  218. * this value as the first entry of the WORK array, and no error
  219. * message related to LWORK is issued by XERBLA.
  220. *
  221. * IWORK (workspace) INTEGER array, dimension (N+6)
  222. * If SENSE = 'E', IWORK is not referenced.
  223. *
  224. * BWORK (workspace) LOGICAL array, dimension (N)
  225. * If SENSE = 'N', BWORK is not referenced.
  226. *
  227. * INFO (output) INTEGER
  228. * = 0: successful exit
  229. * < 0: if INFO = -i, the i-th argument had an illegal value.
  230. * = 1,...,N:
  231. * The QZ iteration failed. No eigenvectors have been
  232. * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
  233. * should be correct for j=INFO+1,...,N.
  234. * > N: =N+1: other than QZ iteration failed in DHGEQZ.
  235. * =N+2: error return from DTGEVC.
  236. *
  237. * Further Details
  238. * ===============
  239. *
  240. * Balancing a matrix pair (A,B) includes, first, permuting rows and
  241. * columns to isolate eigenvalues, second, applying diagonal similarity
  242. * transformation to the rows and columns to make the rows and columns
  243. * as close in norm as possible. The computed reciprocal condition
  244. * numbers correspond to the balanced matrix. Permuting rows and columns
  245. * will not change the condition numbers (in exact arithmetic) but
  246. * diagonal scaling will. For further explanation of balancing, see
  247. * section 4.11.1.2 of LAPACK Users' Guide.
  248. *
  249. * An approximate error bound on the chordal distance between the i-th
  250. * computed generalized eigenvalue w and the corresponding exact
  251. * eigenvalue lambda is
  252. *
  253. * chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
  254. *
  255. * An approximate error bound for the angle between the i-th computed
  256. * eigenvector VL(i) or VR(i) is given by
  257. *
  258. * EPS * norm(ABNRM, BBNRM) / DIF(i).
  259. *
  260. * For further explanation of the reciprocal condition numbers RCONDE
  261. * and RCONDV, see section 4.11 of LAPACK User's Guide.
  262. *
  263. * =====================================================================
  264. *
  265. * .. Parameters ..
  266. DOUBLE PRECISION ZERO, ONE
  267. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  268. * ..
  269. * .. Local Scalars ..
  270. LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
  271. $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
  272. CHARACTER CHTEMP
  273. INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
  274. $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
  275. $ MINWRK, MM
  276. DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
  277. $ SMLNUM, TEMP
  278. * ..
  279. * .. Local Arrays ..
  280. LOGICAL LDUMMA( 1 )
  281. * ..
  282. * .. External Subroutines ..
  283. EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
  284. $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
  285. $ DTGSNA, XERBLA
  286. * ..
  287. * .. External Functions ..
  288. LOGICAL LSAME
  289. INTEGER ILAENV
  290. DOUBLE PRECISION DLAMCH, DLANGE
  291. EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
  292. * ..
  293. * .. Intrinsic Functions ..
  294. INTRINSIC ABS, MAX, SQRT
  295. * ..
  296. * .. Executable Statements ..
  297. *
  298. * Decode the input arguments
  299. *
  300. IF( LSAME( JOBVL, 'N' ) ) THEN
  301. IJOBVL = 1
  302. ILVL = .FALSE.
  303. ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
  304. IJOBVL = 2
  305. ILVL = .TRUE.
  306. ELSE
  307. IJOBVL = -1
  308. ILVL = .FALSE.
  309. END IF
  310. *
  311. IF( LSAME( JOBVR, 'N' ) ) THEN
  312. IJOBVR = 1
  313. ILVR = .FALSE.
  314. ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
  315. IJOBVR = 2
  316. ILVR = .TRUE.
  317. ELSE
  318. IJOBVR = -1
  319. ILVR = .FALSE.
  320. END IF
  321. ILV = ILVL .OR. ILVR
  322. *
  323. NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
  324. WANTSN = LSAME( SENSE, 'N' )
  325. WANTSE = LSAME( SENSE, 'E' )
  326. WANTSV = LSAME( SENSE, 'V' )
  327. WANTSB = LSAME( SENSE, 'B' )
  328. *
  329. * Test the input arguments
  330. *
  331. INFO = 0
  332. LQUERY = ( LWORK.EQ.-1 )
  333. IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
  334. $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
  335. $ THEN
  336. INFO = -1
  337. ELSE IF( IJOBVL.LE.0 ) THEN
  338. INFO = -2
  339. ELSE IF( IJOBVR.LE.0 ) THEN
  340. INFO = -3
  341. ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
  342. $ THEN
  343. INFO = -4
  344. ELSE IF( N.LT.0 ) THEN
  345. INFO = -5
  346. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  347. INFO = -7
  348. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  349. INFO = -9
  350. ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
  351. INFO = -14
  352. ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
  353. INFO = -16
  354. END IF
  355. *
  356. * Compute workspace
  357. * (Note: Comments in the code beginning "Workspace:" describe the
  358. * minimal amount of workspace needed at that point in the code,
  359. * as well as the preferred amount for good performance.
  360. * NB refers to the optimal block size for the immediately
  361. * following subroutine, as returned by ILAENV. The workspace is
  362. * computed assuming ILO = 1 and IHI = N, the worst case.)
  363. *
  364. IF( INFO.EQ.0 ) THEN
  365. IF( N.EQ.0 ) THEN
  366. MINWRK = 1
  367. MAXWRK = 1
  368. ELSE
  369. IF( NOSCL .AND. .NOT.ILV ) THEN
  370. MINWRK = 2*N
  371. ELSE
  372. MINWRK = 6*N
  373. END IF
  374. IF( WANTSE .OR. WANTSB ) THEN
  375. MINWRK = 10*N
  376. END IF
  377. IF( WANTSV .OR. WANTSB ) THEN
  378. MINWRK = MAX( MINWRK, 2*N*( N + 4 ) + 16 )
  379. END IF
  380. MAXWRK = MINWRK
  381. MAXWRK = MAX( MAXWRK,
  382. $ N + N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) )
  383. MAXWRK = MAX( MAXWRK,
  384. $ N + N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) )
  385. IF( ILVL ) THEN
  386. MAXWRK = MAX( MAXWRK, N +
  387. $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, 0 ) )
  388. END IF
  389. END IF
  390. WORK( 1 ) = MAXWRK
  391. *
  392. IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
  393. INFO = -26
  394. END IF
  395. END IF
  396. *
  397. IF( INFO.NE.0 ) THEN
  398. CALL XERBLA( 'DGGEVX', -INFO )
  399. RETURN
  400. ELSE IF( LQUERY ) THEN
  401. RETURN
  402. END IF
  403. *
  404. * Quick return if possible
  405. *
  406. IF( N.EQ.0 )
  407. $ RETURN
  408. *
  409. *
  410. * Get machine constants
  411. *
  412. EPS = DLAMCH( 'P' )
  413. SMLNUM = DLAMCH( 'S' )
  414. BIGNUM = ONE / SMLNUM
  415. CALL DLABAD( SMLNUM, BIGNUM )
  416. SMLNUM = SQRT( SMLNUM ) / EPS
  417. BIGNUM = ONE / SMLNUM
  418. *
  419. * Scale A if max element outside range [SMLNUM,BIGNUM]
  420. *
  421. ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
  422. ILASCL = .FALSE.
  423. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  424. ANRMTO = SMLNUM
  425. ILASCL = .TRUE.
  426. ELSE IF( ANRM.GT.BIGNUM ) THEN
  427. ANRMTO = BIGNUM
  428. ILASCL = .TRUE.
  429. END IF
  430. IF( ILASCL )
  431. $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  432. *
  433. * Scale B if max element outside range [SMLNUM,BIGNUM]
  434. *
  435. BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
  436. ILBSCL = .FALSE.
  437. IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  438. BNRMTO = SMLNUM
  439. ILBSCL = .TRUE.
  440. ELSE IF( BNRM.GT.BIGNUM ) THEN
  441. BNRMTO = BIGNUM
  442. ILBSCL = .TRUE.
  443. END IF
  444. IF( ILBSCL )
  445. $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  446. *
  447. * Permute and/or balance the matrix pair (A,B)
  448. * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
  449. *
  450. CALL DGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
  451. $ WORK, IERR )
  452. *
  453. * Compute ABNRM and BBNRM
  454. *
  455. ABNRM = DLANGE( '1', N, N, A, LDA, WORK( 1 ) )
  456. IF( ILASCL ) THEN
  457. WORK( 1 ) = ABNRM
  458. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, WORK( 1 ), 1,
  459. $ IERR )
  460. ABNRM = WORK( 1 )
  461. END IF
  462. *
  463. BBNRM = DLANGE( '1', N, N, B, LDB, WORK( 1 ) )
  464. IF( ILBSCL ) THEN
  465. WORK( 1 ) = BBNRM
  466. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, WORK( 1 ), 1,
  467. $ IERR )
  468. BBNRM = WORK( 1 )
  469. END IF
  470. *
  471. * Reduce B to triangular form (QR decomposition of B)
  472. * (Workspace: need N, prefer N*NB )
  473. *
  474. IROWS = IHI + 1 - ILO
  475. IF( ILV .OR. .NOT.WANTSN ) THEN
  476. ICOLS = N + 1 - ILO
  477. ELSE
  478. ICOLS = IROWS
  479. END IF
  480. ITAU = 1
  481. IWRK = ITAU + IROWS
  482. CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  483. $ WORK( IWRK ), LWORK+1-IWRK, IERR )
  484. *
  485. * Apply the orthogonal transformation to A
  486. * (Workspace: need N, prefer N*NB)
  487. *
  488. CALL DORMQR( '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 VL and/or VR
  493. * (Workspace: need N, prefer N*NB)
  494. *
  495. IF( ILVL ) THEN
  496. CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
  497. IF( IROWS.GT.1 ) THEN
  498. CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  499. $ VL( ILO+1, ILO ), LDVL )
  500. END IF
  501. CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
  502. $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  503. END IF
  504. *
  505. IF( ILVR )
  506. $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
  507. *
  508. * Reduce to generalized Hessenberg form
  509. * (Workspace: none needed)
  510. *
  511. IF( ILV .OR. .NOT.WANTSN ) THEN
  512. *
  513. * Eigenvectors requested -- work on whole matrix.
  514. *
  515. CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
  516. $ LDVL, VR, LDVR, IERR )
  517. ELSE
  518. CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
  519. $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
  520. END IF
  521. *
  522. * Perform QZ algorithm (Compute eigenvalues, and optionally, the
  523. * Schur forms and Schur vectors)
  524. * (Workspace: need N)
  525. *
  526. IF( ILV .OR. .NOT.WANTSN ) THEN
  527. CHTEMP = 'S'
  528. ELSE
  529. CHTEMP = 'E'
  530. END IF
  531. *
  532. CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
  533. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK,
  534. $ LWORK, IERR )
  535. IF( IERR.NE.0 ) THEN
  536. IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
  537. INFO = IERR
  538. ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
  539. INFO = IERR - N
  540. ELSE
  541. INFO = N + 1
  542. END IF
  543. GO TO 130
  544. END IF
  545. *
  546. * Compute Eigenvectors and estimate condition numbers if desired
  547. * (Workspace: DTGEVC: need 6*N
  548. * DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
  549. * need N otherwise )
  550. *
  551. IF( ILV .OR. .NOT.WANTSN ) THEN
  552. IF( ILV ) THEN
  553. IF( ILVL ) THEN
  554. IF( ILVR ) THEN
  555. CHTEMP = 'B'
  556. ELSE
  557. CHTEMP = 'L'
  558. END IF
  559. ELSE
  560. CHTEMP = 'R'
  561. END IF
  562. *
  563. CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
  564. $ LDVL, VR, LDVR, N, IN, WORK, IERR )
  565. IF( IERR.NE.0 ) THEN
  566. INFO = N + 2
  567. GO TO 130
  568. END IF
  569. END IF
  570. *
  571. IF( .NOT.WANTSN ) THEN
  572. *
  573. * compute eigenvectors (DTGEVC) and estimate condition
  574. * numbers (DTGSNA). Note that the definition of the condition
  575. * number is not invariant under transformation (u,v) to
  576. * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
  577. * Schur form (S,T), Q and Z are orthogonal matrices. In order
  578. * to avoid using extra 2*N*N workspace, we have to recalculate
  579. * eigenvectors and estimate one condition numbers at a time.
  580. *
  581. PAIR = .FALSE.
  582. DO 20 I = 1, N
  583. *
  584. IF( PAIR ) THEN
  585. PAIR = .FALSE.
  586. GO TO 20
  587. END IF
  588. MM = 1
  589. IF( I.LT.N ) THEN
  590. IF( A( I+1, I ).NE.ZERO ) THEN
  591. PAIR = .TRUE.
  592. MM = 2
  593. END IF
  594. END IF
  595. *
  596. DO 10 J = 1, N
  597. BWORK( J ) = .FALSE.
  598. 10 CONTINUE
  599. IF( MM.EQ.1 ) THEN
  600. BWORK( I ) = .TRUE.
  601. ELSE IF( MM.EQ.2 ) THEN
  602. BWORK( I ) = .TRUE.
  603. BWORK( I+1 ) = .TRUE.
  604. END IF
  605. *
  606. IWRK = MM*N + 1
  607. IWRK1 = IWRK + MM*N
  608. *
  609. * Compute a pair of left and right eigenvectors.
  610. * (compute workspace: need up to 4*N + 6*N)
  611. *
  612. IF( WANTSE .OR. WANTSB ) THEN
  613. CALL DTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
  614. $ WORK( 1 ), N, WORK( IWRK ), N, MM, M,
  615. $ WORK( IWRK1 ), IERR )
  616. IF( IERR.NE.0 ) THEN
  617. INFO = N + 2
  618. GO TO 130
  619. END IF
  620. END IF
  621. *
  622. CALL DTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
  623. $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
  624. $ RCONDV( I ), MM, M, WORK( IWRK1 ),
  625. $ LWORK-IWRK1+1, IWORK, IERR )
  626. *
  627. 20 CONTINUE
  628. END IF
  629. END IF
  630. *
  631. * Undo balancing on VL and VR and normalization
  632. * (Workspace: none needed)
  633. *
  634. IF( ILVL ) THEN
  635. CALL DGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
  636. $ LDVL, IERR )
  637. *
  638. DO 70 JC = 1, N
  639. IF( ALPHAI( JC ).LT.ZERO )
  640. $ GO TO 70
  641. TEMP = ZERO
  642. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  643. DO 30 JR = 1, N
  644. TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
  645. 30 CONTINUE
  646. ELSE
  647. DO 40 JR = 1, N
  648. TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
  649. $ ABS( VL( JR, JC+1 ) ) )
  650. 40 CONTINUE
  651. END IF
  652. IF( TEMP.LT.SMLNUM )
  653. $ GO TO 70
  654. TEMP = ONE / TEMP
  655. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  656. DO 50 JR = 1, N
  657. VL( JR, JC ) = VL( JR, JC )*TEMP
  658. 50 CONTINUE
  659. ELSE
  660. DO 60 JR = 1, N
  661. VL( JR, JC ) = VL( JR, JC )*TEMP
  662. VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
  663. 60 CONTINUE
  664. END IF
  665. 70 CONTINUE
  666. END IF
  667. IF( ILVR ) THEN
  668. CALL DGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
  669. $ LDVR, IERR )
  670. DO 120 JC = 1, N
  671. IF( ALPHAI( JC ).LT.ZERO )
  672. $ GO TO 120
  673. TEMP = ZERO
  674. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  675. DO 80 JR = 1, N
  676. TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
  677. 80 CONTINUE
  678. ELSE
  679. DO 90 JR = 1, N
  680. TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
  681. $ ABS( VR( JR, JC+1 ) ) )
  682. 90 CONTINUE
  683. END IF
  684. IF( TEMP.LT.SMLNUM )
  685. $ GO TO 120
  686. TEMP = ONE / TEMP
  687. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  688. DO 100 JR = 1, N
  689. VR( JR, JC ) = VR( JR, JC )*TEMP
  690. 100 CONTINUE
  691. ELSE
  692. DO 110 JR = 1, N
  693. VR( JR, JC ) = VR( JR, JC )*TEMP
  694. VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
  695. 110 CONTINUE
  696. END IF
  697. 120 CONTINUE
  698. END IF
  699. *
  700. * Undo scaling if necessary
  701. *
  702. IF( ILASCL ) THEN
  703. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
  704. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
  705. END IF
  706. *
  707. IF( ILBSCL ) THEN
  708. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  709. END IF
  710. *
  711. 130 CONTINUE
  712. WORK( 1 ) = MAXWRK
  713. *
  714. RETURN
  715. *
  716. * End of DGGEVX
  717. *
  718. END