lapack-netlib/SRC/zgesdd.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 /* Table of constant values */
 514
 515 static doublecomplex c_b1 = {0.,0.};
 516 static doublecomplex c_b2 = {1.,0.};
 517 static integer c_n1 = -1;
 518 static integer c__0 = 0;
 519 static integer c__1 = 1;
 520
 521 /* > \brief \b ZGESDD */
 522
 523 /*  =========== DOCUMENTATION =========== */
 524
 525 /* Online html documentation available at */
 526 /*            http://www.netlib.org/lapack/explore-html/ */
 527
 528 /* > \htmlonly */
 529 /* > Download ZGESDD + dependencies */
 530 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.
 531 f"> */
 532 /* > [TGZ]</a> */
 533 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.
 534 f"> */
 535 /* > [ZIP]</a> */
 536 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.
 537 f"> */
 538 /* > [TXT]</a> */
 539 /* > \endhtmlonly */
 540
 541 /*  Definition: */
 542 /*  =========== */
 543
 544 /*       SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, */
 545 /*                          WORK, LWORK, RWORK, IWORK, INFO ) */
 546
 547 /*       CHARACTER          JOBZ */
 548 /*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N */
 549 /*       INTEGER            IWORK( * ) */
 550 /*       DOUBLE PRECISION   RWORK( * ), S( * ) */
 551 /*       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ), */
 552 /*      $                   WORK( * ) */
 553
 554
 555 /* > \par Purpose: */
 556 /*  ============= */
 557 /* > */
 558 /* > \verbatim */
 559 /* > */
 560 /* > ZGESDD computes the singular value decomposition (SVD) of a complex */
 561 /* > M-by-N matrix A, optionally computing the left and/or right singular */
 562 /* > vectors, by using divide-and-conquer method. The SVD is written */
 563 /* > */
 564 /* >      A = U * SIGMA * conjugate-transpose(V) */
 565 /* > */
 566 /* > where SIGMA is an M-by-N matrix which is zero except for its */
 567 /* > f2cmin(m,n) diagonal elements, U is an M-by-M unitary matrix, and */
 568 /* > V is an N-by-N unitary matrix.  The diagonal elements of SIGMA */
 569 /* > are the singular values of A; they are real and non-negative, and */
 570 /* > are returned in descending order.  The first f2cmin(m,n) columns of */
 571 /* > U and V are the left and right singular vectors of A. */
 572 /* > */
 573 /* > Note that the routine returns VT = V**H, not V. */
 574 /* > */
 575 /* > The divide and conquer algorithm makes very mild assumptions about */
 576 /* > floating point arithmetic. It will work on machines with a guard */
 577 /* > digit in add/subtract, or on those binary machines without guard */
 578 /* > digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
 579 /* > Cray-2. It could conceivably fail on hexadecimal or decimal machines */
 580 /* > without guard digits, but we know of none. */
 581 /* > \endverbatim */
 582
 583 /*  Arguments: */
 584 /*  ========== */
 585
 586 /* > \param[in] JOBZ */
 587 /* > \verbatim */
 588 /* >          JOBZ is CHARACTER*1 */
 589 /* >          Specifies options for computing all or part of the matrix U: */
 590 /* >          = 'A':  all M columns of U and all N rows of V**H are */
 591 /* >                  returned in the arrays U and VT; */
 592 /* >          = 'S':  the first f2cmin(M,N) columns of U and the first */
 593 /* >                  f2cmin(M,N) rows of V**H are returned in the arrays U */
 594 /* >                  and VT; */
 595 /* >          = 'O':  If M >= N, the first N columns of U are overwritten */
 596 /* >                  in the array A and all rows of V**H are returned in */
 597 /* >                  the array VT; */
 598 /* >                  otherwise, all columns of U are returned in the */
 599 /* >                  array U and the first M rows of V**H are overwritten */
 600 /* >                  in the array A; */
 601 /* >          = 'N':  no columns of U or rows of V**H are computed. */
 602 /* > \endverbatim */
 603 /* > */
 604 /* > \param[in] M */
 605 /* > \verbatim */
 606 /* >          M is INTEGER */
 607 /* >          The number of rows of the input matrix A.  M >= 0. */
 608 /* > \endverbatim */
 609 /* > */
 610 /* > \param[in] N */
 611 /* > \verbatim */
 612 /* >          N is INTEGER */
 613 /* >          The number of columns of the input matrix A.  N >= 0. */
 614 /* > \endverbatim */
 615 /* > */
 616 /* > \param[in,out] A */
 617 /* > \verbatim */
 618 /* >          A is COMPLEX*16 array, dimension (LDA,N) */
 619 /* >          On entry, the M-by-N matrix A. */
 620 /* >          On exit, */
 621 /* >          if JOBZ = 'O',  A is overwritten with the first N columns */
 622 /* >                          of U (the left singular vectors, stored */
 623 /* >                          columnwise) if M >= N; */
 624 /* >                          A is overwritten with the first M rows */
 625 /* >                          of V**H (the right singular vectors, stored */
 626 /* >                          rowwise) otherwise. */
 627 /* >          if JOBZ .ne. 'O', the contents of A are destroyed. */
 628 /* > \endverbatim */
 629 /* > */
 630 /* > \param[in] LDA */
 631 /* > \verbatim */
 632 /* >          LDA is INTEGER */
 633 /* >          The leading dimension of the array A.  LDA >= f2cmax(1,M). */
 634 /* > \endverbatim */
 635 /* > */
 636 /* > \param[out] S */
 637 /* > \verbatim */
 638 /* >          S is DOUBLE PRECISION array, dimension (f2cmin(M,N)) */
 639 /* >          The singular values of A, sorted so that S(i) >= S(i+1). */
 640 /* > \endverbatim */
 641 /* > */
 642 /* > \param[out] U */
 643 /* > \verbatim */
 644 /* >          U is COMPLEX*16 array, dimension (LDU,UCOL) */
 645 /* >          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; */
 646 /* >          UCOL = f2cmin(M,N) if JOBZ = 'S'. */
 647 /* >          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M */
 648 /* >          unitary matrix U; */
 649 /* >          if JOBZ = 'S', U contains the first f2cmin(M,N) columns of U */
 650 /* >          (the left singular vectors, stored columnwise); */
 651 /* >          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. */
 652 /* > \endverbatim */
 653 /* > */
 654 /* > \param[in] LDU */
 655 /* > \verbatim */
 656 /* >          LDU is INTEGER */
 657 /* >          The leading dimension of the array U.  LDU >= 1; */
 658 /* >          if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. */
 659 /* > \endverbatim */
 660 /* > */
 661 /* > \param[out] VT */
 662 /* > \verbatim */
 663 /* >          VT is COMPLEX*16 array, dimension (LDVT,N) */
 664 /* >          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the */
 665 /* >          N-by-N unitary matrix V**H; */
 666 /* >          if JOBZ = 'S', VT contains the first f2cmin(M,N) rows of */
 667 /* >          V**H (the right singular vectors, stored rowwise); */
 668 /* >          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. */
 669 /* > \endverbatim */
 670 /* > */
 671 /* > \param[in] LDVT */
 672 /* > \verbatim */
 673 /* >          LDVT is INTEGER */
 674 /* >          The leading dimension of the array VT.  LDVT >= 1; */
 675 /* >          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; */
 676 /* >          if JOBZ = 'S', LDVT >= f2cmin(M,N). */
 677 /* > \endverbatim */
 678 /* > */
 679 /* > \param[out] WORK */
 680 /* > \verbatim */
 681 /* >          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) */
 682 /* >          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 683 /* > \endverbatim */
 684 /* > */
 685 /* > \param[in] LWORK */
 686 /* > \verbatim */
 687 /* >          LWORK is INTEGER */
 688 /* >          The dimension of the array WORK. LWORK >= 1. */
 689 /* >          If LWORK = -1, a workspace query is assumed.  The optimal */
 690 /* >          size for the WORK array is calculated and stored in WORK(1), */
 691 /* >          and no other work except argument checking is performed. */
 692 /* > */
 693 /* >          Let mx = f2cmax(M,N) and mn = f2cmin(M,N). */
 694 /* >          If JOBZ = 'N', LWORK >= 2*mn + mx. */
 695 /* >          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx. */
 696 /* >          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn. */
 697 /* >          If JOBZ = 'A', LWORK >=   mn*mn + 2*mn + mx. */
 698 /* >          These are not tight minimums in all cases; see comments inside code. */
 699 /* >          For good performance, LWORK should generally be larger; */
 700 /* >          a query is recommended. */
 701 /* > \endverbatim */
 702 /* > */
 703 /* > \param[out] RWORK */
 704 /* > \verbatim */
 705 /* >          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) */
 706 /* >          Let mx = f2cmax(M,N) and mn = f2cmin(M,N). */
 707 /* >          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); */
 708 /* >          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn; */
 709 /* >          else              LRWORK >= f2cmax( 5*mn*mn + 5*mn, */
 710 /* >                                           2*mx*mn + 2*mn*mn + mn ). */
 711 /* > \endverbatim */
 712 /* > */
 713 /* > \param[out] IWORK */
 714 /* > \verbatim */
 715 /* >          IWORK is INTEGER array, dimension (8*f2cmin(M,N)) */
 716 /* > \endverbatim */
 717 /* > */
 718 /* > \param[out] INFO */
 719 /* > \verbatim */
 720 /* >          INFO is INTEGER */
 721 /* >          = 0:  successful exit. */
 722 /* >          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 723 /* >          > 0:  The updating process of DBDSDC did not converge. */
 724 /* > \endverbatim */
 725
 726 /*  Authors: */
 727 /*  ======== */
 728
 729 /* > \author Univ. of Tennessee */
 730 /* > \author Univ. of California Berkeley */
 731 /* > \author Univ. of Colorado Denver */
 732 /* > \author NAG Ltd. */
 733
 734 /* > \date June 2016 */
 735
 736 /* > \ingroup complex16GEsing */
 737
 738 /* > \par Contributors: */
 739 /*  ================== */
 740 /* > */
 741 /* >     Ming Gu and Huan Ren, Computer Science Division, University of */
 742 /* >     California at Berkeley, USA */
 743 /* > */
 744 /*  ===================================================================== */
 745 /* Subroutine */ int zgesdd_(char *jobz, integer *m, integer *n,
 746         doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
 747         integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
 748         integer *lwork, doublereal *rwork, integer *iwork, integer *info)
 749 {
 750     /* System generated locals */
 751     integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
 752             i__2, i__3;
 753
 754     /* Local variables */
 755     integer lwork_zgebrd_mm__, lwork_zgebrd_mn__, lwork_zgebrd_nn__,
 756             lwork_zgelqf_mn__, lwork_zgeqrf_mn__;
 757     doublecomplex cdum[1];
 758     integer iscl;
 759     doublereal anrm;
 760     integer idum[1], ierr, itau, lwork_zunglq_mn__, lwork_zunglq_nn__,
 761             lwork_zungqr_mm__, lwork_zungqr_mn__, irvt, lwork_zunmbr_prc_mm__,
 762              lwork_zunmbr_prc_mn__, lwork_zunmbr_prc_nn__,
 763             lwork_zunmbr_qln_mm__, lwork_zunmbr_qln_mn__,
 764             lwork_zunmbr_qln_nn__, i__;
 765     extern logical lsame_(char *, char *);
 766     integer chunk, minmn;
 767     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,
 768             integer *, doublecomplex *, doublecomplex *, integer *,
 769             doublecomplex *, integer *, doublecomplex *, doublecomplex *,
 770             integer *);
 771     integer wrkbl, itaup, itauq;
 772     logical wntqa;
 773     integer nwork;
 774     logical wntqn, wntqo, wntqs;
 775     extern /* Subroutine */ int zlacp2_(char *, integer *, integer *,
 776             doublereal *, integer *, doublecomplex *, integer *);
 777     integer mnthr1, mnthr2, ie;
 778     extern /* Subroutine */ int dbdsdc_(char *, char *, integer *, doublereal
 779             *, doublereal *, doublereal *, integer *, doublereal *, integer *,
 780              doublereal *, integer *, doublereal *, integer *, integer *);
 781     integer il;
 782     extern doublereal dlamch_(char *);
 783     integer ir, iu;
 784     extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
 785             doublereal *, doublereal *, integer *, integer *, doublereal *,
 786             integer *, integer *);
 787     integer lwork_zungbr_p_mn__, lwork_zungbr_p_nn__, lwork_zungbr_q_mn__,
 788             lwork_zungbr_q_mm__;
 789     extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
 790     doublereal bignum;
 791     extern /* Subroutine */ int zgebrd_(integer *, integer *, doublecomplex *,
 792              integer *, doublereal *, doublereal *, doublecomplex *,
 793             doublecomplex *, doublecomplex *, integer *, integer *);
 794     extern logical disnan_(doublereal *);
 795     extern doublereal zlange_(char *, integer *, integer *, doublecomplex *,
 796             integer *, doublereal *);
 797     extern /* Subroutine */ int zgelqf_(integer *, integer *, doublecomplex *,
 798              integer *, doublecomplex *, doublecomplex *, integer *, integer *
 799             ), zlacrm_(integer *, integer *, doublecomplex *, integer *,
 800             doublereal *, integer *, doublecomplex *, integer *, doublereal *)
 801             , zlarcm_(integer *, integer *, doublereal *, integer *,
 802             doublecomplex *, integer *, doublecomplex *, integer *,
 803             doublereal *), zlascl_(char *, integer *, integer *, doublereal *,
 804              doublereal *, integer *, integer *, doublecomplex *, integer *,
 805             integer *), zgeqrf_(integer *, integer *, doublecomplex *,
 806              integer *, doublecomplex *, doublecomplex *, integer *, integer *
 807             );
 808     integer ldwrkl;
 809     extern /* Subroutine */ int zlacpy_(char *, integer *, integer *,
 810             doublecomplex *, integer *, doublecomplex *, integer *),
 811             zlaset_(char *, integer *, integer *, doublecomplex *,
 812             doublecomplex *, doublecomplex *, integer *);
 813     integer ldwrkr, minwrk, ldwrku, maxwrk;
 814     extern /* Subroutine */ int zungbr_(char *, integer *, integer *, integer
 815             *, doublecomplex *, integer *, doublecomplex *, doublecomplex *,
 816             integer *, integer *);
 817     integer ldwkvt;
 818     doublereal smlnum;
 819     logical wntqas;
 820     extern /* Subroutine */ int zunmbr_(char *, char *, char *, integer *,
 821             integer *, integer *, doublecomplex *, integer *, doublecomplex *,
 822              doublecomplex *, integer *, doublecomplex *, integer *, integer *
 823             ), zunglq_(integer *, integer *, integer *
 824             , doublecomplex *, integer *, doublecomplex *, doublecomplex *,
 825             integer *, integer *);
 826     logical lquery;
 827     integer nrwork;
 828     extern /* Subroutine */ int zungqr_(integer *, integer *, integer *,
 829             doublecomplex *, integer *, doublecomplex *, doublecomplex *,
 830             integer *, integer *);
 831     integer blk;
 832     doublereal dum[1], eps;
 833     integer iru, ivt;
 834
 835
 836 /*  -- LAPACK driver routine (version 3.7.0) -- */
 837 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
 838 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
 839 /*     June 2016 */
 840
 841
 842 /*  ===================================================================== */
 843
 844
 845 /*     Test the input arguments */
 846
 847     /* Parameter adjustments */
 848     a_dim1 = *lda;
 849     a_offset = 1 + a_dim1 * 1;
 850     a -= a_offset;
 851     --s;
 852     u_dim1 = *ldu;
 853     u_offset = 1 + u_dim1 * 1;
 854     u -= u_offset;
 855     vt_dim1 = *ldvt;
 856     vt_offset = 1 + vt_dim1 * 1;
 857     vt -= vt_offset;
 858     --work;
 859     --rwork;
 860     --iwork;
 861
 862     /* Function Body */
 863     *info = 0;
 864     minmn = f2cmin(*m,*n);
 865     mnthr1 = (integer) (minmn * 17. / 9.);
 866     mnthr2 = (integer) (minmn * 5. / 3.);
 867     wntqa = lsame_(jobz, "A");
 868     wntqs = lsame_(jobz, "S");
 869     wntqas = wntqa || wntqs;
 870     wntqo = lsame_(jobz, "O");
 871     wntqn = lsame_(jobz, "N");
 872     lquery = *lwork == -1;
 873     minwrk = 1;
 874     maxwrk = 1;
 875
 876     if (! (wntqa || wntqs || wntqo || wntqn)) {
 877         *info = -1;
 878     } else if (*m < 0) {
 879         *info = -2;
 880     } else if (*n < 0) {
 881         *info = -3;
 882     } else if (*lda < f2cmax(1,*m)) {
 883         *info = -5;
 884     } else if (*ldu < 1 || wntqas && *ldu < *m || wntqo && *m < *n && *ldu < *
 885             m) {
 886         *info = -8;
 887     } else if (*ldvt < 1 || wntqa && *ldvt < *n || wntqs && *ldvt < minmn ||
 888             wntqo && *m >= *n && *ldvt < *n) {
 889         *info = -10;
 890     }
 891
 892 /*     Compute workspace */
 893 /*       Note: Comments in the code beginning "Workspace:" describe the */
 894 /*       minimal amount of workspace allocated at that point in the code, */
 895 /*       as well as the preferred amount for good performance. */
 896 /*       CWorkspace refers to complex workspace, and RWorkspace to */
 897 /*       real workspace. NB refers to the optimal block size for the */
 898 /*       immediately following subroutine, as returned by ILAENV.) */
 899
 900     if (*info == 0) {
 901         minwrk = 1;
 902         maxwrk = 1;
 903         if (*m >= *n && minmn > 0) {
 904
 905 /*           There is no complex work space needed for bidiagonal SVD */
 906 /*           The real work space needed for bidiagonal SVD (dbdsdc) is */
 907 /*           BDSPAC = 3*N*N + 4*N for singular values and vectors; */
 908 /*           BDSPAC = 4*N         for singular values only; */
 909 /*           not including e, RU, and RVT matrices. */
 910
 911 /*           Compute space preferred for each routine */
 912             zgebrd_(m, n, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
 913             lwork_zgebrd_mn__ = (integer) cdum[0].r;
 914
 915             zgebrd_(n, n, cdum, n, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
 916             lwork_zgebrd_nn__ = (integer) cdum[0].r;
 917
 918             zgeqrf_(m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
 919             lwork_zgeqrf_mn__ = (integer) cdum[0].r;
 920
 921             zungbr_("P", n, n, n, cdum, n, cdum, cdum, &c_n1, &ierr);
 922             lwork_zungbr_p_nn__ = (integer) cdum[0].r;
 923
 924             zungbr_("Q", m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
 925             lwork_zungbr_q_mm__ = (integer) cdum[0].r;
 926
 927             zungbr_("Q", m, n, n, cdum, m, cdum, cdum, &c_n1, &ierr);
 928             lwork_zungbr_q_mn__ = (integer) cdum[0].r;
 929
 930             zungqr_(m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
 931             lwork_zungqr_mm__ = (integer) cdum[0].r;
 932
 933             zungqr_(m, n, n, cdum, m, cdum, cdum, &c_n1, &ierr);
 934             lwork_zungqr_mn__ = (integer) cdum[0].r;
 935
 936             zunmbr_("P", "R", "C", n, n, n, cdum, n, cdum, cdum, n, cdum, &
 937                     c_n1, &ierr);
 938             lwork_zunmbr_prc_nn__ = (integer) cdum[0].r;
 939
 940             zunmbr_("Q", "L", "N", m, m, n, cdum, m, cdum, cdum, m, cdum, &
 941                     c_n1, &ierr);
 942             lwork_zunmbr_qln_mm__ = (integer) cdum[0].r;
 943
 944             zunmbr_("Q", "L", "N", m, n, n, cdum, m, cdum, cdum, m, cdum, &
 945                     c_n1, &ierr);
 946             lwork_zunmbr_qln_mn__ = (integer) cdum[0].r;
 947
 948             zunmbr_("Q", "L", "N", n, n, n, cdum, n, cdum, cdum, n, cdum, &
 949                     c_n1, &ierr);
 950             lwork_zunmbr_qln_nn__ = (integer) cdum[0].r;
 951
 952             if (*m >= mnthr1) {
 953                 if (wntqn) {
 954
 955 /*                 Path 1 (M >> N, JOBZ='N') */
 956
 957                     maxwrk = *n + lwork_zgeqrf_mn__;
 958 /* Computing MAX */
 959                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zgebrd_nn__;
 960                     maxwrk = f2cmax(i__1,i__2);
 961                     minwrk = *n * 3;
 962                 } else if (wntqo) {
 963
 964 /*                 Path 2 (M >> N, JOBZ='O') */
 965
 966                     wrkbl = *n + lwork_zgeqrf_mn__;
 967 /* Computing MAX */
 968                     i__1 = wrkbl, i__2 = *n + lwork_zungqr_mn__;
 969                     wrkbl = f2cmax(i__1,i__2);
 970 /* Computing MAX */
 971                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zgebrd_nn__;
 972                     wrkbl = f2cmax(i__1,i__2);
 973 /* Computing MAX */
 974                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_qln_nn__;
 975                     wrkbl = f2cmax(i__1,i__2);
 976 /* Computing MAX */
 977                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
 978                     wrkbl = f2cmax(i__1,i__2);
 979                     maxwrk = *m * *n + *n * *n + wrkbl;
 980                     minwrk = (*n << 1) * *n + *n * 3;
 981                 } else if (wntqs) {
 982
 983 /*                 Path 3 (M >> N, JOBZ='S') */
 984
 985                     wrkbl = *n + lwork_zgeqrf_mn__;
 986 /* Computing MAX */
 987                     i__1 = wrkbl, i__2 = *n + lwork_zungqr_mn__;
 988                     wrkbl = f2cmax(i__1,i__2);
 989 /* Computing MAX */
 990                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zgebrd_nn__;
 991                     wrkbl = f2cmax(i__1,i__2);
 992 /* Computing MAX */
 993                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_qln_nn__;
 994                     wrkbl = f2cmax(i__1,i__2);
 995 /* Computing MAX */
 996                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
 997                     wrkbl = f2cmax(i__1,i__2);
 998                     maxwrk = *n * *n + wrkbl;
 999                     minwrk = *n * *n + *n * 3;
1000                 } else if (wntqa) {
1001
1002 /*                 Path 4 (M >> N, JOBZ='A') */
1003
1004                     wrkbl = *n + lwork_zgeqrf_mn__;
1005 /* Computing MAX */
1006                     i__1 = wrkbl, i__2 = *n + lwork_zungqr_mm__;
1007                     wrkbl = f2cmax(i__1,i__2);
1008 /* Computing MAX */
1009                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zgebrd_nn__;
1010                     wrkbl = f2cmax(i__1,i__2);
1011 /* Computing MAX */
1012                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_qln_nn__;
1013                     wrkbl = f2cmax(i__1,i__2);
1014 /* Computing MAX */
1015                     i__1 = wrkbl, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
1016                     wrkbl = f2cmax(i__1,i__2);
1017                     maxwrk = *n * *n + wrkbl;
1018 /* Computing MAX */
1019                     i__1 = *n * 3, i__2 = *n + *m;
1020                     minwrk = *n * *n + f2cmax(i__1,i__2);
1021                 }
1022             } else if (*m >= mnthr2) {
1023
1024 /*              Path 5 (M >> N, but not as much as MNTHR1) */
1025
1026                 maxwrk = (*n << 1) + lwork_zgebrd_mn__;
1027                 minwrk = (*n << 1) + *m;
1028                 if (wntqo) {
1029 /*                 Path 5o (M >> N, JOBZ='O') */
1030 /* Computing MAX */
1031                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_p_nn__;
1032                     maxwrk = f2cmax(i__1,i__2);
1033 /* Computing MAX */
1034                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_q_mn__;
1035                     maxwrk = f2cmax(i__1,i__2);
1036                     maxwrk += *m * *n;
1037                     minwrk += *n * *n;
1038                 } else if (wntqs) {
1039 /*                 Path 5s (M >> N, JOBZ='S') */
1040 /* Computing MAX */
1041                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_p_nn__;
1042                     maxwrk = f2cmax(i__1,i__2);
1043 /* Computing MAX */
1044                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_q_mn__;
1045                     maxwrk = f2cmax(i__1,i__2);
1046                 } else if (wntqa) {
1047 /*                 Path 5a (M >> N, JOBZ='A') */
1048 /* Computing MAX */
1049                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_p_nn__;
1050                     maxwrk = f2cmax(i__1,i__2);
1051 /* Computing MAX */
1052                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zungbr_q_mm__;
1053                     maxwrk = f2cmax(i__1,i__2);
1054                 }
1055             } else {
1056
1057 /*              Path 6 (M >= N, but not much larger) */
1058
1059                 maxwrk = (*n << 1) + lwork_zgebrd_mn__;
1060                 minwrk = (*n << 1) + *m;
1061                 if (wntqo) {
1062 /*                 Path 6o (M >= N, JOBZ='O') */
1063 /* Computing MAX */
1064                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
1065                     maxwrk = f2cmax(i__1,i__2);
1066 /* Computing MAX */
1067                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_qln_mn__;
1068                     maxwrk = f2cmax(i__1,i__2);
1069                     maxwrk += *m * *n;
1070                     minwrk += *n * *n;
1071                 } else if (wntqs) {
1072 /*                 Path 6s (M >= N, JOBZ='S') */
1073 /* Computing MAX */
1074                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_qln_mn__;
1075                     maxwrk = f2cmax(i__1,i__2);
1076 /* Computing MAX */
1077                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
1078                     maxwrk = f2cmax(i__1,i__2);
1079                 } else if (wntqa) {
1080 /*                 Path 6a (M >= N, JOBZ='A') */
1081 /* Computing MAX */
1082                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_qln_mm__;
1083                     maxwrk = f2cmax(i__1,i__2);
1084 /* Computing MAX */
1085                     i__1 = maxwrk, i__2 = (*n << 1) + lwork_zunmbr_prc_nn__;
1086                     maxwrk = f2cmax(i__1,i__2);
1087                 }
1088             }
1089         } else if (minmn > 0) {
1090
1091 /*           There is no complex work space needed for bidiagonal SVD */
1092 /*           The real work space needed for bidiagonal SVD (dbdsdc) is */
1093 /*           BDSPAC = 3*M*M + 4*M for singular values and vectors; */
1094 /*           BDSPAC = 4*M         for singular values only; */
1095 /*           not including e, RU, and RVT matrices. */
1096
1097 /*           Compute space preferred for each routine */
1098             zgebrd_(m, n, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
1099             lwork_zgebrd_mn__ = (integer) cdum[0].r;
1100
1101             zgebrd_(m, m, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
1102             lwork_zgebrd_mm__ = (integer) cdum[0].r;
1103
1104             zgelqf_(m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
1105             lwork_zgelqf_mn__ = (integer) cdum[0].r;
1106
1107             zungbr_("P", m, n, m, cdum, m, cdum, cdum, &c_n1, &ierr);
1108             lwork_zungbr_p_mn__ = (integer) cdum[0].r;
1109
1110             zungbr_("P", n, n, m, cdum, n, cdum, cdum, &c_n1, &ierr);
1111             lwork_zungbr_p_nn__ = (integer) cdum[0].r;
1112
1113             zungbr_("Q", m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
1114             lwork_zungbr_q_mm__ = (integer) cdum[0].r;
1115
1116             zunglq_(m, n, m, cdum, m, cdum, cdum, &c_n1, &ierr);
1117             lwork_zunglq_mn__ = (integer) cdum[0].r;
1118
1119             zunglq_(n, n, m, cdum, n, cdum, cdum, &c_n1, &ierr);
1120             lwork_zunglq_nn__ = (integer) cdum[0].r;
1121
1122             zunmbr_("P", "R", "C", m, m, m, cdum, m, cdum, cdum, m, cdum, &
1123                     c_n1, &ierr);
1124             lwork_zunmbr_prc_mm__ = (integer) cdum[0].r;
1125
1126             zunmbr_("P", "R", "C", m, n, m, cdum, m, cdum, cdum, m, cdum, &
1127                     c_n1, &ierr);
1128             lwork_zunmbr_prc_mn__ = (integer) cdum[0].r;
1129
1130             zunmbr_("P", "R", "C", n, n, m, cdum, n, cdum, cdum, n, cdum, &
1131                     c_n1, &ierr);
1132             lwork_zunmbr_prc_nn__ = (integer) cdum[0].r;
1133
1134             zunmbr_("Q", "L", "N", m, m, m, cdum, m, cdum, cdum, m, cdum, &
1135                     c_n1, &ierr);
1136             lwork_zunmbr_qln_mm__ = (integer) cdum[0].r;
1137
1138             if (*n >= mnthr1) {
1139                 if (wntqn) {
1140
1141 /*                 Path 1t (N >> M, JOBZ='N') */
1142
1143                     maxwrk = *m + lwork_zgelqf_mn__;
1144 /* Computing MAX */
1145                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zgebrd_mm__;
1146                     maxwrk = f2cmax(i__1,i__2);
1147                     minwrk = *m * 3;
1148                 } else if (wntqo) {
1149
1150 /*                 Path 2t (N >> M, JOBZ='O') */
1151
1152                     wrkbl = *m + lwork_zgelqf_mn__;
1153 /* Computing MAX */
1154                     i__1 = wrkbl, i__2 = *m + lwork_zunglq_mn__;
1155                     wrkbl = f2cmax(i__1,i__2);
1156 /* Computing MAX */
1157                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zgebrd_mm__;
1158                     wrkbl = f2cmax(i__1,i__2);
1159 /* Computing MAX */
1160                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1161                     wrkbl = f2cmax(i__1,i__2);
1162 /* Computing MAX */
1163                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_prc_mm__;
1164                     wrkbl = f2cmax(i__1,i__2);
1165                     maxwrk = *m * *n + *m * *m + wrkbl;
1166                     minwrk = (*m << 1) * *m + *m * 3;
1167                 } else if (wntqs) {
1168
1169 /*                 Path 3t (N >> M, JOBZ='S') */
1170
1171                     wrkbl = *m + lwork_zgelqf_mn__;
1172 /* Computing MAX */
1173                     i__1 = wrkbl, i__2 = *m + lwork_zunglq_mn__;
1174                     wrkbl = f2cmax(i__1,i__2);
1175 /* Computing MAX */
1176                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zgebrd_mm__;
1177                     wrkbl = f2cmax(i__1,i__2);
1178 /* Computing MAX */
1179                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1180                     wrkbl = f2cmax(i__1,i__2);
1181 /* Computing MAX */
1182                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_prc_mm__;
1183                     wrkbl = f2cmax(i__1,i__2);
1184                     maxwrk = *m * *m + wrkbl;
1185                     minwrk = *m * *m + *m * 3;
1186                 } else if (wntqa) {
1187
1188 /*                 Path 4t (N >> M, JOBZ='A') */
1189
1190                     wrkbl = *m + lwork_zgelqf_mn__;
1191 /* Computing MAX */
1192                     i__1 = wrkbl, i__2 = *m + lwork_zunglq_nn__;
1193                     wrkbl = f2cmax(i__1,i__2);
1194 /* Computing MAX */
1195                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zgebrd_mm__;
1196                     wrkbl = f2cmax(i__1,i__2);
1197 /* Computing MAX */
1198                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1199                     wrkbl = f2cmax(i__1,i__2);
1200 /* Computing MAX */
1201                     i__1 = wrkbl, i__2 = (*m << 1) + lwork_zunmbr_prc_mm__;
1202                     wrkbl = f2cmax(i__1,i__2);
1203                     maxwrk = *m * *m + wrkbl;
1204 /* Computing MAX */
1205                     i__1 = *m * 3, i__2 = *m + *n;
1206                     minwrk = *m * *m + f2cmax(i__1,i__2);
1207                 }
1208             } else if (*n >= mnthr2) {
1209
1210 /*              Path 5t (N >> M, but not as much as MNTHR1) */
1211
1212                 maxwrk = (*m << 1) + lwork_zgebrd_mn__;
1213                 minwrk = (*m << 1) + *n;
1214                 if (wntqo) {
1215 /*                 Path 5to (N >> M, JOBZ='O') */
1216 /* Computing MAX */
1217                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_q_mm__;
1218                     maxwrk = f2cmax(i__1,i__2);
1219 /* Computing MAX */
1220                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_p_mn__;
1221                     maxwrk = f2cmax(i__1,i__2);
1222                     maxwrk += *m * *n;
1223                     minwrk += *m * *m;
1224                 } else if (wntqs) {
1225 /*                 Path 5ts (N >> M, JOBZ='S') */
1226 /* Computing MAX */
1227                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_q_mm__;
1228                     maxwrk = f2cmax(i__1,i__2);
1229 /* Computing MAX */
1230                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_p_mn__;
1231                     maxwrk = f2cmax(i__1,i__2);
1232                 } else if (wntqa) {
1233 /*                 Path 5ta (N >> M, JOBZ='A') */
1234 /* Computing MAX */
1235                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_q_mm__;
1236                     maxwrk = f2cmax(i__1,i__2);
1237 /* Computing MAX */
1238                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zungbr_p_nn__;
1239                     maxwrk = f2cmax(i__1,i__2);
1240                 }
1241             } else {
1242
1243 /*              Path 6t (N > M, but not much larger) */
1244
1245                 maxwrk = (*m << 1) + lwork_zgebrd_mn__;
1246                 minwrk = (*m << 1) + *n;
1247                 if (wntqo) {
1248 /*                 Path 6to (N > M, JOBZ='O') */
1249 /* Computing MAX */
1250                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1251                     maxwrk = f2cmax(i__1,i__2);
1252 /* Computing MAX */
1253                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_prc_mn__;
1254                     maxwrk = f2cmax(i__1,i__2);
1255                     maxwrk += *m * *n;
1256                     minwrk += *m * *m;
1257                 } else if (wntqs) {
1258 /*                 Path 6ts (N > M, JOBZ='S') */
1259 /* Computing MAX */
1260                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1261                     maxwrk = f2cmax(i__1,i__2);
1262 /* Computing MAX */
1263                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_prc_mn__;
1264                     maxwrk = f2cmax(i__1,i__2);
1265                 } else if (wntqa) {
1266 /*                 Path 6ta (N > M, JOBZ='A') */
1267 /* Computing MAX */
1268                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_qln_mm__;
1269                     maxwrk = f2cmax(i__1,i__2);
1270 /* Computing MAX */
1271                     i__1 = maxwrk, i__2 = (*m << 1) + lwork_zunmbr_prc_nn__;
1272                     maxwrk = f2cmax(i__1,i__2);
1273                 }
1274             }
1275         }
1276         maxwrk = f2cmax(maxwrk,minwrk);
1277     }
1278     if (*info == 0) {
1279         work[1].r = (doublereal) maxwrk, work[1].i = 0.;
1280         if (*lwork < minwrk && ! lquery) {
1281             *info = -12;
1282         }
1283     }
1284
1285     if (*info != 0) {
1286         i__1 = -(*info);
1287         xerbla_("ZGESDD", &i__1, (ftnlen)6);
1288         return 0;
1289     } else if (lquery) {
1290         return 0;
1291     }
1292
1293 /*     Quick return if possible */
1294
1295     if (*m == 0 || *n == 0) {
1296         return 0;
1297     }
1298
1299 /*     Get machine constants */
1300
1301     eps = dlamch_("P");
1302     smlnum = sqrt(dlamch_("S")) / eps;
1303     bignum = 1. / smlnum;
1304
1305 /*     Scale A if f2cmax element outside range [SMLNUM,BIGNUM] */
1306
1307     anrm = zlange_("M", m, n, &a[a_offset], lda, dum);
1308     if (disnan_(&anrm)) {
1309         *info = -4;
1310         return 0;
1311     }
1312     iscl = 0;
1313     if (anrm > 0. && anrm < smlnum) {
1314         iscl = 1;
1315         zlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
1316                 ierr);
1317     } else if (anrm > bignum) {
1318         iscl = 1;
1319         zlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
1320                 ierr);
1321     }
1322
1323     if (*m >= *n) {
1324
1325 /*        A has at least as many rows as columns. If A has sufficiently */
1326 /*        more rows than columns, first reduce using the QR */
1327 /*        decomposition (if sufficient workspace available) */
1328
1329         if (*m >= mnthr1) {
1330
1331             if (wntqn) {
1332
1333 /*              Path 1 (M >> N, JOBZ='N') */
1334 /*              No singular vectors to be computed */
1335
1336                 itau = 1;
1337                 nwork = itau + *n;
1338
1339 /*              Compute A=Q*R */
1340 /*              CWorkspace: need   N [tau] + N    [work] */
1341 /*              CWorkspace: prefer N [tau] + N*NB [work] */
1342 /*              RWorkspace: need   0 */
1343
1344                 i__1 = *lwork - nwork + 1;
1345                 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1346                         i__1, &ierr);
1347
1348 /*              Zero out below R */
1349
1350                 i__1 = *n - 1;
1351                 i__2 = *n - 1;
1352                 zlaset_("L", &i__1, &i__2, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1353                 ie = 1;
1354                 itauq = 1;
1355                 itaup = itauq + *n;
1356                 nwork = itaup + *n;
1357
1358 /*              Bidiagonalize R in A */
1359 /*              CWorkspace: need   2*N [tauq, taup] + N      [work] */
1360 /*              CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work] */
1361 /*              RWorkspace: need   N [e] */
1362
1363                 i__1 = *lwork - nwork + 1;
1364                 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1365                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1366                 nrwork = ie + *n;
1367
1368 /*              Perform bidiagonal SVD, compute singular values only */
1369 /*              CWorkspace: need   0 */
1370 /*              RWorkspace: need   N [e] + BDSPAC */
1371
1372                 dbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1373                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1374
1375             } else if (wntqo) {
1376
1377 /*              Path 2 (M >> N, JOBZ='O') */
1378 /*              N left singular vectors to be overwritten on A and */
1379 /*              N right singular vectors to be computed in VT */
1380
1381                 iu = 1;
1382
1383 /*              WORK(IU) is N by N */
1384
1385                 ldwrku = *n;
1386                 ir = iu + ldwrku * *n;
1387                 if (*lwork >= *m * *n + *n * *n + *n * 3) {
1388
1389 /*                 WORK(IR) is M by N */
1390
1391                     ldwrkr = *m;
1392                 } else {
1393                     ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
1394                 }
1395                 itau = ir + ldwrkr * *n;
1396                 nwork = itau + *n;
1397
1398 /*              Compute A=Q*R */
1399 /*              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work] */
1400 /*              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] */
1401 /*              RWorkspace: need   0 */
1402
1403                 i__1 = *lwork - nwork + 1;
1404                 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1405                         i__1, &ierr);
1406
1407 /*              Copy R to WORK( IR ), zeroing out below it */
1408
1409                 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1410                 i__1 = *n - 1;
1411                 i__2 = *n - 1;
1412                 zlaset_("L", &i__1, &i__2, &c_b1, &c_b1, &work[ir + 1], &
1413                         ldwrkr);
1414
1415 /*              Generate Q in A */
1416 /*              CWorkspace: need   N*N [U] + N*N [R] + N [tau] + N    [work] */
1417 /*              CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] */
1418 /*              RWorkspace: need   0 */
1419
1420                 i__1 = *lwork - nwork + 1;
1421                 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
1422                          &i__1, &ierr);
1423                 ie = 1;
1424                 itauq = itau;
1425                 itaup = itauq + *n;
1426                 nwork = itaup + *n;
1427
1428 /*              Bidiagonalize R in WORK(IR) */
1429 /*              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N      [work] */
1430 /*              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] */
1431 /*              RWorkspace: need   N [e] */
1432
1433                 i__1 = *lwork - nwork + 1;
1434                 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
1435                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1436
1437 /*              Perform bidiagonal SVD, computing left singular vectors */
1438 /*              of R in WORK(IRU) and computing right singular vectors */
1439 /*              of R in WORK(IRVT) */
1440 /*              CWorkspace: need   0 */
1441 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1442
1443                 iru = ie + *n;
1444                 irvt = iru + *n * *n;
1445                 nrwork = irvt + *n * *n;
1446                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1447                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1448                         info);
1449
1450 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1451 /*              Overwrite WORK(IU) by the left singular vectors of R */
1452 /*              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work] */
1453 /*              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1454 /*              RWorkspace: need   0 */
1455
1456                 zlacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
1457                 i__1 = *lwork - nwork + 1;
1458                 zunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
1459                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
1460                         ierr);
1461
1462 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
1463 /*              Overwrite VT by the right singular vectors of R */
1464 /*              CWorkspace: need   N*N [U] + N*N [R] + 2*N [tauq, taup] + N    [work] */
1465 /*              CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1466 /*              RWorkspace: need   0 */
1467
1468                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1469                 i__1 = *lwork - nwork + 1;
1470                 zunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
1471                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1472                         ierr);
1473
1474 /*              Multiply Q in A by left singular vectors of R in */
1475 /*              WORK(IU), storing result in WORK(IR) and copying to A */
1476 /*              CWorkspace: need   N*N [U] + N*N [R] */
1477 /*              CWorkspace: prefer N*N [U] + M*N [R] */
1478 /*              RWorkspace: need   0 */
1479
1480                 i__1 = *m;
1481                 i__2 = ldwrkr;
1482                 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
1483                         i__2) {
1484 /* Computing MIN */
1485                     i__3 = *m - i__ + 1;
1486                     chunk = f2cmin(i__3,ldwrkr);
1487                     zgemm_("N", "N", &chunk, n, n, &c_b2, &a[i__ + a_dim1],
1488                             lda, &work[iu], &ldwrku, &c_b1, &work[ir], &
1489                             ldwrkr);
1490                     zlacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a[i__ +
1491                             a_dim1], lda);
1492 /* L10: */
1493                 }
1494
1495             } else if (wntqs) {
1496
1497 /*              Path 3 (M >> N, JOBZ='S') */
1498 /*              N left singular vectors to be computed in U and */
1499 /*              N right singular vectors to be computed in VT */
1500
1501                 ir = 1;
1502
1503 /*              WORK(IR) is N by N */
1504
1505                 ldwrkr = *n;
1506                 itau = ir + ldwrkr * *n;
1507                 nwork = itau + *n;
1508
1509 /*              Compute A=Q*R */
1510 /*              CWorkspace: need   N*N [R] + N [tau] + N    [work] */
1511 /*              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] */
1512 /*              RWorkspace: need   0 */
1513
1514                 i__2 = *lwork - nwork + 1;
1515                 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1516                         i__2, &ierr);
1517
1518 /*              Copy R to WORK(IR), zeroing out below it */
1519
1520                 zlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1521                 i__2 = *n - 1;
1522                 i__1 = *n - 1;
1523                 zlaset_("L", &i__2, &i__1, &c_b1, &c_b1, &work[ir + 1], &
1524                         ldwrkr);
1525
1526 /*              Generate Q in A */
1527 /*              CWorkspace: need   N*N [R] + N [tau] + N    [work] */
1528 /*              CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] */
1529 /*              RWorkspace: need   0 */
1530
1531                 i__2 = *lwork - nwork + 1;
1532                 zungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
1533                          &i__2, &ierr);
1534                 ie = 1;
1535                 itauq = itau;
1536                 itaup = itauq + *n;
1537                 nwork = itaup + *n;
1538
1539 /*              Bidiagonalize R in WORK(IR) */
1540 /*              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N      [work] */
1541 /*              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] */
1542 /*              RWorkspace: need   N [e] */
1543
1544                 i__2 = *lwork - nwork + 1;
1545                 zgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
1546                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1547
1548 /*              Perform bidiagonal SVD, computing left singular vectors */
1549 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1550 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1551 /*              CWorkspace: need   0 */
1552 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1553
1554                 iru = ie + *n;
1555                 irvt = iru + *n * *n;
1556                 nrwork = irvt + *n * *n;
1557                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1558                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1559                         info);
1560
1561 /*              Copy real matrix RWORK(IRU) to complex matrix U */
1562 /*              Overwrite U by left singular vectors of R */
1563 /*              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work] */
1564 /*              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1565 /*              RWorkspace: need   0 */
1566
1567                 zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
1568                 i__2 = *lwork - nwork + 1;
1569                 zunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
1570                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
1571
1572 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
1573 /*              Overwrite VT by right singular vectors of R */
1574 /*              CWorkspace: need   N*N [R] + 2*N [tauq, taup] + N    [work] */
1575 /*              CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1576 /*              RWorkspace: need   0 */
1577
1578                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1579                 i__2 = *lwork - nwork + 1;
1580                 zunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
1581                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1582                         ierr);
1583
1584 /*              Multiply Q in A by left singular vectors of R in */
1585 /*              WORK(IR), storing result in U */
1586 /*              CWorkspace: need   N*N [R] */
1587 /*              RWorkspace: need   0 */
1588
1589                 zlacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
1590                 zgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &work[ir],
1591                          &ldwrkr, &c_b1, &u[u_offset], ldu);
1592
1593             } else if (wntqa) {
1594
1595 /*              Path 4 (M >> N, JOBZ='A') */
1596 /*              M left singular vectors to be computed in U and */
1597 /*              N right singular vectors to be computed in VT */
1598
1599                 iu = 1;
1600
1601 /*              WORK(IU) is N by N */
1602
1603                 ldwrku = *n;
1604                 itau = iu + ldwrku * *n;
1605                 nwork = itau + *n;
1606
1607 /*              Compute A=Q*R, copying result to U */
1608 /*              CWorkspace: need   N*N [U] + N [tau] + N    [work] */
1609 /*              CWorkspace: prefer N*N [U] + N [tau] + N*NB [work] */
1610 /*              RWorkspace: need   0 */
1611
1612                 i__2 = *lwork - nwork + 1;
1613                 zgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1614                         i__2, &ierr);
1615                 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1616
1617 /*              Generate Q in U */
1618 /*              CWorkspace: need   N*N [U] + N [tau] + M    [work] */
1619 /*              CWorkspace: prefer N*N [U] + N [tau] + M*NB [work] */
1620 /*              RWorkspace: need   0 */
1621
1622                 i__2 = *lwork - nwork + 1;
1623                 zungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
1624                          &i__2, &ierr);
1625
1626 /*              Produce R in A, zeroing out below it */
1627
1628                 i__2 = *n - 1;
1629                 i__1 = *n - 1;
1630                 zlaset_("L", &i__2, &i__1, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1631                 ie = 1;
1632                 itauq = itau;
1633                 itaup = itauq + *n;
1634                 nwork = itaup + *n;
1635
1636 /*              Bidiagonalize R in A */
1637 /*              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N      [work] */
1638 /*              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work] */
1639 /*              RWorkspace: need   N [e] */
1640
1641                 i__2 = *lwork - nwork + 1;
1642                 zgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1643                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1644                 iru = ie + *n;
1645                 irvt = iru + *n * *n;
1646                 nrwork = irvt + *n * *n;
1647
1648 /*              Perform bidiagonal SVD, computing left singular vectors */
1649 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1650 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1651 /*              CWorkspace: need   0 */
1652 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1653
1654                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1655                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1656                         info);
1657
1658 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1659 /*              Overwrite WORK(IU) by left singular vectors of R */
1660 /*              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work] */
1661 /*              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] */
1662 /*              RWorkspace: need   0 */
1663
1664                 zlacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
1665                 i__2 = *lwork - nwork + 1;
1666                 zunmbr_("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
1667                         itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
1668                         ierr);
1669
1670 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
1671 /*              Overwrite VT by right singular vectors of R */
1672 /*              CWorkspace: need   N*N [U] + 2*N [tauq, taup] + N    [work] */
1673 /*              CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] */
1674 /*              RWorkspace: need   0 */
1675
1676                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1677                 i__2 = *lwork - nwork + 1;
1678                 zunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
1679                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1680                         ierr);
1681
1682 /*              Multiply Q in U by left singular vectors of R in */
1683 /*              WORK(IU), storing result in A */
1684 /*              CWorkspace: need   N*N [U] */
1685 /*              RWorkspace: need   0 */
1686
1687                 zgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &work[iu],
1688                          &ldwrku, &c_b1, &a[a_offset], lda);
1689
1690 /*              Copy left singular vectors of A from A to U */
1691
1692                 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1693
1694             }
1695
1696         } else if (*m >= mnthr2) {
1697
1698 /*           MNTHR2 <= M < MNTHR1 */
1699
1700 /*           Path 5 (M >> N, but not as much as MNTHR1) */
1701 /*           Reduce to bidiagonal form without QR decomposition, use */
1702 /*           ZUNGBR and matrix multiplication to compute singular vectors */
1703
1704             ie = 1;
1705             nrwork = ie + *n;
1706             itauq = 1;
1707             itaup = itauq + *n;
1708             nwork = itaup + *n;
1709
1710 /*           Bidiagonalize A */
1711 /*           CWorkspace: need   2*N [tauq, taup] + M        [work] */
1712 /*           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] */
1713 /*           RWorkspace: need   N [e] */
1714
1715             i__2 = *lwork - nwork + 1;
1716             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
1717                     &work[itaup], &work[nwork], &i__2, &ierr);
1718             if (wntqn) {
1719
1720 /*              Path 5n (M >> N, JOBZ='N') */
1721 /*              Compute singular values only */
1722 /*              CWorkspace: need   0 */
1723 /*              RWorkspace: need   N [e] + BDSPAC */
1724
1725                 dbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1726                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1727             } else if (wntqo) {
1728                 iu = nwork;
1729                 iru = nrwork;
1730                 irvt = iru + *n * *n;
1731                 nrwork = irvt + *n * *n;
1732
1733 /*              Path 5o (M >> N, JOBZ='O') */
1734 /*              Copy A to VT, generate P**H */
1735 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
1736 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1737 /*              RWorkspace: need   0 */
1738
1739                 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1740                 i__2 = *lwork - nwork + 1;
1741                 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1742                         work[nwork], &i__2, &ierr);
1743
1744 /*              Generate Q in A */
1745 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
1746 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1747 /*              RWorkspace: need   0 */
1748
1749                 i__2 = *lwork - nwork + 1;
1750                 zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
1751                         nwork], &i__2, &ierr);
1752
1753                 if (*lwork >= *m * *n + *n * 3) {
1754
1755 /*                 WORK( IU ) is M by N */
1756
1757                     ldwrku = *m;
1758                 } else {
1759
1760 /*                 WORK(IU) is LDWRKU by N */
1761
1762                     ldwrku = (*lwork - *n * 3) / *n;
1763                 }
1764                 nwork = iu + ldwrku * *n;
1765
1766 /*              Perform bidiagonal SVD, computing left singular vectors */
1767 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1768 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1769 /*              CWorkspace: need   0 */
1770 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1771
1772                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1773                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1774                         info);
1775
1776 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
1777 /*              storing the result in WORK(IU), copying to VT */
1778 /*              CWorkspace: need   2*N [tauq, taup] + N*N [U] */
1779 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1780
1781                 zlarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &work[iu]
1782                         , &ldwrku, &rwork[nrwork]);
1783                 zlacpy_("F", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);
1784
1785 /*              Multiply Q in A by real matrix RWORK(IRU), storing the */
1786 /*              result in WORK(IU), copying to A */
1787 /*              CWorkspace: need   2*N [tauq, taup] + N*N [U] */
1788 /*              CWorkspace: prefer 2*N [tauq, taup] + M*N [U] */
1789 /*              RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork] */
1790 /*              RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1791
1792                 nrwork = irvt;
1793                 i__2 = *m;
1794                 i__1 = ldwrku;
1795                 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1796                         i__1) {
1797 /* Computing MIN */
1798                     i__3 = *m - i__ + 1;
1799                     chunk = f2cmin(i__3,ldwrku);
1800                     zlacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru], n,
1801                             &work[iu], &ldwrku, &rwork[nrwork]);
1802                     zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
1803                             a_dim1], lda);
1804 /* L20: */
1805                 }
1806
1807             } else if (wntqs) {
1808
1809 /*              Path 5s (M >> N, JOBZ='S') */
1810 /*              Copy A to VT, generate P**H */
1811 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
1812 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1813 /*              RWorkspace: need   0 */
1814
1815                 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1816                 i__1 = *lwork - nwork + 1;
1817                 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1818                         work[nwork], &i__1, &ierr);
1819
1820 /*              Copy A to U, generate Q */
1821 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
1822 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1823 /*              RWorkspace: need   0 */
1824
1825                 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1826                 i__1 = *lwork - nwork + 1;
1827                 zungbr_("Q", m, n, n, &u[u_offset], ldu, &work[itauq], &work[
1828                         nwork], &i__1, &ierr);
1829
1830 /*              Perform bidiagonal SVD, computing left singular vectors */
1831 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1832 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1833 /*              CWorkspace: need   0 */
1834 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1835
1836                 iru = nrwork;
1837                 irvt = iru + *n * *n;
1838                 nrwork = irvt + *n * *n;
1839                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1840                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1841                         info);
1842
1843 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
1844 /*              storing the result in A, copying to VT */
1845 /*              CWorkspace: need   0 */
1846 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1847
1848                 zlarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
1849                         a_offset], lda, &rwork[nrwork]);
1850                 zlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1851
1852 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
1853 /*              result in A, copying to U */
1854 /*              CWorkspace: need   0 */
1855 /*              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1856
1857                 nrwork = irvt;
1858                 zlacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset],
1859                          lda, &rwork[nrwork]);
1860                 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1861             } else {
1862
1863 /*              Path 5a (M >> N, JOBZ='A') */
1864 /*              Copy A to VT, generate P**H */
1865 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
1866 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1867 /*              RWorkspace: need   0 */
1868
1869                 zlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1870                 i__1 = *lwork - nwork + 1;
1871                 zungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1872                         work[nwork], &i__1, &ierr);
1873
1874 /*              Copy A to U, generate Q */
1875 /*              CWorkspace: need   2*N [tauq, taup] + M    [work] */
1876 /*              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] */
1877 /*              RWorkspace: need   0 */
1878
1879                 zlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1880                 i__1 = *lwork - nwork + 1;
1881                 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
1882                         nwork], &i__1, &ierr);
1883
1884 /*              Perform bidiagonal SVD, computing left singular vectors */
1885 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1886 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1887 /*              CWorkspace: need   0 */
1888 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1889
1890                 iru = nrwork;
1891                 irvt = iru + *n * *n;
1892                 nrwork = irvt + *n * *n;
1893                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1894                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1895                         info);
1896
1897 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
1898 /*              storing the result in A, copying to VT */
1899 /*              CWorkspace: need   0 */
1900 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1901
1902                 zlarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
1903                         a_offset], lda, &rwork[nrwork]);
1904                 zlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1905
1906 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
1907 /*              result in A, copying to U */
1908 /*              CWorkspace: need   0 */
1909 /*              RWorkspace: need   N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1910
1911                 nrwork = irvt;
1912                 zlacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset],
1913                          lda, &rwork[nrwork]);
1914                 zlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1915             }
1916
1917         } else {
1918
1919 /*           M .LT. MNTHR2 */
1920
1921 /*           Path 6 (M >= N, but not much larger) */
1922 /*           Reduce to bidiagonal form without QR decomposition */
1923 /*           Use ZUNMBR to compute singular vectors */
1924
1925             ie = 1;
1926             nrwork = ie + *n;
1927             itauq = 1;
1928             itaup = itauq + *n;
1929             nwork = itaup + *n;
1930
1931 /*           Bidiagonalize A */
1932 /*           CWorkspace: need   2*N [tauq, taup] + M        [work] */
1933 /*           CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] */
1934 /*           RWorkspace: need   N [e] */
1935
1936             i__1 = *lwork - nwork + 1;
1937             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
1938                     &work[itaup], &work[nwork], &i__1, &ierr);
1939             if (wntqn) {
1940
1941 /*              Path 6n (M >= N, JOBZ='N') */
1942 /*              Compute singular values only */
1943 /*              CWorkspace: need   0 */
1944 /*              RWorkspace: need   N [e] + BDSPAC */
1945
1946                 dbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1947                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1948             } else if (wntqo) {
1949                 iu = nwork;
1950                 iru = nrwork;
1951                 irvt = iru + *n * *n;
1952                 nrwork = irvt + *n * *n;
1953                 if (*lwork >= *m * *n + *n * 3) {
1954
1955 /*                 WORK( IU ) is M by N */
1956
1957                     ldwrku = *m;
1958                 } else {
1959
1960 /*                 WORK( IU ) is LDWRKU by N */
1961
1962                     ldwrku = (*lwork - *n * 3) / *n;
1963                 }
1964                 nwork = iu + ldwrku * *n;
1965
1966 /*              Path 6o (M >= N, JOBZ='O') */
1967 /*              Perform bidiagonal SVD, computing left singular vectors */
1968 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
1969 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
1970 /*              CWorkspace: need   0 */
1971 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1972
1973                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1974                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1975                         info);
1976
1977 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
1978 /*              Overwrite VT by right singular vectors of A */
1979 /*              CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work] */
1980 /*              CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] */
1981 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] */
1982
1983                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1984                 i__1 = *lwork - nwork + 1;
1985                 zunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
1986                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1987                         ierr);
1988
1989                 if (*lwork >= *m * *n + *n * 3) {
1990
1991 /*                 Path 6o-fast */
1992 /*                 Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1993 /*                 Overwrite WORK(IU) by left singular vectors of A, copying */
1994 /*                 to A */
1995 /*                 CWorkspace: need   2*N [tauq, taup] + M*N [U] + N    [work] */
1996 /*                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work] */
1997 /*                 RWorkspace: need   N [e] + N*N [RU] */
1998
1999                     zlaset_("F", m, n, &c_b1, &c_b1, &work[iu], &ldwrku);
2000                     zlacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
2001                     i__1 = *lwork - nwork + 1;
2002                     zunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
2003                             itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
2004                             ierr);
2005                     zlacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda);
2006                 } else {
2007
2008 /*                 Path 6o-slow */
2009 /*                 Generate Q in A */
2010 /*                 CWorkspace: need   2*N [tauq, taup] + N*N [U] + N    [work] */
2011 /*                 CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] */
2012 /*                 RWorkspace: need   0 */
2013
2014                     i__1 = *lwork - nwork + 1;
2015                     zungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
2016                             work[nwork], &i__1, &ierr);
2017
2018 /*                 Multiply Q in A by real matrix RWORK(IRU), storing the */
2019 /*                 result in WORK(IU), copying to A */
2020 /*                 CWorkspace: need   2*N [tauq, taup] + N*N [U] */
2021 /*                 CWorkspace: prefer 2*N [tauq, taup] + M*N [U] */
2022 /*                 RWorkspace: need   N [e] + N*N [RU] + 2*N*N [rwork] */
2023 /*                 RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
2024
2025                     nrwork = irvt;
2026                     i__1 = *m;
2027                     i__2 = ldwrku;
2028                     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
2029                              i__2) {
2030 /* Computing MIN */
2031                         i__3 = *m - i__ + 1;
2032                         chunk = f2cmin(i__3,ldwrku);
2033                         zlacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru],
2034                                  n, &work[iu], &ldwrku, &rwork[nrwork]);
2035                         zlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
2036                                 a_dim1], lda);
2037 /* L30: */
2038                     }
2039                 }
2040
2041             } else if (wntqs) {
2042
2043 /*              Path 6s (M >= N, JOBZ='S') */
2044 /*              Perform bidiagonal SVD, computing left singular vectors */
2045 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2046 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2047 /*              CWorkspace: need   0 */
2048 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
2049
2050                 iru = nrwork;
2051                 irvt = iru + *n * *n;
2052                 nrwork = irvt + *n * *n;
2053                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
2054                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
2055                         info);
2056
2057 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2058 /*              Overwrite U by left singular vectors of A */
2059 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
2060 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2061 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] */
2062
2063                 zlaset_("F", m, n, &c_b1, &c_b1, &u[u_offset], ldu)
2064                         ;
2065                 zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
2066                 i__2 = *lwork - nwork + 1;
2067                 zunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
2068                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2069
2070 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
2071 /*              Overwrite VT by right singular vectors of A */
2072 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
2073 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2074 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] */
2075
2076                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
2077                 i__2 = *lwork - nwork + 1;
2078                 zunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
2079                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
2080                         ierr);
2081             } else {
2082
2083 /*              Path 6a (M >= N, JOBZ='A') */
2084 /*              Perform bidiagonal SVD, computing left singular vectors */
2085 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2086 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2087 /*              CWorkspace: need   0 */
2088 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
2089
2090                 iru = nrwork;
2091                 irvt = iru + *n * *n;
2092                 nrwork = irvt + *n * *n;
2093                 dbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
2094                         rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
2095                         info);
2096
2097 /*              Set the right corner of U to identity matrix */
2098
2099                 zlaset_("F", m, m, &c_b1, &c_b1, &u[u_offset], ldu)
2100                         ;
2101                 if (*m > *n) {
2102                     i__2 = *m - *n;
2103                     i__1 = *m - *n;
2104                     zlaset_("F", &i__2, &i__1, &c_b1, &c_b2, &u[*n + 1 + (*n
2105                             + 1) * u_dim1], ldu);
2106                 }
2107
2108 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2109 /*              Overwrite U by left singular vectors of A */
2110 /*              CWorkspace: need   2*N [tauq, taup] + M    [work] */
2111 /*              CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] */
2112 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] */
2113
2114                 zlacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
2115                 i__2 = *lwork - nwork + 1;
2116                 zunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2117                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2118
2119 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
2120 /*              Overwrite VT by right singular vectors of A */
2121 /*              CWorkspace: need   2*N [tauq, taup] + N    [work] */
2122 /*              CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2123 /*              RWorkspace: need   N [e] + N*N [RU] + N*N [RVT] */
2124
2125                 zlacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
2126                 i__2 = *lwork - nwork + 1;
2127                 zunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
2128                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
2129                         ierr);
2130             }
2131
2132         }
2133
2134     } else {
2135
2136 /*        A has more columns than rows. If A has sufficiently more */
2137 /*        columns than rows, first reduce using the LQ decomposition (if */
2138 /*        sufficient workspace available) */
2139
2140         if (*n >= mnthr1) {
2141
2142             if (wntqn) {
2143
2144 /*              Path 1t (N >> M, JOBZ='N') */
2145 /*              No singular vectors to be computed */
2146
2147                 itau = 1;
2148                 nwork = itau + *m;
2149
2150 /*              Compute A=L*Q */
2151 /*              CWorkspace: need   M [tau] + M    [work] */
2152 /*              CWorkspace: prefer M [tau] + M*NB [work] */
2153 /*              RWorkspace: need   0 */
2154
2155                 i__2 = *lwork - nwork + 1;
2156                 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2157                         i__2, &ierr);
2158
2159 /*              Zero out above L */
2160
2161                 i__2 = *m - 1;
2162                 i__1 = *m - 1;
2163                 zlaset_("U", &i__2, &i__1, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
2164                         , lda);
2165                 ie = 1;
2166                 itauq = 1;
2167                 itaup = itauq + *m;
2168                 nwork = itaup + *m;
2169
2170 /*              Bidiagonalize L in A */
2171 /*              CWorkspace: need   2*M [tauq, taup] + M      [work] */
2172 /*              CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work] */
2173 /*              RWorkspace: need   M [e] */
2174
2175                 i__2 = *lwork - nwork + 1;
2176                 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2177                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
2178                 nrwork = ie + *m;
2179
2180 /*              Perform bidiagonal SVD, compute singular values only */
2181 /*              CWorkspace: need   0 */
2182 /*              RWorkspace: need   M [e] + BDSPAC */
2183
2184                 dbdsdc_("U", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2185                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2186
2187             } else if (wntqo) {
2188
2189 /*              Path 2t (N >> M, JOBZ='O') */
2190 /*              M right singular vectors to be overwritten on A and */
2191 /*              M left singular vectors to be computed in U */
2192
2193                 ivt = 1;
2194                 ldwkvt = *m;
2195
2196 /*              WORK(IVT) is M by M */
2197
2198                 il = ivt + ldwkvt * *m;
2199                 if (*lwork >= *m * *n + *m * *m + *m * 3) {
2200
2201 /*                 WORK(IL) M by N */
2202
2203                     ldwrkl = *m;
2204                     chunk = *n;
2205                 } else {
2206
2207 /*                 WORK(IL) is M by CHUNK */
2208
2209                     ldwrkl = *m;
2210                     chunk = (*lwork - *m * *m - *m * 3) / *m;
2211                 }
2212                 itau = il + ldwrkl * chunk;
2213                 nwork = itau + *m;
2214
2215 /*              Compute A=L*Q */
2216 /*              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work] */
2217 /*              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] */
2218 /*              RWorkspace: need   0 */
2219
2220                 i__2 = *lwork - nwork + 1;
2221                 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2222                         i__2, &ierr);
2223
2224 /*              Copy L to WORK(IL), zeroing about above it */
2225
2226                 zlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
2227                 i__2 = *m - 1;
2228                 i__1 = *m - 1;
2229                 zlaset_("U", &i__2, &i__1, &c_b1, &c_b1, &work[il + ldwrkl], &
2230                         ldwrkl);
2231
2232 /*              Generate Q in A */
2233 /*              CWorkspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work] */
2234 /*              CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] */
2235 /*              RWorkspace: need   0 */
2236
2237                 i__2 = *lwork - nwork + 1;
2238                 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
2239                          &i__2, &ierr);
2240                 ie = 1;
2241                 itauq = itau;
2242                 itaup = itauq + *m;
2243                 nwork = itaup + *m;
2244
2245 /*              Bidiagonalize L in WORK(IL) */
2246 /*              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M      [work] */
2247 /*              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] */
2248 /*              RWorkspace: need   M [e] */
2249
2250                 i__2 = *lwork - nwork + 1;
2251                 zgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
2252                         itauq], &work[itaup], &work[nwork], &i__2, &ierr);
2253
2254 /*              Perform bidiagonal SVD, computing left singular vectors */
2255 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2256 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2257 /*              CWorkspace: need   0 */
2258 /*              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2259
2260                 iru = ie + *m;
2261                 irvt = iru + *m * *m;
2262                 nrwork = irvt + *m * *m;
2263                 dbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2264                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2265                         info);
2266
2267 /*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
2268 /*              Overwrite WORK(IU) by the left singular vectors of L */
2269 /*              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work] */
2270 /*              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2271 /*              RWorkspace: need   0 */
2272
2273                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2274                 i__2 = *lwork - nwork + 1;
2275                 zunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
2276                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2277
2278 /*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2279 /*              Overwrite WORK(IVT) by the right singular vectors of L */
2280 /*              CWorkspace: need   M*M [VT] + M*M [L] + 2*M [tauq, taup] + M    [work] */
2281 /*              CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2282 /*              RWorkspace: need   0 */
2283
2284                 zlacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2285                 i__2 = *lwork - nwork + 1;
2286                 zunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
2287                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
2288                         ierr);
2289
2290 /*              Multiply right singular vectors of L in WORK(IL) by Q */
2291 /*              in A, storing result in WORK(IL) and copying to A */
2292 /*              CWorkspace: need   M*M [VT] + M*M [L] */
2293 /*              CWorkspace: prefer M*M [VT] + M*N [L] */
2294 /*              RWorkspace: need   0 */
2295
2296                 i__2 = *n;
2297                 i__1 = chunk;
2298                 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2299                         i__1) {
2300 /* Computing MIN */
2301                     i__3 = *n - i__ + 1;
2302                     blk = f2cmin(i__3,chunk);
2303                     zgemm_("N", "N", m, &blk, m, &c_b2, &work[ivt], m, &a[i__
2304                             * a_dim1 + 1], lda, &c_b1, &work[il], &ldwrkl);
2305                     zlacpy_("F", m, &blk, &work[il], &ldwrkl, &a[i__ * a_dim1
2306                             + 1], lda);
2307 /* L40: */
2308                 }
2309
2310             } else if (wntqs) {
2311
2312 /*              Path 3t (N >> M, JOBZ='S') */
2313 /*              M right singular vectors to be computed in VT and */
2314 /*              M left singular vectors to be computed in U */
2315
2316                 il = 1;
2317
2318 /*              WORK(IL) is M by M */
2319
2320                 ldwrkl = *m;
2321                 itau = il + ldwrkl * *m;
2322                 nwork = itau + *m;
2323
2324 /*              Compute A=L*Q */
2325 /*              CWorkspace: need   M*M [L] + M [tau] + M    [work] */
2326 /*              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] */
2327 /*              RWorkspace: need   0 */
2328
2329                 i__1 = *lwork - nwork + 1;
2330                 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2331                         i__1, &ierr);
2332
2333 /*              Copy L to WORK(IL), zeroing out above it */
2334
2335                 zlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
2336                 i__1 = *m - 1;
2337                 i__2 = *m - 1;
2338                 zlaset_("U", &i__1, &i__2, &c_b1, &c_b1, &work[il + ldwrkl], &
2339                         ldwrkl);
2340
2341 /*              Generate Q in A */
2342 /*              CWorkspace: need   M*M [L] + M [tau] + M    [work] */
2343 /*              CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] */
2344 /*              RWorkspace: need   0 */
2345
2346                 i__1 = *lwork - nwork + 1;
2347                 zunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
2348                          &i__1, &ierr);
2349                 ie = 1;
2350                 itauq = itau;
2351                 itaup = itauq + *m;
2352                 nwork = itaup + *m;
2353
2354 /*              Bidiagonalize L in WORK(IL) */
2355 /*              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M      [work] */
2356 /*              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] */
2357 /*              RWorkspace: need   M [e] */
2358
2359                 i__1 = *lwork - nwork + 1;
2360                 zgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
2361                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
2362
2363 /*              Perform bidiagonal SVD, computing left singular vectors */
2364 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2365 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2366 /*              CWorkspace: need   0 */
2367 /*              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2368
2369                 iru = ie + *m;
2370                 irvt = iru + *m * *m;
2371                 nrwork = irvt + *m * *m;
2372                 dbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2373                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2374                         info);
2375
2376 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2377 /*              Overwrite U by left singular vectors of L */
2378 /*              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work] */
2379 /*              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2380 /*              RWorkspace: need   0 */
2381
2382                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2383                 i__1 = *lwork - nwork + 1;
2384                 zunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
2385                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2386
2387 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
2388 /*              Overwrite VT by left singular vectors of L */
2389 /*              CWorkspace: need   M*M [L] + 2*M [tauq, taup] + M    [work] */
2390 /*              CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2391 /*              RWorkspace: need   0 */
2392
2393                 zlacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2394                 i__1 = *lwork - nwork + 1;
2395                 zunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
2396                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2397                         ierr);
2398
2399 /*              Copy VT to WORK(IL), multiply right singular vectors of L */
2400 /*              in WORK(IL) by Q in A, storing result in VT */
2401 /*              CWorkspace: need   M*M [L] */
2402 /*              RWorkspace: need   0 */
2403
2404                 zlacpy_("F", m, m, &vt[vt_offset], ldvt, &work[il], &ldwrkl);
2405                 zgemm_("N", "N", m, n, m, &c_b2, &work[il], &ldwrkl, &a[
2406                         a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
2407
2408             } else if (wntqa) {
2409
2410 /*              Path 4t (N >> M, JOBZ='A') */
2411 /*              N right singular vectors to be computed in VT and */
2412 /*              M left singular vectors to be computed in U */
2413
2414                 ivt = 1;
2415
2416 /*              WORK(IVT) is M by M */
2417
2418                 ldwkvt = *m;
2419                 itau = ivt + ldwkvt * *m;
2420                 nwork = itau + *m;
2421
2422 /*              Compute A=L*Q, copying result to VT */
2423 /*              CWorkspace: need   M*M [VT] + M [tau] + M    [work] */
2424 /*              CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work] */
2425 /*              RWorkspace: need   0 */
2426
2427                 i__1 = *lwork - nwork + 1;
2428                 zgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2429                         i__1, &ierr);
2430                 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2431
2432 /*              Generate Q in VT */
2433 /*              CWorkspace: need   M*M [VT] + M [tau] + N    [work] */
2434 /*              CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work] */
2435 /*              RWorkspace: need   0 */
2436
2437                 i__1 = *lwork - nwork + 1;
2438                 zunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
2439                         nwork], &i__1, &ierr);
2440
2441 /*              Produce L in A, zeroing out above it */
2442
2443                 i__1 = *m - 1;
2444                 i__2 = *m - 1;
2445                 zlaset_("U", &i__1, &i__2, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
2446                         , lda);
2447                 ie = 1;
2448                 itauq = itau;
2449                 itaup = itauq + *m;
2450                 nwork = itaup + *m;
2451
2452 /*              Bidiagonalize L in A */
2453 /*              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M      [work] */
2454 /*              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work] */
2455 /*              RWorkspace: need   M [e] */
2456
2457                 i__1 = *lwork - nwork + 1;
2458                 zgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2459                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
2460
2461 /*              Perform bidiagonal SVD, computing left singular vectors */
2462 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2463 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2464 /*              CWorkspace: need   0 */
2465 /*              RWorkspace: need   M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2466
2467                 iru = ie + *m;
2468                 irvt = iru + *m * *m;
2469                 nrwork = irvt + *m * *m;
2470                 dbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2471                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2472                         info);
2473
2474 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2475 /*              Overwrite U by left singular vectors of L */
2476 /*              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work] */
2477 /*              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] */
2478 /*              RWorkspace: need   0 */
2479
2480                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2481                 i__1 = *lwork - nwork + 1;
2482                 zunmbr_("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
2483                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2484
2485 /*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2486 /*              Overwrite WORK(IVT) by right singular vectors of L */
2487 /*              CWorkspace: need   M*M [VT] + 2*M [tauq, taup] + M    [work] */
2488 /*              CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] */
2489 /*              RWorkspace: need   0 */
2490
2491                 zlacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2492                 i__1 = *lwork - nwork + 1;
2493                 zunmbr_("P", "R", "C", m, m, m, &a[a_offset], lda, &work[
2494                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
2495                         ierr);
2496
2497 /*              Multiply right singular vectors of L in WORK(IVT) by */
2498 /*              Q in VT, storing result in A */
2499 /*              CWorkspace: need   M*M [VT] */
2500 /*              RWorkspace: need   0 */
2501
2502                 zgemm_("N", "N", m, n, m, &c_b2, &work[ivt], &ldwkvt, &vt[
2503                         vt_offset], ldvt, &c_b1, &a[a_offset], lda);
2504
2505 /*              Copy right singular vectors of A from A to VT */
2506
2507                 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2508
2509             }
2510
2511         } else if (*n >= mnthr2) {
2512
2513 /*           MNTHR2 <= N < MNTHR1 */
2514
2515 /*           Path 5t (N >> M, but not as much as MNTHR1) */
2516 /*           Reduce to bidiagonal form without QR decomposition, use */
2517 /*           ZUNGBR and matrix multiplication to compute singular vectors */
2518
2519             ie = 1;
2520             nrwork = ie + *m;
2521             itauq = 1;
2522             itaup = itauq + *m;
2523             nwork = itaup + *m;
2524
2525 /*           Bidiagonalize A */
2526 /*           CWorkspace: need   2*M [tauq, taup] + N        [work] */
2527 /*           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] */
2528 /*           RWorkspace: need   M [e] */
2529
2530             i__1 = *lwork - nwork + 1;
2531             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2532                     &work[itaup], &work[nwork], &i__1, &ierr);
2533
2534             if (wntqn) {
2535
2536 /*              Path 5tn (N >> M, JOBZ='N') */
2537 /*              Compute singular values only */
2538 /*              CWorkspace: need   0 */
2539 /*              RWorkspace: need   M [e] + BDSPAC */
2540
2541                 dbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2542                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2543             } else if (wntqo) {
2544                 irvt = nrwork;
2545                 iru = irvt + *m * *m;
2546                 nrwork = iru + *m * *m;
2547                 ivt = nwork;
2548
2549 /*              Path 5to (N >> M, JOBZ='O') */
2550 /*              Copy A to U, generate Q */
2551 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2552 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2553 /*              RWorkspace: need   0 */
2554
2555                 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2556                 i__1 = *lwork - nwork + 1;
2557                 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2558                         nwork], &i__1, &ierr);
2559
2560 /*              Generate P**H in A */
2561 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2562 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2563 /*              RWorkspace: need   0 */
2564
2565                 i__1 = *lwork - nwork + 1;
2566                 zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
2567                         nwork], &i__1, &ierr);
2568
2569                 ldwkvt = *m;
2570                 if (*lwork >= *m * *n + *m * 3) {
2571
2572 /*                 WORK( IVT ) is M by N */
2573
2574                     nwork = ivt + ldwkvt * *n;
2575                     chunk = *n;
2576                 } else {
2577
2578 /*                 WORK( IVT ) is M by CHUNK */
2579
2580                     chunk = (*lwork - *m * 3) / *m;
2581                     nwork = ivt + ldwkvt * chunk;
2582                 }
2583
2584 /*              Perform bidiagonal SVD, computing left singular vectors */
2585 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2586 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2587 /*              CWorkspace: need   0 */
2588 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2589
2590                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2591                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2592                         info);
2593
2594 /*              Multiply Q in U by real matrix RWORK(IRVT) */
2595 /*              storing the result in WORK(IVT), copying to U */
2596 /*              CWorkspace: need   2*M [tauq, taup] + M*M [VT] */
2597 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2598
2599                 zlacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &work[ivt], &
2600                         ldwkvt, &rwork[nrwork]);
2601                 zlacpy_("F", m, m, &work[ivt], &ldwkvt, &u[u_offset], ldu);
2602
2603 /*              Multiply RWORK(IRVT) by P**H in A, storing the */
2604 /*              result in WORK(IVT), copying to A */
2605 /*              CWorkspace: need   2*M [tauq, taup] + M*M [VT] */
2606 /*              CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] */
2607 /*              RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork] */
2608 /*              RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2609
2610                 nrwork = iru;
2611                 i__1 = *n;
2612                 i__2 = chunk;
2613                 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
2614                         i__2) {
2615 /* Computing MIN */
2616                     i__3 = *n - i__ + 1;
2617                     blk = f2cmin(i__3,chunk);
2618                     zlarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1],
2619                             lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
2620                     zlacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ *
2621                             a_dim1 + 1], lda);
2622 /* L50: */
2623                 }
2624             } else if (wntqs) {
2625
2626 /*              Path 5ts (N >> M, JOBZ='S') */
2627 /*              Copy A to U, generate Q */
2628 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2629 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2630 /*              RWorkspace: need   0 */
2631
2632                 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2633                 i__2 = *lwork - nwork + 1;
2634                 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2635                         nwork], &i__2, &ierr);
2636
2637 /*              Copy A to VT, generate P**H */
2638 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2639 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2640 /*              RWorkspace: need   0 */
2641
2642                 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2643                 i__2 = *lwork - nwork + 1;
2644                 zungbr_("P", m, n, m, &vt[vt_offset], ldvt, &work[itaup], &
2645                         work[nwork], &i__2, &ierr);
2646
2647 /*              Perform bidiagonal SVD, computing left singular vectors */
2648 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2649 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2650 /*              CWorkspace: need   0 */
2651 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2652
2653                 irvt = nrwork;
2654                 iru = irvt + *m * *m;
2655                 nrwork = iru + *m * *m;
2656                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2657                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2658                         info);
2659
2660 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
2661 /*              result in A, copying to U */
2662 /*              CWorkspace: need   0 */
2663 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2664
2665                 zlacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset],
2666                          lda, &rwork[nrwork]);
2667                 zlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2668
2669 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
2670 /*              storing the result in A, copying to VT */
2671 /*              CWorkspace: need   0 */
2672 /*              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2673
2674                 nrwork = iru;
2675                 zlarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
2676                         a_offset], lda, &rwork[nrwork]);
2677                 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2678             } else {
2679
2680 /*              Path 5ta (N >> M, JOBZ='A') */
2681 /*              Copy A to U, generate Q */
2682 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2683 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2684 /*              RWorkspace: need   0 */
2685
2686                 zlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2687                 i__2 = *lwork - nwork + 1;
2688                 zungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2689                         nwork], &i__2, &ierr);
2690
2691 /*              Copy A to VT, generate P**H */
2692 /*              CWorkspace: need   2*M [tauq, taup] + N    [work] */
2693 /*              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] */
2694 /*              RWorkspace: need   0 */
2695
2696                 zlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2697                 i__2 = *lwork - nwork + 1;
2698                 zungbr_("P", n, n, m, &vt[vt_offset], ldvt, &work[itaup], &
2699                         work[nwork], &i__2, &ierr);
2700
2701 /*              Perform bidiagonal SVD, computing left singular vectors */
2702 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2703 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2704 /*              CWorkspace: need   0 */
2705 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2706
2707                 irvt = nrwork;
2708                 iru = irvt + *m * *m;
2709                 nrwork = iru + *m * *m;
2710                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2711                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2712                         info);
2713
2714 /*              Multiply Q in U by real matrix RWORK(IRU), storing the */
2715 /*              result in A, copying to U */
2716 /*              CWorkspace: need   0 */
2717 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2718
2719                 zlacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset],
2720                          lda, &rwork[nrwork]);
2721                 zlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2722
2723 /*              Multiply real matrix RWORK(IRVT) by P**H in VT, */
2724 /*              storing the result in A, copying to VT */
2725 /*              CWorkspace: need   0 */
2726 /*              RWorkspace: need   M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2727
2728                 nrwork = iru;
2729                 zlarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
2730                         a_offset], lda, &rwork[nrwork]);
2731                 zlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2732             }
2733
2734         } else {
2735
2736 /*           N .LT. MNTHR2 */
2737
2738 /*           Path 6t (N > M, but not much larger) */
2739 /*           Reduce to bidiagonal form without LQ decomposition */
2740 /*           Use ZUNMBR to compute singular vectors */
2741
2742             ie = 1;
2743             nrwork = ie + *m;
2744             itauq = 1;
2745             itaup = itauq + *m;
2746             nwork = itaup + *m;
2747
2748 /*           Bidiagonalize A */
2749 /*           CWorkspace: need   2*M [tauq, taup] + N        [work] */
2750 /*           CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] */
2751 /*           RWorkspace: need   M [e] */
2752
2753             i__2 = *lwork - nwork + 1;
2754             zgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2755                     &work[itaup], &work[nwork], &i__2, &ierr);
2756             if (wntqn) {
2757
2758 /*              Path 6tn (N > M, JOBZ='N') */
2759 /*              Compute singular values only */
2760 /*              CWorkspace: need   0 */
2761 /*              RWorkspace: need   M [e] + BDSPAC */
2762
2763                 dbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2764                         c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2765             } else if (wntqo) {
2766 /*              Path 6to (N > M, JOBZ='O') */
2767                 ldwkvt = *m;
2768                 ivt = nwork;
2769                 if (*lwork >= *m * *n + *m * 3) {
2770
2771 /*                 WORK( IVT ) is M by N */
2772
2773                     zlaset_("F", m, n, &c_b1, &c_b1, &work[ivt], &ldwkvt);
2774                     nwork = ivt + ldwkvt * *n;
2775                 } else {
2776
2777 /*                 WORK( IVT ) is M by CHUNK */
2778
2779                     chunk = (*lwork - *m * 3) / *m;
2780                     nwork = ivt + ldwkvt * chunk;
2781                 }
2782
2783 /*              Perform bidiagonal SVD, computing left singular vectors */
2784 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2785 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2786 /*              CWorkspace: need   0 */
2787 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2788
2789                 irvt = nrwork;
2790                 iru = irvt + *m * *m;
2791                 nrwork = iru + *m * *m;
2792                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2793                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2794                         info);
2795
2796 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2797 /*              Overwrite U by left singular vectors of A */
2798 /*              CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work] */
2799 /*              CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] */
2800 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] */
2801
2802                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2803                 i__2 = *lwork - nwork + 1;
2804                 zunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2805                         itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2806
2807                 if (*lwork >= *m * *n + *m * 3) {
2808
2809 /*                 Path 6to-fast */
2810 /*                 Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2811 /*                 Overwrite WORK(IVT) by right singular vectors of A, */
2812 /*                 copying to A */
2813 /*                 CWorkspace: need   2*M [tauq, taup] + M*N [VT] + M    [work] */
2814 /*                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work] */
2815 /*                 RWorkspace: need   M [e] + M*M [RVT] */
2816
2817                     zlacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2818                     i__2 = *lwork - nwork + 1;
2819                     zunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
2820                             itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2,
2821                             &ierr);
2822                     zlacpy_("F", m, n, &work[ivt], &ldwkvt, &a[a_offset], lda);
2823                 } else {
2824
2825 /*                 Path 6to-slow */
2826 /*                 Generate P**H in A */
2827 /*                 CWorkspace: need   2*M [tauq, taup] + M*M [VT] + M    [work] */
2828 /*                 CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] */
2829 /*                 RWorkspace: need   0 */
2830
2831                     i__2 = *lwork - nwork + 1;
2832                     zungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
2833                             work[nwork], &i__2, &ierr);
2834
2835 /*                 Multiply Q in A by real matrix RWORK(IRU), storing the */
2836 /*                 result in WORK(IU), copying to A */
2837 /*                 CWorkspace: need   2*M [tauq, taup] + M*M [VT] */
2838 /*                 CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] */
2839 /*                 RWorkspace: need   M [e] + M*M [RVT] + 2*M*M [rwork] */
2840 /*                 RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2841
2842                     nrwork = iru;
2843                     i__2 = *n;
2844                     i__1 = chunk;
2845                     for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2846                              i__1) {
2847 /* Computing MIN */
2848                         i__3 = *n - i__ + 1;
2849                         blk = f2cmin(i__3,chunk);
2850                         zlarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1]
2851                                 , lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
2852                         zlacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ *
2853                                 a_dim1 + 1], lda);
2854 /* L60: */
2855                     }
2856                 }
2857             } else if (wntqs) {
2858
2859 /*              Path 6ts (N > M, JOBZ='S') */
2860 /*              Perform bidiagonal SVD, computing left singular vectors */
2861 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2862 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2863 /*              CWorkspace: need   0 */
2864 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2865
2866                 irvt = nrwork;
2867                 iru = irvt + *m * *m;
2868                 nrwork = iru + *m * *m;
2869                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2870                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2871                         info);
2872
2873 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2874 /*              Overwrite U by left singular vectors of A */
2875 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2876 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2877 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] */
2878
2879                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2880                 i__1 = *lwork - nwork + 1;
2881                 zunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2882                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2883
2884 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
2885 /*              Overwrite VT by right singular vectors of A */
2886 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2887 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2888 /*              RWorkspace: need   M [e] + M*M [RVT] */
2889
2890                 zlaset_("F", m, n, &c_b1, &c_b1, &vt[vt_offset], ldvt);
2891                 zlacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2892                 i__1 = *lwork - nwork + 1;
2893                 zunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
2894                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2895                         ierr);
2896             } else {
2897
2898 /*              Path 6ta (N > M, JOBZ='A') */
2899 /*              Perform bidiagonal SVD, computing left singular vectors */
2900 /*              of bidiagonal matrix in RWORK(IRU) and computing right */
2901 /*              singular vectors of bidiagonal matrix in RWORK(IRVT) */
2902 /*              CWorkspace: need   0 */
2903 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2904
2905                 irvt = nrwork;
2906                 iru = irvt + *m * *m;
2907                 nrwork = iru + *m * *m;
2908
2909                 dbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2910                         rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2911                         info);
2912
2913 /*              Copy real matrix RWORK(IRU) to complex matrix U */
2914 /*              Overwrite U by left singular vectors of A */
2915 /*              CWorkspace: need   2*M [tauq, taup] + M    [work] */
2916 /*              CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2917 /*              RWorkspace: need   M [e] + M*M [RVT] + M*M [RU] */
2918
2919                 zlacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2920                 i__1 = *lwork - nwork + 1;
2921                 zunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2922                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2923
2924 /*              Set all of VT to identity matrix */
2925
2926                 zlaset_("F", n, n, &c_b1, &c_b2, &vt[vt_offset], ldvt);
2927
2928 /*              Copy real matrix RWORK(IRVT) to complex matrix VT */
2929 /*              Overwrite VT by right singular vectors of A */
2930 /*              CWorkspace: need   2*M [tauq, taup] + N    [work] */
2931 /*              CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] */
2932 /*              RWorkspace: need   M [e] + M*M [RVT] */
2933
2934                 zlacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2935                 i__1 = *lwork - nwork + 1;
2936                 zunmbr_("P", "R", "C", n, n, m, &a[a_offset], lda, &work[
2937                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2938                         ierr);
2939             }
2940
2941         }
2942
2943     }
2944
2945 /*     Undo scaling if necessary */
2946
2947     if (iscl == 1) {
2948         if (anrm > bignum) {
2949             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
2950                     minmn, &ierr);
2951         }
2952         if (*info != 0 && anrm > bignum) {
2953             i__1 = minmn - 1;
2954             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &i__1, &c__1, &rwork[
2955                     ie], &minmn, &ierr);
2956         }
2957         if (anrm < smlnum) {
2958             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
2959                     minmn, &ierr);
2960         }
2961         if (*info != 0 && anrm < smlnum) {
2962             i__1 = minmn - 1;
2963             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__1, &c__1, &rwork[
2964                     ie], &minmn, &ierr);
2965         }
2966     }
2967
2968 /*     Return optimal workspace in WORK(1) */
2969
2970     work[1].r = (doublereal) maxwrk, work[1].i = 0.;
2971
2972     return 0;
2973
2974 /*     End of ZGESDD */
2975
2976 } /* zgesdd_ */
2977