lapack-netlib/SRC/dgesvd.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 integer c__6 = 6;
 516 static integer c__0 = 0;
 517 static integer c__2 = 2;
 518 static integer c_n1 = -1;
 519 static doublereal c_b57 = 0.;
 520 static integer c__1 = 1;
 521 static doublereal c_b79 = 1.;
 522
 523 /* > \brief <b> DGESVD computes the singular value decomposition (SVD) for GE matrices</b> */
 524
 525 /*  =========== DOCUMENTATION =========== */
 526
 527 /* Online html documentation available at */
 528 /*            http://www.netlib.org/lapack/explore-html/ */
 529
 530 /* > \htmlonly */
 531 /* > Download DGESVD + dependencies */
 532 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.
 533 f"> */
 534 /* > [TGZ]</a> */
 535 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.
 536 f"> */
 537 /* > [ZIP]</a> */
 538 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.
 539 f"> */
 540 /* > [TXT]</a> */
 541 /* > \endhtmlonly */
 542
 543 /*  Definition: */
 544 /*  =========== */
 545
 546 /*       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, */
 547 /*                          WORK, LWORK, INFO ) */
 548
 549 /*       CHARACTER          JOBU, JOBVT */
 550 /*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N */
 551 /*       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ), */
 552 /*      $                   VT( LDVT, * ), WORK( * ) */
 553
 554
 555 /* > \par Purpose: */
 556 /*  ============= */
 557 /* > */
 558 /* > \verbatim */
 559 /* > */
 560 /* > DGESVD computes the singular value decomposition (SVD) of a real */
 561 /* > M-by-N matrix A, optionally computing the left and/or right singular */
 562 /* > vectors. The SVD is written */
 563 /* > */
 564 /* >      A = U * SIGMA * 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 orthogonal matrix, and */
 568 /* > V is an N-by-N orthogonal 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 V**T, not V. */
 574 /* > \endverbatim */
 575
 576 /*  Arguments: */
 577 /*  ========== */
 578
 579 /* > \param[in] JOBU */
 580 /* > \verbatim */
 581 /* >          JOBU is CHARACTER*1 */
 582 /* >          Specifies options for computing all or part of the matrix U: */
 583 /* >          = 'A':  all M columns of U are returned in array U: */
 584 /* >          = 'S':  the first f2cmin(m,n) columns of U (the left singular */
 585 /* >                  vectors) are returned in the array U; */
 586 /* >          = 'O':  the first f2cmin(m,n) columns of U (the left singular */
 587 /* >                  vectors) are overwritten on the array A; */
 588 /* >          = 'N':  no columns of U (no left singular vectors) are */
 589 /* >                  computed. */
 590 /* > \endverbatim */
 591 /* > */
 592 /* > \param[in] JOBVT */
 593 /* > \verbatim */
 594 /* >          JOBVT is CHARACTER*1 */
 595 /* >          Specifies options for computing all or part of the matrix */
 596 /* >          V**T: */
 597 /* >          = 'A':  all N rows of V**T are returned in the array VT; */
 598 /* >          = 'S':  the first f2cmin(m,n) rows of V**T (the right singular */
 599 /* >                  vectors) are returned in the array VT; */
 600 /* >          = 'O':  the first f2cmin(m,n) rows of V**T (the right singular */
 601 /* >                  vectors) are overwritten on the array A; */
 602 /* >          = 'N':  no rows of V**T (no right singular vectors) are */
 603 /* >                  computed. */
 604 /* > */
 605 /* >          JOBVT and JOBU cannot both be 'O'. */
 606 /* > \endverbatim */
 607 /* > */
 608 /* > \param[in] M */
 609 /* > \verbatim */
 610 /* >          M is INTEGER */
 611 /* >          The number of rows of the input matrix A.  M >= 0. */
 612 /* > \endverbatim */
 613 /* > */
 614 /* > \param[in] N */
 615 /* > \verbatim */
 616 /* >          N is INTEGER */
 617 /* >          The number of columns of the input matrix A.  N >= 0. */
 618 /* > \endverbatim */
 619 /* > */
 620 /* > \param[in,out] A */
 621 /* > \verbatim */
 622 /* >          A is DOUBLE PRECISION array, dimension (LDA,N) */
 623 /* >          On entry, the M-by-N matrix A. */
 624 /* >          On exit, */
 625 /* >          if JOBU = 'O',  A is overwritten with the first f2cmin(m,n) */
 626 /* >                          columns of U (the left singular vectors, */
 627 /* >                          stored columnwise); */
 628 /* >          if JOBVT = 'O', A is overwritten with the first f2cmin(m,n) */
 629 /* >                          rows of V**T (the right singular vectors, */
 630 /* >                          stored rowwise); */
 631 /* >          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A */
 632 /* >                          are destroyed. */
 633 /* > \endverbatim */
 634 /* > */
 635 /* > \param[in] LDA */
 636 /* > \verbatim */
 637 /* >          LDA is INTEGER */
 638 /* >          The leading dimension of the array A.  LDA >= f2cmax(1,M). */
 639 /* > \endverbatim */
 640 /* > */
 641 /* > \param[out] S */
 642 /* > \verbatim */
 643 /* >          S is DOUBLE PRECISION array, dimension (f2cmin(M,N)) */
 644 /* >          The singular values of A, sorted so that S(i) >= S(i+1). */
 645 /* > \endverbatim */
 646 /* > */
 647 /* > \param[out] U */
 648 /* > \verbatim */
 649 /* >          U is DOUBLE PRECISION array, dimension (LDU,UCOL) */
 650 /* >          (LDU,M) if JOBU = 'A' or (LDU,f2cmin(M,N)) if JOBU = 'S'. */
 651 /* >          If JOBU = 'A', U contains the M-by-M orthogonal matrix U; */
 652 /* >          if JOBU = 'S', U contains the first f2cmin(m,n) columns of U */
 653 /* >          (the left singular vectors, stored columnwise); */
 654 /* >          if JOBU = 'N' or 'O', U is not referenced. */
 655 /* > \endverbatim */
 656 /* > */
 657 /* > \param[in] LDU */
 658 /* > \verbatim */
 659 /* >          LDU is INTEGER */
 660 /* >          The leading dimension of the array U.  LDU >= 1; if */
 661 /* >          JOBU = 'S' or 'A', LDU >= M. */
 662 /* > \endverbatim */
 663 /* > */
 664 /* > \param[out] VT */
 665 /* > \verbatim */
 666 /* >          VT is DOUBLE PRECISION array, dimension (LDVT,N) */
 667 /* >          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix */
 668 /* >          V**T; */
 669 /* >          if JOBVT = 'S', VT contains the first f2cmin(m,n) rows of */
 670 /* >          V**T (the right singular vectors, stored rowwise); */
 671 /* >          if JOBVT = 'N' or 'O', VT is not referenced. */
 672 /* > \endverbatim */
 673 /* > */
 674 /* > \param[in] LDVT */
 675 /* > \verbatim */
 676 /* >          LDVT is INTEGER */
 677 /* >          The leading dimension of the array VT.  LDVT >= 1; if */
 678 /* >          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= f2cmin(M,N). */
 679 /* > \endverbatim */
 680 /* > */
 681 /* > \param[out] WORK */
 682 /* > \verbatim */
 683 /* >          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
 684 /* >          On exit, if INFO = 0, WORK(1) returns the optimal LWORK; */
 685 /* >          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged */
 686 /* >          superdiagonal elements of an upper bidiagonal matrix B */
 687 /* >          whose diagonal is in S (not necessarily sorted). B */
 688 /* >          satisfies A = U * B * VT, so it has the same singular values */
 689 /* >          as A, and singular vectors related by U and VT. */
 690 /* > \endverbatim */
 691 /* > */
 692 /* > \param[in] LWORK */
 693 /* > \verbatim */
 694 /* >          LWORK is INTEGER */
 695 /* >          The dimension of the array WORK. */
 696 /* >          LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code): */
 697 /* >             - PATH 1  (M much larger than N, JOBU='N') */
 698 /* >             - PATH 1t (N much larger than M, JOBVT='N') */
 699 /* >          LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths */
 700 /* >          For good performance, LWORK should generally be larger. */
 701 /* > */
 702 /* >          If LWORK = -1, then a workspace query is assumed; the routine */
 703 /* >          only calculates the optimal size of the WORK array, returns */
 704 /* >          this value as the first entry of the WORK array, and no error */
 705 /* >          message related to LWORK is issued by XERBLA. */
 706 /* > \endverbatim */
 707 /* > */
 708 /* > \param[out] INFO */
 709 /* > \verbatim */
 710 /* >          INFO is INTEGER */
 711 /* >          = 0:  successful exit. */
 712 /* >          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 713 /* >          > 0:  if DBDSQR did not converge, INFO specifies how many */
 714 /* >                superdiagonals of an intermediate bidiagonal form B */
 715 /* >                did not converge to zero. See the description of WORK */
 716 /* >                above for details. */
 717 /* > \endverbatim */
 718
 719 /*  Authors: */
 720 /*  ======== */
 721
 722 /* > \author Univ. of Tennessee */
 723 /* > \author Univ. of California Berkeley */
 724 /* > \author Univ. of Colorado Denver */
 725 /* > \author NAG Ltd. */
 726
 727 /* > \date April 2012 */
 728
 729 /* > \ingroup doubleGEsing */
 730
 731 /*  ===================================================================== */
 732 /* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
 733         doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
 734         ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
 735         integer *info)
 736 {
 737     /* System generated locals */
 738     address a__1[2];
 739     integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1[2],
 740             i__2, i__3, i__4;
 741     char ch__1[2];
 742
 743     /* Local variables */
 744     integer iscl;
 745     doublereal anrm;
 746     integer ierr, itau, ncvt, nrvt, lwork_dgebrd__, lwork_dgelqf__,
 747             lwork_dgeqrf__, i__;
 748     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
 749             integer *, doublereal *, doublereal *, integer *, doublereal *,
 750             integer *, doublereal *, doublereal *, integer *);
 751     extern logical lsame_(char *, char *);
 752     integer chunk, minmn, wrkbl, itaup, itauq, mnthr, iwork;
 753     logical wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
 754     integer ie;
 755     extern /* Subroutine */ int dgebrd_(integer *, integer *, doublereal *,
 756             integer *, doublereal *, doublereal *, doublereal *, doublereal *,
 757              doublereal *, integer *, integer *);
 758     extern doublereal dlamch_(char *), dlange_(char *, integer *,
 759             integer *, doublereal *, integer *, doublereal *);
 760     integer ir, bdspac, iu;
 761     extern /* Subroutine */ int dgelqf_(integer *, integer *, doublereal *,
 762             integer *, doublereal *, doublereal *, integer *, integer *),
 763             dlascl_(char *, integer *, integer *, doublereal *, doublereal *,
 764             integer *, integer *, doublereal *, integer *, integer *),
 765              dgeqrf_(integer *, integer *, doublereal *, integer *,
 766             doublereal *, doublereal *, integer *, integer *), dlacpy_(char *,
 767              integer *, integer *, doublereal *, integer *, doublereal *,
 768             integer *), dlaset_(char *, integer *, integer *,
 769             doublereal *, doublereal *, doublereal *, integer *),
 770             dbdsqr_(char *, integer *, integer *, integer *, integer *,
 771             doublereal *, doublereal *, doublereal *, integer *, doublereal *,
 772              integer *, doublereal *, integer *, doublereal *, integer *), dorgbr_(char *, integer *, integer *, integer *,
 773             doublereal *, integer *, doublereal *, doublereal *, integer *,
 774             integer *);
 775     doublereal bignum;
 776     extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
 777     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
 778             integer *, integer *, ftnlen, ftnlen);
 779     extern /* Subroutine */ int dormbr_(char *, char *, char *, integer *,
 780             integer *, integer *, doublereal *, integer *, doublereal *,
 781             doublereal *, integer *, doublereal *, integer *, integer *), dorglq_(integer *, integer *, integer *,
 782             doublereal *, integer *, doublereal *, doublereal *, integer *,
 783             integer *), dorgqr_(integer *, integer *, integer *, doublereal *,
 784              integer *, doublereal *, doublereal *, integer *, integer *);
 785     integer ldwrkr, minwrk, ldwrku, maxwrk;
 786     doublereal smlnum;
 787     logical lquery, wntuas, wntvas;
 788     integer lwork_dorgbr_p__, lwork_dorgbr_q__, lwork_dorglq_m__,
 789             lwork_dorglq_n__, lwork_dorgqr_m__, lwork_dorgqr_n__, blk, ncu;
 790     doublereal dum[1], eps;
 791     integer nru;
 792
 793
 794 /*  -- LAPACK driver routine (version 3.7.0) -- */
 795 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
 796 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
 797 /*     April 2012 */
 798
 799
 800 /*  ===================================================================== */
 801
 802
 803 /*     Test the input arguments */
 804
 805     /* Parameter adjustments */
 806     a_dim1 = *lda;
 807     a_offset = 1 + a_dim1 * 1;
 808     a -= a_offset;
 809     --s;
 810     u_dim1 = *ldu;
 811     u_offset = 1 + u_dim1 * 1;
 812     u -= u_offset;
 813     vt_dim1 = *ldvt;
 814     vt_offset = 1 + vt_dim1 * 1;
 815     vt -= vt_offset;
 816     --work;
 817
 818     /* Function Body */
 819     *info = 0;
 820     minmn = f2cmin(*m,*n);
 821     wntua = lsame_(jobu, "A");
 822     wntus = lsame_(jobu, "S");
 823     wntuas = wntua || wntus;
 824     wntuo = lsame_(jobu, "O");
 825     wntun = lsame_(jobu, "N");
 826     wntva = lsame_(jobvt, "A");
 827     wntvs = lsame_(jobvt, "S");
 828     wntvas = wntva || wntvs;
 829     wntvo = lsame_(jobvt, "O");
 830     wntvn = lsame_(jobvt, "N");
 831     lquery = *lwork == -1;
 832
 833     if (! (wntua || wntus || wntuo || wntun)) {
 834         *info = -1;
 835     } else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo) {
 836         *info = -2;
 837     } else if (*m < 0) {
 838         *info = -3;
 839     } else if (*n < 0) {
 840         *info = -4;
 841     } else if (*lda < f2cmax(1,*m)) {
 842         *info = -6;
 843     } else if (*ldu < 1 || wntuas && *ldu < *m) {
 844         *info = -9;
 845     } else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn) {
 846         *info = -11;
 847     }
 848
 849 /*     Compute workspace */
 850 /*      (Note: Comments in the code beginning "Workspace:" describe the */
 851 /*       minimal amount of workspace needed at that point in the code, */
 852 /*       as well as the preferred amount for good performance. */
 853 /*       NB refers to the optimal block size for the immediately */
 854 /*       following subroutine, as returned by ILAENV.) */
 855
 856     if (*info == 0) {
 857         minwrk = 1;
 858         maxwrk = 1;
 859         if (*m >= *n && minmn > 0) {
 860
 861 /*           Compute space needed for DBDSQR */
 862
 863 /* Writing concatenation */
 864             i__1[0] = 1, a__1[0] = jobu;
 865             i__1[1] = 1, a__1[1] = jobvt;
 866             s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
 867             mnthr = ilaenv_(&c__6, "DGESVD", ch__1, m, n, &c__0, &c__0, (
 868                     ftnlen)6, (ftnlen)2);
 869             bdspac = *n * 5;
 870 /*           Compute space needed for DGEQRF */
 871             dgeqrf_(m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
 872             lwork_dgeqrf__ = (integer) dum[0];
 873 /*           Compute space needed for DORGQR */
 874             dorgqr_(m, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
 875             lwork_dorgqr_n__ = (integer) dum[0];
 876             dorgqr_(m, m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
 877             lwork_dorgqr_m__ = (integer) dum[0];
 878 /*           Compute space needed for DGEBRD */
 879             dgebrd_(n, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &c_n1,
 880                      &ierr);
 881             lwork_dgebrd__ = (integer) dum[0];
 882 /*           Compute space needed for DORGBR P */
 883             dorgbr_("P", n, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
 884             lwork_dorgbr_p__ = (integer) dum[0];
 885 /*           Compute space needed for DORGBR Q */
 886             dorgbr_("Q", n, n, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
 887             lwork_dorgbr_q__ = (integer) dum[0];
 888
 889             if (*m >= mnthr) {
 890                 if (wntun) {
 891
 892 /*                 Path 1 (M much larger than N, JOBU='N') */
 893
 894                     maxwrk = *n + lwork_dgeqrf__;
 895 /* Computing MAX */
 896                     i__2 = maxwrk, i__3 = *n * 3 + lwork_dgebrd__;
 897                     maxwrk = f2cmax(i__2,i__3);
 898                     if (wntvo || wntvas) {
 899 /* Computing MAX */
 900                         i__2 = maxwrk, i__3 = *n * 3 + lwork_dorgbr_p__;
 901                         maxwrk = f2cmax(i__2,i__3);
 902                     }
 903                     maxwrk = f2cmax(maxwrk,bdspac);
 904 /* Computing MAX */
 905                     i__2 = *n << 2;
 906                     minwrk = f2cmax(i__2,bdspac);
 907                 } else if (wntuo && wntvn) {
 908
 909 /*                 Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
 910
 911                     wrkbl = *n + lwork_dgeqrf__;
 912 /* Computing MAX */
 913                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_n__;
 914                     wrkbl = f2cmax(i__2,i__3);
 915 /* Computing MAX */
 916                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
 917                     wrkbl = f2cmax(i__2,i__3);
 918 /* Computing MAX */
 919                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
 920                     wrkbl = f2cmax(i__2,i__3);
 921                     wrkbl = f2cmax(wrkbl,bdspac);
 922 /* Computing MAX */
 923                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
 924                     maxwrk = f2cmax(i__2,i__3);
 925 /* Computing MAX */
 926                     i__2 = *n * 3 + *m;
 927                     minwrk = f2cmax(i__2,bdspac);
 928                 } else if (wntuo && wntvas) {
 929
 930 /*                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or */
 931 /*                 'A') */
 932
 933                     wrkbl = *n + lwork_dgeqrf__;
 934 /* Computing MAX */
 935                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_n__;
 936                     wrkbl = f2cmax(i__2,i__3);
 937 /* Computing MAX */
 938                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
 939                     wrkbl = f2cmax(i__2,i__3);
 940 /* Computing MAX */
 941                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
 942                     wrkbl = f2cmax(i__2,i__3);
 943 /* Computing MAX */
 944                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_p__;
 945                     wrkbl = f2cmax(i__2,i__3);
 946                     wrkbl = f2cmax(wrkbl,bdspac);
 947 /* Computing MAX */
 948                     i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
 949                     maxwrk = f2cmax(i__2,i__3);
 950 /* Computing MAX */
 951                     i__2 = *n * 3 + *m;
 952                     minwrk = f2cmax(i__2,bdspac);
 953                 } else if (wntus && wntvn) {
 954
 955 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
 956
 957                     wrkbl = *n + lwork_dgeqrf__;
 958 /* Computing MAX */
 959                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_n__;
 960                     wrkbl = f2cmax(i__2,i__3);
 961 /* Computing MAX */
 962                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
 963                     wrkbl = f2cmax(i__2,i__3);
 964 /* Computing MAX */
 965                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
 966                     wrkbl = f2cmax(i__2,i__3);
 967                     wrkbl = f2cmax(wrkbl,bdspac);
 968                     maxwrk = *n * *n + wrkbl;
 969 /* Computing MAX */
 970                     i__2 = *n * 3 + *m;
 971                     minwrk = f2cmax(i__2,bdspac);
 972                 } else if (wntus && wntvo) {
 973
 974 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
 975
 976                     wrkbl = *n + lwork_dgeqrf__;
 977 /* Computing MAX */
 978                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_n__;
 979                     wrkbl = f2cmax(i__2,i__3);
 980 /* Computing MAX */
 981                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
 982                     wrkbl = f2cmax(i__2,i__3);
 983 /* Computing MAX */
 984                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
 985                     wrkbl = f2cmax(i__2,i__3);
 986 /* Computing MAX */
 987                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_p__;
 988                     wrkbl = f2cmax(i__2,i__3);
 989                     wrkbl = f2cmax(wrkbl,bdspac);
 990                     maxwrk = (*n << 1) * *n + wrkbl;
 991 /* Computing MAX */
 992                     i__2 = *n * 3 + *m;
 993                     minwrk = f2cmax(i__2,bdspac);
 994                 } else if (wntus && wntvas) {
 995
 996 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or */
 997 /*                 'A') */
 998
 999                     wrkbl = *n + lwork_dgeqrf__;
1000 /* Computing MAX */
1001                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_n__;
1002                     wrkbl = f2cmax(i__2,i__3);
1003 /* Computing MAX */
1004                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
1005                     wrkbl = f2cmax(i__2,i__3);
1006 /* Computing MAX */
1007                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
1008                     wrkbl = f2cmax(i__2,i__3);
1009 /* Computing MAX */
1010                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_p__;
1011                     wrkbl = f2cmax(i__2,i__3);
1012                     wrkbl = f2cmax(wrkbl,bdspac);
1013                     maxwrk = *n * *n + wrkbl;
1014 /* Computing MAX */
1015                     i__2 = *n * 3 + *m;
1016                     minwrk = f2cmax(i__2,bdspac);
1017                 } else if (wntua && wntvn) {
1018
1019 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
1020
1021                     wrkbl = *n + lwork_dgeqrf__;
1022 /* Computing MAX */
1023                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_m__;
1024                     wrkbl = f2cmax(i__2,i__3);
1025 /* Computing MAX */
1026                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
1027                     wrkbl = f2cmax(i__2,i__3);
1028 /* Computing MAX */
1029                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
1030                     wrkbl = f2cmax(i__2,i__3);
1031                     wrkbl = f2cmax(wrkbl,bdspac);
1032                     maxwrk = *n * *n + wrkbl;
1033 /* Computing MAX */
1034                     i__2 = *n * 3 + *m;
1035                     minwrk = f2cmax(i__2,bdspac);
1036                 } else if (wntua && wntvo) {
1037
1038 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
1039
1040                     wrkbl = *n + lwork_dgeqrf__;
1041 /* Computing MAX */
1042                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_m__;
1043                     wrkbl = f2cmax(i__2,i__3);
1044 /* Computing MAX */
1045                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
1046                     wrkbl = f2cmax(i__2,i__3);
1047 /* Computing MAX */
1048                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
1049                     wrkbl = f2cmax(i__2,i__3);
1050 /* Computing MAX */
1051                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_p__;
1052                     wrkbl = f2cmax(i__2,i__3);
1053                     wrkbl = f2cmax(wrkbl,bdspac);
1054                     maxwrk = (*n << 1) * *n + wrkbl;
1055 /* Computing MAX */
1056                     i__2 = *n * 3 + *m;
1057                     minwrk = f2cmax(i__2,bdspac);
1058                 } else if (wntua && wntvas) {
1059
1060 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or */
1061 /*                 'A') */
1062
1063                     wrkbl = *n + lwork_dgeqrf__;
1064 /* Computing MAX */
1065                     i__2 = wrkbl, i__3 = *n + lwork_dorgqr_m__;
1066                     wrkbl = f2cmax(i__2,i__3);
1067 /* Computing MAX */
1068                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dgebrd__;
1069                     wrkbl = f2cmax(i__2,i__3);
1070 /* Computing MAX */
1071                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_q__;
1072                     wrkbl = f2cmax(i__2,i__3);
1073 /* Computing MAX */
1074                     i__2 = wrkbl, i__3 = *n * 3 + lwork_dorgbr_p__;
1075                     wrkbl = f2cmax(i__2,i__3);
1076                     wrkbl = f2cmax(wrkbl,bdspac);
1077                     maxwrk = *n * *n + wrkbl;
1078 /* Computing MAX */
1079                     i__2 = *n * 3 + *m;
1080                     minwrk = f2cmax(i__2,bdspac);
1081                 }
1082             } else {
1083
1084 /*              Path 10 (M at least N, but not much larger) */
1085
1086                 dgebrd_(m, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &
1087                         c_n1, &ierr);
1088                 lwork_dgebrd__ = (integer) dum[0];
1089                 maxwrk = *n * 3 + lwork_dgebrd__;
1090                 if (wntus || wntuo) {
1091                     dorgbr_("Q", m, n, n, &a[a_offset], lda, dum, dum, &c_n1,
1092                             &ierr);
1093                     lwork_dorgbr_q__ = (integer) dum[0];
1094 /* Computing MAX */
1095                     i__2 = maxwrk, i__3 = *n * 3 + lwork_dorgbr_q__;
1096                     maxwrk = f2cmax(i__2,i__3);
1097                 }
1098                 if (wntua) {
1099                     dorgbr_("Q", m, m, n, &a[a_offset], lda, dum, dum, &c_n1,
1100                             &ierr);
1101                     lwork_dorgbr_q__ = (integer) dum[0];
1102 /* Computing MAX */
1103                     i__2 = maxwrk, i__3 = *n * 3 + lwork_dorgbr_q__;
1104                     maxwrk = f2cmax(i__2,i__3);
1105                 }
1106                 if (! wntvn) {
1107 /* Computing MAX */
1108                     i__2 = maxwrk, i__3 = *n * 3 + lwork_dorgbr_p__;
1109                     maxwrk = f2cmax(i__2,i__3);
1110                 }
1111                 maxwrk = f2cmax(maxwrk,bdspac);
1112 /* Computing MAX */
1113                 i__2 = *n * 3 + *m;
1114                 minwrk = f2cmax(i__2,bdspac);
1115             }
1116         } else if (minmn > 0) {
1117
1118 /*           Compute space needed for DBDSQR */
1119
1120 /* Writing concatenation */
1121             i__1[0] = 1, a__1[0] = jobu;
1122             i__1[1] = 1, a__1[1] = jobvt;
1123             s_cat(ch__1, a__1, i__1, &c__2, (ftnlen)2);
1124             mnthr = ilaenv_(&c__6, "DGESVD", ch__1, m, n, &c__0, &c__0, (
1125                     ftnlen)6, (ftnlen)2);
1126             bdspac = *m * 5;
1127 /*           Compute space needed for DGELQF */
1128             dgelqf_(m, n, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
1129             lwork_dgelqf__ = (integer) dum[0];
1130 /*           Compute space needed for DORGLQ */
1131             dorglq_(n, n, m, dum, n, dum, dum, &c_n1, &ierr);
1132             lwork_dorglq_n__ = (integer) dum[0];
1133             dorglq_(m, n, m, &a[a_offset], lda, dum, dum, &c_n1, &ierr);
1134             lwork_dorglq_m__ = (integer) dum[0];
1135 /*           Compute space needed for DGEBRD */
1136             dgebrd_(m, m, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &c_n1,
1137                      &ierr);
1138             lwork_dgebrd__ = (integer) dum[0];
1139 /*            Compute space needed for DORGBR P */
1140             dorgbr_("P", m, m, m, &a[a_offset], n, dum, dum, &c_n1, &ierr);
1141             lwork_dorgbr_p__ = (integer) dum[0];
1142 /*           Compute space needed for DORGBR Q */
1143             dorgbr_("Q", m, m, m, &a[a_offset], n, dum, dum, &c_n1, &ierr);
1144             lwork_dorgbr_q__ = (integer) dum[0];
1145             if (*n >= mnthr) {
1146                 if (wntvn) {
1147
1148 /*                 Path 1t(N much larger than M, JOBVT='N') */
1149
1150                     maxwrk = *m + lwork_dgelqf__;
1151 /* Computing MAX */
1152                     i__2 = maxwrk, i__3 = *m * 3 + lwork_dgebrd__;
1153                     maxwrk = f2cmax(i__2,i__3);
1154                     if (wntuo || wntuas) {
1155 /* Computing MAX */
1156                         i__2 = maxwrk, i__3 = *m * 3 + lwork_dorgbr_q__;
1157                         maxwrk = f2cmax(i__2,i__3);
1158                     }
1159                     maxwrk = f2cmax(maxwrk,bdspac);
1160 /* Computing MAX */
1161                     i__2 = *m << 2;
1162                     minwrk = f2cmax(i__2,bdspac);
1163                 } else if (wntvo && wntun) {
1164
1165 /*                 Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
1166
1167                     wrkbl = *m + lwork_dgelqf__;
1168 /* Computing MAX */
1169                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_m__;
1170                     wrkbl = f2cmax(i__2,i__3);
1171 /* Computing MAX */
1172                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1173                     wrkbl = f2cmax(i__2,i__3);
1174 /* Computing MAX */
1175                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1176                     wrkbl = f2cmax(i__2,i__3);
1177                     wrkbl = f2cmax(wrkbl,bdspac);
1178 /* Computing MAX */
1179                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
1180                     maxwrk = f2cmax(i__2,i__3);
1181 /* Computing MAX */
1182                     i__2 = *m * 3 + *n;
1183                     minwrk = f2cmax(i__2,bdspac);
1184                 } else if (wntvo && wntuas) {
1185
1186 /*                 Path 3t(N much larger than M, JOBU='S' or 'A', */
1187 /*                 JOBVT='O') */
1188
1189                     wrkbl = *m + lwork_dgelqf__;
1190 /* Computing MAX */
1191                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_m__;
1192                     wrkbl = f2cmax(i__2,i__3);
1193 /* Computing MAX */
1194                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1195                     wrkbl = f2cmax(i__2,i__3);
1196 /* Computing MAX */
1197                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1198                     wrkbl = f2cmax(i__2,i__3);
1199 /* Computing MAX */
1200                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_q__;
1201                     wrkbl = f2cmax(i__2,i__3);
1202                     wrkbl = f2cmax(wrkbl,bdspac);
1203 /* Computing MAX */
1204                     i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
1205                     maxwrk = f2cmax(i__2,i__3);
1206 /* Computing MAX */
1207                     i__2 = *m * 3 + *n;
1208                     minwrk = f2cmax(i__2,bdspac);
1209                 } else if (wntvs && wntun) {
1210
1211 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
1212
1213                     wrkbl = *m + lwork_dgelqf__;
1214 /* Computing MAX */
1215                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_m__;
1216                     wrkbl = f2cmax(i__2,i__3);
1217 /* Computing MAX */
1218                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1219                     wrkbl = f2cmax(i__2,i__3);
1220 /* Computing MAX */
1221                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1222                     wrkbl = f2cmax(i__2,i__3);
1223                     wrkbl = f2cmax(wrkbl,bdspac);
1224                     maxwrk = *m * *m + wrkbl;
1225 /* Computing MAX */
1226                     i__2 = *m * 3 + *n;
1227                     minwrk = f2cmax(i__2,bdspac);
1228                 } else if (wntvs && wntuo) {
1229
1230 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
1231
1232                     wrkbl = *m + lwork_dgelqf__;
1233 /* Computing MAX */
1234                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_m__;
1235                     wrkbl = f2cmax(i__2,i__3);
1236 /* Computing MAX */
1237                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1238                     wrkbl = f2cmax(i__2,i__3);
1239 /* Computing MAX */
1240                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1241                     wrkbl = f2cmax(i__2,i__3);
1242 /* Computing MAX */
1243                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_q__;
1244                     wrkbl = f2cmax(i__2,i__3);
1245                     wrkbl = f2cmax(wrkbl,bdspac);
1246                     maxwrk = (*m << 1) * *m + wrkbl;
1247 /* Computing MAX */
1248                     i__2 = *m * 3 + *n;
1249                     minwrk = f2cmax(i__2,bdspac);
1250                 } else if (wntvs && wntuas) {
1251
1252 /*                 Path 6t(N much larger than M, JOBU='S' or 'A', */
1253 /*                 JOBVT='S') */
1254
1255                     wrkbl = *m + lwork_dgelqf__;
1256 /* Computing MAX */
1257                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_m__;
1258                     wrkbl = f2cmax(i__2,i__3);
1259 /* Computing MAX */
1260                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1261                     wrkbl = f2cmax(i__2,i__3);
1262 /* Computing MAX */
1263                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1264                     wrkbl = f2cmax(i__2,i__3);
1265 /* Computing MAX */
1266                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_q__;
1267                     wrkbl = f2cmax(i__2,i__3);
1268                     wrkbl = f2cmax(wrkbl,bdspac);
1269                     maxwrk = *m * *m + wrkbl;
1270 /* Computing MAX */
1271                     i__2 = *m * 3 + *n;
1272                     minwrk = f2cmax(i__2,bdspac);
1273                 } else if (wntva && wntun) {
1274
1275 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
1276
1277                     wrkbl = *m + lwork_dgelqf__;
1278 /* Computing MAX */
1279                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_n__;
1280                     wrkbl = f2cmax(i__2,i__3);
1281 /* Computing MAX */
1282                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1283                     wrkbl = f2cmax(i__2,i__3);
1284 /* Computing MAX */
1285                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1286                     wrkbl = f2cmax(i__2,i__3);
1287                     wrkbl = f2cmax(wrkbl,bdspac);
1288                     maxwrk = *m * *m + wrkbl;
1289 /* Computing MAX */
1290                     i__2 = *m * 3 + *n;
1291                     minwrk = f2cmax(i__2,bdspac);
1292                 } else if (wntva && wntuo) {
1293
1294 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
1295
1296                     wrkbl = *m + lwork_dgelqf__;
1297 /* Computing MAX */
1298                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_n__;
1299                     wrkbl = f2cmax(i__2,i__3);
1300 /* Computing MAX */
1301                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1302                     wrkbl = f2cmax(i__2,i__3);
1303 /* Computing MAX */
1304                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1305                     wrkbl = f2cmax(i__2,i__3);
1306 /* Computing MAX */
1307                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_q__;
1308                     wrkbl = f2cmax(i__2,i__3);
1309                     wrkbl = f2cmax(wrkbl,bdspac);
1310                     maxwrk = (*m << 1) * *m + wrkbl;
1311 /* Computing MAX */
1312                     i__2 = *m * 3 + *n;
1313                     minwrk = f2cmax(i__2,bdspac);
1314                 } else if (wntva && wntuas) {
1315
1316 /*                 Path 9t(N much larger than M, JOBU='S' or 'A', */
1317 /*                 JOBVT='A') */
1318
1319                     wrkbl = *m + lwork_dgelqf__;
1320 /* Computing MAX */
1321                     i__2 = wrkbl, i__3 = *m + lwork_dorglq_n__;
1322                     wrkbl = f2cmax(i__2,i__3);
1323 /* Computing MAX */
1324                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dgebrd__;
1325                     wrkbl = f2cmax(i__2,i__3);
1326 /* Computing MAX */
1327                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_p__;
1328                     wrkbl = f2cmax(i__2,i__3);
1329 /* Computing MAX */
1330                     i__2 = wrkbl, i__3 = *m * 3 + lwork_dorgbr_q__;
1331                     wrkbl = f2cmax(i__2,i__3);
1332                     wrkbl = f2cmax(wrkbl,bdspac);
1333                     maxwrk = *m * *m + wrkbl;
1334 /* Computing MAX */
1335                     i__2 = *m * 3 + *n;
1336                     minwrk = f2cmax(i__2,bdspac);
1337                 }
1338             } else {
1339
1340 /*              Path 10t(N greater than M, but not much larger) */
1341
1342                 dgebrd_(m, n, &a[a_offset], lda, &s[1], dum, dum, dum, dum, &
1343                         c_n1, &ierr);
1344                 lwork_dgebrd__ = (integer) dum[0];
1345                 maxwrk = *m * 3 + lwork_dgebrd__;
1346                 if (wntvs || wntvo) {
1347 /*                Compute space needed for DORGBR P */
1348                     dorgbr_("P", m, n, m, &a[a_offset], n, dum, dum, &c_n1, &
1349                             ierr);
1350                     lwork_dorgbr_p__ = (integer) dum[0];
1351 /* Computing MAX */
1352                     i__2 = maxwrk, i__3 = *m * 3 + lwork_dorgbr_p__;
1353                     maxwrk = f2cmax(i__2,i__3);
1354                 }
1355                 if (wntva) {
1356                     dorgbr_("P", n, n, m, &a[a_offset], n, dum, dum, &c_n1, &
1357                             ierr);
1358                     lwork_dorgbr_p__ = (integer) dum[0];
1359 /* Computing MAX */
1360                     i__2 = maxwrk, i__3 = *m * 3 + lwork_dorgbr_p__;
1361                     maxwrk = f2cmax(i__2,i__3);
1362                 }
1363                 if (! wntun) {
1364 /* Computing MAX */
1365                     i__2 = maxwrk, i__3 = *m * 3 + lwork_dorgbr_q__;
1366                     maxwrk = f2cmax(i__2,i__3);
1367                 }
1368                 maxwrk = f2cmax(maxwrk,bdspac);
1369 /* Computing MAX */
1370                 i__2 = *m * 3 + *n;
1371                 minwrk = f2cmax(i__2,bdspac);
1372             }
1373         }
1374         maxwrk = f2cmax(maxwrk,minwrk);
1375         work[1] = (doublereal) maxwrk;
1376
1377         if (*lwork < minwrk && ! lquery) {
1378             *info = -13;
1379         }
1380     }
1381
1382     if (*info != 0) {
1383         i__2 = -(*info);
1384         xerbla_("DGESVD", &i__2, (ftnlen)6);
1385         return 0;
1386     } else if (lquery) {
1387         return 0;
1388     }
1389
1390 /*     Quick return if possible */
1391
1392     if (*m == 0 || *n == 0) {
1393         return 0;
1394     }
1395
1396 /*     Get machine constants */
1397
1398     eps = dlamch_("P");
1399     smlnum = sqrt(dlamch_("S")) / eps;
1400     bignum = 1. / smlnum;
1401
1402 /*     Scale A if f2cmax element outside range [SMLNUM,BIGNUM] */
1403
1404     anrm = dlange_("M", m, n, &a[a_offset], lda, dum);
1405     iscl = 0;
1406     if (anrm > 0. && anrm < smlnum) {
1407         iscl = 1;
1408         dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
1409                 ierr);
1410     } else if (anrm > bignum) {
1411         iscl = 1;
1412         dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
1413                 ierr);
1414     }
1415
1416     if (*m >= *n) {
1417
1418 /*        A has at least as many rows as columns. If A has sufficiently */
1419 /*        more rows than columns, first reduce using the QR */
1420 /*        decomposition (if sufficient workspace available) */
1421
1422         if (*m >= mnthr) {
1423
1424             if (wntun) {
1425
1426 /*              Path 1 (M much larger than N, JOBU='N') */
1427 /*              No left singular vectors to be computed */
1428
1429                 itau = 1;
1430                 iwork = itau + *n;
1431
1432 /*              Compute A=Q*R */
1433 /*              (Workspace: need 2*N, prefer N + N*NB) */
1434
1435                 i__2 = *lwork - iwork + 1;
1436                 dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
1437                         i__2, &ierr);
1438
1439 /*              Zero out below R */
1440
1441                 if (*n > 1) {
1442                     i__2 = *n - 1;
1443                     i__3 = *n - 1;
1444                     dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &a[a_dim1 + 2],
1445                              lda);
1446                 }
1447                 ie = 1;
1448                 itauq = ie + *n;
1449                 itaup = itauq + *n;
1450                 iwork = itaup + *n;
1451
1452 /*              Bidiagonalize R in A */
1453 /*              (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
1454
1455                 i__2 = *lwork - iwork + 1;
1456                 dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
1457                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
1458                 ncvt = 0;
1459                 if (wntvo || wntvas) {
1460
1461 /*                 If right singular vectors desired, generate P'. */
1462 /*                 (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
1463
1464                     i__2 = *lwork - iwork + 1;
1465                     dorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &
1466                             work[iwork], &i__2, &ierr);
1467                     ncvt = *n;
1468                 }
1469                 iwork = ie + *n;
1470
1471 /*              Perform bidiagonal QR iteration, computing right */
1472 /*              singular vectors of A in A if desired */
1473 /*              (Workspace: need BDSPAC) */
1474
1475                 dbdsqr_("U", n, &ncvt, &c__0, &c__0, &s[1], &work[ie], &a[
1476                         a_offset], lda, dum, &c__1, dum, &c__1, &work[iwork],
1477                         info);
1478
1479 /*              If right singular vectors desired in VT, copy them there */
1480
1481                 if (wntvas) {
1482                     dlacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset],
1483                             ldvt);
1484                 }
1485
1486             } else if (wntuo && wntvn) {
1487
1488 /*              Path 2 (M much larger than N, JOBU='O', JOBVT='N') */
1489 /*              N left singular vectors to be overwritten on A and */
1490 /*              no right singular vectors to be computed */
1491
1492 /* Computing MAX */
1493                 i__2 = *n << 2;
1494                 if (*lwork >= *n * *n + f2cmax(i__2,bdspac)) {
1495
1496 /*                 Sufficient workspace for a fast algorithm */
1497
1498                     ir = 1;
1499 /* Computing MAX */
1500                     i__2 = wrkbl, i__3 = *lda * *n + *n;
1501                     if (*lwork >= f2cmax(i__2,i__3) + *lda * *n) {
1502
1503 /*                    WORK(IU) is LDA by N, WORK(IR) is LDA by N */
1504
1505                         ldwrku = *lda;
1506                         ldwrkr = *lda;
1507                     } else /* if(complicated condition) */ {
1508 /* Computing MAX */
1509                         i__2 = wrkbl, i__3 = *lda * *n + *n;
1510                         if (*lwork >= f2cmax(i__2,i__3) + *n * *n) {
1511
1512 /*                    WORK(IU) is LDA by N, WORK(IR) is N by N */
1513
1514                             ldwrku = *lda;
1515                             ldwrkr = *n;
1516                         } else {
1517
1518 /*                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N */
1519
1520                             ldwrku = (*lwork - *n * *n - *n) / *n;
1521                             ldwrkr = *n;
1522                         }
1523                     }
1524                     itau = ir + ldwrkr * *n;
1525                     iwork = itau + *n;
1526
1527 /*                 Compute A=Q*R */
1528 /*                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1529
1530                     i__2 = *lwork - iwork + 1;
1531                     dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1532                             , &i__2, &ierr);
1533
1534 /*                 Copy R to WORK(IR) and zero out below it */
1535
1536                     dlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1537                     i__2 = *n - 1;
1538                     i__3 = *n - 1;
1539                     dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[ir + 1],
1540                             &ldwrkr);
1541
1542 /*                 Generate Q in A */
1543 /*                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1544
1545                     i__2 = *lwork - iwork + 1;
1546                     dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1547                             iwork], &i__2, &ierr);
1548                     ie = itau;
1549                     itauq = ie + *n;
1550                     itaup = itauq + *n;
1551                     iwork = itaup + *n;
1552
1553 /*                 Bidiagonalize R in WORK(IR) */
1554 /*                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
1555
1556                     i__2 = *lwork - iwork + 1;
1557                     dgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
1558                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
1559
1560 /*                 Generate left vectors bidiagonalizing R */
1561 /*                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
1562
1563                     i__2 = *lwork - iwork + 1;
1564                     dorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1565                             work[iwork], &i__2, &ierr);
1566                     iwork = ie + *n;
1567
1568 /*                 Perform bidiagonal QR iteration, computing left */
1569 /*                 singular vectors of R in WORK(IR) */
1570 /*                 (Workspace: need N*N + BDSPAC) */
1571
1572                     dbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie], dum, &
1573                             c__1, &work[ir], &ldwrkr, dum, &c__1, &work[iwork]
1574                             , info);
1575                     iu = ie + *n;
1576
1577 /*                 Multiply Q in A by left singular vectors of R in */
1578 /*                 WORK(IR), storing result in WORK(IU) and copying to A */
1579 /*                 (Workspace: need N*N + 2*N, prefer N*N + M*N + N) */
1580
1581                     i__2 = *m;
1582                     i__3 = ldwrku;
1583                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1584                              i__3) {
1585 /* Computing MIN */
1586                         i__4 = *m - i__ + 1;
1587                         chunk = f2cmin(i__4,ldwrku);
1588                         dgemm_("N", "N", &chunk, n, n, &c_b79, &a[i__ +
1589                                 a_dim1], lda, &work[ir], &ldwrkr, &c_b57, &
1590                                 work[iu], &ldwrku);
1591                         dlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
1592                                 a_dim1], lda);
1593 /* L10: */
1594                     }
1595
1596                 } else {
1597
1598 /*                 Insufficient workspace for a fast algorithm */
1599
1600                     ie = 1;
1601                     itauq = ie + *n;
1602                     itaup = itauq + *n;
1603                     iwork = itaup + *n;
1604
1605 /*                 Bidiagonalize A */
1606 /*                 (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) */
1607
1608                     i__3 = *lwork - iwork + 1;
1609                     dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[
1610                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
1611
1612 /*                 Generate left vectors bidiagonalizing A */
1613 /*                 (Workspace: need 4*N, prefer 3*N + N*NB) */
1614
1615                     i__3 = *lwork - iwork + 1;
1616                     dorgbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
1617                             work[iwork], &i__3, &ierr);
1618                     iwork = ie + *n;
1619
1620 /*                 Perform bidiagonal QR iteration, computing left */
1621 /*                 singular vectors of A in A */
1622 /*                 (Workspace: need BDSPAC) */
1623
1624                     dbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie], dum, &
1625                             c__1, &a[a_offset], lda, dum, &c__1, &work[iwork],
1626                              info);
1627
1628                 }
1629
1630             } else if (wntuo && wntvas) {
1631
1632 /*              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') */
1633 /*              N left singular vectors to be overwritten on A and */
1634 /*              N right singular vectors to be computed in VT */
1635
1636 /* Computing MAX */
1637                 i__3 = *n << 2;
1638                 if (*lwork >= *n * *n + f2cmax(i__3,bdspac)) {
1639
1640 /*                 Sufficient workspace for a fast algorithm */
1641
1642                     ir = 1;
1643 /* Computing MAX */
1644                     i__3 = wrkbl, i__2 = *lda * *n + *n;
1645                     if (*lwork >= f2cmax(i__3,i__2) + *lda * *n) {
1646
1647 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1648
1649                         ldwrku = *lda;
1650                         ldwrkr = *lda;
1651                     } else /* if(complicated condition) */ {
1652 /* Computing MAX */
1653                         i__3 = wrkbl, i__2 = *lda * *n + *n;
1654                         if (*lwork >= f2cmax(i__3,i__2) + *n * *n) {
1655
1656 /*                    WORK(IU) is LDA by N and WORK(IR) is N by N */
1657
1658                             ldwrku = *lda;
1659                             ldwrkr = *n;
1660                         } else {
1661
1662 /*                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N */
1663
1664                             ldwrku = (*lwork - *n * *n - *n) / *n;
1665                             ldwrkr = *n;
1666                         }
1667                     }
1668                     itau = ir + ldwrkr * *n;
1669                     iwork = itau + *n;
1670
1671 /*                 Compute A=Q*R */
1672 /*                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1673
1674                     i__3 = *lwork - iwork + 1;
1675                     dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1676                             , &i__3, &ierr);
1677
1678 /*                 Copy R to VT, zeroing out below it */
1679
1680                     dlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1681                             ldvt);
1682                     if (*n > 1) {
1683                         i__3 = *n - 1;
1684                         i__2 = *n - 1;
1685                         dlaset_("L", &i__3, &i__2, &c_b57, &c_b57, &vt[
1686                                 vt_dim1 + 2], ldvt);
1687                     }
1688
1689 /*                 Generate Q in A */
1690 /*                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1691
1692                     i__3 = *lwork - iwork + 1;
1693                     dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1694                             iwork], &i__3, &ierr);
1695                     ie = itau;
1696                     itauq = ie + *n;
1697                     itaup = itauq + *n;
1698                     iwork = itaup + *n;
1699
1700 /*                 Bidiagonalize R in VT, copying result to WORK(IR) */
1701 /*                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
1702
1703                     i__3 = *lwork - iwork + 1;
1704                     dgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &
1705                             work[itauq], &work[itaup], &work[iwork], &i__3, &
1706                             ierr);
1707                     dlacpy_("L", n, n, &vt[vt_offset], ldvt, &work[ir], &
1708                             ldwrkr);
1709
1710 /*                 Generate left vectors bidiagonalizing R in WORK(IR) */
1711 /*                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
1712
1713                     i__3 = *lwork - iwork + 1;
1714                     dorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &
1715                             work[iwork], &i__3, &ierr);
1716
1717 /*                 Generate right vectors bidiagonalizing R in VT */
1718 /*                 (Workspace: need N*N + 4*N-1, prefer N*N + 3*N + (N-1)*NB) */
1719
1720                     i__3 = *lwork - iwork + 1;
1721                     dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1722                             &work[iwork], &i__3, &ierr);
1723                     iwork = ie + *n;
1724
1725 /*                 Perform bidiagonal QR iteration, computing left */
1726 /*                 singular vectors of R in WORK(IR) and computing right */
1727 /*                 singular vectors of R in VT */
1728 /*                 (Workspace: need N*N + BDSPAC) */
1729
1730                     dbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
1731                             vt_offset], ldvt, &work[ir], &ldwrkr, dum, &c__1,
1732                             &work[iwork], info);
1733                     iu = ie + *n;
1734
1735 /*                 Multiply Q in A by left singular vectors of R in */
1736 /*                 WORK(IR), storing result in WORK(IU) and copying to A */
1737 /*                 (Workspace: need N*N + 2*N, prefer N*N + M*N + N) */
1738
1739                     i__3 = *m;
1740                     i__2 = ldwrku;
1741                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
1742                              i__2) {
1743 /* Computing MIN */
1744                         i__4 = *m - i__ + 1;
1745                         chunk = f2cmin(i__4,ldwrku);
1746                         dgemm_("N", "N", &chunk, n, n, &c_b79, &a[i__ +
1747                                 a_dim1], lda, &work[ir], &ldwrkr, &c_b57, &
1748                                 work[iu], &ldwrku);
1749                         dlacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
1750                                 a_dim1], lda);
1751 /* L20: */
1752                     }
1753
1754                 } else {
1755
1756 /*                 Insufficient workspace for a fast algorithm */
1757
1758                     itau = 1;
1759                     iwork = itau + *n;
1760
1761 /*                 Compute A=Q*R */
1762 /*                 (Workspace: need 2*N, prefer N + N*NB) */
1763
1764                     i__2 = *lwork - iwork + 1;
1765                     dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
1766                             , &i__2, &ierr);
1767
1768 /*                 Copy R to VT, zeroing out below it */
1769
1770                     dlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
1771                             ldvt);
1772                     if (*n > 1) {
1773                         i__2 = *n - 1;
1774                         i__3 = *n - 1;
1775                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &vt[
1776                                 vt_dim1 + 2], ldvt);
1777                     }
1778
1779 /*                 Generate Q in A */
1780 /*                 (Workspace: need 2*N, prefer N + N*NB) */
1781
1782                     i__2 = *lwork - iwork + 1;
1783                     dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[
1784                             iwork], &i__2, &ierr);
1785                     ie = itau;
1786                     itauq = ie + *n;
1787                     itaup = itauq + *n;
1788                     iwork = itaup + *n;
1789
1790 /*                 Bidiagonalize R in VT */
1791 /*                 (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
1792
1793                     i__2 = *lwork - iwork + 1;
1794                     dgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &
1795                             work[itauq], &work[itaup], &work[iwork], &i__2, &
1796                             ierr);
1797
1798 /*                 Multiply Q in A by left vectors bidiagonalizing R */
1799 /*                 (Workspace: need 3*N + M, prefer 3*N + M*NB) */
1800
1801                     i__2 = *lwork - iwork + 1;
1802                     dormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &
1803                             work[itauq], &a[a_offset], lda, &work[iwork], &
1804                             i__2, &ierr);
1805
1806 /*                 Generate right vectors bidiagonalizing R in VT */
1807 /*                 (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
1808
1809                     i__2 = *lwork - iwork + 1;
1810                     dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup],
1811                             &work[iwork], &i__2, &ierr);
1812                     iwork = ie + *n;
1813
1814 /*                 Perform bidiagonal QR iteration, computing left */
1815 /*                 singular vectors of A in A and computing right */
1816 /*                 singular vectors of A in VT */
1817 /*                 (Workspace: need BDSPAC) */
1818
1819                     dbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
1820                             vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
1821                             work[iwork], info);
1822
1823                 }
1824
1825             } else if (wntus) {
1826
1827                 if (wntvn) {
1828
1829 /*                 Path 4 (M much larger than N, JOBU='S', JOBVT='N') */
1830 /*                 N left singular vectors to be computed in U and */
1831 /*                 no right singular vectors to be computed */
1832
1833 /* Computing MAX */
1834                     i__2 = *n << 2;
1835                     if (*lwork >= *n * *n + f2cmax(i__2,bdspac)) {
1836
1837 /*                    Sufficient workspace for a fast algorithm */
1838
1839                         ir = 1;
1840                         if (*lwork >= wrkbl + *lda * *n) {
1841
1842 /*                       WORK(IR) is LDA by N */
1843
1844                             ldwrkr = *lda;
1845                         } else {
1846
1847 /*                       WORK(IR) is N by N */
1848
1849                             ldwrkr = *n;
1850                         }
1851                         itau = ir + ldwrkr * *n;
1852                         iwork = itau + *n;
1853
1854 /*                    Compute A=Q*R */
1855 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1856
1857                         i__2 = *lwork - iwork + 1;
1858                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1859                                 iwork], &i__2, &ierr);
1860
1861 /*                    Copy R to WORK(IR), zeroing out below it */
1862
1863                         dlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
1864                                 ldwrkr);
1865                         i__2 = *n - 1;
1866                         i__3 = *n - 1;
1867                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[ir +
1868                                 1], &ldwrkr);
1869
1870 /*                    Generate Q in A */
1871 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
1872
1873                         i__2 = *lwork - iwork + 1;
1874                         dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
1875                                 work[iwork], &i__2, &ierr);
1876                         ie = itau;
1877                         itauq = ie + *n;
1878                         itaup = itauq + *n;
1879                         iwork = itaup + *n;
1880
1881 /*                    Bidiagonalize R in WORK(IR) */
1882 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
1883
1884                         i__2 = *lwork - iwork + 1;
1885                         dgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &
1886                                 work[itauq], &work[itaup], &work[iwork], &
1887                                 i__2, &ierr);
1888
1889 /*                    Generate left vectors bidiagonalizing R in WORK(IR) */
1890 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
1891
1892                         i__2 = *lwork - iwork + 1;
1893                         dorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
1894                                 , &work[iwork], &i__2, &ierr);
1895                         iwork = ie + *n;
1896
1897 /*                    Perform bidiagonal QR iteration, computing left */
1898 /*                    singular vectors of R in WORK(IR) */
1899 /*                    (Workspace: need N*N + BDSPAC) */
1900
1901                         dbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie],
1902                                 dum, &c__1, &work[ir], &ldwrkr, dum, &c__1, &
1903                                 work[iwork], info);
1904
1905 /*                    Multiply Q in A by left singular vectors of R in */
1906 /*                    WORK(IR), storing result in U */
1907 /*                    (Workspace: need N*N) */
1908
1909                         dgemm_("N", "N", m, n, n, &c_b79, &a[a_offset], lda, &
1910                                 work[ir], &ldwrkr, &c_b57, &u[u_offset], ldu);
1911
1912                     } else {
1913
1914 /*                    Insufficient workspace for a fast algorithm */
1915
1916                         itau = 1;
1917                         iwork = itau + *n;
1918
1919 /*                    Compute A=Q*R, copying result to U */
1920 /*                    (Workspace: need 2*N, prefer N + N*NB) */
1921
1922                         i__2 = *lwork - iwork + 1;
1923                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
1924                                 iwork], &i__2, &ierr);
1925                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
1926                                 ldu);
1927
1928 /*                    Generate Q in U */
1929 /*                    (Workspace: need 2*N, prefer N + N*NB) */
1930
1931                         i__2 = *lwork - iwork + 1;
1932                         dorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
1933                                 work[iwork], &i__2, &ierr);
1934                         ie = itau;
1935                         itauq = ie + *n;
1936                         itaup = itauq + *n;
1937                         iwork = itaup + *n;
1938
1939 /*                    Zero out below R in A */
1940
1941                         if (*n > 1) {
1942                             i__2 = *n - 1;
1943                             i__3 = *n - 1;
1944                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &a[
1945                                     a_dim1 + 2], lda);
1946                         }
1947
1948 /*                    Bidiagonalize R in A */
1949 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
1950
1951                         i__2 = *lwork - iwork + 1;
1952                         dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
1953                                 work[itauq], &work[itaup], &work[iwork], &
1954                                 i__2, &ierr);
1955
1956 /*                    Multiply Q in U by left vectors bidiagonalizing R */
1957 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
1958
1959                         i__2 = *lwork - iwork + 1;
1960                         dormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
1961                                 work[itauq], &u[u_offset], ldu, &work[iwork],
1962                                 &i__2, &ierr)
1963                                 ;
1964                         iwork = ie + *n;
1965
1966 /*                    Perform bidiagonal QR iteration, computing left */
1967 /*                    singular vectors of A in U */
1968 /*                    (Workspace: need BDSPAC) */
1969
1970                         dbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie],
1971                                 dum, &c__1, &u[u_offset], ldu, dum, &c__1, &
1972                                 work[iwork], info);
1973
1974                     }
1975
1976                 } else if (wntvo) {
1977
1978 /*                 Path 5 (M much larger than N, JOBU='S', JOBVT='O') */
1979 /*                 N left singular vectors to be computed in U and */
1980 /*                 N right singular vectors to be overwritten on A */
1981
1982 /* Computing MAX */
1983                     i__2 = *n << 2;
1984                     if (*lwork >= (*n << 1) * *n + f2cmax(i__2,bdspac)) {
1985
1986 /*                    Sufficient workspace for a fast algorithm */
1987
1988                         iu = 1;
1989                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
1990
1991 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
1992
1993                             ldwrku = *lda;
1994                             ir = iu + ldwrku * *n;
1995                             ldwrkr = *lda;
1996                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
1997
1998 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
1999
2000                             ldwrku = *lda;
2001                             ir = iu + ldwrku * *n;
2002                             ldwrkr = *n;
2003                         } else {
2004
2005 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
2006
2007                             ldwrku = *n;
2008                             ir = iu + ldwrku * *n;
2009                             ldwrkr = *n;
2010                         }
2011                         itau = ir + ldwrkr * *n;
2012                         iwork = itau + *n;
2013
2014 /*                    Compute A=Q*R */
2015 /*                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) */
2016
2017                         i__2 = *lwork - iwork + 1;
2018                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2019                                 iwork], &i__2, &ierr);
2020
2021 /*                    Copy R to WORK(IU), zeroing out below it */
2022
2023                         dlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2024                                 ldwrku);
2025                         i__2 = *n - 1;
2026                         i__3 = *n - 1;
2027                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
2028                                 1], &ldwrku);
2029
2030 /*                    Generate Q in A */
2031 /*                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) */
2032
2033                         i__2 = *lwork - iwork + 1;
2034                         dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
2035                                 work[iwork], &i__2, &ierr);
2036                         ie = itau;
2037                         itauq = ie + *n;
2038                         itaup = itauq + *n;
2039                         iwork = itaup + *n;
2040
2041 /*                    Bidiagonalize R in WORK(IU), copying result to */
2042 /*                    WORK(IR) */
2043 /*                    (Workspace: need 2*N*N + 4*N, */
2044 /*                                prefer 2*N*N+3*N+2*N*NB) */
2045
2046                         i__2 = *lwork - iwork + 1;
2047                         dgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
2048                                 work[itauq], &work[itaup], &work[iwork], &
2049                                 i__2, &ierr);
2050                         dlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
2051                                 ldwrkr);
2052
2053 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
2054 /*                    (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) */
2055
2056                         i__2 = *lwork - iwork + 1;
2057                         dorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2058                                 , &work[iwork], &i__2, &ierr);
2059
2060 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
2061 /*                    (Workspace: need 2*N*N + 4*N-1, */
2062 /*                                prefer 2*N*N+3*N+(N-1)*NB) */
2063
2064                         i__2 = *lwork - iwork + 1;
2065                         dorgbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
2066                                 , &work[iwork], &i__2, &ierr);
2067                         iwork = ie + *n;
2068
2069 /*                    Perform bidiagonal QR iteration, computing left */
2070 /*                    singular vectors of R in WORK(IU) and computing */
2071 /*                    right singular vectors of R in WORK(IR) */
2072 /*                    (Workspace: need 2*N*N + BDSPAC) */
2073
2074                         dbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &work[
2075                                 ir], &ldwrkr, &work[iu], &ldwrku, dum, &c__1,
2076                                 &work[iwork], info);
2077
2078 /*                    Multiply Q in A by left singular vectors of R in */
2079 /*                    WORK(IU), storing result in U */
2080 /*                    (Workspace: need N*N) */
2081
2082                         dgemm_("N", "N", m, n, n, &c_b79, &a[a_offset], lda, &
2083                                 work[iu], &ldwrku, &c_b57, &u[u_offset], ldu);
2084
2085 /*                    Copy right singular vectors of R to A */
2086 /*                    (Workspace: need N*N) */
2087
2088                         dlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
2089                                 lda);
2090
2091                     } else {
2092
2093 /*                    Insufficient workspace for a fast algorithm */
2094
2095                         itau = 1;
2096                         iwork = itau + *n;
2097
2098 /*                    Compute A=Q*R, copying result to U */
2099 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2100
2101                         i__2 = *lwork - iwork + 1;
2102                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2103                                 iwork], &i__2, &ierr);
2104                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2105                                 ldu);
2106
2107 /*                    Generate Q in U */
2108 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2109
2110                         i__2 = *lwork - iwork + 1;
2111                         dorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
2112                                 work[iwork], &i__2, &ierr);
2113                         ie = itau;
2114                         itauq = ie + *n;
2115                         itaup = itauq + *n;
2116                         iwork = itaup + *n;
2117
2118 /*                    Zero out below R in A */
2119
2120                         if (*n > 1) {
2121                             i__2 = *n - 1;
2122                             i__3 = *n - 1;
2123                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &a[
2124                                     a_dim1 + 2], lda);
2125                         }
2126
2127 /*                    Bidiagonalize R in A */
2128 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
2129
2130                         i__2 = *lwork - iwork + 1;
2131                         dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
2132                                 work[itauq], &work[itaup], &work[iwork], &
2133                                 i__2, &ierr);
2134
2135 /*                    Multiply Q in U by left vectors bidiagonalizing R */
2136 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
2137
2138                         i__2 = *lwork - iwork + 1;
2139                         dormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
2140                                 work[itauq], &u[u_offset], ldu, &work[iwork],
2141                                 &i__2, &ierr)
2142                                 ;
2143
2144 /*                    Generate right vectors bidiagonalizing R in A */
2145 /*                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2146
2147                         i__2 = *lwork - iwork + 1;
2148                         dorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
2149                                  &work[iwork], &i__2, &ierr);
2150                         iwork = ie + *n;
2151
2152 /*                    Perform bidiagonal QR iteration, computing left */
2153 /*                    singular vectors of A in U and computing right */
2154 /*                    singular vectors of A in A */
2155 /*                    (Workspace: need BDSPAC) */
2156
2157                         dbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &a[
2158                                 a_offset], lda, &u[u_offset], ldu, dum, &c__1,
2159                                  &work[iwork], info);
2160
2161                     }
2162
2163                 } else if (wntvas) {
2164
2165 /*                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' */
2166 /*                         or 'A') */
2167 /*                 N left singular vectors to be computed in U and */
2168 /*                 N right singular vectors to be computed in VT */
2169
2170 /* Computing MAX */
2171                     i__2 = *n << 2;
2172                     if (*lwork >= *n * *n + f2cmax(i__2,bdspac)) {
2173
2174 /*                    Sufficient workspace for a fast algorithm */
2175
2176                         iu = 1;
2177                         if (*lwork >= wrkbl + *lda * *n) {
2178
2179 /*                       WORK(IU) is LDA by N */
2180
2181                             ldwrku = *lda;
2182                         } else {
2183
2184 /*                       WORK(IU) is N by N */
2185
2186                             ldwrku = *n;
2187                         }
2188                         itau = iu + ldwrku * *n;
2189                         iwork = itau + *n;
2190
2191 /*                    Compute A=Q*R */
2192 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
2193
2194                         i__2 = *lwork - iwork + 1;
2195                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2196                                 iwork], &i__2, &ierr);
2197
2198 /*                    Copy R to WORK(IU), zeroing out below it */
2199
2200                         dlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2201                                 ldwrku);
2202                         i__2 = *n - 1;
2203                         i__3 = *n - 1;
2204                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
2205                                 1], &ldwrku);
2206
2207 /*                    Generate Q in A */
2208 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
2209
2210                         i__2 = *lwork - iwork + 1;
2211                         dorgqr_(m, n, n, &a[a_offset], lda, &work[itau], &
2212                                 work[iwork], &i__2, &ierr);
2213                         ie = itau;
2214                         itauq = ie + *n;
2215                         itaup = itauq + *n;
2216                         iwork = itaup + *n;
2217
2218 /*                    Bidiagonalize R in WORK(IU), copying result to VT */
2219 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
2220
2221                         i__2 = *lwork - iwork + 1;
2222                         dgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
2223                                 work[itauq], &work[itaup], &work[iwork], &
2224                                 i__2, &ierr);
2225                         dlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
2226                                  ldvt);
2227
2228 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
2229 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
2230
2231                         i__2 = *lwork - iwork + 1;
2232                         dorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2233                                 , &work[iwork], &i__2, &ierr);
2234
2235 /*                    Generate right bidiagonalizing vectors in VT */
2236 /*                    (Workspace: need N*N + 4*N-1, */
2237 /*                                prefer N*N+3*N+(N-1)*NB) */
2238
2239                         i__2 = *lwork - iwork + 1;
2240                         dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2241                                 itaup], &work[iwork], &i__2, &ierr)
2242                                 ;
2243                         iwork = ie + *n;
2244
2245 /*                    Perform bidiagonal QR iteration, computing left */
2246 /*                    singular vectors of R in WORK(IU) and computing */
2247 /*                    right singular vectors of R in VT */
2248 /*                    (Workspace: need N*N + BDSPAC) */
2249
2250                         dbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
2251                                 vt_offset], ldvt, &work[iu], &ldwrku, dum, &
2252                                 c__1, &work[iwork], info);
2253
2254 /*                    Multiply Q in A by left singular vectors of R in */
2255 /*                    WORK(IU), storing result in U */
2256 /*                    (Workspace: need N*N) */
2257
2258                         dgemm_("N", "N", m, n, n, &c_b79, &a[a_offset], lda, &
2259                                 work[iu], &ldwrku, &c_b57, &u[u_offset], ldu);
2260
2261                     } else {
2262
2263 /*                    Insufficient workspace for a fast algorithm */
2264
2265                         itau = 1;
2266                         iwork = itau + *n;
2267
2268 /*                    Compute A=Q*R, copying result to U */
2269 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2270
2271                         i__2 = *lwork - iwork + 1;
2272                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2273                                 iwork], &i__2, &ierr);
2274                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2275                                 ldu);
2276
2277 /*                    Generate Q in U */
2278 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2279
2280                         i__2 = *lwork - iwork + 1;
2281                         dorgqr_(m, n, n, &u[u_offset], ldu, &work[itau], &
2282                                 work[iwork], &i__2, &ierr);
2283
2284 /*                    Copy R to VT, zeroing out below it */
2285
2286                         dlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
2287                                 ldvt);
2288                         if (*n > 1) {
2289                             i__2 = *n - 1;
2290                             i__3 = *n - 1;
2291                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &vt[
2292                                     vt_dim1 + 2], ldvt);
2293                         }
2294                         ie = itau;
2295                         itauq = ie + *n;
2296                         itaup = itauq + *n;
2297                         iwork = itaup + *n;
2298
2299 /*                    Bidiagonalize R in VT */
2300 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
2301
2302                         i__2 = *lwork - iwork + 1;
2303                         dgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie],
2304                                 &work[itauq], &work[itaup], &work[iwork], &
2305                                 i__2, &ierr);
2306
2307 /*                    Multiply Q in U by left bidiagonalizing vectors */
2308 /*                    in VT */
2309 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
2310
2311                         i__2 = *lwork - iwork + 1;
2312                         dormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
2313                                 &work[itauq], &u[u_offset], ldu, &work[iwork],
2314                                  &i__2, &ierr);
2315
2316 /*                    Generate right bidiagonalizing vectors in VT */
2317 /*                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2318
2319                         i__2 = *lwork - iwork + 1;
2320                         dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2321                                 itaup], &work[iwork], &i__2, &ierr)
2322                                 ;
2323                         iwork = ie + *n;
2324
2325 /*                    Perform bidiagonal QR iteration, computing left */
2326 /*                    singular vectors of A in U and computing right */
2327 /*                    singular vectors of A in VT */
2328 /*                    (Workspace: need BDSPAC) */
2329
2330                         dbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
2331                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
2332                                 c__1, &work[iwork], info);
2333
2334                     }
2335
2336                 }
2337
2338             } else if (wntua) {
2339
2340                 if (wntvn) {
2341
2342 /*                 Path 7 (M much larger than N, JOBU='A', JOBVT='N') */
2343 /*                 M left singular vectors to be computed in U and */
2344 /*                 no right singular vectors to be computed */
2345
2346 /* Computing MAX */
2347                     i__2 = *n + *m, i__3 = *n << 2, i__2 = f2cmax(i__2,i__3);
2348                     if (*lwork >= *n * *n + f2cmax(i__2,bdspac)) {
2349
2350 /*                    Sufficient workspace for a fast algorithm */
2351
2352                         ir = 1;
2353                         if (*lwork >= wrkbl + *lda * *n) {
2354
2355 /*                       WORK(IR) is LDA by N */
2356
2357                             ldwrkr = *lda;
2358                         } else {
2359
2360 /*                       WORK(IR) is N by N */
2361
2362                             ldwrkr = *n;
2363                         }
2364                         itau = ir + ldwrkr * *n;
2365                         iwork = itau + *n;
2366
2367 /*                    Compute A=Q*R, copying result to U */
2368 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
2369
2370                         i__2 = *lwork - iwork + 1;
2371                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2372                                 iwork], &i__2, &ierr);
2373                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2374                                 ldu);
2375
2376 /*                    Copy R to WORK(IR), zeroing out below it */
2377
2378                         dlacpy_("U", n, n, &a[a_offset], lda, &work[ir], &
2379                                 ldwrkr);
2380                         i__2 = *n - 1;
2381                         i__3 = *n - 1;
2382                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[ir +
2383                                 1], &ldwrkr);
2384
2385 /*                    Generate Q in U */
2386 /*                    (Workspace: need N*N + N + M, prefer N*N + N + M*NB) */
2387
2388                         i__2 = *lwork - iwork + 1;
2389                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2390                                 work[iwork], &i__2, &ierr);
2391                         ie = itau;
2392                         itauq = ie + *n;
2393                         itaup = itauq + *n;
2394                         iwork = itaup + *n;
2395
2396 /*                    Bidiagonalize R in WORK(IR) */
2397 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
2398
2399                         i__2 = *lwork - iwork + 1;
2400                         dgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &
2401                                 work[itauq], &work[itaup], &work[iwork], &
2402                                 i__2, &ierr);
2403
2404 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
2405 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
2406
2407                         i__2 = *lwork - iwork + 1;
2408                         dorgbr_("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq]
2409                                 , &work[iwork], &i__2, &ierr);
2410                         iwork = ie + *n;
2411
2412 /*                    Perform bidiagonal QR iteration, computing left */
2413 /*                    singular vectors of R in WORK(IR) */
2414 /*                    (Workspace: need N*N + BDSPAC) */
2415
2416                         dbdsqr_("U", n, &c__0, n, &c__0, &s[1], &work[ie],
2417                                 dum, &c__1, &work[ir], &ldwrkr, dum, &c__1, &
2418                                 work[iwork], info);
2419
2420 /*                    Multiply Q in U by left singular vectors of R in */
2421 /*                    WORK(IR), storing result in A */
2422 /*                    (Workspace: need N*N) */
2423
2424                         dgemm_("N", "N", m, n, n, &c_b79, &u[u_offset], ldu, &
2425                                 work[ir], &ldwrkr, &c_b57, &a[a_offset], lda);
2426
2427 /*                    Copy left singular vectors of A from A to U */
2428
2429                         dlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2430                                 ldu);
2431
2432                     } else {
2433
2434 /*                    Insufficient workspace for a fast algorithm */
2435
2436                         itau = 1;
2437                         iwork = itau + *n;
2438
2439 /*                    Compute A=Q*R, copying result to U */
2440 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2441
2442                         i__2 = *lwork - iwork + 1;
2443                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2444                                 iwork], &i__2, &ierr);
2445                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2446                                 ldu);
2447
2448 /*                    Generate Q in U */
2449 /*                    (Workspace: need N + M, prefer N + M*NB) */
2450
2451                         i__2 = *lwork - iwork + 1;
2452                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2453                                 work[iwork], &i__2, &ierr);
2454                         ie = itau;
2455                         itauq = ie + *n;
2456                         itaup = itauq + *n;
2457                         iwork = itaup + *n;
2458
2459 /*                    Zero out below R in A */
2460
2461                         if (*n > 1) {
2462                             i__2 = *n - 1;
2463                             i__3 = *n - 1;
2464                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &a[
2465                                     a_dim1 + 2], lda);
2466                         }
2467
2468 /*                    Bidiagonalize R in A */
2469 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
2470
2471                         i__2 = *lwork - iwork + 1;
2472                         dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
2473                                 work[itauq], &work[itaup], &work[iwork], &
2474                                 i__2, &ierr);
2475
2476 /*                    Multiply Q in U by left bidiagonalizing vectors */
2477 /*                    in A */
2478 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
2479
2480                         i__2 = *lwork - iwork + 1;
2481                         dormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
2482                                 work[itauq], &u[u_offset], ldu, &work[iwork],
2483                                 &i__2, &ierr)
2484                                 ;
2485                         iwork = ie + *n;
2486
2487 /*                    Perform bidiagonal QR iteration, computing left */
2488 /*                    singular vectors of A in U */
2489 /*                    (Workspace: need BDSPAC) */
2490
2491                         dbdsqr_("U", n, &c__0, m, &c__0, &s[1], &work[ie],
2492                                 dum, &c__1, &u[u_offset], ldu, dum, &c__1, &
2493                                 work[iwork], info);
2494
2495                     }
2496
2497                 } else if (wntvo) {
2498
2499 /*                 Path 8 (M much larger than N, JOBU='A', JOBVT='O') */
2500 /*                 M left singular vectors to be computed in U and */
2501 /*                 N right singular vectors to be overwritten on A */
2502
2503 /* Computing MAX */
2504                     i__2 = *n + *m, i__3 = *n << 2, i__2 = f2cmax(i__2,i__3);
2505                     if (*lwork >= (*n << 1) * *n + f2cmax(i__2,bdspac)) {
2506
2507 /*                    Sufficient workspace for a fast algorithm */
2508
2509                         iu = 1;
2510                         if (*lwork >= wrkbl + (*lda << 1) * *n) {
2511
2512 /*                       WORK(IU) is LDA by N and WORK(IR) is LDA by N */
2513
2514                             ldwrku = *lda;
2515                             ir = iu + ldwrku * *n;
2516                             ldwrkr = *lda;
2517                         } else if (*lwork >= wrkbl + (*lda + *n) * *n) {
2518
2519 /*                       WORK(IU) is LDA by N and WORK(IR) is N by N */
2520
2521                             ldwrku = *lda;
2522                             ir = iu + ldwrku * *n;
2523                             ldwrkr = *n;
2524                         } else {
2525
2526 /*                       WORK(IU) is N by N and WORK(IR) is N by N */
2527
2528                             ldwrku = *n;
2529                             ir = iu + ldwrku * *n;
2530                             ldwrkr = *n;
2531                         }
2532                         itau = ir + ldwrkr * *n;
2533                         iwork = itau + *n;
2534
2535 /*                    Compute A=Q*R, copying result to U */
2536 /*                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) */
2537
2538                         i__2 = *lwork - iwork + 1;
2539                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2540                                 iwork], &i__2, &ierr);
2541                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2542                                 ldu);
2543
2544 /*                    Generate Q in U */
2545 /*                    (Workspace: need 2*N*N + N + M, prefer 2*N*N + N + M*NB) */
2546
2547                         i__2 = *lwork - iwork + 1;
2548                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2549                                 work[iwork], &i__2, &ierr);
2550
2551 /*                    Copy R to WORK(IU), zeroing out below it */
2552
2553                         dlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2554                                 ldwrku);
2555                         i__2 = *n - 1;
2556                         i__3 = *n - 1;
2557                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
2558                                 1], &ldwrku);
2559                         ie = itau;
2560                         itauq = ie + *n;
2561                         itaup = itauq + *n;
2562                         iwork = itaup + *n;
2563
2564 /*                    Bidiagonalize R in WORK(IU), copying result to */
2565 /*                    WORK(IR) */
2566 /*                    (Workspace: need 2*N*N + 4*N, */
2567 /*                                prefer 2*N*N+3*N+2*N*NB) */
2568
2569                         i__2 = *lwork - iwork + 1;
2570                         dgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
2571                                 work[itauq], &work[itaup], &work[iwork], &
2572                                 i__2, &ierr);
2573                         dlacpy_("U", n, n, &work[iu], &ldwrku, &work[ir], &
2574                                 ldwrkr);
2575
2576 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
2577 /*                    (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) */
2578
2579                         i__2 = *lwork - iwork + 1;
2580                         dorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2581                                 , &work[iwork], &i__2, &ierr);
2582
2583 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
2584 /*                    (Workspace: need 2*N*N + 4*N-1, */
2585 /*                                prefer 2*N*N+3*N+(N-1)*NB) */
2586
2587                         i__2 = *lwork - iwork + 1;
2588                         dorgbr_("P", n, n, n, &work[ir], &ldwrkr, &work[itaup]
2589                                 , &work[iwork], &i__2, &ierr);
2590                         iwork = ie + *n;
2591
2592 /*                    Perform bidiagonal QR iteration, computing left */
2593 /*                    singular vectors of R in WORK(IU) and computing */
2594 /*                    right singular vectors of R in WORK(IR) */
2595 /*                    (Workspace: need 2*N*N + BDSPAC) */
2596
2597                         dbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &work[
2598                                 ir], &ldwrkr, &work[iu], &ldwrku, dum, &c__1,
2599                                 &work[iwork], info);
2600
2601 /*                    Multiply Q in U by left singular vectors of R in */
2602 /*                    WORK(IU), storing result in A */
2603 /*                    (Workspace: need N*N) */
2604
2605                         dgemm_("N", "N", m, n, n, &c_b79, &u[u_offset], ldu, &
2606                                 work[iu], &ldwrku, &c_b57, &a[a_offset], lda);
2607
2608 /*                    Copy left singular vectors of A from A to U */
2609
2610                         dlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2611                                 ldu);
2612
2613 /*                    Copy right singular vectors of R from WORK(IR) to A */
2614
2615                         dlacpy_("F", n, n, &work[ir], &ldwrkr, &a[a_offset],
2616                                 lda);
2617
2618                     } else {
2619
2620 /*                    Insufficient workspace for a fast algorithm */
2621
2622                         itau = 1;
2623                         iwork = itau + *n;
2624
2625 /*                    Compute A=Q*R, copying result to U */
2626 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2627
2628                         i__2 = *lwork - iwork + 1;
2629                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2630                                 iwork], &i__2, &ierr);
2631                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2632                                 ldu);
2633
2634 /*                    Generate Q in U */
2635 /*                    (Workspace: need N + M, prefer N + M*NB) */
2636
2637                         i__2 = *lwork - iwork + 1;
2638                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2639                                 work[iwork], &i__2, &ierr);
2640                         ie = itau;
2641                         itauq = ie + *n;
2642                         itaup = itauq + *n;
2643                         iwork = itaup + *n;
2644
2645 /*                    Zero out below R in A */
2646
2647                         if (*n > 1) {
2648                             i__2 = *n - 1;
2649                             i__3 = *n - 1;
2650                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &a[
2651                                     a_dim1 + 2], lda);
2652                         }
2653
2654 /*                    Bidiagonalize R in A */
2655 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
2656
2657                         i__2 = *lwork - iwork + 1;
2658                         dgebrd_(n, n, &a[a_offset], lda, &s[1], &work[ie], &
2659                                 work[itauq], &work[itaup], &work[iwork], &
2660                                 i__2, &ierr);
2661
2662 /*                    Multiply Q in U by left bidiagonalizing vectors */
2663 /*                    in A */
2664 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
2665
2666                         i__2 = *lwork - iwork + 1;
2667                         dormbr_("Q", "R", "N", m, n, n, &a[a_offset], lda, &
2668                                 work[itauq], &u[u_offset], ldu, &work[iwork],
2669                                 &i__2, &ierr)
2670                                 ;
2671
2672 /*                    Generate right bidiagonalizing vectors in A */
2673 /*                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2674
2675                         i__2 = *lwork - iwork + 1;
2676                         dorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup],
2677                                  &work[iwork], &i__2, &ierr);
2678                         iwork = ie + *n;
2679
2680 /*                    Perform bidiagonal QR iteration, computing left */
2681 /*                    singular vectors of A in U and computing right */
2682 /*                    singular vectors of A in A */
2683 /*                    (Workspace: need BDSPAC) */
2684
2685                         dbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &a[
2686                                 a_offset], lda, &u[u_offset], ldu, dum, &c__1,
2687                                  &work[iwork], info);
2688
2689                     }
2690
2691                 } else if (wntvas) {
2692
2693 /*                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' */
2694 /*                         or 'A') */
2695 /*                 M left singular vectors to be computed in U and */
2696 /*                 N right singular vectors to be computed in VT */
2697
2698 /* Computing MAX */
2699                     i__2 = *n + *m, i__3 = *n << 2, i__2 = f2cmax(i__2,i__3);
2700                     if (*lwork >= *n * *n + f2cmax(i__2,bdspac)) {
2701
2702 /*                    Sufficient workspace for a fast algorithm */
2703
2704                         iu = 1;
2705                         if (*lwork >= wrkbl + *lda * *n) {
2706
2707 /*                       WORK(IU) is LDA by N */
2708
2709                             ldwrku = *lda;
2710                         } else {
2711
2712 /*                       WORK(IU) is N by N */
2713
2714                             ldwrku = *n;
2715                         }
2716                         itau = iu + ldwrku * *n;
2717                         iwork = itau + *n;
2718
2719 /*                    Compute A=Q*R, copying result to U */
2720 /*                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) */
2721
2722                         i__2 = *lwork - iwork + 1;
2723                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2724                                 iwork], &i__2, &ierr);
2725                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2726                                 ldu);
2727
2728 /*                    Generate Q in U */
2729 /*                    (Workspace: need N*N + N + M, prefer N*N + N + M*NB) */
2730
2731                         i__2 = *lwork - iwork + 1;
2732                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2733                                 work[iwork], &i__2, &ierr);
2734
2735 /*                    Copy R to WORK(IU), zeroing out below it */
2736
2737                         dlacpy_("U", n, n, &a[a_offset], lda, &work[iu], &
2738                                 ldwrku);
2739                         i__2 = *n - 1;
2740                         i__3 = *n - 1;
2741                         dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
2742                                 1], &ldwrku);
2743                         ie = itau;
2744                         itauq = ie + *n;
2745                         itaup = itauq + *n;
2746                         iwork = itaup + *n;
2747
2748 /*                    Bidiagonalize R in WORK(IU), copying result to VT */
2749 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) */
2750
2751                         i__2 = *lwork - iwork + 1;
2752                         dgebrd_(n, n, &work[iu], &ldwrku, &s[1], &work[ie], &
2753                                 work[itauq], &work[itaup], &work[iwork], &
2754                                 i__2, &ierr);
2755                         dlacpy_("U", n, n, &work[iu], &ldwrku, &vt[vt_offset],
2756                                  ldvt);
2757
2758 /*                    Generate left bidiagonalizing vectors in WORK(IU) */
2759 /*                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) */
2760
2761                         i__2 = *lwork - iwork + 1;
2762                         dorgbr_("Q", n, n, n, &work[iu], &ldwrku, &work[itauq]
2763                                 , &work[iwork], &i__2, &ierr);
2764
2765 /*                    Generate right bidiagonalizing vectors in VT */
2766 /*                    (Workspace: need N*N + 4*N-1, */
2767 /*                                prefer N*N+3*N+(N-1)*NB) */
2768
2769                         i__2 = *lwork - iwork + 1;
2770                         dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2771                                 itaup], &work[iwork], &i__2, &ierr)
2772                                 ;
2773                         iwork = ie + *n;
2774
2775 /*                    Perform bidiagonal QR iteration, computing left */
2776 /*                    singular vectors of R in WORK(IU) and computing */
2777 /*                    right singular vectors of R in VT */
2778 /*                    (Workspace: need N*N + BDSPAC) */
2779
2780                         dbdsqr_("U", n, n, n, &c__0, &s[1], &work[ie], &vt[
2781                                 vt_offset], ldvt, &work[iu], &ldwrku, dum, &
2782                                 c__1, &work[iwork], info);
2783
2784 /*                    Multiply Q in U by left singular vectors of R in */
2785 /*                    WORK(IU), storing result in A */
2786 /*                    (Workspace: need N*N) */
2787
2788                         dgemm_("N", "N", m, n, n, &c_b79, &u[u_offset], ldu, &
2789                                 work[iu], &ldwrku, &c_b57, &a[a_offset], lda);
2790
2791 /*                    Copy left singular vectors of A from A to U */
2792
2793                         dlacpy_("F", m, n, &a[a_offset], lda, &u[u_offset],
2794                                 ldu);
2795
2796                     } else {
2797
2798 /*                    Insufficient workspace for a fast algorithm */
2799
2800                         itau = 1;
2801                         iwork = itau + *n;
2802
2803 /*                    Compute A=Q*R, copying result to U */
2804 /*                    (Workspace: need 2*N, prefer N + N*NB) */
2805
2806                         i__2 = *lwork - iwork + 1;
2807                         dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[
2808                                 iwork], &i__2, &ierr);
2809                         dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset],
2810                                 ldu);
2811
2812 /*                    Generate Q in U */
2813 /*                    (Workspace: need N + M, prefer N + M*NB) */
2814
2815                         i__2 = *lwork - iwork + 1;
2816                         dorgqr_(m, m, n, &u[u_offset], ldu, &work[itau], &
2817                                 work[iwork], &i__2, &ierr);
2818
2819 /*                    Copy R from A to VT, zeroing out below it */
2820
2821                         dlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset],
2822                                 ldvt);
2823                         if (*n > 1) {
2824                             i__2 = *n - 1;
2825                             i__3 = *n - 1;
2826                             dlaset_("L", &i__2, &i__3, &c_b57, &c_b57, &vt[
2827                                     vt_dim1 + 2], ldvt);
2828                         }
2829                         ie = itau;
2830                         itauq = ie + *n;
2831                         itaup = itauq + *n;
2832                         iwork = itaup + *n;
2833
2834 /*                    Bidiagonalize R in VT */
2835 /*                    (Workspace: need 4*N, prefer 3*N + 2*N*NB) */
2836
2837                         i__2 = *lwork - iwork + 1;
2838                         dgebrd_(n, n, &vt[vt_offset], ldvt, &s[1], &work[ie],
2839                                 &work[itauq], &work[itaup], &work[iwork], &
2840                                 i__2, &ierr);
2841
2842 /*                    Multiply Q in U by left bidiagonalizing vectors */
2843 /*                    in VT */
2844 /*                    (Workspace: need 3*N + M, prefer 3*N + M*NB) */
2845
2846                         i__2 = *lwork - iwork + 1;
2847                         dormbr_("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt,
2848                                 &work[itauq], &u[u_offset], ldu, &work[iwork],
2849                                  &i__2, &ierr);
2850
2851 /*                    Generate right bidiagonalizing vectors in VT */
2852 /*                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2853
2854                         i__2 = *lwork - iwork + 1;
2855                         dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[
2856                                 itaup], &work[iwork], &i__2, &ierr)
2857                                 ;
2858                         iwork = ie + *n;
2859
2860 /*                    Perform bidiagonal QR iteration, computing left */
2861 /*                    singular vectors of A in U and computing right */
2862 /*                    singular vectors of A in VT */
2863 /*                    (Workspace: need BDSPAC) */
2864
2865                         dbdsqr_("U", n, n, m, &c__0, &s[1], &work[ie], &vt[
2866                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
2867                                 c__1, &work[iwork], info);
2868
2869                     }
2870
2871                 }
2872
2873             }
2874
2875         } else {
2876
2877 /*           M .LT. MNTHR */
2878
2879 /*           Path 10 (M at least N, but not much larger) */
2880 /*           Reduce to bidiagonal form without QR decomposition */
2881
2882             ie = 1;
2883             itauq = ie + *n;
2884             itaup = itauq + *n;
2885             iwork = itaup + *n;
2886
2887 /*           Bidiagonalize A */
2888 /*           (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) */
2889
2890             i__2 = *lwork - iwork + 1;
2891             dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
2892                     work[itaup], &work[iwork], &i__2, &ierr);
2893             if (wntuas) {
2894
2895 /*              If left singular vectors desired in U, copy result to U */
2896 /*              and generate left bidiagonalizing vectors in U */
2897 /*              (Workspace: need 3*N + NCU, prefer 3*N + NCU*NB) */
2898
2899                 dlacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
2900                 if (wntus) {
2901                     ncu = *n;
2902                 }
2903                 if (wntua) {
2904                     ncu = *m;
2905                 }
2906                 i__2 = *lwork - iwork + 1;
2907                 dorgbr_("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], &
2908                         work[iwork], &i__2, &ierr);
2909             }
2910             if (wntvas) {
2911
2912 /*              If right singular vectors desired in VT, copy result to */
2913 /*              VT and generate right bidiagonalizing vectors in VT */
2914 /*              (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2915
2916                 dlacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2917                 i__2 = *lwork - iwork + 1;
2918                 dorgbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
2919                         work[iwork], &i__2, &ierr);
2920             }
2921             if (wntuo) {
2922
2923 /*              If left singular vectors desired in A, generate left */
2924 /*              bidiagonalizing vectors in A */
2925 /*              (Workspace: need 4*N, prefer 3*N + N*NB) */
2926
2927                 i__2 = *lwork - iwork + 1;
2928                 dorgbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
2929                         iwork], &i__2, &ierr);
2930             }
2931             if (wntvo) {
2932
2933 /*              If right singular vectors desired in A, generate right */
2934 /*              bidiagonalizing vectors in A */
2935 /*              (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) */
2936
2937                 i__2 = *lwork - iwork + 1;
2938                 dorgbr_("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[
2939                         iwork], &i__2, &ierr);
2940             }
2941             iwork = ie + *n;
2942             if (wntuas || wntuo) {
2943                 nru = *m;
2944             }
2945             if (wntun) {
2946                 nru = 0;
2947             }
2948             if (wntvas || wntvo) {
2949                 ncvt = *n;
2950             }
2951             if (wntvn) {
2952                 ncvt = 0;
2953             }
2954             if (! wntuo && ! wntvo) {
2955
2956 /*              Perform bidiagonal QR iteration, if desired, computing */
2957 /*              left singular vectors in U and computing right singular */
2958 /*              vectors in VT */
2959 /*              (Workspace: need BDSPAC) */
2960
2961                 dbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
2962                         vt_offset], ldvt, &u[u_offset], ldu, dum, &c__1, &
2963                         work[iwork], info);
2964             } else if (! wntuo && wntvo) {
2965
2966 /*              Perform bidiagonal QR iteration, if desired, computing */
2967 /*              left singular vectors in U and computing right singular */
2968 /*              vectors in A */
2969 /*              (Workspace: need BDSPAC) */
2970
2971                 dbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &a[
2972                         a_offset], lda, &u[u_offset], ldu, dum, &c__1, &work[
2973                         iwork], info);
2974             } else {
2975
2976 /*              Perform bidiagonal QR iteration, if desired, computing */
2977 /*              left singular vectors in A and computing right singular */
2978 /*              vectors in VT */
2979 /*              (Workspace: need BDSPAC) */
2980
2981                 dbdsqr_("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
2982                         vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
2983                         work[iwork], info);
2984             }
2985
2986         }
2987
2988     } else {
2989
2990 /*        A has more columns than rows. If A has sufficiently more */
2991 /*        columns than rows, first reduce using the LQ decomposition (if */
2992 /*        sufficient workspace available) */
2993
2994         if (*n >= mnthr) {
2995
2996             if (wntvn) {
2997
2998 /*              Path 1t(N much larger than M, JOBVT='N') */
2999 /*              No right singular vectors to be computed */
3000
3001                 itau = 1;
3002                 iwork = itau + *m;
3003
3004 /*              Compute A=L*Q */
3005 /*              (Workspace: need 2*M, prefer M + M*NB) */
3006
3007                 i__2 = *lwork - iwork + 1;
3008                 dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork], &
3009                         i__2, &ierr);
3010
3011 /*              Zero out above L */
3012
3013                 i__2 = *m - 1;
3014                 i__3 = *m - 1;
3015                 dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &a[(a_dim1 << 1) +
3016                         1], lda);
3017                 ie = 1;
3018                 itauq = ie + *m;
3019                 itaup = itauq + *m;
3020                 iwork = itaup + *m;
3021
3022 /*              Bidiagonalize L in A */
3023 /*              (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
3024
3025                 i__2 = *lwork - iwork + 1;
3026                 dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
3027                         itauq], &work[itaup], &work[iwork], &i__2, &ierr);
3028                 if (wntuo || wntuas) {
3029
3030 /*                 If left singular vectors desired, generate Q */
3031 /*                 (Workspace: need 4*M, prefer 3*M + M*NB) */
3032
3033                     i__2 = *lwork - iwork + 1;
3034                     dorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq], &
3035                             work[iwork], &i__2, &ierr);
3036                 }
3037                 iwork = ie + *m;
3038                 nru = 0;
3039                 if (wntuo || wntuas) {
3040                     nru = *m;
3041                 }
3042
3043 /*              Perform bidiagonal QR iteration, computing left singular */
3044 /*              vectors of A in A if desired */
3045 /*              (Workspace: need BDSPAC) */
3046
3047                 dbdsqr_("U", m, &c__0, &nru, &c__0, &s[1], &work[ie], dum, &
3048                         c__1, &a[a_offset], lda, dum, &c__1, &work[iwork],
3049                         info);
3050
3051 /*              If left singular vectors desired in U, copy them there */
3052
3053                 if (wntuas) {
3054                     dlacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
3055                 }
3056
3057             } else if (wntvo && wntun) {
3058
3059 /*              Path 2t(N much larger than M, JOBU='N', JOBVT='O') */
3060 /*              M right singular vectors to be overwritten on A and */
3061 /*              no left singular vectors to be computed */
3062
3063 /* Computing MAX */
3064                 i__2 = *m << 2;
3065                 if (*lwork >= *m * *m + f2cmax(i__2,bdspac)) {
3066
3067 /*                 Sufficient workspace for a fast algorithm */
3068
3069                     ir = 1;
3070 /* Computing MAX */
3071                     i__2 = wrkbl, i__3 = *lda * *n + *m;
3072                     if (*lwork >= f2cmax(i__2,i__3) + *lda * *m) {
3073
3074 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
3075
3076                         ldwrku = *lda;
3077                         chunk = *n;
3078                         ldwrkr = *lda;
3079                     } else /* if(complicated condition) */ {
3080 /* Computing MAX */
3081                         i__2 = wrkbl, i__3 = *lda * *n + *m;
3082                         if (*lwork >= f2cmax(i__2,i__3) + *m * *m) {
3083
3084 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
3085
3086                             ldwrku = *lda;
3087                             chunk = *n;
3088                             ldwrkr = *m;
3089                         } else {
3090
3091 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
3092
3093                             ldwrku = *m;
3094                             chunk = (*lwork - *m * *m - *m) / *m;
3095                             ldwrkr = *m;
3096                         }
3097                     }
3098                     itau = ir + ldwrkr * *m;
3099                     iwork = itau + *m;
3100
3101 /*                 Compute A=L*Q */
3102 /*                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3103
3104                     i__2 = *lwork - iwork + 1;
3105                     dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
3106                             , &i__2, &ierr);
3107
3108 /*                 Copy L to WORK(IR) and zero out above it */
3109
3110                     dlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
3111                     i__2 = *m - 1;
3112                     i__3 = *m - 1;
3113                     dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[ir +
3114                             ldwrkr], &ldwrkr);
3115
3116 /*                 Generate Q in A */
3117 /*                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3118
3119                     i__2 = *lwork - iwork + 1;
3120                     dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
3121                             iwork], &i__2, &ierr);
3122                     ie = itau;
3123                     itauq = ie + *m;
3124                     itaup = itauq + *m;
3125                     iwork = itaup + *m;
3126
3127 /*                 Bidiagonalize L in WORK(IR) */
3128 /*                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
3129
3130                     i__2 = *lwork - iwork + 1;
3131                     dgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &work[
3132                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
3133
3134 /*                 Generate right vectors bidiagonalizing L */
3135 /*                 (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) */
3136
3137                     i__2 = *lwork - iwork + 1;
3138                     dorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
3139                             work[iwork], &i__2, &ierr);
3140                     iwork = ie + *m;
3141
3142 /*                 Perform bidiagonal QR iteration, computing right */
3143 /*                 singular vectors of L in WORK(IR) */
3144 /*                 (Workspace: need M*M + BDSPAC) */
3145
3146                     dbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &work[
3147                             ir], &ldwrkr, dum, &c__1, dum, &c__1, &work[iwork]
3148                             , info);
3149                     iu = ie + *m;
3150
3151 /*                 Multiply right singular vectors of L in WORK(IR) by Q */
3152 /*                 in A, storing result in WORK(IU) and copying to A */
3153 /*                 (Workspace: need M*M + 2*M, prefer M*M + M*N + M) */
3154
3155                     i__2 = *n;
3156                     i__3 = chunk;
3157                     for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
3158                              i__3) {
3159 /* Computing MIN */
3160                         i__4 = *n - i__ + 1;
3161                         blk = f2cmin(i__4,chunk);
3162                         dgemm_("N", "N", m, &blk, m, &c_b79, &work[ir], &
3163                                 ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b57, &
3164                                 work[iu], &ldwrku);
3165                         dlacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ *
3166                                 a_dim1 + 1], lda);
3167 /* L30: */
3168                     }
3169
3170                 } else {
3171
3172 /*                 Insufficient workspace for a fast algorithm */
3173
3174                     ie = 1;
3175                     itauq = ie + *m;
3176                     itaup = itauq + *m;
3177                     iwork = itaup + *m;
3178
3179 /*                 Bidiagonalize A */
3180 /*                 (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) */
3181
3182                     i__3 = *lwork - iwork + 1;
3183                     dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[
3184                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
3185
3186 /*                 Generate right vectors bidiagonalizing A */
3187 /*                 (Workspace: need 4*M, prefer 3*M + M*NB) */
3188
3189                     i__3 = *lwork - iwork + 1;
3190                     dorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
3191                             work[iwork], &i__3, &ierr);
3192                     iwork = ie + *m;
3193
3194 /*                 Perform bidiagonal QR iteration, computing right */
3195 /*                 singular vectors of A in A */
3196 /*                 (Workspace: need BDSPAC) */
3197
3198                     dbdsqr_("L", m, n, &c__0, &c__0, &s[1], &work[ie], &a[
3199                             a_offset], lda, dum, &c__1, dum, &c__1, &work[
3200                             iwork], info);
3201
3202                 }
3203
3204             } else if (wntvo && wntuas) {
3205
3206 /*              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') */
3207 /*              M right singular vectors to be overwritten on A and */
3208 /*              M left singular vectors to be computed in U */
3209
3210 /* Computing MAX */
3211                 i__3 = *m << 2;
3212                 if (*lwork >= *m * *m + f2cmax(i__3,bdspac)) {
3213
3214 /*                 Sufficient workspace for a fast algorithm */
3215
3216                     ir = 1;
3217 /* Computing MAX */
3218                     i__3 = wrkbl, i__2 = *lda * *n + *m;
3219                     if (*lwork >= f2cmax(i__3,i__2) + *lda * *m) {
3220
3221 /*                    WORK(IU) is LDA by N and WORK(IR) is LDA by M */
3222
3223                         ldwrku = *lda;
3224                         chunk = *n;
3225                         ldwrkr = *lda;
3226                     } else /* if(complicated condition) */ {
3227 /* Computing MAX */
3228                         i__3 = wrkbl, i__2 = *lda * *n + *m;
3229                         if (*lwork >= f2cmax(i__3,i__2) + *m * *m) {
3230
3231 /*                    WORK(IU) is LDA by N and WORK(IR) is M by M */
3232
3233                             ldwrku = *lda;
3234                             chunk = *n;
3235                             ldwrkr = *m;
3236                         } else {
3237
3238 /*                    WORK(IU) is M by CHUNK and WORK(IR) is M by M */
3239
3240                             ldwrku = *m;
3241                             chunk = (*lwork - *m * *m - *m) / *m;
3242                             ldwrkr = *m;
3243                         }
3244                     }
3245                     itau = ir + ldwrkr * *m;
3246                     iwork = itau + *m;
3247
3248 /*                 Compute A=L*Q */
3249 /*                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3250
3251                     i__3 = *lwork - iwork + 1;
3252                     dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
3253                             , &i__3, &ierr);
3254
3255 /*                 Copy L to U, zeroing about above it */
3256
3257                     dlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
3258                     i__3 = *m - 1;
3259                     i__2 = *m - 1;
3260                     dlaset_("U", &i__3, &i__2, &c_b57, &c_b57, &u[(u_dim1 <<
3261                             1) + 1], ldu);
3262
3263 /*                 Generate Q in A */
3264 /*                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3265
3266                     i__3 = *lwork - iwork + 1;
3267                     dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
3268                             iwork], &i__3, &ierr);
3269                     ie = itau;
3270                     itauq = ie + *m;
3271                     itaup = itauq + *m;
3272                     iwork = itaup + *m;
3273
3274 /*                 Bidiagonalize L in U, copying result to WORK(IR) */
3275 /*                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
3276
3277                     i__3 = *lwork - iwork + 1;
3278                     dgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[
3279                             itauq], &work[itaup], &work[iwork], &i__3, &ierr);
3280                     dlacpy_("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);
3281
3282 /*                 Generate right vectors bidiagonalizing L in WORK(IR) */
3283 /*                 (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) */
3284
3285                     i__3 = *lwork - iwork + 1;
3286                     dorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &
3287                             work[iwork], &i__3, &ierr);
3288
3289 /*                 Generate left vectors bidiagonalizing L in U */
3290 /*                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) */
3291
3292                     i__3 = *lwork - iwork + 1;
3293                     dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
3294                             work[iwork], &i__3, &ierr);
3295                     iwork = ie + *m;
3296
3297 /*                 Perform bidiagonal QR iteration, computing left */
3298 /*                 singular vectors of L in U, and computing right */
3299 /*                 singular vectors of L in WORK(IR) */
3300 /*                 (Workspace: need M*M + BDSPAC) */
3301
3302                     dbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[ir],
3303                             &ldwrkr, &u[u_offset], ldu, dum, &c__1, &work[
3304                             iwork], info);
3305                     iu = ie + *m;
3306
3307 /*                 Multiply right singular vectors of L in WORK(IR) by Q */
3308 /*                 in A, storing result in WORK(IU) and copying to A */
3309 /*                 (Workspace: need M*M + 2*M, prefer M*M + M*N + M)) */
3310
3311                     i__3 = *n;
3312                     i__2 = chunk;
3313                     for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ +=
3314                              i__2) {
3315 /* Computing MIN */
3316                         i__4 = *n - i__ + 1;
3317                         blk = f2cmin(i__4,chunk);
3318                         dgemm_("N", "N", m, &blk, m, &c_b79, &work[ir], &
3319                                 ldwrkr, &a[i__ * a_dim1 + 1], lda, &c_b57, &
3320                                 work[iu], &ldwrku);
3321                         dlacpy_("F", m, &blk, &work[iu], &ldwrku, &a[i__ *
3322                                 a_dim1 + 1], lda);
3323 /* L40: */
3324                     }
3325
3326                 } else {
3327
3328 /*                 Insufficient workspace for a fast algorithm */
3329
3330                     itau = 1;
3331                     iwork = itau + *m;
3332
3333 /*                 Compute A=L*Q */
3334 /*                 (Workspace: need 2*M, prefer M + M*NB) */
3335
3336                     i__2 = *lwork - iwork + 1;
3337                     dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[iwork]
3338                             , &i__2, &ierr);
3339
3340 /*                 Copy L to U, zeroing out above it */
3341
3342                     dlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
3343                     i__2 = *m - 1;
3344                     i__3 = *m - 1;
3345                     dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &u[(u_dim1 <<
3346                             1) + 1], ldu);
3347
3348 /*                 Generate Q in A */
3349 /*                 (Workspace: need 2*M, prefer M + M*NB) */
3350
3351                     i__2 = *lwork - iwork + 1;
3352                     dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[
3353                             iwork], &i__2, &ierr);
3354                     ie = itau;
3355                     itauq = ie + *m;
3356                     itaup = itauq + *m;
3357                     iwork = itaup + *m;
3358
3359 /*                 Bidiagonalize L in U */
3360 /*                 (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
3361
3362                     i__2 = *lwork - iwork + 1;
3363                     dgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[
3364                             itauq], &work[itaup], &work[iwork], &i__2, &ierr);
3365
3366 /*                 Multiply right vectors bidiagonalizing L by Q in A */
3367 /*                 (Workspace: need 3*M + N, prefer 3*M + N*NB) */
3368
3369                     i__2 = *lwork - iwork + 1;
3370                     dormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &work[
3371                             itaup], &a[a_offset], lda, &work[iwork], &i__2, &
3372                             ierr);
3373
3374 /*                 Generate left vectors bidiagonalizing L in U */
3375 /*                 (Workspace: need 4*M, prefer 3*M + M*NB) */
3376
3377                     i__2 = *lwork - iwork + 1;
3378                     dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &
3379                             work[iwork], &i__2, &ierr);
3380                     iwork = ie + *m;
3381
3382 /*                 Perform bidiagonal QR iteration, computing left */
3383 /*                 singular vectors of A in U and computing right */
3384 /*                 singular vectors of A in A */
3385 /*                 (Workspace: need BDSPAC) */
3386
3387                     dbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &a[
3388                             a_offset], lda, &u[u_offset], ldu, dum, &c__1, &
3389                             work[iwork], info);
3390
3391                 }
3392
3393             } else if (wntvs) {
3394
3395                 if (wntun) {
3396
3397 /*                 Path 4t(N much larger than M, JOBU='N', JOBVT='S') */
3398 /*                 M right singular vectors to be computed in VT and */
3399 /*                 no left singular vectors to be computed */
3400
3401 /* Computing MAX */
3402                     i__2 = *m << 2;
3403                     if (*lwork >= *m * *m + f2cmax(i__2,bdspac)) {
3404
3405 /*                    Sufficient workspace for a fast algorithm */
3406
3407                         ir = 1;
3408                         if (*lwork >= wrkbl + *lda * *m) {
3409
3410 /*                       WORK(IR) is LDA by M */
3411
3412                             ldwrkr = *lda;
3413                         } else {
3414
3415 /*                       WORK(IR) is M by M */
3416
3417                             ldwrkr = *m;
3418                         }
3419                         itau = ir + ldwrkr * *m;
3420                         iwork = itau + *m;
3421
3422 /*                    Compute A=L*Q */
3423 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3424
3425                         i__2 = *lwork - iwork + 1;
3426                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3427                                 iwork], &i__2, &ierr);
3428
3429 /*                    Copy L to WORK(IR), zeroing out above it */
3430
3431                         dlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3432                                 ldwrkr);
3433                         i__2 = *m - 1;
3434                         i__3 = *m - 1;
3435                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[ir +
3436                                 ldwrkr], &ldwrkr);
3437
3438 /*                    Generate Q in A */
3439 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3440
3441                         i__2 = *lwork - iwork + 1;
3442                         dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3443                                 work[iwork], &i__2, &ierr);
3444                         ie = itau;
3445                         itauq = ie + *m;
3446                         itaup = itauq + *m;
3447                         iwork = itaup + *m;
3448
3449 /*                    Bidiagonalize L in WORK(IR) */
3450 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
3451
3452                         i__2 = *lwork - iwork + 1;
3453                         dgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &
3454                                 work[itauq], &work[itaup], &work[iwork], &
3455                                 i__2, &ierr);
3456
3457 /*                    Generate right vectors bidiagonalizing L in */
3458 /*                    WORK(IR) */
3459 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) */
3460
3461                         i__2 = *lwork - iwork + 1;
3462                         dorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3463                                 , &work[iwork], &i__2, &ierr);
3464                         iwork = ie + *m;
3465
3466 /*                    Perform bidiagonal QR iteration, computing right */
3467 /*                    singular vectors of L in WORK(IR) */
3468 /*                    (Workspace: need M*M + BDSPAC) */
3469
3470                         dbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &
3471                                 work[ir], &ldwrkr, dum, &c__1, dum, &c__1, &
3472                                 work[iwork], info);
3473
3474 /*                    Multiply right singular vectors of L in WORK(IR) by */
3475 /*                    Q in A, storing result in VT */
3476 /*                    (Workspace: need M*M) */
3477
3478                         dgemm_("N", "N", m, n, m, &c_b79, &work[ir], &ldwrkr,
3479                                 &a[a_offset], lda, &c_b57, &vt[vt_offset],
3480                                 ldvt);
3481
3482                     } else {
3483
3484 /*                    Insufficient workspace for a fast algorithm */
3485
3486                         itau = 1;
3487                         iwork = itau + *m;
3488
3489 /*                    Compute A=L*Q */
3490 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3491
3492                         i__2 = *lwork - iwork + 1;
3493                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3494                                 iwork], &i__2, &ierr);
3495
3496 /*                    Copy result to VT */
3497
3498                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3499                                 ldvt);
3500
3501 /*                    Generate Q in VT */
3502 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3503
3504                         i__2 = *lwork - iwork + 1;
3505                         dorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3506                                 work[iwork], &i__2, &ierr);
3507                         ie = itau;
3508                         itauq = ie + *m;
3509                         itaup = itauq + *m;
3510                         iwork = itaup + *m;
3511
3512 /*                    Zero out above L in A */
3513
3514                         i__2 = *m - 1;
3515                         i__3 = *m - 1;
3516                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &a[(a_dim1
3517                                 << 1) + 1], lda);
3518
3519 /*                    Bidiagonalize L in A */
3520 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
3521
3522                         i__2 = *lwork - iwork + 1;
3523                         dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
3524                                 work[itauq], &work[itaup], &work[iwork], &
3525                                 i__2, &ierr);
3526
3527 /*                    Multiply right vectors bidiagonalizing L by Q in VT */
3528 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
3529
3530                         i__2 = *lwork - iwork + 1;
3531                         dormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
3532                                 work[itaup], &vt[vt_offset], ldvt, &work[
3533                                 iwork], &i__2, &ierr);
3534                         iwork = ie + *m;
3535
3536 /*                    Perform bidiagonal QR iteration, computing right */
3537 /*                    singular vectors of A in VT */
3538 /*                    (Workspace: need BDSPAC) */
3539
3540                         dbdsqr_("U", m, n, &c__0, &c__0, &s[1], &work[ie], &
3541                                 vt[vt_offset], ldvt, dum, &c__1, dum, &c__1, &
3542                                 work[iwork], info);
3543
3544                     }
3545
3546                 } else if (wntuo) {
3547
3548 /*                 Path 5t(N much larger than M, JOBU='O', JOBVT='S') */
3549 /*                 M right singular vectors to be computed in VT and */
3550 /*                 M left singular vectors to be overwritten on A */
3551
3552 /* Computing MAX */
3553                     i__2 = *m << 2;
3554                     if (*lwork >= (*m << 1) * *m + f2cmax(i__2,bdspac)) {
3555
3556 /*                    Sufficient workspace for a fast algorithm */
3557
3558                         iu = 1;
3559                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
3560
3561 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
3562
3563                             ldwrku = *lda;
3564                             ir = iu + ldwrku * *m;
3565                             ldwrkr = *lda;
3566                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
3567
3568 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
3569
3570                             ldwrku = *lda;
3571                             ir = iu + ldwrku * *m;
3572                             ldwrkr = *m;
3573                         } else {
3574
3575 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
3576
3577                             ldwrku = *m;
3578                             ir = iu + ldwrku * *m;
3579                             ldwrkr = *m;
3580                         }
3581                         itau = ir + ldwrkr * *m;
3582                         iwork = itau + *m;
3583
3584 /*                    Compute A=L*Q */
3585 /*                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) */
3586
3587                         i__2 = *lwork - iwork + 1;
3588                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3589                                 iwork], &i__2, &ierr);
3590
3591 /*                    Copy L to WORK(IU), zeroing out below it */
3592
3593                         dlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3594                                 ldwrku);
3595                         i__2 = *m - 1;
3596                         i__3 = *m - 1;
3597                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
3598                                 ldwrku], &ldwrku);
3599
3600 /*                    Generate Q in A */
3601 /*                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) */
3602
3603                         i__2 = *lwork - iwork + 1;
3604                         dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3605                                 work[iwork], &i__2, &ierr);
3606                         ie = itau;
3607                         itauq = ie + *m;
3608                         itaup = itauq + *m;
3609                         iwork = itaup + *m;
3610
3611 /*                    Bidiagonalize L in WORK(IU), copying result to */
3612 /*                    WORK(IR) */
3613 /*                    (Workspace: need 2*M*M + 4*M, */
3614 /*                                prefer 2*M*M+3*M+2*M*NB) */
3615
3616                         i__2 = *lwork - iwork + 1;
3617                         dgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
3618                                 work[itauq], &work[itaup], &work[iwork], &
3619                                 i__2, &ierr);
3620                         dlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
3621                                 ldwrkr);
3622
3623 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
3624 /*                    (Workspace: need 2*M*M + 4*M-1, */
3625 /*                                prefer 2*M*M+3*M+(M-1)*NB) */
3626
3627                         i__2 = *lwork - iwork + 1;
3628                         dorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3629                                 , &work[iwork], &i__2, &ierr);
3630
3631 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
3632 /*                    (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) */
3633
3634                         i__2 = *lwork - iwork + 1;
3635                         dorgbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
3636                                 , &work[iwork], &i__2, &ierr);
3637                         iwork = ie + *m;
3638
3639 /*                    Perform bidiagonal QR iteration, computing left */
3640 /*                    singular vectors of L in WORK(IR) and computing */
3641 /*                    right singular vectors of L in WORK(IU) */
3642 /*                    (Workspace: need 2*M*M + BDSPAC) */
3643
3644                         dbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
3645                                 iu], &ldwrku, &work[ir], &ldwrkr, dum, &c__1,
3646                                 &work[iwork], info);
3647
3648 /*                    Multiply right singular vectors of L in WORK(IU) by */
3649 /*                    Q in A, storing result in VT */
3650 /*                    (Workspace: need M*M) */
3651
3652                         dgemm_("N", "N", m, n, m, &c_b79, &work[iu], &ldwrku,
3653                                 &a[a_offset], lda, &c_b57, &vt[vt_offset],
3654                                 ldvt);
3655
3656 /*                    Copy left singular vectors of L to A */
3657 /*                    (Workspace: need M*M) */
3658
3659                         dlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
3660                                 lda);
3661
3662                     } else {
3663
3664 /*                    Insufficient workspace for a fast algorithm */
3665
3666                         itau = 1;
3667                         iwork = itau + *m;
3668
3669 /*                    Compute A=L*Q, copying result to VT */
3670 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3671
3672                         i__2 = *lwork - iwork + 1;
3673                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3674                                 iwork], &i__2, &ierr);
3675                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3676                                 ldvt);
3677
3678 /*                    Generate Q in VT */
3679 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3680
3681                         i__2 = *lwork - iwork + 1;
3682                         dorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3683                                 work[iwork], &i__2, &ierr);
3684                         ie = itau;
3685                         itauq = ie + *m;
3686                         itaup = itauq + *m;
3687                         iwork = itaup + *m;
3688
3689 /*                    Zero out above L in A */
3690
3691                         i__2 = *m - 1;
3692                         i__3 = *m - 1;
3693                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &a[(a_dim1
3694                                 << 1) + 1], lda);
3695
3696 /*                    Bidiagonalize L in A */
3697 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
3698
3699                         i__2 = *lwork - iwork + 1;
3700                         dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
3701                                 work[itauq], &work[itaup], &work[iwork], &
3702                                 i__2, &ierr);
3703
3704 /*                    Multiply right vectors bidiagonalizing L by Q in VT */
3705 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
3706
3707                         i__2 = *lwork - iwork + 1;
3708                         dormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
3709                                 work[itaup], &vt[vt_offset], ldvt, &work[
3710                                 iwork], &i__2, &ierr);
3711
3712 /*                    Generate left bidiagonalizing vectors of L in A */
3713 /*                    (Workspace: need 4*M, prefer 3*M + M*NB) */
3714
3715                         i__2 = *lwork - iwork + 1;
3716                         dorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
3717                                  &work[iwork], &i__2, &ierr);
3718                         iwork = ie + *m;
3719
3720 /*                    Perform bidiagonal QR iteration, compute left */
3721 /*                    singular vectors of A in A and compute right */
3722 /*                    singular vectors of A in VT */
3723 /*                    (Workspace: need BDSPAC) */
3724
3725                         dbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
3726                                 vt_offset], ldvt, &a[a_offset], lda, dum, &
3727                                 c__1, &work[iwork], info);
3728
3729                     }
3730
3731                 } else if (wntuas) {
3732
3733 /*                 Path 6t(N much larger than M, JOBU='S' or 'A', */
3734 /*                         JOBVT='S') */
3735 /*                 M right singular vectors to be computed in VT and */
3736 /*                 M left singular vectors to be computed in U */
3737
3738 /* Computing MAX */
3739                     i__2 = *m << 2;
3740                     if (*lwork >= *m * *m + f2cmax(i__2,bdspac)) {
3741
3742 /*                    Sufficient workspace for a fast algorithm */
3743
3744                         iu = 1;
3745                         if (*lwork >= wrkbl + *lda * *m) {
3746
3747 /*                       WORK(IU) is LDA by N */
3748
3749                             ldwrku = *lda;
3750                         } else {
3751
3752 /*                       WORK(IU) is LDA by M */
3753
3754                             ldwrku = *m;
3755                         }
3756                         itau = iu + ldwrku * *m;
3757                         iwork = itau + *m;
3758
3759 /*                    Compute A=L*Q */
3760 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3761
3762                         i__2 = *lwork - iwork + 1;
3763                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3764                                 iwork], &i__2, &ierr);
3765
3766 /*                    Copy L to WORK(IU), zeroing out above it */
3767
3768                         dlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
3769                                 ldwrku);
3770                         i__2 = *m - 1;
3771                         i__3 = *m - 1;
3772                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
3773                                 ldwrku], &ldwrku);
3774
3775 /*                    Generate Q in A */
3776 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3777
3778                         i__2 = *lwork - iwork + 1;
3779                         dorglq_(m, n, m, &a[a_offset], lda, &work[itau], &
3780                                 work[iwork], &i__2, &ierr);
3781                         ie = itau;
3782                         itauq = ie + *m;
3783                         itaup = itauq + *m;
3784                         iwork = itaup + *m;
3785
3786 /*                    Bidiagonalize L in WORK(IU), copying result to U */
3787 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
3788
3789                         i__2 = *lwork - iwork + 1;
3790                         dgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
3791                                 work[itauq], &work[itaup], &work[iwork], &
3792                                 i__2, &ierr);
3793                         dlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
3794                                 ldu);
3795
3796 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
3797 /*                    (Workspace: need M*M + 4*M-1, */
3798 /*                                prefer M*M+3*M+(M-1)*NB) */
3799
3800                         i__2 = *lwork - iwork + 1;
3801                         dorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
3802                                 , &work[iwork], &i__2, &ierr);
3803
3804 /*                    Generate left bidiagonalizing vectors in U */
3805 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) */
3806
3807                         i__2 = *lwork - iwork + 1;
3808                         dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3809                                  &work[iwork], &i__2, &ierr);
3810                         iwork = ie + *m;
3811
3812 /*                    Perform bidiagonal QR iteration, computing left */
3813 /*                    singular vectors of L in U and computing right */
3814 /*                    singular vectors of L in WORK(IU) */
3815 /*                    (Workspace: need M*M + BDSPAC) */
3816
3817                         dbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
3818                                 iu], &ldwrku, &u[u_offset], ldu, dum, &c__1, &
3819                                 work[iwork], info);
3820
3821 /*                    Multiply right singular vectors of L in WORK(IU) by */
3822 /*                    Q in A, storing result in VT */
3823 /*                    (Workspace: need M*M) */
3824
3825                         dgemm_("N", "N", m, n, m, &c_b79, &work[iu], &ldwrku,
3826                                 &a[a_offset], lda, &c_b57, &vt[vt_offset],
3827                                 ldvt);
3828
3829                     } else {
3830
3831 /*                    Insufficient workspace for a fast algorithm */
3832
3833                         itau = 1;
3834                         iwork = itau + *m;
3835
3836 /*                    Compute A=L*Q, copying result to VT */
3837 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3838
3839                         i__2 = *lwork - iwork + 1;
3840                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3841                                 iwork], &i__2, &ierr);
3842                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3843                                 ldvt);
3844
3845 /*                    Generate Q in VT */
3846 /*                    (Workspace: need 2*M, prefer M + M*NB) */
3847
3848                         i__2 = *lwork - iwork + 1;
3849                         dorglq_(m, n, m, &vt[vt_offset], ldvt, &work[itau], &
3850                                 work[iwork], &i__2, &ierr);
3851
3852 /*                    Copy L to U, zeroing out above it */
3853
3854                         dlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
3855                                 ldu);
3856                         i__2 = *m - 1;
3857                         i__3 = *m - 1;
3858                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &u[(u_dim1
3859                                 << 1) + 1], ldu);
3860                         ie = itau;
3861                         itauq = ie + *m;
3862                         itaup = itauq + *m;
3863                         iwork = itaup + *m;
3864
3865 /*                    Bidiagonalize L in U */
3866 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
3867
3868                         i__2 = *lwork - iwork + 1;
3869                         dgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &
3870                                 work[itauq], &work[itaup], &work[iwork], &
3871                                 i__2, &ierr);
3872
3873 /*                    Multiply right bidiagonalizing vectors in U by Q */
3874 /*                    in VT */
3875 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
3876
3877                         i__2 = *lwork - iwork + 1;
3878                         dormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &
3879                                 work[itaup], &vt[vt_offset], ldvt, &work[
3880                                 iwork], &i__2, &ierr);
3881
3882 /*                    Generate left bidiagonalizing vectors in U */
3883 /*                    (Workspace: need 4*M, prefer 3*M + M*NB) */
3884
3885                         i__2 = *lwork - iwork + 1;
3886                         dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
3887                                  &work[iwork], &i__2, &ierr);
3888                         iwork = ie + *m;
3889
3890 /*                    Perform bidiagonal QR iteration, computing left */
3891 /*                    singular vectors of A in U and computing right */
3892 /*                    singular vectors of A in VT */
3893 /*                    (Workspace: need BDSPAC) */
3894
3895                         dbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
3896                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
3897                                 c__1, &work[iwork], info);
3898
3899                     }
3900
3901                 }
3902
3903             } else if (wntva) {
3904
3905                 if (wntun) {
3906
3907 /*                 Path 7t(N much larger than M, JOBU='N', JOBVT='A') */
3908 /*                 N right singular vectors to be computed in VT and */
3909 /*                 no left singular vectors to be computed */
3910
3911 /* Computing MAX */
3912                     i__2 = *n + *m, i__3 = *m << 2, i__2 = f2cmax(i__2,i__3);
3913                     if (*lwork >= *m * *m + f2cmax(i__2,bdspac)) {
3914
3915 /*                    Sufficient workspace for a fast algorithm */
3916
3917                         ir = 1;
3918                         if (*lwork >= wrkbl + *lda * *m) {
3919
3920 /*                       WORK(IR) is LDA by M */
3921
3922                             ldwrkr = *lda;
3923                         } else {
3924
3925 /*                       WORK(IR) is M by M */
3926
3927                             ldwrkr = *m;
3928                         }
3929                         itau = ir + ldwrkr * *m;
3930                         iwork = itau + *m;
3931
3932 /*                    Compute A=L*Q, copying result to VT */
3933 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
3934
3935                         i__2 = *lwork - iwork + 1;
3936                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
3937                                 iwork], &i__2, &ierr);
3938                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
3939                                 ldvt);
3940
3941 /*                    Copy L to WORK(IR), zeroing out above it */
3942
3943                         dlacpy_("L", m, m, &a[a_offset], lda, &work[ir], &
3944                                 ldwrkr);
3945                         i__2 = *m - 1;
3946                         i__3 = *m - 1;
3947                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[ir +
3948                                 ldwrkr], &ldwrkr);
3949
3950 /*                    Generate Q in VT */
3951 /*                    (Workspace: need M*M + M + N, prefer M*M + M + N*NB) */
3952
3953                         i__2 = *lwork - iwork + 1;
3954                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
3955                                 work[iwork], &i__2, &ierr);
3956                         ie = itau;
3957                         itauq = ie + *m;
3958                         itaup = itauq + *m;
3959                         iwork = itaup + *m;
3960
3961 /*                    Bidiagonalize L in WORK(IR) */
3962 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
3963
3964                         i__2 = *lwork - iwork + 1;
3965                         dgebrd_(m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &
3966                                 work[itauq], &work[itaup], &work[iwork], &
3967                                 i__2, &ierr);
3968
3969 /*                    Generate right bidiagonalizing vectors in WORK(IR) */
3970 /*                    (Workspace: need M*M + 4*M-1, */
3971 /*                                prefer M*M+3*M+(M-1)*NB) */
3972
3973                         i__2 = *lwork - iwork + 1;
3974                         dorgbr_("P", m, m, m, &work[ir], &ldwrkr, &work[itaup]
3975                                 , &work[iwork], &i__2, &ierr);
3976                         iwork = ie + *m;
3977
3978 /*                    Perform bidiagonal QR iteration, computing right */
3979 /*                    singular vectors of L in WORK(IR) */
3980 /*                    (Workspace: need M*M + BDSPAC) */
3981
3982                         dbdsqr_("U", m, m, &c__0, &c__0, &s[1], &work[ie], &
3983                                 work[ir], &ldwrkr, dum, &c__1, dum, &c__1, &
3984                                 work[iwork], info);
3985
3986 /*                    Multiply right singular vectors of L in WORK(IR) by */
3987 /*                    Q in VT, storing result in A */
3988 /*                    (Workspace: need M*M) */
3989
3990                         dgemm_("N", "N", m, n, m, &c_b79, &work[ir], &ldwrkr,
3991                                 &vt[vt_offset], ldvt, &c_b57, &a[a_offset],
3992                                 lda);
3993
3994 /*                    Copy right singular vectors of A from A to VT */
3995
3996                         dlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
3997                                 ldvt);
3998
3999                     } else {
4000
4001 /*                    Insufficient workspace for a fast algorithm */
4002
4003                         itau = 1;
4004                         iwork = itau + *m;
4005
4006 /*                    Compute A=L*Q, copying result to VT */
4007 /*                    (Workspace: need 2*M, prefer M + M*NB) */
4008
4009                         i__2 = *lwork - iwork + 1;
4010                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
4011                                 iwork], &i__2, &ierr);
4012                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4013                                 ldvt);
4014
4015 /*                    Generate Q in VT */
4016 /*                    (Workspace: need M + N, prefer M + N*NB) */
4017
4018                         i__2 = *lwork - iwork + 1;
4019                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4020                                 work[iwork], &i__2, &ierr);
4021                         ie = itau;
4022                         itauq = ie + *m;
4023                         itaup = itauq + *m;
4024                         iwork = itaup + *m;
4025
4026 /*                    Zero out above L in A */
4027
4028                         i__2 = *m - 1;
4029                         i__3 = *m - 1;
4030                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &a[(a_dim1
4031                                 << 1) + 1], lda);
4032
4033 /*                    Bidiagonalize L in A */
4034 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
4035
4036                         i__2 = *lwork - iwork + 1;
4037                         dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
4038                                 work[itauq], &work[itaup], &work[iwork], &
4039                                 i__2, &ierr);
4040
4041 /*                    Multiply right bidiagonalizing vectors in A by Q */
4042 /*                    in VT */
4043 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
4044
4045                         i__2 = *lwork - iwork + 1;
4046                         dormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
4047                                 work[itaup], &vt[vt_offset], ldvt, &work[
4048                                 iwork], &i__2, &ierr);
4049                         iwork = ie + *m;
4050
4051 /*                    Perform bidiagonal QR iteration, computing right */
4052 /*                    singular vectors of A in VT */
4053 /*                    (Workspace: need BDSPAC) */
4054
4055                         dbdsqr_("U", m, n, &c__0, &c__0, &s[1], &work[ie], &
4056                                 vt[vt_offset], ldvt, dum, &c__1, dum, &c__1, &
4057                                 work[iwork], info);
4058
4059                     }
4060
4061                 } else if (wntuo) {
4062
4063 /*                 Path 8t(N much larger than M, JOBU='O', JOBVT='A') */
4064 /*                 N right singular vectors to be computed in VT and */
4065 /*                 M left singular vectors to be overwritten on A */
4066
4067 /* Computing MAX */
4068                     i__2 = *n + *m, i__3 = *m << 2, i__2 = f2cmax(i__2,i__3);
4069                     if (*lwork >= (*m << 1) * *m + f2cmax(i__2,bdspac)) {
4070
4071 /*                    Sufficient workspace for a fast algorithm */
4072
4073                         iu = 1;
4074                         if (*lwork >= wrkbl + (*lda << 1) * *m) {
4075
4076 /*                       WORK(IU) is LDA by M and WORK(IR) is LDA by M */
4077
4078                             ldwrku = *lda;
4079                             ir = iu + ldwrku * *m;
4080                             ldwrkr = *lda;
4081                         } else if (*lwork >= wrkbl + (*lda + *m) * *m) {
4082
4083 /*                       WORK(IU) is LDA by M and WORK(IR) is M by M */
4084
4085                             ldwrku = *lda;
4086                             ir = iu + ldwrku * *m;
4087                             ldwrkr = *m;
4088                         } else {
4089
4090 /*                       WORK(IU) is M by M and WORK(IR) is M by M */
4091
4092                             ldwrku = *m;
4093                             ir = iu + ldwrku * *m;
4094                             ldwrkr = *m;
4095                         }
4096                         itau = ir + ldwrkr * *m;
4097                         iwork = itau + *m;
4098
4099 /*                    Compute A=L*Q, copying result to VT */
4100 /*                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) */
4101
4102                         i__2 = *lwork - iwork + 1;
4103                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
4104                                 iwork], &i__2, &ierr);
4105                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4106                                 ldvt);
4107
4108 /*                    Generate Q in VT */
4109 /*                    (Workspace: need 2*M*M + M + N, prefer 2*M*M + M + N*NB) */
4110
4111                         i__2 = *lwork - iwork + 1;
4112                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4113                                 work[iwork], &i__2, &ierr);
4114
4115 /*                    Copy L to WORK(IU), zeroing out above it */
4116
4117                         dlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
4118                                 ldwrku);
4119                         i__2 = *m - 1;
4120                         i__3 = *m - 1;
4121                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
4122                                 ldwrku], &ldwrku);
4123                         ie = itau;
4124                         itauq = ie + *m;
4125                         itaup = itauq + *m;
4126                         iwork = itaup + *m;
4127
4128 /*                    Bidiagonalize L in WORK(IU), copying result to */
4129 /*                    WORK(IR) */
4130 /*                    (Workspace: need 2*M*M + 4*M, */
4131 /*                                prefer 2*M*M+3*M+2*M*NB) */
4132
4133                         i__2 = *lwork - iwork + 1;
4134                         dgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
4135                                 work[itauq], &work[itaup], &work[iwork], &
4136                                 i__2, &ierr);
4137                         dlacpy_("L", m, m, &work[iu], &ldwrku, &work[ir], &
4138                                 ldwrkr);
4139
4140 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
4141 /*                    (Workspace: need 2*M*M + 4*M-1, */
4142 /*                                prefer 2*M*M+3*M+(M-1)*NB) */
4143
4144                         i__2 = *lwork - iwork + 1;
4145                         dorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
4146                                 , &work[iwork], &i__2, &ierr);
4147
4148 /*                    Generate left bidiagonalizing vectors in WORK(IR) */
4149 /*                    (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) */
4150
4151                         i__2 = *lwork - iwork + 1;
4152                         dorgbr_("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq]
4153                                 , &work[iwork], &i__2, &ierr);
4154                         iwork = ie + *m;
4155
4156 /*                    Perform bidiagonal QR iteration, computing left */
4157 /*                    singular vectors of L in WORK(IR) and computing */
4158 /*                    right singular vectors of L in WORK(IU) */
4159 /*                    (Workspace: need 2*M*M + BDSPAC) */
4160
4161                         dbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
4162                                 iu], &ldwrku, &work[ir], &ldwrkr, dum, &c__1,
4163                                 &work[iwork], info);
4164
4165 /*                    Multiply right singular vectors of L in WORK(IU) by */
4166 /*                    Q in VT, storing result in A */
4167 /*                    (Workspace: need M*M) */
4168
4169                         dgemm_("N", "N", m, n, m, &c_b79, &work[iu], &ldwrku,
4170                                 &vt[vt_offset], ldvt, &c_b57, &a[a_offset],
4171                                 lda);
4172
4173 /*                    Copy right singular vectors of A from A to VT */
4174
4175                         dlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
4176                                 ldvt);
4177
4178 /*                    Copy left singular vectors of A from WORK(IR) to A */
4179
4180                         dlacpy_("F", m, m, &work[ir], &ldwrkr, &a[a_offset],
4181                                 lda);
4182
4183                     } else {
4184
4185 /*                    Insufficient workspace for a fast algorithm */
4186
4187                         itau = 1;
4188                         iwork = itau + *m;
4189
4190 /*                    Compute A=L*Q, copying result to VT */
4191 /*                    (Workspace: need 2*M, prefer M + M*NB) */
4192
4193                         i__2 = *lwork - iwork + 1;
4194                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
4195                                 iwork], &i__2, &ierr);
4196                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4197                                 ldvt);
4198
4199 /*                    Generate Q in VT */
4200 /*                    (Workspace: need M + N, prefer M + N*NB) */
4201
4202                         i__2 = *lwork - iwork + 1;
4203                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4204                                 work[iwork], &i__2, &ierr);
4205                         ie = itau;
4206                         itauq = ie + *m;
4207                         itaup = itauq + *m;
4208                         iwork = itaup + *m;
4209
4210 /*                    Zero out above L in A */
4211
4212                         i__2 = *m - 1;
4213                         i__3 = *m - 1;
4214                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &a[(a_dim1
4215                                 << 1) + 1], lda);
4216
4217 /*                    Bidiagonalize L in A */
4218 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
4219
4220                         i__2 = *lwork - iwork + 1;
4221                         dgebrd_(m, m, &a[a_offset], lda, &s[1], &work[ie], &
4222                                 work[itauq], &work[itaup], &work[iwork], &
4223                                 i__2, &ierr);
4224
4225 /*                    Multiply right bidiagonalizing vectors in A by Q */
4226 /*                    in VT */
4227 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
4228
4229                         i__2 = *lwork - iwork + 1;
4230                         dormbr_("P", "L", "T", m, n, m, &a[a_offset], lda, &
4231                                 work[itaup], &vt[vt_offset], ldvt, &work[
4232                                 iwork], &i__2, &ierr);
4233
4234 /*                    Generate left bidiagonalizing vectors in A */
4235 /*                    (Workspace: need 4*M, prefer 3*M + M*NB) */
4236
4237                         i__2 = *lwork - iwork + 1;
4238                         dorgbr_("Q", m, m, m, &a[a_offset], lda, &work[itauq],
4239                                  &work[iwork], &i__2, &ierr);
4240                         iwork = ie + *m;
4241
4242 /*                    Perform bidiagonal QR iteration, computing left */
4243 /*                    singular vectors of A in A and computing right */
4244 /*                    singular vectors of A in VT */
4245 /*                    (Workspace: need BDSPAC) */
4246
4247                         dbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
4248                                 vt_offset], ldvt, &a[a_offset], lda, dum, &
4249                                 c__1, &work[iwork], info);
4250
4251                     }
4252
4253                 } else if (wntuas) {
4254
4255 /*                 Path 9t(N much larger than M, JOBU='S' or 'A', */
4256 /*                         JOBVT='A') */
4257 /*                 N right singular vectors to be computed in VT and */
4258 /*                 M left singular vectors to be computed in U */
4259
4260 /* Computing MAX */
4261                     i__2 = *n + *m, i__3 = *m << 2, i__2 = f2cmax(i__2,i__3);
4262                     if (*lwork >= *m * *m + f2cmax(i__2,bdspac)) {
4263
4264 /*                    Sufficient workspace for a fast algorithm */
4265
4266                         iu = 1;
4267                         if (*lwork >= wrkbl + *lda * *m) {
4268
4269 /*                       WORK(IU) is LDA by M */
4270
4271                             ldwrku = *lda;
4272                         } else {
4273
4274 /*                       WORK(IU) is M by M */
4275
4276                             ldwrku = *m;
4277                         }
4278                         itau = iu + ldwrku * *m;
4279                         iwork = itau + *m;
4280
4281 /*                    Compute A=L*Q, copying result to VT */
4282 /*                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) */
4283
4284                         i__2 = *lwork - iwork + 1;
4285                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
4286                                 iwork], &i__2, &ierr);
4287                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4288                                 ldvt);
4289
4290 /*                    Generate Q in VT */
4291 /*                    (Workspace: need M*M + M + N, prefer M*M + M + N*NB) */
4292
4293                         i__2 = *lwork - iwork + 1;
4294                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4295                                 work[iwork], &i__2, &ierr);
4296
4297 /*                    Copy L to WORK(IU), zeroing out above it */
4298
4299                         dlacpy_("L", m, m, &a[a_offset], lda, &work[iu], &
4300                                 ldwrku);
4301                         i__2 = *m - 1;
4302                         i__3 = *m - 1;
4303                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &work[iu +
4304                                 ldwrku], &ldwrku);
4305                         ie = itau;
4306                         itauq = ie + *m;
4307                         itaup = itauq + *m;
4308                         iwork = itaup + *m;
4309
4310 /*                    Bidiagonalize L in WORK(IU), copying result to U */
4311 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) */
4312
4313                         i__2 = *lwork - iwork + 1;
4314                         dgebrd_(m, m, &work[iu], &ldwrku, &s[1], &work[ie], &
4315                                 work[itauq], &work[itaup], &work[iwork], &
4316                                 i__2, &ierr);
4317                         dlacpy_("L", m, m, &work[iu], &ldwrku, &u[u_offset],
4318                                 ldu);
4319
4320 /*                    Generate right bidiagonalizing vectors in WORK(IU) */
4321 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) */
4322
4323                         i__2 = *lwork - iwork + 1;
4324                         dorgbr_("P", m, m, m, &work[iu], &ldwrku, &work[itaup]
4325                                 , &work[iwork], &i__2, &ierr);
4326
4327 /*                    Generate left bidiagonalizing vectors in U */
4328 /*                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) */
4329
4330                         i__2 = *lwork - iwork + 1;
4331                         dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
4332                                  &work[iwork], &i__2, &ierr);
4333                         iwork = ie + *m;
4334
4335 /*                    Perform bidiagonal QR iteration, computing left */
4336 /*                    singular vectors of L in U and computing right */
4337 /*                    singular vectors of L in WORK(IU) */
4338 /*                    (Workspace: need M*M + BDSPAC) */
4339
4340                         dbdsqr_("U", m, m, m, &c__0, &s[1], &work[ie], &work[
4341                                 iu], &ldwrku, &u[u_offset], ldu, dum, &c__1, &
4342                                 work[iwork], info);
4343
4344 /*                    Multiply right singular vectors of L in WORK(IU) by */
4345 /*                    Q in VT, storing result in A */
4346 /*                    (Workspace: need M*M) */
4347
4348                         dgemm_("N", "N", m, n, m, &c_b79, &work[iu], &ldwrku,
4349                                 &vt[vt_offset], ldvt, &c_b57, &a[a_offset],
4350                                 lda);
4351
4352 /*                    Copy right singular vectors of A from A to VT */
4353
4354                         dlacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset],
4355                                 ldvt);
4356
4357                     } else {
4358
4359 /*                    Insufficient workspace for a fast algorithm */
4360
4361                         itau = 1;
4362                         iwork = itau + *m;
4363
4364 /*                    Compute A=L*Q, copying result to VT */
4365 /*                    (Workspace: need 2*M, prefer M + M*NB) */
4366
4367                         i__2 = *lwork - iwork + 1;
4368                         dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[
4369                                 iwork], &i__2, &ierr);
4370                         dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset],
4371                                 ldvt);
4372
4373 /*                    Generate Q in VT */
4374 /*                    (Workspace: need M + N, prefer M + N*NB) */
4375
4376                         i__2 = *lwork - iwork + 1;
4377                         dorglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &
4378                                 work[iwork], &i__2, &ierr);
4379
4380 /*                    Copy L to U, zeroing out above it */
4381
4382                         dlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset],
4383                                 ldu);
4384                         i__2 = *m - 1;
4385                         i__3 = *m - 1;
4386                         dlaset_("U", &i__2, &i__3, &c_b57, &c_b57, &u[(u_dim1
4387                                 << 1) + 1], ldu);
4388                         ie = itau;
4389                         itauq = ie + *m;
4390                         itaup = itauq + *m;
4391                         iwork = itaup + *m;
4392
4393 /*                    Bidiagonalize L in U */
4394 /*                    (Workspace: need 4*M, prefer 3*M + 2*M*NB) */
4395
4396                         i__2 = *lwork - iwork + 1;
4397                         dgebrd_(m, m, &u[u_offset], ldu, &s[1], &work[ie], &
4398                                 work[itauq], &work[itaup], &work[iwork], &
4399                                 i__2, &ierr);
4400
4401 /*                    Multiply right bidiagonalizing vectors in U by Q */
4402 /*                    in VT */
4403 /*                    (Workspace: need 3*M + N, prefer 3*M + N*NB) */
4404
4405                         i__2 = *lwork - iwork + 1;
4406                         dormbr_("P", "L", "T", m, n, m, &u[u_offset], ldu, &
4407                                 work[itaup], &vt[vt_offset], ldvt, &work[
4408                                 iwork], &i__2, &ierr);
4409
4410 /*                    Generate left bidiagonalizing vectors in U */
4411 /*                    (Workspace: need 4*M, prefer 3*M + M*NB) */
4412
4413                         i__2 = *lwork - iwork + 1;
4414                         dorgbr_("Q", m, m, m, &u[u_offset], ldu, &work[itauq],
4415                                  &work[iwork], &i__2, &ierr);
4416                         iwork = ie + *m;
4417
4418 /*                    Perform bidiagonal QR iteration, computing left */
4419 /*                    singular vectors of A in U and computing right */
4420 /*                    singular vectors of A in VT */
4421 /*                    (Workspace: need BDSPAC) */
4422
4423                         dbdsqr_("U", m, n, m, &c__0, &s[1], &work[ie], &vt[
4424                                 vt_offset], ldvt, &u[u_offset], ldu, dum, &
4425                                 c__1, &work[iwork], info);
4426
4427                     }
4428
4429                 }
4430
4431             }
4432
4433         } else {
4434
4435 /*           N .LT. MNTHR */
4436
4437 /*           Path 10t(N greater than M, but not much larger) */
4438 /*           Reduce to bidiagonal form without LQ decomposition */
4439
4440             ie = 1;
4441             itauq = ie + *m;
4442             itaup = itauq + *m;
4443             iwork = itaup + *m;
4444
4445 /*           Bidiagonalize A */
4446 /*           (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) */
4447
4448             i__2 = *lwork - iwork + 1;
4449             dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
4450                     work[itaup], &work[iwork], &i__2, &ierr);
4451             if (wntuas) {
4452
4453 /*              If left singular vectors desired in U, copy result to U */
4454 /*              and generate left bidiagonalizing vectors in U */
4455 /*              (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) */
4456
4457                 dlacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
4458                 i__2 = *lwork - iwork + 1;
4459                 dorgbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
4460                         iwork], &i__2, &ierr);
4461             }
4462             if (wntvas) {
4463
4464 /*              If right singular vectors desired in VT, copy result to */
4465 /*              VT and generate right bidiagonalizing vectors in VT */
4466 /*              (Workspace: need 3*M + NRVT, prefer 3*M + NRVT*NB) */
4467
4468                 dlacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
4469                 if (wntva) {
4470                     nrvt = *n;
4471                 }
4472                 if (wntvs) {
4473                     nrvt = *m;
4474                 }
4475                 i__2 = *lwork - iwork + 1;
4476                 dorgbr_("P", &nrvt, n, m, &vt[vt_offset], ldvt, &work[itaup],
4477                         &work[iwork], &i__2, &ierr);
4478             }
4479             if (wntuo) {
4480
4481 /*              If left singular vectors desired in A, generate left */
4482 /*              bidiagonalizing vectors in A */
4483 /*              (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) */
4484
4485                 i__2 = *lwork - iwork + 1;
4486                 dorgbr_("Q", m, m, n, &a[a_offset], lda, &work[itauq], &work[
4487                         iwork], &i__2, &ierr);
4488             }
4489             if (wntvo) {
4490
4491 /*              If right singular vectors desired in A, generate right */
4492 /*              bidiagonalizing vectors in A */
4493 /*              (Workspace: need 4*M, prefer 3*M + M*NB) */
4494
4495                 i__2 = *lwork - iwork + 1;
4496                 dorgbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
4497                         iwork], &i__2, &ierr);
4498             }
4499             iwork = ie + *m;
4500             if (wntuas || wntuo) {
4501                 nru = *m;
4502             }
4503             if (wntun) {
4504                 nru = 0;
4505             }
4506             if (wntvas || wntvo) {
4507                 ncvt = *n;
4508             }
4509             if (wntvn) {
4510                 ncvt = 0;
4511             }
4512             if (! wntuo && ! wntvo) {
4513
4514 /*              Perform bidiagonal QR iteration, if desired, computing */
4515 /*              left singular vectors in U and computing right singular */
4516 /*              vectors in VT */
4517 /*              (Workspace: need BDSPAC) */
4518
4519                 dbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
4520                         vt_offset], ldvt, &u[u_offset], ldu, dum, &c__1, &
4521                         work[iwork], info);
4522             } else if (! wntuo && wntvo) {
4523
4524 /*              Perform bidiagonal QR iteration, if desired, computing */
4525 /*              left singular vectors in U and computing right singular */
4526 /*              vectors in A */
4527 /*              (Workspace: need BDSPAC) */
4528
4529                 dbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &a[
4530                         a_offset], lda, &u[u_offset], ldu, dum, &c__1, &work[
4531                         iwork], info);
4532             } else {
4533
4534 /*              Perform bidiagonal QR iteration, if desired, computing */
4535 /*              left singular vectors in A and computing right singular */
4536 /*              vectors in VT */
4537 /*              (Workspace: need BDSPAC) */
4538
4539                 dbdsqr_("L", m, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[
4540                         vt_offset], ldvt, &a[a_offset], lda, dum, &c__1, &
4541                         work[iwork], info);
4542             }
4543
4544         }
4545
4546     }
4547
4548 /*     If DBDSQR failed to converge, copy unconverged superdiagonals */
4549 /*     to WORK( 2:MINMN ) */
4550
4551     if (*info != 0) {
4552         if (ie > 2) {
4553             i__2 = minmn - 1;
4554             for (i__ = 1; i__ <= i__2; ++i__) {
4555                 work[i__ + 1] = work[i__ + ie - 1];
4556 /* L50: */
4557             }
4558         }
4559         if (ie < 2) {
4560             for (i__ = minmn - 1; i__ >= 1; --i__) {
4561                 work[i__ + 1] = work[i__ + ie - 1];
4562 /* L60: */
4563             }
4564         }
4565     }
4566
4567 /*     Undo scaling if necessary */
4568
4569     if (iscl == 1) {
4570         if (anrm > bignum) {
4571             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
4572                     minmn, &ierr);
4573         }
4574         if (*info != 0 && anrm > bignum) {
4575             i__2 = minmn - 1;
4576             dlascl_("G", &c__0, &c__0, &bignum, &anrm, &i__2, &c__1, &work[2],
4577                      &minmn, &ierr);
4578         }
4579         if (anrm < smlnum) {
4580             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
4581                     minmn, &ierr);
4582         }
4583         if (*info != 0 && anrm < smlnum) {
4584             i__2 = minmn - 1;
4585             dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__2, &c__1, &work[2],
4586                      &minmn, &ierr);
4587         }
4588     }
4589
4590 /*     Return optimal workspace in WORK(1) */
4591
4592     work[1] = (doublereal) maxwrk;
4593
4594     return 0;
4595
4596 /*     End of DGESVD */
4597
4598 } /* dgesvd_ */
4599