lapack-netlib/SRC/dlag2.c

   1 #include <math.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5 #include <complex.h>
   6 #ifdef complex
   7 #undef complex
   8 #endif
   9 #ifdef I
  10 #undef I
  11 #endif
  12
  13 #if defined(_WIN64)
  14 typedef long long BLASLONG;
  15 typedef unsigned long long BLASULONG;
  16 #else
  17 typedef long BLASLONG;
  18 typedef unsigned long BLASULONG;
  19 #endif
  20
  21 #ifdef LAPACK_ILP64
  22 typedef BLASLONG blasint;
  23 #if defined(_WIN64)
  24 #define blasabs(x) llabs(x)
  25 #else
  26 #define blasabs(x) labs(x)
  27 #endif
  28 #else
  29 typedef int blasint;
  30 #define blasabs(x) abs(x)
  31 #endif
  32
  33 typedef blasint integer;
  34
  35 typedef unsigned int uinteger;
  36 typedef char *address;
  37 typedef short int shortint;
  38 typedef float real;
  39 typedef double doublereal;
  40 typedef struct { real r, i; } complex;
  41 typedef struct { doublereal r, i; } doublecomplex;
  42 #ifdef _MSC_VER
  43 static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
  44 static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
  45 static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
  46 static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
  47 #else
  48 static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
  49 static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
  50 static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
  51 static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
  52 #endif
  53 #define pCf(z) (*_pCf(z))
  54 #define pCd(z) (*_pCd(z))
  55 typedef int logical;
  56 typedef short int shortlogical;
  57 typedef char logical1;
  58 typedef char integer1;
  59
  60 #define TRUE_ (1)
  61 #define FALSE_ (0)
  62
  63 /* Extern is for use with -E */
  64 #ifndef Extern
  65 #define Extern extern
  66 #endif
  67
  68 /* I/O stuff */
  69
  70 typedef int flag;
  71 typedef int ftnlen;
  72 typedef int ftnint;
  73
  74 /*external read, write*/
  75 typedef struct
  76 {       flag cierr;
  77         ftnint ciunit;
  78         flag ciend;
  79         char *cifmt;
  80         ftnint cirec;
  81 } cilist;
  82
  83 /*internal read, write*/
  84 typedef struct
  85 {       flag icierr;
  86         char *iciunit;
  87         flag iciend;
  88         char *icifmt;
  89         ftnint icirlen;
  90         ftnint icirnum;
  91 } icilist;
  92
  93 /*open*/
  94 typedef struct
  95 {       flag oerr;
  96         ftnint ounit;
  97         char *ofnm;
  98         ftnlen ofnmlen;
  99         char *osta;
 100         char *oacc;
 101         char *ofm;
 102         ftnint orl;
 103         char *oblnk;
 104 } olist;
 105
 106 /*close*/
 107 typedef struct
 108 {       flag cerr;
 109         ftnint cunit;
 110         char *csta;
 111 } cllist;
 112
 113 /*rewind, backspace, endfile*/
 114 typedef struct
 115 {       flag aerr;
 116         ftnint aunit;
 117 } alist;
 118
 119 /* inquire */
 120 typedef struct
 121 {       flag inerr;
 122         ftnint inunit;
 123         char *infile;
 124         ftnlen infilen;
 125         ftnint  *inex;  /*parameters in standard's order*/
 126         ftnint  *inopen;
 127         ftnint  *innum;
 128         ftnint  *innamed;
 129         char    *inname;
 130         ftnlen  innamlen;
 131         char    *inacc;
 132         ftnlen  inacclen;
 133         char    *inseq;
 134         ftnlen  inseqlen;
 135         char    *indir;
 136         ftnlen  indirlen;
 137         char    *infmt;
 138         ftnlen  infmtlen;
 139         char    *inform;
 140         ftnint  informlen;
 141         char    *inunf;
 142         ftnlen  inunflen;
 143         ftnint  *inrecl;
 144         ftnint  *innrec;
 145         char    *inblank;
 146         ftnlen  inblanklen;
 147 } inlist;
 148
 149 #define VOID void
 150
 151 union Multitype {       /* for multiple entry points */
 152         integer1 g;
 153         shortint h;
 154         integer i;
 155         /* longint j; */
 156         real r;
 157         doublereal d;
 158         complex c;
 159         doublecomplex z;
 160         };
 161
 162 typedef union Multitype Multitype;
 163
 164 struct Vardesc {        /* for Namelist */
 165         char *name;
 166         char *addr;
 167         ftnlen *dims;
 168         int  type;
 169         };
 170 typedef struct Vardesc Vardesc;
 171
 172 struct Namelist {
 173         char *name;
 174         Vardesc **vars;
 175         int nvars;
 176         };
 177 typedef struct Namelist Namelist;
 178
 179 #define abs(x) ((x) >= 0 ? (x) : -(x))
 180 #define dabs(x) (fabs(x))
 181 #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
 182 #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
 183 #define dmin(a,b) (f2cmin(a,b))
 184 #define dmax(a,b) (f2cmax(a,b))
 185 #define bit_test(a,b)   ((a) >> (b) & 1)
 186 #define bit_clear(a,b)  ((a) & ~((uinteger)1 << (b)))
 187 #define bit_set(a,b)    ((a) |  ((uinteger)1 << (b)))
 188
 189 #define abort_() { sig_die("Fortran abort routine called", 1); }
 190 #define c_abs(z) (cabsf(Cf(z)))
 191 #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
 192 #ifdef _MSC_VER
 193 #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
 194 #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
 195 #else
 196 #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
 197 #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
 198 #endif
 199 #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
 200 #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
 201 #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
 202 //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
 203 #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
 204 #define d_abs(x) (fabs(*(x)))
 205 #define d_acos(x) (acos(*(x)))
 206 #define d_asin(x) (asin(*(x)))
 207 #define d_atan(x) (atan(*(x)))
 208 #define d_atn2(x, y) (atan2(*(x),*(y)))
 209 #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
 210 #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
 211 #define d_cos(x) (cos(*(x)))
 212 #define d_cosh(x) (cosh(*(x)))
 213 #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
 214 #define d_exp(x) (exp(*(x)))
 215 #define d_imag(z) (cimag(Cd(z)))
 216 #define r_imag(z) (cimagf(Cf(z)))
 217 #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
 218 #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
 219 #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
 220 #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
 221 #define d_log(x) (log(*(x)))
 222 #define d_mod(x, y) (fmod(*(x), *(y)))
 223 #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
 224 #define d_nint(x) u_nint(*(x))
 225 #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
 226 #define d_sign(a,b) u_sign(*(a),*(b))
 227 #define r_sign(a,b) u_sign(*(a),*(b))
 228 #define d_sin(x) (sin(*(x)))
 229 #define d_sinh(x) (sinh(*(x)))
 230 #define d_sqrt(x) (sqrt(*(x)))
 231 #define d_tan(x) (tan(*(x)))
 232 #define d_tanh(x) (tanh(*(x)))
 233 #define i_abs(x) abs(*(x))
 234 #define i_dnnt(x) ((integer)u_nint(*(x)))
 235 #define i_len(s, n) (n)
 236 #define i_nint(x) ((integer)u_nint(*(x)))
 237 #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
 238 #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
 239 #define pow_si(B,E) spow_ui(*(B),*(E))
 240 #define pow_ri(B,E) spow_ui(*(B),*(E))
 241 #define pow_di(B,E) dpow_ui(*(B),*(E))
 242 #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
 243 #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
 244 #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
 245 #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++ = ' '; }
 246 #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
 247 #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]; }
 248 #define sig_die(s, kill) { exit(1); }
 249 #define s_stop(s, n) {exit(0);}
 250 static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
 251 #define z_abs(z) (cabs(Cd(z)))
 252 #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
 253 #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
 254 #define myexit_() break;
 255 #define mycycle() continue;
 256 #define myceiling(w) {ceil(w)}
 257 #define myhuge(w) {HUGE_VAL}
 258 //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
 259 #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
 260
 261 /* procedure parameter types for -A and -C++ */
 262
 263 #define F2C_proc_par_types 1
 264 #ifdef __cplusplus
 265 typedef logical (*L_fp)(...);
 266 #else
 267 typedef logical (*L_fp)();
 268 #endif
 269
 270 static float spow_ui(float x, integer n) {
 271         float pow=1.0; unsigned long int u;
 272         if(n != 0) {
 273                 if(n < 0) n = -n, x = 1/x;
 274                 for(u = n; ; ) {
 275                         if(u & 01) pow *= x;
 276                         if(u >>= 1) x *= x;
 277                         else break;
 278                 }
 279         }
 280         return pow;
 281 }
 282 static double dpow_ui(double x, integer n) {
 283         double pow=1.0; unsigned long int u;
 284         if(n != 0) {
 285                 if(n < 0) n = -n, x = 1/x;
 286                 for(u = n; ; ) {
 287                         if(u & 01) pow *= x;
 288                         if(u >>= 1) x *= x;
 289                         else break;
 290                 }
 291         }
 292         return pow;
 293 }
 294 #ifdef _MSC_VER
 295 static _Fcomplex cpow_ui(complex x, integer n) {
 296         complex pow={1.0,0.0}; unsigned long int u;
 297                 if(n != 0) {
 298                 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
 299                 for(u = n; ; ) {
 300                         if(u & 01) pow.r *= x.r, pow.i *= x.i;
 301                         if(u >>= 1) x.r *= x.r, x.i *= x.i;
 302                         else break;
 303                 }
 304         }
 305         _Fcomplex p={pow.r, pow.i};
 306         return p;
 307 }
 308 #else
 309 static _Complex float cpow_ui(_Complex float x, integer n) {
 310         _Complex float pow=1.0; unsigned long int u;
 311         if(n != 0) {
 312                 if(n < 0) n = -n, x = 1/x;
 313                 for(u = n; ; ) {
 314                         if(u & 01) pow *= x;
 315                         if(u >>= 1) x *= x;
 316                         else break;
 317                 }
 318         }
 319         return pow;
 320 }
 321 #endif
 322 #ifdef _MSC_VER
 323 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
 324         _Dcomplex pow={1.0,0.0}; unsigned long int u;
 325         if(n != 0) {
 326                 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
 327                 for(u = n; ; ) {
 328                         if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
 329                         if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
 330                         else break;
 331                 }
 332         }
 333         _Dcomplex p = {pow._Val[0], pow._Val[1]};
 334         return p;
 335 }
 336 #else
 337 static _Complex double zpow_ui(_Complex double x, integer n) {
 338         _Complex double pow=1.0; unsigned long int u;
 339         if(n != 0) {
 340                 if(n < 0) n = -n, x = 1/x;
 341                 for(u = n; ; ) {
 342                         if(u & 01) pow *= x;
 343                         if(u >>= 1) x *= x;
 344                         else break;
 345                 }
 346         }
 347         return pow;
 348 }
 349 #endif
 350 static integer pow_ii(integer x, integer n) {
 351         integer pow; unsigned long int u;
 352         if (n <= 0) {
 353                 if (n == 0 || x == 1) pow = 1;
 354                 else if (x != -1) pow = x == 0 ? 1/x : 0;
 355                 else n = -n;
 356         }
 357         if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
 358                 u = n;
 359                 for(pow = 1; ; ) {
 360                         if(u & 01) pow *= x;
 361                         if(u >>= 1) x *= x;
 362                         else break;
 363                 }
 364         }
 365         return pow;
 366 }
 367 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
 368 {
 369         double m; integer i, mi;
 370         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
 371                 if (w[i-1]>m) mi=i ,m=w[i-1];
 372         return mi-s+1;
 373 }
 374 static integer smaxloc_(float *w, integer s, integer e, integer *n)
 375 {
 376         float m; integer i, mi;
 377         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
 378                 if (w[i-1]>m) mi=i ,m=w[i-1];
 379         return mi-s+1;
 380 }
 381 static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
 382         integer n = *n_, incx = *incx_, incy = *incy_, i;
 383 #ifdef _MSC_VER
 384         _Fcomplex zdotc = {0.0, 0.0};
 385         if (incx == 1 && incy == 1) {
 386                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 387                         zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
 388                         zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
 389                 }
 390         } else {
 391                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 392                         zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
 393                         zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
 394                 }
 395         }
 396         pCf(z) = zdotc;
 397 }
 398 #else
 399         _Complex float zdotc = 0.0;
 400         if (incx == 1 && incy == 1) {
 401                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 402                         zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
 403                 }
 404         } else {
 405                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 406                         zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
 407                 }
 408         }
 409         pCf(z) = zdotc;
 410 }
 411 #endif
 412 static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
 413         integer n = *n_, incx = *incx_, incy = *incy_, i;
 414 #ifdef _MSC_VER
 415         _Dcomplex zdotc = {0.0, 0.0};
 416         if (incx == 1 && incy == 1) {
 417                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 418                         zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
 419                         zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
 420                 }
 421         } else {
 422                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 423                         zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
 424                         zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
 425                 }
 426         }
 427         pCd(z) = zdotc;
 428 }
 429 #else
 430         _Complex double zdotc = 0.0;
 431         if (incx == 1 && incy == 1) {
 432                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 433                         zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
 434                 }
 435         } else {
 436                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 437                         zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
 438                 }
 439         }
 440         pCd(z) = zdotc;
 441 }
 442 #endif
 443 static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
 444         integer n = *n_, incx = *incx_, incy = *incy_, i;
 445 #ifdef _MSC_VER
 446         _Fcomplex zdotc = {0.0, 0.0};
 447         if (incx == 1 && incy == 1) {
 448                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 449                         zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
 450                         zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
 451                 }
 452         } else {
 453                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 454                         zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
 455                         zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
 456                 }
 457         }
 458         pCf(z) = zdotc;
 459 }
 460 #else
 461         _Complex float zdotc = 0.0;
 462         if (incx == 1 && incy == 1) {
 463                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 464                         zdotc += Cf(&x[i]) * Cf(&y[i]);
 465                 }
 466         } else {
 467                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 468                         zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
 469                 }
 470         }
 471         pCf(z) = zdotc;
 472 }
 473 #endif
 474 static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
 475         integer n = *n_, incx = *incx_, incy = *incy_, i;
 476 #ifdef _MSC_VER
 477         _Dcomplex zdotc = {0.0, 0.0};
 478         if (incx == 1 && incy == 1) {
 479                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 480                         zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
 481                         zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
 482                 }
 483         } else {
 484                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 485                         zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
 486                         zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
 487                 }
 488         }
 489         pCd(z) = zdotc;
 490 }
 491 #else
 492         _Complex double zdotc = 0.0;
 493         if (incx == 1 && incy == 1) {
 494                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 495                         zdotc += Cd(&x[i]) * Cd(&y[i]);
 496                 }
 497         } else {
 498                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 499                         zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
 500                 }
 501         }
 502         pCd(z) = zdotc;
 503 }
 504 #endif
 505 /*  -- translated by f2c (version 20000121).
 506    You must link the resulting object file with the libraries:
 507         -lf2c -lm   (in that order)
 508 */
 509
 510
 511
 512
 513 /* > \brief \b DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as nece
 514 ssary to avoid over-/underflow. */
 515
 516 /*  =========== DOCUMENTATION =========== */
 517
 518 /* Online html documentation available at */
 519 /*            http://www.netlib.org/lapack/explore-html/ */
 520
 521 /* > \htmlonly */
 522 /* > Download DLAG2 + dependencies */
 523 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlag2.f
 524 "> */
 525 /* > [TGZ]</a> */
 526 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlag2.f
 527 "> */
 528 /* > [ZIP]</a> */
 529 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlag2.f
 530 "> */
 531 /* > [TXT]</a> */
 532 /* > \endhtmlonly */
 533
 534 /*  Definition: */
 535 /*  =========== */
 536
 537 /*       SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, */
 538 /*                         WR2, WI ) */
 539
 540 /*       INTEGER            LDA, LDB */
 541 /*       DOUBLE PRECISION   SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 */
 542 /*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ) */
 543
 544
 545 /* > \par Purpose: */
 546 /*  ============= */
 547 /* > */
 548 /* > \verbatim */
 549 /* > */
 550 /* > DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue */
 551 /* > problem  A - w B, with scaling as necessary to avoid over-/underflow. */
 552 /* > */
 553 /* > The scaling factor "s" results in a modified eigenvalue equation */
 554 /* > */
 555 /* >     s A - w B */
 556 /* > */
 557 /* > where  s  is a non-negative scaling factor chosen so that  w,  w B, */
 558 /* > and  s A  do not overflow and, if possible, do not underflow, either. */
 559 /* > \endverbatim */
 560
 561 /*  Arguments: */
 562 /*  ========== */
 563
 564 /* > \param[in] A */
 565 /* > \verbatim */
 566 /* >          A is DOUBLE PRECISION array, dimension (LDA, 2) */
 567 /* >          On entry, the 2 x 2 matrix A.  It is assumed that its 1-norm */
 568 /* >          is less than 1/SAFMIN.  Entries less than */
 569 /* >          sqrt(SAFMIN)*norm(A) are subject to being treated as zero. */
 570 /* > \endverbatim */
 571 /* > */
 572 /* > \param[in] LDA */
 573 /* > \verbatim */
 574 /* >          LDA is INTEGER */
 575 /* >          The leading dimension of the array A.  LDA >= 2. */
 576 /* > \endverbatim */
 577 /* > */
 578 /* > \param[in] B */
 579 /* > \verbatim */
 580 /* >          B is DOUBLE PRECISION array, dimension (LDB, 2) */
 581 /* >          On entry, the 2 x 2 upper triangular matrix B.  It is */
 582 /* >          assumed that the one-norm of B is less than 1/SAFMIN.  The */
 583 /* >          diagonals should be at least sqrt(SAFMIN) times the largest */
 584 /* >          element of B (in absolute value); if a diagonal is smaller */
 585 /* >          than that, then  +/- sqrt(SAFMIN) will be used instead of */
 586 /* >          that diagonal. */
 587 /* > \endverbatim */
 588 /* > */
 589 /* > \param[in] LDB */
 590 /* > \verbatim */
 591 /* >          LDB is INTEGER */
 592 /* >          The leading dimension of the array B.  LDB >= 2. */
 593 /* > \endverbatim */
 594 /* > */
 595 /* > \param[in] SAFMIN */
 596 /* > \verbatim */
 597 /* >          SAFMIN is DOUBLE PRECISION */
 598 /* >          The smallest positive number s.t. 1/SAFMIN does not */
 599 /* >          overflow.  (This should always be DLAMCH('S') -- it is an */
 600 /* >          argument in order to avoid having to call DLAMCH frequently.) */
 601 /* > \endverbatim */
 602 /* > */
 603 /* > \param[out] SCALE1 */
 604 /* > \verbatim */
 605 /* >          SCALE1 is DOUBLE PRECISION */
 606 /* >          A scaling factor used to avoid over-/underflow in the */
 607 /* >          eigenvalue equation which defines the first eigenvalue.  If */
 608 /* >          the eigenvalues are complex, then the eigenvalues are */
 609 /* >          ( WR1  +/-  WI i ) / SCALE1  (which may lie outside the */
 610 /* >          exponent range of the machine), SCALE1=SCALE2, and SCALE1 */
 611 /* >          will always be positive.  If the eigenvalues are real, then */
 612 /* >          the first (real) eigenvalue is  WR1 / SCALE1 , but this may */
 613 /* >          overflow or underflow, and in fact, SCALE1 may be zero or */
 614 /* >          less than the underflow threshold if the exact eigenvalue */
 615 /* >          is sufficiently large. */
 616 /* > \endverbatim */
 617 /* > */
 618 /* > \param[out] SCALE2 */
 619 /* > \verbatim */
 620 /* >          SCALE2 is DOUBLE PRECISION */
 621 /* >          A scaling factor used to avoid over-/underflow in the */
 622 /* >          eigenvalue equation which defines the second eigenvalue.  If */
 623 /* >          the eigenvalues are complex, then SCALE2=SCALE1.  If the */
 624 /* >          eigenvalues are real, then the second (real) eigenvalue is */
 625 /* >          WR2 / SCALE2 , but this may overflow or underflow, and in */
 626 /* >          fact, SCALE2 may be zero or less than the underflow */
 627 /* >          threshold if the exact eigenvalue is sufficiently large. */
 628 /* > \endverbatim */
 629 /* > */
 630 /* > \param[out] WR1 */
 631 /* > \verbatim */
 632 /* >          WR1 is DOUBLE PRECISION */
 633 /* >          If the eigenvalue is real, then WR1 is SCALE1 times the */
 634 /* >          eigenvalue closest to the (2,2) element of A B**(-1).  If the */
 635 /* >          eigenvalue is complex, then WR1=WR2 is SCALE1 times the real */
 636 /* >          part of the eigenvalues. */
 637 /* > \endverbatim */
 638 /* > */
 639 /* > \param[out] WR2 */
 640 /* > \verbatim */
 641 /* >          WR2 is DOUBLE PRECISION */
 642 /* >          If the eigenvalue is real, then WR2 is SCALE2 times the */
 643 /* >          other eigenvalue.  If the eigenvalue is complex, then */
 644 /* >          WR1=WR2 is SCALE1 times the real part of the eigenvalues. */
 645 /* > \endverbatim */
 646 /* > */
 647 /* > \param[out] WI */
 648 /* > \verbatim */
 649 /* >          WI is DOUBLE PRECISION */
 650 /* >          If the eigenvalue is real, then WI is zero.  If the */
 651 /* >          eigenvalue is complex, then WI is SCALE1 times the imaginary */
 652 /* >          part of the eigenvalues.  WI will always be non-negative. */
 653 /* > \endverbatim */
 654
 655 /*  Authors: */
 656 /*  ======== */
 657
 658 /* > \author Univ. of Tennessee */
 659 /* > \author Univ. of California Berkeley */
 660 /* > \author Univ. of Colorado Denver */
 661 /* > \author NAG Ltd. */
 662
 663 /* > \date June 2016 */
 664
 665 /* > \ingroup doubleOTHERauxiliary */
 666
 667 /*  ===================================================================== */
 668 /* Subroutine */ int dlag2_(doublereal *a, integer *lda, doublereal *b,
 669         integer *ldb, doublereal *safmin, doublereal *scale1, doublereal *
 670         scale2, doublereal *wr1, doublereal *wr2, doublereal *wi)
 671 {
 672     /* System generated locals */
 673     integer a_dim1, a_offset, b_dim1, b_offset;
 674     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
 675
 676     /* Local variables */
 677     doublereal diff, bmin, wbig, wabs, wdet, r__, binv11, binv22, discr,
 678             anorm, bnorm, bsize, shift, c1, c2, c3, c4, c5, rtmin, rtmax,
 679             wsize, s1, s2, a11, a12, a21, a22, b11, b12, b22, ascale, bscale,
 680             pp, qq, ss, wscale, safmax, wsmall, as11, as12, as22, sum, abi22;
 681
 682
 683 /*  -- LAPACK auxiliary routine (version 3.7.0) -- */
 684 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
 685 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
 686 /*     June 2016 */
 687
 688
 689 /*  ===================================================================== */
 690
 691
 692     /* Parameter adjustments */
 693     a_dim1 = *lda;
 694     a_offset = 1 + a_dim1 * 1;
 695     a -= a_offset;
 696     b_dim1 = *ldb;
 697     b_offset = 1 + b_dim1 * 1;
 698     b -= b_offset;
 699
 700     /* Function Body */
 701     rtmin = sqrt(*safmin);
 702     rtmax = 1. / rtmin;
 703     safmax = 1. / *safmin;
 704
 705 /*     Scale A */
 706
 707 /* Computing MAX */
 708     d__5 = (d__1 = a[a_dim1 + 1], abs(d__1)) + (d__2 = a[a_dim1 + 2], abs(
 709             d__2)), d__6 = (d__3 = a[(a_dim1 << 1) + 1], abs(d__3)) + (d__4 =
 710             a[(a_dim1 << 1) + 2], abs(d__4)), d__5 = f2cmax(d__5,d__6);
 711     anorm = f2cmax(d__5,*safmin);
 712     ascale = 1. / anorm;
 713     a11 = ascale * a[a_dim1 + 1];
 714     a21 = ascale * a[a_dim1 + 2];
 715     a12 = ascale * a[(a_dim1 << 1) + 1];
 716     a22 = ascale * a[(a_dim1 << 1) + 2];
 717
 718 /*     Perturb B if necessary to insure non-singularity */
 719
 720     b11 = b[b_dim1 + 1];
 721     b12 = b[(b_dim1 << 1) + 1];
 722     b22 = b[(b_dim1 << 1) + 2];
 723 /* Computing MAX */
 724     d__1 = abs(b11), d__2 = abs(b12), d__1 = f2cmax(d__1,d__2), d__2 = abs(b22),
 725             d__1 = f2cmax(d__1,d__2);
 726     bmin = rtmin * f2cmax(d__1,rtmin);
 727     if (abs(b11) < bmin) {
 728         b11 = d_sign(&bmin, &b11);
 729     }
 730     if (abs(b22) < bmin) {
 731         b22 = d_sign(&bmin, &b22);
 732     }
 733
 734 /*     Scale B */
 735
 736 /* Computing MAX */
 737     d__1 = abs(b11), d__2 = abs(b12) + abs(b22), d__1 = f2cmax(d__1,d__2);
 738     bnorm = f2cmax(d__1,*safmin);
 739 /* Computing MAX */
 740     d__1 = abs(b11), d__2 = abs(b22);
 741     bsize = f2cmax(d__1,d__2);
 742     bscale = 1. / bsize;
 743     b11 *= bscale;
 744     b12 *= bscale;
 745     b22 *= bscale;
 746
 747 /*     Compute larger eigenvalue by method described by C. van Loan */
 748
 749 /*     ( AS is A shifted by -SHIFT*B ) */
 750
 751     binv11 = 1. / b11;
 752     binv22 = 1. / b22;
 753     s1 = a11 * binv11;
 754     s2 = a22 * binv22;
 755     if (abs(s1) <= abs(s2)) {
 756         as12 = a12 - s1 * b12;
 757         as22 = a22 - s1 * b22;
 758         ss = a21 * (binv11 * binv22);
 759         abi22 = as22 * binv22 - ss * b12;
 760         pp = abi22 * .5;
 761         shift = s1;
 762     } else {
 763         as12 = a12 - s2 * b12;
 764         as11 = a11 - s2 * b11;
 765         ss = a21 * (binv11 * binv22);
 766         abi22 = -ss * b12;
 767         pp = (as11 * binv11 + abi22) * .5;
 768         shift = s2;
 769     }
 770     qq = ss * as12;
 771     if ((d__1 = pp * rtmin, abs(d__1)) >= 1.) {
 772 /* Computing 2nd power */
 773         d__1 = rtmin * pp;
 774         discr = d__1 * d__1 + qq * *safmin;
 775         r__ = sqrt((abs(discr))) * rtmax;
 776     } else {
 777 /* Computing 2nd power */
 778         d__1 = pp;
 779         if (d__1 * d__1 + abs(qq) <= *safmin) {
 780 /* Computing 2nd power */
 781             d__1 = rtmax * pp;
 782             discr = d__1 * d__1 + qq * safmax;
 783             r__ = sqrt((abs(discr))) * rtmin;
 784         } else {
 785 /* Computing 2nd power */
 786             d__1 = pp;
 787             discr = d__1 * d__1 + qq;
 788             r__ = sqrt((abs(discr)));
 789         }
 790     }
 791
 792 /*     Note: the test of R in the following IF is to cover the case when */
 793 /*           DISCR is small and negative and is flushed to zero during */
 794 /*           the calculation of R.  On machines which have a consistent */
 795 /*           flush-to-zero threshold and handle numbers above that */
 796 /*           threshold correctly, it would not be necessary. */
 797
 798     if (discr >= 0. || r__ == 0.) {
 799         sum = pp + d_sign(&r__, &pp);
 800         diff = pp - d_sign(&r__, &pp);
 801         wbig = shift + sum;
 802
 803 /*        Compute smaller eigenvalue */
 804
 805         wsmall = shift + diff;
 806 /* Computing MAX */
 807         d__1 = abs(wsmall);
 808         if (abs(wbig) * .5 > f2cmax(d__1,*safmin)) {
 809             wdet = (a11 * a22 - a12 * a21) * (binv11 * binv22);
 810             wsmall = wdet / wbig;
 811         }
 812
 813 /*        Choose (real) eigenvalue closest to 2,2 element of A*B**(-1) */
 814 /*        for WR1. */
 815
 816         if (pp > abi22) {
 817             *wr1 = f2cmin(wbig,wsmall);
 818             *wr2 = f2cmax(wbig,wsmall);
 819         } else {
 820             *wr1 = f2cmax(wbig,wsmall);
 821             *wr2 = f2cmin(wbig,wsmall);
 822         }
 823         *wi = 0.;
 824     } else {
 825
 826 /*        Complex eigenvalues */
 827
 828         *wr1 = shift + pp;
 829         *wr2 = *wr1;
 830         *wi = r__;
 831     }
 832
 833 /*     Further scaling to avoid underflow and overflow in computing */
 834 /*     SCALE1 and overflow in computing w*B. */
 835
 836 /*     This scale factor (WSCALE) is bounded from above using C1 and C2, */
 837 /*     and from below using C3 and C4. */
 838 /*        C1 implements the condition  s A  must never overflow. */
 839 /*        C2 implements the condition  w B  must never overflow. */
 840 /*        C3, with C2, */
 841 /*           implement the condition that s A - w B must never overflow. */
 842 /*        C4 implements the condition  s    should not underflow. */
 843 /*        C5 implements the condition  f2cmax(s,|w|) should be at least 2. */
 844
 845     c1 = bsize * (*safmin * f2cmax(1.,ascale));
 846     c2 = *safmin * f2cmax(1.,bnorm);
 847     c3 = bsize * *safmin;
 848     if (ascale <= 1. && bsize <= 1.) {
 849 /* Computing MIN */
 850         d__1 = 1., d__2 = ascale / *safmin * bsize;
 851         c4 = f2cmin(d__1,d__2);
 852     } else {
 853         c4 = 1.;
 854     }
 855     if (ascale <= 1. || bsize <= 1.) {
 856 /* Computing MIN */
 857         d__1 = 1., d__2 = ascale * bsize;
 858         c5 = f2cmin(d__1,d__2);
 859     } else {
 860         c5 = 1.;
 861     }
 862
 863 /*     Scale first eigenvalue */
 864
 865     wabs = abs(*wr1) + abs(*wi);
 866 /* Computing MAX */
 867 /* Computing MIN */
 868     d__3 = c4, d__4 = f2cmax(wabs,c5) * .5;
 869     d__1 = f2cmax(*safmin,c1), d__2 = (wabs * c2 + c3) * 1.0000100000000001,
 870             d__1 = f2cmax(d__1,d__2), d__2 = f2cmin(d__3,d__4);
 871     wsize = f2cmax(d__1,d__2);
 872     if (wsize != 1.) {
 873         wscale = 1. / wsize;
 874         if (wsize > 1.) {
 875             *scale1 = f2cmax(ascale,bsize) * wscale * f2cmin(ascale,bsize);
 876         } else {
 877             *scale1 = f2cmin(ascale,bsize) * wscale * f2cmax(ascale,bsize);
 878         }
 879         *wr1 *= wscale;
 880         if (*wi != 0.) {
 881             *wi *= wscale;
 882             *wr2 = *wr1;
 883             *scale2 = *scale1;
 884         }
 885     } else {
 886         *scale1 = ascale * bsize;
 887         *scale2 = *scale1;
 888     }
 889
 890 /*     Scale second eigenvalue (if real) */
 891
 892     if (*wi == 0.) {
 893 /* Computing MAX */
 894 /* Computing MIN */
 895 /* Computing MAX */
 896         d__5 = abs(*wr2);
 897         d__3 = c4, d__4 = f2cmax(d__5,c5) * .5;
 898         d__1 = f2cmax(*safmin,c1), d__2 = (abs(*wr2) * c2 + c3) *
 899                 1.0000100000000001, d__1 = f2cmax(d__1,d__2), d__2 = f2cmin(d__3,
 900                 d__4);
 901         wsize = f2cmax(d__1,d__2);
 902         if (wsize != 1.) {
 903             wscale = 1. / wsize;
 904             if (wsize > 1.) {
 905                 *scale2 = f2cmax(ascale,bsize) * wscale * f2cmin(ascale,bsize);
 906             } else {
 907                 *scale2 = f2cmin(ascale,bsize) * wscale * f2cmax(ascale,bsize);
 908             }
 909             *wr2 *= wscale;
 910         } else {
 911             *scale2 = ascale * bsize;
 912         }
 913     }
 914
 915 /*     End of DLAG2 */
 916
 917     return 0;
 918 } /* dlag2_ */
 919