lapack-netlib/SRC/cgejsv.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
 514 /* Table of constant values */
 515
 516 static complex c_b1 = {0.f,0.f};
 517 static complex c_b2 = {1.f,0.f};
 518 static integer c_n1 = -1;
 519 static integer c__1 = 1;
 520 static integer c__0 = 0;
 521 static real c_b141 = 1.f;
 522 static logical c_false = FALSE_;
 523
 524 /* > \brief \b CGEJSV */
 525
 526 /*  =========== DOCUMENTATION =========== */
 527
 528 /* Online html documentation available at */
 529 /*            http://www.netlib.org/lapack/explore-html/ */
 530
 531 /* > \htmlonly */
 532 /* > Download CGEJSV + dependencies */
 533 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgejsv.
 534 f"> */
 535 /* > [TGZ]</a> */
 536 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgejsv.
 537 f"> */
 538 /* > [ZIP]</a> */
 539 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgejsv.
 540 f"> */
 541 /* > [TXT]</a> */
 542 /* > \endhtmlonly */
 543
 544 /*  Definition: */
 545 /*  =========== */
 546
 547 /*     SUBROUTINE CGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, */
 548 /*                         M, N, A, LDA, SVA, U, LDU, V, LDV, */
 549 /*                         CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) */
 550
 551 /*     IMPLICIT    NONE */
 552 /*     INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N */
 553 /*     COMPLEX     A( LDA, * ),  U( LDU, * ), V( LDV, * ), CWORK( LWORK ) */
 554 /*     REAL        SVA( N ), RWORK( LRWORK ) */
 555 /*     INTEGER     IWORK( * ) */
 556 /*     CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV */
 557
 558
 559 /* > \par Purpose: */
 560 /*  ============= */
 561 /* > */
 562 /* > \verbatim */
 563 /* > */
 564 /* > CGEJSV computes the singular value decomposition (SVD) of a complex M-by-N */
 565 /* > matrix [A], where M >= N. The SVD of [A] is written as */
 566 /* > */
 567 /* >              [A] = [U] * [SIGMA] * [V]^*, */
 568 /* > */
 569 /* > where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N */
 570 /* > diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and */
 571 /* > [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are */
 572 /* > the singular values of [A]. The columns of [U] and [V] are the left and */
 573 /* > the right singular vectors of [A], respectively. The matrices [U] and [V] */
 574 /* > are computed and stored in the arrays U and V, respectively. The diagonal */
 575 /* > of [SIGMA] is computed and stored in the array SVA. */
 576 /* > \endverbatim */
 577 /* > */
 578 /* >  Arguments: */
 579 /* >  ========== */
 580 /* > */
 581 /* > \param[in] JOBA */
 582 /* > \verbatim */
 583 /* >          JOBA is CHARACTER*1 */
 584 /* >         Specifies the level of accuracy: */
 585 /* >       = 'C': This option works well (high relative accuracy) if A = B * D, */
 586 /* >              with well-conditioned B and arbitrary diagonal matrix D. */
 587 /* >              The accuracy cannot be spoiled by COLUMN scaling. The */
 588 /* >              accuracy of the computed output depends on the condition of */
 589 /* >              B, and the procedure aims at the best theoretical accuracy. */
 590 /* >              The relative error max_{i=1:N}|d sigma_i| / sigma_i is */
 591 /* >              bounded by f(M,N)*epsilon* cond(B), independent of D. */
 592 /* >              The input matrix is preprocessed with the QRF with column */
 593 /* >              pivoting. This initial preprocessing and preconditioning by */
 594 /* >              a rank revealing QR factorization is common for all values of */
 595 /* >              JOBA. Additional actions are specified as follows: */
 596 /* >       = 'E': Computation as with 'C' with an additional estimate of the */
 597 /* >              condition number of B. It provides a realistic error bound. */
 598 /* >       = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings */
 599 /* >              D1, D2, and well-conditioned matrix C, this option gives */
 600 /* >              higher accuracy than the 'C' option. If the structure of the */
 601 /* >              input matrix is not known, and relative accuracy is */
 602 /* >              desirable, then this option is advisable. The input matrix A */
 603 /* >              is preprocessed with QR factorization with FULL (row and */
 604 /* >              column) pivoting. */
 605 /* >       = 'G': Computation as with 'F' with an additional estimate of the */
 606 /* >              condition number of B, where A=B*D. If A has heavily weighted */
 607 /* >              rows, then using this condition number gives too pessimistic */
 608 /* >              error bound. */
 609 /* >       = 'A': Small singular values are not well determined by the data */
 610 /* >              and are considered as noisy; the matrix is treated as */
 611 /* >              numerically rank deficient. The error in the computed */
 612 /* >              singular values is bounded by f(m,n)*epsilon*||A||. */
 613 /* >              The computed SVD A = U * S * V^* restores A up to */
 614 /* >              f(m,n)*epsilon*||A||. */
 615 /* >              This gives the procedure the licence to discard (set to zero) */
 616 /* >              all singular values below N*epsilon*||A||. */
 617 /* >       = 'R': Similar as in 'A'. Rank revealing property of the initial */
 618 /* >              QR factorization is used do reveal (using triangular factor) */
 619 /* >              a gap sigma_{r+1} < epsilon * sigma_r in which case the */
 620 /* >              numerical RANK is declared to be r. The SVD is computed with */
 621 /* >              absolute error bounds, but more accurately than with 'A'. */
 622 /* > \endverbatim */
 623 /* > */
 624 /* > \param[in] JOBU */
 625 /* > \verbatim */
 626 /* >          JOBU is CHARACTER*1 */
 627 /* >         Specifies whether to compute the columns of U: */
 628 /* >       = 'U': N columns of U are returned in the array U. */
 629 /* >       = 'F': full set of M left sing. vectors is returned in the array U. */
 630 /* >       = 'W': U may be used as workspace of length M*N. See the description */
 631 /* >              of U. */
 632 /* >       = 'N': U is not computed. */
 633 /* > \endverbatim */
 634 /* > */
 635 /* > \param[in] JOBV */
 636 /* > \verbatim */
 637 /* >          JOBV is CHARACTER*1 */
 638 /* >         Specifies whether to compute the matrix V: */
 639 /* >       = 'V': N columns of V are returned in the array V; Jacobi rotations */
 640 /* >              are not explicitly accumulated. */
 641 /* >       = 'J': N columns of V are returned in the array V, but they are */
 642 /* >              computed as the product of Jacobi rotations, if JOBT = 'N'. */
 643 /* >       = 'W': V may be used as workspace of length N*N. See the description */
 644 /* >              of V. */
 645 /* >       = 'N': V is not computed. */
 646 /* > \endverbatim */
 647 /* > */
 648 /* > \param[in] JOBR */
 649 /* > \verbatim */
 650 /* >          JOBR is CHARACTER*1 */
 651 /* >         Specifies the RANGE for the singular values. Issues the licence to */
 652 /* >         set to zero small positive singular values if they are outside */
 653 /* >         specified range. If A .NE. 0 is scaled so that the largest singular */
 654 /* >         value of c*A is around SQRT(BIG), BIG=SLAMCH('O'), then JOBR issues */
 655 /* >         the licence to kill columns of A whose norm in c*A is less than */
 656 /* >         SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN, */
 657 /* >         where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). */
 658 /* >       = 'N': Do not kill small columns of c*A. This option assumes that */
 659 /* >              BLAS and QR factorizations and triangular solvers are */
 660 /* >              implemented to work in that range. If the condition of A */
 661 /* >              is greater than BIG, use CGESVJ. */
 662 /* >       = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)] */
 663 /* >              (roughly, as described above). This option is recommended. */
 664 /* >                                             =========================== */
 665 /* >         For computing the singular values in the FULL range [SFMIN,BIG] */
 666 /* >         use CGESVJ. */
 667 /* > \endverbatim */
 668 /* > */
 669 /* > \param[in] JOBT */
 670 /* > \verbatim */
 671 /* >          JOBT is CHARACTER*1 */
 672 /* >         If the matrix is square then the procedure may determine to use */
 673 /* >         transposed A if A^* seems to be better with respect to convergence. */
 674 /* >         If the matrix is not square, JOBT is ignored. */
 675 /* >         The decision is based on two values of entropy over the adjoint */
 676 /* >         orbit of A^* * A. See the descriptions of WORK(6) and WORK(7). */
 677 /* >       = 'T': transpose if entropy test indicates possibly faster */
 678 /* >         convergence of Jacobi process if A^* is taken as input. If A is */
 679 /* >         replaced with A^*, then the row pivoting is included automatically. */
 680 /* >       = 'N': do not speculate. */
 681 /* >         The option 'T' can be used to compute only the singular values, or */
 682 /* >         the full SVD (U, SIGMA and V). For only one set of singular vectors */
 683 /* >         (U or V), the caller should provide both U and V, as one of the */
 684 /* >         matrices is used as workspace if the matrix A is transposed. */
 685 /* >         The implementer can easily remove this constraint and make the */
 686 /* >         code more complicated. See the descriptions of U and V. */
 687 /* >         In general, this option is considered experimental, and 'N'; should */
 688 /* >         be preferred. This is subject to changes in the future. */
 689 /* > \endverbatim */
 690 /* > */
 691 /* > \param[in] JOBP */
 692 /* > \verbatim */
 693 /* >          JOBP is CHARACTER*1 */
 694 /* >         Issues the licence to introduce structured perturbations to drown */
 695 /* >         denormalized numbers. This licence should be active if the */
 696 /* >         denormals are poorly implemented, causing slow computation, */
 697 /* >         especially in cases of fast convergence (!). For details see [1,2]. */
 698 /* >         For the sake of simplicity, this perturbations are included only */
 699 /* >         when the full SVD or only the singular values are requested. The */
 700 /* >         implementer/user can easily add the perturbation for the cases of */
 701 /* >         computing one set of singular vectors. */
 702 /* >       = 'P': introduce perturbation */
 703 /* >       = 'N': do not perturb */
 704 /* > \endverbatim */
 705 /* > */
 706 /* > \param[in] M */
 707 /* > \verbatim */
 708 /* >          M is INTEGER */
 709 /* >         The number of rows of the input matrix A.  M >= 0. */
 710 /* > \endverbatim */
 711 /* > */
 712 /* > \param[in] N */
 713 /* > \verbatim */
 714 /* >          N is INTEGER */
 715 /* >         The number of columns of the input matrix A. M >= N >= 0. */
 716 /* > \endverbatim */
 717 /* > */
 718 /* > \param[in,out] A */
 719 /* > \verbatim */
 720 /* >          A is COMPLEX array, dimension (LDA,N) */
 721 /* >          On entry, the M-by-N matrix A. */
 722 /* > \endverbatim */
 723 /* > */
 724 /* > \param[in] LDA */
 725 /* > \verbatim */
 726 /* >          LDA is INTEGER */
 727 /* >          The leading dimension of the array A.  LDA >= f2cmax(1,M). */
 728 /* > \endverbatim */
 729 /* > */
 730 /* > \param[out] SVA */
 731 /* > \verbatim */
 732 /* >          SVA is REAL array, dimension (N) */
 733 /* >          On exit, */
 734 /* >          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the */
 735 /* >            computation SVA contains Euclidean column norms of the */
 736 /* >            iterated matrices in the array A. */
 737 /* >          - For WORK(1) .NE. WORK(2): The singular values of A are */
 738 /* >            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if */
 739 /* >            sigma_max(A) overflows or if small singular values have been */
 740 /* >            saved from underflow by scaling the input matrix A. */
 741 /* >          - If JOBR='R' then some of the singular values may be returned */
 742 /* >            as exact zeros obtained by "set to zero" because they are */
 743 /* >            below the numerical rank threshold or are denormalized numbers. */
 744 /* > \endverbatim */
 745 /* > */
 746 /* > \param[out] U */
 747 /* > \verbatim */
 748 /* >          U is COMPLEX array, dimension ( LDU, N ) or ( LDU, M ) */
 749 /* >          If JOBU = 'U', then U contains on exit the M-by-N matrix of */
 750 /* >                         the left singular vectors. */
 751 /* >          If JOBU = 'F', then U contains on exit the M-by-M matrix of */
 752 /* >                         the left singular vectors, including an ONB */
 753 /* >                         of the orthogonal complement of the Range(A). */
 754 /* >          If JOBU = 'W'  .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N), */
 755 /* >                         then U is used as workspace if the procedure */
 756 /* >                         replaces A with A^*. In that case, [V] is computed */
 757 /* >                         in U as left singular vectors of A^* and then */
 758 /* >                         copied back to the V array. This 'W' option is just */
 759 /* >                         a reminder to the caller that in this case U is */
 760 /* >                         reserved as workspace of length N*N. */
 761 /* >          If JOBU = 'N'  U is not referenced, unless JOBT='T'. */
 762 /* > \endverbatim */
 763 /* > */
 764 /* > \param[in] LDU */
 765 /* > \verbatim */
 766 /* >          LDU is INTEGER */
 767 /* >          The leading dimension of the array U,  LDU >= 1. */
 768 /* >          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M. */
 769 /* > \endverbatim */
 770 /* > */
 771 /* > \param[out] V */
 772 /* > \verbatim */
 773 /* >          V is COMPLEX array, dimension ( LDV, N ) */
 774 /* >          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of */
 775 /* >                         the right singular vectors; */
 776 /* >          If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), */
 777 /* >                         then V is used as workspace if the pprocedure */
 778 /* >                         replaces A with A^*. In that case, [U] is computed */
 779 /* >                         in V as right singular vectors of A^* and then */
 780 /* >                         copied back to the U array. This 'W' option is just */
 781 /* >                         a reminder to the caller that in this case V is */
 782 /* >                         reserved as workspace of length N*N. */
 783 /* >          If JOBV = 'N'  V is not referenced, unless JOBT='T'. */
 784 /* > \endverbatim */
 785 /* > */
 786 /* > \param[in] LDV */
 787 /* > \verbatim */
 788 /* >          LDV is INTEGER */
 789 /* >          The leading dimension of the array V,  LDV >= 1. */
 790 /* >          If JOBV = 'V' or 'J' or 'W', then LDV >= N. */
 791 /* > \endverbatim */
 792 /* > */
 793 /* > \param[out] CWORK */
 794 /* > \verbatim */
 795 /* >          CWORK is COMPLEX array, dimension (MAX(2,LWORK)) */
 796 /* >          If the call to CGEJSV is a workspace query (indicated by LWORK=-1 or */
 797 /* >          LRWORK=-1), then on exit CWORK(1) contains the required length of */
 798 /* >          CWORK for the job parameters used in the call. */
 799 /* > \endverbatim */
 800 /* > */
 801 /* > \param[in] LWORK */
 802 /* > \verbatim */
 803 /* >          LWORK is INTEGER */
 804 /* >          Length of CWORK to confirm proper allocation of workspace. */
 805 /* >          LWORK depends on the job: */
 806 /* > */
 807 /* >          1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and */
 808 /* >            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): */
 809 /* >               LWORK >= 2*N+1. This is the minimal requirement. */
 810 /* >               ->> For optimal performance (blocked code) the optimal value */
 811 /* >               is LWORK >= N + (N+1)*NB. Here NB is the optimal */
 812 /* >               block size for CGEQP3 and CGEQRF. */
 813 /* >               In general, optimal LWORK is computed as */
 814 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3),N+LWORK(CGEQRF), LWORK(CGESVJ)). */
 815 /* >            1.2. .. an estimate of the scaled condition number of A is */
 816 /* >               required (JOBA='E', or 'G'). In this case, LWORK the minimal */
 817 /* >               requirement is LWORK >= N*N + 2*N. */
 818 /* >               ->> For optimal performance (blocked code) the optimal value */
 819 /* >               is LWORK >= f2cmax(N+(N+1)*NB, N*N+2*N)=N**2+2*N. */
 820 /* >               In general, the optimal length LWORK is computed as */
 821 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3),N+LWORK(CGEQRF), LWORK(CGESVJ), */
 822 /* >                            N*N+LWORK(CPOCON)). */
 823 /* >          2. If SIGMA and the right singular vectors are needed (JOBV = 'V'), */
 824 /* >             (JOBU = 'N') */
 825 /* >            2.1   .. no scaled condition estimate requested (JOBE = 'N'): */
 826 /* >            -> the minimal requirement is LWORK >= 3*N. */
 827 /* >            -> For optimal performance, */
 828 /* >               LWORK >= f2cmax(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, */
 829 /* >               where NB is the optimal block size for CGEQP3, CGEQRF, CGELQ, */
 830 /* >               CUNMLQ. In general, the optimal length LWORK is computed as */
 831 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3), N+LWORK(CGESVJ), */
 832 /* >                       N+LWORK(CGELQF), 2*N+LWORK(CGEQRF), N+LWORK(CUNMLQ)). */
 833 /* >            2.2 .. an estimate of the scaled condition number of A is */
 834 /* >               required (JOBA='E', or 'G'). */
 835 /* >            -> the minimal requirement is LWORK >= 3*N. */
 836 /* >            -> For optimal performance, */
 837 /* >               LWORK >= f2cmax(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB, */
 838 /* >               where NB is the optimal block size for CGEQP3, CGEQRF, CGELQ, */
 839 /* >               CUNMLQ. In general, the optimal length LWORK is computed as */
 840 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3), LWORK(CPOCON), N+LWORK(CGESVJ), */
 841 /* >                       N+LWORK(CGELQF), 2*N+LWORK(CGEQRF), N+LWORK(CUNMLQ)). */
 842 /* >          3. If SIGMA and the left singular vectors are needed */
 843 /* >            3.1  .. no scaled condition estimate requested (JOBE = 'N'): */
 844 /* >            -> the minimal requirement is LWORK >= 3*N. */
 845 /* >            -> For optimal performance: */
 846 /* >               if JOBU = 'U' :: LWORK >= f2cmax(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, */
 847 /* >               where NB is the optimal block size for CGEQP3, CGEQRF, CUNMQR. */
 848 /* >               In general, the optimal length LWORK is computed as */
 849 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3), 2*N+LWORK(CGEQRF), N+LWORK(CUNMQR)). */
 850 /* >            3.2  .. an estimate of the scaled condition number of A is */
 851 /* >               required (JOBA='E', or 'G'). */
 852 /* >            -> the minimal requirement is LWORK >= 3*N. */
 853 /* >            -> For optimal performance: */
 854 /* >               if JOBU = 'U' :: LWORK >= f2cmax(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, */
 855 /* >               where NB is the optimal block size for CGEQP3, CGEQRF, CUNMQR. */
 856 /* >               In general, the optimal length LWORK is computed as */
 857 /* >               LWORK >= f2cmax(N+LWORK(CGEQP3),N+LWORK(CPOCON), */
 858 /* >                        2*N+LWORK(CGEQRF), N+LWORK(CUNMQR)). */
 859 /* > */
 860 /* >          4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and */
 861 /* >            4.1. if JOBV = 'V' */
 862 /* >               the minimal requirement is LWORK >= 5*N+2*N*N. */
 863 /* >            4.2. if JOBV = 'J' the minimal requirement is */
 864 /* >               LWORK >= 4*N+N*N. */
 865 /* >            In both cases, the allocated CWORK can accommodate blocked runs */
 866 /* >            of CGEQP3, CGEQRF, CGELQF, CUNMQR, CUNMLQ. */
 867 /* > */
 868 /* >          If the call to CGEJSV is a workspace query (indicated by LWORK=-1 or */
 869 /* >          LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the */
 870 /* >          minimal length of CWORK for the job parameters used in the call. */
 871 /* > \endverbatim */
 872 /* > */
 873 /* > \param[out] RWORK */
 874 /* > \verbatim */
 875 /* >          RWORK is REAL array, dimension (MAX(7,LWORK)) */
 876 /* >          On exit, */
 877 /* >          RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1) */
 878 /* >                    such that SCALE*SVA(1:N) are the computed singular values */
 879 /* >                    of A. (See the description of SVA().) */
 880 /* >          RWORK(2) = See the description of RWORK(1). */
 881 /* >          RWORK(3) = SCONDA is an estimate for the condition number of */
 882 /* >                    column equilibrated A. (If JOBA = 'E' or 'G') */
 883 /* >                    SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). */
 884 /* >                    It is computed using SPOCON. It holds */
 885 /* >                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA */
 886 /* >                    where R is the triangular factor from the QRF of A. */
 887 /* >                    However, if R is truncated and the numerical rank is */
 888 /* >                    determined to be strictly smaller than N, SCONDA is */
 889 /* >                    returned as -1, thus indicating that the smallest */
 890 /* >                    singular values might be lost. */
 891 /* > */
 892 /* >          If full SVD is needed, the following two condition numbers are */
 893 /* >          useful for the analysis of the algorithm. They are provied for */
 894 /* >          a developer/implementer who is familiar with the details of */
 895 /* >          the method. */
 896 /* > */
 897 /* >          RWORK(4) = an estimate of the scaled condition number of the */
 898 /* >                    triangular factor in the first QR factorization. */
 899 /* >          RWORK(5) = an estimate of the scaled condition number of the */
 900 /* >                    triangular factor in the second QR factorization. */
 901 /* >          The following two parameters are computed if JOBT = 'T'. */
 902 /* >          They are provided for a developer/implementer who is familiar */
 903 /* >          with the details of the method. */
 904 /* >          RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy */
 905 /* >                    of diag(A^* * A) / Trace(A^* * A) taken as point in the */
 906 /* >                    probability simplex. */
 907 /* >          RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).) */
 908 /* >          If the call to CGEJSV is a workspace query (indicated by LWORK=-1 or */
 909 /* >          LRWORK=-1), then on exit RWORK(1) contains the required length of */
 910 /* >          RWORK for the job parameters used in the call. */
 911 /* > \endverbatim */
 912 /* > */
 913 /* > \param[in] LRWORK */
 914 /* > \verbatim */
 915 /* >          LRWORK is INTEGER */
 916 /* >          Length of RWORK to confirm proper allocation of workspace. */
 917 /* >          LRWORK depends on the job: */
 918 /* > */
 919 /* >       1. If only the singular values are requested i.e. if */
 920 /* >          LSAME(JOBU,'N') .AND. LSAME(JOBV,'N') */
 921 /* >          then: */
 922 /* >          1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), */
 923 /* >               then: LRWORK = f2cmax( 7, 2 * M ). */
 924 /* >          1.2. Otherwise, LRWORK  = f2cmax( 7,  N ). */
 925 /* >       2. If singular values with the right singular vectors are requested */
 926 /* >          i.e. if */
 927 /* >          (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND. */
 928 /* >          .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) */
 929 /* >          then: */
 930 /* >          2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), */
 931 /* >          then LRWORK = f2cmax( 7, 2 * M ). */
 932 /* >          2.2. Otherwise, LRWORK  = f2cmax( 7,  N ). */
 933 /* >       3. If singular values with the left singular vectors are requested, i.e. if */
 934 /* >          (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. */
 935 /* >          .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) */
 936 /* >          then: */
 937 /* >          3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), */
 938 /* >          then LRWORK = f2cmax( 7, 2 * M ). */
 939 /* >          3.2. Otherwise, LRWORK  = f2cmax( 7,  N ). */
 940 /* >       4. If singular values with both the left and the right singular vectors */
 941 /* >          are requested, i.e. if */
 942 /* >          (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. */
 943 /* >          (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) */
 944 /* >          then: */
 945 /* >          4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), */
 946 /* >          then LRWORK = f2cmax( 7, 2 * M ). */
 947 /* >          4.2. Otherwise, LRWORK  = f2cmax( 7, N ). */
 948 /* > */
 949 /* >          If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and */
 950 /* >          the length of RWORK is returned in RWORK(1). */
 951 /* > \endverbatim */
 952 /* > */
 953 /* > \param[out] IWORK */
 954 /* > \verbatim */
 955 /* >          IWORK is INTEGER array, of dimension at least 4, that further depends */
 956 /* >          on the job: */
 957 /* > */
 958 /* >          1. If only the singular values are requested then: */
 959 /* >             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) */
 960 /* >             then the length of IWORK is N+M; otherwise the length of IWORK is N. */
 961 /* >          2. If the singular values and the right singular vectors are requested then: */
 962 /* >             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) */
 963 /* >             then the length of IWORK is N+M; otherwise the length of IWORK is N. */
 964 /* >          3. If the singular values and the left singular vectors are requested then: */
 965 /* >             If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) */
 966 /* >             then the length of IWORK is N+M; otherwise the length of IWORK is N. */
 967 /* >          4. If the singular values with both the left and the right singular vectors */
 968 /* >             are requested, then: */
 969 /* >             4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows: */
 970 /* >                  If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) */
 971 /* >                  then the length of IWORK is N+M; otherwise the length of IWORK is N. */
 972 /* >             4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows: */
 973 /* >                  If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) */
 974 /* >                  then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N. */
 975 /* > */
 976 /* >          On exit, */
 977 /* >          IWORK(1) = the numerical rank determined after the initial */
 978 /* >                     QR factorization with pivoting. See the descriptions */
 979 /* >                     of JOBA and JOBR. */
 980 /* >          IWORK(2) = the number of the computed nonzero singular values */
 981 /* >          IWORK(3) = if nonzero, a warning message: */
 982 /* >                     If IWORK(3) = 1 then some of the column norms of A */
 983 /* >                     were denormalized floats. The requested high accuracy */
 984 /* >                     is not warranted by the data. */
 985 /* >          IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to */
 986 /* >                     do the job as specified by the JOB parameters. */
 987 /* >          If the call to CGEJSV is a workspace query (indicated by LWORK = -1 and */
 988 /* >          LRWORK = -1), then on exit IWORK(1) contains the required length of */
 989 /* >          IWORK for the job parameters used in the call. */
 990 /* > \endverbatim */
 991 /* > */
 992 /* > \param[out] INFO */
 993 /* > \verbatim */
 994 /* >          INFO is INTEGER */
 995 /* >           < 0:  if INFO = -i, then the i-th argument had an illegal value. */
 996 /* >           = 0:  successful exit; */
 997 /* >           > 0:  CGEJSV  did not converge in the maximal allowed number */
 998 /* >                 of sweeps. The computed values may be inaccurate. */
 999 /* > \endverbatim */
