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