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

/contrib/clapack-3.2.1-CMAKE/SRC/dggevx.c

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