PageRenderTime 37ms CodeModel.GetById 6ms RepoModel.GetById 1ms app.codeStats 0ms

/SRC/dggevx.f

https://bitbucket.org/iricpt/lapack
FORTRAN Legacy | 865 lines | 320 code | 0 blank | 545 comment | 0 complexity | d3ac8227e2ab12bf6f5559183441d0e9 MD5 | raw file
Possible License(s): BSD-3-Clause
  1. *> \brief <b> DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download DGGEVX + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggevx.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggevx.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggevx.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
  22. * ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
  23. * IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
  24. * RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
  25. *
  26. * .. Scalar Arguments ..
  27. * CHARACTER BALANC, JOBVL, JOBVR, SENSE
  28. * INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  29. * DOUBLE PRECISION ABNRM, BBNRM
  30. * ..
  31. * .. Array Arguments ..
  32. * LOGICAL BWORK( * )
  33. * INTEGER IWORK( * )
  34. * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  35. * $ B( LDB, * ), BETA( * ), LSCALE( * ),
  36. * $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
  37. * $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
  38. * ..
  39. *
  40. *
  41. *> \par Purpose:
  42. * =============
  43. *>
  44. *> \verbatim
  45. *>
  46. *> DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
  47. *> the generalized eigenvalues, and optionally, the left and/or right
  48. *> generalized eigenvectors.
  49. *>
  50. *> Optionally also, it computes a balancing transformation to improve
  51. *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
  52. *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
  53. *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
  54. *> right eigenvectors (RCONDV).
  55. *>
  56. *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
  57. *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
  58. *> singular. It is usually represented as the pair (alpha,beta), as
  59. *> there is a reasonable interpretation for beta=0, and even for both
  60. *> being zero.
  61. *>
  62. *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
  63. *> of (A,B) satisfies
  64. *>
  65. *> A * v(j) = lambda(j) * B * v(j) .
  66. *>
  67. *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
  68. *> of (A,B) satisfies
  69. *>
  70. *> u(j)**H * A = lambda(j) * u(j)**H * B.
  71. *>
  72. *> where u(j)**H is the conjugate-transpose of u(j).
  73. *>
  74. *> \endverbatim
  75. *
  76. * Arguments:
  77. * ==========
  78. *
  79. *> \param[in] BALANC
  80. *> \verbatim
  81. *> BALANC is CHARACTER*1
  82. *> Specifies the balance option to be performed.
  83. *> = 'N': do not diagonally scale or permute;
  84. *> = 'P': permute only;
  85. *> = 'S': scale only;
  86. *> = 'B': both permute and scale.
  87. *> Computed reciprocal condition numbers will be for the
  88. *> matrices after permuting and/or balancing. Permuting does
  89. *> not change condition numbers (in exact arithmetic), but
  90. *> balancing does.
  91. *> \endverbatim
  92. *>
  93. *> \param[in] JOBVL
  94. *> \verbatim
  95. *> JOBVL is CHARACTER*1
  96. *> = 'N': do not compute the left generalized eigenvectors;
  97. *> = 'V': compute the left generalized eigenvectors.
  98. *> \endverbatim
  99. *>
  100. *> \param[in] JOBVR
  101. *> \verbatim
  102. *> JOBVR is CHARACTER*1
  103. *> = 'N': do not compute the right generalized eigenvectors;
  104. *> = 'V': compute the right generalized eigenvectors.
  105. *> \endverbatim
  106. *>
  107. *> \param[in] SENSE
  108. *> \verbatim
  109. *> SENSE is CHARACTER*1
  110. *> Determines which reciprocal condition numbers are computed.
  111. *> = 'N': none are computed;
  112. *> = 'E': computed for eigenvalues only;
  113. *> = 'V': computed for eigenvectors only;
  114. *> = 'B': computed for eigenvalues and eigenvectors.
  115. *> \endverbatim
  116. *>
  117. *> \param[in] N
  118. *> \verbatim
  119. *> N is INTEGER
  120. *> The order of the matrices A, B, VL, and VR. N >= 0.
  121. *> \endverbatim
  122. *>
  123. *> \param[in,out] A
  124. *> \verbatim
  125. *> A is DOUBLE PRECISION array, dimension (LDA, N)
  126. *> On entry, the matrix A in the pair (A,B).
  127. *> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
  128. *> or both, then A contains the first part of the real Schur
  129. *> form of the "balanced" versions of the input A and B.
  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 DOUBLE PRECISION array, dimension (LDB, N)
  141. *> On entry, the matrix B in the pair (A,B).
  142. *> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
  143. *> or both, then B contains the second part of the real Schur
  144. *> form of the "balanced" versions of the input A and B.
  145. *> \endverbatim
  146. *>
  147. *> \param[in] LDB
  148. *> \verbatim
  149. *> LDB is INTEGER
  150. *> The leading dimension of B. LDB >= max(1,N).
  151. *> \endverbatim
  152. *>
  153. *> \param[out] ALPHAR
  154. *> \verbatim
  155. *> ALPHAR is DOUBLE PRECISION array, dimension (N)
  156. *> \endverbatim
  157. *>
  158. *> \param[out] ALPHAI
  159. *> \verbatim
  160. *> ALPHAI is DOUBLE PRECISION array, dimension (N)
  161. *> \endverbatim
  162. *>
  163. *> \param[out] BETA
  164. *> \verbatim
  165. *> BETA is DOUBLE PRECISION array, dimension (N)
  166. *> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
  167. *> be the generalized eigenvalues. If ALPHAI(j) is zero, then
  168. *> the j-th eigenvalue is real; if positive, then the j-th and
  169. *> (j+1)-st eigenvalues are a complex conjugate pair, with
  170. *> ALPHAI(j+1) negative.
  171. *>
  172. *> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
  173. *> may easily over- or underflow, and BETA(j) may even be zero.
  174. *> Thus, the user should avoid naively computing the ratio
  175. *> ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
  176. *> than and usually comparable with norm(A) in magnitude, and
  177. *> BETA always less than and usually comparable with norm(B).
  178. *> \endverbatim
  179. *>
  180. *> \param[out] VL
  181. *> \verbatim
  182. *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
  183. *> If JOBVL = 'V', the left eigenvectors u(j) are stored one
  184. *> after another in the columns of VL, in the same order as
  185. *> their eigenvalues. If the j-th eigenvalue is real, then
  186. *> u(j) = VL(:,j), the j-th column of VL. If the j-th and
  187. *> (j+1)-th eigenvalues form a complex conjugate pair, then
  188. *> u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
  189. *> Each eigenvector will be scaled so the largest component have
  190. *> abs(real part) + abs(imag. part) = 1.
  191. *> Not referenced if JOBVL = 'N'.
  192. *> \endverbatim
  193. *>
  194. *> \param[in] LDVL
  195. *> \verbatim
  196. *> LDVL is INTEGER
  197. *> The leading dimension of the matrix VL. LDVL >= 1, and
  198. *> if JOBVL = 'V', LDVL >= N.
  199. *> \endverbatim
  200. *>
  201. *> \param[out] VR
  202. *> \verbatim
  203. *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
  204. *> If JOBVR = 'V', the right eigenvectors v(j) are stored one
  205. *> after another in the columns of VR, in the same order as
  206. *> their eigenvalues. If the j-th eigenvalue is real, then
  207. *> v(j) = VR(:,j), the j-th column of VR. If the j-th and
  208. *> (j+1)-th eigenvalues form a complex conjugate pair, then
  209. *> v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
  210. *> Each eigenvector will be scaled so the largest component have
  211. *> abs(real part) + abs(imag. part) = 1.
  212. *> Not referenced if JOBVR = 'N'.
  213. *> \endverbatim
  214. *>
  215. *> \param[in] LDVR
  216. *> \verbatim
  217. *> LDVR is INTEGER
  218. *> The leading dimension of the matrix VR. LDVR >= 1, and
  219. *> if JOBVR = 'V', LDVR >= N.
  220. *> \endverbatim
  221. *>
  222. *> \param[out] ILO
  223. *> \verbatim
  224. *> ILO is INTEGER
  225. *> \endverbatim
  226. *>
  227. *> \param[out] IHI
  228. *> \verbatim
  229. *> IHI is INTEGER
  230. *> ILO and IHI are integer values such that on exit
  231. *> A(i,j) = 0 and B(i,j) = 0 if i > j and
  232. *> j = 1,...,ILO-1 or i = IHI+1,...,N.
  233. *> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
  234. *> \endverbatim
  235. *>
  236. *> \param[out] LSCALE
  237. *> \verbatim
  238. *> LSCALE is DOUBLE PRECISION array, dimension (N)
  239. *> Details of the permutations and scaling factors applied
  240. *> to the left side of A and B. If PL(j) is the index of the
  241. *> row interchanged with row j, and DL(j) is the scaling
  242. *> factor applied to row j, then
  243. *> LSCALE(j) = PL(j) for j = 1,...,ILO-1
  244. *> = DL(j) for j = ILO,...,IHI
  245. *> = PL(j) for j = IHI+1,...,N.
  246. *> The order in which the interchanges are made is N to IHI+1,
  247. *> then 1 to ILO-1.
  248. *> \endverbatim
  249. *>
  250. *> \param[out] RSCALE
  251. *> \verbatim
  252. *> RSCALE is DOUBLE PRECISION array, dimension (N)
  253. *> Details of the permutations and scaling factors applied
  254. *> to the right side of A and B. If PR(j) is the index of the
  255. *> column interchanged with column j, and DR(j) is the scaling
  256. *> factor applied to column j, then
  257. *> RSCALE(j) = PR(j) for j = 1,...,ILO-1
  258. *> = DR(j) for j = ILO,...,IHI
  259. *> = PR(j) for j = IHI+1,...,N
  260. *> The order in which the interchanges are made is N to IHI+1,
  261. *> then 1 to ILO-1.
  262. *> \endverbatim
  263. *>
  264. *> \param[out] ABNRM
  265. *> \verbatim
  266. *> ABNRM is DOUBLE PRECISION
  267. *> The one-norm of the balanced matrix A.
  268. *> \endverbatim
  269. *>
  270. *> \param[out] BBNRM
  271. *> \verbatim
  272. *> BBNRM is DOUBLE PRECISION
  273. *> The one-norm of the balanced matrix B.
  274. *> \endverbatim
  275. *>
  276. *> \param[out] RCONDE
  277. *> \verbatim
  278. *> RCONDE is DOUBLE PRECISION array, dimension (N)
  279. *> If SENSE = 'E' or 'B', the reciprocal condition numbers of
  280. *> the eigenvalues, stored in consecutive elements of the array.
  281. *> For a complex conjugate pair of eigenvalues two consecutive
  282. *> elements of RCONDE are set to the same value. Thus RCONDE(j),
  283. *> RCONDV(j), and the j-th columns of VL and VR all correspond
  284. *> to the j-th eigenpair.
  285. *> If SENSE = 'N or 'V', RCONDE is not referenced.
  286. *> \endverbatim
  287. *>
  288. *> \param[out] RCONDV
  289. *> \verbatim
  290. *> RCONDV is DOUBLE PRECISION array, dimension (N)
  291. *> If SENSE = 'V' or 'B', the estimated reciprocal condition
  292. *> numbers of the eigenvectors, stored in consecutive elements
  293. *> of the array. For a complex eigenvector two consecutive
  294. *> elements of RCONDV are set to the same value. If the
  295. *> eigenvalues cannot be reordered to compute RCONDV(j),
  296. *> RCONDV(j) is set to 0; this can only occur when the true
  297. *> value would be very small anyway.
  298. *> If SENSE = 'N' or 'E', RCONDV is not referenced.
  299. *> \endverbatim
  300. *>
  301. *> \param[out] WORK
  302. *> \verbatim
  303. *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  304. *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  305. *> \endverbatim
  306. *>
  307. *> \param[in] LWORK
  308. *> \verbatim
  309. *> LWORK is INTEGER
  310. *> The dimension of the array WORK. LWORK >= max(1,2*N).
  311. *> If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
  312. *> LWORK >= max(1,6*N).
  313. *> If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
  314. *> If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
  315. *>
  316. *> If LWORK = -1, then a workspace query is assumed; the routine
  317. *> only calculates the optimal size of the WORK array, returns
  318. *> this value as the first entry of the WORK array, and no error
  319. *> message related to LWORK is issued by XERBLA.
  320. *> \endverbatim
  321. *>
  322. *> \param[out] IWORK
  323. *> \verbatim
  324. *> IWORK is INTEGER array, dimension (N+6)
  325. *> If SENSE = 'E', IWORK is not referenced.
  326. *> \endverbatim
  327. *>
  328. *> \param[out] BWORK
  329. *> \verbatim
  330. *> BWORK is LOGICAL array, dimension (N)
  331. *> If SENSE = 'N', BWORK is not referenced.
  332. *> \endverbatim
  333. *>
  334. *> \param[out] INFO
  335. *> \verbatim
  336. *> INFO is INTEGER
  337. *> = 0: successful exit
  338. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  339. *> = 1,...,N:
  340. *> The QZ iteration failed. No eigenvectors have been
  341. *> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
  342. *> should be correct for j=INFO+1,...,N.
  343. *> > N: =N+1: other than QZ iteration failed in DHGEQZ.
  344. *> =N+2: error return from DTGEVC.
  345. *> \endverbatim
  346. *
  347. * Authors:
  348. * ========
  349. *
  350. *> \author Univ. of Tennessee
  351. *> \author Univ. of California Berkeley
  352. *> \author Univ. of Colorado Denver
  353. *> \author NAG Ltd.
  354. *
  355. *> \ingroup doubleGEeigen
  356. *
  357. *> \par Further Details:
  358. * =====================
  359. *>
  360. *> \verbatim
  361. *>
  362. *> Balancing a matrix pair (A,B) includes, first, permuting rows and
  363. *> columns to isolate eigenvalues, second, applying diagonal similarity
  364. *> transformation to the rows and columns to make the rows and columns
  365. *> as close in norm as possible. The computed reciprocal condition
  366. *> numbers correspond to the balanced matrix. Permuting rows and columns
  367. *> will not change the condition numbers (in exact arithmetic) but
  368. *> diagonal scaling will. For further explanation of balancing, see
  369. *> section 4.11.1.2 of LAPACK Users' Guide.
  370. *>
  371. *> An approximate error bound on the chordal distance between the i-th
  372. *> computed generalized eigenvalue w and the corresponding exact
  373. *> eigenvalue lambda is
  374. *>
  375. *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
  376. *>
  377. *> An approximate error bound for the angle between the i-th computed
  378. *> eigenvector VL(i) or VR(i) is given by
  379. *>
  380. *> EPS * norm(ABNRM, BBNRM) / DIF(i).
  381. *>
  382. *> For further explanation of the reciprocal condition numbers RCONDE
  383. *> and RCONDV, see section 4.11 of LAPACK User's Guide.
  384. *> \endverbatim
  385. *>
  386. * =====================================================================
  387. SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
  388. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
  389. $ IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
  390. $ RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
  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. DOUBLE PRECISION ABNRM, BBNRM
  400. * ..
  401. * .. Array Arguments ..
  402. LOGICAL BWORK( * )
  403. INTEGER IWORK( * )
  404. DOUBLE PRECISION 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. DOUBLE PRECISION ZERO, ONE
  414. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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. DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
  424. $ SMLNUM, TEMP
  425. * ..
  426. * .. Local Arrays ..
  427. LOGICAL LDUMMA( 1 )
  428. * ..
  429. * .. External Subroutines ..
  430. EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
  431. $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
  432. $ DTGSNA, XERBLA
  433. * ..
  434. * .. External Functions ..
  435. LOGICAL LSAME
  436. INTEGER ILAENV
  437. DOUBLE PRECISION DLAMCH, DLANGE
  438. EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
  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.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
  481. $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
  482. $ THEN
  483. INFO = -1
  484. ELSE IF( IJOBVL.LE.0 ) THEN
  485. INFO = -2
  486. ELSE IF( IJOBVR.LE.0 ) THEN
  487. INFO = -3
  488. ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
  489. $ THEN
  490. INFO = -4
  491. ELSE IF( N.LT.0 ) THEN
  492. INFO = -5
  493. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  494. INFO = -7
  495. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  496. INFO = -9
  497. ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
  498. INFO = -14
  499. ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
  500. INFO = -16
  501. END IF
  502. *
  503. * Compute workspace
  504. * (Note: Comments in the code beginning "Workspace:" describe the
  505. * minimal amount of workspace needed at that point in the code,
  506. * as well as the preferred amount for good performance.
  507. * NB refers to the optimal block size for the immediately
  508. * following subroutine, as returned by ILAENV. The workspace is
  509. * computed assuming ILO = 1 and IHI = N, the worst case.)
  510. *
  511. IF( INFO.EQ.0 ) THEN
  512. IF( N.EQ.0 ) THEN
  513. MINWRK = 1
  514. MAXWRK = 1
  515. ELSE
  516. IF( NOSCL .AND. .NOT.ILV ) THEN
  517. MINWRK = 2*N
  518. ELSE
  519. MINWRK = 6*N
  520. END IF
  521. IF( WANTSE .OR. WANTSB ) THEN
  522. MINWRK = 10*N
  523. END IF
  524. IF( WANTSV .OR. WANTSB ) THEN
  525. MINWRK = MAX( MINWRK, 2*N*( N + 4 ) + 16 )
  526. END IF
  527. MAXWRK = MINWRK
  528. MAXWRK = MAX( MAXWRK,
  529. $ N + N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) )
  530. MAXWRK = MAX( MAXWRK,
  531. $ N + N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) )
  532. IF( ILVL ) THEN
  533. MAXWRK = MAX( MAXWRK, N +
  534. $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, 0 ) )
  535. END IF
  536. END IF
  537. WORK( 1 ) = MAXWRK
  538. *
  539. IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
  540. INFO = -26
  541. END IF
  542. END IF
  543. *
  544. IF( INFO.NE.0 ) THEN
  545. CALL XERBLA( 'DGGEVX', -INFO )
  546. RETURN
  547. ELSE IF( LQUERY ) THEN
  548. RETURN
  549. END IF
  550. *
  551. * Quick return if possible
  552. *
  553. IF( N.EQ.0 )
  554. $ RETURN
  555. *
  556. *
  557. * Get machine constants
  558. *
  559. EPS = DLAMCH( 'P' )
  560. SMLNUM = DLAMCH( 'S' )
  561. BIGNUM = ONE / SMLNUM
  562. CALL DLABAD( SMLNUM, BIGNUM )
  563. SMLNUM = SQRT( SMLNUM ) / EPS
  564. BIGNUM = ONE / SMLNUM
  565. *
  566. * Scale A if max element outside range [SMLNUM,BIGNUM]
  567. *
  568. ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
  569. ILASCL = .FALSE.
  570. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  571. ANRMTO = SMLNUM
  572. ILASCL = .TRUE.
  573. ELSE IF( ANRM.GT.BIGNUM ) THEN
  574. ANRMTO = BIGNUM
  575. ILASCL = .TRUE.
  576. END IF
  577. IF( ILASCL )
  578. $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  579. *
  580. * Scale B if max element outside range [SMLNUM,BIGNUM]
  581. *
  582. BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
  583. ILBSCL = .FALSE.
  584. IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  585. BNRMTO = SMLNUM
  586. ILBSCL = .TRUE.
  587. ELSE IF( BNRM.GT.BIGNUM ) THEN
  588. BNRMTO = BIGNUM
  589. ILBSCL = .TRUE.
  590. END IF
  591. IF( ILBSCL )
  592. $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  593. *
  594. * Permute and/or balance the matrix pair (A,B)
  595. * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
  596. *
  597. CALL DGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
  598. $ WORK, IERR )
  599. *
  600. * Compute ABNRM and BBNRM
  601. *
  602. ABNRM = DLANGE( '1', N, N, A, LDA, WORK( 1 ) )
  603. IF( ILASCL ) THEN
  604. WORK( 1 ) = ABNRM
  605. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, WORK( 1 ), 1,
  606. $ IERR )
  607. ABNRM = WORK( 1 )
  608. END IF
  609. *
  610. BBNRM = DLANGE( '1', N, N, B, LDB, WORK( 1 ) )
  611. IF( ILBSCL ) THEN
  612. WORK( 1 ) = BBNRM
  613. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, WORK( 1 ), 1,
  614. $ IERR )
  615. BBNRM = WORK( 1 )
  616. END IF
  617. *
  618. * Reduce B to triangular form (QR decomposition of B)
  619. * (Workspace: need N, prefer N*NB )
  620. *
  621. IROWS = IHI + 1 - ILO
  622. IF( ILV .OR. .NOT.WANTSN ) THEN
  623. ICOLS = N + 1 - ILO
  624. ELSE
  625. ICOLS = IROWS
  626. END IF
  627. ITAU = 1
  628. IWRK = ITAU + IROWS
  629. CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  630. $ WORK( IWRK ), LWORK+1-IWRK, IERR )
  631. *
  632. * Apply the orthogonal transformation to A
  633. * (Workspace: need N, prefer N*NB)
  634. *
  635. CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
  636. $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
  637. $ LWORK+1-IWRK, IERR )
  638. *
  639. * Initialize VL and/or VR
  640. * (Workspace: need N, prefer N*NB)
  641. *
  642. IF( ILVL ) THEN
  643. CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
  644. IF( IROWS.GT.1 ) THEN
  645. CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  646. $ VL( ILO+1, ILO ), LDVL )
  647. END IF
  648. CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
  649. $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  650. END IF
  651. *
  652. IF( ILVR )
  653. $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
  654. *
  655. * Reduce to generalized Hessenberg form
  656. * (Workspace: none needed)
  657. *
  658. IF( ILV .OR. .NOT.WANTSN ) THEN
  659. *
  660. * Eigenvectors requested -- work on whole matrix.
  661. *
  662. CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
  663. $ LDVL, VR, LDVR, IERR )
  664. ELSE
  665. CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
  666. $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
  667. END IF
  668. *
  669. * Perform QZ algorithm (Compute eigenvalues, and optionally, the
  670. * Schur forms and Schur vectors)
  671. * (Workspace: need N)
  672. *
  673. IF( ILV .OR. .NOT.WANTSN ) THEN
  674. CHTEMP = 'S'
  675. ELSE
  676. CHTEMP = 'E'
  677. END IF
  678. *
  679. CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
  680. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK,
  681. $ LWORK, IERR )
  682. IF( IERR.NE.0 ) THEN
  683. IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
  684. INFO = IERR
  685. ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
  686. INFO = IERR - N
  687. ELSE
  688. INFO = N + 1
  689. END IF
  690. GO TO 130
  691. END IF
  692. *
  693. * Compute Eigenvectors and estimate condition numbers if desired
  694. * (Workspace: DTGEVC: need 6*N
  695. * DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
  696. * need N otherwise )
  697. *
  698. IF( ILV .OR. .NOT.WANTSN ) THEN
  699. IF( ILV ) THEN
  700. IF( ILVL ) THEN
  701. IF( ILVR ) THEN
  702. CHTEMP = 'B'
  703. ELSE
  704. CHTEMP = 'L'
  705. END IF
  706. ELSE
  707. CHTEMP = 'R'
  708. END IF
  709. *
  710. CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
  711. $ LDVL, VR, LDVR, N, IN, WORK, IERR )
  712. IF( IERR.NE.0 ) THEN
  713. INFO = N + 2
  714. GO TO 130
  715. END IF
  716. END IF
  717. *
  718. IF( .NOT.WANTSN ) THEN
  719. *
  720. * compute eigenvectors (DTGEVC) and estimate condition
  721. * numbers (DTGSNA). Note that the definition of the condition
  722. * number is not invariant under transformation (u,v) to
  723. * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
  724. * Schur form (S,T), Q and Z are orthogonal matrices. In order
  725. * to avoid using extra 2*N*N workspace, we have to recalculate
  726. * eigenvectors and estimate one condition numbers at a time.
  727. *
  728. PAIR = .FALSE.
  729. DO 20 I = 1, N
  730. *
  731. IF( PAIR ) THEN
  732. PAIR = .FALSE.
  733. GO TO 20
  734. END IF
  735. MM = 1
  736. IF( I.LT.N ) THEN
  737. IF( A( I+1, I ).NE.ZERO ) THEN
  738. PAIR = .TRUE.
  739. MM = 2
  740. END IF
  741. END IF
  742. *
  743. DO 10 J = 1, N
  744. BWORK( J ) = .FALSE.
  745. 10 CONTINUE
  746. IF( MM.EQ.1 ) THEN
  747. BWORK( I ) = .TRUE.
  748. ELSE IF( MM.EQ.2 ) THEN
  749. BWORK( I ) = .TRUE.
  750. BWORK( I+1 ) = .TRUE.
  751. END IF
  752. *
  753. IWRK = MM*N + 1
  754. IWRK1 = IWRK + MM*N
  755. *
  756. * Compute a pair of left and right eigenvectors.
  757. * (compute workspace: need up to 4*N + 6*N)
  758. *
  759. IF( WANTSE .OR. WANTSB ) THEN
  760. CALL DTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
  761. $ WORK( 1 ), N, WORK( IWRK ), N, MM, M,
  762. $ WORK( IWRK1 ), IERR )
  763. IF( IERR.NE.0 ) THEN
  764. INFO = N + 2
  765. GO TO 130
  766. END IF
  767. END IF
  768. *
  769. CALL DTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
  770. $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
  771. $ RCONDV( I ), MM, M, WORK( IWRK1 ),
  772. $ LWORK-IWRK1+1, IWORK, IERR )
  773. *
  774. 20 CONTINUE
  775. END IF
  776. END IF
  777. *
  778. * Undo balancing on VL and VR and normalization
  779. * (Workspace: none needed)
  780. *
  781. IF( ILVL ) THEN
  782. CALL DGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
  783. $ LDVL, IERR )
  784. *
  785. DO 70 JC = 1, N
  786. IF( ALPHAI( JC ).LT.ZERO )
  787. $ GO TO 70
  788. TEMP = ZERO
  789. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  790. DO 30 JR = 1, N
  791. TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
  792. 30 CONTINUE
  793. ELSE
  794. DO 40 JR = 1, N
  795. TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
  796. $ ABS( VL( JR, JC+1 ) ) )
  797. 40 CONTINUE
  798. END IF
  799. IF( TEMP.LT.SMLNUM )
  800. $ GO TO 70
  801. TEMP = ONE / TEMP
  802. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  803. DO 50 JR = 1, N
  804. VL( JR, JC ) = VL( JR, JC )*TEMP
  805. 50 CONTINUE
  806. ELSE
  807. DO 60 JR = 1, N
  808. VL( JR, JC ) = VL( JR, JC )*TEMP
  809. VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
  810. 60 CONTINUE
  811. END IF
  812. 70 CONTINUE
  813. END IF
  814. IF( ILVR ) THEN
  815. CALL DGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
  816. $ LDVR, IERR )
  817. DO 120 JC = 1, N
  818. IF( ALPHAI( JC ).LT.ZERO )
  819. $ GO TO 120
  820. TEMP = ZERO
  821. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  822. DO 80 JR = 1, N
  823. TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
  824. 80 CONTINUE
  825. ELSE
  826. DO 90 JR = 1, N
  827. TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
  828. $ ABS( VR( JR, JC+1 ) ) )
  829. 90 CONTINUE
  830. END IF
  831. IF( TEMP.LT.SMLNUM )
  832. $ GO TO 120
  833. TEMP = ONE / TEMP
  834. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  835. DO 100 JR = 1, N
  836. VR( JR, JC ) = VR( JR, JC )*TEMP
  837. 100 CONTINUE
  838. ELSE
  839. DO 110 JR = 1, N
  840. VR( JR, JC ) = VR( JR, JC )*TEMP
  841. VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
  842. 110 CONTINUE
  843. END IF
  844. 120 CONTINUE
  845. END IF
  846. *
  847. * Undo scaling if necessary
  848. *
  849. 130 CONTINUE
  850. *
  851. IF( ILASCL ) THEN
  852. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
  853. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
  854. END IF
  855. *
  856. IF( ILBSCL ) THEN
  857. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  858. END IF
  859. *
  860. WORK( 1 ) = MAXWRK
  861. RETURN
  862. *
  863. * End of DGGEVX
  864. *
  865. END