PageRenderTime 60ms CodeModel.GetById 26ms RepoModel.GetById 1ms app.codeStats 0ms

/src/LAPACK/dggevx.f

http://github.com/exciting/exciting
FORTRAN Legacy | 868 lines | 320 code | 0 blank | 548 comment | 0 complexity | 21ef4595454b1ccfd94de4d339afa7db MD5 | raw file
Possible License(s): LGPL-2.1, BSD-3-Clause, LGPL-3.0
  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. *> \date November 2011
  356. *
  357. *> \ingroup doubleGEeigen
  358. *
  359. *> \par Further Details:
  360. * =====================
  361. *>
  362. *> \verbatim
  363. *>
  364. *> Balancing a matrix pair (A,B) includes, first, permuting rows and
  365. *> columns to isolate eigenvalues, second, applying diagonal similarity
  366. *> transformation to the rows and columns to make the rows and columns
  367. *> as close in norm as possible. The computed reciprocal condition
  368. *> numbers correspond to the balanced matrix. Permuting rows and columns
  369. *> will not change the condition numbers (in exact arithmetic) but
  370. *> diagonal scaling will. For further explanation of balancing, see
  371. *> section 4.11.1.2 of LAPACK Users' Guide.
  372. *>
  373. *> An approximate error bound on the chordal distance between the i-th
  374. *> computed generalized eigenvalue w and the corresponding exact
  375. *> eigenvalue lambda is
  376. *>
  377. *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
  378. *>
  379. *> An approximate error bound for the angle between the i-th computed
  380. *> eigenvector VL(i) or VR(i) is given by
  381. *>
  382. *> EPS * norm(ABNRM, BBNRM) / DIF(i).
  383. *>
  384. *> For further explanation of the reciprocal condition numbers RCONDE
  385. *> and RCONDV, see section 4.11 of LAPACK User's Guide.
  386. *> \endverbatim
  387. *>
  388. * =====================================================================
  389. SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
  390. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
  391. $ IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
  392. $ RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
  393. *
  394. * -- LAPACK driver routine (version 3.4.0) --
  395. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  396. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  397. * November 2011
  398. *
  399. * .. Scalar Arguments ..
  400. CHARACTER BALANC, JOBVL, JOBVR, SENSE
  401. INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  402. DOUBLE PRECISION ABNRM, BBNRM
  403. * ..
  404. * .. Array Arguments ..
  405. LOGICAL BWORK( * )
  406. INTEGER IWORK( * )
  407. DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  408. $ B( LDB, * ), BETA( * ), LSCALE( * ),
  409. $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
  410. $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
  411. * ..
  412. *
  413. * =====================================================================
  414. *
  415. * .. Parameters ..
  416. DOUBLE PRECISION ZERO, ONE
  417. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  418. * ..
  419. * .. Local Scalars ..
  420. LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
  421. $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
  422. CHARACTER CHTEMP
  423. INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
  424. $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
  425. $ MINWRK, MM
  426. DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
  427. $ SMLNUM, TEMP
  428. * ..
  429. * .. Local Arrays ..
  430. LOGICAL LDUMMA( 1 )
  431. * ..
  432. * .. External Subroutines ..
  433. EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
  434. $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
  435. $ DTGSNA, XERBLA
  436. * ..
  437. * .. External Functions ..
  438. LOGICAL LSAME
  439. INTEGER ILAENV
  440. DOUBLE PRECISION DLAMCH, DLANGE
  441. EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
  442. * ..
  443. * .. Intrinsic Functions ..
  444. INTRINSIC ABS, MAX, SQRT
  445. * ..
  446. * .. Executable Statements ..
  447. *
  448. * Decode the input arguments
  449. *
  450. IF( LSAME( JOBVL, 'N' ) ) THEN
  451. IJOBVL = 1
  452. ILVL = .FALSE.
  453. ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
  454. IJOBVL = 2
  455. ILVL = .TRUE.
  456. ELSE
  457. IJOBVL = -1
  458. ILVL = .FALSE.
  459. END IF
  460. *
  461. IF( LSAME( JOBVR, 'N' ) ) THEN
  462. IJOBVR = 1
  463. ILVR = .FALSE.
  464. ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
  465. IJOBVR = 2
  466. ILVR = .TRUE.
  467. ELSE
  468. IJOBVR = -1
  469. ILVR = .FALSE.
  470. END IF
  471. ILV = ILVL .OR. ILVR
  472. *
  473. NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
  474. WANTSN = LSAME( SENSE, 'N' )
  475. WANTSE = LSAME( SENSE, 'E' )
  476. WANTSV = LSAME( SENSE, 'V' )
  477. WANTSB = LSAME( SENSE, 'B' )
  478. *
  479. * Test the input arguments
  480. *
  481. INFO = 0
  482. LQUERY = ( LWORK.EQ.-1 )
  483. IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
  484. $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
  485. $ THEN
  486. INFO = -1
  487. ELSE IF( IJOBVL.LE.0 ) THEN
  488. INFO = -2
  489. ELSE IF( IJOBVR.LE.0 ) THEN
  490. INFO = -3
  491. ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
  492. $ THEN
  493. INFO = -4
  494. ELSE IF( N.LT.0 ) THEN
  495. INFO = -5
  496. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  497. INFO = -7
  498. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  499. INFO = -9
  500. ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
  501. INFO = -14
  502. ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
  503. INFO = -16
  504. END IF
  505. *
  506. * Compute workspace
  507. * (Note: Comments in the code beginning "Workspace:" describe the
  508. * minimal amount of workspace needed at that point in the code,
  509. * as well as the preferred amount for good performance.
  510. * NB refers to the optimal block size for the immediately
  511. * following subroutine, as returned by ILAENV. The workspace is
  512. * computed assuming ILO = 1 and IHI = N, the worst case.)
  513. *
  514. IF( INFO.EQ.0 ) THEN
  515. IF( N.EQ.0 ) THEN
  516. MINWRK = 1
  517. MAXWRK = 1
  518. ELSE
  519. IF( NOSCL .AND. .NOT.ILV ) THEN
  520. MINWRK = 2*N
  521. ELSE
  522. MINWRK = 6*N
  523. END IF
  524. IF( WANTSE .OR. WANTSB ) THEN
  525. MINWRK = 10*N
  526. END IF
  527. IF( WANTSV .OR. WANTSB ) THEN
  528. MINWRK = MAX( MINWRK, 2*N*( N + 4 ) + 16 )
  529. END IF
  530. MAXWRK = MINWRK
  531. MAXWRK = MAX( MAXWRK,
  532. $ N + N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) )
  533. MAXWRK = MAX( MAXWRK,
  534. $ N + N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) )
  535. IF( ILVL ) THEN
  536. MAXWRK = MAX( MAXWRK, N +
  537. $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, 0 ) )
  538. END IF
  539. END IF
  540. WORK( 1 ) = MAXWRK
  541. *
  542. IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
  543. INFO = -26
  544. END IF
  545. END IF
  546. *
  547. IF( INFO.NE.0 ) THEN
  548. CALL XERBLA( 'DGGEVX', -INFO )
  549. RETURN
  550. ELSE IF( LQUERY ) THEN
  551. RETURN
  552. END IF
  553. *
  554. * Quick return if possible
  555. *
  556. IF( N.EQ.0 )
  557. $ RETURN
  558. *
  559. *
  560. * Get machine constants
  561. *
  562. EPS = DLAMCH( 'P' )
  563. SMLNUM = DLAMCH( 'S' )
  564. BIGNUM = ONE / SMLNUM
  565. CALL DLABAD( SMLNUM, BIGNUM )
  566. SMLNUM = SQRT( SMLNUM ) / EPS
  567. BIGNUM = ONE / SMLNUM
  568. *
  569. * Scale A if max element outside range [SMLNUM,BIGNUM]
  570. *
  571. ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
  572. ILASCL = .FALSE.
  573. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  574. ANRMTO = SMLNUM
  575. ILASCL = .TRUE.
  576. ELSE IF( ANRM.GT.BIGNUM ) THEN
  577. ANRMTO = BIGNUM
  578. ILASCL = .TRUE.
  579. END IF
  580. IF( ILASCL )
  581. $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  582. *
  583. * Scale B if max element outside range [SMLNUM,BIGNUM]
  584. *
  585. BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
  586. ILBSCL = .FALSE.
  587. IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  588. BNRMTO = SMLNUM
  589. ILBSCL = .TRUE.
  590. ELSE IF( BNRM.GT.BIGNUM ) THEN
  591. BNRMTO = BIGNUM
  592. ILBSCL = .TRUE.
  593. END IF
  594. IF( ILBSCL )
  595. $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  596. *
  597. * Permute and/or balance the matrix pair (A,B)
  598. * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
  599. *
  600. CALL DGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
  601. $ WORK, IERR )
  602. *
  603. * Compute ABNRM and BBNRM
  604. *
  605. ABNRM = DLANGE( '1', N, N, A, LDA, WORK( 1 ) )
  606. IF( ILASCL ) THEN
  607. WORK( 1 ) = ABNRM
  608. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, WORK( 1 ), 1,
  609. $ IERR )
  610. ABNRM = WORK( 1 )
  611. END IF
  612. *
  613. BBNRM = DLANGE( '1', N, N, B, LDB, WORK( 1 ) )
  614. IF( ILBSCL ) THEN
  615. WORK( 1 ) = BBNRM
  616. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, WORK( 1 ), 1,
  617. $ IERR )
  618. BBNRM = WORK( 1 )
  619. END IF
  620. *
  621. * Reduce B to triangular form (QR decomposition of B)
  622. * (Workspace: need N, prefer N*NB )
  623. *
  624. IROWS = IHI + 1 - ILO
  625. IF( ILV .OR. .NOT.WANTSN ) THEN
  626. ICOLS = N + 1 - ILO
  627. ELSE
  628. ICOLS = IROWS
  629. END IF
  630. ITAU = 1
  631. IWRK = ITAU + IROWS
  632. CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  633. $ WORK( IWRK ), LWORK+1-IWRK, IERR )
  634. *
  635. * Apply the orthogonal transformation to A
  636. * (Workspace: need N, prefer N*NB)
  637. *
  638. CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
  639. $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
  640. $ LWORK+1-IWRK, IERR )
  641. *
  642. * Initialize VL and/or VR
  643. * (Workspace: need N, prefer N*NB)
  644. *
  645. IF( ILVL ) THEN
  646. CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
  647. IF( IROWS.GT.1 ) THEN
  648. CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  649. $ VL( ILO+1, ILO ), LDVL )
  650. END IF
  651. CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
  652. $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  653. END IF
  654. *
  655. IF( ILVR )
  656. $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
  657. *
  658. * Reduce to generalized Hessenberg form
  659. * (Workspace: none needed)
  660. *
  661. IF( ILV .OR. .NOT.WANTSN ) THEN
  662. *
  663. * Eigenvectors requested -- work on whole matrix.
  664. *
  665. CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
  666. $ LDVL, VR, LDVR, IERR )
  667. ELSE
  668. CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
  669. $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
  670. END IF
  671. *
  672. * Perform QZ algorithm (Compute eigenvalues, and optionally, the
  673. * Schur forms and Schur vectors)
  674. * (Workspace: need N)
  675. *
  676. IF( ILV .OR. .NOT.WANTSN ) THEN
  677. CHTEMP = 'S'
  678. ELSE
  679. CHTEMP = 'E'
  680. END IF
  681. *
  682. CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
  683. $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK,
  684. $ LWORK, IERR )
  685. IF( IERR.NE.0 ) THEN
  686. IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
  687. INFO = IERR
  688. ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
  689. INFO = IERR - N
  690. ELSE
  691. INFO = N + 1
  692. END IF
  693. GO TO 130
  694. END IF
  695. *
  696. * Compute Eigenvectors and estimate condition numbers if desired
  697. * (Workspace: DTGEVC: need 6*N
  698. * DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
  699. * need N otherwise )
  700. *
  701. IF( ILV .OR. .NOT.WANTSN ) THEN
  702. IF( ILV ) THEN
  703. IF( ILVL ) THEN
  704. IF( ILVR ) THEN
  705. CHTEMP = 'B'
  706. ELSE
  707. CHTEMP = 'L'
  708. END IF
  709. ELSE
  710. CHTEMP = 'R'
  711. END IF
  712. *
  713. CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
  714. $ LDVL, VR, LDVR, N, IN, WORK, IERR )
  715. IF( IERR.NE.0 ) THEN
  716. INFO = N + 2
  717. GO TO 130
  718. END IF
  719. END IF
  720. *
  721. IF( .NOT.WANTSN ) THEN
  722. *
  723. * compute eigenvectors (DTGEVC) and estimate condition
  724. * numbers (DTGSNA). Note that the definition of the condition
  725. * number is not invariant under transformation (u,v) to
  726. * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
  727. * Schur form (S,T), Q and Z are orthogonal matrices. In order
  728. * to avoid using extra 2*N*N workspace, we have to recalculate
  729. * eigenvectors and estimate one condition numbers at a time.
  730. *
  731. PAIR = .FALSE.
  732. DO 20 I = 1, N
  733. *
  734. IF( PAIR ) THEN
  735. PAIR = .FALSE.
  736. GO TO 20
  737. END IF
  738. MM = 1
  739. IF( I.LT.N ) THEN
  740. IF( A( I+1, I ).NE.ZERO ) THEN
  741. PAIR = .TRUE.
  742. MM = 2
  743. END IF
  744. END IF
  745. *
  746. DO 10 J = 1, N
  747. BWORK( J ) = .FALSE.
  748. 10 CONTINUE
  749. IF( MM.EQ.1 ) THEN
  750. BWORK( I ) = .TRUE.
  751. ELSE IF( MM.EQ.2 ) THEN
  752. BWORK( I ) = .TRUE.
  753. BWORK( I+1 ) = .TRUE.
  754. END IF
  755. *
  756. IWRK = MM*N + 1
  757. IWRK1 = IWRK + MM*N
  758. *
  759. * Compute a pair of left and right eigenvectors.
  760. * (compute workspace: need up to 4*N + 6*N)
  761. *
  762. IF( WANTSE .OR. WANTSB ) THEN
  763. CALL DTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
  764. $ WORK( 1 ), N, WORK( IWRK ), N, MM, M,
  765. $ WORK( IWRK1 ), IERR )
  766. IF( IERR.NE.0 ) THEN
  767. INFO = N + 2
  768. GO TO 130
  769. END IF
  770. END IF
  771. *
  772. CALL DTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
  773. $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
  774. $ RCONDV( I ), MM, M, WORK( IWRK1 ),
  775. $ LWORK-IWRK1+1, IWORK, IERR )
  776. *
  777. 20 CONTINUE
  778. END IF
  779. END IF
  780. *
  781. * Undo balancing on VL and VR and normalization
  782. * (Workspace: none needed)
  783. *
  784. IF( ILVL ) THEN
  785. CALL DGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
  786. $ LDVL, IERR )
  787. *
  788. DO 70 JC = 1, N
  789. IF( ALPHAI( JC ).LT.ZERO )
  790. $ GO TO 70
  791. TEMP = ZERO
  792. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  793. DO 30 JR = 1, N
  794. TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
  795. 30 CONTINUE
  796. ELSE
  797. DO 40 JR = 1, N
  798. TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
  799. $ ABS( VL( JR, JC+1 ) ) )
  800. 40 CONTINUE
  801. END IF
  802. IF( TEMP.LT.SMLNUM )
  803. $ GO TO 70
  804. TEMP = ONE / TEMP
  805. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  806. DO 50 JR = 1, N
  807. VL( JR, JC ) = VL( JR, JC )*TEMP
  808. 50 CONTINUE
  809. ELSE
  810. DO 60 JR = 1, N
  811. VL( JR, JC ) = VL( JR, JC )*TEMP
  812. VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
  813. 60 CONTINUE
  814. END IF
  815. 70 CONTINUE
  816. END IF
  817. IF( ILVR ) THEN
  818. CALL DGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
  819. $ LDVR, IERR )
  820. DO 120 JC = 1, N
  821. IF( ALPHAI( JC ).LT.ZERO )
  822. $ GO TO 120
  823. TEMP = ZERO
  824. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  825. DO 80 JR = 1, N
  826. TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
  827. 80 CONTINUE
  828. ELSE
  829. DO 90 JR = 1, N
  830. TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
  831. $ ABS( VR( JR, JC+1 ) ) )
  832. 90 CONTINUE
  833. END IF
  834. IF( TEMP.LT.SMLNUM )
  835. $ GO TO 120
  836. TEMP = ONE / TEMP
  837. IF( ALPHAI( JC ).EQ.ZERO ) THEN
  838. DO 100 JR = 1, N
  839. VR( JR, JC ) = VR( JR, JC )*TEMP
  840. 100 CONTINUE
  841. ELSE
  842. DO 110 JR = 1, N
  843. VR( JR, JC ) = VR( JR, JC )*TEMP
  844. VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
  845. 110 CONTINUE
  846. END IF
  847. 120 CONTINUE
  848. END IF
  849. *
  850. * Undo scaling if necessary
  851. *
  852. IF( ILASCL ) THEN
  853. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
  854. CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
  855. END IF
  856. *
  857. IF( ILBSCL ) THEN
  858. CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  859. END IF
  860. *
  861. 130 CONTINUE
  862. WORK( 1 ) = MAXWRK
  863. *
  864. RETURN
  865. *
  866. * End of DGGEVX
  867. *
  868. END