dgeesx.c 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650
  1. /* dgeesx.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 integer c_n1 = -1;
  17. /* Subroutine */ int _starpu_dgeesx_(char *jobvs, char *sort, L_fp select, char *
  18. sense, integer *n, doublereal *a, integer *lda, integer *sdim,
  19. doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs,
  20. doublereal *rconde, doublereal *rcondv, doublereal *work, integer *
  21. lwork, integer *iwork, integer *liwork, logical *bwork, integer *info)
  22. {
  23. /* System generated locals */
  24. integer a_dim1, a_offset, vs_dim1, vs_offset, i__1, i__2, i__3;
  25. /* Builtin functions */
  26. double sqrt(doublereal);
  27. /* Local variables */
  28. integer i__, i1, i2, ip, ihi, ilo;
  29. doublereal dum[1], eps;
  30. integer ibal;
  31. doublereal anrm;
  32. integer ierr, itau, iwrk, lwrk, inxt, icond, ieval;
  33. extern logical _starpu_lsame_(char *, char *);
  34. extern /* Subroutine */ int _starpu_dcopy_(integer *, doublereal *, integer *,
  35. doublereal *, integer *), _starpu_dswap_(integer *, doublereal *, integer
  36. *, doublereal *, integer *);
  37. logical cursl;
  38. integer liwrk;
  39. extern /* Subroutine */ int _starpu_dlabad_(doublereal *, doublereal *), _starpu_dgebak_(
  40. char *, char *, integer *, integer *, integer *, doublereal *,
  41. integer *, doublereal *, integer *, integer *),
  42. _starpu_dgebal_(char *, integer *, doublereal *, integer *, integer *,
  43. integer *, doublereal *, integer *);
  44. logical lst2sl, scalea;
  45. extern doublereal _starpu_dlamch_(char *);
  46. doublereal cscale;
  47. extern doublereal _starpu_dlange_(char *, integer *, integer *, doublereal *,
  48. integer *, doublereal *);
  49. extern /* Subroutine */ int _starpu_dgehrd_(integer *, integer *, integer *,
  50. doublereal *, integer *, doublereal *, doublereal *, integer *,
  51. integer *), _starpu_dlascl_(char *, integer *, integer *, doublereal *,
  52. doublereal *, integer *, integer *, doublereal *, integer *,
  53. integer *), _starpu_dlacpy_(char *, integer *, integer *,
  54. doublereal *, integer *, doublereal *, integer *),
  55. _starpu_xerbla_(char *, integer *);
  56. extern integer _starpu_ilaenv_(integer *, char *, char *, integer *, integer *,
  57. integer *, integer *);
  58. doublereal bignum;
  59. extern /* Subroutine */ int _starpu_dorghr_(integer *, integer *, integer *,
  60. doublereal *, integer *, doublereal *, doublereal *, integer *,
  61. integer *), _starpu_dhseqr_(char *, char *, integer *, integer *, integer
  62. *, doublereal *, integer *, doublereal *, doublereal *,
  63. doublereal *, integer *, doublereal *, integer *, integer *);
  64. logical wantsb;
  65. extern /* Subroutine */ int _starpu_dtrsen_(char *, char *, logical *, integer *,
  66. doublereal *, integer *, doublereal *, integer *, doublereal *,
  67. doublereal *, integer *, doublereal *, doublereal *, doublereal *,
  68. integer *, integer *, integer *, integer *);
  69. logical wantse, lastsl;
  70. integer minwrk, maxwrk;
  71. logical wantsn;
  72. doublereal smlnum;
  73. integer hswork;
  74. logical wantst, lquery, wantsv, wantvs;
  75. /* -- LAPACK driver routine (version 3.2) -- */
  76. /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  77. /* November 2006 */
  78. /* .. Scalar Arguments .. */
  79. /* .. */
  80. /* .. Array Arguments .. */
  81. /* .. */
  82. /* .. Function Arguments .. */
  83. /* .. */
  84. /* Purpose */
  85. /* ======= */
  86. /* DGEESX computes for an N-by-N real nonsymmetric matrix A, the */
  87. /* eigenvalues, the real Schur form T, and, optionally, the matrix of */
  88. /* Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T). */
  89. /* Optionally, it also orders the eigenvalues on the diagonal of the */
  90. /* real Schur form so that selected eigenvalues are at the top left; */
  91. /* computes a reciprocal condition number for the average of the */
  92. /* selected eigenvalues (RCONDE); and computes a reciprocal condition */
  93. /* number for the right invariant subspace corresponding to the */
  94. /* selected eigenvalues (RCONDV). The leading columns of Z form an */
  95. /* orthonormal basis for this invariant subspace. */
  96. /* For further explanation of the reciprocal condition numbers RCONDE */
  97. /* and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where */
  98. /* these quantities are called s and sep respectively). */
  99. /* A real matrix is in real Schur form if it is upper quasi-triangular */
  100. /* with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in */
  101. /* the form */
  102. /* [ a b ] */
  103. /* [ c a ] */
  104. /* where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). */
  105. /* Arguments */
  106. /* ========= */
  107. /* JOBVS (input) CHARACTER*1 */
  108. /* = 'N': Schur vectors are not computed; */
  109. /* = 'V': Schur vectors are computed. */
  110. /* SORT (input) CHARACTER*1 */
  111. /* Specifies whether or not to order the eigenvalues on the */
  112. /* diagonal of the Schur form. */
  113. /* = 'N': Eigenvalues are not ordered; */
  114. /* = 'S': Eigenvalues are ordered (see SELECT). */
  115. /* SELECT (external procedure) LOGICAL FUNCTION of two DOUBLE PRECISION arguments */
  116. /* SELECT must be declared EXTERNAL in the calling subroutine. */
  117. /* If SORT = 'S', SELECT is used to select eigenvalues to sort */
  118. /* to the top left of the Schur form. */
  119. /* If SORT = 'N', SELECT is not referenced. */
  120. /* An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if */
  121. /* SELECT(WR(j),WI(j)) is true; i.e., if either one of a */
  122. /* complex conjugate pair of eigenvalues is selected, then both */
  123. /* are. Note that a selected complex eigenvalue may no longer */
  124. /* satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since */
  125. /* ordering may change the value of complex eigenvalues */
  126. /* (especially if the eigenvalue is ill-conditioned); in this */
  127. /* case INFO may be set to N+3 (see INFO below). */
  128. /* SENSE (input) CHARACTER*1 */
  129. /* Determines which reciprocal condition numbers are computed. */
  130. /* = 'N': None are computed; */
  131. /* = 'E': Computed for average of selected eigenvalues only; */
  132. /* = 'V': Computed for selected right invariant subspace only; */
  133. /* = 'B': Computed for both. */
  134. /* If SENSE = 'E', 'V' or 'B', SORT must equal 'S'. */
  135. /* N (input) INTEGER */
  136. /* The order of the matrix A. N >= 0. */
  137. /* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
  138. /* On entry, the N-by-N matrix A. */
  139. /* On exit, A is overwritten by its real Schur form T. */
  140. /* LDA (input) INTEGER */
  141. /* The leading dimension of the array A. LDA >= max(1,N). */
  142. /* SDIM (output) INTEGER */
  143. /* If SORT = 'N', SDIM = 0. */
  144. /* If SORT = 'S', SDIM = number of eigenvalues (after sorting) */
  145. /* for which SELECT is true. (Complex conjugate */
  146. /* pairs for which SELECT is true for either */
  147. /* eigenvalue count as 2.) */
  148. /* WR (output) DOUBLE PRECISION array, dimension (N) */
  149. /* WI (output) DOUBLE PRECISION array, dimension (N) */
  150. /* WR and WI contain the real and imaginary parts, respectively, */
  151. /* of the computed eigenvalues, in the same order that they */
  152. /* appear on the diagonal of the output Schur form T. Complex */
  153. /* conjugate pairs of eigenvalues appear consecutively with the */
  154. /* eigenvalue having the positive imaginary part first. */
  155. /* VS (output) DOUBLE PRECISION array, dimension (LDVS,N) */
  156. /* If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur */
  157. /* vectors. */
  158. /* If JOBVS = 'N', VS is not referenced. */
  159. /* LDVS (input) INTEGER */
  160. /* The leading dimension of the array VS. LDVS >= 1, and if */
  161. /* JOBVS = 'V', LDVS >= N. */
  162. /* RCONDE (output) DOUBLE PRECISION */
  163. /* If SENSE = 'E' or 'B', RCONDE contains the reciprocal */
  164. /* condition number for the average of the selected eigenvalues. */
  165. /* Not referenced if SENSE = 'N' or 'V'. */
  166. /* RCONDV (output) DOUBLE PRECISION */
  167. /* If SENSE = 'V' or 'B', RCONDV contains the reciprocal */
  168. /* condition number for the selected right invariant subspace. */
  169. /* Not referenced if SENSE = 'N' or 'E'. */
  170. /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
  171. /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  172. /* LWORK (input) INTEGER */
  173. /* The dimension of the array WORK. LWORK >= max(1,3*N). */
  174. /* Also, if SENSE = 'E' or 'V' or 'B', */
  175. /* LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of */
  176. /* selected eigenvalues computed by this routine. Note that */
  177. /* N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only */
  178. /* returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or */
  179. /* 'B' this may not be large enough. */
  180. /* For good performance, LWORK must generally be larger. */
  181. /* If LWORK = -1, then a workspace query is assumed; the routine */
  182. /* only calculates upper bounds on the optimal sizes of the */
  183. /* arrays WORK and IWORK, returns these values as the first */
  184. /* entries of the WORK and IWORK arrays, and no error messages */
  185. /* related to LWORK or LIWORK are issued by XERBLA. */
  186. /* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */
  187. /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
  188. /* LIWORK (input) INTEGER */
  189. /* The dimension of the array IWORK. */
  190. /* LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM). */
  191. /* Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is */
  192. /* only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this */
  193. /* may not be large enough. */
  194. /* If LIWORK = -1, then a workspace query is assumed; the */
  195. /* routine only calculates upper bounds on the optimal sizes of */
  196. /* the arrays WORK and IWORK, returns these values as the first */
  197. /* entries of the WORK and IWORK arrays, and no error messages */
  198. /* related to LWORK or LIWORK are issued by XERBLA. */
  199. /* BWORK (workspace) LOGICAL array, dimension (N) */
  200. /* Not referenced if SORT = 'N'. */
  201. /* INFO (output) INTEGER */
  202. /* = 0: successful exit */
  203. /* < 0: if INFO = -i, the i-th argument had an illegal value. */
  204. /* > 0: if INFO = i, and i is */
  205. /* <= N: the QR algorithm failed to compute all the */
  206. /* eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI */
  207. /* contain those eigenvalues which have converged; if */
  208. /* JOBVS = 'V', VS contains the transformation which */
  209. /* reduces A to its partially converged Schur form. */
  210. /* = N+1: the eigenvalues could not be reordered because some */
  211. /* eigenvalues were too close to separate (the problem */
  212. /* is very ill-conditioned); */
  213. /* = N+2: after reordering, roundoff changed values of some */
  214. /* complex eigenvalues so that leading eigenvalues in */
  215. /* the Schur form no longer satisfy SELECT=.TRUE. This */
  216. /* could also be caused by underflow due to scaling. */
  217. /* ===================================================================== */
  218. /* .. Parameters .. */
  219. /* .. */
  220. /* .. Local Scalars .. */
  221. /* .. */
  222. /* .. Local Arrays .. */
  223. /* .. */
  224. /* .. External Subroutines .. */
  225. /* .. */
  226. /* .. External Functions .. */
  227. /* .. */
  228. /* .. Intrinsic Functions .. */
  229. /* .. */
  230. /* .. Executable Statements .. */
  231. /* Test the input arguments */
  232. /* Parameter adjustments */
  233. a_dim1 = *lda;
  234. a_offset = 1 + a_dim1;
  235. a -= a_offset;
  236. --wr;
  237. --wi;
  238. vs_dim1 = *ldvs;
  239. vs_offset = 1 + vs_dim1;
  240. vs -= vs_offset;
  241. --work;
  242. --iwork;
  243. --bwork;
  244. /* Function Body */
  245. *info = 0;
  246. wantvs = _starpu_lsame_(jobvs, "V");
  247. wantst = _starpu_lsame_(sort, "S");
  248. wantsn = _starpu_lsame_(sense, "N");
  249. wantse = _starpu_lsame_(sense, "E");
  250. wantsv = _starpu_lsame_(sense, "V");
  251. wantsb = _starpu_lsame_(sense, "B");
  252. lquery = *lwork == -1 || *liwork == -1;
  253. if (! wantvs && ! _starpu_lsame_(jobvs, "N")) {
  254. *info = -1;
  255. } else if (! wantst && ! _starpu_lsame_(sort, "N")) {
  256. *info = -2;
  257. } else if (! (wantsn || wantse || wantsv || wantsb) || ! wantst && !
  258. wantsn) {
  259. *info = -4;
  260. } else if (*n < 0) {
  261. *info = -5;
  262. } else if (*lda < max(1,*n)) {
  263. *info = -7;
  264. } else if (*ldvs < 1 || wantvs && *ldvs < *n) {
  265. *info = -12;
  266. }
  267. /* Compute workspace */
  268. /* (Note: Comments in the code beginning "RWorkspace:" describe the */
  269. /* minimal amount of real workspace needed at that point in the */
  270. /* code, as well as the preferred amount for good performance. */
  271. /* IWorkspace refers to integer workspace. */
  272. /* NB refers to the optimal block size for the immediately */
  273. /* following subroutine, as returned by ILAENV. */
  274. /* HSWORK refers to the workspace preferred by DHSEQR, as */
  275. /* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, */
  276. /* the worst case. */
  277. /* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed */
  278. /* depends on SDIM, which is computed by the routine DTRSEN later */
  279. /* in the code.) */
  280. if (*info == 0) {
  281. liwrk = 1;
  282. if (*n == 0) {
  283. minwrk = 1;
  284. lwrk = 1;
  285. } else {
  286. maxwrk = (*n << 1) + *n * _starpu_ilaenv_(&c__1, "DGEHRD", " ", n, &c__1,
  287. n, &c__0);
  288. minwrk = *n * 3;
  289. _starpu_dhseqr_("S", jobvs, n, &c__1, n, &a[a_offset], lda, &wr[1], &wi[1]
  290. , &vs[vs_offset], ldvs, &work[1], &c_n1, &ieval);
  291. hswork = (integer) work[1];
  292. if (! wantvs) {
  293. /* Computing MAX */
  294. i__1 = maxwrk, i__2 = *n + hswork;
  295. maxwrk = max(i__1,i__2);
  296. } else {
  297. /* Computing MAX */
  298. i__1 = maxwrk, i__2 = (*n << 1) + (*n - 1) * _starpu_ilaenv_(&c__1,
  299. "DORGHR", " ", n, &c__1, n, &c_n1);
  300. maxwrk = max(i__1,i__2);
  301. /* Computing MAX */
  302. i__1 = maxwrk, i__2 = *n + hswork;
  303. maxwrk = max(i__1,i__2);
  304. }
  305. lwrk = maxwrk;
  306. if (! wantsn) {
  307. /* Computing MAX */
  308. i__1 = lwrk, i__2 = *n + *n * *n / 2;
  309. lwrk = max(i__1,i__2);
  310. }
  311. if (wantsv || wantsb) {
  312. liwrk = *n * *n / 4;
  313. }
  314. }
  315. iwork[1] = liwrk;
  316. work[1] = (doublereal) lwrk;
  317. if (*lwork < minwrk && ! lquery) {
  318. *info = -16;
  319. } else if (*liwork < 1 && ! lquery) {
  320. *info = -18;
  321. }
  322. }
  323. if (*info != 0) {
  324. i__1 = -(*info);
  325. _starpu_xerbla_("DGEESX", &i__1);
  326. return 0;
  327. }
  328. /* Quick return if possible */
  329. if (*n == 0) {
  330. *sdim = 0;
  331. return 0;
  332. }
  333. /* Get machine constants */
  334. eps = _starpu_dlamch_("P");
  335. smlnum = _starpu_dlamch_("S");
  336. bignum = 1. / smlnum;
  337. _starpu_dlabad_(&smlnum, &bignum);
  338. smlnum = sqrt(smlnum) / eps;
  339. bignum = 1. / smlnum;
  340. /* Scale A if max element outside range [SMLNUM,BIGNUM] */
  341. anrm = _starpu_dlange_("M", n, n, &a[a_offset], lda, dum);
  342. scalea = FALSE_;
  343. if (anrm > 0. && anrm < smlnum) {
  344. scalea = TRUE_;
  345. cscale = smlnum;
  346. } else if (anrm > bignum) {
  347. scalea = TRUE_;
  348. cscale = bignum;
  349. }
  350. if (scalea) {
  351. _starpu_dlascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &
  352. ierr);
  353. }
  354. /* Permute the matrix to make it more nearly triangular */
  355. /* (RWorkspace: need N) */
  356. ibal = 1;
  357. _starpu_dgebal_("P", n, &a[a_offset], lda, &ilo, &ihi, &work[ibal], &ierr);
  358. /* Reduce to upper Hessenberg form */
  359. /* (RWorkspace: need 3*N, prefer 2*N+N*NB) */
  360. itau = *n + ibal;
  361. iwrk = *n + itau;
  362. i__1 = *lwork - iwrk + 1;
  363. _starpu_dgehrd_(n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1,
  364. &ierr);
  365. if (wantvs) {
  366. /* Copy Householder vectors to VS */
  367. _starpu_dlacpy_("L", n, n, &a[a_offset], lda, &vs[vs_offset], ldvs)
  368. ;
  369. /* Generate orthogonal matrix in VS */
  370. /* (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) */
  371. i__1 = *lwork - iwrk + 1;
  372. _starpu_dorghr_(n, &ilo, &ihi, &vs[vs_offset], ldvs, &work[itau], &work[iwrk],
  373. &i__1, &ierr);
  374. }
  375. *sdim = 0;
  376. /* Perform QR iteration, accumulating Schur vectors in VS if desired */
  377. /* (RWorkspace: need N+1, prefer N+HSWORK (see comments) ) */
  378. iwrk = itau;
  379. i__1 = *lwork - iwrk + 1;
  380. _starpu_dhseqr_("S", jobvs, n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vs[
  381. vs_offset], ldvs, &work[iwrk], &i__1, &ieval);
  382. if (ieval > 0) {
  383. *info = ieval;
  384. }
  385. /* Sort eigenvalues if desired */
  386. if (wantst && *info == 0) {
  387. if (scalea) {
  388. _starpu_dlascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &wr[1], n, &
  389. ierr);
  390. _starpu_dlascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &wi[1], n, &
  391. ierr);
  392. }
  393. i__1 = *n;
  394. for (i__ = 1; i__ <= i__1; ++i__) {
  395. bwork[i__] = (*select)(&wr[i__], &wi[i__]);
  396. /* L10: */
  397. }
  398. /* Reorder eigenvalues, transform Schur vectors, and compute */
  399. /* reciprocal condition numbers */
  400. /* (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM) */
  401. /* otherwise, need N ) */
  402. /* (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM) */
  403. /* otherwise, need 0 ) */
  404. i__1 = *lwork - iwrk + 1;
  405. _starpu_dtrsen_(sense, jobvs, &bwork[1], n, &a[a_offset], lda, &vs[vs_offset],
  406. ldvs, &wr[1], &wi[1], sdim, rconde, rcondv, &work[iwrk], &
  407. i__1, &iwork[1], liwork, &icond);
  408. if (! wantsn) {
  409. /* Computing MAX */
  410. i__1 = maxwrk, i__2 = *n + (*sdim << 1) * (*n - *sdim);
  411. maxwrk = max(i__1,i__2);
  412. }
  413. if (icond == -15) {
  414. /* Not enough real workspace */
  415. *info = -16;
  416. } else if (icond == -17) {
  417. /* Not enough integer workspace */
  418. *info = -18;
  419. } else if (icond > 0) {
  420. /* DTRSEN failed to reorder or to restore standard Schur form */
  421. *info = icond + *n;
  422. }
  423. }
  424. if (wantvs) {
  425. /* Undo balancing */
  426. /* (RWorkspace: need N) */
  427. _starpu_dgebak_("P", "R", n, &ilo, &ihi, &work[ibal], n, &vs[vs_offset], ldvs,
  428. &ierr);
  429. }
  430. if (scalea) {
  431. /* Undo scaling for the Schur form of A */
  432. _starpu_dlascl_("H", &c__0, &c__0, &cscale, &anrm, n, n, &a[a_offset], lda, &
  433. ierr);
  434. i__1 = *lda + 1;
  435. _starpu_dcopy_(n, &a[a_offset], &i__1, &wr[1], &c__1);
  436. if ((wantsv || wantsb) && *info == 0) {
  437. dum[0] = *rcondv;
  438. _starpu_dlascl_("G", &c__0, &c__0, &cscale, &anrm, &c__1, &c__1, dum, &
  439. c__1, &ierr);
  440. *rcondv = dum[0];
  441. }
  442. if (cscale == smlnum) {
  443. /* If scaling back towards underflow, adjust WI if an */
  444. /* offdiagonal element of a 2-by-2 block in the Schur form */
  445. /* underflows. */
  446. if (ieval > 0) {
  447. i1 = ieval + 1;
  448. i2 = ihi - 1;
  449. i__1 = ilo - 1;
  450. _starpu_dlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[
  451. 1], n, &ierr);
  452. } else if (wantst) {
  453. i1 = 1;
  454. i2 = *n - 1;
  455. } else {
  456. i1 = ilo;
  457. i2 = ihi - 1;
  458. }
  459. inxt = i1 - 1;
  460. i__1 = i2;
  461. for (i__ = i1; i__ <= i__1; ++i__) {
  462. if (i__ < inxt) {
  463. goto L20;
  464. }
  465. if (wi[i__] == 0.) {
  466. inxt = i__ + 1;
  467. } else {
  468. if (a[i__ + 1 + i__ * a_dim1] == 0.) {
  469. wi[i__] = 0.;
  470. wi[i__ + 1] = 0.;
  471. } else if (a[i__ + 1 + i__ * a_dim1] != 0. && a[i__ + (
  472. i__ + 1) * a_dim1] == 0.) {
  473. wi[i__] = 0.;
  474. wi[i__ + 1] = 0.;
  475. if (i__ > 1) {
  476. i__2 = i__ - 1;
  477. _starpu_dswap_(&i__2, &a[i__ * a_dim1 + 1], &c__1, &a[(
  478. i__ + 1) * a_dim1 + 1], &c__1);
  479. }
  480. if (*n > i__ + 1) {
  481. i__2 = *n - i__ - 1;
  482. _starpu_dswap_(&i__2, &a[i__ + (i__ + 2) * a_dim1], lda, &
  483. a[i__ + 1 + (i__ + 2) * a_dim1], lda);
  484. }
  485. _starpu_dswap_(n, &vs[i__ * vs_dim1 + 1], &c__1, &vs[(i__ + 1)
  486. * vs_dim1 + 1], &c__1);
  487. a[i__ + (i__ + 1) * a_dim1] = a[i__ + 1 + i__ *
  488. a_dim1];
  489. a[i__ + 1 + i__ * a_dim1] = 0.;
  490. }
  491. inxt = i__ + 2;
  492. }
  493. L20:
  494. ;
  495. }
  496. }
  497. i__1 = *n - ieval;
  498. /* Computing MAX */
  499. i__3 = *n - ieval;
  500. i__2 = max(i__3,1);
  501. _starpu_dlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[ieval +
  502. 1], &i__2, &ierr);
  503. }
  504. if (wantst && *info == 0) {
  505. /* Check if reordering successful */
  506. lastsl = TRUE_;
  507. lst2sl = TRUE_;
  508. *sdim = 0;
  509. ip = 0;
  510. i__1 = *n;
  511. for (i__ = 1; i__ <= i__1; ++i__) {
  512. cursl = (*select)(&wr[i__], &wi[i__]);
  513. if (wi[i__] == 0.) {
  514. if (cursl) {
  515. ++(*sdim);
  516. }
  517. ip = 0;
  518. if (cursl && ! lastsl) {
  519. *info = *n + 2;
  520. }
  521. } else {
  522. if (ip == 1) {
  523. /* Last eigenvalue of conjugate pair */
  524. cursl = cursl || lastsl;
  525. lastsl = cursl;
  526. if (cursl) {
  527. *sdim += 2;
  528. }
  529. ip = -1;
  530. if (cursl && ! lst2sl) {
  531. *info = *n + 2;
  532. }
  533. } else {
  534. /* First eigenvalue of conjugate pair */
  535. ip = 1;
  536. }
  537. }
  538. lst2sl = lastsl;
  539. lastsl = cursl;
  540. /* L30: */
  541. }
  542. }
  543. work[1] = (doublereal) maxwrk;
  544. if (wantsv || wantsb) {
  545. /* Computing MAX */
  546. i__1 = 1, i__2 = *sdim * (*n - *sdim);
  547. iwork[1] = max(i__1,i__2);
  548. } else {
  549. iwork[1] = 1;
  550. }
  551. return 0;
  552. /* End of DGEESX */
  553. } /* _starpu_dgeesx_ */