1000
1001 /*  Authors: */
1002 /*  ======== */
1003
1004 /* > \author Univ. of Tennessee */
1005 /* > \author Univ. of California Berkeley */
1006 /* > \author Univ. of Colorado Denver */
1007 /* > \author NAG Ltd. */
1008
1009 /* > \date June 2016 */
1010
1011 /* > \ingroup complexGEsing */
1012
1013 /* > \par Further Details: */
1014 /*  ===================== */
1015 /* > */
1016 /* > \verbatim */
1017 /* >  CGEJSV implements a preconditioned Jacobi SVD algorithm. It uses CGEQP3, */
1018 /* >  CGEQRF, and CGELQF as preprocessors and preconditioners. Optionally, an */
1019 /* >  additional row pivoting can be used as a preprocessor, which in some */
1020 /* >  cases results in much higher accuracy. An example is matrix A with the */
1021 /* >  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned */
1022 /* >  diagonal matrices and C is well-conditioned matrix. In that case, complete */
1023 /* >  pivoting in the first QR factorizations provides accuracy dependent on the */
1024 /* >  condition number of C, and independent of D1, D2. Such higher accuracy is */
1025 /* >  not completely understood theoretically, but it works well in practice. */
1026 /* >  Further, if A can be written as A = B*D, with well-conditioned B and some */
1027 /* >  diagonal D, then the high accuracy is guaranteed, both theoretically and */
1028 /* >  in software, independent of D. For more details see [1], [2]. */
1029 /* >     The computational range for the singular values can be the full range */
1030 /* >  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS */
1031 /* >  & LAPACK routines called by CGEJSV are implemented to work in that range. */
1032 /* >  If that is not the case, then the restriction for safe computation with */
1033 /* >  the singular values in the range of normalized IEEE numbers is that the */
1034 /* >  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not */
1035 /* >  overflow. This code (CGEJSV) is best used in this restricted range, */
1036 /* >  meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are */
1037 /* >  returned as zeros. See JOBR for details on this. */
1038 /* >     Further, this implementation is somewhat slower than the one described */
1039 /* >  in [1,2] due to replacement of some non-LAPACK components, and because */
1040 /* >  the choice of some tuning parameters in the iterative part (CGESVJ) is */
1041 /* >  left to the implementer on a particular machine. */
1042 /* >     The rank revealing QR factorization (in this code: CGEQP3) should be */
1043 /* >  implemented as in [3]. We have a new version of CGEQP3 under development */
1044 /* >  that is more robust than the current one in LAPACK, with a cleaner cut in */
1045 /* >  rank deficient cases. It will be available in the SIGMA library [4]. */
1046 /* >  If M is much larger than N, it is obvious that the initial QRF with */
1047 /* >  column pivoting can be preprocessed by the QRF without pivoting. That */
1048 /* >  well known trick is not used in CGEJSV because in some cases heavy row */
1049 /* >  weighting can be treated with complete pivoting. The overhead in cases */
1050 /* >  M much larger than N is then only due to pivoting, but the benefits in */
1051 /* >  terms of accuracy have prevailed. The implementer/user can incorporate */
1052 /* >  this extra QRF step easily. The implementer can also improve data movement */
1053 /* >  (matrix transpose, matrix copy, matrix transposed copy) - this */
1054 /* >  implementation of CGEJSV uses only the simplest, naive data movement. */
1055 /* > \endverbatim */
1056
1057 /* > \par Contributor: */
1058 /*  ================== */
1059 /* > */
1060 /* >  Zlatko Drmac (Zagreb, Croatia) */
1061
1062 /* > \par References: */
1063 /*  ================ */
1064 /* > */
1065 /* > \verbatim */
1066 /* > */
1067 /* > [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. */
1068 /* >     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. */
1069 /* >     LAPACK Working note 169. */
1070 /* > [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. */
1071 /* >     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. */
1072 /* >     LAPACK Working note 170. */
1073 /* > [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR */
1074 /* >     factorization software - a case study. */
1075 /* >     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. */
1076 /* >     LAPACK Working note 176. */
1077 /* > [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, */
1078 /* >     QSVD, (H,K)-SVD computations. */
1079 /* >     Department of Mathematics, University of Zagreb, 2008, 2016. */
1080 /* > \endverbatim */
1081
1082 /* >  \par Bugs, examples and comments: */
1083 /*   ================================= */
1084 /* > */
1085 /* >  Please report all bugs and send interesting examples and/or comments to */
1086 /* >  drmac@math.hr. Thank you. */
1087 /* > */
1088 /*  ===================================================================== */
1089 /* Subroutine */ int cgejsv_(char *joba, char *jobu, char *jobv, char *jobr,
1090         char *jobt, char *jobp, integer *m, integer *n, complex *a, integer *
1091         lda, real *sva, complex *u, integer *ldu, complex *v, integer *ldv,
1092         complex *cwork, integer *lwork, real *rwork, integer *lrwork, integer
1093         *iwork, integer *info)
1094 {
1095     /* System generated locals */
1096     integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
1097             i__3, i__4, i__5, i__6, i__7, i__8, i__9, i__10, i__11;
1098     real r__1, r__2, r__3;
1099     complex q__1;
1100
1101     /* Local variables */
1102     integer lwrk_cunmqr__;
1103     logical defr;
1104     real aapp, aaqq;
1105     logical kill;
1106     integer ierr, lwrk_cgeqp3n__;
1107     real temp1;
1108     integer lwunmqrm, lwrk_cgesvju__, lwrk_cgesvjv__, lwqp3, lwrk_cunmqrm__,
1109             p, q;
1110     logical jracc;
1111     extern logical lsame_(char *, char *);
1112     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
1113     complex ctemp;
1114     real entra, small;
1115     integer iwoff;
1116     real sfmin;
1117     logical lsvec;
1118     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
1119             complex *, integer *), cswap_(integer *, complex *, integer *,
1120             complex *, integer *);
1121     real epsln;
1122     logical rsvec;
1123     integer lwcon, lwlqf;
1124     extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *,
1125             integer *, integer *, complex *, complex *, integer *, complex *,
1126             integer *);
1127     integer lwqrf, n1;
1128     logical l2aber;
1129     extern /* Subroutine */ int cgeqp3_(integer *, integer *, complex *,
1130             integer *, integer *, complex *, complex *, integer *, real *,
1131             integer *);
1132     real condr1, condr2, uscal1, uscal2;
1133     logical l2kill, l2rank, l2tran;
1134     extern real scnrm2_(integer *, complex *, integer *);
1135     logical l2pert;
1136     integer lrwqp3;
1137     extern /* Subroutine */ int clacgv_(integer *, complex *, integer *);
1138     integer nr;
1139     extern /* Subroutine */ int cgelqf_(integer *, integer *, complex *,
1140             integer *, complex *, complex *, integer *, integer *);
1141     extern integer icamax_(integer *, complex *, integer *);
1142     extern /* Subroutine */ int clascl_(char *, integer *, integer *, real *,
1143             real *, integer *, integer *, complex *, integer *, integer *);
1144     real scalem, sconda;
1145     logical goscal;
1146     real aatmin;
1147     extern real slamch_(char *);
1148     real aatmax;
1149     extern /* Subroutine */ int cgeqrf_(integer *, integer *, complex *,
1150             integer *, complex *, complex *, integer *, integer *), clacpy_(
1151             char *, integer *, integer *, complex *, integer *, complex *,
1152             integer *), clapmr_(logical *, integer *, integer *,
1153             complex *, integer *, integer *);
1154     logical noscal;
1155     extern /* Subroutine */ int claset_(char *, integer *, integer *, complex
1156             *, complex *, complex *, integer *);
1157     extern integer isamax_(integer *, real *, integer *);
1158     extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
1159             real *, integer *, integer *, real *, integer *, integer *), cpocon_(char *, integer *, complex *, integer *, real *,
1160             real *, complex *, real *, integer *), csscal_(integer *,
1161             real *, complex *, integer *), classq_(integer *, complex *,
1162             integer *, real *, real *), xerbla_(char *, integer *, ftnlen),
1163             cgesvj_(char *, char *, char *, integer *, integer *, complex *,
1164             integer *, real *, integer *, complex *, integer *, complex *,
1165             integer *, real *, integer *, integer *),
1166             claswp_(integer *, complex *, integer *, integer *, integer *,
1167             integer *, integer *);
1168     real entrat;
1169     logical almort;
1170     complex cdummy[1];
1171     extern /* Subroutine */ int cungqr_(integer *, integer *, integer *,
1172             complex *, integer *, complex *, complex *, integer *, integer *);
1173     real maxprj;
1174     extern /* Subroutine */ int cunmlq_(char *, char *, integer *, integer *,
1175             integer *, complex *, integer *, complex *, complex *, integer *,
1176             complex *, integer *, integer *);
1177     logical errest;
1178     integer lrwcon;
1179     extern /* Subroutine */ int slassq_(integer *, real *, integer *, real *,
1180             real *);
1181     logical transp;
1182     integer minwrk, lwsvdj;
1183     extern /* Subroutine */ int cunmqr_(char *, char *, integer *, integer *,
1184             integer *, complex *, integer *, complex *, complex *, integer *,
1185             complex *, integer *, integer *);
1186     real rdummy[1];
1187     logical lquery, rowpiv;
1188     integer optwrk;
1189     real big;
1190     integer lwrk_cgeqp3__;
1191     real cond_ok__, xsc, big1;
1192     integer warning, numrank, lwrk_cgelqf__, miniwrk, lwrk_cgeqrf__, minrwrk,
1193             lrwsvdj, lwunmlq, lwsvdjv, lwrk_cgesvj__, lwunmqr, lwrk_cunmlq__;
1194
1195
1196 /*  -- LAPACK computational routine (version 3.7.1) -- */
1197 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
1198 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
1199 /*     June 2017 */
1200
1201
1202 /*  =========================================================================== */
1203
1204
1205
1206
1207
1208 /*     Test the input arguments */
1209
1210     /* Parameter adjustments */
1211     --sva;
1212     a_dim1 = *lda;
1213     a_offset = 1 + a_dim1 * 1;
1214     a -= a_offset;
1215     u_dim1 = *ldu;
1216     u_offset = 1 + u_dim1 * 1;
1217     u -= u_offset;
1218     v_dim1 = *ldv;
1219     v_offset = 1 + v_dim1 * 1;
1220     v -= v_offset;
1221     --cwork;
1222     --rwork;
1223     --iwork;
1224
1225     /* Function Body */
1226     lsvec = lsame_(jobu, "U") || lsame_(jobu, "F");
1227     jracc = lsame_(jobv, "J");
1228     rsvec = lsame_(jobv, "V") || jracc;
1229     rowpiv = lsame_(joba, "F") || lsame_(joba, "G");
1230     l2rank = lsame_(joba, "R");
1231     l2aber = lsame_(joba, "A");
1232     errest = lsame_(joba, "E") || lsame_(joba, "G");
1233     l2tran = lsame_(jobt, "T") && *m == *n;
1234     l2kill = lsame_(jobr, "R");
1235     defr = lsame_(jobr, "N");
1236     l2pert = lsame_(jobp, "P");
1237
1238     lquery = *lwork == -1 || *lrwork == -1;
1239
1240     if (! (rowpiv || l2rank || l2aber || errest || lsame_(joba, "C"))) {
1241         *info = -1;
1242     } else if (! (lsvec || lsame_(jobu, "N") || lsame_(
1243             jobu, "W") && rsvec && l2tran)) {
1244         *info = -2;
1245     } else if (! (rsvec || lsame_(jobv, "N") || lsame_(
1246             jobv, "W") && lsvec && l2tran)) {
1247         *info = -3;
1248     } else if (! (l2kill || defr)) {
1249         *info = -4;
1250     } else if (! (lsame_(jobt, "T") || lsame_(jobt,
1251             "N"))) {
1252         *info = -5;
1253     } else if (! (l2pert || lsame_(jobp, "N"))) {
1254         *info = -6;
1255     } else if (*m < 0) {
1256         *info = -7;
1257     } else if (*n < 0 || *n > *m) {
1258         *info = -8;
1259     } else if (*lda < *m) {
1260         *info = -10;
1261     } else if (lsvec && *ldu < *m) {
1262         *info = -13;
1263     } else if (rsvec && *ldv < *n) {
1264         *info = -15;
1265     } else {
1266 /*        #:) */
1267         *info = 0;
1268     }
1269
1270     if (*info == 0) {
1271 /*         [[The expressions for computing the minimal and the optimal */
1272 /*         values of LCWORK, LRWORK are written with a lot of redundancy and */
1273 /*         can be simplified. However, this verbose form is useful for */
1274 /*         maintenance and modifications of the code.]] */
1275
1276 /*         CGEQRF of an N x N matrix, CGELQF of an N x N matrix, */
1277 /*         CUNMLQ for computing N x N matrix, CUNMQR for computing N x N */
1278 /*         matrix, CUNMQR for computing M x N matrix, respectively. */
1279         lwqp3 = *n + 1;
1280         lwqrf = f2cmax(1,*n);
1281         lwlqf = f2cmax(1,*n);
1282         lwunmlq = f2cmax(1,*n);
1283         lwunmqr = f2cmax(1,*n);
1284         lwunmqrm = f2cmax(1,*m);
1285         lwcon = *n << 1;
1286 /*         without and with explicit accumulation of Jacobi rotations */
1287 /* Computing MAX */
1288         i__1 = *n << 1;
1289         lwsvdj = f2cmax(i__1,1);
1290 /* Computing MAX */
1291         i__1 = *n << 1;
1292         lwsvdjv = f2cmax(i__1,1);
1293         lrwqp3 = *n << 1;
1294         lrwcon = *n;
1295         lrwsvdj = *n;
1296         if (lquery) {
1297             cgeqp3_(m, n, &a[a_offset], lda, &iwork[1], cdummy, cdummy, &c_n1,
1298                      rdummy, &ierr);
1299             lwrk_cgeqp3__ = cdummy[0].r;
1300             cgeqrf_(n, n, &a[a_offset], lda, cdummy, cdummy, &c_n1, &ierr);
1301             lwrk_cgeqrf__ = cdummy[0].r;
1302             cgelqf_(n, n, &a[a_offset], lda, cdummy, cdummy, &c_n1, &ierr);
1303             lwrk_cgelqf__ = cdummy[0].r;
1304         }
1305         minwrk = 2;
1306         optwrk = 2;
1307         miniwrk = *n;
1308         if (! (lsvec || rsvec)) {
1309 /*             only the singular values are requested */
1310             if (errest) {
1311 /* Computing MAX */
1312 /* Computing 2nd power */
1313                 i__3 = *n;
1314                 i__1 = *n + lwqp3, i__2 = i__3 * i__3 + lwcon, i__1 = f2cmax(
1315                         i__1,i__2), i__2 = *n + lwqrf, i__1 = f2cmax(i__1,i__2);
1316                 minwrk = f2cmax(i__1,lwsvdj);
1317             } else {
1318 /* Computing MAX */
1319                 i__1 = *n + lwqp3, i__2 = *n + lwqrf, i__1 = f2cmax(i__1,i__2);
1320                 minwrk = f2cmax(i__1,lwsvdj);
1321             }
1322             if (lquery) {
1323                 cgesvj_("L", "N", "N", n, n, &a[a_offset], lda, &sva[1], n, &
1324                         v[v_offset], ldv, cdummy, &c_n1, rdummy, &c_n1, &ierr);
1325                 lwrk_cgesvj__ = cdummy[0].r;
1326                 if (errest) {
1327 /* Computing MAX */
1328 /* Computing 2nd power */
1329                     i__3 = *n;
1330                     i__1 = *n + lwrk_cgeqp3__, i__2 = i__3 * i__3 + lwcon,
1331                             i__1 = f2cmax(i__1,i__2), i__2 = *n + lwrk_cgeqrf__,
1332                             i__1 = f2cmax(i__1,i__2);
1333                     optwrk = f2cmax(i__1,lwrk_cgesvj__);
1334                 } else {
1335 /* Computing MAX */
1336                     i__1 = *n + lwrk_cgeqp3__, i__2 = *n + lwrk_cgeqrf__,
1337                             i__1 = f2cmax(i__1,i__2);
1338                     optwrk = f2cmax(i__1,lwrk_cgesvj__);
1339                 }
1340             }
1341             if (l2tran || rowpiv) {
1342                 if (errest) {
1343 /* Computing MAX */
1344                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1345                             f2cmax(i__1,lrwqp3), i__1 = f2cmax(i__1,lrwcon);
1346                     minrwrk = f2cmax(i__1,lrwsvdj);
1347                 } else {
1348 /* Computing MAX */
1349                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1350                             f2cmax(i__1,lrwqp3);
1351                     minrwrk = f2cmax(i__1,lrwsvdj);
1352                 }
1353             } else {
1354                 if (errest) {
1355 /* Computing MAX */
1356                     i__1 = f2cmax(7,lrwqp3), i__1 = f2cmax(i__1,lrwcon);
1357                     minrwrk = f2cmax(i__1,lrwsvdj);
1358                 } else {
1359 /* Computing MAX */
1360                     i__1 = f2cmax(7,lrwqp3);
1361                     minrwrk = f2cmax(i__1,lrwsvdj);
1362                 }
1363             }
1364             if (rowpiv || l2tran) {
1365                 miniwrk += *m;
1366             }
1367         } else if (rsvec && ! lsvec) {
1368 /*            singular values and the right singular vectors are requested */
1369             if (errest) {
1370 /* Computing MAX */
1371                 i__1 = *n + lwqp3, i__1 = f2cmax(i__1,lwcon), i__1 = f2cmax(i__1,
1372                         lwsvdj), i__2 = *n + lwlqf, i__1 = f2cmax(i__1,i__2),
1373                         i__2 = (*n << 1) + lwqrf, i__1 = f2cmax(i__1,i__2), i__2
1374                         = *n + lwsvdj, i__1 = f2cmax(i__1,i__2), i__2 = *n +
1375                         lwunmlq;
1376                 minwrk = f2cmax(i__1,i__2);
1377             } else {
1378 /* Computing MAX */
1379                 i__1 = *n + lwqp3, i__1 = f2cmax(i__1,lwsvdj), i__2 = *n + lwlqf,
1380                          i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + lwqrf,
1381                         i__1 = f2cmax(i__1,i__2), i__2 = *n + lwsvdj, i__1 = f2cmax(
1382                         i__1,i__2), i__2 = *n + lwunmlq;
1383                 minwrk = f2cmax(i__1,i__2);
1384             }
1385             if (lquery) {
1386                 cgesvj_("L", "U", "N", n, n, &u[u_offset], ldu, &sva[1], n, &
1387                         a[a_offset], lda, cdummy, &c_n1, rdummy, &c_n1, &ierr);
1388                 lwrk_cgesvj__ = cdummy[0].r;
1389                 cunmlq_("L", "C", n, n, n, &a[a_offset], lda, cdummy, &v[
1390                         v_offset], ldv, cdummy, &c_n1, &ierr);
1391                 lwrk_cunmlq__ = cdummy[0].r;
1392                 if (errest) {
1393 /* Computing MAX */
1394                     i__1 = *n + lwrk_cgeqp3__, i__1 = f2cmax(i__1,lwcon), i__1 =
1395                             f2cmax(i__1,lwrk_cgesvj__), i__2 = *n +
1396                             lwrk_cgelqf__, i__1 = f2cmax(i__1,i__2), i__2 = (*n
1397                             << 1) + lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2),
1398                             i__2 = *n + lwrk_cgesvj__, i__1 = f2cmax(i__1,i__2),
1399                             i__2 = *n + lwrk_cunmlq__;
1400                     optwrk = f2cmax(i__1,i__2);
1401                 } else {
1402 /* Computing MAX */
1403                     i__1 = *n + lwrk_cgeqp3__, i__1 = f2cmax(i__1,lwrk_cgesvj__),
1404                              i__2 = *n + lwrk_cgelqf__, i__1 = f2cmax(i__1,i__2),
1405                              i__2 = (*n << 1) + lwrk_cgeqrf__, i__1 = f2cmax(
1406                             i__1,i__2), i__2 = *n + lwrk_cgesvj__, i__1 = f2cmax(
1407                             i__1,i__2), i__2 = *n + lwrk_cunmlq__;
1408                     optwrk = f2cmax(i__1,i__2);
1409                 }
1410             }
1411             if (l2tran || rowpiv) {
1412                 if (errest) {
1413 /* Computing MAX */
1414                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1415                             f2cmax(i__1,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1416                     minrwrk = f2cmax(i__1,lrwcon);
1417                 } else {
1418 /* Computing MAX */
1419                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1420                             f2cmax(i__1,lrwqp3);
1421                     minrwrk = f2cmax(i__1,lrwsvdj);
1422                 }
1423             } else {
1424                 if (errest) {
1425 /* Computing MAX */
1426                     i__1 = f2cmax(7,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1427                     minrwrk = f2cmax(i__1,lrwcon);
1428                 } else {
1429 /* Computing MAX */
1430                     i__1 = f2cmax(7,lrwqp3);
1431                     minrwrk = f2cmax(i__1,lrwsvdj);
1432                 }
1433             }
1434             if (rowpiv || l2tran) {
1435                 miniwrk += *m;
1436             }
1437         } else if (lsvec && ! rsvec) {
1438 /*            singular values and the left singular vectors are requested */
1439             if (errest) {
1440 /* Computing MAX */
1441                 i__1 = f2cmax(lwqp3,lwcon), i__2 = *n + lwqrf, i__1 = f2cmax(i__1,
1442                         i__2), i__1 = f2cmax(i__1,lwsvdj);
1443                 minwrk = *n + f2cmax(i__1,lwunmqrm);
1444             } else {
1445 /* Computing MAX */
1446                 i__1 = lwqp3, i__2 = *n + lwqrf, i__1 = f2cmax(i__1,i__2), i__1 =
1447                          f2cmax(i__1,lwsvdj);
1448                 minwrk = *n + f2cmax(i__1,lwunmqrm);
1449             }
1450             if (lquery) {
1451                 cgesvj_("L", "U", "N", n, n, &u[u_offset], ldu, &sva[1], n, &
1452                         a[a_offset], lda, cdummy, &c_n1, rdummy, &c_n1, &ierr);
1453                 lwrk_cgesvj__ = cdummy[0].r;
1454                 cunmqr_("L", "N", m, n, n, &a[a_offset], lda, cdummy, &u[
1455                         u_offset], ldu, cdummy, &c_n1, &ierr);
1456                 lwrk_cunmqrm__ = cdummy[0].r;
1457                 if (errest) {
1458 /* Computing MAX */
1459                     i__1 = f2cmax(lwrk_cgeqp3__,lwcon), i__2 = *n +
1460                             lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2), i__1 = f2cmax(
1461                             i__1,lwrk_cgesvj__);
1462                     optwrk = *n + f2cmax(i__1,lwrk_cunmqrm__);
1463                 } else {
1464 /* Computing MAX */
1465                     i__1 = lwrk_cgeqp3__, i__2 = *n + lwrk_cgeqrf__, i__1 =
1466                             f2cmax(i__1,i__2), i__1 = f2cmax(i__1,lwrk_cgesvj__);
1467                     optwrk = *n + f2cmax(i__1,lwrk_cunmqrm__);
1468                 }
1469             }
1470             if (l2tran || rowpiv) {
1471                 if (errest) {
1472 /* Computing MAX */
1473                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1474                             f2cmax(i__1,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1475                     minrwrk = f2cmax(i__1,lrwcon);
1476                 } else {
1477 /* Computing MAX */
1478                     i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 =
1479                             f2cmax(i__1,lrwqp3);
1480                     minrwrk = f2cmax(i__1,lrwsvdj);
1481                 }
1482             } else {
1483                 if (errest) {
1484 /* Computing MAX */
1485                     i__1 = f2cmax(7,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1486                     minrwrk = f2cmax(i__1,lrwcon);
1487                 } else {
1488 /* Computing MAX */
1489                     i__1 = f2cmax(7,lrwqp3);
1490                     minrwrk = f2cmax(i__1,lrwsvdj);
1491                 }
1492             }
1493             if (rowpiv || l2tran) {
1494                 miniwrk += *m;
1495             }
1496         } else {
1497 /*            full SVD is requested */
1498             if (! jracc) {
1499                 if (errest) {
1500 /* Computing MAX */
1501 /* Computing 2nd power */
1502                     i__3 = *n;
1503 /* Computing 2nd power */
1504                     i__4 = *n;
1505 /* Computing 2nd power */
1506                     i__5 = *n;
1507 /* Computing 2nd power */
1508                     i__6 = *n;
1509 /* Computing 2nd power */
1510                     i__7 = *n;
1511 /* Computing 2nd power */
1512                     i__8 = *n;
1513 /* Computing 2nd power */
1514                     i__9 = *n;
1515 /* Computing 2nd power */
1516                     i__10 = *n;
1517 /* Computing 2nd power */
1518                     i__11 = *n;
1519                     i__1 = *n + lwqp3, i__2 = *n + lwcon, i__1 = f2cmax(i__1,
1520                             i__2), i__2 = (*n << 1) + i__3 * i__3 + lwcon,
1521                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + lwqrf,
1522                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + lwqp3,
1523                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + i__4 *
1524                             i__4 + *n + lwlqf, i__1 = f2cmax(i__1,i__2), i__2 = (
1525                             *n << 1) + i__5 * i__5 + *n + i__6 * i__6 + lwcon,
1526                              i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + i__7 *
1527                             i__7 + *n + lwsvdj, i__1 = f2cmax(i__1,i__2), i__2 =
1528                             (*n << 1) + i__8 * i__8 + *n + lwsvdjv, i__1 =
1529                             f2cmax(i__1,i__2), i__2 = (*n << 1) + i__9 * i__9 + *
1530                             n + lwunmqr, i__1 = f2cmax(i__1,i__2), i__2 = (*n <<
1531                             1) + i__10 * i__10 + *n + lwunmlq, i__1 = f2cmax(
1532                             i__1,i__2), i__2 = *n + i__11 * i__11 + lwsvdj,
1533                             i__1 = f2cmax(i__1,i__2), i__2 = *n + lwunmqrm;
1534                     minwrk = f2cmax(i__1,i__2);
1535                 } else {
1536 /* Computing MAX */
1537 /* Computing 2nd power */
1538                     i__3 = *n;
1539 /* Computing 2nd power */
1540                     i__4 = *n;
1541 /* Computing 2nd power */
1542                     i__5 = *n;
1543 /* Computing 2nd power */
1544                     i__6 = *n;
1545 /* Computing 2nd power */
1546                     i__7 = *n;
1547 /* Computing 2nd power */
1548                     i__8 = *n;
1549 /* Computing 2nd power */
1550                     i__9 = *n;
1551 /* Computing 2nd power */
1552                     i__10 = *n;
1553 /* Computing 2nd power */
1554                     i__11 = *n;
1555                     i__1 = *n + lwqp3, i__2 = (*n << 1) + i__3 * i__3 + lwcon,
1556                              i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + lwqrf,
1557                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + lwqp3,
1558                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + i__4 *
1559                             i__4 + *n + lwlqf, i__1 = f2cmax(i__1,i__2), i__2 = (
1560                             *n << 1) + i__5 * i__5 + *n + i__6 * i__6 + lwcon,
1561                              i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + i__7 *
1562                             i__7 + *n + lwsvdj, i__1 = f2cmax(i__1,i__2), i__2 =
1563                             (*n << 1) + i__8 * i__8 + *n + lwsvdjv, i__1 =
1564                             f2cmax(i__1,i__2), i__2 = (*n << 1) + i__9 * i__9 + *
1565                             n + lwunmqr, i__1 = f2cmax(i__1,i__2), i__2 = (*n <<
1566                             1) + i__10 * i__10 + *n + lwunmlq, i__1 = f2cmax(
1567                             i__1,i__2), i__2 = *n + i__11 * i__11 + lwsvdj,
1568                             i__1 = f2cmax(i__1,i__2), i__2 = *n + lwunmqrm;
1569                     minwrk = f2cmax(i__1,i__2);
1570                 }
1571                 miniwrk += *n;
1572                 if (rowpiv || l2tran) {
1573                     miniwrk += *m;
1574                 }
1575             } else {
1576                 if (errest) {
1577 /* Computing MAX */
1578 /* Computing 2nd power */
1579                     i__3 = *n;
1580 /* Computing 2nd power */
1581                     i__4 = *n;
1582                     i__1 = *n + lwqp3, i__2 = *n + lwcon, i__1 = f2cmax(i__1,
1583                             i__2), i__2 = (*n << 1) + lwqrf, i__1 = f2cmax(i__1,
1584                             i__2), i__2 = (*n << 1) + i__3 * i__3 + lwsvdjv,
1585                             i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) + i__4 *
1586                             i__4 + *n + lwunmqr, i__1 = f2cmax(i__1,i__2), i__2 =
1587                              *n + lwunmqrm;
1588                     minwrk = f2cmax(i__1,i__2);
1589                 } else {
1590 /* Computing MAX */
1591 /* Computing 2nd power */
1592                     i__3 = *n;
1593 /* Computing 2nd power */
1594                     i__4 = *n;
1595                     i__1 = *n + lwqp3, i__2 = (*n << 1) + lwqrf, i__1 = f2cmax(
1596                             i__1,i__2), i__2 = (*n << 1) + i__3 * i__3 +
1597                             lwsvdjv, i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1)
1598                             + i__4 * i__4 + *n + lwunmqr, i__1 = f2cmax(i__1,
1599                             i__2), i__2 = *n + lwunmqrm;
1600                     minwrk = f2cmax(i__1,i__2);
1601                 }
1602                 if (rowpiv || l2tran) {
1603                     miniwrk += *m;
1604                 }
1605             }
1606             if (lquery) {
1607                 cunmqr_("L", "N", m, n, n, &a[a_offset], lda, cdummy, &u[
1608                         u_offset], ldu, cdummy, &c_n1, &ierr);
1609                 lwrk_cunmqrm__ = cdummy[0].r;
1610                 cunmqr_("L", "N", n, n, n, &a[a_offset], lda, cdummy, &u[
1611                         u_offset], ldu, cdummy, &c_n1, &ierr);
1612                 lwrk_cunmqr__ = cdummy[0].r;
1613                 if (! jracc) {
1614                     cgeqp3_(n, n, &a[a_offset], lda, &iwork[1], cdummy,
1615                             cdummy, &c_n1, rdummy, &ierr);
1616                     lwrk_cgeqp3n__ = cdummy[0].r;
1617                     cgesvj_("L", "U", "N", n, n, &u[u_offset], ldu, &sva[1],
1618                             n, &v[v_offset], ldv, cdummy, &c_n1, rdummy, &
1619                             c_n1, &ierr);
1620                     lwrk_cgesvj__ = cdummy[0].r;
1621                     cgesvj_("U", "U", "N", n, n, &u[u_offset], ldu, &sva[1],
1622                             n, &v[v_offset], ldv, cdummy, &c_n1, rdummy, &
1623                             c_n1, &ierr);
1624                     lwrk_cgesvju__ = cdummy[0].r;
1625                     cgesvj_("L", "U", "V", n, n, &u[u_offset], ldu, &sva[1],
1626                             n, &v[v_offset], ldv, cdummy, &c_n1, rdummy, &
1627                             c_n1, &ierr);
1628                     lwrk_cgesvjv__ = cdummy[0].r;
1629                     cunmlq_("L", "C", n, n, n, &a[a_offset], lda, cdummy, &v[
1630                             v_offset], ldv, cdummy, &c_n1, &ierr);
1631                     lwrk_cunmlq__ = cdummy[0].r;
1632                     if (errest) {
1633 /* Computing MAX */
1634 /* Computing 2nd power */
1635                         i__3 = *n;
1636 /* Computing 2nd power */
1637                         i__4 = *n;
1638 /* Computing 2nd power */
1639                         i__5 = *n;
1640 /* Computing 2nd power */
1641                         i__6 = *n;
1642 /* Computing 2nd power */
1643                         i__7 = *n;
1644 /* Computing 2nd power */
1645                         i__8 = *n;
1646 /* Computing 2nd power */
1647                         i__9 = *n;
1648 /* Computing 2nd power */
1649                         i__10 = *n;
1650 /* Computing 2nd power */
1651                         i__11 = *n;
1652                         i__1 = *n + lwrk_cgeqp3__, i__2 = *n + lwcon, i__1 =
1653                                 f2cmax(i__1,i__2), i__2 = (*n << 1) + i__3 *
1654                                 i__3 + lwcon, i__1 = f2cmax(i__1,i__2), i__2 = (*
1655                                 n << 1) + lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2)
1656                                 , i__2 = (*n << 1) + lwrk_cgeqp3n__, i__1 =
1657                                 f2cmax(i__1,i__2), i__2 = (*n << 1) + i__4 *
1658                                 i__4 + *n + lwrk_cgelqf__, i__1 = f2cmax(i__1,
1659                                 i__2), i__2 = (*n << 1) + i__5 * i__5 + *n +
1660                                 i__6 * i__6 + lwcon, i__1 = f2cmax(i__1,i__2),
1661                                 i__2 = (*n << 1) + i__7 * i__7 + *n +
1662                                 lwrk_cgesvj__, i__1 = f2cmax(i__1,i__2), i__2 = (
1663                                 *n << 1) + i__8 * i__8 + *n + lwrk_cgesvjv__,
1664                                 i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) +
1665                                 i__9 * i__9 + *n + lwrk_cunmqr__, i__1 = f2cmax(
1666                                 i__1,i__2), i__2 = (*n << 1) + i__10 * i__10
1667                                 + *n + lwrk_cunmlq__, i__1 = f2cmax(i__1,i__2),
1668                                 i__2 = *n + i__11 * i__11 + lwrk_cgesvju__,
1669                                 i__1 = f2cmax(i__1,i__2), i__2 = *n +
1670                                 lwrk_cunmqrm__;
1671                         optwrk = f2cmax(i__1,i__2);
1672                     } else {
1673 /* Computing MAX */
1674 /* Computing 2nd power */
1675                         i__3 = *n;
1676 /* Computing 2nd power */
1677                         i__4 = *n;
1678 /* Computing 2nd power */
1679                         i__5 = *n;
1680 /* Computing 2nd power */
1681                         i__6 = *n;
1682 /* Computing 2nd power */
1683                         i__7 = *n;
1684 /* Computing 2nd power */
1685                         i__8 = *n;
1686 /* Computing 2nd power */
1687                         i__9 = *n;
1688 /* Computing 2nd power */
1689                         i__10 = *n;
1690 /* Computing 2nd power */
1691                         i__11 = *n;
1692                         i__1 = *n + lwrk_cgeqp3__, i__2 = (*n << 1) + i__3 *
1693                                 i__3 + lwcon, i__1 = f2cmax(i__1,i__2), i__2 = (*
1694                                 n << 1) + lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2)
1695                                 , i__2 = (*n << 1) + lwrk_cgeqp3n__, i__1 =
1696                                 f2cmax(i__1,i__2), i__2 = (*n << 1) + i__4 *
1697                                 i__4 + *n + lwrk_cgelqf__, i__1 = f2cmax(i__1,
1698                                 i__2), i__2 = (*n << 1) + i__5 * i__5 + *n +
1699                                 i__6 * i__6 + lwcon, i__1 = f2cmax(i__1,i__2),
1700                                 i__2 = (*n << 1) + i__7 * i__7 + *n +
1701                                 lwrk_cgesvj__, i__1 = f2cmax(i__1,i__2), i__2 = (
1702                                 *n << 1) + i__8 * i__8 + *n + lwrk_cgesvjv__,
1703                                 i__1 = f2cmax(i__1,i__2), i__2 = (*n << 1) +
1704                                 i__9 * i__9 + *n + lwrk_cunmqr__, i__1 = f2cmax(
1705                                 i__1,i__2), i__2 = (*n << 1) + i__10 * i__10
1706                                 + *n + lwrk_cunmlq__, i__1 = f2cmax(i__1,i__2),
1707                                 i__2 = *n + i__11 * i__11 + lwrk_cgesvju__,
1708                                 i__1 = f2cmax(i__1,i__2), i__2 = *n +
1709                                 lwrk_cunmqrm__;
1710                         optwrk = f2cmax(i__1,i__2);
1711                     }
1712                 } else {
1713                     cgesvj_("L", "U", "V", n, n, &u[u_offset], ldu, &sva[1],
1714                             n, &v[v_offset], ldv, cdummy, &c_n1, rdummy, &
1715                             c_n1, &ierr);
1716                     lwrk_cgesvjv__ = cdummy[0].r;
1717                     cunmqr_("L", "N", n, n, n, cdummy, n, cdummy, &v[v_offset]
1718                             , ldv, cdummy, &c_n1, &ierr)
1719                             ;
1720                     lwrk_cunmqr__ = cdummy[0].r;
1721                     cunmqr_("L", "N", m, n, n, &a[a_offset], lda, cdummy, &u[
1722                             u_offset], ldu, cdummy, &c_n1, &ierr);
1723                     lwrk_cunmqrm__ = cdummy[0].r;
1724                     if (errest) {
1725 /* Computing MAX */
1726 /* Computing 2nd power */
1727                         i__3 = *n;
1728 /* Computing 2nd power */
1729                         i__4 = *n;
1730 /* Computing 2nd power */
1731                         i__5 = *n;
1732                         i__1 = *n + lwrk_cgeqp3__, i__2 = *n + lwcon, i__1 =
1733                                 f2cmax(i__1,i__2), i__2 = (*n << 1) +
1734                                 lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2), i__2 = (
1735                                 *n << 1) + i__3 * i__3, i__1 = f2cmax(i__1,i__2),
1736                                  i__2 = (*n << 1) + i__4 * i__4 +
1737                                 lwrk_cgesvjv__, i__1 = f2cmax(i__1,i__2), i__2 =
1738                                 (*n << 1) + i__5 * i__5 + *n + lwrk_cunmqr__,
1739                                 i__1 = f2cmax(i__1,i__2), i__2 = *n +
1740                                 lwrk_cunmqrm__;
1741                         optwrk = f2cmax(i__1,i__2);
1742                     } else {
1743 /* Computing MAX */
1744 /* Computing 2nd power */
1745                         i__3 = *n;
1746 /* Computing 2nd power */
1747                         i__4 = *n;
1748 /* Computing 2nd power */
1749                         i__5 = *n;
1750                         i__1 = *n + lwrk_cgeqp3__, i__2 = (*n << 1) +
1751                                 lwrk_cgeqrf__, i__1 = f2cmax(i__1,i__2), i__2 = (
1752                                 *n << 1) + i__3 * i__3, i__1 = f2cmax(i__1,i__2),
1753                                  i__2 = (*n << 1) + i__4 * i__4 +
1754                                 lwrk_cgesvjv__, i__1 = f2cmax(i__1,i__2), i__2 =
1755                                 (*n << 1) + i__5 * i__5 + *n + lwrk_cunmqr__,
1756                                 i__1 = f2cmax(i__1,i__2), i__2 = *n +
1757                                 lwrk_cunmqrm__;
1758                         optwrk = f2cmax(i__1,i__2);
1759                     }
1760                 }
1761             }
1762             if (l2tran || rowpiv) {
1763 /* Computing MAX */
1764                 i__1 = 7, i__2 = *m << 1, i__1 = f2cmax(i__1,i__2), i__1 = f2cmax(
1765                         i__1,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1766                 minrwrk = f2cmax(i__1,lrwcon);
1767             } else {
1768 /* Computing MAX */
1769                 i__1 = f2cmax(7,lrwqp3), i__1 = f2cmax(i__1,lrwsvdj);
1770                 minrwrk = f2cmax(i__1,lrwcon);
1771             }
1772         }
1773         minwrk = f2cmax(2,minwrk);
1774         optwrk = f2cmax(optwrk,minwrk);
1775         if (*lwork < minwrk && ! lquery) {
1776             *info = -17;
1777         }
1778         if (*lrwork < minrwrk && ! lquery) {
1779             *info = -19;
1780         }
1781     }
1782
1783     if (*info != 0) {
1784 /*       #:( */
1785         i__1 = -(*info);
1786         xerbla_("CGEJSV", &i__1, (ftnlen)6);
1787         return 0;
1788     } else if (lquery) {
1789         cwork[1].r = (real) optwrk, cwork[1].i = 0.f;
1790         cwork[2].r = (real) minwrk, cwork[2].i = 0.f;
1791         rwork[1] = (real) minrwrk;
1792         iwork[1] = f2cmax(4,miniwrk);
1793         return 0;
1794     }
1795
1796 /*     Quick return for void matrix (Y3K safe) */
1797 /* #:) */
1798     if (*m == 0 || *n == 0) {
1799         iwork[1] = 0;
1800         iwork[2] = 0;
1801         iwork[3] = 0;
1802         iwork[4] = 0;
1803         rwork[1] = 0.f;
1804         rwork[2] = 0.f;
1805         rwork[3] = 0.f;
1806         rwork[4] = 0.f;
1807         rwork[5] = 0.f;
1808         rwork[6] = 0.f;
1809         rwork[7] = 0.f;
1810         return 0;
1811     }
1812
1813 /*     Determine whether the matrix U should be M x N or M x M */
1814
1815     if (lsvec) {
1816         n1 = *n;
1817         if (lsame_(jobu, "F")) {
1818             n1 = *m;
1819         }
1820     }
1821
1822 /*     Set numerical parameters */
1823
1824 /* !    NOTE: Make sure SLAMCH() does not fail on the target architecture. */
1825
1826     epsln = slamch_("Epsilon");
1827     sfmin = slamch_("SafeMinimum");
1828     small = sfmin / epsln;
1829     big = slamch_("O");
1830 /*     BIG   = ONE / SFMIN */
1831
1832 /*     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N */
1833
1834 /* (!)  If necessary, scale SVA() to protect the largest norm from */
1835 /*     overflow. It is possible that this scaling pushes the smallest */
1836 /*     column norm left from the underflow threshold (extreme case). */
1837
1838     scalem = 1.f / sqrt((real) (*m) * (real) (*n));
1839     noscal = TRUE_;
1840     goscal = TRUE_;
1841     i__1 = *n;
1842     for (p = 1; p <= i__1; ++p) {
1843         aapp = 0.f;
1844         aaqq = 1.f;
1845         classq_(m, &a[p * a_dim1 + 1], &c__1, &aapp, &aaqq);
1846         if (aapp > big) {
1847             *info = -9;
1848             i__2 = -(*info);
1849             xerbla_("CGEJSV", &i__2, (ftnlen)6);
1850             return 0;
1851         }
1852         aaqq = sqrt(aaqq);
1853         if (aapp < big / aaqq && noscal) {
1854             sva[p] = aapp * aaqq;
1855         } else {
1856             noscal = FALSE_;
1857             sva[p] = aapp * (aaqq * scalem);
1858             if (goscal) {
1859                 goscal = FALSE_;
1860                 i__2 = p - 1;
1861                 sscal_(&i__2, &scalem, &sva[1], &c__1);
1862             }
1863         }
1864 /* L1874: */
1865     }
1866
1867     if (noscal) {
1868         scalem = 1.f;
1869     }
1870
1871     aapp = 0.f;
1872     aaqq = big;
1873     i__1 = *n;
1874     for (p = 1; p <= i__1; ++p) {
1875 /* Computing MAX */
1876         r__1 = aapp, r__2 = sva[p];
1877         aapp = f2cmax(r__1,r__2);
1878         if (sva[p] != 0.f) {
1879 /* Computing MIN */
1880             r__1 = aaqq, r__2 = sva[p];
1881             aaqq = f2cmin(r__1,r__2);
1882         }
1883 /* L4781: */
1884     }
1885
1886 /*     Quick return for zero M x N matrix */
1887 /* #:) */
1888     if (aapp == 0.f) {
1889         if (lsvec) {
1890             claset_("G", m, &n1, &c_b1, &c_b2, &u[u_offset], ldu);
1891         }
1892         if (rsvec) {
1893             claset_("G", n, n, &c_b1, &c_b2, &v[v_offset], ldv);
1894         }
1895         rwork[1] = 1.f;
1896         rwork[2] = 1.f;
1897         if (errest) {
1898             rwork[3] = 1.f;
1899         }
1900         if (lsvec && rsvec) {
1901             rwork[4] = 1.f;
1902             rwork[5] = 1.f;
1903         }
1904         if (l2tran) {
1905             rwork[6] = 0.f;
1906             rwork[7] = 0.f;
1907         }
1908         iwork[1] = 0;
1909         iwork[2] = 0;
1910         iwork[3] = 0;
1911         iwork[4] = -1;
1912         return 0;
1913     }
1914
1915 /*     Issue warning if denormalized column norms detected. Override the */
1916 /*     high relative accuracy request. Issue licence to kill nonzero columns */
1917 /*     (set them to zero) whose norm is less than sigma_max / BIG (roughly). */
1918 /* #:( */
1919     warning = 0;
1920     if (aaqq <= sfmin) {
1921         l2rank = TRUE_;
1922         l2kill = TRUE_;
1923         warning = 1;
1924     }
1925
1926 /*     Quick return for one-column matrix */
1927 /* #:) */
1928     if (*n == 1) {
1929
1930         if (lsvec) {
1931             clascl_("G", &c__0, &c__0, &sva[1], &scalem, m, &c__1, &a[a_dim1
1932                     + 1], lda, &ierr);
1933             clacpy_("A", m, &c__1, &a[a_offset], lda, &u[u_offset], ldu);
1934 /*           computing all M left singular vectors of the M x 1 matrix */
1935             if (n1 != *n) {
1936                 i__1 = *lwork - *n;
1937                 cgeqrf_(m, n, &u[u_offset], ldu, &cwork[1], &cwork[*n + 1], &
1938                         i__1, &ierr);
1939                 i__1 = *lwork - *n;
1940                 cungqr_(m, &n1, &c__1, &u[u_offset], ldu, &cwork[1], &cwork[*
1941                         n + 1], &i__1, &ierr);
1942                 ccopy_(m, &a[a_dim1 + 1], &c__1, &u[u_dim1 + 1], &c__1);
1943             }
1944         }
1945         if (rsvec) {
1946             i__1 = v_dim1 + 1;
1947             v[i__1].r = 1.f, v[i__1].i = 0.f;
1948         }
1949         if (sva[1] < big * scalem) {
1950             sva[1] /= scalem;
1951             scalem = 1.f;
1952         }
1953         rwork[1] = 1.f / scalem;
1954         rwork[2] = 1.f;
1955         if (sva[1] != 0.f) {
1956             iwork[1] = 1;
1957             if (sva[1] / scalem >= sfmin) {
1958                 iwork[2] = 1;
1959             } else {
1960                 iwork[2] = 0;
1961             }
1962         } else {
1963             iwork[1] = 0;
1964             iwork[2] = 0;
1965         }
1966         iwork[3] = 0;
1967         iwork[4] = -1;
1968         if (errest) {
1969             rwork[3] = 1.f;
1970         }
1971         if (lsvec && rsvec) {
1972             rwork[4] = 1.f;
1973             rwork[5] = 1.f;
1974         }
1975         if (l2tran) {
1976             rwork[6] = 0.f;
1977             rwork[7] = 0.f;
1978         }
1979         return 0;
1980
1981     }
1982
1983     transp = FALSE_;
1984
1985     aatmax = -1.f;
1986     aatmin = big;
1987     if (rowpiv || l2tran) {
1988
1989 /*     Compute the row norms, needed to determine row pivoting sequence */
1990 /*     (in the case of heavily row weighted A, row pivoting is strongly */
1991 /*     advised) and to collect information needed to compare the */
1992 /*     structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.). */
1993
1994         if (l2tran) {
1995             i__1 = *m;
1996             for (p = 1; p <= i__1; ++p) {
1997                 xsc = 0.f;
1998                 temp1 = 1.f;
1999                 classq_(n, &a[p + a_dim1], lda, &xsc, &temp1);
2000 /*              CLASSQ gets both the ell_2 and the ell_infinity norm */
2001 /*              in one pass through the vector */
2002                 rwork[*m + p] = xsc * scalem;
2003                 rwork[p] = xsc * (scalem * sqrt(temp1));
2004 /* Computing MAX */
2005                 r__1 = aatmax, r__2 = rwork[p];
2006                 aatmax = f2cmax(r__1,r__2);
2007                 if (rwork[p] != 0.f) {
2008 /* Computing MIN */
2009                     r__1 = aatmin, r__2 = rwork[p];
2010                     aatmin = f2cmin(r__1,r__2);
2011                 }
2012 /* L1950: */
2013             }
2014         } else {
2015             i__1 = *m;
2016             for (p = 1; p <= i__1; ++p) {
2017                 rwork[*m + p] = scalem * c_abs(&a[p + icamax_(n, &a[p +
2018                         a_dim1], lda) * a_dim1]);
2019 /* Computing MAX */
2020                 r__1 = aatmax, r__2 = rwork[*m + p];
2021                 aatmax = f2cmax(r__1,r__2);
2022 /* Computing MIN */
2023                 r__1 = aatmin, r__2 = rwork[*m + p];
2024                 aatmin = f2cmin(r__1,r__2);
2025 /* L1904: */
2026             }
2027         }
2028
2029     }
2030
2031 /*     For square matrix A try to determine whether A^*  would be better */
2032 /*     input for the preconditioned Jacobi SVD, with faster convergence. */
2033 /*     The decision is based on an O(N) function of the vector of column */
2034 /*     and row norms of A, based on the Shannon entropy. This should give */
2035 /*     the right choice in most cases when the difference actually matters. */
2036 /*     It may fail and pick the slower converging side. */
2037
2038     entra = 0.f;
2039     entrat = 0.f;
2040     if (l2tran) {
2041
2042         xsc = 0.f;
2043         temp1 = 1.f;
2044         slassq_(n, &sva[1], &c__1, &xsc, &temp1);
2045         temp1 = 1.f / temp1;
2046
2047         entra = 0.f;
2048         i__1 = *n;
2049         for (p = 1; p <= i__1; ++p) {
2050 /* Computing 2nd power */
2051             r__1 = sva[p] / xsc;
2052             big1 = r__1 * r__1 * temp1;
2053             if (big1 != 0.f) {
2054                 entra += big1 * log(big1);
2055             }
2056 /* L1113: */
2057         }
2058         entra = -entra / log((real) (*n));
2059
2060 /*        Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex. */
2061 /*        It is derived from the diagonal of  A^* * A.  Do the same with the */
2062 /*        diagonal of A * A^*, compute the entropy of the corresponding */
2063 /*        probability distribution. Note that A * A^* and A^* * A have the */
2064 /*        same trace. */
2065
2066         entrat = 0.f;
2067         i__1 = *m;
2068         for (p = 1; p <= i__1; ++p) {
2069 /* Computing 2nd power */
2070             r__1 = rwork[p] / xsc;
2071             big1 = r__1 * r__1 * temp1;
2072             if (big1 != 0.f) {
2073                 entrat += big1 * log(big1);
2074             }
2075 /* L1114: */
2076         }
2077         entrat = -entrat / log((real) (*m));
2078
2079 /*        Analyze the entropies and decide A or A^*. Smaller entropy */
2080 /*        usually means better input for the algorithm. */
2081
2082         transp = entrat < entra;
2083
2084 /*        If A^* is better than A, take the adjoint of A. This is allowed */
2085 /*        only for square matrices, M=N. */
2086         if (transp) {
2087 /*           In an optimal implementation, this trivial transpose */
2088 /*           should be replaced with faster transpose. */
2089             i__1 = *n - 1;
2090             for (p = 1; p <= i__1; ++p) {
2091                 i__2 = p + p * a_dim1;
2092                 r_cnjg(&q__1, &a[p + p * a_dim1]);
2093                 a[i__2].r = q__1.r, a[i__2].i = q__1.i;
2094                 i__2 = *n;
2095                 for (q = p + 1; q <= i__2; ++q) {
2096                     r_cnjg(&q__1, &a[q + p * a_dim1]);
2097                     ctemp.r = q__1.r, ctemp.i = q__1.i;
2098                     i__3 = q + p * a_dim1;
2099                     r_cnjg(&q__1, &a[p + q * a_dim1]);
2100                     a[i__3].r = q__1.r, a[i__3].i = q__1.i;
2101                     i__3 = p + q * a_dim1;
2102                     a[i__3].r = ctemp.r, a[i__3].i = ctemp.i;
2103 /* L1116: */
2104                 }
2105 /* L1115: */
2106             }
2107             i__1 = *n + *n * a_dim1;
2108             r_cnjg(&q__1, &a[*n + *n * a_dim1]);
2109             a[i__1].r = q__1.r, a[i__1].i = q__1.i;
2110             i__1 = *n;
2111             for (p = 1; p <= i__1; ++p) {
2112                 rwork[*m + p] = sva[p];
2113                 sva[p] = rwork[p];
2114 /*              previously computed row 2-norms are now column 2-norms */
2115 /*              of the transposed matrix */
2116 /* L1117: */
2117             }
2118             temp1 = aapp;
2119             aapp = aatmax;
2120             aatmax = temp1;
2121             temp1 = aaqq;
2122             aaqq = aatmin;
2123             aatmin = temp1;
2124             kill = lsvec;
2125             lsvec = rsvec;
2126             rsvec = kill;
2127             if (lsvec) {
2128                 n1 = *n;
2129             }
2130
2131             rowpiv = TRUE_;
2132         }
2133
2134     }
2135 /*     END IF L2TRAN */
2136
2137 /*     Scale the matrix so that its maximal singular value remains less */
2138 /*     than SQRT(BIG) -- the matrix is scaled so that its maximal column */
2139 /*     has Euclidean norm equal to SQRT(BIG/N). The only reason to keep */
2140 /*     SQRT(BIG) instead of BIG is the fact that CGEJSV uses LAPACK and */
2141 /*     BLAS routines that, in some implementations, are not capable of */
2142 /*     working in the full interval [SFMIN,BIG] and that they may provoke */
2143 /*     overflows in the intermediate results. If the singular values spread */
2144 /*     from SFMIN to BIG, then CGESVJ will compute them. So, in that case, */
2145 /*     one should use CGESVJ instead of CGEJSV. */
2146     big1 = sqrt(big);
2147     temp1 = sqrt(big / (real) (*n));
2148 /*     >> for future updates: allow bigger range, i.e. the largest column */
2149 /*     will be allowed up to BIG/N and CGESVJ will do the rest. However, for */
2150 /*     this all other (LAPACK) components must allow such a range. */
2151 /*     TEMP1  = BIG/REAL(N) */
2152 /*     TEMP1  = BIG * EPSLN  this should 'almost' work with current LAPACK components */
2153     slascl_("G", &c__0, &c__0, &aapp, &temp1, n, &c__1, &sva[1], n, &ierr);
2154     if (aaqq > aapp * sfmin) {
2155         aaqq = aaqq / aapp * temp1;
2156     } else {
2157         aaqq = aaqq * temp1 / aapp;
2158     }
2159     temp1 *= scalem;
2160     clascl_("G", &c__0, &c__0, &aapp, &temp1, m, n, &a[a_offset], lda, &ierr);
2161
2162 /*     To undo scaling at the end of this procedure, multiply the */
2163 /*     computed singular values with USCAL2 / USCAL1. */
2164
2165     uscal1 = temp1;
2166     uscal2 = aapp;
2167
2168     if (l2kill) {
2169 /*        L2KILL enforces computation of nonzero singular values in */
2170 /*        the restricted range of condition number of the initial A, */
2171 /*        sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN). */
2172         xsc = sqrt(sfmin);
2173     } else {
2174         xsc = small;
2175
2176 /*        Now, if the condition number of A is too big, */
2177 /*        sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN, */
2178 /*        as a precaution measure, the full SVD is computed using CGESVJ */
2179 /*        with accumulated Jacobi rotations. This provides numerically */
2180 /*        more robust computation, at the cost of slightly increased run */
2181 /*        time. Depending on the concrete implementation of BLAS and LAPACK */
2182 /*        (i.e. how they behave in presence of extreme ill-conditioning) the */
2183 /*        implementor may decide to remove this switch. */
2184         if (aaqq < sqrt(sfmin) && lsvec && rsvec) {
2185             jracc = TRUE_;
2186         }
2187
2188     }
2189     if (aaqq < xsc) {
2190         i__1 = *n;
2191         for (p = 1; p <= i__1; ++p) {
2192             if (sva[p] < xsc) {
2193                 claset_("A", m, &c__1, &c_b1, &c_b1, &a[p * a_dim1 + 1], lda);
2194                 sva[p] = 0.f;
2195             }
2196 /* L700: */
2197         }
2198     }
2199
2200 /*     Preconditioning using QR factorization with pivoting */
2201
2202     if (rowpiv) {
2203 /*        Optional row permutation (Bjoerck row pivoting): */
2204 /*        A result by Cox and Higham shows that the Bjoerck's */
2205 /*        row pivoting combined with standard column pivoting */
2206 /*        has similar effect as Powell-Reid complete pivoting. */
2207 /*        The ell-infinity norms of A are made nonincreasing. */
2208         if (lsvec && rsvec && ! jracc) {
2209             iwoff = *n << 1;
2210         } else {
2211             iwoff = *n;
2212         }
2213         i__1 = *m - 1;
2214         for (p = 1; p <= i__1; ++p) {
2215             i__2 = *m - p + 1;
2216             q = isamax_(&i__2, &rwork[*m + p], &c__1) + p - 1;
2217             iwork[iwoff + p] = q;
2218             if (p != q) {
2219                 temp1 = rwork[*m + p];
2220                 rwork[*m + p] = rwork[*m + q];
2221                 rwork[*m + q] = temp1;
2222             }
2223 /* L1952: */
2224         }
2225         i__1 = *m - 1;
2226         claswp_(n, &a[a_offset], lda, &c__1, &i__1, &iwork[iwoff + 1], &c__1);
2227     }
2228
2229 /*     End of the preparation phase (scaling, optional sorting and */
2230 /*     transposing, optional flushing of small columns). */
2231
2232 /*     Preconditioning */
2233
2234 /*     If the full SVD is needed, the right singular vectors are computed */
2235 /*     from a matrix equation, and for that we need theoretical analysis */
2236 /*     of the Businger-Golub pivoting. So we use CGEQP3 as the first RR QRF. */
2237 /*     In all other cases the first RR QRF can be chosen by other criteria */
2238 /*     (eg speed by replacing global with restricted window pivoting, such */
2239 /*     as in xGEQPX from TOMS # 782). Good results will be obtained using */
2240 /*     xGEQPX with properly (!) chosen numerical parameters. */
2241 /*     Any improvement of CGEQP3 improves overal performance of CGEJSV. */
2242
2243 /*     A * P1 = Q1 * [ R1^* 0]^*: */
2244     i__1 = *n;
2245     for (p = 1; p <= i__1; ++p) {
2246         iwork[p] = 0;
2247 /* L1963: */
2248     }
2249     i__1 = *lwork - *n;
2250     cgeqp3_(m, n, &a[a_offset], lda, &iwork[1], &cwork[1], &cwork[*n + 1], &
2251             i__1, &rwork[1], &ierr);
2252
2253 /*     The upper triangular matrix R1 from the first QRF is inspected for */
2254 /*     rank deficiency and possibilities for deflation, or possible */
2255 /*     ill-conditioning. Depending on the user specified flag L2RANK, */
2256 /*     the procedure explores possibilities to reduce the numerical */
2257 /*     rank by inspecting the computed upper triangular factor. If */
2258 /*     L2RANK or L2ABER are up, then CGEJSV will compute the SVD of */
2259 /*     A + dA, where ||dA|| <= f(M,N)*EPSLN. */
2260
2261     nr = 1;
2262     if (l2aber) {
2263 /*        Standard absolute error bound suffices. All sigma_i with */
2264 /*        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an */
2265 /*        aggressive enforcement of lower numerical rank by introducing a */
2266 /*        backward error of the order of N*EPSLN*||A||. */
2267         temp1 = sqrt((real) (*n)) * epsln;
2268         i__1 = *n;
2269         for (p = 2; p <= i__1; ++p) {
2270             if (c_abs(&a[p + p * a_dim1]) >= temp1 * c_abs(&a[a_dim1 + 1])) {
2271                 ++nr;
2272             } else {
2273                 goto L3002;
2274             }
2275 /* L3001: */
2276         }
2277 L3002:
2278         ;
2279     } else if (l2rank) {
2280 /*        Sudden drop on the diagonal of R1 is used as the criterion for */
2281 /*        close-to-rank-deficient. */
2282         temp1 = sqrt(sfmin);
2283         i__1 = *n;
2284         for (p = 2; p <= i__1; ++p) {
2285             if (c_abs(&a[p + p * a_dim1]) < epsln * c_abs(&a[p - 1 + (p - 1) *
2286                      a_dim1]) || c_abs(&a[p + p * a_dim1]) < small || l2kill
2287                     && c_abs(&a[p + p * a_dim1]) < temp1) {
2288                 goto L3402;
2289             }
2290             ++nr;
2291 /* L3401: */
2292         }
2293 L3402:
2294
2295         ;
2296     } else {
2297 /*        The goal is high relative accuracy. However, if the matrix */
2298 /*        has high scaled condition number the relative accuracy is in */
2299 /*        general not feasible. Later on, a condition number estimator */
2300 /*        will be deployed to estimate the scaled condition number. */
2301 /*        Here we just remove the underflowed part of the triangular */
2302 /*        factor. This prevents the situation in which the code is */
2303 /*        working hard to get the accuracy not warranted by the data. */
2304         temp1 = sqrt(sfmin);
2305         i__1 = *n;
2306         for (p = 2; p <= i__1; ++p) {
2307             if (c_abs(&a[p + p * a_dim1]) < small || l2kill && c_abs(&a[p + p
2308                     * a_dim1]) < temp1) {
2309                 goto L3302;
2310             }
2311             ++nr;
2312 /* L3301: */
2313         }
2314 L3302:
2315
2316         ;
2317     }
2318
2319     almort = FALSE_;
2320     if (nr == *n) {
2321         maxprj = 1.f;
2322         i__1 = *n;
2323         for (p = 2; p <= i__1; ++p) {
2324             temp1 = c_abs(&a[p + p * a_dim1]) / sva[iwork[p]];
2325             maxprj = f2cmin(maxprj,temp1);
2326 /* L3051: */
2327         }
2328 /* Computing 2nd power */
2329         r__1 = maxprj;
2330         if (r__1 * r__1 >= 1.f - (real) (*n) * epsln) {
2331             almort = TRUE_;
2332         }
2333     }
2334
2335
2336     sconda = -1.f;
2337     condr1 = -1.f;
2338     condr2 = -1.f;
2339
2340     if (errest) {
2341         if (*n == nr) {
2342             if (rsvec) {
2343                 clacpy_("U", n, n, &a[a_offset], lda, &v[v_offset], ldv);
2344                 i__1 = *n;
2345                 for (p = 1; p <= i__1; ++p) {
2346                     temp1 = sva[iwork[p]];
2347                     r__1 = 1.f / temp1;
2348                     csscal_(&p, &r__1, &v[p * v_dim1 + 1], &c__1);
2349 /* L3053: */
2350                 }
2351                 if (lsvec) {
2352                     cpocon_("U", n, &v[v_offset], ldv, &c_b141, &temp1, &
2353                             cwork[*n + 1], &rwork[1], &ierr);
2354                 } else {
2355                     cpocon_("U", n, &v[v_offset], ldv, &c_b141, &temp1, &
2356                             cwork[1], &rwork[1], &ierr);
2357                 }
2358
2359             } else if (lsvec) {
2360                 clacpy_("U", n, n, &a[a_offset], lda, &u[u_offset], ldu);
2361                 i__1 = *n;
2362                 for (p = 1; p <= i__1; ++p) {
2363                     temp1 = sva[iwork[p]];
2364                     r__1 = 1.f / temp1;
2365                     csscal_(&p, &r__1, &u[p * u_dim1 + 1], &c__1);
2366 /* L3054: */
2367                 }
2368                 cpocon_("U", n, &u[u_offset], ldu, &c_b141, &temp1, &cwork[*n
2369                         + 1], &rwork[1], &ierr);
2370             } else {
2371                 clacpy_("U", n, n, &a[a_offset], lda, &cwork[1], n)
2372                         ;
2373 /* []            CALL CLACPY( 'U', N, N, A, LDA, CWORK(N+1), N ) */
2374 /*              Change: here index shifted by N to the left, CWORK(1:N) */
2375 /*              not needed for SIGMA only computation */
2376                 i__1 = *n;
2377                 for (p = 1; p <= i__1; ++p) {
2378                     temp1 = sva[iwork[p]];
2379 /* []               CALL CSSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 ) */
2380                     r__1 = 1.f / temp1;
2381                     csscal_(&p, &r__1, &cwork[(p - 1) * *n + 1], &c__1);
2382 /* L3052: */
2383                 }
2384 /* []               CALL CPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1, */
2385 /* []     $              CWORK(N+N*N+1), RWORK, IERR ) */
2386                 cpocon_("U", n, &cwork[1], n, &c_b141, &temp1, &cwork[*n * *n
2387                         + 1], &rwork[1], &ierr);
2388
2389             }
2390             if (temp1 != 0.f) {
2391                 sconda = 1.f / sqrt(temp1);
2392             } else {
2393                 sconda = -1.f;
2394             }
2395 /*           SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). */
2396 /*           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA */
2397         } else {
2398             sconda = -1.f;
2399         }
2400     }
2401
2402     c_div(&q__1, &a[a_dim1 + 1], &a[nr + nr * a_dim1]);
2403     l2pert = l2pert && c_abs(&q__1) > sqrt(big1);
2404 /*     If there is no violent scaling, artificial perturbation is not needed. */
2405
2406 /*     Phase 3: */
2407
2408     if (! (rsvec || lsvec)) {
2409
2410 /*         Singular Values only */
2411
2412 /* Computing MIN */
2413         i__2 = *n - 1;
2414         i__1 = f2cmin(i__2,nr);
2415         for (p = 1; p <= i__1; ++p) {
2416             i__2 = *n - p;
2417             ccopy_(&i__2, &a[p + (p + 1) * a_dim1], lda, &a[p + 1 + p *
2418                     a_dim1], &c__1);
2419             i__2 = *n - p + 1;
2420             clacgv_(&i__2, &a[p + p * a_dim1], &c__1);
2421 /* L1946: */
2422         }
2423         if (nr == *n) {
2424             i__1 = *n + *n * a_dim1;
2425             r_cnjg(&q__1, &a[*n + *n * a_dim1]);
2426             a[i__1].r = q__1.r, a[i__1].i = q__1.i;
2427         }
2428
2429 /*        The following two DO-loops introduce small relative perturbation */
2430 /*        into the strict upper triangle of the lower triangular matrix. */
2431 /*        Small entries below the main diagonal are also changed. */
2432 /*        This modification is useful if the computing environment does not */
2433 /*        provide/allow FLUSH TO ZERO underflow, for it prevents many */
2434 /*        annoying denormalized numbers in case of strongly scaled matrices. */
2435 /*        The perturbation is structured so that it does not introduce any */
2436 /*        new perturbation of the singular values, and it does not destroy */
2437 /*        the job done by the preconditioner. */
2438 /*        The licence for this perturbation is in the variable L2PERT, which */
2439 /*        should be .FALSE. if FLUSH TO ZERO underflow is active. */
2440
2441         if (! almort) {
2442
2443             if (l2pert) {
2444 /*              XSC = SQRT(SMALL) */
2445                 xsc = epsln / (real) (*n);
2446                 i__1 = nr;
2447                 for (q = 1; q <= i__1; ++q) {
2448                     r__1 = xsc * c_abs(&a[q + q * a_dim1]);
2449                     q__1.r = r__1, q__1.i = 0.f;
2450                     ctemp.r = q__1.r, ctemp.i = q__1.i;
2451                     i__2 = *n;
2452                     for (p = 1; p <= i__2; ++p) {
2453                         if (p > q && c_abs(&a[p + q * a_dim1]) <= temp1 || p <
2454                                  q) {
2455                             i__3 = p + q * a_dim1;
2456                             a[i__3].r = ctemp.r, a[i__3].i = ctemp.i;
2457                         }
2458 /*     $                     A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) */
2459 /* L4949: */
2460                     }
2461 /* L4947: */
2462                 }
2463             } else {
2464                 i__1 = nr - 1;
2465                 i__2 = nr - 1;
2466                 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
2467                         , lda);
2468             }
2469
2470
2471             i__1 = *lwork - *n;
2472             cgeqrf_(n, &nr, &a[a_offset], lda, &cwork[1], &cwork[*n + 1], &
2473                     i__1, &ierr);
2474
2475             i__1 = nr - 1;
2476             for (p = 1; p <= i__1; ++p) {
2477                 i__2 = nr - p;
2478                 ccopy_(&i__2, &a[p + (p + 1) * a_dim1], lda, &a[p + 1 + p *
2479                         a_dim1], &c__1);
2480                 i__2 = nr - p + 1;
2481                 clacgv_(&i__2, &a[p + p * a_dim1], &c__1);
2482 /* L1948: */
2483             }
2484
2485         }
2486
2487 /*           Row-cyclic Jacobi SVD algorithm with column pivoting */
2488
2489 /*           to drown denormals */
2490         if (l2pert) {
2491 /*              XSC = SQRT(SMALL) */
2492             xsc = epsln / (real) (*n);
2493             i__1 = nr;
2494             for (q = 1; q <= i__1; ++q) {
2495                 r__1 = xsc * c_abs(&a[q + q * a_dim1]);
2496                 q__1.r = r__1, q__1.i = 0.f;
2497                 ctemp.r = q__1.r, ctemp.i = q__1.i;
2498                 i__2 = nr;
2499                 for (p = 1; p <= i__2; ++p) {
2500                     if (p > q && c_abs(&a[p + q * a_dim1]) <= temp1 || p < q)
2501                             {
2502                         i__3 = p + q * a_dim1;
2503                         a[i__3].r = ctemp.r, a[i__3].i = ctemp.i;
2504                     }
2505 /*     $                   A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) */
2506 /* L1949: */
2507                 }
2508 /* L1947: */
2509             }
2510         } else {
2511             i__1 = nr - 1;
2512             i__2 = nr - 1;
2513             claset_("U", &i__1, &i__2, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1],
2514                     lda);
2515         }
2516
2517 /*           triangular matrix (plus perturbation which is ignored in */
2518 /*           the part which destroys triangular form (confusing?!)) */
2519
2520         cgesvj_("L", "N", "N", &nr, &nr, &a[a_offset], lda, &sva[1], n, &v[
2521                 v_offset], ldv, &cwork[1], lwork, &rwork[1], lrwork, info);
2522
2523         scalem = rwork[1];
2524         numrank = i_nint(&rwork[2]);
2525
2526
2527     } else if (rsvec && ! lsvec && ! jracc || jracc && ! lsvec && nr != *n) {
2528
2529 /*        -> Singular Values and Right Singular Vectors <- */
2530
2531         if (almort) {
2532
2533             i__1 = nr;
2534             for (p = 1; p <= i__1; ++p) {
2535                 i__2 = *n - p + 1;
2536                 ccopy_(&i__2, &a[p + p * a_dim1], lda, &v[p + p * v_dim1], &
2537                         c__1);
2538                 i__2 = *n - p + 1;
2539                 clacgv_(&i__2, &v[p + p * v_dim1], &c__1);
2540 /* L1998: */
2541             }
2542             i__1 = nr - 1;
2543             i__2 = nr - 1;
2544             claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1) + 1],
2545                     ldv);
2546
2547             cgesvj_("L", "U", "N", n, &nr, &v[v_offset], ldv, &sva[1], &nr, &
2548                     a[a_offset], lda, &cwork[1], lwork, &rwork[1], lrwork,
2549                     info);
2550             scalem = rwork[1];
2551             numrank = i_nint(&rwork[2]);
2552         } else {
2553
2554 /*        accumulated product of Jacobi rotations, three are perfect ) */
2555
2556             i__1 = nr - 1;
2557             i__2 = nr - 1;
2558             claset_("L", &i__1, &i__2, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
2559             i__1 = *lwork - *n;
2560             cgelqf_(&nr, n, &a[a_offset], lda, &cwork[1], &cwork[*n + 1], &
2561                     i__1, &ierr);
2562             clacpy_("L", &nr, &nr, &a[a_offset], lda, &v[v_offset], ldv);
2563             i__1 = nr - 1;
2564             i__2 = nr - 1;
2565             claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1) + 1],
2566                     ldv);
2567             i__1 = *lwork - (*n << 1);
2568             cgeqrf_(&nr, &nr, &v[v_offset], ldv, &cwork[*n + 1], &cwork[(*n <<
2569                      1) + 1], &i__1, &ierr);
2570             i__1 = nr;
2571             for (p = 1; p <= i__1; ++p) {
2572                 i__2 = nr - p + 1;
2573                 ccopy_(&i__2, &v[p + p * v_dim1], ldv, &v[p + p * v_dim1], &
2574                         c__1);
2575                 i__2 = nr - p + 1;
2576                 clacgv_(&i__2, &v[p + p * v_dim1], &c__1);
2577 /* L8998: */
2578             }
2579             i__1 = nr - 1;
2580             i__2 = nr - 1;
2581             claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1) + 1],
2582                     ldv);
2583
2584             i__1 = *lwork - *n;
2585             cgesvj_("L", "U", "N", &nr, &nr, &v[v_offset], ldv, &sva[1], &nr,
2586                     &u[u_offset], ldu, &cwork[*n + 1], &i__1, &rwork[1],
2587                     lrwork, info);
2588             scalem = rwork[1];
2589             numrank = i_nint(&rwork[2]);
2590             if (nr < *n) {
2591                 i__1 = *n - nr;
2592                 claset_("A", &i__1, &nr, &c_b1, &c_b1, &v[nr + 1 + v_dim1],
2593                         ldv);
2594                 i__1 = *n - nr;
2595                 claset_("A", &nr, &i__1, &c_b1, &c_b1, &v[(nr + 1) * v_dim1 +
2596                         1], ldv);
2597                 i__1 = *n - nr;
2598                 i__2 = *n - nr;
2599                 claset_("A", &i__1, &i__2, &c_b1, &c_b2, &v[nr + 1 + (nr + 1)
2600                         * v_dim1], ldv);
2601             }
2602
2603             i__1 = *lwork - *n;
2604             cunmlq_("L", "C", n, n, &nr, &a[a_offset], lda, &cwork[1], &v[
2605                     v_offset], ldv, &cwork[*n + 1], &i__1, &ierr);
2606
2607         }
2608 /*         DO 8991 p = 1, N */
2609 /*            CALL CCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA ) */
2610 /* 8991    CONTINUE */
2611 /*         CALL CLACPY( 'All', N, N, A, LDA, V, LDV ) */
2612         clapmr_(&c_false, n, n, &v[v_offset], ldv, &iwork[1]);
2613
2614         if (transp) {
2615             clacpy_("A", n, n, &v[v_offset], ldv, &u[u_offset], ldu);
2616         }
2617
2618     } else if (jracc && ! lsvec && nr == *n) {
2619
2620         i__1 = *n - 1;
2621         i__2 = *n - 1;
2622         claset_("L", &i__1, &i__2, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
2623
2624         cgesvj_("U", "N", "V", n, n, &a[a_offset], lda, &sva[1], n, &v[
2625                 v_offset], ldv, &cwork[1], lwork, &rwork[1], lrwork, info);
2626         scalem = rwork[1];
2627         numrank = i_nint(&rwork[2]);
2628         clapmr_(&c_false, n, n, &v[v_offset], ldv, &iwork[1]);
2629
2630     } else if (lsvec && ! rsvec) {
2631
2632
2633 /*        Jacobi rotations in the Jacobi iterations. */
2634         i__1 = nr;
2635         for (p = 1; p <= i__1; ++p) {
2636             i__2 = *n - p + 1;
2637             ccopy_(&i__2, &a[p + p * a_dim1], lda, &u[p + p * u_dim1], &c__1);
2638             i__2 = *n - p + 1;
2639             clacgv_(&i__2, &u[p + p * u_dim1], &c__1);
2640 /* L1965: */
2641         }
2642         i__1 = nr - 1;
2643         i__2 = nr - 1;
2644         claset_("U", &i__1, &i__2, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
2645
2646         i__1 = *lwork - (*n << 1);
2647         cgeqrf_(n, &nr, &u[u_offset], ldu, &cwork[*n + 1], &cwork[(*n << 1) +
2648                 1], &i__1, &ierr);
2649
2650         i__1 = nr - 1;
2651         for (p = 1; p <= i__1; ++p) {
2652             i__2 = nr - p;
2653             ccopy_(&i__2, &u[p + (p + 1) * u_dim1], ldu, &u[p + 1 + p *
2654                     u_dim1], &c__1);
2655             i__2 = *n - p + 1;
2656             clacgv_(&i__2, &u[p + p * u_dim1], &c__1);
2657 /* L1967: */
2658         }
2659         i__1 = nr - 1;
2660         i__2 = nr - 1;
2661         claset_("U", &i__1, &i__2, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1], ldu);
2662
2663         i__1 = *lwork - *n;
2664         cgesvj_("L", "U", "N", &nr, &nr, &u[u_offset], ldu, &sva[1], &nr, &a[
2665                 a_offset], lda, &cwork[*n + 1], &i__1, &rwork[1], lrwork,
2666                 info);
2667         scalem = rwork[1];
2668         numrank = i_nint(&rwork[2]);
2669
2670         if (nr < *m) {
2671             i__1 = *m - nr;
2672             claset_("A", &i__1, &nr, &c_b1, &c_b1, &u[nr + 1 + u_dim1], ldu);
2673             if (nr < n1) {
2674                 i__1 = n1 - nr;
2675                 claset_("A", &nr, &i__1, &c_b1, &c_b1, &u[(nr + 1) * u_dim1 +
2676                         1], ldu);
2677                 i__1 = *m - nr;
2678                 i__2 = n1 - nr;
2679                 claset_("A", &i__1, &i__2, &c_b1, &c_b2, &u[nr + 1 + (nr + 1)
2680                         * u_dim1], ldu);
2681             }
2682         }
2683
2684         i__1 = *lwork - *n;
2685         cunmqr_("L", "N", m, &n1, n, &a[a_offset], lda, &cwork[1], &u[
2686                 u_offset], ldu, &cwork[*n + 1], &i__1, &ierr);
2687
2688         if (rowpiv) {
2689             i__1 = *m - 1;
2690             claswp_(&n1, &u[u_offset], ldu, &c__1, &i__1, &iwork[iwoff + 1], &
2691                     c_n1);
2692         }
2693
2694         i__1 = n1;
2695         for (p = 1; p <= i__1; ++p) {
2696             xsc = 1.f / scnrm2_(m, &u[p * u_dim1 + 1], &c__1);
2697             csscal_(m, &xsc, &u[p * u_dim1 + 1], &c__1);
2698 /* L1974: */
2699         }
2700
2701         if (transp) {
2702             clacpy_("A", n, n, &u[u_offset], ldu, &v[v_offset], ldv);
2703         }
2704
2705     } else {
2706
2707
2708         if (! jracc) {
2709
2710             if (! almort) {
2711
2712 /*           Second Preconditioning Step (QRF [with pivoting]) */
2713 /*           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is */
2714 /*           equivalent to an LQF CALL. Since in many libraries the QRF */
2715 /*           seems to be better optimized than the LQF, we do explicit */
2716 /*           transpose and use the QRF. This is subject to changes in an */
2717 /*           optimized implementation of CGEJSV. */
2718
2719                 i__1 = nr;
2720                 for (p = 1; p <= i__1; ++p) {
2721                     i__2 = *n - p + 1;
2722                     ccopy_(&i__2, &a[p + p * a_dim1], lda, &v[p + p * v_dim1],
2723                              &c__1);
2724                     i__2 = *n - p + 1;
2725                     clacgv_(&i__2, &v[p + p * v_dim1], &c__1);
2726 /* L1968: */
2727                 }
2728
2729 /*           denormals in the second QR factorization, where they are */
2730 /*           as good as zeros. This is done to avoid painfully slow */
2731 /*           computation with denormals. The relative size of the perturbation */
2732 /*           is a parameter that can be changed by the implementer. */
2733 /*           This perturbation device will be obsolete on machines with */
2734 /*           properly implemented arithmetic. */
2735 /*           To switch it off, set L2PERT=.FALSE. To remove it from  the */
2736 /*           code, remove the action under L2PERT=.TRUE., leave the ELSE part. */
2737 /*           The following two loops should be blocked and fused with the */
2738 /*           transposed copy above. */
2739
2740                 if (l2pert) {
2741                     xsc = sqrt(small);
2742                     i__1 = nr;
2743                     for (q = 1; q <= i__1; ++q) {
2744                         r__1 = xsc * c_abs(&v[q + q * v_dim1]);
2745                         q__1.r = r__1, q__1.i = 0.f;
2746                         ctemp.r = q__1.r, ctemp.i = q__1.i;
2747                         i__2 = *n;
2748                         for (p = 1; p <= i__2; ++p) {
2749                             if (p > q && c_abs(&v[p + q * v_dim1]) <= temp1 ||
2750                                      p < q) {
2751                                 i__3 = p + q * v_dim1;
2752                                 v[i__3].r = ctemp.r, v[i__3].i = ctemp.i;
2753                             }
2754 /*     $                   V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) */
2755                             if (p < q) {
2756                                 i__3 = p + q * v_dim1;
2757                                 i__4 = p + q * v_dim1;
2758                                 q__1.r = -v[i__4].r, q__1.i = -v[i__4].i;
2759                                 v[i__3].r = q__1.r, v[i__3].i = q__1.i;
2760                             }
2761 /* L2968: */
2762                         }
2763 /* L2969: */
2764                     }
2765                 } else {
2766                     i__1 = nr - 1;
2767                     i__2 = nr - 1;
2768                     claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1)
2769                             + 1], ldv);
2770                 }
2771
2772 /*           Estimate the row scaled condition number of R1 */
2773 /*           (If R1 is rectangular, N > NR, then the condition number */
2774 /*           of the leading NR x NR submatrix is estimated.) */
2775
2776                 clacpy_("L", &nr, &nr, &v[v_offset], ldv, &cwork[(*n << 1) +
2777                         1], &nr);
2778                 i__1 = nr;
2779                 for (p = 1; p <= i__1; ++p) {
2780                     i__2 = nr - p + 1;
2781                     temp1 = scnrm2_(&i__2, &cwork[(*n << 1) + (p - 1) * nr +
2782                             p], &c__1);
2783                     i__2 = nr - p + 1;
2784                     r__1 = 1.f / temp1;
2785                     csscal_(&i__2, &r__1, &cwork[(*n << 1) + (p - 1) * nr + p]
2786                             , &c__1);
2787 /* L3950: */
2788                 }
2789                 cpocon_("L", &nr, &cwork[(*n << 1) + 1], &nr, &c_b141, &temp1,
2790                          &cwork[(*n << 1) + nr * nr + 1], &rwork[1], &ierr);
2791                 condr1 = 1.f / sqrt(temp1);
2792 /*           R1 is OK for inverse <=> CONDR1 .LT. REAL(N) */
2793 /*           more conservative    <=> CONDR1 .LT. SQRT(REAL(N)) */
2794
2795                 cond_ok__ = sqrt(sqrt((real) nr));
2796 /* [TP]       COND_OK is a tuning parameter. */
2797
2798                 if (condr1 < cond_ok__) {
2799 /*              implementation, this QRF should be implemented as the QRF */
2800 /*              of a lower triangular matrix. */
2801 /*              R1^* = Q2 * R2 */
2802                     i__1 = *lwork - (*n << 1);
2803                     cgeqrf_(n, &nr, &v[v_offset], ldv, &cwork[*n + 1], &cwork[
2804                             (*n << 1) + 1], &i__1, &ierr);
2805
2806                     if (l2pert) {
2807                         xsc = sqrt(small) / epsln;
2808                         i__1 = nr;
2809                         for (p = 2; p <= i__1; ++p) {
2810                             i__2 = p - 1;
2811                             for (q = 1; q <= i__2; ++q) {
2812 /* Computing MIN */
2813                                 r__2 = c_abs(&v[p + p * v_dim1]), r__3 =
2814                                         c_abs(&v[q + q * v_dim1]);
2815                                 r__1 = xsc * f2cmin(r__2,r__3);
2816                                 q__1.r = r__1, q__1.i = 0.f;
2817                                 ctemp.r = q__1.r, ctemp.i = q__1.i;
2818                                 if (c_abs(&v[q + p * v_dim1]) <= temp1) {
2819                                     i__3 = q + p * v_dim1;
2820                                     v[i__3].r = ctemp.r, v[i__3].i = ctemp.i;
2821                                 }
2822 /*     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) */
2823 /* L3958: */
2824                             }
2825 /* L3959: */
2826                         }
2827                     }
2828
2829                     if (nr != *n) {
2830                         clacpy_("A", n, &nr, &v[v_offset], ldv, &cwork[(*n <<
2831                                 1) + 1], n);
2832                     }
2833
2834                     i__1 = nr - 1;
2835                     for (p = 1; p <= i__1; ++p) {
2836                         i__2 = nr - p;
2837                         ccopy_(&i__2, &v[p + (p + 1) * v_dim1], ldv, &v[p + 1
2838                                 + p * v_dim1], &c__1);
2839                         i__2 = nr - p + 1;
2840                         clacgv_(&i__2, &v[p + p * v_dim1], &c__1);
2841 /* L1969: */
2842                     }
2843                     i__1 = nr + nr * v_dim1;
2844                     r_cnjg(&q__1, &v[nr + nr * v_dim1]);
2845                     v[i__1].r = q__1.r, v[i__1].i = q__1.i;
2846
2847                     condr2 = condr1;
2848
2849                 } else {
2850
2851 /*              Note that windowed pivoting would be equally good */
2852 /*              numerically, and more run-time efficient. So, in */
2853 /*              an optimal implementation, the next call to CGEQP3 */
2854 /*              should be replaced with eg. CALL CGEQPX (ACM TOMS #782) */
2855 /*              with properly (carefully) chosen parameters. */
2856
2857 /*              R1^* * P2 = Q2 * R2 */
2858                     i__1 = nr;
2859                     for (p = 1; p <= i__1; ++p) {
2860                         iwork[*n + p] = 0;
2861 /* L3003: */
2862                     }
2863                     i__1 = *lwork - (*n << 1);
2864                     cgeqp3_(n, &nr, &v[v_offset], ldv, &iwork[*n + 1], &cwork[
2865                             *n + 1], &cwork[(*n << 1) + 1], &i__1, &rwork[1],
2866                             &ierr);
2867 /* *               CALL CGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), */
2868 /* *     $              LWORK-2*N, IERR ) */
2869                     if (l2pert) {
2870                         xsc = sqrt(small);
2871                         i__1 = nr;
2872                         for (p = 2; p <= i__1; ++p) {
2873                             i__2 = p - 1;
2874                             for (q = 1; q <= i__2; ++q) {
2875 /* Computing MIN */
2876                                 r__2 = c_abs(&v[p + p * v_dim1]), r__3 =
2877                                         c_abs(&v[q + q * v_dim1]);
2878                                 r__1 = xsc * f2cmin(r__2,r__3);
2879                                 q__1.r = r__1, q__1.i = 0.f;
2880                                 ctemp.r = q__1.r, ctemp.i = q__1.i;
2881                                 if (c_abs(&v[q + p * v_dim1]) <= temp1) {
2882                                     i__3 = q + p * v_dim1;
2883                                     v[i__3].r = ctemp.r, v[i__3].i = ctemp.i;
2884                                 }
2885 /*     $                     V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) */
2886 /* L3968: */
2887                             }
2888 /* L3969: */
2889                         }
2890                     }
2891
2892                     clacpy_("A", n, &nr, &v[v_offset], ldv, &cwork[(*n << 1)
2893                             + 1], n);
2894
2895                     if (l2pert) {
2896                         xsc = sqrt(small);
2897                         i__1 = nr;
2898                         for (p = 2; p <= i__1; ++p) {
2899                             i__2 = p - 1;
2900                             for (q = 1; q <= i__2; ++q) {
2901 /* Computing MIN */
2902                                 r__2 = c_abs(&v[p + p * v_dim1]), r__3 =
2903                                         c_abs(&v[q + q * v_dim1]);
2904                                 r__1 = xsc * f2cmin(r__2,r__3);
2905                                 q__1.r = r__1, q__1.i = 0.f;
2906                                 ctemp.r = q__1.r, ctemp.i = q__1.i;
2907 /*                        V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) ) */
2908                                 i__3 = p + q * v_dim1;
2909                                 q__1.r = -ctemp.r, q__1.i = -ctemp.i;
2910                                 v[i__3].r = q__1.r, v[i__3].i = q__1.i;
2911 /* L8971: */
2912                             }
2913 /* L8970: */
2914                         }
2915                     } else {
2916                         i__1 = nr - 1;
2917                         i__2 = nr - 1;
2918                         claset_("L", &i__1, &i__2, &c_b1, &c_b1, &v[v_dim1 +
2919                                 2], ldv);
2920                     }
2921 /*              Now, compute R2 = L3 * Q3, the LQ factorization. */
2922                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
2923                     cgelqf_(&nr, &nr, &v[v_offset], ldv, &cwork[(*n << 1) + *
2924                             n * nr + 1], &cwork[(*n << 1) + *n * nr + nr + 1],
2925                              &i__1, &ierr);
2926                     clacpy_("L", &nr, &nr, &v[v_offset], ldv, &cwork[(*n << 1)
2927                              + *n * nr + nr + 1], &nr);
2928                     i__1 = nr;
2929                     for (p = 1; p <= i__1; ++p) {
2930                         temp1 = scnrm2_(&p, &cwork[(*n << 1) + *n * nr + nr +
2931                                 p], &nr);
2932                         r__1 = 1.f / temp1;
2933                         csscal_(&p, &r__1, &cwork[(*n << 1) + *n * nr + nr +
2934                                 p], &nr);
2935 /* L4950: */
2936                     }
2937                     cpocon_("L", &nr, &cwork[(*n << 1) + *n * nr + nr + 1], &
2938                             nr, &c_b141, &temp1, &cwork[(*n << 1) + *n * nr +
2939                             nr + nr * nr + 1], &rwork[1], &ierr);
2940                     condr2 = 1.f / sqrt(temp1);
2941
2942
2943                     if (condr2 >= cond_ok__) {
2944 /*                 (this overwrites the copy of R2, as it will not be */
2945 /*                 needed in this branch, but it does not overwritte the */
2946 /*                 Huseholder vectors of Q2.). */
2947                         clacpy_("U", &nr, &nr, &v[v_offset], ldv, &cwork[(*n
2948                                 << 1) + 1], n);
2949 /*                 WORK(2*N+N*NR+1:2*N+N*NR+N) */
2950                     }
2951
2952                 }
2953
2954                 if (l2pert) {
2955                     xsc = sqrt(small);
2956                     i__1 = nr;
2957                     for (q = 2; q <= i__1; ++q) {
2958                         i__2 = q + q * v_dim1;
2959                         q__1.r = xsc * v[i__2].r, q__1.i = xsc * v[i__2].i;
2960                         ctemp.r = q__1.r, ctemp.i = q__1.i;
2961                         i__2 = q - 1;
2962                         for (p = 1; p <= i__2; ++p) {
2963 /*                     V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) ) */
2964                             i__3 = p + q * v_dim1;
2965                             q__1.r = -ctemp.r, q__1.i = -ctemp.i;
2966                             v[i__3].r = q__1.r, v[i__3].i = q__1.i;
2967 /* L4969: */
2968                         }
2969 /* L4968: */
2970                     }
2971                 } else {
2972                     i__1 = nr - 1;
2973                     i__2 = nr - 1;
2974                     claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1)
2975                             + 1], ldv);
2976                 }
2977
2978 /*        Second preconditioning finished; continue with Jacobi SVD */
2979 /*        The input matrix is lower trinagular. */
2980
2981 /*        Recover the right singular vectors as solution of a well */
2982 /*        conditioned triangular matrix equation. */
2983
2984                 if (condr1 < cond_ok__) {
2985
2986                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
2987                     cgesvj_("L", "U", "N", &nr, &nr, &v[v_offset], ldv, &sva[
2988                             1], &nr, &u[u_offset], ldu, &cwork[(*n << 1) + *n
2989                             * nr + nr + 1], &i__1, &rwork[1], lrwork, info);
2990                     scalem = rwork[1];
2991                     numrank = i_nint(&rwork[2]);
2992                     i__1 = nr;
2993                     for (p = 1; p <= i__1; ++p) {
2994                         ccopy_(&nr, &v[p * v_dim1 + 1], &c__1, &u[p * u_dim1
2995                                 + 1], &c__1);
2996                         csscal_(&nr, &sva[p], &v[p * v_dim1 + 1], &c__1);
2997 /* L3970: */
2998                     }
2999
3000                     if (nr == *n) {
3001 /* :))             .. best case, R1 is inverted. The solution of this matrix */
3002 /*                 equation is Q2*V2 = the product of the Jacobi rotations */
3003 /*                 used in CGESVJ, premultiplied with the orthogonal matrix */
3004 /*                 from the second QR factorization. */
3005                         ctrsm_("L", "U", "N", "N", &nr, &nr, &c_b2, &a[
3006                                 a_offset], lda, &v[v_offset], ldv);
3007                     } else {
3008 /*                 is inverted to get the product of the Jacobi rotations */
3009 /*                 used in CGESVJ. The Q-factor from the second QR */
3010 /*                 factorization is then built in explicitly. */
3011                         ctrsm_("L", "U", "C", "N", &nr, &nr, &c_b2, &cwork[(*
3012                                 n << 1) + 1], n, &v[v_offset], ldv);
3013                         if (nr < *n) {
3014                             i__1 = *n - nr;
3015                             claset_("A", &i__1, &nr, &c_b1, &c_b1, &v[nr + 1
3016                                     + v_dim1], ldv);
3017                             i__1 = *n - nr;
3018                             claset_("A", &nr, &i__1, &c_b1, &c_b1, &v[(nr + 1)
3019                                      * v_dim1 + 1], ldv);
3020                             i__1 = *n - nr;
3021                             i__2 = *n - nr;
3022                             claset_("A", &i__1, &i__2, &c_b1, &c_b2, &v[nr +
3023                                     1 + (nr + 1) * v_dim1], ldv);
3024                         }
3025                         i__1 = *lwork - (*n << 1) - *n * nr - nr;
3026                         cunmqr_("L", "N", n, n, &nr, &cwork[(*n << 1) + 1], n,
3027                                  &cwork[*n + 1], &v[v_offset], ldv, &cwork[(*
3028                                 n << 1) + *n * nr + nr + 1], &i__1, &ierr);
3029                     }
3030
3031                 } else if (condr2 < cond_ok__) {
3032
3033 /*              The matrix R2 is inverted. The solution of the matrix equation */
3034 /*              is Q3^* * V3 = the product of the Jacobi rotations (appplied to */
3035 /*              the lower triangular L3 from the LQ factorization of */
3036 /*              R2=L3*Q3), pre-multiplied with the transposed Q3. */
3037                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
3038                     cgesvj_("L", "U", "N", &nr, &nr, &v[v_offset], ldv, &sva[
3039                             1], &nr, &u[u_offset], ldu, &cwork[(*n << 1) + *n
3040                             * nr + nr + 1], &i__1, &rwork[1], lrwork, info);
3041                     scalem = rwork[1];
3042                     numrank = i_nint(&rwork[2]);
3043                     i__1 = nr;
3044                     for (p = 1; p <= i__1; ++p) {
3045                         ccopy_(&nr, &v[p * v_dim1 + 1], &c__1, &u[p * u_dim1
3046                                 + 1], &c__1);
3047                         csscal_(&nr, &sva[p], &u[p * u_dim1 + 1], &c__1);
3048 /* L3870: */
3049                     }
3050                     ctrsm_("L", "U", "N", "N", &nr, &nr, &c_b2, &cwork[(*n <<
3051                             1) + 1], n, &u[u_offset], ldu);
3052                     i__1 = nr;
3053                     for (q = 1; q <= i__1; ++q) {
3054                         i__2 = nr;
3055                         for (p = 1; p <= i__2; ++p) {
3056                             i__3 = (*n << 1) + *n * nr + nr + iwork[*n + p];
3057                             i__4 = p + q * u_dim1;
3058                             cwork[i__3].r = u[i__4].r, cwork[i__3].i = u[i__4]
3059                                     .i;
3060 /* L872: */
3061                         }
3062                         i__2 = nr;
3063                         for (p = 1; p <= i__2; ++p) {
3064                             i__3 = p + q * u_dim1;
3065                             i__4 = (*n << 1) + *n * nr + nr + p;
3066                             u[i__3].r = cwork[i__4].r, u[i__3].i = cwork[i__4]
3067                                     .i;
3068 /* L874: */
3069                         }
3070 /* L873: */
3071                     }
3072                     if (nr < *n) {
3073                         i__1 = *n - nr;
3074                         claset_("A", &i__1, &nr, &c_b1, &c_b1, &v[nr + 1 +
3075                                 v_dim1], ldv);
3076                         i__1 = *n - nr;
3077                         claset_("A", &nr, &i__1, &c_b1, &c_b1, &v[(nr + 1) *
3078                                 v_dim1 + 1], ldv);
3079                         i__1 = *n - nr;
3080                         i__2 = *n - nr;
3081                         claset_("A", &i__1, &i__2, &c_b1, &c_b2, &v[nr + 1 + (
3082                                 nr + 1) * v_dim1], ldv);
3083                     }
3084                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
3085                     cunmqr_("L", "N", n, n, &nr, &cwork[(*n << 1) + 1], n, &
3086                             cwork[*n + 1], &v[v_offset], ldv, &cwork[(*n << 1)
3087                              + *n * nr + nr + 1], &i__1, &ierr);
3088                 } else {
3089 /*              Last line of defense. */
3090 /* #:(          This is a rather pathological case: no scaled condition */
3091 /*              improvement after two pivoted QR factorizations. Other */
3092 /*              possibility is that the rank revealing QR factorization */
3093 /*              or the condition estimator has failed, or the COND_OK */
3094 /*              is set very close to ONE (which is unnecessary). Normally, */
3095 /*              this branch should never be executed, but in rare cases of */
3096 /*              failure of the RRQR or condition estimator, the last line of */
3097 /*              defense ensures that CGEJSV completes the task. */
3098 /*              Compute the full SVD of L3 using CGESVJ with explicit */
3099 /*              accumulation of Jacobi rotations. */
3100                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
3101                     cgesvj_("L", "U", "V", &nr, &nr, &v[v_offset], ldv, &sva[
3102                             1], &nr, &u[u_offset], ldu, &cwork[(*n << 1) + *n
3103                             * nr + nr + 1], &i__1, &rwork[1], lrwork, info);
3104                     scalem = rwork[1];
3105                     numrank = i_nint(&rwork[2]);
3106                     if (nr < *n) {
3107                         i__1 = *n - nr;
3108                         claset_("A", &i__1, &nr, &c_b1, &c_b1, &v[nr + 1 +
3109                                 v_dim1], ldv);
3110                         i__1 = *n - nr;
3111                         claset_("A", &nr, &i__1, &c_b1, &c_b1, &v[(nr + 1) *
3112                                 v_dim1 + 1], ldv);
3113                         i__1 = *n - nr;
3114                         i__2 = *n - nr;
3115                         claset_("A", &i__1, &i__2, &c_b1, &c_b2, &v[nr + 1 + (
3116                                 nr + 1) * v_dim1], ldv);
3117                     }
3118                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
3119                     cunmqr_("L", "N", n, n, &nr, &cwork[(*n << 1) + 1], n, &
3120                             cwork[*n + 1], &v[v_offset], ldv, &cwork[(*n << 1)
3121                              + *n * nr + nr + 1], &i__1, &ierr);
3122
3123                     i__1 = *lwork - (*n << 1) - *n * nr - nr;
3124                     cunmlq_("L", "C", &nr, &nr, &nr, &cwork[(*n << 1) + 1], n,
3125                              &cwork[(*n << 1) + *n * nr + 1], &u[u_offset],
3126                             ldu, &cwork[(*n << 1) + *n * nr + nr + 1], &i__1,
3127                             &ierr);
3128                     i__1 = nr;
3129                     for (q = 1; q <= i__1; ++q) {
3130                         i__2 = nr;
3131                         for (p = 1; p <= i__2; ++p) {
3132                             i__3 = (*n << 1) + *n * nr + nr + iwork[*n + p];
3133                             i__4 = p + q * u_dim1;
3134                             cwork[i__3].r = u[i__4].r, cwork[i__3].i = u[i__4]
3135                                     .i;
3136 /* L772: */
3137                         }
3138                         i__2 = nr;
3139                         for (p = 1; p <= i__2; ++p) {
3140                             i__3 = p + q * u_dim1;
3141                             i__4 = (*n << 1) + *n * nr + nr + p;
3142                             u[i__3].r = cwork[i__4].r, u[i__3].i = cwork[i__4]
3143                                     .i;
3144 /* L774: */
3145                         }
3146 /* L773: */
3147                     }
3148
3149                 }
3150
3151 /*           Permute the rows of V using the (column) permutation from the */
3152 /*           first QRF. Also, scale the columns to make them unit in */
3153 /*           Euclidean norm. This applies to all cases. */
3154
3155                 temp1 = sqrt((real) (*n)) * epsln;
3156                 i__1 = *n;
3157                 for (q = 1; q <= i__1; ++q) {
3158                     i__2 = *n;
3159                     for (p = 1; p <= i__2; ++p) {
3160                         i__3 = (*n << 1) + *n * nr + nr + iwork[p];
3161                         i__4 = p + q * v_dim1;
3162                         cwork[i__3].r = v[i__4].r, cwork[i__3].i = v[i__4].i;
3163 /* L972: */
3164                     }
3165                     i__2 = *n;
3166                     for (p = 1; p <= i__2; ++p) {
3167                         i__3 = p + q * v_dim1;
3168                         i__4 = (*n << 1) + *n * nr + nr + p;
3169                         v[i__3].r = cwork[i__4].r, v[i__3].i = cwork[i__4].i;
3170 /* L973: */
3171                     }
3172                     xsc = 1.f / scnrm2_(n, &v[q * v_dim1 + 1], &c__1);
3173                     if (xsc < 1.f - temp1 || xsc > temp1 + 1.f) {
3174                         csscal_(n, &xsc, &v[q * v_dim1 + 1], &c__1);
3175                     }
3176 /* L1972: */
3177                 }
3178 /*           At this moment, V contains the right singular vectors of A. */
3179 /*           Next, assemble the left singular vector matrix U (M x N). */
3180                 if (nr < *m) {
3181                     i__1 = *m - nr;
3182                     claset_("A", &i__1, &nr, &c_b1, &c_b1, &u[nr + 1 + u_dim1]
3183                             , ldu);
3184                     if (nr < n1) {
3185                         i__1 = n1 - nr;
3186                         claset_("A", &nr, &i__1, &c_b1, &c_b1, &u[(nr + 1) *
3187                                 u_dim1 + 1], ldu);
3188                         i__1 = *m - nr;
3189                         i__2 = n1 - nr;
3190                         claset_("A", &i__1, &i__2, &c_b1, &c_b2, &u[nr + 1 + (
3191                                 nr + 1) * u_dim1], ldu);
3192                     }
3193                 }
3194
3195 /*           The Q matrix from the first QRF is built into the left singular */
3196 /*           matrix U. This applies to all cases. */
3197
3198                 i__1 = *lwork - *n;
3199                 cunmqr_("L", "N", m, &n1, n, &a[a_offset], lda, &cwork[1], &u[
3200                         u_offset], ldu, &cwork[*n + 1], &i__1, &ierr);
3201 /*           The columns of U are normalized. The cost is O(M*N) flops. */
3202                 temp1 = sqrt((real) (*m)) * epsln;
3203                 i__1 = nr;
3204                 for (p = 1; p <= i__1; ++p) {
3205                     xsc = 1.f / scnrm2_(m, &u[p * u_dim1 + 1], &c__1);
3206                     if (xsc < 1.f - temp1 || xsc > temp1 + 1.f) {
3207                         csscal_(m, &xsc, &u[p * u_dim1 + 1], &c__1);
3208                     }
3209 /* L1973: */
3210                 }
3211
3212 /*           If the initial QRF is computed with row pivoting, the left */
3213 /*           singular vectors must be adjusted. */
3214
3215                 if (rowpiv) {
3216                     i__1 = *m - 1;
3217                     claswp_(&n1, &u[u_offset], ldu, &c__1, &i__1, &iwork[
3218                             iwoff + 1], &c_n1);
3219                 }
3220
3221             } else {
3222
3223 /*        the second QRF is not needed */
3224
3225                 clacpy_("U", n, n, &a[a_offset], lda, &cwork[*n + 1], n);
3226                 if (l2pert) {
3227                     xsc = sqrt(small);
3228                     i__1 = *n;
3229                     for (p = 2; p <= i__1; ++p) {
3230                         i__2 = *n + (p - 1) * *n + p;
3231                         q__1.r = xsc * cwork[i__2].r, q__1.i = xsc * cwork[
3232                                 i__2].i;
3233                         ctemp.r = q__1.r, ctemp.i = q__1.i;
3234                         i__2 = p - 1;
3235                         for (q = 1; q <= i__2; ++q) {
3236 /*                     CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) / */
3237 /*     $                                        ABS(CWORK(N+(p-1)*N+q)) ) */
3238                             i__3 = *n + (q - 1) * *n + p;
3239                             q__1.r = -ctemp.r, q__1.i = -ctemp.i;
3240                             cwork[i__3].r = q__1.r, cwork[i__3].i = q__1.i;
3241 /* L5971: */
3242                         }
3243 /* L5970: */
3244                     }
3245                 } else {
3246                     i__1 = *n - 1;
3247                     i__2 = *n - 1;
3248                     claset_("L", &i__1, &i__2, &c_b1, &c_b1, &cwork[*n + 2],
3249                             n);
3250                 }
3251
3252                 i__1 = *lwork - *n - *n * *n;
3253                 cgesvj_("U", "U", "N", n, n, &cwork[*n + 1], n, &sva[1], n, &
3254                         u[u_offset], ldu, &cwork[*n + *n * *n + 1], &i__1, &
3255                         rwork[1], lrwork, info);
3256
3257                 scalem = rwork[1];
3258                 numrank = i_nint(&rwork[2]);
3259                 i__1 = *n;
3260                 for (p = 1; p <= i__1; ++p) {
3261                     ccopy_(n, &cwork[*n + (p - 1) * *n + 1], &c__1, &u[p *
3262                             u_dim1 + 1], &c__1);
3263                     csscal_(n, &sva[p], &cwork[*n + (p - 1) * *n + 1], &c__1);
3264 /* L6970: */
3265                 }
3266
3267                 ctrsm_("L", "U", "N", "N", n, n, &c_b2, &a[a_offset], lda, &
3268                         cwork[*n + 1], n);
3269                 i__1 = *n;
3270                 for (p = 1; p <= i__1; ++p) {
3271                     ccopy_(n, &cwork[*n + p], n, &v[iwork[p] + v_dim1], ldv);
3272 /* L6972: */
3273                 }
3274                 temp1 = sqrt((real) (*n)) * epsln;
3275                 i__1 = *n;
3276                 for (p = 1; p <= i__1; ++p) {
3277                     xsc = 1.f / scnrm2_(n, &v[p * v_dim1 + 1], &c__1);
3278                     if (xsc < 1.f - temp1 || xsc > temp1 + 1.f) {
3279                         csscal_(n, &xsc, &v[p * v_dim1 + 1], &c__1);
3280                     }
3281 /* L6971: */
3282                 }
3283
3284 /*           Assemble the left singular vector matrix U (M x N). */
3285
3286                 if (*n < *m) {
3287                     i__1 = *m - *n;
3288                     claset_("A", &i__1, n, &c_b1, &c_b1, &u[*n + 1 + u_dim1],
3289                             ldu);
3290                     if (*n < n1) {
3291                         i__1 = n1 - *n;
3292                         claset_("A", n, &i__1, &c_b1, &c_b1, &u[(*n + 1) *
3293                                 u_dim1 + 1], ldu);
3294                         i__1 = *m - *n;
3295                         i__2 = n1 - *n;
3296                         claset_("A", &i__1, &i__2, &c_b1, &c_b2, &u[*n + 1 + (
3297                                 *n + 1) * u_dim1], ldu);
3298                     }
3299                 }
3300                 i__1 = *lwork - *n;
3301                 cunmqr_("L", "N", m, &n1, n, &a[a_offset], lda, &cwork[1], &u[
3302                         u_offset], ldu, &cwork[*n + 1], &i__1, &ierr);
3303                 temp1 = sqrt((real) (*m)) * epsln;
3304                 i__1 = n1;
3305                 for (p = 1; p <= i__1; ++p) {
3306                     xsc = 1.f / scnrm2_(m, &u[p * u_dim1 + 1], &c__1);
3307                     if (xsc < 1.f - temp1 || xsc > temp1 + 1.f) {
3308                         csscal_(m, &xsc, &u[p * u_dim1 + 1], &c__1);
3309                     }
3310 /* L6973: */
3311                 }
3312
3313                 if (rowpiv) {
3314                     i__1 = *m - 1;
3315                     claswp_(&n1, &u[u_offset], ldu, &c__1, &i__1, &iwork[
3316                             iwoff + 1], &c_n1);
3317                 }
3318
3319             }
3320
3321 /*        end of the  >> almost orthogonal case <<  in the full SVD */
3322
3323         } else {
3324
3325 /*        This branch deploys a preconditioned Jacobi SVD with explicitly */
3326 /*        accumulated rotations. It is included as optional, mainly for */
3327 /*        experimental purposes. It does perform well, and can also be used. */
3328 /*        In this implementation, this branch will be automatically activated */
3329 /*        if the  condition number sigma_max(A) / sigma_min(A) is predicted */
3330 /*        to be greater than the overflow threshold. This is because the */
3331 /*        a posteriori computation of the singular vectors assumes robust */
3332 /*        implementation of BLAS and some LAPACK procedures, capable of working */
3333 /*        in presence of extreme values, e.g. when the singular values spread from */
3334 /*        the underflow to the overflow threshold. */
3335
3336             i__1 = nr;
3337             for (p = 1; p <= i__1; ++p) {
3338                 i__2 = *n - p + 1;
3339                 ccopy_(&i__2, &a[p + p * a_dim1], lda, &v[p + p * v_dim1], &
3340                         c__1);
3341                 i__2 = *n - p + 1;
3342                 clacgv_(&i__2, &v[p + p * v_dim1], &c__1);
3343 /* L7968: */
3344             }
3345
3346             if (l2pert) {
3347                 xsc = sqrt(small / epsln);
3348                 i__1 = nr;
3349                 for (q = 1; q <= i__1; ++q) {
3350                     r__1 = xsc * c_abs(&v[q + q * v_dim1]);
3351                     q__1.r = r__1, q__1.i = 0.f;
3352                     ctemp.r = q__1.r, ctemp.i = q__1.i;
3353                     i__2 = *n;
3354                     for (p = 1; p <= i__2; ++p) {
3355                         if (p > q && c_abs(&v[p + q * v_dim1]) <= temp1 || p <
3356                                  q) {
3357                             i__3 = p + q * v_dim1;
3358                             v[i__3].r = ctemp.r, v[i__3].i = ctemp.i;
3359                         }
3360 /*     $                V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) */
3361                         if (p < q) {
3362                             i__3 = p + q * v_dim1;
3363                             i__4 = p + q * v_dim1;
3364                             q__1.r = -v[i__4].r, q__1.i = -v[i__4].i;
3365                             v[i__3].r = q__1.r, v[i__3].i = q__1.i;
3366                         }
3367 /* L5968: */
3368                     }
3369 /* L5969: */
3370                 }
3371             } else {
3372                 i__1 = nr - 1;
3373                 i__2 = nr - 1;
3374                 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &v[(v_dim1 << 1) + 1]
3375                         , ldv);
3376             }
3377             i__1 = *lwork - (*n << 1);
3378             cgeqrf_(n, &nr, &v[v_offset], ldv, &cwork[*n + 1], &cwork[(*n <<
3379                     1) + 1], &i__1, &ierr);
3380             clacpy_("L", n, &nr, &v[v_offset], ldv, &cwork[(*n << 1) + 1], n);
3381
3382             i__1 = nr;
3383             for (p = 1; p <= i__1; ++p) {
3384                 i__2 = nr - p + 1;
3385                 ccopy_(&i__2, &v[p + p * v_dim1], ldv, &u[p + p * u_dim1], &
3386                         c__1);
3387                 i__2 = nr - p + 1;
3388                 clacgv_(&i__2, &u[p + p * u_dim1], &c__1);
3389 /* L7969: */
3390             }
3391             if (l2pert) {
3392                 xsc = sqrt(small / epsln);
3393                 i__1 = nr;
3394                 for (q = 2; q <= i__1; ++q) {
3395                     i__2 = q - 1;
3396                     for (p = 1; p <= i__2; ++p) {
3397 /* Computing MIN */
3398                         r__2 = c_abs(&u[p + p * u_dim1]), r__3 = c_abs(&u[q +
3399                                 q * u_dim1]);
3400                         r__1 = xsc * f2cmin(r__2,r__3);
3401                         q__1.r = r__1, q__1.i = 0.f;
3402                         ctemp.r = q__1.r, ctemp.i = q__1.i;
3403 /*                  U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) ) */
3404                         i__3 = p + q * u_dim1;
3405                         q__1.r = -ctemp.r, q__1.i = -ctemp.i;
3406                         u[i__3].r = q__1.r, u[i__3].i = q__1.i;
3407 /* L9971: */
3408                     }
3409 /* L9970: */
3410                 }
3411             } else {
3412                 i__1 = nr - 1;
3413                 i__2 = nr - 1;
3414                 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &u[(u_dim1 << 1) + 1]
3415                         , ldu);
3416             }
3417             i__1 = *lwork - (*n << 1) - *n * nr;
3418             cgesvj_("L", "U", "V", &nr, &nr, &u[u_offset], ldu, &sva[1], n, &
3419                     v[v_offset], ldv, &cwork[(*n << 1) + *n * nr + 1], &i__1,
3420                     &rwork[1], lrwork, info);
3421             scalem = rwork[1];
3422             numrank = i_nint(&rwork[2]);
3423             if (nr < *n) {
3424                 i__1 = *n - nr;
3425                 claset_("A", &i__1, &nr, &c_b1, &c_b1, &v[nr + 1 + v_dim1],
3426                         ldv);
3427                 i__1 = *n - nr;
3428                 claset_("A", &nr, &i__1, &c_b1, &c_b1, &v[(nr + 1) * v_dim1 +
3429                         1], ldv);
3430                 i__1 = *n - nr;
3431                 i__2 = *n - nr;
3432                 claset_("A", &i__1, &i__2, &c_b1, &c_b2, &v[nr + 1 + (nr + 1)
3433                         * v_dim1], ldv);
3434             }
3435             i__1 = *lwork - (*n << 1) - *n * nr - nr;
3436             cunmqr_("L", "N", n, n, &nr, &cwork[(*n << 1) + 1], n, &cwork[*n
3437                     + 1], &v[v_offset], ldv, &cwork[(*n << 1) + *n * nr + nr
3438                     + 1], &i__1, &ierr);
3439
3440 /*           Permute the rows of V using the (column) permutation from the */
3441 /*           first QRF. Also, scale the columns to make them unit in */
3442 /*           Euclidean norm. This applies to all cases. */
3443
3444             temp1 = sqrt((real) (*n)) * epsln;
3445             i__1 = *n;
3446             for (q = 1; q <= i__1; ++q) {
3447                 i__2 = *n;
3448                 for (p = 1; p <= i__2; ++p) {
3449                     i__3 = (*n << 1) + *n * nr + nr + iwork[p];
3450                     i__4 = p + q * v_dim1;
3451                     cwork[i__3].r = v[i__4].r, cwork[i__3].i = v[i__4].i;
3452 /* L8972: */
3453                 }
3454                 i__2 = *n;
3455                 for (p = 1; p <= i__2; ++p) {
3456                     i__3 = p + q * v_dim1;
3457                     i__4 = (*n << 1) + *n * nr + nr + p;
3458                     v[i__3].r = cwork[i__4].r, v[i__3].i = cwork[i__4].i;
3459 /* L8973: */
3460                 }
3461                 xsc = 1.f / scnrm2_(n, &v[q * v_dim1 + 1], &c__1);
3462                 if (xsc < 1.f - temp1 || xsc > temp1 + 1.f) {
3463                     csscal_(n, &xsc, &v[q * v_dim1 + 1], &c__1);
3464                 }
3465 /* L7972: */
3466             }
3467
3468 /*           At this moment, V contains the right singular vectors of A. */
3469 /*           Next, assemble the left singular vector matrix U (M x N). */
3470
3471             if (nr < *m) {
3472                 i__1 = *m - nr;
3473                 claset_("A", &i__1, &nr, &c_b1, &c_b1, &u[nr + 1 + u_dim1],
3474                         ldu);
3475                 if (nr < n1) {
3476                     i__1 = n1 - nr;
3477                     claset_("A", &nr, &i__1, &c_b1, &c_b1, &u[(nr + 1) *
3478                             u_dim1 + 1], ldu);
3479                     i__1 = *m - nr;
3480                     i__2 = n1 - nr;
3481                     claset_("A", &i__1, &i__2, &c_b1, &c_b2, &u[nr + 1 + (nr
3482                             + 1) * u_dim1], ldu);
3483                 }
3484             }
3485
3486             i__1 = *lwork - *n;
3487             cunmqr_("L", "N", m, &n1, n, &a[a_offset], lda, &cwork[1], &u[
3488                     u_offset], ldu, &cwork[*n + 1], &i__1, &ierr);
3489
3490             if (rowpiv) {
3491                 i__1 = *m - 1;
3492                 claswp_(&n1, &u[u_offset], ldu, &c__1, &i__1, &iwork[iwoff +
3493                         1], &c_n1);
3494             }
3495
3496
3497         }
3498         if (transp) {
3499             i__1 = *n;
3500             for (p = 1; p <= i__1; ++p) {
3501                 cswap_(n, &u[p * u_dim1 + 1], &c__1, &v[p * v_dim1 + 1], &
3502                         c__1);
3503 /* L6974: */
3504             }
3505         }
3506
3507     }
3508 /*     end of the full SVD */
3509
3510 /*     Undo scaling, if necessary (and possible) */
3511
3512     if (uscal2 <= big / sva[1] * uscal1) {
3513         slascl_("G", &c__0, &c__0, &uscal1, &uscal2, &nr, &c__1, &sva[1], n, &
3514                 ierr);
3515         uscal1 = 1.f;
3516         uscal2 = 1.f;
3517     }
3518
3519     if (nr < *n) {
3520         i__1 = *n;
3521         for (p = nr + 1; p <= i__1; ++p) {
3522             sva[p] = 0.f;
3523 /* L3004: */
3524         }
3525     }
3526
3527     rwork[1] = uscal2 * scalem;
3528     rwork[2] = uscal1;
3529     if (errest) {
3530         rwork[3] = sconda;
3531     }
3532     if (lsvec && rsvec) {
3533         rwork[4] = condr1;
3534         rwork[5] = condr2;
3535     }
3536     if (l2tran) {
3537         rwork[6] = entra;
3538         rwork[7] = entrat;
3539     }
3540
3541     iwork[1] = nr;
3542     iwork[2] = numrank;
3543     iwork[3] = warning;
3544     if (transp) {
3545         iwork[4] = 1;
3546     } else {
3547         iwork[4] = -1;
3548     }
3549
3550     return 0;
3551 } /* cgejsv_ */
3552