lapack-netlib/SRC/iparam2stage.c

   1 #include <math.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5 #include <ctype.h>
   6 #include <complex.h>
   7 #ifdef complex
   8 #undef complex
   9 #endif
  10 #ifdef I
  11 #undef I
  12 #endif
  13
  14 #if defined(_WIN64)
  15 typedef long long BLASLONG;
  16 typedef unsigned long long BLASULONG;
  17 #else
  18 typedef long BLASLONG;
  19 typedef unsigned long BLASULONG;
  20 #endif
  21
  22 #ifdef LAPACK_ILP64
  23 typedef BLASLONG blasint;
  24 #if defined(_WIN64)
  25 #define blasabs(x) llabs(x)
  26 #else
  27 #define blasabs(x) labs(x)
  28 #endif
  29 #else
  30 typedef int blasint;
  31 #define blasabs(x) abs(x)
  32 #endif
  33
  34 typedef blasint integer;
  35
  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 #ifdef _MSC_VER
  44 static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
  45 static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
  46 static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
  47 static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
  48 #else
  49 static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
  50 static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
  51 static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
  52 static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
  53 #endif
  54 #define pCf(z) (*_pCf(z))
  55 #define pCd(z) (*_pCd(z))
  56 typedef int logical;
  57 typedef short int shortlogical;
  58 typedef char logical1;
  59 typedef char integer1;
  60
  61 #define TRUE_ (1)
  62 #define FALSE_ (0)
  63
  64 /* Extern is for use with -E */
  65 #ifndef Extern
  66 #define Extern extern
  67 #endif
  68
  69 /* I/O stuff */
  70
  71 typedef int flag;
  72 typedef int ftnlen;
  73 typedef int ftnint;
  74
  75 /*external read, write*/
  76 typedef struct
  77 {       flag cierr;
  78         ftnint ciunit;
  79         flag ciend;
  80         char *cifmt;
  81         ftnint cirec;
  82 } cilist;
  83
  84 /*internal read, write*/
  85 typedef struct
  86 {       flag icierr;
  87         char *iciunit;
  88         flag iciend;
  89         char *icifmt;
  90         ftnint icirlen;
  91         ftnint icirnum;
  92 } icilist;
  93
  94 /*open*/
  95 typedef struct
  96 {       flag oerr;
  97         ftnint ounit;
  98         char *ofnm;
  99         ftnlen ofnmlen;
 100         char *osta;
 101         char *oacc;
 102         char *ofm;
 103         ftnint orl;
 104         char *oblnk;
 105 } olist;
 106
 107 /*close*/
 108 typedef struct
 109 {       flag cerr;
 110         ftnint cunit;
 111         char *csta;
 112 } cllist;
 113
 114 /*rewind, backspace, endfile*/
 115 typedef struct
 116 {       flag aerr;
 117         ftnint aunit;
 118 } alist;
 119
 120 /* inquire */
 121 typedef struct
 122 {       flag inerr;
 123         ftnint inunit;
 124         char *infile;
 125         ftnlen infilen;
 126         ftnint  *inex;  /*parameters in standard's order*/
 127         ftnint  *inopen;
 128         ftnint  *innum;
 129         ftnint  *innamed;
 130         char    *inname;
 131         ftnlen  innamlen;
 132         char    *inacc;
 133         ftnlen  inacclen;
 134         char    *inseq;
 135         ftnlen  inseqlen;
 136         char    *indir;
 137         ftnlen  indirlen;
 138         char    *infmt;
 139         ftnlen  infmtlen;
 140         char    *inform;
 141         ftnint  informlen;
 142         char    *inunf;
 143         ftnlen  inunflen;
 144         ftnint  *inrecl;
 145         ftnint  *innrec;
 146         char    *inblank;
 147         ftnlen  inblanklen;
 148 } inlist;
 149
 150 #define VOID void
 151
 152 union Multitype {       /* for multiple entry points */
 153         integer1 g;
 154         shortint h;
 155         integer i;
 156         /* longint j; */
 157         real r;
 158         doublereal d;
 159         complex c;
 160         doublecomplex z;
 161         };
 162
 163 typedef union Multitype Multitype;
 164
 165 struct Vardesc {        /* for Namelist */
 166         char *name;
 167         char *addr;
 168         ftnlen *dims;
 169         int  type;
 170         };
 171 typedef struct Vardesc Vardesc;
 172
 173 struct Namelist {
 174         char *name;
 175         Vardesc **vars;
 176         int nvars;
 177         };
 178 typedef struct Namelist Namelist;
 179
 180 #define abs(x) ((x) >= 0 ? (x) : -(x))
 181 #define dabs(x) (fabs(x))
 182 #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
 183 #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
 184 #define dmin(a,b) (f2cmin(a,b))
 185 #define dmax(a,b) (f2cmax(a,b))
 186 #define bit_test(a,b)   ((a) >> (b) & 1)
 187 #define bit_clear(a,b)  ((a) & ~((uinteger)1 << (b)))
 188 #define bit_set(a,b)    ((a) |  ((uinteger)1 << (b)))
 189
 190 #define abort_() { sig_die("Fortran abort routine called", 1); }
 191 #define c_abs(z) (cabsf(Cf(z)))
 192 #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
 193 #ifdef _MSC_VER
 194 #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]);}
 195 #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]);}
 196 #else
 197 #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
 198 #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
 199 #endif
 200 #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
 201 #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
 202 #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
 203 //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
 204 #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
 205 #define d_abs(x) (fabs(*(x)))
 206 #define d_acos(x) (acos(*(x)))
 207 #define d_asin(x) (asin(*(x)))
 208 #define d_atan(x) (atan(*(x)))
 209 #define d_atn2(x, y) (atan2(*(x),*(y)))
 210 #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
 211 #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
 212 #define d_cos(x) (cos(*(x)))
 213 #define d_cosh(x) (cosh(*(x)))
 214 #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
 215 #define d_exp(x) (exp(*(x)))
 216 #define d_imag(z) (cimag(Cd(z)))
 217 #define r_imag(z) (cimagf(Cf(z)))
 218 #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
 219 #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
 220 #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
 221 #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
 222 #define d_log(x) (log(*(x)))
 223 #define d_mod(x, y) (fmod(*(x), *(y)))
 224 #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
 225 #define d_nint(x) u_nint(*(x))
 226 #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
 227 #define d_sign(a,b) u_sign(*(a),*(b))
 228 #define r_sign(a,b) u_sign(*(a),*(b))
 229 #define d_sin(x) (sin(*(x)))
 230 #define d_sinh(x) (sinh(*(x)))
 231 #define d_sqrt(x) (sqrt(*(x)))
 232 #define d_tan(x) (tan(*(x)))
 233 #define d_tanh(x) (tanh(*(x)))
 234 #define i_abs(x) abs(*(x))
 235 #define i_dnnt(x) ((integer)u_nint(*(x)))
 236 #define i_len(s, n) (n)
 237 #define i_nint(x) ((integer)u_nint(*(x)))
 238 #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
 239 #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
 240 #define pow_si(B,E) spow_ui(*(B),*(E))
 241 #define pow_ri(B,E) spow_ui(*(B),*(E))
 242 #define pow_di(B,E) dpow_ui(*(B),*(E))
 243 #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
 244 #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
 245 #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
 246 #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++ = ' '; }
 247 #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
 248 #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]; }
 249 #define sig_die(s, kill) { exit(1); }
 250 #define s_stop(s, n) {exit(0);}
 251 static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
 252 #define z_abs(z) (cabs(Cd(z)))
 253 #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
 254 #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
 255 #define myexit_() break;
 256 #define mycycle() continue;
 257 #define myceiling(w) {ceil(w)}
 258 #define myhuge(w) {HUGE_VAL}
 259 //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
 260 #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
 261
 262 /* procedure parameter types for -A and -C++ */
 263
 264 #define F2C_proc_par_types 1
 265 #ifdef __cplusplus
 266 typedef logical (*L_fp)(...);
 267 #else
 268 typedef logical (*L_fp)();
 269 #endif
 270
 271 static float spow_ui(float x, integer n) {
 272         float pow=1.0; unsigned long int u;
 273         if(n != 0) {
 274                 if(n < 0) n = -n, x = 1/x;
 275                 for(u = n; ; ) {
 276                         if(u & 01) pow *= x;
 277                         if(u >>= 1) x *= x;
 278                         else break;
 279                 }
 280         }
 281         return pow;
 282 }
 283 static double dpow_ui(double x, integer n) {
 284         double pow=1.0; unsigned long int u;
 285         if(n != 0) {
 286                 if(n < 0) n = -n, x = 1/x;
 287                 for(u = n; ; ) {
 288                         if(u & 01) pow *= x;
 289                         if(u >>= 1) x *= x;
 290                         else break;
 291                 }
 292         }
 293         return pow;
 294 }
 295 #ifdef _MSC_VER
 296 static _Fcomplex cpow_ui(complex x, integer n) {
 297         complex pow={1.0,0.0}; unsigned long int u;
 298                 if(n != 0) {
 299                 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
 300                 for(u = n; ; ) {
 301                         if(u & 01) pow.r *= x.r, pow.i *= x.i;
 302                         if(u >>= 1) x.r *= x.r, x.i *= x.i;
 303                         else break;
 304                 }
 305         }
 306         _Fcomplex p={pow.r, pow.i};
 307         return p;
 308 }
 309 #else
 310 static _Complex float cpow_ui(_Complex float x, integer n) {
 311         _Complex float pow=1.0; unsigned long int u;
 312         if(n != 0) {
 313                 if(n < 0) n = -n, x = 1/x;
 314                 for(u = n; ; ) {
 315                         if(u & 01) pow *= x;
 316                         if(u >>= 1) x *= x;
 317                         else break;
 318                 }
 319         }
 320         return pow;
 321 }
 322 #endif
 323 #ifdef _MSC_VER
 324 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
 325         _Dcomplex pow={1.0,0.0}; unsigned long int u;
 326         if(n != 0) {
 327                 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
 328                 for(u = n; ; ) {
 329                         if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
 330                         if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
 331                         else break;
 332                 }
 333         }
 334         _Dcomplex p = {pow._Val[0], pow._Val[1]};
 335         return p;
 336 }
 337 #else
 338 static _Complex double zpow_ui(_Complex double x, integer n) {
 339         _Complex double pow=1.0; unsigned long int u;
 340         if(n != 0) {
 341                 if(n < 0) n = -n, x = 1/x;
 342                 for(u = n; ; ) {
 343                         if(u & 01) pow *= x;
 344                         if(u >>= 1) x *= x;
 345                         else break;
 346                 }
 347         }
 348         return pow;
 349 }
 350 #endif
 351 static integer pow_ii(integer x, integer n) {
 352         integer pow; unsigned long int u;
 353         if (n <= 0) {
 354                 if (n == 0 || x == 1) pow = 1;
 355                 else if (x != -1) pow = x == 0 ? 1/x : 0;
 356                 else n = -n;
 357         }
 358         if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
 359                 u = n;
 360                 for(pow = 1; ; ) {
 361                         if(u & 01) pow *= x;
 362                         if(u >>= 1) x *= x;
 363                         else break;
 364                 }
 365         }
 366         return pow;
 367 }
 368 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
 369 {
 370         double m; integer i, mi;
 371         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
 372                 if (w[i-1]>m) mi=i ,m=w[i-1];
 373         return mi-s+1;
 374 }
 375 static integer smaxloc_(float *w, integer s, integer e, integer *n)
 376 {
 377         float m; integer i, mi;
 378         for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
 379                 if (w[i-1]>m) mi=i ,m=w[i-1];
 380         return mi-s+1;
 381 }
 382 static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
 383         integer n = *n_, incx = *incx_, incy = *incy_, i;
 384 #ifdef _MSC_VER
 385         _Fcomplex zdotc = {0.0, 0.0};
 386         if (incx == 1 && incy == 1) {
 387                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 388                         zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
 389                         zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
 390                 }
 391         } else {
 392                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 393                         zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
 394                         zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
 395                 }
 396         }
 397         pCf(z) = zdotc;
 398 }
 399 #else
 400         _Complex float zdotc = 0.0;
 401         if (incx == 1 && incy == 1) {
 402                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 403                         zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
 404                 }
 405         } else {
 406                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 407                         zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
 408                 }
 409         }
 410         pCf(z) = zdotc;
 411 }
 412 #endif
 413 static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
 414         integer n = *n_, incx = *incx_, incy = *incy_, i;
 415 #ifdef _MSC_VER
 416         _Dcomplex zdotc = {0.0, 0.0};
 417         if (incx == 1 && incy == 1) {
 418                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 419                         zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
 420                         zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
 421                 }
 422         } else {
 423                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 424                         zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
 425                         zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
 426                 }
 427         }
 428         pCd(z) = zdotc;
 429 }
 430 #else
 431         _Complex double zdotc = 0.0;
 432         if (incx == 1 && incy == 1) {
 433                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 434                         zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
 435                 }
 436         } else {
 437                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 438                         zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
 439                 }
 440         }
 441         pCd(z) = zdotc;
 442 }
 443 #endif
 444 static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
 445         integer n = *n_, incx = *incx_, incy = *incy_, i;
 446 #ifdef _MSC_VER
 447         _Fcomplex zdotc = {0.0, 0.0};
 448         if (incx == 1 && incy == 1) {
 449                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 450                         zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
 451                         zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
 452                 }
 453         } else {
 454                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 455                         zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
 456                         zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
 457                 }
 458         }
 459         pCf(z) = zdotc;
 460 }
 461 #else
 462         _Complex float zdotc = 0.0;
 463         if (incx == 1 && incy == 1) {
 464                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 465                         zdotc += Cf(&x[i]) * Cf(&y[i]);
 466                 }
 467         } else {
 468                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 469                         zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
 470                 }
 471         }
 472         pCf(z) = zdotc;
 473 }
 474 #endif
 475 static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
 476         integer n = *n_, incx = *incx_, incy = *incy_, i;
 477 #ifdef _MSC_VER
 478         _Dcomplex zdotc = {0.0, 0.0};
 479         if (incx == 1 && incy == 1) {
 480                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 481                         zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
 482                         zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
 483                 }
 484         } else {
 485                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 486                         zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
 487                         zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
 488                 }
 489         }
 490         pCd(z) = zdotc;
 491 }
 492 #else
 493         _Complex double zdotc = 0.0;
 494         if (incx == 1 && incy == 1) {
 495                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 496                         zdotc += Cd(&x[i]) * Cd(&y[i]);
 497                 }
 498         } else {
 499                 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
 500                         zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
 501                 }
 502         }
 503         pCd(z) = zdotc;
 504 }
 505 #endif
 506 /*  -- translated by f2c (version 20000121).
 507    You must link the resulting object file with the libraries:
 508         -lf2c -lm   (in that order)
 509 */
 510
 511
 512
 513
 514 /* Table of constant values */
 515
 516 static integer c__1 = 1;
 517 static integer c_n1 = -1;
 518
 519 /* > \brief \b IPARAM2STAGE */
 520
 521 /*  =========== DOCUMENTATION =========== */
 522
 523 /* Online html documentation available at */
 524 /*            http://www.netlib.org/lapack/explore-html/ */
 525
 526 /* > \htmlonly */
 527 /* > Download IPARAM2STAGE + dependencies */
 528 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/iparam2
 529 stage.F"> */
 530 /* > [TGZ]</a> */
 531 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/iparam2
 532 stage.F"> */
 533 /* > [ZIP]</a> */
 534 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/iparam2
 535 stage.F"> */
 536 /* > [TXT]</a> */
 537 /* > \endhtmlonly */
 538
 539 /*  Definition: */
 540 /*  =========== */
 541
 542 /*       INTEGER FUNCTION IPARAM2STAGE( ISPEC, NAME, OPTS, */
 543 /*                                    NI, NBI, IBI, NXI ) */
 544 /*       #if defined(_OPENMP) */
 545 /*           use omp_lib */
 546 /*       #endif */
 547 /*       IMPLICIT NONE */
 548
 549 /*       CHARACTER*( * )    NAME, OPTS */
 550 /*       INTEGER            ISPEC, NI, NBI, IBI, NXI */
 551
 552 /* > \par Purpose: */
 553 /*  ============= */
 554 /* > */
 555 /* > \verbatim */
 556 /* > */
 557 /* >      This program sets problem and machine dependent parameters */
 558 /* >      useful for xHETRD_2STAGE, xHETRD_HE2HB, xHETRD_HB2ST, */
 559 /* >      xGEBRD_2STAGE, xGEBRD_GE2GB, xGEBRD_GB2BD */
 560 /* >      and related subroutines for eigenvalue problems. */
 561 /* >      It is called whenever ILAENV is called with 17 <= ISPEC <= 21. */
 562 /* >      It is called whenever ILAENV2STAGE is called with 1 <= ISPEC <= 5 */
 563 /* >      with a direct conversion ISPEC + 16. */
 564 /* > \endverbatim */
 565
 566 /*  Arguments: */
 567 /*  ========== */
 568
 569 /* > \param[in] ISPEC */
 570 /* > \verbatim */
 571 /* >          ISPEC is integer scalar */
 572 /* >              ISPEC specifies which tunable parameter IPARAM2STAGE should */
 573 /* >              return. */
 574 /* > */
 575 /* >              ISPEC=17: the optimal blocksize nb for the reduction to */
 576 /* >                        BAND */
 577 /* > */
 578 /* >              ISPEC=18: the optimal blocksize ib for the eigenvectors */
 579 /* >                        singular vectors update routine */
 580 /* > */
 581 /* >              ISPEC=19: The length of the array that store the Housholder */
 582 /* >                        representation for the second stage */
 583 /* >                        Band to Tridiagonal or Bidiagonal */
 584 /* > */
 585 /* >              ISPEC=20: The workspace needed for the routine in input. */
 586 /* > */
 587 /* >              ISPEC=21: For future release. */
 588 /* > \endverbatim */
 589 /* > */
 590 /* > \param[in] NAME */
 591 /* > \verbatim */
 592 /* >          NAME is character string */
 593 /* >               Name of the calling subroutine */
 594 /* > \endverbatim */
 595 /* > */
 596 /* > \param[in] OPTS */
 597 /* > \verbatim */
 598 /* >          OPTS is CHARACTER*(*) */
 599 /* >          The character options to the subroutine NAME, concatenated */
 600 /* >          into a single character string.  For example, UPLO = 'U', */
 601 /* >          TRANS = 'T', and DIAG = 'N' for a triangular routine would */
 602 /* >          be specified as OPTS = 'UTN'. */
 603 /* > \endverbatim */
 604 /* > */
 605 /* > \param[in] NI */
 606 /* > \verbatim */
 607 /* >          NI is INTEGER which is the size of the matrix */
 608 /* > \endverbatim */
 609 /* > */
 610 /* > \param[in] NBI */
 611 /* > \verbatim */
 612 /* >          NBI is INTEGER which is the used in the reduciton, */
 613 /* >          (e.g., the size of the band), needed to compute workspace */
 614 /* >          and LHOUS2. */
 615 /* > \endverbatim */
 616 /* > */
 617 /* > \param[in] IBI */
 618 /* > \verbatim */
 619 /* >          IBI is INTEGER which represent the IB of the reduciton, */
 620 /* >          needed to compute workspace and LHOUS2. */
 621 /* > \endverbatim */
 622 /* > */
 623 /* > \param[in] NXI */
 624 /* > \verbatim */
 625 /* >          NXI is INTEGER needed in the future release. */
 626 /* > \endverbatim */
 627
 628 /*  Authors: */
 629 /*  ======== */
 630
 631 /* > \author Univ. of Tennessee */
 632 /* > \author Univ. of California Berkeley */
 633 /* > \author Univ. of Colorado Denver */
 634 /* > \author NAG Ltd. */
 635
 636 /* > \date June 2016 */
 637
 638 /* > \ingroup auxOTHERauxiliary */
 639
 640 /* > \par Further Details: */
 641 /*  ===================== */
 642 /* > */
 643 /* > \verbatim */
 644 /* > */
 645 /* >  Implemented by Azzam Haidar. */
 646 /* > */
 647 /* >  All detail are available on technical report, SC11, SC13 papers. */
 648 /* > */
 649 /* >  Azzam Haidar, Hatem Ltaief, and Jack Dongarra. */
 650 /* >  Parallel reduction to condensed forms for symmetric eigenvalue problems */
 651 /* >  using aggregated fine-grained and memory-aware kernels. In Proceedings */
 652 /* >  of 2011 International Conference for High Performance Computing, */
 653 /* >  Networking, Storage and Analysis (SC '11), New York, NY, USA, */
 654 /* >  Article 8 , 11 pages. */
 655 /* >  http://doi.acm.org/10.1145/2063384.2063394 */
 656 /* > */
 657 /* >  A. Haidar, J. Kurzak, P. Luszczek, 2013. */
 658 /* >  An improved parallel singular value algorithm and its implementation */
 659 /* >  for multicore hardware, In Proceedings of 2013 International Conference */
 660 /* >  for High Performance Computing, Networking, Storage and Analysis (SC '13). */
 661 /* >  Denver, Colorado, USA, 2013. */
 662 /* >  Article 90, 12 pages. */
 663 /* >  http://doi.acm.org/10.1145/2503210.2503292 */
 664 /* > */
 665 /* >  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra. */
 666 /* >  A novel hybrid CPU-GPU generalized eigensolver for electronic structure */
 667 /* >  calculations based on fine-grained memory aware tasks. */
 668 /* >  International Journal of High Performance Computing Applications. */
 669 /* >  Volume 28 Issue 2, Pages 196-209, May 2014. */
 670 /* >  http://hpc.sagepub.com/content/28/2/196 */
 671 /* > */
 672 /* > \endverbatim */
 673 /* > */
 674 /*  ===================================================================== */
 675 integer iparam2stage_(integer *ispec, char *name__, char *opts, integer *ni,
 676         integer *nbi, integer *ibi, integer *nxi)
 677 {
 678     /* System generated locals */
 679     integer ret_val, i__1, i__2, i__3;
 680
 681     /* Local variables */
 682     char algo[4], prec[1], stag[6], vect[1];
 683     integer nthreads, i__;
 684     logical cprec, rprec;
 685     integer lhous, lwork, factoptnb, ib, ic, kd, iz;
 686     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
 687             integer *, integer *, ftnlen, ftnlen);
 688     char subnam[14];
 689     integer lqoptnb, qroptnb;
 690     integer name_len;
 691
 692
 693 /*  -- LAPACK auxiliary routine (version 3.8.0) -- */
 694 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
 695 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
 696 /*     June 2016 */
 697
 698
 699 /*  ================================================================ */
 700
 701 /*     Invalid value for ISPEC */
 702
 703     if (*ispec < 17 || *ispec > 21) {
 704         ret_val = -1;
 705         return ret_val;
 706     }
 707
 708 /*     Get the number of threads */
 709
 710     nthreads = 1;
 711 /*      WRITE(*,*) 'IPARAM VOICI NTHREADS ISPEC ',NTHREADS, ISPEC */
 712
 713     if (*ispec != 19) {
 714
 715 /*        Convert NAME to upper case if the first character is lower case. */
 716
 717         ret_val = -1;
 718
 719 //      s_copy(subnam, name__, (ftnlen)12, name_len);
 720 strncpy(subnam,name__,13);
 721 subnam[13]='\0';
 722 for (int i=0;i<13;i++) subnam[i]=toupper(subnam[i]);
 723         //fprintf(stderr,"iparam2stage, name__ gelesen #%s#\n",name__);
 724 //fprintf(stderr,"iparam2stage, subnam gelesen #%s#\n",subnam);
 725
 726 #if 0
 727
 728         ic = *(unsigned char *)subnam;
 729         iz = 'Z';
 730         if (iz == 90 || iz == 122) {
 731
 732 /*           ASCII character set */
 733
 734             if (ic >= 97 && ic <= 122) {
 735                 *(unsigned char *)subnam = (char) (ic - 32);
 736                 for (i__ = 2; i__ <= 12; ++i__) {
 737                     ic = *(unsigned char *)&subnam[i__ - 1];
 738                     if (ic >= 97 && ic <= 122) {
 739                         *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
 740                     }
 741 /* L100: */
 742                 }
 743             }
 744
 745         } else if (iz == 233 || iz == 169) {
 746
 747 /*           EBCDIC character set */
 748
 749             if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162
 750                     && ic <= 169) {
 751                 *(unsigned char *)subnam = (char) (ic + 64);
 752                 for (i__ = 2; i__ <= 12; ++i__) {
 753                     ic = *(unsigned char *)&subnam[i__ - 1];
 754                     if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 ||
 755                             ic >= 162 && ic <= 169) {
 756                         *(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64);
 757                     }
 758 /* L110: */
 759                 }
 760             }
 761
 762         } else if (iz == 218 || iz == 250) {
 763
 764 /*           Prime machines:  ASCII+128 */
 765
 766             if (ic >= 225 && ic <= 250) {
 767                 *(unsigned char *)subnam = (char) (ic - 32);
 768                 for (i__ = 2; i__ <= 12; ++i__) {
 769                     ic = *(unsigned char *)&subnam[i__ - 1];
 770                     if (ic >= 225 && ic <= 250) {
 771                         *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
 772                     }
 773 /* L120: */
 774                 }
 775             }
 776         }
 777 #endif
 778
 779 //fprintf(stderr,"iparam2stage, subnam gross #%s#\n",subnam);
 780
 781 //      *(unsigned char *)prec = *(unsigned char *)subnam;
 782         strncpy(prec,subnam,1);
 783         strncpy(algo, subnam+3,3);
 784         algo[3]='\0';
 785         strncpy(stag, subnam+7,5);
 786         stag[5]='\0';
 787 //      s_copy(algo, subnam + 3, (ftnlen)3, (ftnlen)3);
 788 //      s_copy(stag, subnam + 7, (ftnlen)5, (ftnlen)5);
 789         rprec = *(unsigned char *)prec == 'S' || *(unsigned char *)prec ==
 790                 'D';
 791         cprec = *(unsigned char *)prec == 'C' || *(unsigned char *)prec ==
 792                 'Z';
 793
 794 /*        Invalid value for PRECISION */
 795 //fprintf(stderr," prec %s algo %s stag %s\n",prec,algo,stag);
 796         if (! (rprec || cprec)) {
 797             ret_val = -1;
 798             return ret_val;
 799         }
 800     }
 801 /*      WRITE(*,*),'RPREC,CPREC ',RPREC,CPREC, */
 802 /*     $           '   ALGO ',ALGO,'    STAGE ',STAG */
 803
 804
 805     if (*ispec == 17 || *ispec == 18) {
 806 //fprintf(stderr,"iparam2stage spec 17/18");
 807 /*     ISPEC = 17, 18:  block size KD, IB */
 808 /*     Could be also dependent from N but for now it */
 809 /*     depend only on sequential or parallel */
 810
 811         if (nthreads > 4) {
 812             if (cprec) {
 813                 kd = 128;
 814                 ib = 32;
 815             } else {
 816                 kd = 160;
 817                 ib = 40;
 818             }
 819         } else if (nthreads > 1) {
 820             if (cprec) {
 821                 kd = 64;
 822                 ib = 32;
 823             } else {
 824                 kd = 64;
 825                 ib = 32;
 826             }
 827         } else {
 828             if (cprec) {
 829                 kd = 16;
 830                 ib = 16;
 831             } else {
 832                 kd = 32;
 833                 ib = 16;
 834             }
 835         }
 836         if (*ispec == 17) {
 837             ret_val = kd;
 838         }
 839         if (*ispec == 18) {
 840             ret_val = ib;
 841         }
 842
 843     } else if (*ispec == 19) {
 844 //fprintf(stderr,"iparam2stage spec 19\n");
 845 /*     ISPEC = 19: */
 846 /*     LHOUS length of the Houselholder representation */
 847 /*     matrix (V,T) of the second stage. should be >= 1. */
 848
 849 /*     Will add the VECT OPTION HERE next release */
 850         *(unsigned char *)vect = *(unsigned char *)opts;
 851         if (*(unsigned char *)vect == 'N') {
 852 /* Computing MAX */
 853             i__1 = 1, i__2 = *ni << 2;
 854             lhous = f2cmax(i__1,i__2);
 855         } else {
 856 /*           This is not correct, it need to call the ALGO and the stage2 */
 857 /* Computing MAX */
 858             i__1 = 1, i__2 = *ni << 2;
 859             lhous = f2cmax(i__1,i__2) + *ibi;
 860         }
 861         if (lhous >= 0) {
 862             ret_val = lhous;
 863         } else {
 864             ret_val = -1;
 865         }
 866
 867     } else if (*ispec == 20) {
 868 //fprintf(stderr,"iparam2stage spec 20\n");
 869 /*     ISPEC = 20: (21 for future use) */
 870 /*     LWORK length of the workspace for */
 871 /*     either or both stages for TRD and BRD. should be >= 1. */
 872 /*     TRD: */
 873 /*     TRD_stage 1: = LT + LW + LS1 + LS2 */
 874 /*                  = LDT*KD + N*KD + N*MAX(KD,FACTOPTNB) + LDS2*KD */
 875 /*                    where LDT=LDS2=KD */
 876 /*                  = N*KD + N*f2cmax(KD,FACTOPTNB) + 2*KD*KD */
 877 /*     TRD_stage 2: = (2NB+1)*N + KD*NTHREADS */
 878 /*     TRD_both   : = f2cmax(stage1,stage2) + AB ( AB=(KD+1)*N ) */
 879 /*                  = N*KD + N*f2cmax(KD+1,FACTOPTNB) */
 880 /*                    + f2cmax(2*KD*KD, KD*NTHREADS) */
 881 /*                    + (KD+1)*N */
 882         lwork = -1;
 883         char *subnam=malloc(7*sizeof(char));
 884                         strncpy(subnam,prec,1);
 885                         sprintf(subnam+1,"GEQRF\0");
 886 //      *(unsigned char *)subnam = *(unsigned char *)prec;
 887 //      s_copy(subnam + 1, "GEQRF", (ftnlen)5, (ftnlen)5);
 888         qroptnb = ilaenv_(&c__1, subnam, " ", ni, nbi, &c_n1, &c_n1, (ftnlen)
 889                 12, (ftnlen)1);
 890         sprintf(subnam+1,"GELQF\0");
 891         s_copy(subnam + 1, "GELQF", (ftnlen)5, (ftnlen)5);
 892         lqoptnb = ilaenv_(&c__1, subnam, " ", nbi, ni, &c_n1, &c_n1, (ftnlen)
 893                 12, (ftnlen)1);
 894 /*        Could be QR or LQ for TRD and the f2cmax for BRD */
 895         factoptnb = f2cmax(qroptnb,lqoptnb);
 896         if (s_cmp(algo, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
 897             if (s_cmp(stag, "2STAG", (ftnlen)5, (ftnlen)5) == 0) {
 898 /* Computing MAX */
 899                 i__1 = *nbi + 1;
 900 /* Computing MAX */
 901                 i__2 = (*nbi << 1) * *nbi, i__3 = *nbi * nthreads;
 902                 lwork = *ni * *nbi + *ni * f2cmax(i__1,factoptnb) + f2cmax(i__2,
 903                         i__3) + (*nbi + 1) * *ni;
 904             } else if (s_cmp(stag, "HE2HB", (ftnlen)5, (ftnlen)5) == 0 ||
 905                     s_cmp(stag, "SY2SB", (ftnlen)5, (ftnlen)5) == 0) {
 906                 lwork = *ni * *nbi + *ni * f2cmax(*nbi,factoptnb) + (*nbi << 1) *
 907                          *nbi;
 908             } else if (s_cmp(stag, "HB2ST", (ftnlen)5, (ftnlen)5) == 0 ||
 909                     s_cmp(stag, "SB2ST", (ftnlen)5, (ftnlen)5) == 0) {
 910                 lwork = ((*nbi << 1) + 1) * *ni + *nbi * nthreads;
 911             }
 912         } else if (s_cmp(algo, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
 913             if (s_cmp(stag, "2STAG", (ftnlen)5, (ftnlen)5) == 0) {
 914 /* Computing MAX */
 915                 i__1 = *nbi + 1;
 916 /* Computing MAX */
 917                 i__2 = (*nbi << 1) * *nbi, i__3 = *nbi * nthreads;
 918                 lwork = (*ni << 1) * *nbi + *ni * f2cmax(i__1,factoptnb) + f2cmax(
 919                         i__2,i__3) + (*nbi + 1) * *ni;
 920             } else if (s_cmp(stag, "GE2GB", (ftnlen)5, (ftnlen)5) == 0) {
 921                 lwork = *ni * *nbi + *ni * f2cmax(*nbi,factoptnb) + (*nbi << 1) *
 922                          *nbi;
 923             } else if (s_cmp(stag, "GB2BD", (ftnlen)5, (ftnlen)5) == 0) {
 924                 lwork = (*nbi * 3 + 1) * *ni + *nbi * nthreads;
 925             }
 926         }
 927         lwork = f2cmax(1,lwork);
 928         if (lwork > 0) {
 929             ret_val = lwork;
 930         } else {
 931             ret_val = -1;
 932         }
 933
 934     } else if (*ispec == 21) {
 935 //fprintf(stderr,"iparam2stage spec 21\n");
 936 /*     ISPEC = 21 for future use */
 937         ret_val = *nxi;
 938     }
 939
 940 /*     ==== End of IPARAM2STAGE ==== */
 941
 942     return ret_val;
 943 } /* iparam2stage_ */
 944