You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

sgegv.c 39 kB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316
  1. /* f2c.h -- Standard Fortran to C header file */
  2. /** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed."
  3. - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */
  4. #ifndef F2C_INCLUDE
  5. #define F2C_INCLUDE
  6. #include <math.h>
  7. #include <stdlib.h>
  8. #include <string.h>
  9. #include <stdio.h>
  10. #include <complex.h>
  11. #ifdef complex
  12. #undef complex
  13. #endif
  14. #ifdef I
  15. #undef I
  16. #endif
  17. #if defined(_WIN64)
  18. typedef long long BLASLONG;
  19. typedef unsigned long long BLASULONG;
  20. #else
  21. typedef long BLASLONG;
  22. typedef unsigned long BLASULONG;
  23. #endif
  24. #ifdef LAPACK_ILP64
  25. typedef BLASLONG blasint;
  26. #if defined(_WIN64)
  27. #define blasabs(x) llabs(x)
  28. #else
  29. #define blasabs(x) labs(x)
  30. #endif
  31. #else
  32. typedef int blasint;
  33. #define blasabs(x) abs(x)
  34. #endif
  35. typedef blasint integer;
  36. typedef unsigned int uinteger;
  37. typedef char *address;
  38. typedef short int shortint;
  39. typedef float real;
  40. typedef double doublereal;
  41. typedef struct { real r, i; } complex;
  42. typedef struct { doublereal r, i; } doublecomplex;
  43. static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
  44. static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
  45. static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
  46. static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
  47. #define pCf(z) (*_pCf(z))
  48. #define pCd(z) (*_pCd(z))
  49. typedef int logical;
  50. typedef short int shortlogical;
  51. typedef char logical1;
  52. typedef char integer1;
  53. #define TRUE_ (1)
  54. #define FALSE_ (0)
  55. /* Extern is for use with -E */
  56. #ifndef Extern
  57. #define Extern extern
  58. #endif
  59. /* I/O stuff */
  60. typedef int flag;
  61. typedef int ftnlen;
  62. typedef int ftnint;
  63. /*external read, write*/
  64. typedef struct
  65. { flag cierr;
  66. ftnint ciunit;
  67. flag ciend;
  68. char *cifmt;
  69. ftnint cirec;
  70. } cilist;
  71. /*internal read, write*/
  72. typedef struct
  73. { flag icierr;
  74. char *iciunit;
  75. flag iciend;
  76. char *icifmt;
  77. ftnint icirlen;
  78. ftnint icirnum;
  79. } icilist;
  80. /*open*/
  81. typedef struct
  82. { flag oerr;
  83. ftnint ounit;
  84. char *ofnm;
  85. ftnlen ofnmlen;
  86. char *osta;
  87. char *oacc;
  88. char *ofm;
  89. ftnint orl;
  90. char *oblnk;
  91. } olist;
  92. /*close*/
  93. typedef struct
  94. { flag cerr;
  95. ftnint cunit;
  96. char *csta;
  97. } cllist;
  98. /*rewind, backspace, endfile*/
  99. typedef struct
  100. { flag aerr;
  101. ftnint aunit;
  102. } alist;
  103. /* inquire */
  104. typedef struct
  105. { flag inerr;
  106. ftnint inunit;
  107. char *infile;
  108. ftnlen infilen;
  109. ftnint *inex; /*parameters in standard's order*/
  110. ftnint *inopen;
  111. ftnint *innum;
  112. ftnint *innamed;
  113. char *inname;
  114. ftnlen innamlen;
  115. char *inacc;
  116. ftnlen inacclen;
  117. char *inseq;
  118. ftnlen inseqlen;
  119. char *indir;
  120. ftnlen indirlen;
  121. char *infmt;
  122. ftnlen infmtlen;
  123. char *inform;
  124. ftnint informlen;
  125. char *inunf;
  126. ftnlen inunflen;
  127. ftnint *inrecl;
  128. ftnint *innrec;
  129. char *inblank;
  130. ftnlen inblanklen;
  131. } inlist;
  132. #define VOID void
  133. union Multitype { /* for multiple entry points */
  134. integer1 g;
  135. shortint h;
  136. integer i;
  137. /* longint j; */
  138. real r;
  139. doublereal d;
  140. complex c;
  141. doublecomplex z;
  142. };
  143. typedef union Multitype Multitype;
  144. struct Vardesc { /* for Namelist */
  145. char *name;
  146. char *addr;
  147. ftnlen *dims;
  148. int type;
  149. };
  150. typedef struct Vardesc Vardesc;
  151. struct Namelist {
  152. char *name;
  153. Vardesc **vars;
  154. int nvars;
  155. };
  156. typedef struct Namelist Namelist;
  157. #define abs(x) ((x) >= 0 ? (x) : -(x))
  158. #define dabs(x) (fabs(x))
  159. #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
  160. #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
  161. #define dmin(a,b) (f2cmin(a,b))
  162. #define dmax(a,b) (f2cmax(a,b))
  163. #define bit_test(a,b) ((a) >> (b) & 1)
  164. #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
  165. #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
  166. #define abort_() { sig_die("Fortran abort routine called", 1); }
  167. #define c_abs(z) (cabsf(Cf(z)))
  168. #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
  169. #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
  170. #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
  171. #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
  172. #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
  173. #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
  174. //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
  175. #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
  176. #define d_abs(x) (fabs(*(x)))
  177. #define d_acos(x) (acos(*(x)))
  178. #define d_asin(x) (asin(*(x)))
  179. #define d_atan(x) (atan(*(x)))
  180. #define d_atn2(x, y) (atan2(*(x),*(y)))
  181. #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
  182. #define r_cnjg(R, Z) { pCf(R) = conj(Cf(Z)); }
  183. #define d_cos(x) (cos(*(x)))
  184. #define d_cosh(x) (cosh(*(x)))
  185. #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
  186. #define d_exp(x) (exp(*(x)))
  187. #define d_imag(z) (cimag(Cd(z)))
  188. #define r_imag(z) (cimag(Cf(z)))
  189. #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
  190. #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
  191. #define d_log(x) (log(*(x)))
  192. #define d_mod(x, y) (fmod(*(x), *(y)))
  193. #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
  194. #define d_nint(x) u_nint(*(x))
  195. #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
  196. #define d_sign(a,b) u_sign(*(a),*(b))
  197. #define d_sin(x) (sin(*(x)))
  198. #define d_sinh(x) (sinh(*(x)))
  199. #define d_sqrt(x) (sqrt(*(x)))
  200. #define d_tan(x) (tan(*(x)))
  201. #define d_tanh(x) (tanh(*(x)))
  202. #define i_abs(x) abs(*(x))
  203. #define i_dnnt(x) ((integer)u_nint(*(x)))
  204. #define i_len(s, n) (n)
  205. #define i_nint(x) ((integer)u_nint(*(x)))
  206. #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
  207. #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
  208. #define pow_si(B,E) spow_ui(*(B),*(E))
  209. #define pow_di(B,E) dpow_ui(*(B),*(E))
  210. #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
  211. #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
  212. #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
  213. #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
  214. #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
  215. #define sig_die(s, kill) { exit(1); }
  216. #define s_stop(s, n) {exit(0);}
  217. static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
  218. #define z_abs(z) (cabs(Cd(z)))
  219. #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
  220. #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
  221. #define myexit_() break;
  222. #define mycycle() continue;
  223. #define myceiling(w) {ceil(w)}
  224. #define myhuge(w) {HUGE_VAL}
  225. #define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
  226. /* procedure parameter types for -A and -C++ */
  227. #define F2C_proc_par_types 1
  228. #ifdef __cplusplus
  229. typedef logical (*L_fp)(...);
  230. #else
  231. typedef logical (*L_fp)();
  232. #endif
  233. static float spow_ui(float x, integer n) {
  234. float pow=1.0; unsigned long int u;
  235. if(n != 0) {
  236. if(n < 0) n = -n, x = 1/x;
  237. for(u = n; ; ) {
  238. if(u & 01) pow *= x;
  239. if(u >>= 1) x *= x;
  240. else break;
  241. }
  242. }
  243. return pow;
  244. }
  245. static double dpow_ui(double x, integer n) {
  246. double pow=1.0; unsigned long int u;
  247. if(n != 0) {
  248. if(n < 0) n = -n, x = 1/x;
  249. for(u = n; ; ) {
  250. if(u & 01) pow *= x;
  251. if(u >>= 1) x *= x;
  252. else break;
  253. }
  254. }
  255. return pow;
  256. }
  257. static _Complex float cpow_ui(_Complex float x, integer n) {
  258. _Complex float pow=1.0; unsigned long int u;
  259. if(n != 0) {
  260. if(n < 0) n = -n, x = 1/x;
  261. for(u = n; ; ) {
  262. if(u & 01) pow *= x;
  263. if(u >>= 1) x *= x;
  264. else break;
  265. }
  266. }
  267. return pow;
  268. }
  269. static _Complex double zpow_ui(_Complex double x, integer n) {
  270. _Complex double pow=1.0; unsigned long int u;
  271. if(n != 0) {
  272. if(n < 0) n = -n, x = 1/x;
  273. for(u = n; ; ) {
  274. if(u & 01) pow *= x;
  275. if(u >>= 1) x *= x;
  276. else break;
  277. }
  278. }
  279. return pow;
  280. }
  281. static integer pow_ii(integer x, integer n) {
  282. integer pow; unsigned long int u;
  283. if (n <= 0) {
  284. if (n == 0 || x == 1) pow = 1;
  285. else if (x != -1) pow = x == 0 ? 1/x : 0;
  286. else n = -n;
  287. }
  288. if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
  289. u = n;
  290. for(pow = 1; ; ) {
  291. if(u & 01) pow *= x;
  292. if(u >>= 1) x *= x;
  293. else break;
  294. }
  295. }
  296. return pow;
  297. }
  298. static integer dmaxloc_(double *w, integer s, integer e, integer *n)
  299. {
  300. double m; integer i, mi;
  301. for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
  302. if (w[i-1]>m) mi=i ,m=w[i-1];
  303. return mi-s+1;
  304. }
  305. static integer smaxloc_(float *w, integer s, integer e, integer *n)
  306. {
  307. float m; integer i, mi;
  308. for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
  309. if (w[i-1]>m) mi=i ,m=w[i-1];
  310. return mi-s+1;
  311. }
  312. static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
  313. integer n = *n_, incx = *incx_, incy = *incy_, i;
  314. _Complex float zdotc = 0.0;
  315. if (incx == 1 && incy == 1) {
  316. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  317. zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
  318. }
  319. } else {
  320. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  321. zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
  322. }
  323. }
  324. pCf(z) = zdotc;
  325. }
  326. static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
  327. integer n = *n_, incx = *incx_, incy = *incy_, i;
  328. _Complex double zdotc = 0.0;
  329. if (incx == 1 && incy == 1) {
  330. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  331. zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
  332. }
  333. } else {
  334. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  335. zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
  336. }
  337. }
  338. pCd(z) = zdotc;
  339. }
  340. static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
  341. integer n = *n_, incx = *incx_, incy = *incy_, i;
  342. _Complex float zdotc = 0.0;
  343. if (incx == 1 && incy == 1) {
  344. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  345. zdotc += Cf(&x[i]) * Cf(&y[i]);
  346. }
  347. } else {
  348. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  349. zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
  350. }
  351. }
  352. pCf(z) = zdotc;
  353. }
  354. static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
  355. integer n = *n_, incx = *incx_, incy = *incy_, i;
  356. _Complex double zdotc = 0.0;
  357. if (incx == 1 && incy == 1) {
  358. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  359. zdotc += Cd(&x[i]) * Cd(&y[i]);
  360. }
  361. } else {
  362. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  363. zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
  364. }
  365. }
  366. pCd(z) = zdotc;
  367. }
  368. #endif
  369. /* -- translated by f2c (version 20000121).
  370. You must link the resulting object file with the libraries:
  371. -lf2c -lm (in that order)
  372. */
  373. /* Table of constant values */
  374. static integer c__1 = 1;
  375. static integer c_n1 = -1;
  376. static real c_b27 = 1.f;
  377. static real c_b38 = 0.f;
  378. /* > \brief <b> SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE mat
  379. rices</b> */
  380. /* =========== DOCUMENTATION =========== */
  381. /* Online html documentation available at */
  382. /* http://www.netlib.org/lapack/explore-html/ */
  383. /* > \htmlonly */
  384. /* > Download SGEGV + dependencies */
  385. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgegv.f
  386. "> */
  387. /* > [TGZ]</a> */
  388. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgegv.f
  389. "> */
  390. /* > [ZIP]</a> */
  391. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgegv.f
  392. "> */
  393. /* > [TXT]</a> */
  394. /* > \endhtmlonly */
  395. /* Definition: */
  396. /* =========== */
  397. /* SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, */
  398. /* BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) */
  399. /* CHARACTER JOBVL, JOBVR */
  400. /* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N */
  401. /* REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ), */
  402. /* $ B( LDB, * ), BETA( * ), VL( LDVL, * ), */
  403. /* $ VR( LDVR, * ), WORK( * ) */
  404. /* > \par Purpose: */
  405. /* ============= */
  406. /* > */
  407. /* > \verbatim */
  408. /* > */
  409. /* > This routine is deprecated and has been replaced by routine SGGEV. */
  410. /* > */
  411. /* > SGEGV computes the eigenvalues and, optionally, the left and/or right */
  412. /* > eigenvectors of a real matrix pair (A,B). */
  413. /* > Given two square matrices A and B, */
  414. /* > the generalized nonsymmetric eigenvalue problem (GNEP) is to find the */
  415. /* > eigenvalues lambda and corresponding (non-zero) eigenvectors x such */
  416. /* > that */
  417. /* > */
  418. /* > A*x = lambda*B*x. */
  419. /* > */
  420. /* > An alternate form is to find the eigenvalues mu and corresponding */
  421. /* > eigenvectors y such that */
  422. /* > */
  423. /* > mu*A*y = B*y. */
  424. /* > */
  425. /* > These two forms are equivalent with mu = 1/lambda and x = y if */
  426. /* > neither lambda nor mu is zero. In order to deal with the case that */
  427. /* > lambda or mu is zero or small, two values alpha and beta are returned */
  428. /* > for each eigenvalue, such that lambda = alpha/beta and */
  429. /* > mu = beta/alpha. */
  430. /* > */
  431. /* > The vectors x and y in the above equations are right eigenvectors of */
  432. /* > the matrix pair (A,B). Vectors u and v satisfying */
  433. /* > */
  434. /* > u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B */
  435. /* > */
  436. /* > are left eigenvectors of (A,B). */
  437. /* > */
  438. /* > Note: this routine performs "full balancing" on A and B */
  439. /* > \endverbatim */
  440. /* Arguments: */
  441. /* ========== */
  442. /* > \param[in] JOBVL */
  443. /* > \verbatim */
  444. /* > JOBVL is CHARACTER*1 */
  445. /* > = 'N': do not compute the left generalized eigenvectors; */
  446. /* > = 'V': compute the left generalized eigenvectors (returned */
  447. /* > in VL). */
  448. /* > \endverbatim */
  449. /* > */
  450. /* > \param[in] JOBVR */
  451. /* > \verbatim */
  452. /* > JOBVR is CHARACTER*1 */
  453. /* > = 'N': do not compute the right generalized eigenvectors; */
  454. /* > = 'V': compute the right generalized eigenvectors (returned */
  455. /* > in VR). */
  456. /* > \endverbatim */
  457. /* > */
  458. /* > \param[in] N */
  459. /* > \verbatim */
  460. /* > N is INTEGER */
  461. /* > The order of the matrices A, B, VL, and VR. N >= 0. */
  462. /* > \endverbatim */
  463. /* > */
  464. /* > \param[in,out] A */
  465. /* > \verbatim */
  466. /* > A is REAL array, dimension (LDA, N) */
  467. /* > On entry, the matrix A. */
  468. /* > If JOBVL = 'V' or JOBVR = 'V', then on exit A */
  469. /* > contains the real Schur form of A from the generalized Schur */
  470. /* > factorization of the pair (A,B) after balancing. */
  471. /* > If no eigenvectors were computed, then only the diagonal */
  472. /* > blocks from the Schur form will be correct. See SGGHRD and */
  473. /* > SHGEQZ for details. */
  474. /* > \endverbatim */
  475. /* > */
  476. /* > \param[in] LDA */
  477. /* > \verbatim */
  478. /* > LDA is INTEGER */
  479. /* > The leading dimension of A. LDA >= f2cmax(1,N). */
  480. /* > \endverbatim */
  481. /* > */
  482. /* > \param[in,out] B */
  483. /* > \verbatim */
  484. /* > B is REAL array, dimension (LDB, N) */
  485. /* > On entry, the matrix B. */
  486. /* > If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the */
  487. /* > upper triangular matrix obtained from B in the generalized */
  488. /* > Schur factorization of the pair (A,B) after balancing. */
  489. /* > If no eigenvectors were computed, then only those elements of */
  490. /* > B corresponding to the diagonal blocks from the Schur form of */
  491. /* > A will be correct. See SGGHRD and SHGEQZ for details. */
  492. /* > \endverbatim */
  493. /* > */
  494. /* > \param[in] LDB */
  495. /* > \verbatim */
  496. /* > LDB is INTEGER */
  497. /* > The leading dimension of B. LDB >= f2cmax(1,N). */
  498. /* > \endverbatim */
  499. /* > */
  500. /* > \param[out] ALPHAR */
  501. /* > \verbatim */
  502. /* > ALPHAR is REAL array, dimension (N) */
  503. /* > The real parts of each scalar alpha defining an eigenvalue of */
  504. /* > GNEP. */
  505. /* > \endverbatim */
  506. /* > */
  507. /* > \param[out] ALPHAI */
  508. /* > \verbatim */
  509. /* > ALPHAI is REAL array, dimension (N) */
  510. /* > The imaginary parts of each scalar alpha defining an */
  511. /* > eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th */
  512. /* > eigenvalue is real; if positive, then the j-th and */
  513. /* > (j+1)-st eigenvalues are a complex conjugate pair, with */
  514. /* > ALPHAI(j+1) = -ALPHAI(j). */
  515. /* > \endverbatim */
  516. /* > */
  517. /* > \param[out] BETA */
  518. /* > \verbatim */
  519. /* > BETA is REAL array, dimension (N) */
  520. /* > The scalars beta that define the eigenvalues of GNEP. */
  521. /* > */
  522. /* > Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and */
  523. /* > beta = BETA(j) represent the j-th eigenvalue of the matrix */
  524. /* > pair (A,B), in one of the forms lambda = alpha/beta or */
  525. /* > mu = beta/alpha. Since either lambda or mu may overflow, */
  526. /* > they should not, in general, be computed. */
  527. /* > \endverbatim */
  528. /* > */
  529. /* > \param[out] VL */
  530. /* > \verbatim */
  531. /* > VL is REAL array, dimension (LDVL,N) */
  532. /* > If JOBVL = 'V', the left eigenvectors u(j) are stored */
  533. /* > in the columns of VL, in the same order as their eigenvalues. */
  534. /* > If the j-th eigenvalue is real, then u(j) = VL(:,j). */
  535. /* > If the j-th and (j+1)-st eigenvalues form a complex conjugate */
  536. /* > pair, then */
  537. /* > u(j) = VL(:,j) + i*VL(:,j+1) */
  538. /* > and */
  539. /* > u(j+1) = VL(:,j) - i*VL(:,j+1). */
  540. /* > */
  541. /* > Each eigenvector is scaled so that its largest component has */
  542. /* > abs(real part) + abs(imag. part) = 1, except for eigenvectors */
  543. /* > corresponding to an eigenvalue with alpha = beta = 0, which */
  544. /* > are set to zero. */
  545. /* > Not referenced if JOBVL = 'N'. */
  546. /* > \endverbatim */
  547. /* > */
  548. /* > \param[in] LDVL */
  549. /* > \verbatim */
  550. /* > LDVL is INTEGER */
  551. /* > The leading dimension of the matrix VL. LDVL >= 1, and */
  552. /* > if JOBVL = 'V', LDVL >= N. */
  553. /* > \endverbatim */
  554. /* > */
  555. /* > \param[out] VR */
  556. /* > \verbatim */
  557. /* > VR is REAL array, dimension (LDVR,N) */
  558. /* > If JOBVR = 'V', the right eigenvectors x(j) are stored */
  559. /* > in the columns of VR, in the same order as their eigenvalues. */
  560. /* > If the j-th eigenvalue is real, then x(j) = VR(:,j). */
  561. /* > If the j-th and (j+1)-st eigenvalues form a complex conjugate */
  562. /* > pair, then */
  563. /* > x(j) = VR(:,j) + i*VR(:,j+1) */
  564. /* > and */
  565. /* > x(j+1) = VR(:,j) - i*VR(:,j+1). */
  566. /* > */
  567. /* > Each eigenvector is scaled so that its largest component has */
  568. /* > abs(real part) + abs(imag. part) = 1, except for eigenvalues */
  569. /* > corresponding to an eigenvalue with alpha = beta = 0, which */
  570. /* > are set to zero. */
  571. /* > Not referenced if JOBVR = 'N'. */
  572. /* > \endverbatim */
  573. /* > */
  574. /* > \param[in] LDVR */
  575. /* > \verbatim */
  576. /* > LDVR is INTEGER */
  577. /* > The leading dimension of the matrix VR. LDVR >= 1, and */
  578. /* > if JOBVR = 'V', LDVR >= N. */
  579. /* > \endverbatim */
  580. /* > */
  581. /* > \param[out] WORK */
  582. /* > \verbatim */
  583. /* > WORK is REAL array, dimension (MAX(1,LWORK)) */
  584. /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  585. /* > \endverbatim */
  586. /* > */
  587. /* > \param[in] LWORK */
  588. /* > \verbatim */
  589. /* > LWORK is INTEGER */
  590. /* > The dimension of the array WORK. LWORK >= f2cmax(1,8*N). */
  591. /* > For good performance, LWORK must generally be larger. */
  592. /* > To compute the optimal value of LWORK, call ILAENV to get */
  593. /* > blocksizes (for SGEQRF, SORMQR, and SORGQR.) Then compute: */
  594. /* > NB -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR; */
  595. /* > The optimal LWORK is: */
  596. /* > 2*N + MAX( 6*N, N*(NB+1) ). */
  597. /* > */
  598. /* > If LWORK = -1, then a workspace query is assumed; the routine */
  599. /* > only calculates the optimal size of the WORK array, returns */
  600. /* > this value as the first entry of the WORK array, and no error */
  601. /* > message related to LWORK is issued by XERBLA. */
  602. /* > \endverbatim */
  603. /* > */
  604. /* > \param[out] INFO */
  605. /* > \verbatim */
  606. /* > INFO is INTEGER */
  607. /* > = 0: successful exit */
  608. /* > < 0: if INFO = -i, the i-th argument had an illegal value. */
  609. /* > = 1,...,N: */
  610. /* > The QZ iteration failed. No eigenvectors have been */
  611. /* > calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) */
  612. /* > should be correct for j=INFO+1,...,N. */
  613. /* > > N: errors that usually indicate LAPACK problems: */
  614. /* > =N+1: error return from SGGBAL */
  615. /* > =N+2: error return from SGEQRF */
  616. /* > =N+3: error return from SORMQR */
  617. /* > =N+4: error return from SORGQR */
  618. /* > =N+5: error return from SGGHRD */
  619. /* > =N+6: error return from SHGEQZ (other than failed */
  620. /* > iteration) */
  621. /* > =N+7: error return from STGEVC */
  622. /* > =N+8: error return from SGGBAK (computing VL) */
  623. /* > =N+9: error return from SGGBAK (computing VR) */
  624. /* > =N+10: error return from SLASCL (various calls) */
  625. /* > \endverbatim */
  626. /* Authors: */
  627. /* ======== */
  628. /* > \author Univ. of Tennessee */
  629. /* > \author Univ. of California Berkeley */
  630. /* > \author Univ. of Colorado Denver */
  631. /* > \author NAG Ltd. */
  632. /* > \date December 2016 */
  633. /* > \ingroup realGEeigen */
  634. /* > \par Further Details: */
  635. /* ===================== */
  636. /* > */
  637. /* > \verbatim */
  638. /* > */
  639. /* > Balancing */
  640. /* > --------- */
  641. /* > */
  642. /* > This driver calls SGGBAL to both permute and scale rows and columns */
  643. /* > of A and B. The permutations PL and PR are chosen so that PL*A*PR */
  644. /* > and PL*B*R will be upper triangular except for the diagonal blocks */
  645. /* > A(i:j,i:j) and B(i:j,i:j), with i and j as close together as */
  646. /* > possible. The diagonal scaling matrices DL and DR are chosen so */
  647. /* > that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to */
  648. /* > one (except for the elements that start out zero.) */
  649. /* > */
  650. /* > After the eigenvalues and eigenvectors of the balanced matrices */
  651. /* > have been computed, SGGBAK transforms the eigenvectors back to what */
  652. /* > they would have been (in perfect arithmetic) if they had not been */
  653. /* > balanced. */
  654. /* > */
  655. /* > Contents of A and B on Exit */
  656. /* > -------- -- - --- - -- ---- */
  657. /* > */
  658. /* > If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or */
  659. /* > both), then on exit the arrays A and B will contain the real Schur */
  660. /* > form[*] of the "balanced" versions of A and B. If no eigenvectors */
  661. /* > are computed, then only the diagonal blocks will be correct. */
  662. /* > */
  663. /* > [*] See SHGEQZ, SGEGS, or read the book "Matrix Computations", */
  664. /* > by Golub & van Loan, pub. by Johns Hopkins U. Press. */
  665. /* > \endverbatim */
  666. /* > */
  667. /* ===================================================================== */
  668. /* Subroutine */ int sgegv_(char *jobvl, char *jobvr, integer *n, real *a,
  669. integer *lda, real *b, integer *ldb, real *alphar, real *alphai, real
  670. *beta, real *vl, integer *ldvl, real *vr, integer *ldvr, real *work,
  671. integer *lwork, integer *info)
  672. {
  673. /* System generated locals */
  674. integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
  675. vr_offset, i__1, i__2;
  676. real r__1, r__2, r__3, r__4;
  677. /* Local variables */
  678. real absb, anrm, bnrm;
  679. integer itau;
  680. real temp;
  681. logical ilvl, ilvr;
  682. integer lopt;
  683. real anrm1, anrm2, bnrm1, bnrm2, absai, scale, absar, sbeta;
  684. extern logical lsame_(char *, char *);
  685. integer ileft, iinfo, icols, iwork, irows, jc, nb, in, jr;
  686. real salfai;
  687. extern /* Subroutine */ int sggbak_(char *, char *, integer *, integer *,
  688. integer *, real *, real *, integer *, real *, integer *, integer *
  689. ), sggbal_(char *, integer *, real *, integer *,
  690. real *, integer *, integer *, integer *, real *, real *, real *,
  691. integer *);
  692. real salfar;
  693. extern real slamch_(char *), slange_(char *, integer *, integer *,
  694. real *, integer *, real *);
  695. real safmin;
  696. extern /* Subroutine */ int sgghrd_(char *, char *, integer *, integer *,
  697. integer *, real *, integer *, real *, integer *, real *, integer *
  698. , real *, integer *, integer *);
  699. real safmax;
  700. char chtemp[1];
  701. logical ldumma[1];
  702. extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
  703. real *, integer *, integer *, real *, integer *, integer *), xerbla_(char *, integer *);
  704. extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  705. integer *, integer *, ftnlen, ftnlen);
  706. integer ijobvl, iright;
  707. logical ilimit;
  708. extern /* Subroutine */ int sgeqrf_(integer *, integer *, real *, integer
  709. *, real *, real *, integer *, integer *);
  710. integer ijobvr;
  711. extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *,
  712. integer *, real *, integer *), slaset_(char *, integer *,
  713. integer *, real *, real *, real *, integer *), stgevc_(
  714. char *, char *, logical *, integer *, real *, integer *, real *,
  715. integer *, real *, integer *, real *, integer *, integer *,
  716. integer *, real *, integer *);
  717. real onepls;
  718. integer lwkmin, nb1, nb2, nb3;
  719. extern /* Subroutine */ int shgeqz_(char *, char *, char *, integer *,
  720. integer *, integer *, real *, integer *, real *, integer *, real *
  721. , real *, real *, real *, integer *, real *, integer *, real *,
  722. integer *, integer *), sorgqr_(integer *,
  723. integer *, integer *, real *, integer *, real *, real *, integer *
  724. , integer *);
  725. integer lwkopt;
  726. logical lquery;
  727. extern /* Subroutine */ int sormqr_(char *, char *, integer *, integer *,
  728. integer *, real *, integer *, real *, real *, integer *, real *,
  729. integer *, integer *);
  730. integer ihi, ilo;
  731. real eps;
  732. logical ilv;
  733. /* -- LAPACK driver routine (version 3.7.0) -- */
  734. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  735. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  736. /* December 2016 */
  737. /* ===================================================================== */
  738. /* Decode the input arguments */
  739. /* Parameter adjustments */
  740. a_dim1 = *lda;
  741. a_offset = 1 + a_dim1 * 1;
  742. a -= a_offset;
  743. b_dim1 = *ldb;
  744. b_offset = 1 + b_dim1 * 1;
  745. b -= b_offset;
  746. --alphar;
  747. --alphai;
  748. --beta;
  749. vl_dim1 = *ldvl;
  750. vl_offset = 1 + vl_dim1 * 1;
  751. vl -= vl_offset;
  752. vr_dim1 = *ldvr;
  753. vr_offset = 1 + vr_dim1 * 1;
  754. vr -= vr_offset;
  755. --work;
  756. /* Function Body */
  757. if (lsame_(jobvl, "N")) {
  758. ijobvl = 1;
  759. ilvl = FALSE_;
  760. } else if (lsame_(jobvl, "V")) {
  761. ijobvl = 2;
  762. ilvl = TRUE_;
  763. } else {
  764. ijobvl = -1;
  765. ilvl = FALSE_;
  766. }
  767. if (lsame_(jobvr, "N")) {
  768. ijobvr = 1;
  769. ilvr = FALSE_;
  770. } else if (lsame_(jobvr, "V")) {
  771. ijobvr = 2;
  772. ilvr = TRUE_;
  773. } else {
  774. ijobvr = -1;
  775. ilvr = FALSE_;
  776. }
  777. ilv = ilvl || ilvr;
  778. /* Test the input arguments */
  779. /* Computing MAX */
  780. i__1 = *n << 3;
  781. lwkmin = f2cmax(i__1,1);
  782. lwkopt = lwkmin;
  783. work[1] = (real) lwkopt;
  784. lquery = *lwork == -1;
  785. *info = 0;
  786. if (ijobvl <= 0) {
  787. *info = -1;
  788. } else if (ijobvr <= 0) {
  789. *info = -2;
  790. } else if (*n < 0) {
  791. *info = -3;
  792. } else if (*lda < f2cmax(1,*n)) {
  793. *info = -5;
  794. } else if (*ldb < f2cmax(1,*n)) {
  795. *info = -7;
  796. } else if (*ldvl < 1 || ilvl && *ldvl < *n) {
  797. *info = -12;
  798. } else if (*ldvr < 1 || ilvr && *ldvr < *n) {
  799. *info = -14;
  800. } else if (*lwork < lwkmin && ! lquery) {
  801. *info = -16;
  802. }
  803. if (*info == 0) {
  804. nb1 = ilaenv_(&c__1, "SGEQRF", " ", n, n, &c_n1, &c_n1, (ftnlen)6, (
  805. ftnlen)1);
  806. nb2 = ilaenv_(&c__1, "SORMQR", " ", n, n, n, &c_n1, (ftnlen)6, (
  807. ftnlen)1);
  808. nb3 = ilaenv_(&c__1, "SORGQR", " ", n, n, n, &c_n1, (ftnlen)6, (
  809. ftnlen)1);
  810. /* Computing MAX */
  811. i__1 = f2cmax(nb1,nb2);
  812. nb = f2cmax(i__1,nb3);
  813. /* Computing MAX */
  814. i__1 = *n * 6, i__2 = *n * (nb + 1);
  815. lopt = (*n << 1) + f2cmax(i__1,i__2);
  816. work[1] = (real) lopt;
  817. }
  818. if (*info != 0) {
  819. i__1 = -(*info);
  820. xerbla_("SGEGV ", &i__1);
  821. return 0;
  822. } else if (lquery) {
  823. return 0;
  824. }
  825. /* Quick return if possible */
  826. if (*n == 0) {
  827. return 0;
  828. }
  829. /* Get machine constants */
  830. eps = slamch_("E") * slamch_("B");
  831. safmin = slamch_("S");
  832. safmin += safmin;
  833. safmax = 1.f / safmin;
  834. onepls = eps * 4 + 1.f;
  835. /* Scale A */
  836. anrm = slange_("M", n, n, &a[a_offset], lda, &work[1]);
  837. anrm1 = anrm;
  838. anrm2 = 1.f;
  839. if (anrm < 1.f) {
  840. if (safmax * anrm < 1.f) {
  841. anrm1 = safmin;
  842. anrm2 = safmax * anrm;
  843. }
  844. }
  845. if (anrm > 0.f) {
  846. slascl_("G", &c_n1, &c_n1, &anrm, &c_b27, n, n, &a[a_offset], lda, &
  847. iinfo);
  848. if (iinfo != 0) {
  849. *info = *n + 10;
  850. return 0;
  851. }
  852. }
  853. /* Scale B */
  854. bnrm = slange_("M", n, n, &b[b_offset], ldb, &work[1]);
  855. bnrm1 = bnrm;
  856. bnrm2 = 1.f;
  857. if (bnrm < 1.f) {
  858. if (safmax * bnrm < 1.f) {
  859. bnrm1 = safmin;
  860. bnrm2 = safmax * bnrm;
  861. }
  862. }
  863. if (bnrm > 0.f) {
  864. slascl_("G", &c_n1, &c_n1, &bnrm, &c_b27, n, n, &b[b_offset], ldb, &
  865. iinfo);
  866. if (iinfo != 0) {
  867. *info = *n + 10;
  868. return 0;
  869. }
  870. }
  871. /* Permute the matrix to make it more nearly triangular */
  872. /* Workspace layout: (8*N words -- "work" requires 6*N words) */
  873. /* left_permutation, right_permutation, work... */
  874. ileft = 1;
  875. iright = *n + 1;
  876. iwork = iright + *n;
  877. sggbal_("P", n, &a[a_offset], lda, &b[b_offset], ldb, &ilo, &ihi, &work[
  878. ileft], &work[iright], &work[iwork], &iinfo);
  879. if (iinfo != 0) {
  880. *info = *n + 1;
  881. goto L120;
  882. }
  883. /* Reduce B to triangular form, and initialize VL and/or VR */
  884. /* Workspace layout: ("work..." must have at least N words) */
  885. /* left_permutation, right_permutation, tau, work... */
  886. irows = ihi + 1 - ilo;
  887. if (ilv) {
  888. icols = *n + 1 - ilo;
  889. } else {
  890. icols = irows;
  891. }
  892. itau = iwork;
  893. iwork = itau + irows;
  894. i__1 = *lwork + 1 - iwork;
  895. sgeqrf_(&irows, &icols, &b[ilo + ilo * b_dim1], ldb, &work[itau], &work[
  896. iwork], &i__1, &iinfo);
  897. if (iinfo >= 0) {
  898. /* Computing MAX */
  899. i__1 = lwkopt, i__2 = (integer) work[iwork] + iwork - 1;
  900. lwkopt = f2cmax(i__1,i__2);
  901. }
  902. if (iinfo != 0) {
  903. *info = *n + 2;
  904. goto L120;
  905. }
  906. i__1 = *lwork + 1 - iwork;
  907. sormqr_("L", "T", &irows, &icols, &irows, &b[ilo + ilo * b_dim1], ldb, &
  908. work[itau], &a[ilo + ilo * a_dim1], lda, &work[iwork], &i__1, &
  909. iinfo);
  910. if (iinfo >= 0) {
  911. /* Computing MAX */
  912. i__1 = lwkopt, i__2 = (integer) work[iwork] + iwork - 1;
  913. lwkopt = f2cmax(i__1,i__2);
  914. }
  915. if (iinfo != 0) {
  916. *info = *n + 3;
  917. goto L120;
  918. }
  919. if (ilvl) {
  920. slaset_("Full", n, n, &c_b38, &c_b27, &vl[vl_offset], ldvl)
  921. ;
  922. i__1 = irows - 1;
  923. i__2 = irows - 1;
  924. slacpy_("L", &i__1, &i__2, &b[ilo + 1 + ilo * b_dim1], ldb, &vl[ilo +
  925. 1 + ilo * vl_dim1], ldvl);
  926. i__1 = *lwork + 1 - iwork;
  927. sorgqr_(&irows, &irows, &irows, &vl[ilo + ilo * vl_dim1], ldvl, &work[
  928. itau], &work[iwork], &i__1, &iinfo);
  929. if (iinfo >= 0) {
  930. /* Computing MAX */
  931. i__1 = lwkopt, i__2 = (integer) work[iwork] + iwork - 1;
  932. lwkopt = f2cmax(i__1,i__2);
  933. }
  934. if (iinfo != 0) {
  935. *info = *n + 4;
  936. goto L120;
  937. }
  938. }
  939. if (ilvr) {
  940. slaset_("Full", n, n, &c_b38, &c_b27, &vr[vr_offset], ldvr)
  941. ;
  942. }
  943. /* Reduce to generalized Hessenberg form */
  944. if (ilv) {
  945. /* Eigenvectors requested -- work on whole matrix. */
  946. sgghrd_(jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset],
  947. ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &iinfo);
  948. } else {
  949. sgghrd_("N", "N", &irows, &c__1, &irows, &a[ilo + ilo * a_dim1], lda,
  950. &b[ilo + ilo * b_dim1], ldb, &vl[vl_offset], ldvl, &vr[
  951. vr_offset], ldvr, &iinfo);
  952. }
  953. if (iinfo != 0) {
  954. *info = *n + 5;
  955. goto L120;
  956. }
  957. /* Perform QZ algorithm */
  958. /* Workspace layout: ("work..." must have at least 1 word) */
  959. /* left_permutation, right_permutation, work... */
  960. iwork = itau;
  961. if (ilv) {
  962. *(unsigned char *)chtemp = 'S';
  963. } else {
  964. *(unsigned char *)chtemp = 'E';
  965. }
  966. i__1 = *lwork + 1 - iwork;
  967. shgeqz_(chtemp, jobvl, jobvr, n, &ilo, &ihi, &a[a_offset], lda, &b[
  968. b_offset], ldb, &alphar[1], &alphai[1], &beta[1], &vl[vl_offset],
  969. ldvl, &vr[vr_offset], ldvr, &work[iwork], &i__1, &iinfo);
  970. if (iinfo >= 0) {
  971. /* Computing MAX */
  972. i__1 = lwkopt, i__2 = (integer) work[iwork] + iwork - 1;
  973. lwkopt = f2cmax(i__1,i__2);
  974. }
  975. if (iinfo != 0) {
  976. if (iinfo > 0 && iinfo <= *n) {
  977. *info = iinfo;
  978. } else if (iinfo > *n && iinfo <= *n << 1) {
  979. *info = iinfo - *n;
  980. } else {
  981. *info = *n + 6;
  982. }
  983. goto L120;
  984. }
  985. if (ilv) {
  986. /* Compute Eigenvectors (STGEVC requires 6*N words of workspace) */
  987. if (ilvl) {
  988. if (ilvr) {
  989. *(unsigned char *)chtemp = 'B';
  990. } else {
  991. *(unsigned char *)chtemp = 'L';
  992. }
  993. } else {
  994. *(unsigned char *)chtemp = 'R';
  995. }
  996. stgevc_(chtemp, "B", ldumma, n, &a[a_offset], lda, &b[b_offset], ldb,
  997. &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, n, &in, &work[
  998. iwork], &iinfo);
  999. if (iinfo != 0) {
  1000. *info = *n + 7;
  1001. goto L120;
  1002. }
  1003. /* Undo balancing on VL and VR, rescale */
  1004. if (ilvl) {
  1005. sggbak_("P", "L", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
  1006. vl[vl_offset], ldvl, &iinfo);
  1007. if (iinfo != 0) {
  1008. *info = *n + 8;
  1009. goto L120;
  1010. }
  1011. i__1 = *n;
  1012. for (jc = 1; jc <= i__1; ++jc) {
  1013. if (alphai[jc] < 0.f) {
  1014. goto L50;
  1015. }
  1016. temp = 0.f;
  1017. if (alphai[jc] == 0.f) {
  1018. i__2 = *n;
  1019. for (jr = 1; jr <= i__2; ++jr) {
  1020. /* Computing MAX */
  1021. r__2 = temp, r__3 = (r__1 = vl[jr + jc * vl_dim1],
  1022. abs(r__1));
  1023. temp = f2cmax(r__2,r__3);
  1024. /* L10: */
  1025. }
  1026. } else {
  1027. i__2 = *n;
  1028. for (jr = 1; jr <= i__2; ++jr) {
  1029. /* Computing MAX */
  1030. r__3 = temp, r__4 = (r__1 = vl[jr + jc * vl_dim1],
  1031. abs(r__1)) + (r__2 = vl[jr + (jc + 1) *
  1032. vl_dim1], abs(r__2));
  1033. temp = f2cmax(r__3,r__4);
  1034. /* L20: */
  1035. }
  1036. }
  1037. if (temp < safmin) {
  1038. goto L50;
  1039. }
  1040. temp = 1.f / temp;
  1041. if (alphai[jc] == 0.f) {
  1042. i__2 = *n;
  1043. for (jr = 1; jr <= i__2; ++jr) {
  1044. vl[jr + jc * vl_dim1] *= temp;
  1045. /* L30: */
  1046. }
  1047. } else {
  1048. i__2 = *n;
  1049. for (jr = 1; jr <= i__2; ++jr) {
  1050. vl[jr + jc * vl_dim1] *= temp;
  1051. vl[jr + (jc + 1) * vl_dim1] *= temp;
  1052. /* L40: */
  1053. }
  1054. }
  1055. L50:
  1056. ;
  1057. }
  1058. }
  1059. if (ilvr) {
  1060. sggbak_("P", "R", n, &ilo, &ihi, &work[ileft], &work[iright], n, &
  1061. vr[vr_offset], ldvr, &iinfo);
  1062. if (iinfo != 0) {
  1063. *info = *n + 9;
  1064. goto L120;
  1065. }
  1066. i__1 = *n;
  1067. for (jc = 1; jc <= i__1; ++jc) {
  1068. if (alphai[jc] < 0.f) {
  1069. goto L100;
  1070. }
  1071. temp = 0.f;
  1072. if (alphai[jc] == 0.f) {
  1073. i__2 = *n;
  1074. for (jr = 1; jr <= i__2; ++jr) {
  1075. /* Computing MAX */
  1076. r__2 = temp, r__3 = (r__1 = vr[jr + jc * vr_dim1],
  1077. abs(r__1));
  1078. temp = f2cmax(r__2,r__3);
  1079. /* L60: */
  1080. }
  1081. } else {
  1082. i__2 = *n;
  1083. for (jr = 1; jr <= i__2; ++jr) {
  1084. /* Computing MAX */
  1085. r__3 = temp, r__4 = (r__1 = vr[jr + jc * vr_dim1],
  1086. abs(r__1)) + (r__2 = vr[jr + (jc + 1) *
  1087. vr_dim1], abs(r__2));
  1088. temp = f2cmax(r__3,r__4);
  1089. /* L70: */
  1090. }
  1091. }
  1092. if (temp < safmin) {
  1093. goto L100;
  1094. }
  1095. temp = 1.f / temp;
  1096. if (alphai[jc] == 0.f) {
  1097. i__2 = *n;
  1098. for (jr = 1; jr <= i__2; ++jr) {
  1099. vr[jr + jc * vr_dim1] *= temp;
  1100. /* L80: */
  1101. }
  1102. } else {
  1103. i__2 = *n;
  1104. for (jr = 1; jr <= i__2; ++jr) {
  1105. vr[jr + jc * vr_dim1] *= temp;
  1106. vr[jr + (jc + 1) * vr_dim1] *= temp;
  1107. /* L90: */
  1108. }
  1109. }
  1110. L100:
  1111. ;
  1112. }
  1113. }
  1114. /* End of eigenvector calculation */
  1115. }
  1116. /* Undo scaling in alpha, beta */
  1117. /* Note: this does not give the alpha and beta for the unscaled */
  1118. /* problem. */
  1119. /* Un-scaling is limited to avoid underflow in alpha and beta */
  1120. /* if they are significant. */
  1121. i__1 = *n;
  1122. for (jc = 1; jc <= i__1; ++jc) {
  1123. absar = (r__1 = alphar[jc], abs(r__1));
  1124. absai = (r__1 = alphai[jc], abs(r__1));
  1125. absb = (r__1 = beta[jc], abs(r__1));
  1126. salfar = anrm * alphar[jc];
  1127. salfai = anrm * alphai[jc];
  1128. sbeta = bnrm * beta[jc];
  1129. ilimit = FALSE_;
  1130. scale = 1.f;
  1131. /* Check for significant underflow in ALPHAI */
  1132. /* Computing MAX */
  1133. r__1 = safmin, r__2 = eps * absar, r__1 = f2cmax(r__1,r__2), r__2 = eps *
  1134. absb;
  1135. if (abs(salfai) < safmin && absai >= f2cmax(r__1,r__2)) {
  1136. ilimit = TRUE_;
  1137. /* Computing MAX */
  1138. r__1 = onepls * safmin, r__2 = anrm2 * absai;
  1139. scale = onepls * safmin / anrm1 / f2cmax(r__1,r__2);
  1140. } else if (salfai == 0.f) {
  1141. /* If insignificant underflow in ALPHAI, then make the */
  1142. /* conjugate eigenvalue real. */
  1143. if (alphai[jc] < 0.f && jc > 1) {
  1144. alphai[jc - 1] = 0.f;
  1145. } else if (alphai[jc] > 0.f && jc < *n) {
  1146. alphai[jc + 1] = 0.f;
  1147. }
  1148. }
  1149. /* Check for significant underflow in ALPHAR */
  1150. /* Computing MAX */
  1151. r__1 = safmin, r__2 = eps * absai, r__1 = f2cmax(r__1,r__2), r__2 = eps *
  1152. absb;
  1153. if (abs(salfar) < safmin && absar >= f2cmax(r__1,r__2)) {
  1154. ilimit = TRUE_;
  1155. /* Computing MAX */
  1156. /* Computing MAX */
  1157. r__3 = onepls * safmin, r__4 = anrm2 * absar;
  1158. r__1 = scale, r__2 = onepls * safmin / anrm1 / f2cmax(r__3,r__4);
  1159. scale = f2cmax(r__1,r__2);
  1160. }
  1161. /* Check for significant underflow in BETA */
  1162. /* Computing MAX */
  1163. r__1 = safmin, r__2 = eps * absar, r__1 = f2cmax(r__1,r__2), r__2 = eps *
  1164. absai;
  1165. if (abs(sbeta) < safmin && absb >= f2cmax(r__1,r__2)) {
  1166. ilimit = TRUE_;
  1167. /* Computing MAX */
  1168. /* Computing MAX */
  1169. r__3 = onepls * safmin, r__4 = bnrm2 * absb;
  1170. r__1 = scale, r__2 = onepls * safmin / bnrm1 / f2cmax(r__3,r__4);
  1171. scale = f2cmax(r__1,r__2);
  1172. }
  1173. /* Check for possible overflow when limiting scaling */
  1174. if (ilimit) {
  1175. /* Computing MAX */
  1176. r__1 = abs(salfar), r__2 = abs(salfai), r__1 = f2cmax(r__1,r__2),
  1177. r__2 = abs(sbeta);
  1178. temp = scale * safmin * f2cmax(r__1,r__2);
  1179. if (temp > 1.f) {
  1180. scale /= temp;
  1181. }
  1182. if (scale < 1.f) {
  1183. ilimit = FALSE_;
  1184. }
  1185. }
  1186. /* Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary. */
  1187. if (ilimit) {
  1188. salfar = scale * alphar[jc] * anrm;
  1189. salfai = scale * alphai[jc] * anrm;
  1190. sbeta = scale * beta[jc] * bnrm;
  1191. }
  1192. alphar[jc] = salfar;
  1193. alphai[jc] = salfai;
  1194. beta[jc] = sbeta;
  1195. /* L110: */
  1196. }
  1197. L120:
  1198. work[1] = (real) lwkopt;
  1199. return 0;
  1200. /* End of SGEGV */
  1201. } /* sgegv_ */