PageRenderTime 48ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 0ms

/structural/supporting_matl_misc_and_old/support/converted/lapack/dggevx.c

https://bitbucket.org/slawton/windturbinemdo-stevebitb
C | 1027 lines | 467 code | 89 blank | 471 comment | 146 complexity | 6099438bd0495c04ba64d7fcc986d7d5 MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause
  1. /* dggevx.f -- translated by f2c (version 20000121).
  2. You must link the resulting object file with the libraries:
  3. -lf2c -lm (in that order)
  4. */
  5. #include "f2c.h"
  6. /* Table of constant values */
  7. static integer c__1 = 1;
  8. static integer c__0 = 0;
  9. static doublereal c_b59 = 0.;
  10. static doublereal c_b60 = 1.;
  11. /* > \brief <b> DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE mat
  12. rices</b> */
  13. /* =========== DOCUMENTATION =========== */
  14. /* Online html documentation available at */
  15. /* http://www.netlib.org/lapack/explore-html/ */
  16. /* > \htmlonly */
  17. /* > Download DGGEVX + dependencies */
  18. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggevx.
  19. f"> */
  20. /* > [TGZ]</a> */
  21. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggevx.
  22. f"> */
  23. /* > [ZIP]</a> */
  24. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggevx.
  25. f"> */
  26. /* > [TXT]</a> */
  27. /* > \endhtmlonly */
  28. /* Definition: */
  29. /* =========== */
  30. /* SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, */
  31. /* ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO, */
  32. /* IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, */
  33. /* RCONDV, WORK, LWORK, IWORK, BWORK, INFO ) */
  34. /* .. Scalar Arguments .. */
  35. /* CHARACTER BALANC, JOBVL, JOBVR, SENSE */
  36. /* INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N */
  37. /* DOUBLE PRECISION ABNRM, BBNRM */
  38. /* .. */
  39. /* .. Array Arguments .. */
  40. /* LOGICAL BWORK( * ) */
  41. /* INTEGER IWORK( * ) */
  42. /* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), */
  43. /* $ B( LDB, * ), BETA( * ), LSCALE( * ), */
  44. /* $ RCONDE( * ), RCONDV( * ), RSCALE( * ), */
  45. /* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * ) */
  46. /* .. */
  47. /* > \par Purpose: */
  48. /* ============= */
  49. /* > */
  50. /* > \verbatim */
  51. /* > */
  52. /* > DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B) */
  53. /* > the generalized eigenvalues, and optionally, the left and/or right */
  54. /* > generalized eigenvectors. */
  55. /* > */
  56. /* > Optionally also, it computes a balancing transformation to improve */
  57. /* > the conditioning of the eigenvalues and eigenvectors (ILO, IHI, */
  58. /* > LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for */
  59. /* > the eigenvalues (RCONDE), and reciprocal condition numbers for the */
  60. /* > right eigenvectors (RCONDV). */
  61. /* > */
  62. /* > A generalized eigenvalue for a pair of matrices (A,B) is a scalar */
  63. /* > lambda or a ratio alpha/beta = lambda, such that A - lambda*B is */
  64. /* > singular. It is usually represented as the pair (alpha,beta), as */
  65. /* > there is a reasonable interpretation for beta=0, and even for both */
  66. /* > being zero. */
  67. /* > */
  68. /* > The right eigenvector v(j) corresponding to the eigenvalue lambda(j) */
  69. /* > of (A,B) satisfies */
  70. /* > */
  71. /* > A * v(j) = lambda(j) * B * v(j) . */
  72. /* > */
  73. /* > The left eigenvector u(j) corresponding to the eigenvalue lambda(j) */
  74. /* > of (A,B) satisfies */
  75. /* > */
  76. /* > u(j)**H * A = lambda(j) * u(j)**H * B. */
  77. /* > */
  78. /* > where u(j)**H is the conjugate-transpose of u(j). */
  79. /* > */
  80. /* > \endverbatim */
  81. /* Arguments: */
  82. /* ========== */
  83. /* > \param[in] BALANC */
  84. /* > \verbatim */
  85. /* > BALANC is CHARACTER*1 */
  86. /* > Specifies the balance option to be performed. */
  87. /* > = 'N': do not diagonally scale or permute; */
  88. /* > = 'P': permute only; */
  89. /* > = 'S': scale only; */
  90. /* > = 'B': both permute and scale. */
  91. /* > Computed reciprocal condition numbers will be for the */
  92. /* > matrices after permuting and/or balancing. Permuting does */
  93. /* > not change condition numbers (in exact arithmetic), but */
  94. /* > balancing does. */
  95. /* > \endverbatim */
  96. /* > */
  97. /* > \param[in] JOBVL */
  98. /* > \verbatim */
  99. /* > JOBVL is CHARACTER*1 */
  100. /* > = 'N': do not compute the left generalized eigenvectors; */
  101. /* > = 'V': compute the left generalized eigenvectors. */
  102. /* > \endverbatim */
  103. /* > */
  104. /* > \param[in] JOBVR */
  105. /* > \verbatim */
  106. /* > JOBVR is CHARACTER*1 */
  107. /* > = 'N': do not compute the right generalized eigenvectors; */
  108. /* > = 'V': compute the right generalized eigenvectors. */
  109. /* > \endverbatim */
  110. /* > */
  111. /* > \param[in] SENSE */
  112. /* > \verbatim */
  113. /* > SENSE is CHARACTER*1 */
  114. /* > Determines which reciprocal condition numbers are computed. */
  115. /* > = 'N': none are computed; */
  116. /* > = 'E': computed for eigenvalues only; */
  117. /* > = 'V': computed for eigenvectors only; */
  118. /* > = 'B': computed for eigenvalues and eigenvectors. */
  119. /* > \endverbatim */
  120. /* > */
  121. /* > \param[in] N */
  122. /* > \verbatim */
  123. /* > N is INTEGER */
  124. /* > The order of the matrices A, B, VL, and VR. N >= 0. */
  125. /* > \endverbatim */
  126. /* > */
  127. /* > \param[in,out] A */
  128. /* > \verbatim */
  129. /* > A is DOUBLE PRECISION array, dimension (LDA, N) */
  130. /* > On entry, the matrix A in the pair (A,B). */
  131. /* > On exit, A has been overwritten. If JOBVL='V' or JOBVR='V' */
  132. /* > or both, then A contains the first part of the real Schur */
  133. /* > form of the "balanced" versions of the input A and B. */
  134. /* > \endverbatim */
  135. /* > */
  136. /* > \param[in] LDA */
  137. /* > \verbatim */
  138. /* > LDA is INTEGER */
  139. /* > The leading dimension of A. LDA >= max(1,N). */
  140. /* > \endverbatim */
  141. /* > */
  142. /* > \param[in,out] B */
  143. /* > \verbatim */
  144. /* > B is DOUBLE PRECISION array, dimension (LDB, N) */
  145. /* > On entry, the matrix B in the pair (A,B). */
  146. /* > On exit, B has been overwritten. If JOBVL='V' or JOBVR='V' */
  147. /* > or both, then B contains the second part of the real Schur */
  148. /* > form of the "balanced" versions of the input A and B. */
  149. /* > \endverbatim */
  150. /* > */
  151. /* > \param[in] LDB */
  152. /* > \verbatim */
  153. /* > LDB is INTEGER */
  154. /* > The leading dimension of B. LDB >= max(1,N). */
  155. /* > \endverbatim */
  156. /* > */
  157. /* > \param[out] ALPHAR */
  158. /* > \verbatim */
  159. /* > ALPHAR is DOUBLE PRECISION array, dimension (N) */
  160. /* > \endverbatim */
  161. /* > */
  162. /* > \param[out] ALPHAI */
  163. /* > \verbatim */
  164. /* > ALPHAI is DOUBLE PRECISION array, dimension (N) */
  165. /* > \endverbatim */
  166. /* > */
  167. /* > \param[out] BETA */
  168. /* > \verbatim */
  169. /* > BETA is DOUBLE PRECISION array, dimension (N) */
  170. /* > On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will */
  171. /* > be the generalized eigenvalues. If ALPHAI(j) is zero, then */
  172. /* > the j-th eigenvalue is real; if positive, then the j-th and */
  173. /* > (j+1)-st eigenvalues are a complex conjugate pair, with */
  174. /* > ALPHAI(j+1) negative. */
  175. /* > */
  176. /* > Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) */
  177. /* > may easily over- or underflow, and BETA(j) may even be zero. */
  178. /* > Thus, the user should avoid naively computing the ratio */
  179. /* > ALPHA/BETA. However, ALPHAR and ALPHAI will be always less */
  180. /* > than and usually comparable with norm(A) in magnitude, and */
  181. /* > BETA always less than and usually comparable with norm(B). */
  182. /* > \endverbatim */
  183. /* > */
  184. /* > \param[out] VL */
  185. /* > \verbatim */
  186. /* > VL is DOUBLE PRECISION array, dimension (LDVL,N) */
  187. /* > If JOBVL = 'V', the left eigenvectors u(j) are stored one */
  188. /* > after another in the columns of VL, in the same order as */
  189. /* > their eigenvalues. If the j-th eigenvalue is real, then */
  190. /* > u(j) = VL(:,j), the j-th column of VL. If the j-th and */
  191. /* > (j+1)-th eigenvalues form a complex conjugate pair, then */
  192. /* > u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). */
  193. /* > Each eigenvector will be scaled so the largest component have */
  194. /* > abs(real part) + abs(imag. part) = 1. */
  195. /* > Not referenced if JOBVL = 'N'. */
  196. /* > \endverbatim */
  197. /* > */
  198. /* > \param[in] LDVL */
  199. /* > \verbatim */
  200. /* > LDVL is INTEGER */
  201. /* > The leading dimension of the matrix VL. LDVL >= 1, and */
  202. /* > if JOBVL = 'V', LDVL >= N. */
  203. /* > \endverbatim */
  204. /* > */
  205. /* > \param[out] VR */
  206. /* > \verbatim */
  207. /* > VR is DOUBLE PRECISION array, dimension (LDVR,N) */
  208. /* > If JOBVR = 'V', the right eigenvectors v(j) are stored one */
  209. /* > after another in the columns of VR, in the same order as */
  210. /* > their eigenvalues. If the j-th eigenvalue is real, then */
  211. /* > v(j) = VR(:,j), the j-th column of VR. If the j-th and */
  212. /* > (j+1)-th eigenvalues form a complex conjugate pair, then */
  213. /* > v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). */
  214. /* > Each eigenvector will be scaled so the largest component have */
  215. /* > abs(real part) + abs(imag. part) = 1. */
  216. /* > Not referenced if JOBVR = 'N'. */
  217. /* > \endverbatim */
  218. /* > */
  219. /* > \param[in] LDVR */
  220. /* > \verbatim */
  221. /* > LDVR is INTEGER */
  222. /* > The leading dimension of the matrix VR. LDVR >= 1, and */
  223. /* > if JOBVR = 'V', LDVR >= N. */
  224. /* > \endverbatim */
  225. /* > */
  226. /* > \param[out] ILO */
  227. /* > \verbatim */
  228. /* > ILO is INTEGER */
  229. /* > \endverbatim */
  230. /* > */
  231. /* > \param[out] IHI */
  232. /* > \verbatim */
  233. /* > IHI is INTEGER */
  234. /* > ILO and IHI are integer values such that on exit */
  235. /* > A(i,j) = 0 and B(i,j) = 0 if i > j and */
  236. /* > j = 1,...,ILO-1 or i = IHI+1,...,N. */
  237. /* > If BALANC = 'N' or 'S', ILO = 1 and IHI = N. */
  238. /* > \endverbatim */
  239. /* > */
  240. /* > \param[out] LSCALE */
  241. /* > \verbatim */
  242. /* > LSCALE is DOUBLE PRECISION array, dimension (N) */
  243. /* > Details of the permutations and scaling factors applied */
  244. /* > to the left side of A and B. If PL(j) is the index of the */
  245. /* > row interchanged with row j, and DL(j) is the scaling */
  246. /* > factor applied to row j, then */
  247. /* > LSCALE(j) = PL(j) for j = 1,...,ILO-1 */
  248. /* > = DL(j) for j = ILO,...,IHI */
  249. /* > = PL(j) for j = IHI+1,...,N. */
  250. /* > The order in which the interchanges are made is N to IHI+1, */
  251. /* > then 1 to ILO-1. */
  252. /* > \endverbatim */
  253. /* > */
  254. /* > \param[out] RSCALE */
  255. /* > \verbatim */
  256. /* > RSCALE is DOUBLE PRECISION array, dimension (N) */
  257. /* > Details of the permutations and scaling factors applied */
  258. /* > to the right side of A and B. If PR(j) is the index of the */
  259. /* > column interchanged with column j, and DR(j) is the scaling */
  260. /* > factor applied to column j, then */
  261. /* > RSCALE(j) = PR(j) for j = 1,...,ILO-1 */
  262. /* > = DR(j) for j = ILO,...,IHI */
  263. /* > = PR(j) for j = IHI+1,...,N */
  264. /* > The order in which the interchanges are made is N to IHI+1, */
  265. /* > then 1 to ILO-1. */
  266. /* > \endverbatim */
  267. /* > */
  268. /* > \param[out] ABNRM */
  269. /* > \verbatim */
  270. /* > ABNRM is DOUBLE PRECISION */
  271. /* > The one-norm of the balanced matrix A. */
  272. /* > \endverbatim */
  273. /* > */
  274. /* > \param[out] BBNRM */
  275. /* > \verbatim */
  276. /* > BBNRM is DOUBLE PRECISION */
  277. /* > The one-norm of the balanced matrix B. */
  278. /* > \endverbatim */
  279. /* > */
  280. /* > \param[out] RCONDE */
  281. /* > \verbatim */
  282. /* > RCONDE is DOUBLE PRECISION array, dimension (N) */
  283. /* > If SENSE = 'E' or 'B', the reciprocal condition numbers of */
  284. /* > the eigenvalues, stored in consecutive elements of the array. */
  285. /* > For a complex conjugate pair of eigenvalues two consecutive */
  286. /* > elements of RCONDE are set to the same value. Thus RCONDE(j), */
  287. /* > RCONDV(j), and the j-th columns of VL and VR all correspond */
  288. /* > to the j-th eigenpair. */
  289. /* > If SENSE = 'N or 'V', RCONDE is not referenced. */
  290. /* > \endverbatim */
  291. /* > */
  292. /* > \param[out] RCONDV */
  293. /* > \verbatim */
  294. /* > RCONDV is DOUBLE PRECISION array, dimension (N) */
  295. /* > If SENSE = 'V' or 'B', the estimated reciprocal condition */
  296. /* > numbers of the eigenvectors, stored in consecutive elements */
  297. /* > of the array. For a complex eigenvector two consecutive */
  298. /* > elements of RCONDV are set to the same value. If the */
  299. /* > eigenvalues cannot be reordered to compute RCONDV(j), */
  300. /* > RCONDV(j) is set to 0; this can only occur when the true */
  301. /* > value would be very small anyway. */
  302. /* > If SENSE = 'N' or 'E', RCONDV is not referenced. */
  303. /* > \endverbatim */
  304. /* > */
  305. /* > \param[out] WORK */
  306. /* > \verbatim */
  307. /* > WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
  308. /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  309. /* > \endverbatim */
  310. /* > */
  311. /* > \param[in] LWORK */
  312. /* > \verbatim */
  313. /* > LWORK is INTEGER */
  314. /* > The dimension of the array WORK. LWORK >= max(1,2*N). */
  315. /* > If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V', */
  316. /* > LWORK >= max(1,6*N). */
  317. /* > If SENSE = 'E' or 'B', LWORK >= max(1,10*N). */
  318. /* > If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16. */
  319. /* > */
  320. /* > If LWORK = -1, then a workspace query is assumed; the routine */
  321. /* > only calculates the optimal size of the WORK array, returns */
  322. /* > this value as the first entry of the WORK array, and no error */
  323. /* > message related to LWORK is issued by XERBLA. */
  324. /* > \endverbatim */
  325. /* > */
  326. /* > \param[out] IWORK */
  327. /* > \verbatim */
  328. /* > IWORK is INTEGER array, dimension (N+6) */
  329. /* > If SENSE = 'E', IWORK is not referenced. */
  330. /* > \endverbatim */
  331. /* > */
  332. /* > \param[out] BWORK */
  333. /* > \verbatim */
  334. /* > BWORK is LOGICAL array, dimension (N) */
  335. /* > If SENSE = 'N', BWORK is not referenced. */
  336. /* > \endverbatim */
  337. /* > */
  338. /* > \param[out] INFO */
  339. /* > \verbatim */
  340. /* > INFO is INTEGER */
  341. /* > = 0: successful exit */
  342. /* > < 0: if INFO = -i, the i-th argument had an illegal value. */
  343. /* > = 1,...,N: */
  344. /* > The QZ iteration failed. No eigenvectors have been */
  345. /* > calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) */
  346. /* > should be correct for j=INFO+1,...,N. */
  347. /* > > N: =N+1: other than QZ iteration failed in DHGEQZ. */
  348. /* > =N+2: error return from DTGEVC. */
  349. /* > \endverbatim */
  350. /* Authors: */
  351. /* ======== */
  352. /* > \author Univ. of Tennessee */
  353. /* > \author Univ. of California Berkeley */
  354. /* > \author Univ. of Colorado Denver */
  355. /* > \author NAG Ltd. */
  356. /* > \date November 2011 */
  357. /* > \ingroup doubleGEeigen */
  358. /* > \par Further Details: */
  359. /* ===================== */
  360. /* > */
  361. /* > \verbatim */
  362. /* > */
  363. /* > Balancing a matrix pair (A,B) includes, first, permuting rows and */
  364. /* > columns to isolate eigenvalues, second, applying diagonal similarity */
  365. /* > transformation to the rows and columns to make the rows and columns */
  366. /* > as close in norm as possible. The computed reciprocal condition */
  367. /* > numbers correspond to the balanced matrix. Permuting rows and columns */
  368. /* > will not change the condition numbers (in exact arithmetic) but */
  369. /* > diagonal scaling will. For further explanation of balancing, see */
  370. /* > section 4.11.1.2 of LAPACK Users' Guide. */
  371. /* > */
  372. /* > An approximate error bound on the chordal distance between the i-th */
  373. /* > computed generalized eigenvalue w and the corresponding exact */
  374. /* > eigenvalue lambda is */
  375. /* > */
  376. /* > chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I) */
  377. /* > */
  378. /* > An approximate error bound for the angle between the i-th computed */
  379. /* > eigenvector VL(i) or VR(i) is given by */
  380. /* > */
  381. /* > EPS * norm(ABNRM, BBNRM) / DIF(i). */
  382. /* > */
  383. /* > For further explanation of the reciprocal condition numbers RCONDE */
  384. /* > and RCONDV, see section 4.11 of LAPACK User's Guide. */
  385. /* > \endverbatim */
  386. /* > */
  387. /* ===================================================================== */
  388. /* Subroutine */ int dggevx_(balanc, jobvl, jobvr, sense, n, a, lda, b, ldb,
  389. alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale,
  390. abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info,
  391. balanc_len, jobvl_len, jobvr_len, sense_len)
  392. char *balanc, *jobvl, *jobvr, *sense;
  393. integer *n;
  394. doublereal *a;
  395. integer *lda;
  396. doublereal *b;
  397. integer *ldb;
  398. doublereal *alphar, *alphai, *beta, *vl;
  399. integer *ldvl;
  400. doublereal *vr;
  401. integer *ldvr, *ilo, *ihi;
  402. doublereal *lscale, *rscale, *abnrm, *bbnrm, *rconde, *rcondv, *work;
  403. integer *lwork, *iwork;
  404. logical *bwork;
  405. integer *info;
  406. ftnlen balanc_len;
  407. ftnlen jobvl_len;
  408. ftnlen jobvr_len;
  409. ftnlen sense_len;
  410. {
  411. /* System generated locals */
  412. integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
  413. vr_offset, i__1, i__2;
  414. doublereal d__1, d__2, d__3, d__4;
  415. /* Builtin functions */
  416. double sqrt();
  417. /* Local variables */
  418. static logical pair;
  419. static doublereal anrm, bnrm;
  420. static integer ierr, itau;
  421. static doublereal temp;
  422. static logical ilvl, ilvr;
  423. static integer iwrk, iwrk1, i__, j, m;
  424. extern logical lsame_();
  425. static integer icols;
  426. static logical noscl;
  427. static integer irows;
  428. extern /* Subroutine */ int dlabad_();
  429. static integer jc;
  430. extern /* Subroutine */ int dggbak_(), dggbal_();
  431. static integer in;
  432. extern doublereal dlamch_();
  433. static integer mm;
  434. extern doublereal dlange_();
  435. static integer jr;
  436. extern /* Subroutine */ int dgghrd_(), dlascl_();
  437. static logical ilascl, ilbscl;
  438. extern /* Subroutine */ int dgeqrf_(), dlacpy_();
  439. static logical ldumma[1];
  440. static char chtemp[1];
  441. static doublereal bignum;
  442. extern /* Subroutine */ int dhgeqz_(), dlaset_();
  443. static integer ijobvl;
  444. extern /* Subroutine */ int dtgevc_(), dtgsna_(), xerbla_();
  445. extern integer ilaenv_();
  446. static integer ijobvr;
  447. static logical wantsb;
  448. extern /* Subroutine */ int dorgqr_();
  449. static doublereal anrmto;
  450. static logical wantse;
  451. static doublereal bnrmto;
  452. extern /* Subroutine */ int dormqr_();
  453. static integer minwrk, maxwrk;
  454. static logical wantsn;
  455. static doublereal smlnum;
  456. static logical lquery, wantsv;
  457. static doublereal eps;
  458. static logical ilv;
  459. /* -- LAPACK driver routine (version 3.4.0) -- */
  460. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  461. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  462. /* November 2011 */
  463. /* .. Scalar Arguments .. */
  464. /* .. */
  465. /* .. Array Arguments .. */
  466. /* .. */
  467. /* ===================================================================== */
  468. /* .. Parameters .. */
  469. /* .. */
  470. /* .. Local Scalars .. */
  471. /* .. */
  472. /* .. Local Arrays .. */
  473. /* .. */
  474. /* .. External Subroutines .. */
  475. /* .. */
  476. /* .. External Functions .. */
  477. /* .. */
  478. /* .. Intrinsic Functions .. */
  479. /* .. */
  480. /* .. Executable Statements .. */
  481. /* Decode the input arguments */
  482. /* Parameter adjustments */
  483. a_dim1 = *lda;
  484. a_offset = 1 + a_dim1 * 1;
  485. a -= a_offset;
  486. b_dim1 = *ldb;
  487. b_offset = 1 + b_dim1 * 1;
  488. b -= b_offset;
  489. --alphar;
  490. --alphai;
  491. --beta;
  492. vl_dim1 = *ldvl;
  493. vl_offset = 1 + vl_dim1 * 1;
  494. vl -= vl_offset;
  495. vr_dim1 = *ldvr;
  496. vr_offset = 1 + vr_dim1 * 1;
  497. vr -= vr_offset;
  498. --lscale;
  499. --rscale;
  500. --rconde;
  501. --rcondv;
  502. --work;
  503. --iwork;
  504. --bwork;
  505. /* Function Body */
  506. if (lsame_(jobvl, "N", (ftnlen)1, (ftnlen)1)) {
  507. ijobvl = 1;
  508. ilvl = FALSE_;
  509. } else if (lsame_(jobvl, "V", (ftnlen)1, (ftnlen)1)) {
  510. ijobvl = 2;
  511. ilvl = TRUE_;
  512. } else {
  513. ijobvl = -1;
  514. ilvl = FALSE_;
  515. }
  516. if (lsame_(jobvr, "N", (ftnlen)1, (ftnlen)1)) {
  517. ijobvr = 1;
  518. ilvr = FALSE_;
  519. } else if (lsame_(jobvr, "V", (ftnlen)1, (ftnlen)1)) {
  520. ijobvr = 2;
  521. ilvr = TRUE_;
  522. } else {
  523. ijobvr = -1;
  524. ilvr = FALSE_;
  525. }
  526. ilv = ilvl || ilvr;
  527. noscl = lsame_(balanc, "N", (ftnlen)1, (ftnlen)1) || lsame_(balanc, "P", (
  528. ftnlen)1, (ftnlen)1);
  529. wantsn = lsame_(sense, "N", (ftnlen)1, (ftnlen)1);
  530. wantse = lsame_(sense, "E", (ftnlen)1, (ftnlen)1);
  531. wantsv = lsame_(sense, "V", (ftnlen)1, (ftnlen)1);
  532. wantsb = lsame_(sense, "B", (ftnlen)1, (ftnlen)1);
  533. /* Test the input arguments */
  534. *info = 0;
  535. lquery = *lwork == -1;
  536. if (! (lsame_(balanc, "N", (ftnlen)1, (ftnlen)1) || lsame_(balanc, "S", (
  537. ftnlen)1, (ftnlen)1) || lsame_(balanc, "P", (ftnlen)1, (ftnlen)1)
  538. || lsame_(balanc, "B", (ftnlen)1, (ftnlen)1))) {
  539. *info = -1;
  540. } else if (ijobvl <= 0) {
  541. *info = -2;
  542. } else if (ijobvr <= 0) {
  543. *info = -3;
  544. } else if (! (wantsn || wantse || wantsb || wantsv)) {
  545. *info = -4;
  546. } else if (*n < 0) {
  547. *info = -5;
  548. } else if (*lda < max(1,*n)) {
  549. *info = -7;
  550. } else if (*ldb < max(1,*n)) {
  551. *info = -9;
  552. } else if (*ldvl < 1 || ilvl && *ldvl < *n) {
  553. *info = -14;
  554. } else if (*ldvr < 1 || ilvr && *ldvr < *n) {
  555. *info = -16;
  556. }
  557. /* Compute workspace */
  558. /* (Note: Comments in the code beginning "Workspace:" describe the */
  559. /* minimal amount of workspace needed at that point in the code, */
  560. /* as well as the preferred amount for good performance. */
  561. /* NB refers to the optimal block size for the immediately */
  562. /* following subroutine, as returned by ILAENV. The workspace is */
  563. /* computed assuming ILO = 1 and IHI = N, the worst case.) */
  564. if (*info == 0) {
  565. if (*n == 0) {
  566. minwrk = 1;
  567. maxwrk = 1;
  568. } else {
  569. if (noscl && ! ilv) {
  570. minwrk = *n << 1;
  571. } else {
  572. minwrk = *n * 6;
  573. }
  574. if (wantse || wantsb) {
  575. minwrk = *n * 10;
  576. }
  577. if (wantsv || wantsb) {
  578. /* Computing MAX */
  579. i__1 = minwrk, i__2 = (*n << 1) * (*n + 4) + 16;
  580. minwrk = max(i__1,i__2);
  581. }
  582. maxwrk = minwrk;
  583. /* Computing MAX */
  584. i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", n, &
  585. c__1, n, &c__0, (ftnlen)6, (ftnlen)1);
  586. maxwrk = max(i__1,i__2);
  587. /* Computing MAX */
  588. i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "DORMQR", " ", n, &
  589. c__1, n, &c__0, (ftnlen)6, (ftnlen)1);
  590. maxwrk = max(i__1,i__2);
  591. if (ilvl) {
  592. /* Computing MAX */
  593. i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "DORGQR",
  594. " ", n, &c__1, n, &c__0, (ftnlen)6, (ftnlen)1);
  595. maxwrk = max(i__1,i__2);
  596. }
  597. }
  598. work[1] = (doublereal) maxwrk;
  599. if (*lwork < minwrk && ! lquery) {
  600. *info = -26;
  601. }
  602. }
  603. if (*info != 0) {
  604. i__1 = -(*info);
  605. xerbla_("DGGEVX", &i__1, (ftnlen)6);
  606. return 0;
  607. } else if (lquery) {
  608. return 0;
  609. }
  610. /* Quick return if possible */
  611. if (*n == 0) {
  612. return 0;
  613. }
  614. /* Get machine constants */
  615. eps = dlamch_("P", (ftnlen)1);
  616. smlnum = dlamch_("S", (ftnlen)1);
  617. bignum = 1. / smlnum;
  618. dlabad_(&smlnum, &bignum);
  619. smlnum = sqrt(smlnum) / eps;
  620. bignum = 1. / smlnum;
  621. /* Scale A if max element outside range [SMLNUM,BIGNUM] */
  622. anrm = dlange_("M", n, n, &a[a_offset], lda, &work[1], (ftnlen)1);
  623. ilascl = FALSE_;
  624. if (anrm > 0. && anrm < smlnum) {
  625. anrmto = smlnum;
  626. ilascl = TRUE_;
  627. } else if (anrm > bignum) {
  628. anrmto = bignum;
  629. ilascl = TRUE_;
  630. }
  631. if (ilascl) {
  632. dlascl_("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
  633. ierr, (ftnlen)1);
  634. }
  635. /* Scale B if max element outside range [SMLNUM,BIGNUM] */
  636. bnrm = dlange_("M", n, n, &b[b_offset], ldb, &work[1], (ftnlen)1);
  637. ilbscl = FALSE_;
  638. if (bnrm > 0. && bnrm < smlnum) {
  639. bnrmto = smlnum;
  640. ilbscl = TRUE_;
  641. } else if (bnrm > bignum) {
  642. bnrmto = bignum;
  643. ilbscl = TRUE_;
  644. }
  645. if (ilbscl) {
  646. dlascl_("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
  647. ierr, (ftnlen)1);
  648. }
  649. /* Permute and/or balance the matrix pair (A,B) */
  650. /* (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise) */
  651. dggbal_(balanc, n, &a[a_offset], lda, &b[b_offset], ldb, ilo, ihi, &
  652. lscale[1], &rscale[1], &work[1], &ierr, (ftnlen)1);
  653. /* Compute ABNRM and BBNRM */
  654. *abnrm = dlange_("1", n, n, &a[a_offset], lda, &work[1], (ftnlen)1);
  655. if (ilascl) {
  656. work[1] = *abnrm;
  657. dlascl_("G", &c__0, &c__0, &anrmto, &anrm, &c__1, &c__1, &work[1], &
  658. c__1, &ierr, (ftnlen)1);
  659. *abnrm = work[1];
  660. }
  661. *bbnrm = dlange_("1", n, n, &b[b_offset], ldb, &work[1], (ftnlen)1);
  662. if (ilbscl) {
  663. work[1] = *bbnrm;
  664. dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, &c__1, &c__1, &work[1], &
  665. c__1, &ierr, (ftnlen)1);
  666. *bbnrm = work[1];
  667. }
  668. /* Reduce B to triangular form (QR decomposition of B) */
  669. /* (Workspace: need N, prefer N*NB ) */
  670. irows = *ihi + 1 - *ilo;
  671. if (ilv || ! wantsn) {
  672. icols = *n + 1 - *ilo;
  673. } else {
  674. icols = irows;
  675. }
  676. itau = 1;
  677. iwrk = itau + irows;
  678. i__1 = *lwork + 1 - iwrk;
  679. dgeqrf_(&irows, &icols, &b[*ilo + *ilo * b_dim1], ldb, &work[itau], &work[
  680. iwrk], &i__1, &ierr);
  681. /* Apply the orthogonal transformation to A */
  682. /* (Workspace: need N, prefer N*NB) */
  683. i__1 = *lwork + 1 - iwrk;
  684. dormqr_("L", "T", &irows, &icols, &irows, &b[*ilo + *ilo * b_dim1], ldb, &
  685. work[itau], &a[*ilo + *ilo * a_dim1], lda, &work[iwrk], &i__1, &
  686. ierr, (ftnlen)1, (ftnlen)1);
  687. /* Initialize VL and/or VR */
  688. /* (Workspace: need N, prefer N*NB) */
  689. if (ilvl) {
  690. dlaset_("Full", n, n, &c_b59, &c_b60, &vl[vl_offset], ldvl, (ftnlen)4)
  691. ;
  692. if (irows > 1) {
  693. i__1 = irows - 1;
  694. i__2 = irows - 1;
  695. dlacpy_("L", &i__1, &i__2, &b[*ilo + 1 + *ilo * b_dim1], ldb, &vl[
  696. *ilo + 1 + *ilo * vl_dim1], ldvl, (ftnlen)1);
  697. }
  698. i__1 = *lwork + 1 - iwrk;
  699. dorgqr_(&irows, &irows, &irows, &vl[*ilo + *ilo * vl_dim1], ldvl, &
  700. work[itau], &work[iwrk], &i__1, &ierr);
  701. }
  702. if (ilvr) {
  703. dlaset_("Full", n, n, &c_b59, &c_b60, &vr[vr_offset], ldvr, (ftnlen)4)
  704. ;
  705. }
  706. /* Reduce to generalized Hessenberg form */
  707. /* (Workspace: none needed) */
  708. if (ilv || ! wantsn) {
  709. /* Eigenvectors requested -- work on whole matrix. */
  710. dgghrd_(jobvl, jobvr, n, ilo, ihi, &a[a_offset], lda, &b[b_offset],
  711. ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ierr, (
  712. ftnlen)1, (ftnlen)1);
  713. } else {
  714. dgghrd_("N", "N", &irows, &c__1, &irows, &a[*ilo + *ilo * a_dim1],
  715. lda, &b[*ilo + *ilo * b_dim1], ldb, &vl[vl_offset], ldvl, &vr[
  716. vr_offset], ldvr, &ierr, (ftnlen)1, (ftnlen)1);
  717. }
  718. /* Perform QZ algorithm (Compute eigenvalues, and optionally, the */
  719. /* Schur forms and Schur vectors) */
  720. /* (Workspace: need N) */
  721. if (ilv || ! wantsn) {
  722. *(unsigned char *)chtemp = 'S';
  723. } else {
  724. *(unsigned char *)chtemp = 'E';
  725. }
  726. dhgeqz_(chtemp, jobvl, jobvr, n, ilo, ihi, &a[a_offset], lda, &b[b_offset]
  727. , ldb, &alphar[1], &alphai[1], &beta[1], &vl[vl_offset], ldvl, &
  728. vr[vr_offset], ldvr, &work[1], lwork, &ierr, (ftnlen)1, (ftnlen)1,
  729. (ftnlen)1);
  730. if (ierr != 0) {
  731. if (ierr > 0 && ierr <= *n) {
  732. *info = ierr;
  733. } else if (ierr > *n && ierr <= *n << 1) {
  734. *info = ierr - *n;
  735. } else {
  736. *info = *n + 1;
  737. }
  738. goto L130;
  739. }
  740. /* Compute Eigenvectors and estimate condition numbers if desired */
  741. /* (Workspace: DTGEVC: need 6*N */
  742. /* DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B', */
  743. /* need N otherwise ) */
  744. if (ilv || ! wantsn) {
  745. if (ilv) {
  746. if (ilvl) {
  747. if (ilvr) {
  748. *(unsigned char *)chtemp = 'B';
  749. } else {
  750. *(unsigned char *)chtemp = 'L';
  751. }
  752. } else {
  753. *(unsigned char *)chtemp = 'R';
  754. }
  755. dtgevc_(chtemp, "B", ldumma, n, &a[a_offset], lda, &b[b_offset],
  756. ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, n, &in, &
  757. work[1], &ierr, (ftnlen)1, (ftnlen)1);
  758. if (ierr != 0) {
  759. *info = *n + 2;
  760. goto L130;
  761. }
  762. }
  763. if (! wantsn) {
  764. /* compute eigenvectors (DTGEVC) and estimate condition */
  765. /* numbers (DTGSNA). Note that the definition of the condition */
  766. /* number is not invariant under transformation (u,v) to */
  767. /* (Q*u, Z*v), where (u,v) are eigenvectors of the generalized */
  768. /* Schur form (S,T), Q and Z are orthogonal matrices. In order */
  769. /* to avoid using extra 2*N*N workspace, we have to recalculate */
  770. /* eigenvectors and estimate one condition numbers at a time. */
  771. pair = FALSE_;
  772. i__1 = *n;
  773. for (i__ = 1; i__ <= i__1; ++i__) {
  774. if (pair) {
  775. pair = FALSE_;
  776. goto L20;
  777. }
  778. mm = 1;
  779. if (i__ < *n) {
  780. if (a[i__ + 1 + i__ * a_dim1] != 0.) {
  781. pair = TRUE_;
  782. mm = 2;
  783. }
  784. }
  785. i__2 = *n;
  786. for (j = 1; j <= i__2; ++j) {
  787. bwork[j] = FALSE_;
  788. /* L10: */
  789. }
  790. if (mm == 1) {
  791. bwork[i__] = TRUE_;
  792. } else if (mm == 2) {
  793. bwork[i__] = TRUE_;
  794. bwork[i__ + 1] = TRUE_;
  795. }
  796. iwrk = mm * *n + 1;
  797. iwrk1 = iwrk + mm * *n;
  798. /* Compute a pair of left and right eigenvectors. */
  799. /* (compute workspace: need up to 4*N + 6*N) */
  800. if (wantse || wantsb) {
  801. dtgevc_("B", "S", &bwork[1], n, &a[a_offset], lda, &b[
  802. b_offset], ldb, &work[1], n, &work[iwrk], n, &mm,
  803. &m, &work[iwrk1], &ierr, (ftnlen)1, (ftnlen)1);
  804. if (ierr != 0) {
  805. *info = *n + 2;
  806. goto L130;
  807. }
  808. }
  809. i__2 = *lwork - iwrk1 + 1;
  810. dtgsna_(sense, "S", &bwork[1], n, &a[a_offset], lda, &b[
  811. b_offset], ldb, &work[1], n, &work[iwrk], n, &rconde[
  812. i__], &rcondv[i__], &mm, &m, &work[iwrk1], &i__2, &
  813. iwork[1], &ierr, (ftnlen)1, (ftnlen)1);
  814. L20:
  815. ;
  816. }
  817. }
  818. }
  819. /* Undo balancing on VL and VR and normalization */
  820. /* (Workspace: none needed) */
  821. if (ilvl) {
  822. dggbak_(balanc, "L", n, ilo, ihi, &lscale[1], &rscale[1], n, &vl[
  823. vl_offset], ldvl, &ierr, (ftnlen)1, (ftnlen)1);
  824. i__1 = *n;
  825. for (jc = 1; jc <= i__1; ++jc) {
  826. if (alphai[jc] < 0.) {
  827. goto L70;
  828. }
  829. temp = 0.;
  830. if (alphai[jc] == 0.) {
  831. i__2 = *n;
  832. for (jr = 1; jr <= i__2; ++jr) {
  833. /* Computing MAX */
  834. d__2 = temp, d__3 = (d__1 = vl[jr + jc * vl_dim1], abs(
  835. d__1));
  836. temp = max(d__2,d__3);
  837. /* L30: */
  838. }
  839. } else {
  840. i__2 = *n;
  841. for (jr = 1; jr <= i__2; ++jr) {
  842. /* Computing MAX */
  843. d__3 = temp, d__4 = (d__1 = vl[jr + jc * vl_dim1], abs(
  844. d__1)) + (d__2 = vl[jr + (jc + 1) * vl_dim1], abs(
  845. d__2));
  846. temp = max(d__3,d__4);
  847. /* L40: */
  848. }
  849. }
  850. if (temp < smlnum) {
  851. goto L70;
  852. }
  853. temp = 1. / temp;
  854. if (alphai[jc] == 0.) {
  855. i__2 = *n;
  856. for (jr = 1; jr <= i__2; ++jr) {
  857. vl[jr + jc * vl_dim1] *= temp;
  858. /* L50: */
  859. }
  860. } else {
  861. i__2 = *n;
  862. for (jr = 1; jr <= i__2; ++jr) {
  863. vl[jr + jc * vl_dim1] *= temp;
  864. vl[jr + (jc + 1) * vl_dim1] *= temp;
  865. /* L60: */
  866. }
  867. }
  868. L70:
  869. ;
  870. }
  871. }
  872. if (ilvr) {
  873. dggbak_(balanc, "R", n, ilo, ihi, &lscale[1], &rscale[1], n, &vr[
  874. vr_offset], ldvr, &ierr, (ftnlen)1, (ftnlen)1);
  875. i__1 = *n;
  876. for (jc = 1; jc <= i__1; ++jc) {
  877. if (alphai[jc] < 0.) {
  878. goto L120;
  879. }
  880. temp = 0.;
  881. if (alphai[jc] == 0.) {
  882. i__2 = *n;
  883. for (jr = 1; jr <= i__2; ++jr) {
  884. /* Computing MAX */
  885. d__2 = temp, d__3 = (d__1 = vr[jr + jc * vr_dim1], abs(
  886. d__1));
  887. temp = max(d__2,d__3);
  888. /* L80: */
  889. }
  890. } else {
  891. i__2 = *n;
  892. for (jr = 1; jr <= i__2; ++jr) {
  893. /* Computing MAX */
  894. d__3 = temp, d__4 = (d__1 = vr[jr + jc * vr_dim1], abs(
  895. d__1)) + (d__2 = vr[jr + (jc + 1) * vr_dim1], abs(
  896. d__2));
  897. temp = max(d__3,d__4);
  898. /* L90: */
  899. }
  900. }
  901. if (temp < smlnum) {
  902. goto L120;
  903. }
  904. temp = 1. / temp;
  905. if (alphai[jc] == 0.) {
  906. i__2 = *n;
  907. for (jr = 1; jr <= i__2; ++jr) {
  908. vr[jr + jc * vr_dim1] *= temp;
  909. /* L100: */
  910. }
  911. } else {
  912. i__2 = *n;
  913. for (jr = 1; jr <= i__2; ++jr) {
  914. vr[jr + jc * vr_dim1] *= temp;
  915. vr[jr + (jc + 1) * vr_dim1] *= temp;
  916. /* L110: */
  917. }
  918. }
  919. L120:
  920. ;
  921. }
  922. }
  923. /* Undo scaling if necessary */
  924. if (ilascl) {
  925. dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], n, &
  926. ierr, (ftnlen)1);
  927. dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], n, &
  928. ierr, (ftnlen)1);
  929. }
  930. if (ilbscl) {
  931. dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
  932. ierr, (ftnlen)1);
  933. }
  934. L130:
  935. work[1] = (doublereal) maxwrk;
  936. return 0;
  937. /* End of DGGEVX */
  938. } /* dggevx_ */