PageRenderTime 54ms CodeModel.GetById 27ms RepoModel.GetById 1ms app.codeStats 0ms

/lib/clapack/SRC/dggevx.c

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