14 typedef long long BLASLONG;
15 typedef unsigned long long BLASULONG;
17 typedef long BLASLONG;
18 typedef unsigned long BLASULONG;
22 typedef BLASLONG blasint;
24 #define blasabs(x) llabs(x)
26 #define blasabs(x) labs(x)
30 #define blasabs(x) abs(x)
33 typedef blasint integer;
35 typedef unsigned int uinteger;
36 typedef char *address;
37 typedef short int shortint;
39 typedef double doublereal;
40 typedef struct { real r, i; } complex;
41 typedef struct { doublereal r, i; } doublecomplex;
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;}
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;}
53 #define pCf(z) (*_pCf(z))
54 #define pCd(z) (*_pCd(z))
56 typedef short int shortlogical;
57 typedef char logical1;
58 typedef char integer1;
63 /* Extern is for use with -E */
74 /*external read, write*/
83 /*internal read, write*/
113 /*rewind, backspace, endfile*/
125 ftnint *inex; /*parameters in standard's order*/
151 union Multitype { /* for multiple entry points */
162 typedef union Multitype Multitype;
164 struct Vardesc { /* for Namelist */
170 typedef struct Vardesc Vardesc;
177 typedef struct Namelist Namelist;
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)))
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)); }
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]);}
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);}
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)}
261 /* procedure parameter types for -A and -C++ */
263 #define F2C_proc_par_types 1
265 typedef logical (*L_fp)(...);
267 typedef logical (*L_fp)();
270 static float spow_ui(float x, integer n) {
271 float pow=1.0; unsigned long int u;
273 if(n < 0) n = -n, x = 1/x;
282 static double dpow_ui(double x, integer n) {
283 double pow=1.0; unsigned long int u;
285 if(n < 0) n = -n, x = 1/x;
295 static _Fcomplex cpow_ui(complex x, integer n) {
296 complex pow={1.0,0.0}; unsigned long int u;
298 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
300 if(u & 01) pow.r *= x.r, pow.i *= x.i;
301 if(u >>= 1) x.r *= x.r, x.i *= x.i;
305 _Fcomplex p={pow.r, pow.i};
309 static _Complex float cpow_ui(_Complex float x, integer n) {
310 _Complex float pow=1.0; unsigned long int u;
312 if(n < 0) n = -n, x = 1/x;
323 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
324 _Dcomplex pow={1.0,0.0}; unsigned long int u;
326 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
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];
333 _Dcomplex p = {pow._Val[0], pow._Val[1]};
337 static _Complex double zpow_ui(_Complex double x, integer n) {
338 _Complex double pow=1.0; unsigned long int u;
340 if(n < 0) n = -n, x = 1/x;
350 static integer pow_ii(integer x, integer n) {
351 integer pow; unsigned long int u;
353 if (n == 0 || x == 1) pow = 1;
354 else if (x != -1) pow = x == 0 ? 1/x : 0;
357 if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
367 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
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];
374 static integer smaxloc_(float *w, integer s, integer e, integer *n)
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];
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;
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];
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];
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]);
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]);
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;
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];
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];
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]);
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]);
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;
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];
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];
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]);
467 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
468 zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
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;
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];
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];
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]);
498 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
499 zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
505 /* -- translated by f2c (version 20000121).
506 You must link the resulting object file with the libraries:
507 -lf2c -lm (in that order)
513 /* Table of constant values */
515 static complex c_b1 = {0.f,0.f};
516 static complex c_b2 = {1.f,0.f};
517 static integer c_n1 = -1;
518 static integer c__0 = 0;
519 static integer c__1 = 1;
521 /* > \brief \b CGESDD */
523 /* =========== DOCUMENTATION =========== */
525 /* Online html documentation available at */
526 /* http://www.netlib.org/lapack/explore-html/ */
529 /* > Download CGESDD + dependencies */
530 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesdd.
533 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesdd.
536 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesdd.
544 /* SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, */
545 /* WORK, LWORK, RWORK, IWORK, INFO ) */
548 /* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N */
549 /* INTEGER IWORK( * ) */
550 /* REAL RWORK( * ), S( * ) */
551 /* COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), */
555 /* > \par Purpose: */
560 /* > CGESDD computes the singular value decomposition (SVD) of a complex */
561 /* > M-by-N matrix A, optionally computing the left and/or right singular */
562 /* > vectors, by using divide-and-conquer method. The SVD is written */
564 /* > A = U * SIGMA * conjugate-transpose(V) */
566 /* > where SIGMA is an M-by-N matrix which is zero except for its */
567 /* > f2cmin(m,n) diagonal elements, U is an M-by-M unitary matrix, and */
568 /* > V is an N-by-N unitary matrix. The diagonal elements of SIGMA */
569 /* > are the singular values of A; they are real and non-negative, and */
570 /* > are returned in descending order. The first f2cmin(m,n) columns of */
571 /* > U and V are the left and right singular vectors of A. */
573 /* > Note that the routine returns VT = V**H, not V. */
575 /* > The divide and conquer algorithm makes very mild assumptions about */
576 /* > floating point arithmetic. It will work on machines with a guard */
577 /* > digit in add/subtract, or on those binary machines without guard */
578 /* > digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */
579 /* > Cray-2. It could conceivably fail on hexadecimal or decimal machines */
580 /* > without guard digits, but we know of none. */
586 /* > \param[in] JOBZ */
588 /* > JOBZ is CHARACTER*1 */
589 /* > Specifies options for computing all or part of the matrix U: */
590 /* > = 'A': all M columns of U and all N rows of V**H are */
591 /* > returned in the arrays U and VT; */
592 /* > = 'S': the first f2cmin(M,N) columns of U and the first */
593 /* > f2cmin(M,N) rows of V**H are returned in the arrays U */
595 /* > = 'O': If M >= N, the first N columns of U are overwritten */
596 /* > in the array A and all rows of V**H are returned in */
597 /* > the array VT; */
598 /* > otherwise, all columns of U are returned in the */
599 /* > array U and the first M rows of V**H are overwritten */
600 /* > in the array A; */
601 /* > = 'N': no columns of U or rows of V**H are computed. */
607 /* > The number of rows of the input matrix A. M >= 0. */
613 /* > The number of columns of the input matrix A. N >= 0. */
616 /* > \param[in,out] A */
618 /* > A is COMPLEX array, dimension (LDA,N) */
619 /* > On entry, the M-by-N matrix A. */
621 /* > if JOBZ = 'O', A is overwritten with the first N columns */
622 /* > of U (the left singular vectors, stored */
623 /* > columnwise) if M >= N; */
624 /* > A is overwritten with the first M rows */
625 /* > of V**H (the right singular vectors, stored */
626 /* > rowwise) otherwise. */
627 /* > if JOBZ .ne. 'O', the contents of A are destroyed. */
630 /* > \param[in] LDA */
632 /* > LDA is INTEGER */
633 /* > The leading dimension of the array A. LDA >= f2cmax(1,M). */
636 /* > \param[out] S */
638 /* > S is REAL array, dimension (f2cmin(M,N)) */
639 /* > The singular values of A, sorted so that S(i) >= S(i+1). */
642 /* > \param[out] U */
644 /* > U is COMPLEX array, dimension (LDU,UCOL) */
645 /* > UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; */
646 /* > UCOL = f2cmin(M,N) if JOBZ = 'S'. */
647 /* > If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M */
648 /* > unitary matrix U; */
649 /* > if JOBZ = 'S', U contains the first f2cmin(M,N) columns of U */
650 /* > (the left singular vectors, stored columnwise); */
651 /* > if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. */
654 /* > \param[in] LDU */
656 /* > LDU is INTEGER */
657 /* > The leading dimension of the array U. LDU >= 1; */
658 /* > if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. */
661 /* > \param[out] VT */
663 /* > VT is COMPLEX array, dimension (LDVT,N) */
664 /* > If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the */
665 /* > N-by-N unitary matrix V**H; */
666 /* > if JOBZ = 'S', VT contains the first f2cmin(M,N) rows of */
667 /* > V**H (the right singular vectors, stored rowwise); */
668 /* > if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. */
671 /* > \param[in] LDVT */
673 /* > LDVT is INTEGER */
674 /* > The leading dimension of the array VT. LDVT >= 1; */
675 /* > if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; */
676 /* > if JOBZ = 'S', LDVT >= f2cmin(M,N). */
679 /* > \param[out] WORK */
681 /* > WORK is COMPLEX array, dimension (MAX(1,LWORK)) */
682 /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
685 /* > \param[in] LWORK */
687 /* > LWORK is INTEGER */
688 /* > The dimension of the array WORK. LWORK >= 1. */
689 /* > If LWORK = -1, a workspace query is assumed. The optimal */
690 /* > size for the WORK array is calculated and stored in WORK(1), */
691 /* > and no other work except argument checking is performed. */
693 /* > Let mx = f2cmax(M,N) and mn = f2cmin(M,N). */
694 /* > If JOBZ = 'N', LWORK >= 2*mn + mx. */
695 /* > If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx. */
696 /* > If JOBZ = 'S', LWORK >= mn*mn + 3*mn. */
697 /* > If JOBZ = 'A', LWORK >= mn*mn + 2*mn + mx. */
698 /* > These are not tight minimums in all cases; see comments inside code. */
699 /* > For good performance, LWORK should generally be larger; */
700 /* > a query is recommended. */
703 /* > \param[out] RWORK */
705 /* > RWORK is REAL array, dimension (MAX(1,LRWORK)) */
706 /* > Let mx = f2cmax(M,N) and mn = f2cmin(M,N). */
707 /* > If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); */
708 /* > else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn; */
709 /* > else LRWORK >= f2cmax( 5*mn*mn + 5*mn, */
710 /* > 2*mx*mn + 2*mn*mn + mn ). */
713 /* > \param[out] IWORK */
715 /* > IWORK is INTEGER array, dimension (8*f2cmin(M,N)) */
718 /* > \param[out] INFO */
720 /* > INFO is INTEGER */
721 /* > = 0: successful exit. */
722 /* > < 0: if INFO = -i, the i-th argument had an illegal value. */
723 /* > > 0: The updating process of SBDSDC did not converge. */
729 /* > \author Univ. of Tennessee */
730 /* > \author Univ. of California Berkeley */
731 /* > \author Univ. of Colorado Denver */
732 /* > \author NAG Ltd. */
734 /* > \date June 2016 */
736 /* > \ingroup complexGEsing */
738 /* > \par Contributors: */
739 /* ================== */
741 /* > Ming Gu and Huan Ren, Computer Science Division, University of */
742 /* > California at Berkeley, USA */
744 /* ===================================================================== */
745 /* Subroutine */ int cgesdd_(char *jobz, integer *m, integer *n, complex *a,
746 integer *lda, real *s, complex *u, integer *ldu, complex *vt, integer
747 *ldvt, complex *work, integer *lwork, real *rwork, integer *iwork,
750 /* System generated locals */
751 integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
754 /* Local variables */
755 integer lwork_cunglq_mn__, lwork_cunglq_nn__, lwork_cungqr_mm__,
758 integer iscl, lwork_cunmbr_prc_mm__, lwork_cunmbr_prc_mn__,
759 lwork_cunmbr_prc_nn__;
761 integer ierr, itau, lwork_cunmbr_qln_mm__, lwork_cunmbr_qln_mn__,
762 lwork_cunmbr_qln_nn__, idum[1], irvt, i__;
763 extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *,
764 integer *, complex *, complex *, integer *, complex *, integer *,
765 complex *, complex *, integer *);
766 extern logical lsame_(char *, char *);
767 integer chunk, minmn, wrkbl, itaup, itauq;
770 extern /* Subroutine */ int clacp2_(char *, integer *, integer *, real *,
771 integer *, complex *, integer *);
772 logical wntqn, wntqo, wntqs;
773 integer mnthr1, mnthr2, ie, lwork_cungbr_p_mn__, il, lwork_cungbr_p_nn__,
774 lwork_cungbr_q_mn__, lwork_cungbr_q_mm__;
775 extern /* Subroutine */ int cgebrd_(integer *, integer *, complex *,
776 integer *, real *, real *, complex *, complex *, complex *,
777 integer *, integer *);
779 extern real clange_(char *, integer *, integer *, complex *, integer *,
782 extern /* Subroutine */ int cgelqf_(integer *, integer *, complex *,
783 integer *, complex *, complex *, integer *, integer *), clacrm_(
784 integer *, integer *, complex *, integer *, real *, integer *,
785 complex *, integer *, real *), clarcm_(integer *, integer *, real
786 *, integer *, complex *, integer *, complex *, integer *, real *),
787 clascl_(char *, integer *, integer *, real *, real *, integer *,
788 integer *, complex *, integer *, integer *), sbdsdc_(char
789 *, char *, integer *, real *, real *, real *, integer *, real *,
790 integer *, real *, integer *, real *, integer *, integer *), cgeqrf_(integer *, integer *, complex *, integer
791 *, complex *, complex *, integer *, integer *);
792 extern real slamch_(char *);
793 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex
794 *, integer *, complex *, integer *), claset_(char *,
795 integer *, integer *, complex *, complex *, complex *, integer *), xerbla_(char *, integer *, ftnlen), cungbr_(char *,
796 integer *, integer *, integer *, complex *, integer *, complex *,
797 complex *, integer *, integer *);
799 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
800 real *, integer *, integer *, real *, integer *, integer *), cunmbr_(char *, char *, char *, integer *, integer *,
801 integer *, complex *, integer *, complex *, complex *, integer *,
802 complex *, integer *, integer *), cunglq_(
803 integer *, integer *, integer *, complex *, integer *, complex *,
804 complex *, integer *, integer *);
805 extern logical sisnan_(real *);
807 extern /* Subroutine */ int cungqr_(integer *, integer *, integer *,
808 complex *, integer *, complex *, complex *, integer *, integer *);
809 integer ldwrkr, minwrk, ldwrku, maxwrk, ldwkvt;
811 logical wntqas, lquery;
814 integer iru, ivt, lwork_cgebrd_mm__, lwork_cgebrd_mn__, lwork_cgebrd_nn__,
815 lwork_cgelqf_mn__, lwork_cgeqrf_mn__;
818 /* -- LAPACK driver routine (version 3.7.0) -- */
819 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
820 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
824 /* ===================================================================== */
827 /* Test the input arguments */
829 /* Parameter adjustments */
831 a_offset = 1 + a_dim1 * 1;
835 u_offset = 1 + u_dim1 * 1;
838 vt_offset = 1 + vt_dim1 * 1;
846 minmn = f2cmin(*m,*n);
847 mnthr1 = (integer) (minmn * 17.f / 9.f);
848 mnthr2 = (integer) (minmn * 5.f / 3.f);
849 wntqa = lsame_(jobz, "A");
850 wntqs = lsame_(jobz, "S");
851 wntqas = wntqa || wntqs;
852 wntqo = lsame_(jobz, "O");
853 wntqn = lsame_(jobz, "N");
854 lquery = *lwork == -1;
858 if (! (wntqa || wntqs || wntqo || wntqn)) {
864 } else if (*lda < f2cmax(1,*m)) {
866 } else if (*ldu < 1 || wntqas && *ldu < *m || wntqo && *m < *n && *ldu < *
869 } else if (*ldvt < 1 || wntqa && *ldvt < *n || wntqs && *ldvt < minmn ||
870 wntqo && *m >= *n && *ldvt < *n) {
874 /* Compute workspace */
875 /* Note: Comments in the code beginning "Workspace:" describe the */
876 /* minimal amount of workspace allocated at that point in the code, */
877 /* as well as the preferred amount for good performance. */
878 /* CWorkspace refers to complex workspace, and RWorkspace to */
879 /* real workspace. NB refers to the optimal block size for the */
880 /* immediately following subroutine, as returned by ILAENV.) */
885 if (*m >= *n && minmn > 0) {
887 /* There is no complex work space needed for bidiagonal SVD */
888 /* The real work space needed for bidiagonal SVD (sbdsdc) is */
889 /* BDSPAC = 3*N*N + 4*N for singular values and vectors; */
890 /* BDSPAC = 4*N for singular values only; */
891 /* not including e, RU, and RVT matrices. */
893 /* Compute space preferred for each routine */
894 cgebrd_(m, n, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
895 lwork_cgebrd_mn__ = (integer) cdum[0].r;
897 cgebrd_(n, n, cdum, n, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
898 lwork_cgebrd_nn__ = (integer) cdum[0].r;
900 cgeqrf_(m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
901 lwork_cgeqrf_mn__ = (integer) cdum[0].r;
903 cungbr_("P", n, n, n, cdum, n, cdum, cdum, &c_n1, &ierr);
904 lwork_cungbr_p_nn__ = (integer) cdum[0].r;
906 cungbr_("Q", m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
907 lwork_cungbr_q_mm__ = (integer) cdum[0].r;
909 cungbr_("Q", m, n, n, cdum, m, cdum, cdum, &c_n1, &ierr);
910 lwork_cungbr_q_mn__ = (integer) cdum[0].r;
912 cungqr_(m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
913 lwork_cungqr_mm__ = (integer) cdum[0].r;
915 cungqr_(m, n, n, cdum, m, cdum, cdum, &c_n1, &ierr);
916 lwork_cungqr_mn__ = (integer) cdum[0].r;
918 cunmbr_("P", "R", "C", n, n, n, cdum, n, cdum, cdum, n, cdum, &
920 lwork_cunmbr_prc_nn__ = (integer) cdum[0].r;
922 cunmbr_("Q", "L", "N", m, m, n, cdum, m, cdum, cdum, m, cdum, &
924 lwork_cunmbr_qln_mm__ = (integer) cdum[0].r;
926 cunmbr_("Q", "L", "N", m, n, n, cdum, m, cdum, cdum, m, cdum, &
928 lwork_cunmbr_qln_mn__ = (integer) cdum[0].r;
930 cunmbr_("Q", "L", "N", n, n, n, cdum, n, cdum, cdum, n, cdum, &
932 lwork_cunmbr_qln_nn__ = (integer) cdum[0].r;
937 /* Path 1 (M >> N, JOBZ='N') */
939 maxwrk = *n + lwork_cgeqrf_mn__;
941 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cgebrd_nn__;
942 maxwrk = f2cmax(i__1,i__2);
946 /* Path 2 (M >> N, JOBZ='O') */
948 wrkbl = *n + lwork_cgeqrf_mn__;
950 i__1 = wrkbl, i__2 = *n + lwork_cungqr_mn__;
951 wrkbl = f2cmax(i__1,i__2);
953 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cgebrd_nn__;
954 wrkbl = f2cmax(i__1,i__2);
956 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_qln_nn__;
957 wrkbl = f2cmax(i__1,i__2);
959 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
960 wrkbl = f2cmax(i__1,i__2);
961 maxwrk = *m * *n + *n * *n + wrkbl;
962 minwrk = (*n << 1) * *n + *n * 3;
965 /* Path 3 (M >> N, JOBZ='S') */
967 wrkbl = *n + lwork_cgeqrf_mn__;
969 i__1 = wrkbl, i__2 = *n + lwork_cungqr_mn__;
970 wrkbl = f2cmax(i__1,i__2);
972 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cgebrd_nn__;
973 wrkbl = f2cmax(i__1,i__2);
975 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_qln_nn__;
976 wrkbl = f2cmax(i__1,i__2);
978 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
979 wrkbl = f2cmax(i__1,i__2);
980 maxwrk = *n * *n + wrkbl;
981 minwrk = *n * *n + *n * 3;
984 /* Path 4 (M >> N, JOBZ='A') */
986 wrkbl = *n + lwork_cgeqrf_mn__;
988 i__1 = wrkbl, i__2 = *n + lwork_cungqr_mm__;
989 wrkbl = f2cmax(i__1,i__2);
991 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cgebrd_nn__;
992 wrkbl = f2cmax(i__1,i__2);
994 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_qln_nn__;
995 wrkbl = f2cmax(i__1,i__2);
997 i__1 = wrkbl, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
998 wrkbl = f2cmax(i__1,i__2);
999 maxwrk = *n * *n + wrkbl;
1001 i__1 = *n * 3, i__2 = *n + *m;
1002 minwrk = *n * *n + f2cmax(i__1,i__2);
1004 } else if (*m >= mnthr2) {
1006 /* Path 5 (M >> N, but not as much as MNTHR1) */
1008 maxwrk = (*n << 1) + lwork_cgebrd_mn__;
1009 minwrk = (*n << 1) + *m;
1011 /* Path 5o (M >> N, JOBZ='O') */
1013 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_p_nn__;
1014 maxwrk = f2cmax(i__1,i__2);
1016 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_q_mn__;
1017 maxwrk = f2cmax(i__1,i__2);
1021 /* Path 5s (M >> N, JOBZ='S') */
1023 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_p_nn__;
1024 maxwrk = f2cmax(i__1,i__2);
1026 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_q_mn__;
1027 maxwrk = f2cmax(i__1,i__2);
1029 /* Path 5a (M >> N, JOBZ='A') */
1031 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_p_nn__;
1032 maxwrk = f2cmax(i__1,i__2);
1034 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cungbr_q_mm__;
1035 maxwrk = f2cmax(i__1,i__2);
1039 /* Path 6 (M >= N, but not much larger) */
1041 maxwrk = (*n << 1) + lwork_cgebrd_mn__;
1042 minwrk = (*n << 1) + *m;
1044 /* Path 6o (M >= N, JOBZ='O') */
1046 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
1047 maxwrk = f2cmax(i__1,i__2);
1049 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_qln_mn__;
1050 maxwrk = f2cmax(i__1,i__2);
1054 /* Path 6s (M >= N, JOBZ='S') */
1056 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_qln_mn__;
1057 maxwrk = f2cmax(i__1,i__2);
1059 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
1060 maxwrk = f2cmax(i__1,i__2);
1062 /* Path 6a (M >= N, JOBZ='A') */
1064 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_qln_mm__;
1065 maxwrk = f2cmax(i__1,i__2);
1067 i__1 = maxwrk, i__2 = (*n << 1) + lwork_cunmbr_prc_nn__;
1068 maxwrk = f2cmax(i__1,i__2);
1071 } else if (minmn > 0) {
1073 /* There is no complex work space needed for bidiagonal SVD */
1074 /* The real work space needed for bidiagonal SVD (sbdsdc) is */
1075 /* BDSPAC = 3*M*M + 4*M for singular values and vectors; */
1076 /* BDSPAC = 4*M for singular values only; */
1077 /* not including e, RU, and RVT matrices. */
1079 /* Compute space preferred for each routine */
1080 cgebrd_(m, n, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
1081 lwork_cgebrd_mn__ = (integer) cdum[0].r;
1083 cgebrd_(m, m, cdum, m, dum, dum, cdum, cdum, cdum, &c_n1, &ierr);
1084 lwork_cgebrd_mm__ = (integer) cdum[0].r;
1086 cgelqf_(m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
1087 lwork_cgelqf_mn__ = (integer) cdum[0].r;
1089 cungbr_("P", m, n, m, cdum, m, cdum, cdum, &c_n1, &ierr);
1090 lwork_cungbr_p_mn__ = (integer) cdum[0].r;
1092 cungbr_("P", n, n, m, cdum, n, cdum, cdum, &c_n1, &ierr);
1093 lwork_cungbr_p_nn__ = (integer) cdum[0].r;
1095 cungbr_("Q", m, m, n, cdum, m, cdum, cdum, &c_n1, &ierr);
1096 lwork_cungbr_q_mm__ = (integer) cdum[0].r;
1098 cunglq_(m, n, m, cdum, m, cdum, cdum, &c_n1, &ierr);
1099 lwork_cunglq_mn__ = (integer) cdum[0].r;
1101 cunglq_(n, n, m, cdum, n, cdum, cdum, &c_n1, &ierr);
1102 lwork_cunglq_nn__ = (integer) cdum[0].r;
1104 cunmbr_("P", "R", "C", m, m, m, cdum, m, cdum, cdum, m, cdum, &
1106 lwork_cunmbr_prc_mm__ = (integer) cdum[0].r;
1108 cunmbr_("P", "R", "C", m, n, m, cdum, m, cdum, cdum, m, cdum, &
1110 lwork_cunmbr_prc_mn__ = (integer) cdum[0].r;
1112 cunmbr_("P", "R", "C", n, n, m, cdum, n, cdum, cdum, n, cdum, &
1114 lwork_cunmbr_prc_nn__ = (integer) cdum[0].r;
1116 cunmbr_("Q", "L", "N", m, m, m, cdum, m, cdum, cdum, m, cdum, &
1118 lwork_cunmbr_qln_mm__ = (integer) cdum[0].r;
1123 /* Path 1t (N >> M, JOBZ='N') */
1125 maxwrk = *m + lwork_cgelqf_mn__;
1127 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cgebrd_mm__;
1128 maxwrk = f2cmax(i__1,i__2);
1132 /* Path 2t (N >> M, JOBZ='O') */
1134 wrkbl = *m + lwork_cgelqf_mn__;
1136 i__1 = wrkbl, i__2 = *m + lwork_cunglq_mn__;
1137 wrkbl = f2cmax(i__1,i__2);
1139 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cgebrd_mm__;
1140 wrkbl = f2cmax(i__1,i__2);
1142 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1143 wrkbl = f2cmax(i__1,i__2);
1145 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_prc_mm__;
1146 wrkbl = f2cmax(i__1,i__2);
1147 maxwrk = *m * *n + *m * *m + wrkbl;
1148 minwrk = (*m << 1) * *m + *m * 3;
1151 /* Path 3t (N >> M, JOBZ='S') */
1153 wrkbl = *m + lwork_cgelqf_mn__;
1155 i__1 = wrkbl, i__2 = *m + lwork_cunglq_mn__;
1156 wrkbl = f2cmax(i__1,i__2);
1158 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cgebrd_mm__;
1159 wrkbl = f2cmax(i__1,i__2);
1161 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1162 wrkbl = f2cmax(i__1,i__2);
1164 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_prc_mm__;
1165 wrkbl = f2cmax(i__1,i__2);
1166 maxwrk = *m * *m + wrkbl;
1167 minwrk = *m * *m + *m * 3;
1170 /* Path 4t (N >> M, JOBZ='A') */
1172 wrkbl = *m + lwork_cgelqf_mn__;
1174 i__1 = wrkbl, i__2 = *m + lwork_cunglq_nn__;
1175 wrkbl = f2cmax(i__1,i__2);
1177 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cgebrd_mm__;
1178 wrkbl = f2cmax(i__1,i__2);
1180 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1181 wrkbl = f2cmax(i__1,i__2);
1183 i__1 = wrkbl, i__2 = (*m << 1) + lwork_cunmbr_prc_mm__;
1184 wrkbl = f2cmax(i__1,i__2);
1185 maxwrk = *m * *m + wrkbl;
1187 i__1 = *m * 3, i__2 = *m + *n;
1188 minwrk = *m * *m + f2cmax(i__1,i__2);
1190 } else if (*n >= mnthr2) {
1192 /* Path 5t (N >> M, but not as much as MNTHR1) */
1194 maxwrk = (*m << 1) + lwork_cgebrd_mn__;
1195 minwrk = (*m << 1) + *n;
1197 /* Path 5to (N >> M, JOBZ='O') */
1199 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_q_mm__;
1200 maxwrk = f2cmax(i__1,i__2);
1202 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_p_mn__;
1203 maxwrk = f2cmax(i__1,i__2);
1207 /* Path 5ts (N >> M, JOBZ='S') */
1209 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_q_mm__;
1210 maxwrk = f2cmax(i__1,i__2);
1212 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_p_mn__;
1213 maxwrk = f2cmax(i__1,i__2);
1215 /* Path 5ta (N >> M, JOBZ='A') */
1217 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_q_mm__;
1218 maxwrk = f2cmax(i__1,i__2);
1220 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cungbr_p_nn__;
1221 maxwrk = f2cmax(i__1,i__2);
1225 /* Path 6t (N > M, but not much larger) */
1227 maxwrk = (*m << 1) + lwork_cgebrd_mn__;
1228 minwrk = (*m << 1) + *n;
1230 /* Path 6to (N > M, JOBZ='O') */
1232 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1233 maxwrk = f2cmax(i__1,i__2);
1235 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_prc_mn__;
1236 maxwrk = f2cmax(i__1,i__2);
1240 /* Path 6ts (N > M, JOBZ='S') */
1242 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1243 maxwrk = f2cmax(i__1,i__2);
1245 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_prc_mn__;
1246 maxwrk = f2cmax(i__1,i__2);
1248 /* Path 6ta (N > M, JOBZ='A') */
1250 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_qln_mm__;
1251 maxwrk = f2cmax(i__1,i__2);
1253 i__1 = maxwrk, i__2 = (*m << 1) + lwork_cunmbr_prc_nn__;
1254 maxwrk = f2cmax(i__1,i__2);
1258 maxwrk = f2cmax(maxwrk,minwrk);
1261 work[1].r = (real) maxwrk, work[1].i = 0.f;
1262 if (*lwork < minwrk && ! lquery) {
1269 xerbla_("CGESDD", &i__1, (ftnlen)6);
1271 } else if (lquery) {
1275 /* Quick return if possible */
1277 if (*m == 0 || *n == 0) {
1281 /* Get machine constants */
1284 smlnum = sqrt(slamch_("S")) / eps;
1285 bignum = 1.f / smlnum;
1287 /* Scale A if f2cmax element outside range [SMLNUM,BIGNUM] */
1289 anrm = clange_("M", m, n, &a[a_offset], lda, dum);
1290 if (sisnan_(&anrm)) {
1295 if (anrm > 0.f && anrm < smlnum) {
1297 clascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &
1299 } else if (anrm > bignum) {
1301 clascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &
1307 /* A has at least as many rows as columns. If A has sufficiently */
1308 /* more rows than columns, first reduce using the QR */
1309 /* decomposition (if sufficient workspace available) */
1315 /* Path 1 (M >> N, JOBZ='N') */
1316 /* No singular vectors to be computed */
1322 /* CWorkspace: need N [tau] + N [work] */
1323 /* CWorkspace: prefer N [tau] + N*NB [work] */
1324 /* RWorkspace: need 0 */
1326 i__1 = *lwork - nwork + 1;
1327 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1330 /* Zero out below R */
1334 claset_("L", &i__1, &i__2, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1340 /* Bidiagonalize R in A */
1341 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1342 /* CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work] */
1343 /* RWorkspace: need N [e] */
1345 i__1 = *lwork - nwork + 1;
1346 cgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1347 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1350 /* Perform bidiagonal SVD, compute singular values only */
1351 /* CWorkspace: need 0 */
1352 /* RWorkspace: need N [e] + BDSPAC */
1354 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1355 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1359 /* Path 2 (M >> N, JOBZ='O') */
1360 /* N left singular vectors to be overwritten on A and */
1361 /* N right singular vectors to be computed in VT */
1365 /* WORK(IU) is N by N */
1368 ir = iu + ldwrku * *n;
1369 if (*lwork >= *m * *n + *n * *n + *n * 3) {
1371 /* WORK(IR) is M by N */
1375 ldwrkr = (*lwork - *n * *n - *n * 3) / *n;
1377 itau = ir + ldwrkr * *n;
1381 /* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work] */
1382 /* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] */
1383 /* RWorkspace: need 0 */
1385 i__1 = *lwork - nwork + 1;
1386 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1389 /* Copy R to WORK( IR ), zeroing out below it */
1391 clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1394 claset_("L", &i__1, &i__2, &c_b1, &c_b1, &work[ir + 1], &
1397 /* Generate Q in A */
1398 /* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work] */
1399 /* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] */
1400 /* RWorkspace: need 0 */
1402 i__1 = *lwork - nwork + 1;
1403 cungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
1410 /* Bidiagonalize R in WORK(IR) */
1411 /* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] */
1412 /* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] */
1413 /* RWorkspace: need N [e] */
1415 i__1 = *lwork - nwork + 1;
1416 cgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
1417 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
1419 /* Perform bidiagonal SVD, computing left singular vectors */
1420 /* of R in WORK(IRU) and computing right singular vectors */
1421 /* of R in WORK(IRVT) */
1422 /* CWorkspace: need 0 */
1423 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1426 irvt = iru + *n * *n;
1427 nrwork = irvt + *n * *n;
1428 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1429 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1432 /* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1433 /* Overwrite WORK(IU) by the left singular vectors of R */
1434 /* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] */
1435 /* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1436 /* RWorkspace: need 0 */
1438 clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
1439 i__1 = *lwork - nwork + 1;
1440 cunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
1441 itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
1444 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
1445 /* Overwrite VT by the right singular vectors of R */
1446 /* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] */
1447 /* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1448 /* RWorkspace: need 0 */
1450 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1451 i__1 = *lwork - nwork + 1;
1452 cunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
1453 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1456 /* Multiply Q in A by left singular vectors of R in */
1457 /* WORK(IU), storing result in WORK(IR) and copying to A */
1458 /* CWorkspace: need N*N [U] + N*N [R] */
1459 /* CWorkspace: prefer N*N [U] + M*N [R] */
1460 /* RWorkspace: need 0 */
1464 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
1467 i__3 = *m - i__ + 1;
1468 chunk = f2cmin(i__3,ldwrkr);
1469 cgemm_("N", "N", &chunk, n, n, &c_b2, &a[i__ + a_dim1],
1470 lda, &work[iu], &ldwrku, &c_b1, &work[ir], &
1472 clacpy_("F", &chunk, n, &work[ir], &ldwrkr, &a[i__ +
1479 /* Path 3 (M >> N, JOBZ='S') */
1480 /* N left singular vectors to be computed in U and */
1481 /* N right singular vectors to be computed in VT */
1485 /* WORK(IR) is N by N */
1488 itau = ir + ldwrkr * *n;
1492 /* CWorkspace: need N*N [R] + N [tau] + N [work] */
1493 /* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] */
1494 /* RWorkspace: need 0 */
1496 i__2 = *lwork - nwork + 1;
1497 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1500 /* Copy R to WORK(IR), zeroing out below it */
1502 clacpy_("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
1505 claset_("L", &i__2, &i__1, &c_b1, &c_b1, &work[ir + 1], &
1508 /* Generate Q in A */
1509 /* CWorkspace: need N*N [R] + N [tau] + N [work] */
1510 /* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] */
1511 /* RWorkspace: need 0 */
1513 i__2 = *lwork - nwork + 1;
1514 cungqr_(m, n, n, &a[a_offset], lda, &work[itau], &work[nwork],
1521 /* Bidiagonalize R in WORK(IR) */
1522 /* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] */
1523 /* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] */
1524 /* RWorkspace: need N [e] */
1526 i__2 = *lwork - nwork + 1;
1527 cgebrd_(n, n, &work[ir], &ldwrkr, &s[1], &rwork[ie], &work[
1528 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1530 /* Perform bidiagonal SVD, computing left singular vectors */
1531 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1532 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1533 /* CWorkspace: need 0 */
1534 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1537 irvt = iru + *n * *n;
1538 nrwork = irvt + *n * *n;
1539 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1540 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1543 /* Copy real matrix RWORK(IRU) to complex matrix U */
1544 /* Overwrite U by left singular vectors of R */
1545 /* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] */
1546 /* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1547 /* RWorkspace: need 0 */
1549 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
1550 i__2 = *lwork - nwork + 1;
1551 cunmbr_("Q", "L", "N", n, n, n, &work[ir], &ldwrkr, &work[
1552 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
1554 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
1555 /* Overwrite VT by right singular vectors of R */
1556 /* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] */
1557 /* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] */
1558 /* RWorkspace: need 0 */
1560 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1561 i__2 = *lwork - nwork + 1;
1562 cunmbr_("P", "R", "C", n, n, n, &work[ir], &ldwrkr, &work[
1563 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1566 /* Multiply Q in A by left singular vectors of R in */
1567 /* WORK(IR), storing result in U */
1568 /* CWorkspace: need N*N [R] */
1569 /* RWorkspace: need 0 */
1571 clacpy_("F", n, n, &u[u_offset], ldu, &work[ir], &ldwrkr);
1572 cgemm_("N", "N", m, n, n, &c_b2, &a[a_offset], lda, &work[ir],
1573 &ldwrkr, &c_b1, &u[u_offset], ldu);
1577 /* Path 4 (M >> N, JOBZ='A') */
1578 /* M left singular vectors to be computed in U and */
1579 /* N right singular vectors to be computed in VT */
1583 /* WORK(IU) is N by N */
1586 itau = iu + ldwrku * *n;
1589 /* Compute A=Q*R, copying result to U */
1590 /* CWorkspace: need N*N [U] + N [tau] + N [work] */
1591 /* CWorkspace: prefer N*N [U] + N [tau] + N*NB [work] */
1592 /* RWorkspace: need 0 */
1594 i__2 = *lwork - nwork + 1;
1595 cgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
1597 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1599 /* Generate Q in U */
1600 /* CWorkspace: need N*N [U] + N [tau] + M [work] */
1601 /* CWorkspace: prefer N*N [U] + N [tau] + M*NB [work] */
1602 /* RWorkspace: need 0 */
1604 i__2 = *lwork - nwork + 1;
1605 cungqr_(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
1608 /* Produce R in A, zeroing out below it */
1612 claset_("L", &i__2, &i__1, &c_b1, &c_b1, &a[a_dim1 + 2], lda);
1618 /* Bidiagonalize R in A */
1619 /* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] */
1620 /* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work] */
1621 /* RWorkspace: need N [e] */
1623 i__2 = *lwork - nwork + 1;
1624 cgebrd_(n, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[
1625 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
1627 irvt = iru + *n * *n;
1628 nrwork = irvt + *n * *n;
1630 /* Perform bidiagonal SVD, computing left singular vectors */
1631 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1632 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1633 /* CWorkspace: need 0 */
1634 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1636 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1637 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1640 /* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1641 /* Overwrite WORK(IU) by left singular vectors of R */
1642 /* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] */
1643 /* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] */
1644 /* RWorkspace: need 0 */
1646 clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
1647 i__2 = *lwork - nwork + 1;
1648 cunmbr_("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
1649 itauq], &work[iu], &ldwrku, &work[nwork], &i__2, &
1652 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
1653 /* Overwrite VT by right singular vectors of R */
1654 /* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] */
1655 /* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] */
1656 /* RWorkspace: need 0 */
1658 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1659 i__2 = *lwork - nwork + 1;
1660 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
1661 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
1664 /* Multiply Q in U by left singular vectors of R in */
1665 /* WORK(IU), storing result in A */
1666 /* CWorkspace: need N*N [U] */
1667 /* RWorkspace: need 0 */
1669 cgemm_("N", "N", m, n, n, &c_b2, &u[u_offset], ldu, &work[iu],
1670 &ldwrku, &c_b1, &a[a_offset], lda);
1672 /* Copy left singular vectors of A from A to U */
1674 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1678 } else if (*m >= mnthr2) {
1680 /* MNTHR2 <= M < MNTHR1 */
1682 /* Path 5 (M >> N, but not as much as MNTHR1) */
1683 /* Reduce to bidiagonal form without QR decomposition, use */
1684 /* CUNGBR and matrix multiplication to compute singular vectors */
1692 /* Bidiagonalize A */
1693 /* CWorkspace: need 2*N [tauq, taup] + M [work] */
1694 /* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] */
1695 /* RWorkspace: need N [e] */
1697 i__2 = *lwork - nwork + 1;
1698 cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
1699 &work[itaup], &work[nwork], &i__2, &ierr);
1702 /* Path 5n (M >> N, JOBZ='N') */
1703 /* Compute singular values only */
1704 /* CWorkspace: need 0 */
1705 /* RWorkspace: need N [e] + BDSPAC */
1707 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1708 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1712 irvt = iru + *n * *n;
1713 nrwork = irvt + *n * *n;
1715 /* Path 5o (M >> N, JOBZ='O') */
1716 /* Copy A to VT, generate P**H */
1717 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1718 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1719 /* RWorkspace: need 0 */
1721 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1722 i__2 = *lwork - nwork + 1;
1723 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1724 work[nwork], &i__2, &ierr);
1726 /* Generate Q in A */
1727 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1728 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1729 /* RWorkspace: need 0 */
1731 i__2 = *lwork - nwork + 1;
1732 cungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[
1733 nwork], &i__2, &ierr);
1735 if (*lwork >= *m * *n + *n * 3) {
1737 /* WORK( IU ) is M by N */
1742 /* WORK(IU) is LDWRKU by N */
1744 ldwrku = (*lwork - *n * 3) / *n;
1746 nwork = iu + ldwrku * *n;
1748 /* Perform bidiagonal SVD, computing left singular vectors */
1749 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1750 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1751 /* CWorkspace: need 0 */
1752 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1754 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1755 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1758 /* Multiply real matrix RWORK(IRVT) by P**H in VT, */
1759 /* storing the result in WORK(IU), copying to VT */
1760 /* CWorkspace: need 2*N [tauq, taup] + N*N [U] */
1761 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1763 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &work[iu]
1764 , &ldwrku, &rwork[nrwork]);
1765 clacpy_("F", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);
1767 /* Multiply Q in A by real matrix RWORK(IRU), storing the */
1768 /* result in WORK(IU), copying to A */
1769 /* CWorkspace: need 2*N [tauq, taup] + N*N [U] */
1770 /* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] */
1771 /* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork] */
1772 /* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1777 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
1780 i__3 = *m - i__ + 1;
1781 chunk = f2cmin(i__3,ldwrku);
1782 clacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru], n,
1783 &work[iu], &ldwrku, &rwork[nrwork]);
1784 clacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
1791 /* Path 5s (M >> N, JOBZ='S') */
1792 /* Copy A to VT, generate P**H */
1793 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1794 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1795 /* RWorkspace: need 0 */
1797 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1798 i__1 = *lwork - nwork + 1;
1799 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1800 work[nwork], &i__1, &ierr);
1802 /* Copy A to U, generate Q */
1803 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1804 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1805 /* RWorkspace: need 0 */
1807 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1808 i__1 = *lwork - nwork + 1;
1809 cungbr_("Q", m, n, n, &u[u_offset], ldu, &work[itauq], &work[
1810 nwork], &i__1, &ierr);
1812 /* Perform bidiagonal SVD, computing left singular vectors */
1813 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1814 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1815 /* CWorkspace: need 0 */
1816 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1819 irvt = iru + *n * *n;
1820 nrwork = irvt + *n * *n;
1821 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1822 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1825 /* Multiply real matrix RWORK(IRVT) by P**H in VT, */
1826 /* storing the result in A, copying to VT */
1827 /* CWorkspace: need 0 */
1828 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1830 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
1831 a_offset], lda, &rwork[nrwork]);
1832 clacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1834 /* Multiply Q in U by real matrix RWORK(IRU), storing the */
1835 /* result in A, copying to U */
1836 /* CWorkspace: need 0 */
1837 /* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1840 clacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset],
1841 lda, &rwork[nrwork]);
1842 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1845 /* Path 5a (M >> N, JOBZ='A') */
1846 /* Copy A to VT, generate P**H */
1847 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
1848 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
1849 /* RWorkspace: need 0 */
1851 clacpy_("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1852 i__1 = *lwork - nwork + 1;
1853 cungbr_("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &
1854 work[nwork], &i__1, &ierr);
1856 /* Copy A to U, generate Q */
1857 /* CWorkspace: need 2*N [tauq, taup] + M [work] */
1858 /* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] */
1859 /* RWorkspace: need 0 */
1861 clacpy_("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1862 i__1 = *lwork - nwork + 1;
1863 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
1864 nwork], &i__1, &ierr);
1866 /* Perform bidiagonal SVD, computing left singular vectors */
1867 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1868 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1869 /* CWorkspace: need 0 */
1870 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1873 irvt = iru + *n * *n;
1874 nrwork = irvt + *n * *n;
1875 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1876 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1879 /* Multiply real matrix RWORK(IRVT) by P**H in VT, */
1880 /* storing the result in A, copying to VT */
1881 /* CWorkspace: need 0 */
1882 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] */
1884 clarcm_(n, n, &rwork[irvt], n, &vt[vt_offset], ldvt, &a[
1885 a_offset], lda, &rwork[nrwork]);
1886 clacpy_("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
1888 /* Multiply Q in U by real matrix RWORK(IRU), storing the */
1889 /* result in A, copying to U */
1890 /* CWorkspace: need 0 */
1891 /* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
1894 clacrm_(m, n, &u[u_offset], ldu, &rwork[iru], n, &a[a_offset],
1895 lda, &rwork[nrwork]);
1896 clacpy_("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
1903 /* Path 6 (M >= N, but not much larger) */
1904 /* Reduce to bidiagonal form without QR decomposition */
1905 /* Use CUNMBR to compute singular vectors */
1913 /* Bidiagonalize A */
1914 /* CWorkspace: need 2*N [tauq, taup] + M [work] */
1915 /* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] */
1916 /* RWorkspace: need N [e] */
1918 i__1 = *lwork - nwork + 1;
1919 cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
1920 &work[itaup], &work[nwork], &i__1, &ierr);
1923 /* Path 6n (M >= N, JOBZ='N') */
1924 /* Compute singular values only */
1925 /* CWorkspace: need 0 */
1926 /* RWorkspace: need N [e] + BDSPAC */
1928 sbdsdc_("U", "N", n, &s[1], &rwork[ie], dum, &c__1, dum, &
1929 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
1933 irvt = iru + *n * *n;
1934 nrwork = irvt + *n * *n;
1935 if (*lwork >= *m * *n + *n * 3) {
1937 /* WORK( IU ) is M by N */
1942 /* WORK( IU ) is LDWRKU by N */
1944 ldwrku = (*lwork - *n * 3) / *n;
1946 nwork = iu + ldwrku * *n;
1948 /* Path 6o (M >= N, JOBZ='O') */
1949 /* Perform bidiagonal SVD, computing left singular vectors */
1950 /* of bidiagonal matrix in RWORK(IRU) and computing right */
1951 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
1952 /* CWorkspace: need 0 */
1953 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
1955 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
1956 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
1959 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
1960 /* Overwrite VT by right singular vectors of A */
1961 /* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work] */
1962 /* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] */
1963 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] */
1965 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
1966 i__1 = *lwork - nwork + 1;
1967 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
1968 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
1971 if (*lwork >= *m * *n + *n * 3) {
1974 /* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
1975 /* Overwrite WORK(IU) by left singular vectors of A, copying */
1977 /* CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work] */
1978 /* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work] */
1979 /* RWorkspace: need N [e] + N*N [RU] */
1981 claset_("F", m, n, &c_b1, &c_b1, &work[iu], &ldwrku);
1982 clacp2_("F", n, n, &rwork[iru], n, &work[iu], &ldwrku);
1983 i__1 = *lwork - nwork + 1;
1984 cunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
1985 itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
1987 clacpy_("F", m, n, &work[iu], &ldwrku, &a[a_offset], lda);
1991 /* Generate Q in A */
1992 /* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work] */
1993 /* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] */
1994 /* RWorkspace: need 0 */
1996 i__1 = *lwork - nwork + 1;
1997 cungbr_("Q", m, n, n, &a[a_offset], lda, &work[itauq], &
1998 work[nwork], &i__1, &ierr);
2000 /* Multiply Q in A by real matrix RWORK(IRU), storing the */
2001 /* result in WORK(IU), copying to A */
2002 /* CWorkspace: need 2*N [tauq, taup] + N*N [U] */
2003 /* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] */
2004 /* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork] */
2005 /* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here */
2010 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
2013 i__3 = *m - i__ + 1;
2014 chunk = f2cmin(i__3,ldwrku);
2015 clacrm_(&chunk, n, &a[i__ + a_dim1], lda, &rwork[iru],
2016 n, &work[iu], &ldwrku, &rwork[nrwork]);
2017 clacpy_("F", &chunk, n, &work[iu], &ldwrku, &a[i__ +
2025 /* Path 6s (M >= N, JOBZ='S') */
2026 /* Perform bidiagonal SVD, computing left singular vectors */
2027 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2028 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2029 /* CWorkspace: need 0 */
2030 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
2033 irvt = iru + *n * *n;
2034 nrwork = irvt + *n * *n;
2035 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
2036 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
2039 /* Copy real matrix RWORK(IRU) to complex matrix U */
2040 /* Overwrite U by left singular vectors of A */
2041 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
2042 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2043 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] */
2045 claset_("F", m, n, &c_b1, &c_b1, &u[u_offset], ldu)
2047 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
2048 i__2 = *lwork - nwork + 1;
2049 cunmbr_("Q", "L", "N", m, n, n, &a[a_offset], lda, &work[
2050 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2052 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
2053 /* Overwrite VT by right singular vectors of A */
2054 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
2055 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2056 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] */
2058 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
2059 i__2 = *lwork - nwork + 1;
2060 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
2061 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
2065 /* Path 6a (M >= N, JOBZ='A') */
2066 /* Perform bidiagonal SVD, computing left singular vectors */
2067 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2068 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2069 /* CWorkspace: need 0 */
2070 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC */
2073 irvt = iru + *n * *n;
2074 nrwork = irvt + *n * *n;
2075 sbdsdc_("U", "I", n, &s[1], &rwork[ie], &rwork[iru], n, &
2076 rwork[irvt], n, dum, idum, &rwork[nrwork], &iwork[1],
2079 /* Set the right corner of U to identity matrix */
2081 claset_("F", m, m, &c_b1, &c_b1, &u[u_offset], ldu)
2086 claset_("F", &i__2, &i__1, &c_b1, &c_b2, &u[*n + 1 + (*n
2087 + 1) * u_dim1], ldu);
2090 /* Copy real matrix RWORK(IRU) to complex matrix U */
2091 /* Overwrite U by left singular vectors of A */
2092 /* CWorkspace: need 2*N [tauq, taup] + M [work] */
2093 /* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] */
2094 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] */
2096 clacp2_("F", n, n, &rwork[iru], n, &u[u_offset], ldu);
2097 i__2 = *lwork - nwork + 1;
2098 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2099 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2101 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
2102 /* Overwrite VT by right singular vectors of A */
2103 /* CWorkspace: need 2*N [tauq, taup] + N [work] */
2104 /* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] */
2105 /* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] */
2107 clacp2_("F", n, n, &rwork[irvt], n, &vt[vt_offset], ldvt);
2108 i__2 = *lwork - nwork + 1;
2109 cunmbr_("P", "R", "C", n, n, n, &a[a_offset], lda, &work[
2110 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__2, &
2118 /* A has more columns than rows. If A has sufficiently more */
2119 /* columns than rows, first reduce using the LQ decomposition (if */
2120 /* sufficient workspace available) */
2126 /* Path 1t (N >> M, JOBZ='N') */
2127 /* No singular vectors to be computed */
2133 /* CWorkspace: need M [tau] + M [work] */
2134 /* CWorkspace: prefer M [tau] + M*NB [work] */
2135 /* RWorkspace: need 0 */
2137 i__2 = *lwork - nwork + 1;
2138 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2141 /* Zero out above L */
2145 claset_("U", &i__2, &i__1, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
2152 /* Bidiagonalize L in A */
2153 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2154 /* CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work] */
2155 /* RWorkspace: need M [e] */
2157 i__2 = *lwork - nwork + 1;
2158 cgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2159 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
2162 /* Perform bidiagonal SVD, compute singular values only */
2163 /* CWorkspace: need 0 */
2164 /* RWorkspace: need M [e] + BDSPAC */
2166 sbdsdc_("U", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2167 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2171 /* Path 2t (N >> M, JOBZ='O') */
2172 /* M right singular vectors to be overwritten on A and */
2173 /* M left singular vectors to be computed in U */
2178 /* WORK(IVT) is M by M */
2180 il = ivt + ldwkvt * *m;
2181 if (*lwork >= *m * *n + *m * *m + *m * 3) {
2183 /* WORK(IL) M by N */
2189 /* WORK(IL) is M by CHUNK */
2192 chunk = (*lwork - *m * *m - *m * 3) / *m;
2194 itau = il + ldwrkl * chunk;
2198 /* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work] */
2199 /* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] */
2200 /* RWorkspace: need 0 */
2202 i__2 = *lwork - nwork + 1;
2203 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2206 /* Copy L to WORK(IL), zeroing about above it */
2208 clacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
2211 claset_("U", &i__2, &i__1, &c_b1, &c_b1, &work[il + ldwrkl], &
2214 /* Generate Q in A */
2215 /* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work] */
2216 /* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] */
2217 /* RWorkspace: need 0 */
2219 i__2 = *lwork - nwork + 1;
2220 cunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
2227 /* Bidiagonalize L in WORK(IL) */
2228 /* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] */
2229 /* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] */
2230 /* RWorkspace: need M [e] */
2232 i__2 = *lwork - nwork + 1;
2233 cgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
2234 itauq], &work[itaup], &work[nwork], &i__2, &ierr);
2236 /* Perform bidiagonal SVD, computing left singular vectors */
2237 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2238 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2239 /* CWorkspace: need 0 */
2240 /* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2243 irvt = iru + *m * *m;
2244 nrwork = irvt + *m * *m;
2245 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2246 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2249 /* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) */
2250 /* Overwrite WORK(IU) by the left singular vectors of L */
2251 /* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] */
2252 /* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2253 /* RWorkspace: need 0 */
2255 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2256 i__2 = *lwork - nwork + 1;
2257 cunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
2258 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2260 /* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2261 /* Overwrite WORK(IVT) by the right singular vectors of L */
2262 /* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] */
2263 /* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2264 /* RWorkspace: need 0 */
2266 clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2267 i__2 = *lwork - nwork + 1;
2268 cunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
2269 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2, &
2272 /* Multiply right singular vectors of L in WORK(IL) by Q */
2273 /* in A, storing result in WORK(IL) and copying to A */
2274 /* CWorkspace: need M*M [VT] + M*M [L] */
2275 /* CWorkspace: prefer M*M [VT] + M*N [L] */
2276 /* RWorkspace: need 0 */
2280 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2283 i__3 = *n - i__ + 1;
2284 blk = f2cmin(i__3,chunk);
2285 cgemm_("N", "N", m, &blk, m, &c_b2, &work[ivt], m, &a[i__
2286 * a_dim1 + 1], lda, &c_b1, &work[il], &ldwrkl);
2287 clacpy_("F", m, &blk, &work[il], &ldwrkl, &a[i__ * a_dim1
2294 /* Path 3t (N >> M, JOBZ='S') */
2295 /* M right singular vectors to be computed in VT and */
2296 /* M left singular vectors to be computed in U */
2300 /* WORK(IL) is M by M */
2303 itau = il + ldwrkl * *m;
2307 /* CWorkspace: need M*M [L] + M [tau] + M [work] */
2308 /* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] */
2309 /* RWorkspace: need 0 */
2311 i__1 = *lwork - nwork + 1;
2312 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2315 /* Copy L to WORK(IL), zeroing out above it */
2317 clacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwrkl);
2320 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &work[il + ldwrkl], &
2323 /* Generate Q in A */
2324 /* CWorkspace: need M*M [L] + M [tau] + M [work] */
2325 /* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] */
2326 /* RWorkspace: need 0 */
2328 i__1 = *lwork - nwork + 1;
2329 cunglq_(m, n, m, &a[a_offset], lda, &work[itau], &work[nwork],
2336 /* Bidiagonalize L in WORK(IL) */
2337 /* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] */
2338 /* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] */
2339 /* RWorkspace: need M [e] */
2341 i__1 = *lwork - nwork + 1;
2342 cgebrd_(m, m, &work[il], &ldwrkl, &s[1], &rwork[ie], &work[
2343 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
2345 /* Perform bidiagonal SVD, computing left singular vectors */
2346 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2347 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2348 /* CWorkspace: need 0 */
2349 /* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2352 irvt = iru + *m * *m;
2353 nrwork = irvt + *m * *m;
2354 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2355 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2358 /* Copy real matrix RWORK(IRU) to complex matrix U */
2359 /* Overwrite U by left singular vectors of L */
2360 /* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] */
2361 /* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2362 /* RWorkspace: need 0 */
2364 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2365 i__1 = *lwork - nwork + 1;
2366 cunmbr_("Q", "L", "N", m, m, m, &work[il], &ldwrkl, &work[
2367 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2369 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
2370 /* Overwrite VT by left singular vectors of L */
2371 /* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] */
2372 /* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] */
2373 /* RWorkspace: need 0 */
2375 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2376 i__1 = *lwork - nwork + 1;
2377 cunmbr_("P", "R", "C", m, m, m, &work[il], &ldwrkl, &work[
2378 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2381 /* Copy VT to WORK(IL), multiply right singular vectors of L */
2382 /* in WORK(IL) by Q in A, storing result in VT */
2383 /* CWorkspace: need M*M [L] */
2384 /* RWorkspace: need 0 */
2386 clacpy_("F", m, m, &vt[vt_offset], ldvt, &work[il], &ldwrkl);
2387 cgemm_("N", "N", m, n, m, &c_b2, &work[il], &ldwrkl, &a[
2388 a_offset], lda, &c_b1, &vt[vt_offset], ldvt);
2392 /* Path 4t (N >> M, JOBZ='A') */
2393 /* N right singular vectors to be computed in VT and */
2394 /* M left singular vectors to be computed in U */
2398 /* WORK(IVT) is M by M */
2401 itau = ivt + ldwkvt * *m;
2404 /* Compute A=L*Q, copying result to VT */
2405 /* CWorkspace: need M*M [VT] + M [tau] + M [work] */
2406 /* CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work] */
2407 /* RWorkspace: need 0 */
2409 i__1 = *lwork - nwork + 1;
2410 cgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
2412 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2414 /* Generate Q in VT */
2415 /* CWorkspace: need M*M [VT] + M [tau] + N [work] */
2416 /* CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work] */
2417 /* RWorkspace: need 0 */
2419 i__1 = *lwork - nwork + 1;
2420 cunglq_(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
2421 nwork], &i__1, &ierr);
2423 /* Produce L in A, zeroing out above it */
2427 claset_("U", &i__1, &i__2, &c_b1, &c_b1, &a[(a_dim1 << 1) + 1]
2434 /* Bidiagonalize L in A */
2435 /* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] */
2436 /* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work] */
2437 /* RWorkspace: need M [e] */
2439 i__1 = *lwork - nwork + 1;
2440 cgebrd_(m, m, &a[a_offset], lda, &s[1], &rwork[ie], &work[
2441 itauq], &work[itaup], &work[nwork], &i__1, &ierr);
2443 /* Perform bidiagonal SVD, computing left singular vectors */
2444 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2445 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2446 /* CWorkspace: need 0 */
2447 /* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC */
2450 irvt = iru + *m * *m;
2451 nrwork = irvt + *m * *m;
2452 sbdsdc_("U", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2453 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2456 /* Copy real matrix RWORK(IRU) to complex matrix U */
2457 /* Overwrite U by left singular vectors of L */
2458 /* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] */
2459 /* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] */
2460 /* RWorkspace: need 0 */
2462 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2463 i__1 = *lwork - nwork + 1;
2464 cunmbr_("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
2465 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2467 /* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2468 /* Overwrite WORK(IVT) by right singular vectors of L */
2469 /* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] */
2470 /* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] */
2471 /* RWorkspace: need 0 */
2473 clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2474 i__1 = *lwork - nwork + 1;
2475 cunmbr_("P", "R", "C", m, m, m, &a[a_offset], lda, &work[
2476 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
2479 /* Multiply right singular vectors of L in WORK(IVT) by */
2480 /* Q in VT, storing result in A */
2481 /* CWorkspace: need M*M [VT] */
2482 /* RWorkspace: need 0 */
2484 cgemm_("N", "N", m, n, m, &c_b2, &work[ivt], &ldwkvt, &vt[
2485 vt_offset], ldvt, &c_b1, &a[a_offset], lda);
2487 /* Copy right singular vectors of A from A to VT */
2489 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2493 } else if (*n >= mnthr2) {
2495 /* MNTHR2 <= N < MNTHR1 */
2497 /* Path 5t (N >> M, but not as much as MNTHR1) */
2498 /* Reduce to bidiagonal form without QR decomposition, use */
2499 /* CUNGBR and matrix multiplication to compute singular vectors */
2507 /* Bidiagonalize A */
2508 /* CWorkspace: need 2*M [tauq, taup] + N [work] */
2509 /* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] */
2510 /* RWorkspace: need M [e] */
2512 i__1 = *lwork - nwork + 1;
2513 cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2514 &work[itaup], &work[nwork], &i__1, &ierr);
2518 /* Path 5tn (N >> M, JOBZ='N') */
2519 /* Compute singular values only */
2520 /* CWorkspace: need 0 */
2521 /* RWorkspace: need M [e] + BDSPAC */
2523 sbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2524 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2527 iru = irvt + *m * *m;
2528 nrwork = iru + *m * *m;
2531 /* Path 5to (N >> M, JOBZ='O') */
2532 /* Copy A to U, generate Q */
2533 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2534 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2535 /* RWorkspace: need 0 */
2537 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2538 i__1 = *lwork - nwork + 1;
2539 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2540 nwork], &i__1, &ierr);
2542 /* Generate P**H in A */
2543 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2544 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2545 /* RWorkspace: need 0 */
2547 i__1 = *lwork - nwork + 1;
2548 cungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[
2549 nwork], &i__1, &ierr);
2552 if (*lwork >= *m * *n + *m * 3) {
2554 /* WORK( IVT ) is M by N */
2556 nwork = ivt + ldwkvt * *n;
2560 /* WORK( IVT ) is M by CHUNK */
2562 chunk = (*lwork - *m * 3) / *m;
2563 nwork = ivt + ldwkvt * chunk;
2566 /* Perform bidiagonal SVD, computing left singular vectors */
2567 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2568 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2569 /* CWorkspace: need 0 */
2570 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2572 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2573 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2576 /* Multiply Q in U by real matrix RWORK(IRVT) */
2577 /* storing the result in WORK(IVT), copying to U */
2578 /* CWorkspace: need 2*M [tauq, taup] + M*M [VT] */
2579 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2581 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &work[ivt], &
2582 ldwkvt, &rwork[nrwork]);
2583 clacpy_("F", m, m, &work[ivt], &ldwkvt, &u[u_offset], ldu);
2585 /* Multiply RWORK(IRVT) by P**H in A, storing the */
2586 /* result in WORK(IVT), copying to A */
2587 /* CWorkspace: need 2*M [tauq, taup] + M*M [VT] */
2588 /* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] */
2589 /* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork] */
2590 /* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2595 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ +=
2598 i__3 = *n - i__ + 1;
2599 blk = f2cmin(i__3,chunk);
2600 clarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1],
2601 lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
2602 clacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ *
2608 /* Path 5ts (N >> M, JOBZ='S') */
2609 /* Copy A to U, generate Q */
2610 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2611 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2612 /* RWorkspace: need 0 */
2614 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2615 i__2 = *lwork - nwork + 1;
2616 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2617 nwork], &i__2, &ierr);
2619 /* Copy A to VT, generate P**H */
2620 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2621 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2622 /* RWorkspace: need 0 */
2624 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2625 i__2 = *lwork - nwork + 1;
2626 cungbr_("P", m, n, m, &vt[vt_offset], ldvt, &work[itaup], &
2627 work[nwork], &i__2, &ierr);
2629 /* Perform bidiagonal SVD, computing left singular vectors */
2630 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2631 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2632 /* CWorkspace: need 0 */
2633 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2636 iru = irvt + *m * *m;
2637 nrwork = iru + *m * *m;
2638 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2639 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2642 /* Multiply Q in U by real matrix RWORK(IRU), storing the */
2643 /* result in A, copying to U */
2644 /* CWorkspace: need 0 */
2645 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2647 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset],
2648 lda, &rwork[nrwork]);
2649 clacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2651 /* Multiply real matrix RWORK(IRVT) by P**H in VT, */
2652 /* storing the result in A, copying to VT */
2653 /* CWorkspace: need 0 */
2654 /* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2657 clarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
2658 a_offset], lda, &rwork[nrwork]);
2659 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2662 /* Path 5ta (N >> M, JOBZ='A') */
2663 /* Copy A to U, generate Q */
2664 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2665 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2666 /* RWorkspace: need 0 */
2668 clacpy_("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2669 i__2 = *lwork - nwork + 1;
2670 cungbr_("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[
2671 nwork], &i__2, &ierr);
2673 /* Copy A to VT, generate P**H */
2674 /* CWorkspace: need 2*M [tauq, taup] + N [work] */
2675 /* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] */
2676 /* RWorkspace: need 0 */
2678 clacpy_("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2679 i__2 = *lwork - nwork + 1;
2680 cungbr_("P", n, n, m, &vt[vt_offset], ldvt, &work[itaup], &
2681 work[nwork], &i__2, &ierr);
2683 /* Perform bidiagonal SVD, computing left singular vectors */
2684 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2685 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2686 /* CWorkspace: need 0 */
2687 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2690 iru = irvt + *m * *m;
2691 nrwork = iru + *m * *m;
2692 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2693 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2696 /* Multiply Q in U by real matrix RWORK(IRU), storing the */
2697 /* result in A, copying to U */
2698 /* CWorkspace: need 0 */
2699 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] */
2701 clacrm_(m, m, &u[u_offset], ldu, &rwork[iru], m, &a[a_offset],
2702 lda, &rwork[nrwork]);
2703 clacpy_("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
2705 /* Multiply real matrix RWORK(IRVT) by P**H in VT, */
2706 /* storing the result in A, copying to VT */
2707 /* CWorkspace: need 0 */
2708 /* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2711 clarcm_(m, n, &rwork[irvt], m, &vt[vt_offset], ldvt, &a[
2712 a_offset], lda, &rwork[nrwork]);
2713 clacpy_("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
2720 /* Path 6t (N > M, but not much larger) */
2721 /* Reduce to bidiagonal form without LQ decomposition */
2722 /* Use CUNMBR to compute singular vectors */
2730 /* Bidiagonalize A */
2731 /* CWorkspace: need 2*M [tauq, taup] + N [work] */
2732 /* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] */
2733 /* RWorkspace: need M [e] */
2735 i__2 = *lwork - nwork + 1;
2736 cgebrd_(m, n, &a[a_offset], lda, &s[1], &rwork[ie], &work[itauq],
2737 &work[itaup], &work[nwork], &i__2, &ierr);
2740 /* Path 6tn (N > M, JOBZ='N') */
2741 /* Compute singular values only */
2742 /* CWorkspace: need 0 */
2743 /* RWorkspace: need M [e] + BDSPAC */
2745 sbdsdc_("L", "N", m, &s[1], &rwork[ie], dum, &c__1, dum, &
2746 c__1, dum, idum, &rwork[nrwork], &iwork[1], info);
2748 /* Path 6to (N > M, JOBZ='O') */
2751 if (*lwork >= *m * *n + *m * 3) {
2753 /* WORK( IVT ) is M by N */
2755 claset_("F", m, n, &c_b1, &c_b1, &work[ivt], &ldwkvt);
2756 nwork = ivt + ldwkvt * *n;
2759 /* WORK( IVT ) is M by CHUNK */
2761 chunk = (*lwork - *m * 3) / *m;
2762 nwork = ivt + ldwkvt * chunk;
2765 /* Perform bidiagonal SVD, computing left singular vectors */
2766 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2767 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2768 /* CWorkspace: need 0 */
2769 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2772 iru = irvt + *m * *m;
2773 nrwork = iru + *m * *m;
2774 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2775 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2778 /* Copy real matrix RWORK(IRU) to complex matrix U */
2779 /* Overwrite U by left singular vectors of A */
2780 /* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work] */
2781 /* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] */
2782 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] */
2784 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2785 i__2 = *lwork - nwork + 1;
2786 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2787 itauq], &u[u_offset], ldu, &work[nwork], &i__2, &ierr);
2789 if (*lwork >= *m * *n + *m * 3) {
2792 /* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) */
2793 /* Overwrite WORK(IVT) by right singular vectors of A, */
2795 /* CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work] */
2796 /* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work] */
2797 /* RWorkspace: need M [e] + M*M [RVT] */
2799 clacp2_("F", m, m, &rwork[irvt], m, &work[ivt], &ldwkvt);
2800 i__2 = *lwork - nwork + 1;
2801 cunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
2802 itaup], &work[ivt], &ldwkvt, &work[nwork], &i__2,
2804 clacpy_("F", m, n, &work[ivt], &ldwkvt, &a[a_offset], lda);
2808 /* Generate P**H in A */
2809 /* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work] */
2810 /* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] */
2811 /* RWorkspace: need 0 */
2813 i__2 = *lwork - nwork + 1;
2814 cungbr_("P", m, n, m, &a[a_offset], lda, &work[itaup], &
2815 work[nwork], &i__2, &ierr);
2817 /* Multiply Q in A by real matrix RWORK(IRU), storing the */
2818 /* result in WORK(IU), copying to A */
2819 /* CWorkspace: need 2*M [tauq, taup] + M*M [VT] */
2820 /* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] */
2821 /* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork] */
2822 /* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here */
2827 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ +=
2830 i__3 = *n - i__ + 1;
2831 blk = f2cmin(i__3,chunk);
2832 clarcm_(m, &blk, &rwork[irvt], m, &a[i__ * a_dim1 + 1]
2833 , lda, &work[ivt], &ldwkvt, &rwork[nrwork]);
2834 clacpy_("F", m, &blk, &work[ivt], &ldwkvt, &a[i__ *
2841 /* Path 6ts (N > M, JOBZ='S') */
2842 /* Perform bidiagonal SVD, computing left singular vectors */
2843 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2844 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2845 /* CWorkspace: need 0 */
2846 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2849 iru = irvt + *m * *m;
2850 nrwork = iru + *m * *m;
2851 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2852 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2855 /* Copy real matrix RWORK(IRU) to complex matrix U */
2856 /* Overwrite U by left singular vectors of A */
2857 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2858 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2859 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] */
2861 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2862 i__1 = *lwork - nwork + 1;
2863 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2864 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2866 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
2867 /* Overwrite VT by right singular vectors of A */
2868 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2869 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2870 /* RWorkspace: need M [e] + M*M [RVT] */
2872 claset_("F", m, n, &c_b1, &c_b1, &vt[vt_offset], ldvt);
2873 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2874 i__1 = *lwork - nwork + 1;
2875 cunmbr_("P", "R", "C", m, n, m, &a[a_offset], lda, &work[
2876 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2880 /* Path 6ta (N > M, JOBZ='A') */
2881 /* Perform bidiagonal SVD, computing left singular vectors */
2882 /* of bidiagonal matrix in RWORK(IRU) and computing right */
2883 /* singular vectors of bidiagonal matrix in RWORK(IRVT) */
2884 /* CWorkspace: need 0 */
2885 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC */
2888 iru = irvt + *m * *m;
2889 nrwork = iru + *m * *m;
2891 sbdsdc_("L", "I", m, &s[1], &rwork[ie], &rwork[iru], m, &
2892 rwork[irvt], m, dum, idum, &rwork[nrwork], &iwork[1],
2895 /* Copy real matrix RWORK(IRU) to complex matrix U */
2896 /* Overwrite U by left singular vectors of A */
2897 /* CWorkspace: need 2*M [tauq, taup] + M [work] */
2898 /* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] */
2899 /* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] */
2901 clacp2_("F", m, m, &rwork[iru], m, &u[u_offset], ldu);
2902 i__1 = *lwork - nwork + 1;
2903 cunmbr_("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
2904 itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
2906 /* Set all of VT to identity matrix */
2908 claset_("F", n, n, &c_b1, &c_b2, &vt[vt_offset], ldvt);
2910 /* Copy real matrix RWORK(IRVT) to complex matrix VT */
2911 /* Overwrite VT by right singular vectors of A */
2912 /* CWorkspace: need 2*M [tauq, taup] + N [work] */
2913 /* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] */
2914 /* RWorkspace: need M [e] + M*M [RVT] */
2916 clacp2_("F", m, m, &rwork[irvt], m, &vt[vt_offset], ldvt);
2917 i__1 = *lwork - nwork + 1;
2918 cunmbr_("P", "R", "C", n, n, m, &a[a_offset], lda, &work[
2919 itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
2927 /* Undo scaling if necessary */
2930 if (anrm > bignum) {
2931 slascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
2934 if (*info != 0 && anrm > bignum) {
2936 slascl_("G", &c__0, &c__0, &bignum, &anrm, &i__1, &c__1, &rwork[
2937 ie], &minmn, &ierr);
2939 if (anrm < smlnum) {
2940 slascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
2943 if (*info != 0 && anrm < smlnum) {
2945 slascl_("G", &c__0, &c__0, &smlnum, &anrm, &i__1, &c__1, &rwork[
2946 ie], &minmn, &ierr);
2950 /* Return optimal workspace in WORK(1) */
2952 work[1].r = (real) maxwrk, work[1].i = 0.f;