praat/dwsys/NUMclapack.cpp

/* NUMclapack.c */

/* -- LAPACK driver routines (version 3.0) -- Univ. of Tennessee, Univ. of
   California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab,
   and Rice University October 31, 1999 -- translated by f2c (version
   19990503).

	Adapted by David Weenink 20021201

 djmw 20030205 Latest modification
*/
/* #include "blaswrap.h" */
#include "NUMf2c.h"
#include "NUMclapack.h"
#include "NUMcblas.h"
#include "NUM2.h"
#include "melder.h"
#define FALSE  0
#define TRUE  1

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;
static integer c_n1 = -1;
static integer c__2 = 2;
static integer c__3 = 3;
static integer c__4 = 4;
static integer c__6 = 6;
static integer c__10 = 10;
static integer c__11 = 11;
static double c_b15 = -.125;
static double c_b49 = 1.;
static double c_b72 = -1.;
static double c_b74 = 0.;
static double c_b108 = 1.;
static double c_b416 = 0.;
static double c_b438 = 1.;

#define MAX(m,n) ((m) > (n) ? (m) : (n))
#define MIN(m,n) ((m) < (n) ? (m) : (n))

#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]

/* --------------------------------------------------- */

#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1]

int NUMlapack_dbdsqr (const char *uplo, integer *n, integer *ncvt, integer *nru, integer *ncc, double *d__, double *e, double *vt,
                      integer *ldvt, double *u, integer *ldu, double *c__, integer *ldc, double *work, integer *info) {
	/* System generated locals */
	integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
	double d__1, d__2, d__3, d__4;

	/* Local variables */
	static double abse;
	static integer idir;
	static double abss;
	static integer oldm;
	static double cosl;
	static integer isub, iter;
	static double unfl, sinl, cosr, smin, smax, sinr;
	static double f, g, h__;
	static integer i__, j, m;
	static double r__;
	static double oldcs;
	static integer oldll;
	static double shift, sigmn, oldsn;
	static integer maxit;
	static double sminl, sigmx;
	static integer lower;
	static double cs;
	static integer ll;
	static double sn, mu;
	static double sminoa, thresh;
	static integer rotate;
	static double sminlo;
	static integer nm1;
	static double tolmul;
	static integer nm12, nm13, lll;
	static double eps, sll, tol;

	/* Parameter adjustments */
	--d__;
	--e;
	vt_dim1 = *ldvt;
	vt_offset = 1 + vt_dim1 * 1;
	vt -= vt_offset;
	u_dim1 = *ldu;
	u_offset = 1 + u_dim1 * 1;
	u -= u_offset;
	c_dim1 = *ldc;
	c_offset = 1 + c_dim1 * 1;
	c__ -= c_offset;
	--work;

	/* Function Body */
	*info = 0;
	lower = lsame_ (uplo, "L");
	if (!lsame_ (uplo, "U") && !lower) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*ncvt < 0) {
		*info = -3;
	} else if (*nru < 0) {
		*info = -4;
	} else if (*ncc < 0) {
		*info = -5;
	} else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < MAX (1, *n)) {
		*info = -9;
	} else if (*ldu < MAX (1, *nru)) {
		*info = -11;
	} else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < MAX (1, *n)) {
		*info = -13;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DBDSQR", &i__1);
		return 0;
	}
	if (*n == 0) {
		return 0;
	}
	if (*n == 1) {
		goto L160;
	}

	/* ROTATE is true if any singular vectors desired, false otherwise */

	rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;

	/* If no singular vectors desired, use qd algorithm */

	if (!rotate) {
		NUMlapack_dlasq1 (n, &d__[1], &e[1], &work[1], info);
		return 0;
	}

	nm1 = *n - 1;
	nm12 = nm1 + nm1;
	nm13 = nm12 + nm1;
	idir = 0;

	/* Get machine constants */

	eps = NUMblas_dlamch ("Epsilon");
	unfl = NUMblas_dlamch ("Safe minimum");

	/* If matrix lower bidiagonal, rotate to be upper bidiagonal by applying
	   Givens rotations on the left */

	if (lower) {
		i__1 = *n - 1;
		for (i__ = 1; i__ <= i__1; ++i__) {
			NUMlapack_dlartg (&d__[i__], &e[i__], &cs, &sn, &r__);
			d__[i__] = r__;
			e[i__] = sn * d__[i__ + 1];
			d__[i__ + 1] = cs * d__[i__ + 1];
			work[i__] = cs;
			work[nm1 + i__] = sn;
			/* L10: */
		}

		/* Update singular vectors if desired */

		if (*nru > 0) {
			NUMlapack_dlasr ("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset], ldu);
		}
		if (*ncc > 0) {
			NUMlapack_dlasr ("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset], ldc);
		}
	}

	/* Compute singular values to relative accuracy TOL (By setting TOL to be
	   negative, algorithm will compute singular values to absolute accuracy
	   ABS(TOL)*norm(input matrix))

	   Computing MAX Computing MIN */
	d__3 = 100., d__4 = pow (eps, c_b15);
	d__1 = 10., d__2 = MIN (d__3, d__4);
	tolmul = MAX (d__1, d__2);
	tol = tolmul * eps;

	/* Compute approximate maximum, minimum singular values */

	smax = 0.;
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
		/* Computing MAX */
		d__2 = smax, d__3 = (d__1 = d__[i__], fabs (d__1));
		smax = MAX (d__2, d__3);
		/* L20: */
	}
	i__1 = *n - 1;
	for (i__ = 1; i__ <= i__1; ++i__) {
		/* Computing MAX */
		d__2 = smax, d__3 = (d__1 = e[i__], fabs (d__1));
		smax = MAX (d__2, d__3);
		/* L30: */
	}
	sminl = 0.;
	if (tol >= 0.) {

		/* Relative accuracy desired */

		sminoa = fabs (d__[1]);
		if (sminoa == 0.) {
			goto L50;
		}
		mu = sminoa;
		i__1 = *n;
		for (i__ = 2; i__ <= i__1; ++i__) {
			mu = (d__2 = d__[i__], fabs (d__2)) * (mu / (mu + (d__1 = e[i__ - 1], fabs (d__1))));
			sminoa = MIN (sminoa, mu);
			if (sminoa == 0.) {
				goto L50;
			}
			/* L40: */
		}
L50:
		sminoa /= sqrt ( (double) (*n));
		/* Computing MAX */
		d__1 = tol * sminoa, d__2 = *n * 6 * *n * unfl;
		thresh = MAX (d__1, d__2);
	} else {

		/* Absolute accuracy desired

		   Computing MAX */
		d__1 = fabs (tol) * smax, d__2 = *n * 6 * *n * unfl;
		thresh = MAX (d__1, d__2);
	}

	/* Prepare for main iteration loop for the singular values (MAXIT is the
	   maximum number of passes through the inner loop permitted before
	   nonconvergence signalled.) */

	maxit = *n * 6 * *n;
	iter = 0;
	oldll = -1;
	oldm = -1;

	/* M points to last element of unconverged part of matrix */

	m = *n;

	/* Begin main iteration loop */

L60:

	/* Check for convergence or exceeding iteration count */

	if (m <= 1) {
		goto L160;
	}
	if (iter > maxit) {
		goto L200;
	}

	/* Find diagonal block of matrix to work on */

	if (tol < 0. && (d__1 = d__[m], fabs (d__1)) <= thresh) {
		d__[m] = 0.;
	}
	smax = (d__1 = d__[m], fabs (d__1));
	smin = smax;
	i__1 = m - 1;
	for (lll = 1; lll <= i__1; ++lll) {
		ll = m - lll;
		abss = (d__1 = d__[ll], fabs (d__1));
		abse = (d__1 = e[ll], fabs (d__1));
		if (tol < 0. && abss <= thresh) {
			d__[ll] = 0.;
		}
		if (abse <= thresh) {
			goto L80;
		}
		smin = MIN (smin, abss);
		/* Computing MAX */
		d__1 = MAX (smax, abss);
		smax = MAX (d__1, abse);
		/* L70: */
	}
	ll = 0;
	goto L90;
L80:
	e[ll] = 0.;

	/* Matrix splits since E(LL) = 0 */

	if (ll == m - 1) {

		/* Convergence of bottom singular value, return to top of loop */

		--m;
		goto L60;
	}
L90:
	++ll;

	/* E(LL) through E(M-1) are nonzero, E(LL-1) is zero */

	if (ll == m - 1) {

		/* 2 by 2 block, handle separately */

		NUMlapack_dlasv2 (&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr, &sinl, &cosl);
		d__[m - 1] = sigmx;
		e[m - 1] = 0.;
		d__[m] = sigmn;

		/* Compute singular vectors, if desired */

		if (*ncvt > 0) {
			NUMblas_drot (ncvt, &vt_ref (m - 1, 1), ldvt, &vt_ref (m, 1), ldvt, &cosr, &sinr);
		}
		if (*nru > 0) {
			NUMblas_drot (nru, &u_ref (1, m - 1), &c__1, &u_ref (1, m), &c__1, &cosl, &sinl);
		}
		if (*ncc > 0) {
			NUMblas_drot (ncc, &c___ref (m - 1, 1), ldc, &c___ref (m, 1), ldc, &cosl, &sinl);
		}
		m += -2;
		goto L60;
	}

	/* If working on new submatrix, choose shift direction (from larger end
	   diagonal element towards smaller) */

	if (ll > oldm || m < oldll) {
		if ( (d__1 = d__[ll], fabs (d__1)) >= (d__2 = d__[m], fabs (d__2))) {

			/* Chase bulge from top (big end) to bottom (small end) */

			idir = 1;
		} else {

			/* Chase bulge from bottom (big end) to top (small end) */

			idir = 2;
		}
	}

	/* Apply convergence tests */

	if (idir == 1) {

		/* Run convergence test in forward direction First apply standard
		   test to bottom of matrix */

		if ( (d__2 = e[m - 1], fabs (d__2)) <= fabs (tol) * (d__1 = d__[m], fabs (d__1)) || tol < 0. &&
		        (d__3 = e[m - 1], fabs (d__3)) <= thresh) {
			e[m - 1] = 0.;
			goto L60;
		}

		if (tol >= 0.) {

			/* If relative accuracy desired, apply convergence criterion
			   forward */

			mu = (d__1 = d__[ll], fabs (d__1));
			sminl = mu;
			i__1 = m - 1;
			for (lll = ll; lll <= i__1; ++lll) {
				if ( (d__1 = e[lll], fabs (d__1)) <= tol * mu) {
					e[lll] = 0.;
					goto L60;
				}
				sminlo = sminl;
				mu = (d__2 = d__[lll + 1], fabs (d__2)) * (mu / (mu + (d__1 = e[lll], fabs (d__1))));
				sminl = MIN (sminl, mu);
				/* L100: */
			}
		}

	} else {

		/* Run convergence test in backward direction First apply standard
		   test to top of matrix */

		if ( (d__2 = e[ll], fabs (d__2)) <= fabs (tol) * (d__1 = d__[ll], fabs (d__1)) || tol < 0. &&
		        (d__3 = e[ll], fabs (d__3)) <= thresh) {
			e[ll] = 0.;
			goto L60;
		}

		if (tol >= 0.) {

			/* If relative accuracy desired, apply convergence criterion
			   backward */

			mu = (d__1 = d__[m], fabs (d__1));
			sminl = mu;
			i__1 = ll;
			for (lll = m - 1; lll >= i__1; --lll) {
				if ( (d__1 = e[lll], fabs (d__1)) <= tol * mu) {
					e[lll] = 0.;
					goto L60;
				}
				sminlo = sminl;
				mu = (d__2 = d__[lll], fabs (d__2)) * (mu / (mu + (d__1 = e[lll], fabs (d__1))));
				sminl = MIN (sminl, mu);
				/* L110: */
			}
		}
	}
	oldll = ll;
	oldm = m;

	/* Compute shift.  First, test if shifting would ruin relative accuracy,
	   and if so set the shift to zero.

	   Computing MAX */
	d__1 = eps, d__2 = tol * .01;
	if (tol >= 0. && *n * tol * (sminl / smax) <= MAX (d__1, d__2)) {

		/* Use a zero shift to avoid loss of relative accuracy */

		shift = 0.;
	} else {

		/* Compute the shift from 2-by-2 block at end of matrix */

		if (idir == 1) {
			sll = (d__1 = d__[ll], fabs (d__1));
			NUMlapack_dlas2 (&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
		} else {
			sll = (d__1 = d__[m], fabs (d__1));
			NUMlapack_dlas2 (&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
		}

		/* Test if shift negligible, and if so set to zero */

		if (sll > 0.) {
			/* Computing 2nd power */
			d__1 = shift / sll;
			if (d__1 * d__1 < eps) {
				shift = 0.;
			}
		}
	}

	/* Increment iteration count */

	iter = iter + m - ll;

	/* If SHIFT = 0, do simplified QR iteration */

	if (shift == 0.) {
		if (idir == 1) {

			/* Chase bulge from top to bottom Save cosines and sines for
			   later singular vector updates */

			cs = 1.;
			oldcs = 1.;
			i__1 = m - 1;
			for (i__ = ll; i__ <= i__1; ++i__) {
				d__1 = d__[i__] * cs;
				NUMlapack_dlartg (&d__1, &e[i__], &cs, &sn, &r__);
				if (i__ > ll) {
					e[i__ - 1] = oldsn * r__;
				}
				d__1 = oldcs * r__;
				d__2 = d__[i__ + 1] * sn;
				NUMlapack_dlartg (&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
				work[i__ - ll + 1] = cs;
				work[i__ - ll + 1 + nm1] = sn;
				work[i__ - ll + 1 + nm12] = oldcs;
				work[i__ - ll + 1 + nm13] = oldsn;
				/* L120: */
			}
			h__ = d__[m] * cs;
			d__[m] = h__ * oldcs;
			e[m - 1] = h__ * oldsn;

			/* Update singular vectors */

			if (*ncvt > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt_ref (ll, 1), ldvt);
			}
			if (*nru > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u_ref (1, ll),
				                 ldu);
			}
			if (*ncc > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c___ref (ll, 1),
				                 ldc);
			}

			/* Test convergence */

			if ( (d__1 = e[m - 1], fabs (d__1)) <= thresh) {
				e[m - 1] = 0.;
			}

		} else {

			/* Chase bulge from bottom to top Save cosines and sines for
			   later singular vector updates */

			cs = 1.;
			oldcs = 1.;
			i__1 = ll + 1;
			for (i__ = m; i__ >= i__1; --i__) {
				d__1 = d__[i__] * cs;
				NUMlapack_dlartg (&d__1, &e[i__ - 1], &cs, &sn, &r__);
				if (i__ < m) {
					e[i__] = oldsn * r__;
				}
				d__1 = oldcs * r__;
				d__2 = d__[i__ - 1] * sn;
				NUMlapack_dlartg (&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
				work[i__ - ll] = cs;
				work[i__ - ll + nm1] = -sn;
				work[i__ - ll + nm12] = oldcs;
				work[i__ - ll + nm13] = -oldsn;
				/* L130: */
			}
			h__ = d__[ll] * cs;
			d__[ll] = h__ * oldcs;
			e[ll] = h__ * oldsn;

			/* Update singular vectors */

			if (*ncvt > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[nm13 + 1], &vt_ref (ll, 1),
				                 ldvt);
			}
			if (*nru > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u_ref (1, ll), ldu);
			}
			if (*ncc > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c___ref (ll, 1), ldc);
			}

			/* Test convergence */

			if ( (d__1 = e[ll], fabs (d__1)) <= thresh) {
				e[ll] = 0.;
			}
		}
	} else {

		/* Use nonzero shift */

		if (idir == 1) {

			/* Chase bulge from top to bottom Save cosines and sines for
			   later singular vector updates */

			f = ( (d__1 = d__[ll], fabs (d__1)) - shift) * (d_sign (&c_b49, &d__[ll]) + shift / d__[ll]);
			g = e[ll];
			i__1 = m - 1;
			for (i__ = ll; i__ <= i__1; ++i__) {
				NUMlapack_dlartg (&f, &g, &cosr, &sinr, &r__);
				if (i__ > ll) {
					e[i__ - 1] = r__;
				}
				f = cosr * d__[i__] + sinr * e[i__];
				e[i__] = cosr * e[i__] - sinr * d__[i__];
				g = sinr * d__[i__ + 1];
				d__[i__ + 1] = cosr * d__[i__ + 1];
				NUMlapack_dlartg (&f, &g, &cosl, &sinl, &r__);
				d__[i__] = r__;
				f = cosl * e[i__] + sinl * d__[i__ + 1];
				d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
				if (i__ < m - 1) {
					g = sinl * e[i__ + 1];
					e[i__ + 1] = cosl * e[i__ + 1];
				}
				work[i__ - ll + 1] = cosr;
				work[i__ - ll + 1 + nm1] = sinr;
				work[i__ - ll + 1 + nm12] = cosl;
				work[i__ - ll + 1 + nm13] = sinl;
				/* L140: */
			}
			e[m - 1] = f;

			/* Update singular vectors */

			if (*ncvt > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt_ref (ll, 1), ldvt);
			}
			if (*nru > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u_ref (1, ll),
				                 ldu);
			}
			if (*ncc > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c___ref (ll, 1),
				                 ldc);
			}

			/* Test convergence */

			if ( (d__1 = e[m - 1], fabs (d__1)) <= thresh) {
				e[m - 1] = 0.;
			}

		} else {

			/* Chase bulge from bottom to top Save cosines and sines for
			   later singular vector updates */

			f = ( (d__1 = d__[m], fabs (d__1)) - shift) * (d_sign (&c_b49, &d__[m]) + shift / d__[m]);
			g = e[m - 1];
			i__1 = ll + 1;
			for (i__ = m; i__ >= i__1; --i__) {
				NUMlapack_dlartg (&f, &g, &cosr, &sinr, &r__);
				if (i__ < m) {
					e[i__] = r__;
				}
				f = cosr * d__[i__] + sinr * e[i__ - 1];
				e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
				g = sinr * d__[i__ - 1];
				d__[i__ - 1] = cosr * d__[i__ - 1];
				NUMlapack_dlartg (&f, &g, &cosl, &sinl, &r__);
				d__[i__] = r__;
				f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
				d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
				if (i__ > ll + 1) {
					g = sinl * e[i__ - 2];
					e[i__ - 2] = cosl * e[i__ - 2];
				}
				work[i__ - ll] = cosr;
				work[i__ - ll + nm1] = -sinr;
				work[i__ - ll + nm12] = cosl;
				work[i__ - ll + nm13] = -sinl;
				/* L150: */
			}
			e[ll] = f;

			/* Test convergence */

			if ( (d__1 = e[ll], fabs (d__1)) <= thresh) {
				e[ll] = 0.;
			}

			/* Update singular vectors if desired */

			if (*ncvt > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[nm13 + 1], &vt_ref (ll, 1),
				                 ldvt);
			}
			if (*nru > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u_ref (1, ll), ldu);
			}
			if (*ncc > 0) {
				i__1 = m - ll + 1;
				NUMlapack_dlasr ("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c___ref (ll, 1), ldc);
			}
		}
	}

	/* QR iteration finished, go back and check convergence */

	goto L60;

	/* All singular values converged, so make them positive */

L160:
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
		if (d__[i__] < 0.) {
			d__[i__] = -d__[i__];

			/* Change sign of singular vectors, if desired */

			if (*ncvt > 0) {
				NUMblas_dscal (ncvt, &c_b72, &vt_ref (i__, 1), ldvt);
			}
		}
		/* L170: */
	}

	/* Sort the singular values into decreasing order (insertion sort on
	   singular values, but only one transposition per singular vector) */

	i__1 = *n - 1;
	for (i__ = 1; i__ <= i__1; ++i__) {

		/* Scan for smallest D(I) */

		isub = 1;
		smin = d__[1];
		i__2 = *n + 1 - i__;
		for (j = 2; j <= i__2; ++j) {
			if (d__[j] <= smin) {
				isub = j;
				smin = d__[j];
			}
			/* L180: */
		}
		if (isub != *n + 1 - i__) {

			/* Swap singular values and vectors */

			d__[isub] = d__[*n + 1 - i__];
			d__[*n + 1 - i__] = smin;
			if (*ncvt > 0) {
				NUMblas_dswap (ncvt, &vt_ref (isub, 1), ldvt, &vt_ref (*n + 1 - i__, 1), ldvt);
			}
			if (*nru > 0) {
				NUMblas_dswap (nru, &u_ref (1, isub), &c__1, &u_ref (1, *n + 1 - i__), &c__1);
			}
			if (*ncc > 0) {
				NUMblas_dswap (ncc, &c___ref (isub, 1), ldc, &c___ref (*n + 1 - i__, 1), ldc);
			}
		}
		/* L190: */
	}
	goto L220;

	/* Maximum number of iterations exceeded, failure to converge */

L200:
	*info = 0;
	i__1 = *n - 1;
	for (i__ = 1; i__ <= i__1; ++i__) {
		if (e[i__] != 0.) {
			++ (*info);
		}
		/* L210: */
	}
L220:
	return 0;
}								/* NUMlapack_dbdsqr */

#undef vt_ref
#undef u_ref


int NUMlapack_dgebd2 (integer *m, integer *n, double *a, integer *lda, double *d__, double *e, double *tauq,
                      double *taup, double *work, integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3, i__4;

	/* Local variables */
	static integer i__;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--d__;
	--e;
	--tauq;
	--taup;
	--work;

	/* Function Body */
	*info = 0;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	}
	if (*info < 0) {
		i__1 = - (*info);
		xerbla_ ("DGEBD2", &i__1);
		return 0;
	}

	if (*m >= *n) {

		/* Reduce to upper bidiagonal form */

		i__1 = *n;
		for (i__ = 1; i__ <= i__1; ++i__) {

			/* Generate elementary reflector H(i) to annihilate A(i+1:m,i)

			   Computing MIN */
			i__2 = i__ + 1;
			i__3 = *m - i__ + 1;
			NUMlapack_dlarfg (&i__3, &a_ref (i__, i__), &a_ref (MIN (i__2, *m), i__), &c__1, &tauq[i__]);
			d__[i__] = a_ref (i__, i__);
			a_ref (i__, i__) = 1.;

			/* Apply H(i) to A(i:m,i+1:n) from the left */

			i__2 = *m - i__ + 1;
			i__3 = *n - i__;
			NUMlapack_dlarf ("Left", &i__2, &i__3, &a_ref (i__, i__), &c__1, &tauq[i__], &a_ref (i__, i__ + 1),
			                 lda, &work[1]);
			a_ref (i__, i__) = d__[i__];

			if (i__ < *n) {

				/* Generate elementary reflector G(i) to annihilate
				   A(i,i+2:n)

				   Computing MIN */
				i__2 = i__ + 2;
				i__3 = *n - i__;
				NUMlapack_dlarfg (&i__3, &a_ref (i__, i__ + 1), &a_ref (i__, MIN (i__2, *n)), lda, &taup[i__]);
				e[i__] = a_ref (i__, i__ + 1);
				a_ref (i__, i__ + 1) = 1.;

				/* Apply G(i) to A(i+1:m,i+1:n) from the right */

				i__2 = *m - i__;
				i__3 = *n - i__;
				NUMlapack_dlarf ("Right", &i__2, &i__3, &a_ref (i__, i__ + 1), lda, &taup[i__], &a_ref (i__ + 1,
				                 i__ + 1), lda, &work[1]);
				a_ref (i__, i__ + 1) = e[i__];
			} else {
				taup[i__] = 0.;
			}
			/* L10: */
		}
	} else {

		/* Reduce to lower bidiagonal form */

		i__1 = *m;
		for (i__ = 1; i__ <= i__1; ++i__) {

			/* Generate elementary reflector G(i) to annihilate A(i,i+1:n)

			   Computing MIN */
			i__2 = i__ + 1;
			i__3 = *n - i__ + 1;
			NUMlapack_dlarfg (&i__3, &a_ref (i__, i__), &a_ref (i__, MIN (i__2, *n)), lda, &taup[i__]);
			d__[i__] = a_ref (i__, i__);
			a_ref (i__, i__) = 1.;

			/* Apply G(i) to A(i+1:m,i:n) from the right

			   Computing MIN */
			i__2 = i__ + 1;
			i__3 = *m - i__;
			i__4 = *n - i__ + 1;
			NUMlapack_dlarf ("Right", &i__3, &i__4, &a_ref (i__, i__), lda, &taup[i__], &a_ref (MIN (i__2, *m),
			                 i__), lda, &work[1]);
			a_ref (i__, i__) = d__[i__];

			if (i__ < *m) {

				/* Generate elementary reflector H(i) to annihilate
				   A(i+2:m,i)

				   Computing MIN */
				i__2 = i__ + 2;
				i__3 = *m - i__;
				NUMlapack_dlarfg (&i__3, &a_ref (i__ + 1, i__), &a_ref (MIN (i__2, *m), i__), &c__1,
				                  &tauq[i__]);
				e[i__] = a_ref (i__ + 1, i__);
				a_ref (i__ + 1, i__) = 1.;

				/* Apply H(i) to A(i+1:m,i+1:n) from the left */

				i__2 = *m - i__;
				i__3 = *n - i__;
				NUMlapack_dlarf ("Left", &i__2, &i__3, &a_ref (i__ + 1, i__), &c__1, &tauq[i__],
				                 &a_ref (i__ + 1, i__ + 1), lda, &work[1]);
				a_ref (i__ + 1, i__) = e[i__];
			} else {
				tauq[i__] = 0.;
			}
			/* L20: */
		}
	}
	return 0;
}								/* NUMlapack_dgebd2 */


int NUMlapack_dgebak (const char *job, const char *side, integer *n, integer *ilo, integer *ihi, double *scale, integer *m,
                      double *v, integer *ldv, integer *info) {
	/* System generated locals */
	integer v_dim1, v_offset, i__1;

	/* Local variables */
	static integer i__, k;
	static double s;
	static int leftv;
	static integer ii;
	static int rightv;

#define v_ref(a_1,a_2) v[(a_2)*v_dim1 + a_1]

	--scale;
	v_dim1 = *ldv;
	v_offset = 1 + v_dim1 * 1;
	v -= v_offset;

	/* Function Body */
	rightv = lsame_ (side, "R");
	leftv = lsame_ (side, "L");

	*info = 0;
	if (!lsame_ (job, "N") && !lsame_ (job, "P") && !lsame_ (job, "S") && !lsame_ (job, "B")) {
		*info = -1;
	} else if (!rightv && !leftv) {
		*info = -2;
	} else if (*n < 0) {
		*info = -3;
	} else if (*ilo < 1 || *ilo > MAX (1, *n)) {
		*info = -4;
	} else if (*ihi < MIN (*ilo, *n) || *ihi > *n) {
		*info = -5;
	} else if (*m < 0) {
		*info = -7;
	} else if (*ldv < MAX (1, *n)) {
		*info = -9;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("NUMlapack_dgebak ", &i__1);
		return 0;
	}

	/* Quick return if possible */

	if (*n == 0) {
		return 0;
	}
	if (*m == 0) {
		return 0;
	}
	if (lsame_ (job, "N")) {
		return 0;
	}

	if (*ilo == *ihi) {
		goto L30;
	}

	/* Backward balance */

	if (lsame_ (job, "S") || lsame_ (job, "B")) {

		if (rightv) {
			i__1 = *ihi;
			for (i__ = *ilo; i__ <= i__1; ++i__) {
				s = scale[i__];
				NUMblas_dscal (m, &s, &v_ref (i__, 1), ldv);
				/* L10: */
			}
		}

		if (leftv) {
			i__1 = *ihi;
			for (i__ = *ilo; i__ <= i__1; ++i__) {
				s = 1. / scale[i__];
				NUMblas_dscal (m, &s, &v_ref (i__, 1), ldv);
				/* L20: */
			}
		}

	}

	/* Backward permutation

	   For I = ILO-1 step -1 until 1, IHI+1 step 1 until N do -- */

L30:
	if (lsame_ (job, "P") || lsame_ (job, "B")) {
		if (rightv) {
			i__1 = *n;
			for (ii = 1; ii <= i__1; ++ii) {
				i__ = ii;
				if (i__ >= *ilo && i__ <= *ihi) {
					goto L40;
				}
				if (i__ < *ilo) {
					i__ = *ilo - ii;
				}
				k = (integer) scale[i__];
				if (k == i__) {
					goto L40;
				}
				NUMblas_dswap (m, &v_ref (i__, 1), ldv, &v_ref (k, 1), ldv);
L40:
				;
			}
		}

		if (leftv) {
			i__1 = *n;
			for (ii = 1; ii <= i__1; ++ii) {
				i__ = ii;
				if (i__ >= *ilo && i__ <= *ihi) {
					goto L50;
				}
				if (i__ < *ilo) {
					i__ = *ilo - ii;
				}
				k = (integer) scale[i__];
				if (k == i__) {
					goto L50;
				}
				NUMblas_dswap (m, &v_ref (i__, 1), ldv, &v_ref (k, 1), ldv);
L50:
				;
			}
		}
	}

	return 0;
}								/* NUMlapack_dgebak */

#undef v_ref

int NUMlapack_dgebal (const char *job, integer *n, double *a, integer *lda, integer *ilo, integer *ihi, double *scale,
                      integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2;
	double d__1, d__2;

	/* Local variables */
	static integer iexc;
	static double c__, f, g;
	static integer i__, j, k, l, m;
	static double r__, s;
	static double sfmin1, sfmin2, sfmax1, sfmax2, ca, ra;
	static int noconv;
	static integer ica, ira;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--scale;

	/* Function Body */
	*info = 0;
	if (!lsame_ (job, "N") && !lsame_ (job, "P") && !lsame_ (job, "S") && !lsame_ (job, "B")) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *n)) {
		*info = -4;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("NUMlapack_dgebal ", &i__1);
		return 0;
	}

	k = 1;
	l = *n;

	if (*n == 0) {
		goto L210;
	}

	if (lsame_ (job, "N")) {
		i__1 = *n;
		for (i__ = 1; i__ <= i__1; ++i__) {
			scale[i__] = 1.;
			/* L10: */
		}
		goto L210;
	}

	if (lsame_ (job, "S")) {
		goto L120;
	}

	/* Permutation to isolate eigenvalues if possible */

	goto L50;

	/* Row and column exchange. */

L20:
	scale[m] = (double) j;
	if (j == m) {
		goto L30;
	}

	NUMblas_dswap (&l, &a_ref (1, j), &c__1, &a_ref (1, m), &c__1);
	i__1 = *n - k + 1;
	NUMblas_dswap (&i__1, &a_ref (j, k), lda, &a_ref (m, k), lda);

L30:
	switch (iexc) {
		case 1:
			goto L40;
		case 2:
			goto L80;
	}

	/* Search for rows isolating an eigenvalue and push them down. */

L40:
	if (l == 1) {
		goto L210;
	}
	--l;

L50:
	for (j = l; j >= 1; --j) {

		i__1 = l;
		for (i__ = 1; i__ <= i__1; ++i__) {
			if (i__ == j) {
				goto L60;
			}
			if (a_ref (j, i__) != 0.) {
				goto L70;
			}
L60:
			;
		}

		m = l;
		iexc = 1;
		goto L20;
L70:
		;
	}

	goto L90;

	/* Search for columns isolating an eigenvalue and push them left. */

L80:
	++k;

L90:
	i__1 = l;
	for (j = k; j <= i__1; ++j) {

		i__2 = l;
		for (i__ = k; i__ <= i__2; ++i__) {
			if (i__ == j) {
				goto L100;
			}
			if (a_ref (i__, j) != 0.) {
				goto L110;
			}
L100:
			;
		}

		m = k;
		iexc = 2;
		goto L20;
L110:
		;
	}

L120:
	i__1 = l;
	for (i__ = k; i__ <= i__1; ++i__) {
		scale[i__] = 1.;
		/* L130: */
	}

	if (lsame_ (job, "P")) {
		goto L210;
	}

	/* Balance the submatrix in rows K to L.

	   Iterative loop for norm reduction */

	sfmin1 = NUMblas_dlamch ("S") / NUMblas_dlamch ("P");
	sfmax1 = 1. / sfmin1;
	sfmin2 = sfmin1 * 8.;
	sfmax2 = 1. / sfmin2;
L140:
	noconv = FALSE;

	i__1 = l;
	for (i__ = k; i__ <= i__1; ++i__) {
		c__ = 0.;
		r__ = 0.;

		i__2 = l;
		for (j = k; j <= i__2; ++j) {
			if (j == i__) {
				goto L150;
			}
			c__ += (d__1 = a_ref (j, i__), fabs (d__1));
			r__ += (d__1 = a_ref (i__, j), fabs (d__1));
L150:
			;
		}
		ica = NUMblas_idamax (&l, &a_ref (1, i__), &c__1);
		ca = (d__1 = a_ref (ica, i__), fabs (d__1));
		i__2 = *n - k + 1;
		ira = NUMblas_idamax (&i__2, &a_ref (i__, k), lda);
		ra = (d__1 = a_ref (i__, ira + k - 1), fabs (d__1));

		/* Guard against zero C or R due to underflow. */

		if (c__ == 0. || r__ == 0.) {
			goto L200;
		}
		g = r__ / 8.;
		f = 1.;
		s = c__ + r__;
L160:
		/* Computing MAX */
		d__1 = MAX (f, c__);
		/* Computing MIN */
		d__2 = MIN (r__, g);
		if (c__ >= g || MAX (d__1, ca) >= sfmax2 || MIN (d__2, ra) <= sfmin2) {
			goto L170;
		}
		f *= 8.;
		c__ *= 8.;
		ca *= 8.;
		r__ /= 8.;
		g /= 8.;
		ra /= 8.;
		goto L160;

L170:
		g = c__ / 8.;
L180:
		/* Computing MIN */
		d__1 = MIN (f, c__), d__1 = MIN (d__1, g);
		if (g < r__ || MAX (r__, ra) >= sfmax2 || MIN (d__1, ca) <= sfmin2) {
			goto L190;
		}
		f /= 8.;
		c__ /= 8.;
		g /= 8.;
		ca /= 8.;
		r__ *= 8.;
		ra *= 8.;
		goto L180;

		/* Now balance. */

L190:
		if (c__ + r__ >= s * .95) {
			goto L200;
		}
		if (f < 1. && scale[i__] < 1.) {
			if (f * scale[i__] <= sfmin1) {
				goto L200;
			}
		}
		if (f > 1. && scale[i__] > 1.) {
			if (scale[i__] >= sfmax1 / f) {
				goto L200;
			}
		}
		g = 1. / f;
		scale[i__] *= f;
		noconv = TRUE;

		i__2 = *n - k + 1;
		NUMblas_dscal (&i__2, &g, &a_ref (i__, k), lda);
		NUMblas_dscal (&l, &f, &a_ref (1, i__), &c__1);

L200:
		;
	}

	if (noconv) {
		goto L140;
	}

L210:
	*ilo = k;
	*ihi = l;

	return 0;
}								/* NUMlapack_dgebal */


int NUMlapack_dgebrd (integer *m, integer *n, double *a, integer *lda, double *d__, double *e, double *tauq,
                      double *taup, double *work, integer *lwork, integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;
	static integer c_n1 = -1;
	static integer c__3 = 3;
	static integer c__2 = 2;
	static double c_b21 = -1.;
	static double c_b22 = 1.;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3, i__4;

	/* Local variables */
	static integer i__, j;
	static integer nbmin, iinfo, minmn;
	static integer nb;
	static integer nx;
	static double ws;
	static integer ldwrkx, ldwrky, lwkopt;
	static integer lquery;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--d__;
	--e;
	--tauq;
	--taup;
	--work;

	/* Function Body */
	*info = 0;
	/* Computing MAX */
	i__1 = 1, i__2 = NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
	nb = MAX (i__1, i__2);
	lwkopt = (*m + *n) * nb;
	work[1] = (double) lwkopt;
	lquery = *lwork == -1;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	} else {					/* if(complicated condition) */
		/* Computing MAX */
		i__1 = MAX (1, *m);
		if (*lwork < MAX (i__1, *n) && !lquery) {
			*info = -10;
		}
	}
	if (*info < 0) {
		i__1 = - (*info);
		xerbla_ ("DGEBRD", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	minmn = MIN (*m, *n);
	if (minmn == 0) {
		work[1] = 1.;
		return 0;
	}

	ws = (double) MAX (*m, *n);
	ldwrkx = *m;
	ldwrky = *n;

	if (nb > 1 && nb < minmn) {

		/* Set the crossover point NX.

		   Computing MAX */
		i__1 = nb, i__2 = NUMlapack_ilaenv (&c__3, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
		nx = MAX (i__1, i__2);

		/* Determine when to switch from blocked to unblocked code. */

		if (nx < minmn) {
			ws = (double) ( (*m + *n) * nb);
			if ( (double) (*lwork) < ws) {

				/* Not enough work space for the optimal NB, consider using a
				   smaller block size. */

				nbmin = NUMlapack_ilaenv (&c__2, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
				if (*lwork >= (*m + *n) * nbmin) {
					nb = *lwork / (*m + *n);
				} else {
					nb = 1;
					nx = minmn;
				}
			}
		}
	} else {
		nx = minmn;
	}

	i__1 = minmn - nx;
	i__2 = nb;
	for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {

		/* Reduce rows and columns i:i+nb-1 to bidiagonal form and return the
		   matrices X and Y which are needed to update the unreduced part of
		   the matrix */

		i__3 = *m - i__ + 1;
		i__4 = *n - i__ + 1;
		NUMlapack_dlabrd (&i__3, &i__4, &nb, &a_ref (i__, i__), lda, &d__[i__], &e[i__], &tauq[i__], &taup[i__],
		                  &work[1], &ldwrkx, &work[ldwrkx * nb + 1], &ldwrky);

		/* Update the trailing submatrix A(i+nb:m,i+nb:n), using an update of
		   the form A := A - V*Y' - X*U' */

		i__3 = *m - i__ - nb + 1;
		i__4 = *n - i__ - nb + 1;
		NUMblas_dgemm ("No transpose", "Transpose", &i__3, &i__4, &nb, &c_b21, &a_ref (i__ + nb, i__), lda,
		               &work[ldwrkx * nb + nb + 1], &ldwrky, &c_b22, &a_ref (i__ + nb, i__ + nb), lda);
		i__3 = *m - i__ - nb + 1;
		i__4 = *n - i__ - nb + 1;
		NUMblas_dgemm ("No transpose", "No transpose", &i__3, &i__4, &nb, &c_b21, &work[nb + 1], &ldwrkx,
		               &a_ref (i__, i__ + nb), lda, &c_b22, &a_ref (i__ + nb, i__ + nb), lda);

		/* Copy diagonal and off-diagonal elements of B back into A */

		if (*m >= *n) {
			i__3 = i__ + nb - 1;
			for (j = i__; j <= i__3; ++j) {
				a_ref (j, j) = d__[j];
				a_ref (j, j + 1) = e[j];
				/* L10: */
			}
		} else {
			i__3 = i__ + nb - 1;
			for (j = i__; j <= i__3; ++j) {
				a_ref (j, j) = d__[j];
				a_ref (j + 1, j) = e[j];
				/* L20: */
			}
		}
		/* L30: */
	}

	/* Use unblocked code to reduce the remainder of the matrix */

	i__2 = *m - i__ + 1;
	i__1 = *n - i__ + 1;
	NUMlapack_dgebd2 (&i__2, &i__1, &a_ref (i__, i__), lda, &d__[i__], &e[i__], &tauq[i__], &taup[i__],
	                  &work[1], &iinfo);
	work[1] = ws;
	return 0;
}								/* NUMlapack_dgebrd */


int NUMlapack_dgeev (const char *jobvl, const char *jobvr, integer *n, double *a, integer *lda, double *wr, double *wi,
                     double *vl, integer *ldvl, double *vr, integer *ldvr, double *work, integer *lwork, integer *info) {
	/* Table of constant values */
	static integer c__8 = 8;
	static integer c_n1 = -1;

	/* System generated locals */
	integer a_dim1, a_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3, i__4;
	double d__1, d__2;

	/* Local variables */
	static integer ibal;
	static char side[1];
	static integer maxb;
	static double anrm;
	static integer ierr, itau;
	static integer iwrk, nout;
	static integer i__, k;
	static double r__;
	static double cs;
	static int scalea;
	static double cscale;
	static double sn;
	static int select[1];
	static double bignum;
	static integer minwrk, maxwrk;
	static int wantvl;
	static double smlnum;
	static integer hswork;
	static int lquery, wantvr;
	static integer ihi;
	static double scl;
	static integer ilo;
	static double dum[1], eps;

#define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
#define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--wr;
	--wi;
	vl_dim1 = *ldvl;
	vl_offset = 1 + vl_dim1 * 1;
	vl -= vl_offset;
	vr_dim1 = *ldvr;
	vr_offset = 1 + vr_dim1 * 1;
	vr -= vr_offset;
	--work;

	/* Function Body */
	*info = 0;
	lquery = *lwork == -1;
	wantvl = lsame_ (jobvl, "V");
	wantvr = lsame_ (jobvr, "V");
	if (!wantvl && !lsame_ (jobvl, "N")) {
		*info = -1;
	} else if (!wantvr && !lsame_ (jobvr, "N")) {
		*info = -2;
	} else if (*n < 0) {
		*info = -3;
	} else if (*lda < MAX (1, *n)) {
		*info = -5;
	} else if (*ldvl < 1 || wantvl && *ldvl < *n) {
		*info = -9;
	} else if (*ldvr < 1 || wantvr && *ldvr < *n) {
		*info = -11;
	}

	/* Compute workspace (Note: Comments in the code beginning "Workspace:"
	   describe the minimal amount of workspace needed at that point in the
	   code, as well as the preferred amount for good performance. NB refers
	   to the optimal block size for the immediately following subroutine, as
	   returned by ILAENV. HSWORK refers to the workspace preferred by
	   NUMlapack_dhseqr , as calculated below. HSWORK is computed assuming
	   ILO=1 and IHI=N, the worst case.) */

	minwrk = 1;
	if (*info == 0 && (*lwork >= 1 || lquery)) {
		maxwrk =
		    (*n << 1) + *n * NUMlapack_ilaenv (&c__1, "NUMlapack_dgehrd ", " ", n, &c__1, n, &c__0, 6, 1);
		if (!wantvl && !wantvr) {
			/* Computing MAX */
			i__1 = 1, i__2 = *n * 3;
			minwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = NUMlapack_ilaenv (&c__8, "NUMlapack_dhseqr ", "EN", n, &c__1, n, &c_n1, 6, 2);
			maxb = MAX (i__1, 2);
			/* Computing MIN Computing MAX */
			i__3 = 2, i__4 = NUMlapack_ilaenv (&c__4, "NUMlapack_dhseqr ", "EN", n, &c__1, n, &c_n1, 6, 2);
			i__1 = MIN (maxb, *n), i__2 = MAX (i__3, i__4);
			k = MIN (i__1, i__2);
			/* Computing MAX */
			i__1 = k * (k + 2), i__2 = *n << 1;
			hswork = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 = *n + 1, i__1 = MAX (i__1, i__2), i__2 = *n + hswork;
			maxwrk = MAX (i__1, i__2);
		} else {
			/* Computing MAX */
			i__1 = 1, i__2 = *n << 2;
			minwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 =
			           (*n << 1) + (*n - 1) * NUMlapack_ilaenv (&c__1, "DOR" "GHR", " ", n, &c__1, n, &c_n1, 6, 1);
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = NUMlapack_ilaenv (&c__8, "NUMlapack_dhseqr ", "SV", n, &c__1, n, &c_n1, 6, 2);
			maxb = MAX (i__1, 2);
			/* Computing MIN Computing MAX */
			i__3 = 2, i__4 = NUMlapack_ilaenv (&c__4, "NUMlapack_dhseqr ", "SV", n, &c__1, n, &c_n1, 6, 2);
			i__1 = MIN (maxb, *n), i__2 = MAX (i__3, i__4);
			k = MIN (i__1, i__2);
			/* Computing MAX */
			i__1 = k * (k + 2), i__2 = *n << 1;
			hswork = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 = *n + 1, i__1 = MAX (i__1, i__2), i__2 = *n + hswork;
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 = *n << 2;
			maxwrk = MAX (i__1, i__2);
		}
		work[1] = (double) maxwrk;
	}
	if (*lwork < minwrk && !lquery) {
		*info = -13;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("NUMlapack_dgeev ", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	if (*n == 0) {
		return 0;
	}

	/* Get machine constants */

	eps = NUMblas_dlamch ("P");
	smlnum = NUMblas_dlamch ("S");
	bignum = 1. / smlnum;
	NUMlapack_dlabad (&smlnum, &bignum);
	smlnum = sqrt (smlnum) / eps;
	bignum = 1. / smlnum;

	/* Scale A if max element outside range [SMLNUM,BIGNUM] */

	anrm = NUMlapack_dlange ("M", n, n, &a[a_offset], lda, dum);
	scalea = FALSE;
	if (anrm > 0. && anrm < smlnum) {
		scalea = TRUE;
		cscale = smlnum;
	} else if (anrm > bignum) {
		scalea = TRUE;
		cscale = bignum;
	}
	if (scalea) {
		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &ierr);
	}

	/* Balance the matrix (Workspace: need N) */

	ibal = 1;
	NUMlapack_dgebal ("B", n, &a[a_offset], lda, &ilo, &ihi, &work[ibal], &ierr);

	/* Reduce to upper Hessenberg form (Workspace: need 3*N, prefer 2*N+N*NB)
	 */

	itau = ibal + *n;
	iwrk = itau + *n;
	i__1 = *lwork - iwrk + 1;
	NUMlapack_dgehrd (n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1, &ierr);

	if (wantvl) {

		/* Want left eigenvectors Copy Householder vectors to VL */

		* (unsigned char *) side = 'L';
		NUMlapack_dlacpy ("L", n, n, &a[a_offset], lda, &vl[vl_offset], ldvl);

		/* Generate orthogonal matrix in VL (Workspace: need 3*N-1, prefer
		   2*N+(N-1)*NB) */

		i__1 = *lwork - iwrk + 1;
		NUMlapack_dorghr (n, &ilo, &ihi, &vl[vl_offset], ldvl, &work[itau], &work[iwrk], &i__1, &ierr);

		/* Perform QR iteration, accumulating Schur vectors in VL (Workspace:
		   need N+1, prefer N+HSWORK (see comments) ) */

		iwrk = itau;
		i__1 = *lwork - iwrk + 1;
		NUMlapack_dhseqr ("S", "V", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vl[vl_offset], ldvl,
		                  &work[iwrk], &i__1, info);

		if (wantvr) {

			/* Want left and right eigenvectors Copy Schur vectors to VR */

			* (unsigned char *) side = 'B';
			NUMlapack_dlacpy ("F", n, n, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr);
		}

	} else if (wantvr) {

		/* Want right eigenvectors Copy Householder vectors to VR */

		* (unsigned char *) side = 'R';
		NUMlapack_dlacpy ("L", n, n, &a[a_offset], lda, &vr[vr_offset], ldvr);

		/* Generate orthogonal matrix in VR (Workspace: need 3*N-1, prefer
		   2*N+(N-1)*NB) */

		i__1 = *lwork - iwrk + 1;
		NUMlapack_dorghr (n, &ilo, &ihi, &vr[vr_offset], ldvr, &work[itau], &work[iwrk], &i__1, &ierr);

		/* Perform QR iteration, accumulating Schur vectors in VR (Workspace:
		   need N+1, prefer N+HSWORK (see comments) ) */

		iwrk = itau;
		i__1 = *lwork - iwrk + 1;
		NUMlapack_dhseqr ("S", "V", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset], ldvr,
		                  &work[iwrk], &i__1, info);

	} else {

		/* Compute eigenvalues only (Workspace: need N+1, prefer N+HSWORK
		   (see comments) ) */

		iwrk = itau;
		i__1 = *lwork - iwrk + 1;
		NUMlapack_dhseqr ("E", "N", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset], ldvr,
		                  &work[iwrk], &i__1, info);
	}

	/* If INFO > 0 from NUMlapack_dhseqr , then quit */

	if (*info > 0) {
		goto L50;
	}

	if (wantvl || wantvr) {

		/* Compute left and/or right eigenvectors (Workspace: need 4*N) */

		NUMlapack_dtrevc (side, "B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl, &vr[vr_offset],
		                  ldvr, n, &nout, &work[iwrk], &ierr);
	}

	if (wantvl) {

		/* Undo balancing of left eigenvectors (Workspace: need N) */

		NUMlapack_dgebak ("B", "L", n, &ilo, &ihi, &work[ibal], n, &vl[vl_offset], ldvl, &ierr);

		/* Normalize left eigenvectors and make largest component real */

		i__1 = *n;
		for (i__ = 1; i__ <= i__1; ++i__) {
			if (wi[i__] == 0.) {
				scl = 1. / NUMblas_dnrm2 (n, &vl_ref (1, i__), &c__1);
				NUMblas_dscal (n, &scl, &vl_ref (1, i__), &c__1);
			} else if (wi[i__] > 0.) {
				d__1 = NUMblas_dnrm2 (n, &vl_ref (1, i__), &c__1);
				d__2 = NUMblas_dnrm2 (n, &vl_ref (1, i__ + 1), &c__1);
				scl = 1. / NUMlapack_dlapy2 (&d__1, &d__2);
				NUMblas_dscal (n, &scl, &vl_ref (1, i__), &c__1);
				NUMblas_dscal (n, &scl, &vl_ref (1, i__ + 1), &c__1);
				i__2 = *n;
				for (k = 1; k <= i__2; ++k) {
					/* Computing 2nd power */
					d__1 = vl_ref (k, i__);
					/* Computing 2nd power */
					d__2 = vl_ref (k, i__ + 1);
					work[iwrk + k - 1] = d__1 * d__1 + d__2 * d__2;
					/* L10: */
				}
				k = NUMblas_idamax (n, &work[iwrk], &c__1);
				NUMlapack_dlartg (&vl_ref (k, i__), &vl_ref (k, i__ + 1), &cs, &sn, &r__);
				NUMblas_drot (n, &vl_ref (1, i__), &c__1, &vl_ref (1, i__ + 1), &c__1, &cs, &sn);
				vl_ref (k, i__ + 1) = 0.;
			}
			/* L20: */
		}
	}

	if (wantvr) {

		/* Undo balancing of right eigenvectors (Workspace: need N) */

		NUMlapack_dgebak ("B", "R", n, &ilo, &ihi, &work[ibal], n, &vr[vr_offset], ldvr, &ierr);

		/* Normalize right eigenvectors and make largest component real */

		i__1 = *n;
		for (i__ = 1; i__ <= i__1; ++i__) {
			if (wi[i__] == 0.) {
				scl = 1. / NUMblas_dnrm2 (n, &vr_ref (1, i__), &c__1);
				NUMblas_dscal (n, &scl, &vr_ref (1, i__), &c__1);
			} else if (wi[i__] > 0.) {
				d__1 = NUMblas_dnrm2 (n, &vr_ref (1, i__), &c__1);
				d__2 = NUMblas_dnrm2 (n, &vr_ref (1, i__ + 1), &c__1);
				scl = 1. / NUMlapack_dlapy2 (&d__1, &d__2);
				NUMblas_dscal (n, &scl, &vr_ref (1, i__), &c__1);
				NUMblas_dscal (n, &scl, &vr_ref (1, i__ + 1), &c__1);
				i__2 = *n;
				for (k = 1; k <= i__2; ++k) {
					/* Computing 2nd power */
					d__1 = vr_ref (k, i__);
					/* Computing 2nd power */
					d__2 = vr_ref (k, i__ + 1);
					work[iwrk + k - 1] = d__1 * d__1 + d__2 * d__2;
					/* L30: */
				}
				k = NUMblas_idamax (n, &work[iwrk], &c__1);
				NUMlapack_dlartg (&vr_ref (k, i__), &vr_ref (k, i__ + 1), &cs, &sn, &r__);
				NUMblas_drot (n, &vr_ref (1, i__), &c__1, &vr_ref (1, i__ + 1), &c__1, &cs, &sn);
				vr_ref (k, i__ + 1) = 0.;
			}
			/* L40: */
		}
	}

	/* Undo scaling if necessary */

L50:
	if (scalea) {
		i__1 = *n - *info;
		/* Computing MAX */
		i__3 = *n - *info;
		i__2 = MAX (i__3, 1);
		NUMlapack_dlascl ("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wr[*info + 1], &i__2, &ierr);
		i__1 = *n - *info;
		/* Computing MAX */
		i__3 = *n - *info;
		i__2 = MAX (i__3, 1);
		NUMlapack_dlascl ("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[*info + 1], &i__2, &ierr);
		if (*info > 0) {
			i__1 = ilo - 1;
			NUMlapack_dlascl ("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wr[1], n, &ierr);
			i__1 = ilo - 1;
			NUMlapack_dlascl ("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[1], n, &ierr);
		}
	}

	work[1] = (double) maxwrk;
	return 0;
}								/* NUMlapack_dgeev */

#undef vr_ref
#undef vl_ref


int NUMlapack_dgehd2 (integer *n, integer *ilo, integer *ihi, double *a, integer *lda, double *tau, double *work,
                      integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3;

	/* Local variables */
	static integer i__;
	static double aii;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	if (*n < 0) {
		*info = -1;
	} else if (*ilo < 1 || *ilo > MAX (1, *n)) {
		*info = -2;
	} else if (*ihi < MIN (*ilo, *n) || *ihi > *n) {
		*info = -3;
	} else if (*lda < MAX (1, *n)) {
		*info = -5;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("NUMlapack_dgehd2 ", &i__1);
		return 0;
	}

	i__1 = *ihi - 1;
	for (i__ = *ilo; i__ <= i__1; ++i__) {

		/* Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)

		   Computing MIN */
		i__2 = i__ + 2;
		i__3 = *ihi - i__;
		NUMlapack_dlarfg (&i__3, &a_ref (i__ + 1, i__), &a_ref (MIN (i__2, *n), i__), &c__1, &tau[i__]);
		aii = a_ref (i__ + 1, i__);
		a_ref (i__ + 1, i__) = 1.;

		/* Apply H(i) to A(1:ihi,i+1:ihi) from the right */

		i__2 = *ihi - i__;
		NUMlapack_dlarf ("Right", ihi, &i__2, &a_ref (i__ + 1, i__), &c__1, &tau[i__], &a_ref (1, i__ + 1),
		                 lda, &work[1]);

		/* Apply H(i) to A(i+1:ihi,i+1:n) from the left */

		i__2 = *ihi - i__;
		i__3 = *n - i__;
		NUMlapack_dlarf ("Left", &i__2, &i__3, &a_ref (i__ + 1, i__), &c__1, &tau[i__], &a_ref (i__ + 1,
		                 i__ + 1), lda, &work[1]);

		a_ref (i__ + 1, i__) = aii;
		/* L10: */
	}

	return 0;
}								/* NUMlapack_dgehd2 */


int NUMlapack_dgehrd (integer *n, integer *ilo, integer *ihi, double *a, integer *lda, double *tau, double *work,
                      integer *lwork, integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;
	static integer c_n1 = -1;
	static integer c__3 = 3;
	static integer c__2 = 2;
	static integer c__65 = 65;
	static double c_b25 = -1.;
	static double c_b26 = 1.;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3, i__4;

	/* Local variables */
	static integer i__;
	static double t[4160] /* was [65][64] */ ;
	static integer nbmin, iinfo;
	static integer ib;
	static double ei;
	static integer nb, nh;
	static integer nx;
	static integer ldwork, lwkopt;
	static int lquery;
	static integer iws;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	/* Computing MIN */
	i__1 = 64, i__2 = NUMlapack_ilaenv (&c__1, "NUMlapack_dgehrd ", " ", n, ilo, ihi, &c_n1, 6, 1);
	nb = MIN (i__1, i__2);
	lwkopt = *n * nb;
	work[1] = (double) lwkopt;
	lquery = *lwork == -1;
	if (*n < 0) {
		*info = -1;
	} else if (*ilo < 1 || *ilo > MAX (1, *n)) {
		*info = -2;
	} else if (*ihi < MIN (*ilo, *n) || *ihi > *n) {
		*info = -3;
	} else if (*lda < MAX (1, *n)) {
		*info = -5;
	} else if (*lwork < MAX (1, *n) && !lquery) {
		*info = -8;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("NUMlapack_dgehrd ", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Set elements 1:ILO-1 and IHI:N-1 of TAU to zero */

	i__1 = *ilo - 1;
	for (i__ = 1; i__ <= i__1; ++i__) {
		tau[i__] = 0.;
		/* L10: */
	}
	i__1 = *n - 1;
	for (i__ = MAX (1, *ihi); i__ <= i__1; ++i__) {
		tau[i__] = 0.;
		/* L20: */
	}

	/* Quick return if possible */

	nh = *ihi - *ilo + 1;
	if (nh <= 1) {
		work[1] = 1.;
		return 0;
	}

	/* Determine the block size.

	   Computing MIN */
	i__1 = 64, i__2 = NUMlapack_ilaenv (&c__1, "NUMlapack_dgehrd ", " ", n, ilo, ihi, &c_n1, 6, 1);
	nb = MIN (i__1, i__2);
	nbmin = 2;
	iws = 1;
	if (nb > 1 && nb < nh) {

		/* Determine when to cross over from blocked to unblocked code (last
		   block is always handled by unblocked code).

		   Computing MAX */
		i__1 = nb, i__2 = NUMlapack_ilaenv (&c__3, "NUMlapack_dgehrd ", " ", n, ilo, ihi, &c_n1, 6, 1);
		nx = MAX (i__1, i__2);
		if (nx < nh) {

			/* Determine if workspace is large enough for blocked code. */

			iws = *n * nb;
			if (*lwork < iws) {

				/* Not enough workspace to use optimal NB: determine the
				   minimum value of NB, and reduce NB or force use of
				   unblocked code.

				   Computing MAX */
				i__1 = 2, i__2 =
				           NUMlapack_ilaenv (&c__2, "NUMlapack_dgehrd ", " ", n, ilo, ihi, &c_n1, 6, 1);
				nbmin = MAX (i__1, i__2);
				if (*lwork >= *n * nbmin) {
					nb = *lwork / *n;
				} else {
					nb = 1;
				}
			}
		}
	}
	ldwork = *n;

	if (nb < nbmin || nb >= nh) {

		/* Use unblocked code below */

		i__ = *ilo;

	} else {

		/* Use blocked code */

		i__1 = *ihi - 1 - nx;
		i__2 = nb;
		for (i__ = *ilo; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
			/* Computing MIN */
			i__3 = nb, i__4 = *ihi - i__;
			ib = MIN (i__3, i__4);

			/* Reduce columns i:i+ib-1 to Hessenberg form, returning the
			   matrices V and T of the block reflector H = I - V*T*V' which
			   performs the reduction, and also the matrix Y = A*V*T */

			NUMlapack_dlahrd (ihi, &i__, &ib, &a_ref (1, i__), lda, &tau[i__], t, &c__65, &work[1], &ldwork);

			/* Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
			   right, computing A := A - Y * V'. V(i+ib,ib-1) should be set to
			   1. */

			ei = a_ref (i__ + ib, i__ + ib - 1);
			a_ref (i__ + ib, i__ + ib - 1) = 1.;
			i__3 = *ihi - i__ - ib + 1;
			NUMblas_dgemm ("No transpose", "Transpose", ihi, &i__3, &ib, &c_b25, &work[1], &ldwork,
			               &a_ref (i__ + ib, i__), lda, &c_b26, &a_ref (1, i__ + ib), lda);
			a_ref (i__ + ib, i__ + ib - 1) = ei;

			/* Apply the block reflector H to A(i+1:ihi,i+ib:n) from the left
			 */

			i__3 = *ihi - i__;
			i__4 = *n - i__ - ib + 1;
			NUMlapack_dlarfb ("Left", "Transpose", "Forward", "Columnwise", &i__3, &i__4, &ib,
			                  &a_ref (i__ + 1, i__), lda, t, &c__65, &a_ref (i__ + 1, i__ + ib), lda, &work[1], &ldwork);
			/* L30: */
		}
	}

	/* Use unblocked code to reduce the rest of the matrix */

	NUMlapack_dgehd2 (n, &i__, ihi, &a[a_offset], lda, &tau[1], &work[1], &iinfo);
	work[1] = (double) iws;

	return 0;
}								/* NUMlapack_dgehrd */

int NUMlapack_dgelq2 (integer *m, integer *n, double *a, integer *lda, double *tau, double *work, integer *info) {
	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3;

	/* Local variables */
	static integer i__, k;
	static double aii;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGELQ2", &i__1);
		return 0;
	}

	k = MIN (*m, *n);

	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {

		/* Generate elementary reflector H(i) to annihilate A(i,i+1:n)

		   Computing MIN */
		i__2 = i__ + 1;
		i__3 = *n - i__ + 1;
		NUMlapack_dlarfg (&i__3, &a_ref (i__, i__), &a_ref (i__, MIN (i__2, *n)), lda, &tau[i__]);
		if (i__ < *m) {

			/* Apply H(i) to A(i+1:m,i:n) from the right */

			aii = a_ref (i__, i__);
			a_ref (i__, i__) = 1.;
			i__2 = *m - i__;
			i__3 = *n - i__ + 1;
			NUMlapack_dlarf ("Right", &i__2, &i__3, &a_ref (i__, i__), lda, &tau[i__], &a_ref (i__ + 1, i__),
			                 lda, &work[1]);
			a_ref (i__, i__) = aii;
		}
		/* L10: */
	}
	return 0;
}								/* NUMlapack_dgelq2 */

int NUMlapack_dgelqf (integer *m, integer *n, double *a, integer *lda, double *tau, double *work, integer *lwork,
                      integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;
	static integer c_n1 = -1;
	static integer c__3 = 3;
	static integer c__2 = 2;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3, i__4;

	/* Local variables */
	static integer i__, k, nbmin, iinfo;
	static integer ib, nb;
	static integer nx;
	static integer ldwork, lwkopt;
	static integer lquery;
	static integer iws;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	nb = NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
	lwkopt = *m * nb;
	work[1] = (double) lwkopt;
	lquery = *lwork == -1;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	} else if (*lwork < MAX (1, *m) && !lquery) {
		*info = -7;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGELQF", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	k = MIN (*m, *n);
	if (k == 0) {
		work[1] = 1.;
		return 0;
	}

	nbmin = 2;
	nx = 0;
	iws = *m;
	if (nb > 1 && nb < k) {

		/* Determine when to cross over from blocked to unblocked code.

		   Computing MAX */
		i__1 = 0, i__2 = NUMlapack_ilaenv (&c__3, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
		nx = MAX (i__1, i__2);
		if (nx < k) {

			/* Determine if workspace is large enough for blocked code. */

			ldwork = *m;
			iws = ldwork * nb;
			if (*lwork < iws) {

				/* Not enough workspace to use optimal NB: reduce NB and
				   determine the minimum value of NB. */

				nb = *lwork / ldwork;
				/* Computing MAX */
				i__1 = 2, i__2 = NUMlapack_ilaenv (&c__2, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
				nbmin = MAX (i__1, i__2);
			}
		}
	}

	if (nb >= nbmin && nb < k && nx < k) {

		/* Use blocked code initially */

		i__1 = k - nx;
		i__2 = nb;
		for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
			/* Computing MIN */
			i__3 = k - i__ + 1;
			ib = MIN (i__3, nb);

			/* Compute the LQ factorization of the current block
			   A(i:i+ib-1,i:n) */

			i__3 = *n - i__ + 1;
			NUMlapack_dgelq2 (&ib, &i__3, &a_ref (i__, i__), lda, &tau[i__], &work[1], &iinfo);
			if (i__ + ib <= *m) {

				/* Form the triangular factor of the block reflector H = H(i)
				   H(i+1) . . . H(i+ib-1) */

				i__3 = *n - i__ + 1;
				NUMlapack_dlarft ("Forward", "Rowwise", &i__3, &ib, &a_ref (i__, i__), lda, &tau[i__], &work[1],
				                  &ldwork);

				/* Apply H to A(i+ib:m,i:n) from the right */

				i__3 = *m - i__ - ib + 1;
				i__4 = *n - i__ + 1;
				NUMlapack_dlarfb ("Right", "No transpose", "Forward", "Rowwise", &i__3, &i__4, &ib, &a_ref (i__,
				                  i__), lda, &work[1], &ldwork, &a_ref (i__ + ib, i__), lda, &work[ib + 1], &ldwork);
			}
			/* L10: */
		}
	} else {
		i__ = 1;
	}

	/* Use unblocked code to factor the last or only block. */

	if (i__ <= k) {
		i__2 = *m - i__ + 1;
		i__1 = *n - i__ + 1;
		NUMlapack_dgelq2 (&i__2, &i__1, &a_ref (i__, i__), lda, &tau[i__], &work[1], &iinfo);
	}

	work[1] = (double) iws;
	return 0;
}								/* NUMlapack_dgelqf */

#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]

int NUMlapack_dgelss (integer *m, integer *n, integer *nrhs, double *a, integer *lda, double *b, integer *ldb, double *s,
                      double *rcond, integer *rank, double *work, integer *lwork, integer *info) {
	/* System generated locals */
	integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
	double d__1;

	/* Local variables */
	static double anrm, bnrm;
	static integer itau;
	static double vdum[1];
	static integer i__;
	static integer iascl, ibscl;
	static integer chunk;
	static double sfmin;
	static integer minmn;
	static integer maxmn, itaup, itauq, mnthr, iwork;
	static integer bl, ie, il;
	static integer mm;
	static integer bdspac;
	static double bignum;
	static integer ldwork;
	static integer minwrk, maxwrk;
	static double smlnum;
	static integer lquery;
	static double eps, thr;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	b_dim1 = *ldb;
	b_offset = 1 + b_dim1 * 1;
	b -= b_offset;
	--s;
	--work;

	/* Function Body */
	*info = 0;
	minmn = MIN (*m, *n);
	maxmn = MAX (*m, *n);
	mnthr = NUMlapack_ilaenv (&c__6, "DGELSS", " ", m, n, nrhs, &c_n1, 6, 1);
	lquery = *lwork == -1;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*nrhs < 0) {
		*info = -3;
	} else if (*lda < MAX (1, *m)) {
		*info = -5;
	} else if (*ldb < MAX (1, maxmn)) {
		*info = -7;
	}

	/* Compute workspace (Note: Comments in the code beginning "Workspace:"
	   describe the minimal amount of workspace needed at that point in the
	   code, as well as the preferred amount for good performance. NB refers
	   to the optimal block size for the immediately following subroutine, as
	   returned by ILAENV.) */

	minwrk = 1;
	if (*info == 0 && (*lwork >= 1 || lquery)) {
		maxwrk = 0;
		mm = *m;
		if (*m >= *n && *m >= mnthr) {

			/* Path 1a - overdetermined, with many more rows than columns */

			mm = *n;
			/* Computing MAX */
			i__1 = maxwrk, i__2 = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 =
			           *n + *nrhs * NUMlapack_ilaenv (&c__1, "DORMQR", "LT", m, nrhs, n, &c_n1, 6, 2);
			maxwrk = MAX (i__1, i__2);
		}
		if (*m >= *n) {

			/* Path 1 - overdetermined or exactly determined

			   Compute workspace needed for DBDSQR

			   Computing MAX */
			i__1 = 1, i__2 = *n * 5;
			bdspac = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 =
			           *n * 3 + (mm + *n) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", &mm, n, &c_n1, &c_n1, 6, 1);
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 =
			           *n * 3 + *nrhs * NUMlapack_ilaenv (&c__1, "DORMBR", "QLT", &mm, nrhs, n, &c_n1, 6, 3);
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = maxwrk, i__2 =
			           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
			maxwrk = MAX (i__1, i__2);
			maxwrk = MAX (maxwrk, bdspac);
			/* Computing MAX */
			i__1 = maxwrk, i__2 = *n * *nrhs;
			maxwrk = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = *n * 3 + mm, i__2 = *n * 3 + *nrhs, i__1 = MAX (i__1, i__2);
			minwrk = MAX (i__1, bdspac);
			maxwrk = MAX (minwrk, maxwrk);
		}
		if (*n > *m) {

			/* Compute workspace needed for DBDSQR

			   Computing MAX */
			i__1 = 1, i__2 = *m * 5;
			bdspac = MAX (i__1, i__2);
			/* Computing MAX */
			i__1 = *m * 3 + *nrhs, i__2 = *m * 3 + *n, i__1 = MAX (i__1, i__2);
			minwrk = MAX (i__1, bdspac);
			if (*n >= mnthr) {

				/* Path 2a - underdetermined, with many more columns than
				   rows */

				maxwrk = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m * *m + (*m << 2) + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1,
				                   &c_n1, 6, 1);
				maxwrk = MAX (i__1, i__2);
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m * *m + (*m << 2) + *nrhs * NUMlapack_ilaenv (&c__1, "DORMBR", "QLT", m, nrhs, m, &c_n1,
				                   6, 3);
				maxwrk = MAX (i__1, i__2);
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m * *m + (*m << 2) + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6,
				                   1);
				maxwrk = MAX (i__1, i__2);
				/* Computing MAX */
				i__1 = maxwrk, i__2 = *m * *m + *m + bdspac;
				maxwrk = MAX (i__1, i__2);
				if (*nrhs > 1) {
					/* Computing MAX */
					i__1 = maxwrk, i__2 = *m * *m + *m + *m * *nrhs;
					maxwrk = MAX (i__1, i__2);
				} else {
					/* Computing MAX */
					i__1 = maxwrk, i__2 = *m * *m + (*m << 1);
					maxwrk = MAX (i__1, i__2);
				}
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m + *nrhs * NUMlapack_ilaenv (&c__1, "DORMLQ", "LT", n, nrhs, m, &c_n1, 6, 2);
				maxwrk = MAX (i__1, i__2);
			} else {

				/* Path 2 - underdetermined */

				maxwrk = *m * 3 + (*n + *m) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m * 3 + *nrhs * NUMlapack_ilaenv (&c__1, "DORMBR", "QLT", m, nrhs, m, &c_n1, 6, 3);
				maxwrk = MAX (i__1, i__2);
				/* Computing MAX */
				i__1 = maxwrk, i__2 =
				           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, n, m, &c_n1, 6, 1);
				maxwrk = MAX (i__1, i__2);
				maxwrk = MAX (maxwrk, bdspac);
				/* Computing MAX */
				i__1 = maxwrk, i__2 = *n * *nrhs;
				maxwrk = MAX (i__1, i__2);
			}
		}
		maxwrk = MAX (minwrk, maxwrk);
		work[1] = (double) maxwrk;
	}

	minwrk = MAX (minwrk, 1);
	if (*lwork < minwrk && !lquery) {
		*info = -12;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGELSS", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	if (*m == 0 || *n == 0) {
		*rank = 0;
		return 0;
	}

	/* Get machine parameters */

	eps = NUMblas_dlamch ("P");
	sfmin = NUMblas_dlamch ("S");
	smlnum = sfmin / eps;
	bignum = 1. / smlnum;
	NUMlapack_dlabad (&smlnum, &bignum);

	/* Scale A if max element outside range [SMLNUM,BIGNUM] */

	anrm = NUMlapack_dlange ("M", m, n, &a[a_offset], lda, &work[1]);
	iascl = 0;
	if (anrm > 0. && anrm < smlnum) {

		/* Scale matrix norm up to SMLNUM */

		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, info);
		iascl = 1;
	} else if (anrm > bignum) {

		/* Scale matrix norm down to BIGNUM */

		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, info);
		iascl = 2;
	} else if (anrm == 0.) {

		/* Matrix all zero. Return zero solution. */

		i__1 = MAX (*m, *n);
		NUMlapack_dlaset ("F", &i__1, nrhs, &c_b74, &c_b74, &b[b_offset], ldb);
		NUMlapack_dlaset ("F", &minmn, &c__1, &c_b74, &c_b74, &s[1], &c__1);
		*rank = 0;
		goto L70;
	}

	/* Scale B if max element outside range [SMLNUM,BIGNUM] */

	bnrm = NUMlapack_dlange ("M", m, nrhs, &b[b_offset], ldb, &work[1]);
	ibscl = 0;
	if (bnrm > 0. && bnrm < smlnum) {

		/* Scale matrix norm up to SMLNUM */

		NUMlapack_dlascl ("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, info);
		ibscl = 1;
	} else if (bnrm > bignum) {

		/* Scale matrix norm down to BIGNUM */

		NUMlapack_dlascl ("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, info);
		ibscl = 2;
	}

	/* Overdetermined case */

	if (*m >= *n) {

		/* Path 1 - overdetermined or exactly determined */

		mm = *m;
		if (*m >= mnthr) {

			/* Path 1a - overdetermined, with many more rows than columns */

			mm = *n;
			itau = 1;
			iwork = itau + *n;

			/* Compute A=Q*R (Workspace: need 2*N, prefer N+N*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__1, info);

			/* Multiply B by transpose(Q) (Workspace: need N+NRHS, prefer
			   N+NRHS*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dormqr ("L", "T", m, nrhs, n, &a[a_offset], lda, &work[itau], &b[b_offset], ldb,
			                  &work[iwork], &i__1, info);

			/* Zero out below R */

			if (*n > 1) {
				i__1 = *n - 1;
				i__2 = *n - 1;
				NUMlapack_dlaset ("L", &i__1, &i__2, &c_b74, &c_b74, &a_ref (2, 1), lda);
			}
		}

		ie = 1;
		itauq = ie + *n;
		itaup = itauq + *n;
		iwork = itaup + *n;

		/* Bidiagonalize R in A (Workspace: need 3*N+MM, prefer
		   3*N+(MM+N)*NB) */

		i__1 = *lwork - iwork + 1;
		NUMlapack_dgebrd (&mm, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup], &work[iwork],
		                  &i__1, info);

		/* Multiply B by transpose of left bidiagonalizing vectors of R
		   (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) */

		i__1 = *lwork - iwork + 1;
		NUMlapack_dormbr ("Q", "L", "T", &mm, nrhs, n, &a[a_offset], lda, &work[itauq], &b[b_offset], ldb,
		                  &work[iwork], &i__1, info);

		/* Generate right bidiagonalizing vectors of R in A (Workspace: need
		   4*N-1, prefer 3*N+(N-1)*NB) */

		i__1 = *lwork - iwork + 1;
		NUMlapack_dorgbr ("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__1, info);
		iwork = ie + *n;

		/* Perform bidiagonal QR iteration multiply B by transpose of left
		   singular vectors compute right singular vectors in A (Workspace:
		   need BDSPAC) */

		NUMlapack_dbdsqr ("U", n, n, &c__0, nrhs, &s[1], &work[ie], &a[a_offset], lda, vdum, &c__1,
		                  &b[b_offset], ldb, &work[iwork], info);
		if (*info != 0) {
			goto L70;
		}

		/* Multiply B by reciprocals of singular values

		   Computing MAX */
		d__1 = *rcond * s[1];
		thr = MAX (d__1, sfmin);
		if (*rcond < 0.) {
			/* Computing MAX */
			d__1 = eps * s[1];
			thr = MAX (d__1, sfmin);
		}
		*rank = 0;
		i__1 = *n;
		for (i__ = 1; i__ <= i__1; ++i__) {
			if (s[i__] > thr) {
				NUMlapack_drscl (nrhs, &s[i__], &b_ref (i__, 1), ldb);
				++ (*rank);
			} else {
				NUMlapack_dlaset ("F", &c__1, nrhs, &c_b74, &c_b74, &b_ref (i__, 1), ldb);
			}
			/* L10: */
		}

		/* Multiply B by right singular vectors (Workspace: need N, prefer
		   N*NRHS) */

		if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
			NUMblas_dgemm ("T", "N", n, nrhs, n, &c_b108, &a[a_offset], lda, &b[b_offset], ldb, &c_b74, &work[1],
			               ldb);
			NUMlapack_dlacpy ("G", n, nrhs, &work[1], ldb, &b[b_offset], ldb);
		} else if (*nrhs > 1) {
			chunk = *lwork / *n;
			i__1 = *nrhs;
			i__2 = chunk;
			for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
				/* Computing MIN */
				i__3 = *nrhs - i__ + 1;
				bl = MIN (i__3, chunk);
				NUMblas_dgemm ("T", "N", n, &bl, n, &c_b108, &a[a_offset], lda, &b_ref (1, i__), ldb, &c_b74,
				               &work[1], n);
				NUMlapack_dlacpy ("G", n, &bl, &work[1], n, &b_ref (1, i__), ldb);
				/* L20: */
			}
		} else {
			NUMblas_dgemv ("T", n, n, &c_b108, &a[a_offset], lda, &b[b_offset], &c__1, &c_b74, &work[1], &c__1);
			NUMblas_dcopy (n, &work[1], &c__1, &b[b_offset], &c__1);
		}

	} else {					/* if(complicated condition) */
		/* Computing MAX */
		i__2 = *m, i__1 = (*m << 1) - 4, i__2 = MAX (i__2, i__1), i__2 = MAX (i__2, *nrhs), i__1 =
		        *n - *m * 3;
		if (*n >= mnthr && *lwork >= (*m << 2) + *m * *m + MAX (i__2, i__1)) {

			/* Path 2a - underdetermined, with many more columns than rows
			   and sufficient workspace for an efficient algorithm */

			ldwork = *m;
			/* Computing MAX Computing MAX */
			i__3 = *m, i__4 = (*m << 1) - 4, i__3 = MAX (i__3, i__4), i__3 = MAX (i__3, *nrhs), i__4 =
			        *n - *m * 3;
			i__2 = (*m << 2) + *m * *lda + MAX (i__3, i__4), i__1 = *m * *lda + *m + *m * *nrhs;
			if (*lwork >= MAX (i__2, i__1)) {
				ldwork = *lda;
			}
			itau = 1;
			iwork = *m + 1;

			/* Compute A=L*Q (Workspace: need 2*M, prefer M+M*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, info);
			il = iwork;

			/* Copy L to WORK(IL), zeroing out above it */

			NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[il], &ldwork);
			i__2 = *m - 1;
			i__1 = *m - 1;
			NUMlapack_dlaset ("U", &i__2, &i__1, &c_b74, &c_b74, &work[il + ldwork], &ldwork);
			ie = il + ldwork * *m;
			itauq = ie + *m;
			itaup = itauq + *m;
			iwork = itaup + *m;

			/* Bidiagonalize L in WORK(IL) (Workspace: need M*M+5*M, prefer
			   M*M+4*M+2*M*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dgebrd (m, m, &work[il], &ldwork, &s[1], &work[ie], &work[itauq], &work[itaup],
			                  &work[iwork], &i__2, info);

			/* Multiply B by transpose of left bidiagonalizing vectors of L
			   (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dormbr ("Q", "L", "T", m, nrhs, m, &work[il], &ldwork, &work[itauq], &b[b_offset], ldb,
			                  &work[iwork], &i__2, info);

			/* Generate right bidiagonalizing vectors of R in WORK(IL)
			   (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dorgbr ("P", m, m, m, &work[il], &ldwork, &work[itaup], &work[iwork], &i__2, info);
			iwork = ie + *m;

			/* Perform bidiagonal QR iteration, computing right singular
			   vectors of L in WORK(IL) and multiplying B by transpose of
			   left singular vectors (Workspace: need M*M+M+BDSPAC) */

			NUMlapack_dbdsqr ("U", m, m, &c__0, nrhs, &s[1], &work[ie], &work[il], &ldwork, &a[a_offset], lda,
			                  &b[b_offset], ldb, &work[iwork], info);
			if (*info != 0) {
				goto L70;
			}

			/* Multiply B by reciprocals of singular values

			   Computing MAX */
			d__1 = *rcond * s[1];
			thr = MAX (d__1, sfmin);
			if (*rcond < 0.) {
				/* Computing MAX */
				d__1 = eps * s[1];
				thr = MAX (d__1, sfmin);
			}
			*rank = 0;
			i__2 = *m;
			for (i__ = 1; i__ <= i__2; ++i__) {
				if (s[i__] > thr) {
					NUMlapack_drscl (nrhs, &s[i__], &b_ref (i__, 1), ldb);
					++ (*rank);
				} else {
					NUMlapack_dlaset ("F", &c__1, nrhs, &c_b74, &c_b74, &b_ref (i__, 1), ldb);
				}
				/* L30: */
			}
			iwork = ie;

			/* Multiply B by right singular vectors of L in WORK(IL)
			   (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS) */

			if (*lwork >= *ldb * *nrhs + iwork - 1 && *nrhs > 1) {
				NUMblas_dgemm ("T", "N", m, nrhs, m, &c_b108, &work[il], &ldwork, &b[b_offset], ldb, &c_b74,
				               &work[iwork], ldb);
				NUMlapack_dlacpy ("G", m, nrhs, &work[iwork], ldb, &b[b_offset], ldb);
			} else if (*nrhs > 1) {
				chunk = (*lwork - iwork + 1) / *m;
				i__2 = *nrhs;
				i__1 = chunk;
				for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
					/* Computing MIN */
					i__3 = *nrhs - i__ + 1;
					bl = MIN (i__3, chunk);
					NUMblas_dgemm ("T", "N", m, &bl, m, &c_b108, &work[il], &ldwork, &b_ref (1, i__), ldb, &c_b74,
					               &work[iwork], n);
					NUMlapack_dlacpy ("G", m, &bl, &work[iwork], n, &b_ref (1, i__), ldb);
					/* L40: */
				}
			} else {
				NUMblas_dgemv ("T", m, m, &c_b108, &work[il], &ldwork, &b_ref (1, 1), &c__1, &c_b74, &work[iwork],
				               &c__1);
				NUMblas_dcopy (m, &work[iwork], &c__1, &b_ref (1, 1), &c__1);
			}

			/* Zero out below first M rows of B */

			i__1 = *n - *m;
			NUMlapack_dlaset ("F", &i__1, nrhs, &c_b74, &c_b74, &b_ref (*m + 1, 1), ldb);
			iwork = itau + *m;

			/* Multiply transpose(Q) by B (Workspace: need M+NRHS, prefer
			   M+NRHS*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dormlq ("L", "T", n, nrhs, m, &a[a_offset], lda, &work[itau], &b[b_offset], ldb,
			                  &work[iwork], &i__1, info);

		} else {

			/* Path 2 - remaining underdetermined cases */

			ie = 1;
			itauq = ie + *m;
			itaup = itauq + *m;
			iwork = itaup + *m;

			/* Bidiagonalize A (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dgebrd (m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
			                  &work[iwork], &i__1, info);

			/* Multiply B by transpose of left bidiagonalizing vectors
			   (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dormbr ("Q", "L", "T", m, nrhs, n, &a[a_offset], lda, &work[itauq], &b[b_offset], ldb,
			                  &work[iwork], &i__1, info);

			/* Generate right bidiagonalizing vectors in A (Workspace: need
			   4*M, prefer 3*M+M*NB) */

			i__1 = *lwork - iwork + 1;
			NUMlapack_dorgbr ("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[iwork], &i__1, info);
			iwork = ie + *m;

			/* Perform bidiagonal QR iteration, computing right singular
			   vectors of A in A and multiplying B by transpose of left
			   singular vectors (Workspace: need BDSPAC) */

			NUMlapack_dbdsqr ("L", m, n, &c__0, nrhs, &s[1], &work[ie], &a[a_offset], lda, vdum, &c__1,
			                  &b[b_offset], ldb, &work[iwork], info);
			if (*info != 0) {
				goto L70;
			}

			/* Multiply B by reciprocals of singular values

			   Computing MAX */
			d__1 = *rcond * s[1];
			thr = MAX (d__1, sfmin);
			if (*rcond < 0.) {
				/* Computing MAX */
				d__1 = eps * s[1];
				thr = MAX (d__1, sfmin);
			}
			*rank = 0;
			i__1 = *m;
			for (i__ = 1; i__ <= i__1; ++i__) {
				if (s[i__] > thr) {
					NUMlapack_drscl (nrhs, &s[i__], &b_ref (i__, 1), ldb);
					++ (*rank);
				} else {
					NUMlapack_dlaset ("F", &c__1, nrhs, &c_b74, &c_b74, &b_ref (i__, 1), ldb);
				}
				/* L50: */
			}

			/* Multiply B by right singular vectors of A (Workspace: need N,
			   prefer N*NRHS) */

			if (*lwork >= *ldb * *nrhs && *nrhs > 1) {
				NUMblas_dgemm ("T", "N", n, nrhs, m, &c_b108, &a[a_offset], lda, &b[b_offset], ldb, &c_b74, &work[1],
				               ldb);
				NUMlapack_dlacpy ("F", n, nrhs, &work[1], ldb, &b[b_offset], ldb);
			} else if (*nrhs > 1) {
				chunk = *lwork / *n;
				i__1 = *nrhs;
				i__2 = chunk;
				for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
					/* Computing MIN */
					i__3 = *nrhs - i__ + 1;
					bl = MIN (i__3, chunk);
					NUMblas_dgemm ("T", "N", n, &bl, m, &c_b108, &a[a_offset], lda, &b_ref (1, i__), ldb, &c_b74,
					               &work[1], n);
					NUMlapack_dlacpy ("F", n, &bl, &work[1], n, &b_ref (1, i__), ldb);
					/* L60: */
				}
			} else {
				NUMblas_dgemv ("T", m, n, &c_b108, &a[a_offset], lda, &b[b_offset], &c__1, &c_b74, &work[1], &c__1);
				NUMblas_dcopy (n, &work[1], &c__1, &b[b_offset], &c__1);
			}
		}
	}

	/* Undo scaling */

	if (iascl == 1) {
		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, info);
		NUMlapack_dlascl ("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &minmn, info);
	} else if (iascl == 2) {
		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, info);
		NUMlapack_dlascl ("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &minmn, info);
	}
	if (ibscl == 1) {
		NUMlapack_dlascl ("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, info);
	} else if (ibscl == 2) {
		NUMlapack_dlascl ("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, info);
	}

L70:
	work[1] = (double) maxwrk;
	return 0;
}								/* NUMlapack_dgelss */

#undef b_ref

int NUMlapack_dgeqpf (integer *m, integer *n, double *a, integer *lda, integer *jpvt, double *tau, double *work, integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3;
	double d__1, d__2;

	/* Local variables */
	static double temp;
	static double temp2;
	static integer i__, j;
	static integer itemp;
	static integer ma, mn;
	static double aii;
	static integer pvt;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--jpvt;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGEQPF", &i__1);
		return 0;
	}

	mn = MIN (*m, *n);

	/* Move initial columns up front */

	itemp = 1;
	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
		if (jpvt[i__] != 0) {
			if (i__ != itemp) {
				NUMblas_dswap (m, &a_ref (1, i__), &c__1, &a_ref (1, itemp), &c__1);
				jpvt[i__] = jpvt[itemp];
				jpvt[itemp] = i__;
			} else {
				jpvt[i__] = i__;
			}
			++itemp;
		} else {
			jpvt[i__] = i__;
		}
		/* L10: */
	}
	--itemp;

	/* Compute the QR factorization and update remaining columns */

	if (itemp > 0) {
		ma = MIN (itemp, *m);
		NUMlapack_dgeqr2 (m, &ma, &a[a_offset], lda, &tau[1], &work[1], info);
		if (ma < *n) {
			i__1 = *n - ma;
			NUMlapack_dorm2r ("Left", "Transpose", m, &i__1, &ma, &a[a_offset], lda, &tau[1], &a_ref (1,
			                  ma + 1), lda, &work[1], info);
		}
	}

	if (itemp < mn) {

		/* Initialize partial column norms. The first n elements of work
		   store the exact column norms. */

		i__1 = *n;
		for (i__ = itemp + 1; i__ <= i__1; ++i__) {
			i__2 = *m - itemp;
			work[i__] = NUMblas_dnrm2 (&i__2, &a_ref (itemp + 1, i__), &c__1);
			work[*n + i__] = work[i__];
			/* L20: */
		}

		/* Compute factorization */

		i__1 = mn;
		for (i__ = itemp + 1; i__ <= i__1; ++i__) {

			/* Determine ith pivot column and swap if necessary */

			i__2 = *n - i__ + 1;
			pvt = i__ - 1 + NUMblas_idamax (&i__2, &work[i__], &c__1);

			if (pvt != i__) {
				NUMblas_dswap (m, &a_ref (1, pvt), &c__1, &a_ref (1, i__), &c__1);
				itemp = jpvt[pvt];
				jpvt[pvt] = jpvt[i__];
				jpvt[i__] = itemp;
				work[pvt] = work[i__];
				work[*n + pvt] = work[*n + i__];
			}

			/* Generate elementary reflector H(i) */

			if (i__ < *m) {
				i__2 = *m - i__ + 1;
				NUMlapack_dlarfg (&i__2, &a_ref (i__, i__), &a_ref (i__ + 1, i__), &c__1, &tau[i__]);
			} else {
				NUMlapack_dlarfg (&c__1, &a_ref (*m, *m), &a_ref (*m, *m), &c__1, &tau[*m]);
			}

			if (i__ < *n) {

				/* Apply H(i) to A(i:m,i+1:n) from the left */

				aii = a_ref (i__, i__);
				a_ref (i__, i__) = 1.;
				i__2 = *m - i__ + 1;
				i__3 = *n - i__;
				NUMlapack_dlarf ("LEFT", &i__2, &i__3, &a_ref (i__, i__), &c__1, &tau[i__], &a_ref (i__,
				                 i__ + 1), lda, &work[ (*n << 1) + 1]);
				a_ref (i__, i__) = aii;
			}

			/* Update partial column norms */

			i__2 = *n;
			for (j = i__ + 1; j <= i__2; ++j) {
				if (work[j] != 0.) {
					/* Computing 2nd power */
					d__2 = (d__1 = a_ref (i__, j), fabs (d__1)) / work[j];
					temp = 1. - d__2 * d__2;
					temp = MAX (temp, 0.);
					/* Computing 2nd power */
					d__1 = work[j] / work[*n + j];
					temp2 = temp * .05 * (d__1 * d__1) + 1.;
					if (temp2 == 1.) {
						if (*m - i__ > 0) {
							i__3 = *m - i__;
							work[j] = NUMblas_dnrm2 (&i__3, &a_ref (i__ + 1, j), &c__1);
							work[*n + j] = work[j];
						} else {
							work[j] = 0.;
							work[*n + j] = 0.;
						}
					} else {
						work[j] *= sqrt (temp);
					}
				}
				/* L30: */
			}

			/* L40: */
		}
	}
	return 0;
}								/* NUMlapack_dgeqpf */

int NUMlapack_dgeqr2 (integer *m, integer *n, double *a, integer *lda, double *tau, double *work, integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3;

	/* Local variables */
	static integer i__, k;
	static double aii;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGEQR2", &i__1);
		return 0;
	}

	k = MIN (*m, *n);

	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {

		/* Generate elementary reflector H(i) to annihilate A(i+1:m,i)

		   Computing MIN */
		i__2 = i__ + 1;
		i__3 = *m - i__ + 1;
		NUMlapack_dlarfg (&i__3, &a_ref (i__, i__), &a_ref (MIN (i__2, *m), i__), &c__1, &tau[i__]);
		if (i__ < *n) {

			/* Apply H(i) to A(i:m,i+1:n) from the left */

			aii = a_ref (i__, i__);
			a_ref (i__, i__) = 1.;
			i__2 = *m - i__ + 1;
			i__3 = *n - i__;
			NUMlapack_dlarf ("Left", &i__2, &i__3, &a_ref (i__, i__), &c__1, &tau[i__], &a_ref (i__, i__ + 1),
			                 lda, &work[1]);
			a_ref (i__, i__) = aii;
		}
		/* L10: */
	}
	return 0;
}								/* NUMlapack_dgeqr2 */

int NUMlapack_dgeqrf (integer *m, integer *n, double *a, integer *lda, double *tau, double *work, integer *lwork,
                      integer *info) {
	/* Table of constant values */
	static integer c__1 = 1;
	static integer c_n1 = -1;
	static integer c__3 = 3;
	static integer c__2 = 2;

	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2, i__3, i__4;

	/* Local variables */
	static integer i__, k, nbmin, iinfo;
	static integer ib, nb;
	static integer nx;
	static integer ldwork, lwkopt;
	static integer lquery;
	static integer iws;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	nb = NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
	lwkopt = *n * nb;
	work[1] = (double) lwkopt;
	lquery = *lwork == -1;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	} else if (*lwork < MAX (1, *n) && !lquery) {
		*info = -7;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGEQRF", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	k = MIN (*m, *n);
	if (k == 0) {
		work[1] = 1.;
		return 0;
	}

	nbmin = 2;
	nx = 0;
	iws = *n;
	if (nb > 1 && nb < k) {

		/* Determine when to cross over from blocked to unblocked code.

		   Computing MAX */
		i__1 = 0, i__2 = NUMlapack_ilaenv (&c__3, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
		nx = MAX (i__1, i__2);
		if (nx < k) {

			/* Determine if workspace is large enough for blocked code. */

			ldwork = *n;
			iws = ldwork * nb;
			if (*lwork < iws) {

				/* Not enough workspace to use optimal NB: reduce NB and
				   determine the minimum value of NB. */

				nb = *lwork / ldwork;
				/* Computing MAX */
				i__1 = 2, i__2 = NUMlapack_ilaenv (&c__2, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
				nbmin = MAX (i__1, i__2);
			}
		}
	}

	if (nb >= nbmin && nb < k && nx < k) {

		/* Use blocked code initially */

		i__1 = k - nx;
		i__2 = nb;
		for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
			/* Computing MIN */
			i__3 = k - i__ + 1;
			ib = MIN (i__3, nb);

			/* Compute the QR factorization of the current block
			   A(i:m,i:i+ib-1) */

			i__3 = *m - i__ + 1;
			NUMlapack_dgeqr2 (&i__3, &ib, &a_ref (i__, i__), lda, &tau[i__], &work[1], &iinfo);
			if (i__ + ib <= *n) {

				/* Form the triangular factor of the block reflector H = H(i)
				   H(i+1) . . . H(i+ib-1) */

				i__3 = *m - i__ + 1;
				NUMlapack_dlarft ("Forward", "Columnwise", &i__3, &ib, &a_ref (i__, i__), lda, &tau[i__],
				                  &work[1], &ldwork);

				/* Apply H' to A(i:m,i+ib:n) from the left */

				i__3 = *m - i__ + 1;
				i__4 = *n - i__ - ib + 1;
				NUMlapack_dlarfb ("Left", "Transpose", "Forward", "Columnwise", &i__3, &i__4, &ib, &a_ref (i__,
				                  i__), lda, &work[1], &ldwork, &a_ref (i__, i__ + ib), lda, &work[ib + 1], &ldwork);
			}
			/* L10: */
		}
	} else {
		i__ = 1;
	}

	/* Use unblocked code to factor the last or only block. */

	if (i__ <= k) {
		i__2 = *m - i__ + 1;
		i__1 = *n - i__ + 1;
		NUMlapack_dgeqr2 (&i__2, &i__1, &a_ref (i__, i__), lda, &tau[i__], &work[1], &iinfo);
	}

	work[1] = (double) iws;
	return 0;
}								/* NUMlapack_dgeqrf */

int NUMlapack_dgerq2 (integer *m, integer *n, double *a, integer *lda, double *tau, double *work, integer *info) {
	/* System generated locals */
	integer a_dim1, a_offset, i__1, i__2;

	/* Local variables */
	static integer i__, k;
	static double aii;

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--tau;
	--work;

	/* Function Body */
	*info = 0;
	if (*m < 0) {
		*info = -1;
	} else if (*n < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *m)) {
		*info = -4;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGERQ2", &i__1);
		return 0;
	}

	k = MIN (*m, *n);

	for (i__ = k; i__ >= 1; --i__) {

		/* Generate elementary reflector H(i) to annihilate
		   A(m-k+i,1:n-k+i-1) */

		i__1 = *n - k + i__;
		NUMlapack_dlarfg (&i__1, &a_ref (*m - k + i__, *n - k + i__), &a_ref (*m - k + i__, 1), lda, &tau[i__]);

		/* Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right */

		aii = a_ref (*m - k + i__, *n - k + i__);
		a_ref (*m - k + i__, *n - k + i__) = 1.;
		i__1 = *m - k + i__ - 1;
		i__2 = *n - k + i__;
		NUMlapack_dlarf ("Right", &i__1, &i__2, &a_ref (*m - k + i__, 1), lda, &tau[i__], &a[a_offset], lda,
		                 &work[1]);
		a_ref (*m - k + i__, *n - k + i__) = aii;
		/* L10: */
	}
	return 0;
}								/* NUMlapack_dgerq2 */

int NUMlapack_dgesv (integer *n, integer *nrhs, double *a, integer *lda, integer *ipiv, double *b, integer *ldb, integer *info) {
	/* System generated locals */
	integer a_dim1, a_offset, b_dim1, b_offset, i__1;

	/* Local variables */

	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--ipiv;
	b_dim1 = *ldb;
	b_offset = 1 + b_dim1 * 1;
	b -= b_offset;

	/* Function Body */
	*info = 0;
	if (*n < 0) {
		*info = -1;
	} else if (*nrhs < 0) {
		*info = -2;
	} else if (*lda < MAX (1, *n)) {
		*info = -4;
	} else if (*ldb < MAX (1, *n)) {
		*info = -7;
	}
	if (*info != 0) {
		i__1 = - (*info);
		xerbla_ ("DGESV ", &i__1);
		return 0;
	}

	/* Compute the LU factorization of A. */

	NUMlapack_dgetrf (n, n, &a[a_offset], lda, &ipiv[1], info);
	if (*info == 0) {

		/* Solve the system A*X = B, overwriting B with X. */

		NUMlapack_dgetrs ("No transpose", n, nrhs, &a[a_offset], lda, &ipiv[1], &b[b_offset], ldb, info);
	}
	return 0;
}								/* NUMlapack_dgesv */

#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1]

int NUMlapack_dgesvd (const char *jobu, const char *jobvt, integer *m, integer *n, double *a, integer *lda, double *s, double *u,
                      integer *ldu, double *vt, integer *ldvt, double *work, integer *lwork, integer *info) {
	/* System generated locals */
	const char *a__1[2];
	integer a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1[2], i__2, i__3, i__4;
	char ch__1[2];

	/* Local variables */
	static integer iscl;
	static double anrm;
	static integer ierr, itau, ncvt, nrvt, i__;
	static integer chunk, minmn, wrkbl, itaup, itauq, mnthr, iwork;
	static integer wntua, wntva, wntun, wntuo, wntvn, wntvo, wntus, wntvs;
	static integer ie;
	static integer ir, bdspac, iu;
	static double bignum;
	static integer ldwrkr, minwrk, ldwrku, maxwrk;
	static double smlnum;
	static integer lquery, wntuas, wntvas;
	static integer blk, ncu;
	static double dum[1], eps;
	static integer nru;

	/* Parameter adjustments */
	a_dim1 = *lda;
	a_offset = 1 + a_dim1 * 1;
	a -= a_offset;
	--s;
	u_dim1 = *ldu;
	u_offset = 1 + u_dim1 * 1;
	u -= u_offset;
	vt_dim1 = *ldvt;
	vt_offset = 1 + vt_dim1 * 1;
	vt -= vt_offset;
	--work;

	/* Function Body */
	*info = 0;
	minmn = MIN (*m, *n);
	/* Writing concatenation */
	i__1[0] = 1, a__1[0] = jobu;
	i__1[1] = 1, a__1[1] = jobvt;
	s_cat (ch__1, a__1, i__1, &c__2, 2);
	mnthr = NUMlapack_ilaenv (&c__6, "DGESVD", ch__1, m, n, &c__0, &c__0, 6, 2);
	wntua = lsame_ (jobu, "A");
	wntus = lsame_ (jobu, "S");
	wntuas = wntua || wntus;
	wntuo = lsame_ (jobu, "O");
	wntun = lsame_ (jobu, "N");
	wntva = lsame_ (jobvt, "A");
	wntvs = lsame_ (jobvt, "S");
	wntvas = wntva || wntvs;
	wntvo = lsame_ (jobvt, "O");
	wntvn = lsame_ (jobvt, "N");
	minwrk = 1;
	lquery = *lwork == -1;

	if (! (wntua || wntus || wntuo || wntun)) {
		*info = -1;
	} else if (! (wntva || wntvs || wntvo || wntvn) || wntvo && wntuo) {
		*info = -2;
	} else if (*m < 0) {
		*info = -3;
	} else if (*n < 0) {
		*info = -4;
	} else if (*lda < MAX (1, *m)) {
		*info = -6;
	} else if (*ldu < 1 || wntuas && *ldu < *m) {
		*info = -9;
	} else if (*ldvt < 1 || wntva && *ldvt < *n || wntvs && *ldvt < minmn) {
		*info = -11;
	}

	/* Compute workspace (Note: Comments in the code beginning "Workspace:"
	   describe the minimal amount of workspace needed at that point in the
	   code, as well as the preferred amount for good performance. NB refers
	   to the optimal block size for the immediately following subroutine, as
	   returned by ILAENV.) */

	if (*info == 0 && (*lwork >= 1 || lquery) && *m > 0 && *n > 0) {
		if (*m >= *n) {

			/* Compute space needed for DBDSQR */

			bdspac = *n * 5;
			if (*m >= mnthr) {
				if (wntun) {

					/* Path 1 (M much larger than N, JOBU='N') */

					maxwrk = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
					if (wntvo || wntvas) {
						/* Computing MAX */
						i__2 = maxwrk, i__3 =
						           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
						maxwrk = MAX (i__2, i__3);
					}
					maxwrk = MAX (maxwrk, bdspac);
					/* Computing MAX */
					i__2 = *n << 2;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntuo && wntvn) {

					/* Path 2 (M much larger than N, JOBU='O', JOBVT='N') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *n * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					/* Computing MAX */
					i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
					maxwrk = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntuo && wntvas) {

					/* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
					   'A') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *n * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					/* Computing MAX */
					i__2 = *n * *n + wrkbl, i__3 = *n * *n + *m * *n + *n;
					maxwrk = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntus && wntvn) {

					/* Path 4 (M much larger than N, JOBU='S', JOBVT='N') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *n * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *n * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntus && wntvo) {

					/* Path 5 (M much larger than N, JOBU='S', JOBVT='O') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *n * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = (*n << 1) * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntus && wntvas) {

					/* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
					   'A') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *n * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *n * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntua && wntvn) {

					/* Path 7 (M much larger than N, JOBU='A', JOBVT='N') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *m * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, m, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *n * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntua && wntvo) {

					/* Path 8 (M much larger than N, JOBU='A', JOBVT='O') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *m * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, m, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = (*n << 1) * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntua && wntvas) {

					/* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
					   'A') */

					wrkbl = *n + *n * NUMlapack_ilaenv (&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n + *m * NUMlapack_ilaenv (&c__1, "DORGQR", " ", m, m, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", n, n, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *n * *n + wrkbl;
					/* Computing MAX */
					i__2 = *n * 3 + *m;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				}
			} else {

				/* Path 10 (M at least N, but not much larger) */

				maxwrk = *n * 3 + (*m + *n) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
				if (wntus || wntuo) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *n * 3 + *n * NUMlapack_ilaenv (&c__1, "DORG" "BR", "Q", m, n, n, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				if (wntua) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *n * 3 + *m * NUMlapack_ilaenv (&c__1, "DORG" "BR", "Q", m, m, n, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				if (!wntvn) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *n * 3 + (*n - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", n, n, n, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				maxwrk = MAX (maxwrk, bdspac);
				/* Computing MAX */
				i__2 = *n * 3 + *m;
				minwrk = MAX (i__2, bdspac);
				maxwrk = MAX (maxwrk, minwrk);
			}
		} else {

			/* Compute space needed for DBDSQR */

			bdspac = *m * 5;
			if (*n >= mnthr) {
				if (wntvn) {

					/* Path 1t(N much larger than M, JOBVT='N') */

					maxwrk = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
					if (wntuo || wntuas) {
						/* Computing MAX */
						i__2 = maxwrk, i__3 =
						           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
						maxwrk = MAX (i__2, i__3);
					}
					maxwrk = MAX (maxwrk, bdspac);
					/* Computing MAX */
					i__2 = *m << 2;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntvo && wntun) {

					/* Path 2t(N much larger than M, JOBU='N', JOBVT='O') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *m * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", m, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					/* Computing MAX */
					i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
					maxwrk = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntvo && wntuas) {

					/* Path 3t(N much larger than M, JOBU='S' or 'A',
					   JOBVT='O') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *m * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", m, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					/* Computing MAX */
					i__2 = *m * *m + wrkbl, i__3 = *m * *m + *m * *n + *m;
					maxwrk = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntvs && wntun) {

					/* Path 4t(N much larger than M, JOBU='N', JOBVT='S') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *m * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", m, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *m * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntvs && wntuo) {

					/* Path 5t(N much larger than M, JOBU='O', JOBVT='S') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *m * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", m, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = (*m << 1) * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntvs && wntuas) {

					/* Path 6t(N much larger than M, JOBU='S' or 'A',
					   JOBVT='S') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *m * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", m, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *m * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntva && wntun) {

					/* Path 7t(N much larger than M, JOBU='N', JOBVT='A') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *n * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", n, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *m * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntva && wntuo) {

					/* Path 8t(N much larger than M, JOBU='O', JOBVT='A') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *n * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", n, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = (*m << 1) * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				} else if (wntva && wntuas) {

					/* Path 9t(N much larger than M, JOBU='S' or 'A',
					   JOBVT='A') */

					wrkbl = *m + *m * NUMlapack_ilaenv (&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1, 6, 1);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m + *n * NUMlapack_ilaenv (&c__1, "DORGLQ", " ", n, n, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m << 1) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "P", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					/* Computing MAX */
					i__2 = wrkbl, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					wrkbl = MAX (i__2, i__3);
					wrkbl = MAX (wrkbl, bdspac);
					maxwrk = *m * *m + wrkbl;
					/* Computing MAX */
					i__2 = *m * 3 + *n;
					minwrk = MAX (i__2, bdspac);
					maxwrk = MAX (maxwrk, minwrk);
				}
			} else {

				/* Path 10t(N greater than M, but not much larger) */

				maxwrk = *m * 3 + (*m + *n) * NUMlapack_ilaenv (&c__1, "DGEBRD", " ", m, n, &c_n1, &c_n1, 6, 1);
				if (wntvs || wntvo) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *m * 3 + *m * NUMlapack_ilaenv (&c__1, "DORG" "BR", "P", m, n, m, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				if (wntva) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *m * 3 + *n * NUMlapack_ilaenv (&c__1, "DORG" "BR", "P", n, n, m, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				if (!wntun) {
					/* Computing MAX */
					i__2 = maxwrk, i__3 =
					           *m * 3 + (*m - 1) * NUMlapack_ilaenv (&c__1, "DORGBR", "Q", m, m, m, &c_n1, 6, 1);
					maxwrk = MAX (i__2, i__3);
				}
				maxwrk = MAX (maxwrk, bdspac);
				/* Computing MAX */
				i__2 = *m * 3 + *n;
				minwrk = MAX (i__2, bdspac);
				maxwrk = MAX (maxwrk, minwrk);
			}
		}
		work[1] = (double) maxwrk;
	}

	if (*lwork < minwrk && !lquery) {
		*info = -13;
	}
	if (*info != 0) {
		i__2 = - (*info);
		xerbla_ ("DGESVD", &i__2);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/* Quick return if possible */

	if (*m == 0 || *n == 0) {
		if (*lwork >= 1) {
			work[1] = 1.;
		}
		return 0;
	}

	/* Get machine constants */

	eps = NUMblas_dlamch ("P");
	smlnum = sqrt (NUMblas_dlamch ("S")) / eps;
	bignum = 1. / smlnum;

	/* Scale A if max element outside range [SMLNUM,BIGNUM] */

	anrm = NUMlapack_dlange ("M", m, n, &a[a_offset], lda, dum);
	iscl = 0;
	if (anrm > 0. && anrm < smlnum) {
		iscl = 1;
		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, &ierr);
	} else if (anrm > bignum) {
		iscl = 1;
		NUMlapack_dlascl ("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, &ierr);
	}

	if (*m >= *n) {

		/* A has at least as many rows as columns. If A has sufficiently more
		   rows than columns, first reduce using the QR decomposition (if
		   sufficient workspace available) */

		if (*m >= mnthr) {

			if (wntun) {

				/* Path 1 (M much larger than N, JOBU='N') No left singular
				   vectors to be computed */

				itau = 1;
				iwork = itau + *n;

				/* Compute A=Q*R (Workspace: need 2*N, prefer N+N*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

				/* Zero out below R */

				i__2 = *n - 1;
				i__3 = *n - 1;
				NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &a_ref (2, 1), lda);
				ie = 1;
				itauq = ie + *n;
				itaup = itauq + *n;
				iwork = itaup + *n;

				/* Bidiagonalize R in A (Workspace: need 4*N, prefer
				   3*N+2*N*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dgebrd (n, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
				                  &work[iwork], &i__2, &ierr);
				ncvt = 0;
				if (wntvo || wntvas) {

					/* If right singular vectors desired, generate P'.
					   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2,
					                  &ierr);
					ncvt = *n;
				}
				iwork = ie + *n;

				/* Perform bidiagonal QR iteration, computing right singular
				   vectors of A in A if desired (Workspace: need BDSPAC) */

				NUMlapack_dbdsqr ("U", n, &ncvt, &c__0, &c__0, &s[1], &work[ie], &a[a_offset], lda, dum, &c__1,
				                  dum, &c__1, &work[iwork], info);

				/* If right singular vectors desired in VT, copy them there */

				if (wntvas) {
					NUMlapack_dlacpy ("F", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
				}

			} else if (wntuo && wntvn) {

				/* Path 2 (M much larger than N, JOBU='O', JOBVT='N') N left
				   singular vectors to be overwritten on A and no right
				   singular vectors to be computed

				   Computing MAX */
				i__2 = *n << 2;
				if (*lwork >= *n * *n + MAX (i__2, bdspac)) {

					/* Sufficient workspace for a fast algorithm */

					ir = 1;
					/* Computing MAX */
					i__2 = wrkbl, i__3 = *lda * *n + *n;
					if (*lwork >= MAX (i__2, i__3) + *lda * *n) {

						/* WORK(IU) is LDA by N, WORK(IR) is LDA by N */

						ldwrku = *lda;
						ldwrkr = *lda;
					} else {	/* if(complicated condition) */
						/* Computing MAX */
						i__2 = wrkbl, i__3 = *lda * *n + *n;
						if (*lwork >= MAX (i__2, i__3) + *n * *n) {

							/* WORK(IU) is LDA by N, WORK(IR) is N by N */

							ldwrku = *lda;
							ldwrkr = *n;
						} else {

							/* WORK(IU) is LDWRKU by N, WORK(IR) is N by N */

							ldwrku = (*lwork - *n * *n - *n) / *n;
							ldwrkr = *n;
						}
					}
					itau = ir + ldwrkr * *n;
					iwork = itau + *n;

					/* Compute A=Q*R (Workspace: need N*N+2*N, prefer
					   N*N+N+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

					/* Copy R to WORK(IR) and zero out below it */

					NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
					i__2 = *n - 1;
					i__3 = *n - 1;
					NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[ir + 1], &ldwrkr);

					/* Generate Q in A (Workspace: need N*N+2*N, prefer
					   N*N+N+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
					ie = itau;
					itauq = ie + *n;
					itaup = itauq + *n;
					iwork = itaup + *n;

					/* Bidiagonalize R in WORK(IR) (Workspace: need N*N+4*N,
					   prefer N*N+3*N+2*N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgebrd (n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__2, &ierr);

					/* Generate left vectors bidiagonalizing R (Workspace:
					   need N*N+4*N, prefer N*N+3*N+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__2,
					                  &ierr);
					iwork = ie + *n;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of R in WORK(IR) (Workspace: need
					   N*N+BDSPAC) */

					NUMlapack_dbdsqr ("U", n, &c__0, n, &c__0, &s[1], &work[ie], dum, &c__1, &work[ir], &ldwrkr,
					                  dum, &c__1, &work[iwork], info);
					iu = ie + *n;

					/* Multiply Q in A by left singular vectors of R in
					   WORK(IR), storing result in WORK(IU) and copying to A
					   (Workspace: need N*N+2*N, prefer N*N+M*N+N) */

					i__2 = *m;
					i__3 = ldwrku;
					for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__3) {
						/* Computing MIN */
						i__4 = *m - i__ + 1;
						chunk = MIN (i__4, ldwrku);
						NUMblas_dgemm ("N", "N", &chunk, n, n, &c_b438, &a_ref (i__, 1), lda, &work[ir], &ldwrkr,
						               &c_b416, &work[iu], &ldwrku);
						NUMlapack_dlacpy ("F", &chunk, n, &work[iu], &ldwrku, &a_ref (i__, 1), lda);
						/* L10: */
					}

				} else {

					/* Insufficient workspace for a fast algorithm */

					ie = 1;
					itauq = ie + *n;
					itaup = itauq + *n;
					iwork = itaup + *n;

					/* Bidiagonalize A (Workspace: need 3*N+M, prefer
					   3*N+(M+N)*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgebrd (m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__3, &ierr);

					/* Generate left vectors bidiagonalizing A (Workspace:
					   need 4*N, prefer 3*N+N*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[iwork], &i__3,
					                  &ierr);
					iwork = ie + *n;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of A in A (Workspace: need BDSPAC) */

					NUMlapack_dbdsqr ("U", n, &c__0, m, &c__0, &s[1], &work[ie], dum, &c__1, &a[a_offset], lda,
					                  dum, &c__1, &work[iwork], info);

				}

			} else if (wntuo && wntvas) {

				/* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
				   N left singular vectors to be overwritten on A and N right
				   singular vectors to be computed in VT

				   Computing MAX */
				i__3 = *n << 2;
				if (*lwork >= *n * *n + MAX (i__3, bdspac)) {

					/* Sufficient workspace for a fast algorithm */

					ir = 1;
					/* Computing MAX */
					i__3 = wrkbl, i__2 = *lda * *n + *n;
					if (*lwork >= MAX (i__3, i__2) + *lda * *n) {

						/* WORK(IU) is LDA by N and WORK(IR) is LDA by N */

						ldwrku = *lda;
						ldwrkr = *lda;
					} else {	/* if(complicated condition) */
						/* Computing MAX */
						i__3 = wrkbl, i__2 = *lda * *n + *n;
						if (*lwork >= MAX (i__3, i__2) + *n * *n) {

							/* WORK(IU) is LDA by N and WORK(IR) is N by N */

							ldwrku = *lda;
							ldwrkr = *n;
						} else {

							/* WORK(IU) is LDWRKU by N and WORK(IR) is N by N
							 */

							ldwrku = (*lwork - *n * *n - *n) / *n;
							ldwrkr = *n;
						}
					}
					itau = ir + ldwrkr * *n;
					iwork = itau + *n;

					/* Compute A=Q*R (Workspace: need N*N+2*N, prefer
					   N*N+N+N*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__3, &ierr);

					/* Copy R to VT, zeroing out below it */

					NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
					i__3 = *n - 1;
					i__2 = *n - 1;
					NUMlapack_dlaset ("L", &i__3, &i__2, &c_b416, &c_b416, &vt_ref (2, 1), ldvt);

					/* Generate Q in A (Workspace: need N*N+2*N, prefer
					   N*N+N+N*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__3, &ierr);
					ie = itau;
					itauq = ie + *n;
					itaup = itauq + *n;
					iwork = itaup + *n;

					/* Bidiagonalize R in VT, copying result to WORK(IR)
					   (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgebrd (n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__3, &ierr);
					NUMlapack_dlacpy ("L", n, n, &vt[vt_offset], ldvt, &work[ir], &ldwrkr);

					/* Generate left vectors bidiagonalizing R in WORK(IR)
					   (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__3,
					                  &ierr);

					/* Generate right vectors bidiagonalizing R in VT
					   (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__3,
					                  &ierr);
					iwork = ie + *n;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of R in WORK(IR) and computing right
					   singular vectors of R in VT (Workspace: need
					   N*N+BDSPAC) */

					NUMlapack_dbdsqr ("U", n, n, n, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt, &work[ir],
					                  &ldwrkr, dum, &c__1, &work[iwork], info);
					iu = ie + *n;

					/* Multiply Q in A by left singular vectors of R in
					   WORK(IR), storing result in WORK(IU) and copying to A
					   (Workspace: need N*N+2*N, prefer N*N+M*N+N) */

					i__3 = *m;
					i__2 = ldwrku;
					for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ += i__2) {
						/* Computing MIN */
						i__4 = *m - i__ + 1;
						chunk = MIN (i__4, ldwrku);
						NUMblas_dgemm ("N", "N", &chunk, n, n, &c_b438, &a_ref (i__, 1), lda, &work[ir], &ldwrkr,
						               &c_b416, &work[iu], &ldwrku);
						NUMlapack_dlacpy ("F", &chunk, n, &work[iu], &ldwrku, &a_ref (i__, 1), lda);
						/* L20: */
					}

				} else {

					/* Insufficient workspace for a fast algorithm */

					itau = 1;
					iwork = itau + *n;

					/* Compute A=Q*R (Workspace: need 2*N, prefer N+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

					/* Copy R to VT, zeroing out below it */

					NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
					i__2 = *n - 1;
					i__3 = *n - 1;
					NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &vt_ref (2, 1), ldvt);

					/* Generate Q in A (Workspace: need 2*N, prefer N+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
					ie = itau;
					itauq = ie + *n;
					itaup = itauq + *n;
					iwork = itaup + *n;

					/* Bidiagonalize R in VT (Workspace: need 4*N, prefer
					   3*N+2*N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgebrd (n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__2, &ierr);

					/* Multiply Q in A by left vectors bidiagonalizing R
					   (Workspace: need 3*N+M, prefer 3*N+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dormbr ("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &work[itauq], &a[a_offset],
					                  lda, &work[iwork], &i__2, &ierr);

					/* Generate right vectors bidiagonalizing R in VT
					   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2,
					                  &ierr);
					iwork = ie + *n;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of A in A and computing right
					   singular vectors of A in VT (Workspace: need BDSPAC) */

					NUMlapack_dbdsqr ("U", n, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt, &a[a_offset],
					                  lda, dum, &c__1, &work[iwork], info);

				}

			} else if (wntus) {

				if (wntvn) {

					/* Path 4 (M much larger than N, JOBU='S', JOBVT='N') N
					   left singular vectors to be computed in U and no right
					   singular vectors to be computed

					   Computing MAX */
					i__2 = *n << 2;
					if (*lwork >= *n * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						ir = 1;
						if (*lwork >= wrkbl + *lda * *n) {

							/* WORK(IR) is LDA by N */

							ldwrkr = *lda;
						} else {

							/* WORK(IR) is N by N */

							ldwrkr = *n;
						}
						itau = ir + ldwrkr * *n;
						iwork = itau + *n;

						/* Compute A=Q*R (Workspace: need N*N+2*N, prefer
						   N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to WORK(IR), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[ir + 1], &ldwrkr);

						/* Generate Q in A (Workspace: need N*N+2*N, prefer
						   N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IR) (Workspace: need
						   N*N+4*N, prefer N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Generate left vectors bidiagonalizing R in
						   WORK(IR) (Workspace: need N*N+4*N, prefer
						   N*N+3*N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IR) (Workspace: need
						   N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, &c__0, n, &c__0, &s[1], &work[ie], dum, &c__1, &work[ir],
						                  &ldwrkr, dum, &c__1, &work[iwork], info);

						/* Multiply Q in A by left singular vectors of R in
						   WORK(IR), storing result in U (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &a[a_offset], lda, &work[ir], &ldwrkr, &c_b416,
						               &u[u_offset], ldu);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need 2*N, prefer
						   N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Zero out below R in A */

						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &a_ref (2, 1), lda);

						/* Bidiagonalize R in A (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left vectors bidiagonalizing R
						   (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &a[a_offset], lda, &work[itauq], &u[u_offset],
						                  ldu, &work[iwork], &i__2, &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, &c__0, m, &c__0, &s[1], &work[ie], dum, &c__1, &u[u_offset],
						                  ldu, dum, &c__1, &work[iwork], info);

					}

				} else if (wntvo) {

					/* Path 5 (M much larger than N, JOBU='S', JOBVT='O') N
					   left singular vectors to be computed in U and N right
					   singular vectors to be overwritten on A

					   Computing MAX */
					i__2 = *n << 2;
					if (*lwork >= (*n << 1) * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + (*lda << 1) * *n) {

							/* WORK(IU) is LDA by N and WORK(IR) is LDA by N */

							ldwrku = *lda;
							ir = iu + ldwrku * *n;
							ldwrkr = *lda;
						} else if (*lwork >= wrkbl + (*lda + *n) * *n) {

							/* WORK(IU) is LDA by N and WORK(IR) is N by N */

							ldwrku = *lda;
							ir = iu + ldwrku * *n;
							ldwrkr = *n;
						} else {

							/* WORK(IU) is N by N and WORK(IR) is N by N */

							ldwrku = *n;
							ir = iu + ldwrku * *n;
							ldwrkr = *n;
						}
						itau = ir + ldwrkr * *n;
						iwork = itau + *n;

						/* Compute A=Q*R (Workspace: need 2*N*N+2*N, prefer
						   2*N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to WORK(IU), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[iu + 1], &ldwrku);

						/* Generate Q in A (Workspace: need 2*N*N+2*N, prefer
						   2*N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IU), copying result to
						   WORK(IR) (Workspace: need 2*N*N+4*N, prefer
						   2*N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", n, n, &work[iu], &ldwrku, &work[ir], &ldwrkr);

						/* Generate left bidiagonalizing vectors in WORK(IU)
						   (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[iu], &ldwrku, &work[itauq], &work[iwork], &i__2,
						                  &ierr);

						/* Generate right bidiagonalizing vectors in WORK(IR)
						   (Workspace: need 2*N*N+4*N-1, prefer
						   2*N*N+3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IU) and computing
						   right singular vectors of R in WORK(IR)
						   (Workspace: need 2*N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, n, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, &work[iu],
						                  &ldwrku, dum, &c__1, &work[iwork], info);

						/* Multiply Q in A by left singular vectors of R in
						   WORK(IU), storing result in U (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &a[a_offset], lda, &work[iu], &ldwrku, &c_b416,
						               &u[u_offset], ldu);

						/* Copy right singular vectors of R to A (Workspace:
						   need N*N) */

						NUMlapack_dlacpy ("F", n, n, &work[ir], &ldwrkr, &a[a_offset], lda);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need 2*N, prefer
						   N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Zero out below R in A */

						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &a_ref (2, 1), lda);

						/* Bidiagonalize R in A (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left vectors bidiagonalizing R
						   (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &a[a_offset], lda, &work[itauq], &u[u_offset],
						                  ldu, &work[iwork], &i__2, &ierr);

						/* Generate right vectors bidiagonalizing R in A
						   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in A (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, m, &c__0, &s[1], &work[ie], &a[a_offset], lda,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				} else if (wntvas) {

					/* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
					   'A') N left singular vectors to be computed in U and N
					   right singular vectors to be computed in VT

					   Computing MAX */
					i__2 = *n << 2;
					if (*lwork >= *n * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + *lda * *n) {

							/* WORK(IU) is LDA by N */

							ldwrku = *lda;
						} else {

							/* WORK(IU) is N by N */

							ldwrku = *n;
						}
						itau = iu + ldwrku * *n;
						iwork = itau + *n;

						/* Compute A=Q*R (Workspace: need N*N+2*N, prefer
						   N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to WORK(IU), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[iu + 1], &ldwrku);

						/* Generate Q in A (Workspace: need N*N+2*N, prefer
						   N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IU), copying result to VT
						   (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);

						/* Generate left bidiagonalizing vectors in WORK(IU)
						   (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[iu], &ldwrku, &work[itauq], &work[iwork], &i__2,
						                  &ierr);

						/* Generate right bidiagonalizing vectors in VT
						   (Workspace: need N*N+4*N-1, prefer
						   N*N+3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IU) and computing
						   right singular vectors of R in VT (Workspace: need
						   N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, n, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &work[iu], &ldwrku, dum, &c__1, &work[iwork], info);

						/* Multiply Q in A by left singular vectors of R in
						   WORK(IU), storing result in U (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &a[a_offset], lda, &work[iu], &ldwrku, &c_b416,
						               &u[u_offset], ldu);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need 2*N, prefer
						   N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, n, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to VT, zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &vt_ref (2, 1), ldvt);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in VT (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left bidiagonalizing vectors in
						   VT (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &work[itauq],
						                  &u[u_offset], ldu, &work[iwork], &i__2, &ierr);

						/* Generate right bidiagonalizing vectors in VT
						   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				}

			} else if (wntua) {

				if (wntvn) {

					/* Path 7 (M much larger than N, JOBU='A', JOBVT='N') M
					   left singular vectors to be computed in U and no right
					   singular vectors to be computed

					   Computing MAX */
					i__2 = *n + *m, i__3 = *n << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= *n * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						ir = 1;
						if (*lwork >= wrkbl + *lda * *n) {

							/* WORK(IR) is LDA by N */

							ldwrkr = *lda;
						} else {

							/* WORK(IR) is N by N */

							ldwrkr = *n;
						}
						itau = ir + ldwrkr * *n;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need N*N+2*N, prefer N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Copy R to WORK(IR), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[ir], &ldwrkr);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[ir + 1], &ldwrkr);

						/* Generate Q in U (Workspace: need N*N+N+M, prefer
						   N*N+N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IR) (Workspace: need
						   N*N+4*N, prefer N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Generate left bidiagonalizing vectors in WORK(IR)
						   (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IR) (Workspace: need
						   N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, &c__0, n, &c__0, &s[1], &work[ie], dum, &c__1, &work[ir],
						                  &ldwrkr, dum, &c__1, &work[iwork], info);

						/* Multiply Q in U by left singular vectors of R in
						   WORK(IR), storing result in A (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &u[u_offset], ldu, &work[ir], &ldwrkr, &c_b416,
						               &a[a_offset], lda);

						/* Copy left singular vectors of A from A to U */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need N+M, prefer
						   N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Zero out below R in A */

						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &a_ref (2, 1), lda);

						/* Bidiagonalize R in A (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left bidiagonalizing vectors in
						   A (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &a[a_offset], lda, &work[itauq], &u[u_offset],
						                  ldu, &work[iwork], &i__2, &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, &c__0, m, &c__0, &s[1], &work[ie], dum, &c__1, &u[u_offset],
						                  ldu, dum, &c__1, &work[iwork], info);

					}

				} else if (wntvo) {

					/* Path 8 (M much larger than N, JOBU='A', JOBVT='O') M
					   left singular vectors to be computed in U and N right
					   singular vectors to be overwritten on A

					   Computing MAX */
					i__2 = *n + *m, i__3 = *n << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= (*n << 1) * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + (*lda << 1) * *n) {

							/* WORK(IU) is LDA by N and WORK(IR) is LDA by N */

							ldwrku = *lda;
							ir = iu + ldwrku * *n;
							ldwrkr = *lda;
						} else if (*lwork >= wrkbl + (*lda + *n) * *n) {

							/* WORK(IU) is LDA by N and WORK(IR) is N by N */

							ldwrku = *lda;
							ir = iu + ldwrku * *n;
							ldwrkr = *n;
						} else {

							/* WORK(IU) is N by N and WORK(IR) is N by N */

							ldwrku = *n;
							ir = iu + ldwrku * *n;
							ldwrkr = *n;
						}
						itau = ir + ldwrkr * *n;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N*N+2*N, prefer 2*N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need 2*N*N+N+M, prefer
						   2*N*N+N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to WORK(IU), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[iu + 1], &ldwrku);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IU), copying result to
						   WORK(IR) (Workspace: need 2*N*N+4*N, prefer
						   2*N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", n, n, &work[iu], &ldwrku, &work[ir], &ldwrkr);

						/* Generate left bidiagonalizing vectors in WORK(IU)
						   (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[iu], &ldwrku, &work[itauq], &work[iwork], &i__2,
						                  &ierr);

						/* Generate right bidiagonalizing vectors in WORK(IR)
						   (Workspace: need 2*N*N+4*N-1, prefer
						   2*N*N+3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IU) and computing
						   right singular vectors of R in WORK(IR)
						   (Workspace: need 2*N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, n, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, &work[iu],
						                  &ldwrku, dum, &c__1, &work[iwork], info);

						/* Multiply Q in U by left singular vectors of R in
						   WORK(IU), storing result in A (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &u[u_offset], ldu, &work[iu], &ldwrku, &c_b416,
						               &a[a_offset], lda);

						/* Copy left singular vectors of A from A to U */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Copy right singular vectors of R from WORK(IR) to
						   A */

						NUMlapack_dlacpy ("F", n, n, &work[ir], &ldwrkr, &a[a_offset], lda);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need N+M, prefer
						   N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Zero out below R in A */

						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &a_ref (2, 1), lda);

						/* Bidiagonalize R in A (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left bidiagonalizing vectors in
						   A (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &a[a_offset], lda, &work[itauq], &u[u_offset],
						                  ldu, &work[iwork], &i__2, &ierr);

						/* Generate right bidiagonalizing vectors in A
						   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in A (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, m, &c__0, &s[1], &work[ie], &a[a_offset], lda,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				} else if (wntvas) {

					/* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
					   'A') M left singular vectors to be computed in U and N
					   right singular vectors to be computed in VT

					   Computing MAX */
					i__2 = *n + *m, i__3 = *n << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= *n * *n + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + *lda * *n) {

							/* WORK(IU) is LDA by N */

							ldwrku = *lda;
						} else {

							/* WORK(IU) is N by N */

							ldwrku = *n;
						}
						itau = iu + ldwrku * *n;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need N*N+2*N, prefer N*N+N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need N*N+N+M, prefer
						   N*N+N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R to WORK(IU), zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &work[iu + 1], &ldwrku);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in WORK(IU), copying result to VT
						   (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", n, n, &work[iu], &ldwrku, &vt[vt_offset], ldvt);

						/* Generate left bidiagonalizing vectors in WORK(IU)
						   (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", n, n, n, &work[iu], &ldwrku, &work[itauq], &work[iwork], &i__2,
						                  &ierr);

						/* Generate right bidiagonalizing vectors in VT
						   (Workspace: need N*N+4*N-1, prefer
						   N*N+3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of R in WORK(IU) and computing
						   right singular vectors of R in VT (Workspace: need
						   N*N+BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, n, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &work[iu], &ldwrku, dum, &c__1, &work[iwork], info);

						/* Multiply Q in U by left singular vectors of R in
						   WORK(IU), storing result in A (Workspace: need
						   N*N) */

						NUMblas_dgemm ("N", "N", m, n, n, &c_b438, &u[u_offset], ldu, &work[iu], &ldwrku, &c_b416,
						               &a[a_offset], lda);

						/* Copy left singular vectors of A from A to U */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *n;

						/* Compute A=Q*R, copying result to U (Workspace:
						   need 2*N, prefer N+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgeqrf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);

						/* Generate Q in U (Workspace: need N+M, prefer
						   N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgqr (m, m, n, &u[u_offset], ldu, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy R from A to VT, zeroing out below it */

						NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
						i__2 = *n - 1;
						i__3 = *n - 1;
						NUMlapack_dlaset ("L", &i__2, &i__3, &c_b416, &c_b416, &vt_ref (2, 1), ldvt);
						ie = itau;
						itauq = ie + *n;
						itaup = itauq + *n;
						iwork = itaup + *n;

						/* Bidiagonalize R in VT (Workspace: need 4*N, prefer
						   3*N+2*N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (n, n, &vt[vt_offset], ldvt, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Multiply Q in U by left bidiagonalizing vectors in
						   VT (Workspace: need 3*N+M, prefer 3*N+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("Q", "R", "N", m, n, n, &vt[vt_offset], ldvt, &work[itauq],
						                  &u[u_offset], ldu, &work[iwork], &i__2, &ierr);

						/* Generate right bidiagonalizing vectors in VT
						   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *n;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", n, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				}

			}

		} else {

			/* M .LT. MNTHR

			   Path 10 (M at least N, but not much larger) Reduce to
			   bidiagonal form without QR decomposition */

			ie = 1;
			itauq = ie + *n;
			itaup = itauq + *n;
			iwork = itaup + *n;

			/* Bidiagonalize A (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dgebrd (m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
			                  &work[iwork], &i__2, &ierr);
			if (wntuas) {

				/* If left singular vectors desired in U, copy result to U
				   and generate left bidiagonalizing vectors in U (Workspace:
				   need 3*N+NCU, prefer 3*N+NCU*NB) */

				NUMlapack_dlacpy ("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
				if (wntus) {
					ncu = *n;
				}
				if (wntua) {
					ncu = *m;
				}
				i__2 = *lwork - iwork + 1;
				NUMlapack_dorgbr ("Q", m, &ncu, n, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
			}
			if (wntvas) {

				/* If right singular vectors desired in VT, copy result to VT
				   and generate right bidiagonalizing vectors in VT
				   (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) */

				NUMlapack_dlacpy ("U", n, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
				i__2 = *lwork - iwork + 1;
				NUMlapack_dorgbr ("P", n, n, n, &vt[vt_offset], ldvt, &work[itaup], &work[iwork], &i__2, &ierr);
			}
			if (wntuo) {

				/* If left singular vectors desired in A, generate left
				   bidiagonalizing vectors in A (Workspace: need 4*N, prefer
				   3*N+N*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dorgbr ("Q", m, n, n, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2, &ierr);
			}
			if (wntvo) {

				/* If right singular vectors desired in A, generate right
				   bidiagonalizing vectors in A (Workspace: need 4*N-1,
				   prefer 3*N+(N-1)*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dorgbr ("P", n, n, n, &a[a_offset], lda, &work[itaup], &work[iwork], &i__2, &ierr);
			}
			iwork = ie + *n;
			if (wntuas || wntuo) {
				nru = *m;
			}
			if (wntun) {
				nru = 0;
			}
			if (wntvas || wntvo) {
				ncvt = *n;
			}
			if (wntvn) {
				ncvt = 0;
			}
			if (!wntuo && !wntvo) {

				/* Perform bidiagonal QR iteration, if desired, computing
				   left singular vectors in U and computing right singular
				   vectors in VT (Workspace: need BDSPAC) */

				NUMlapack_dbdsqr ("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
				                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);
			} else if (!wntuo && wntvo) {

				/* Perform bidiagonal QR iteration, if desired, computing
				   left singular vectors in U and computing right singular
				   vectors in A (Workspace: need BDSPAC) */

				NUMlapack_dbdsqr ("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &a[a_offset], lda, &u[u_offset],
				                  ldu, dum, &c__1, &work[iwork], info);
			} else {

				/* Perform bidiagonal QR iteration, if desired, computing
				   left singular vectors in A and computing right singular
				   vectors in VT (Workspace: need BDSPAC) */

				NUMlapack_dbdsqr ("U", n, &ncvt, &nru, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
				                  &a[a_offset], lda, dum, &c__1, &work[iwork], info);
			}

		}

	} else {

		/* A has more columns than rows. If A has sufficiently more columns
		   than rows, first reduce using the LQ decomposition (if sufficient
		   workspace available) */

		if (*n >= mnthr) {

			if (wntvn) {

				/* Path 1t(N much larger than M, JOBVT='N') No right singular
				   vectors to be computed */

				itau = 1;
				iwork = itau + *m;

				/* Compute A=L*Q (Workspace: need 2*M, prefer M+M*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

				/* Zero out above L */

				i__2 = *m - 1;
				i__3 = *m - 1;
				NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &a_ref (1, 2), lda);
				ie = 1;
				itauq = ie + *m;
				itaup = itauq + *m;
				iwork = itaup + *m;

				/* Bidiagonalize L in A (Workspace: need 4*M, prefer
				   3*M+2*M*NB) */

				i__2 = *lwork - iwork + 1;
				NUMlapack_dgebrd (m, m, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
				                  &work[iwork], &i__2, &ierr);
				if (wntuo || wntuas) {

					/* If left singular vectors desired, generate Q
					   (Workspace: need 4*M, prefer 3*M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", m, m, m, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2,
					                  &ierr);
				}
				iwork = ie + *m;
				nru = 0;
				if (wntuo || wntuas) {
					nru = *m;
				}

				/* Perform bidiagonal QR iteration, computing left singular
				   vectors of A in A if desired (Workspace: need BDSPAC) */

				NUMlapack_dbdsqr ("U", m, &c__0, &nru, &c__0, &s[1], &work[ie], dum, &c__1, &a[a_offset], lda,
				                  dum, &c__1, &work[iwork], info);

				/* If left singular vectors desired in U, copy them there */

				if (wntuas) {
					NUMlapack_dlacpy ("F", m, m, &a[a_offset], lda, &u[u_offset], ldu);
				}

			} else if (wntvo && wntun) {

				/* Path 2t(N much larger than M, JOBU='N', JOBVT='O') M right
				   singular vectors to be overwritten on A and no left
				   singular vectors to be computed

				   Computing MAX */
				i__2 = *m << 2;
				if (*lwork >= *m * *m + MAX (i__2, bdspac)) {

					/* Sufficient workspace for a fast algorithm */

					ir = 1;
					/* Computing MAX */
					i__2 = wrkbl, i__3 = *lda * *n + *m;
					if (*lwork >= MAX (i__2, i__3) + *lda * *m) {

						/* WORK(IU) is LDA by N and WORK(IR) is LDA by M */

						ldwrku = *lda;
						chunk = *n;
						ldwrkr = *lda;
					} else {	/* if(complicated condition) */
						/* Computing MAX */
						i__2 = wrkbl, i__3 = *lda * *n + *m;
						if (*lwork >= MAX (i__2, i__3) + *m * *m) {

							/* WORK(IU) is LDA by N and WORK(IR) is M by M */

							ldwrku = *lda;
							chunk = *n;
							ldwrkr = *m;
						} else {

							/* WORK(IU) is M by CHUNK and WORK(IR) is M by M */

							ldwrku = *m;
							chunk = (*lwork - *m * *m - *m) / *m;
							ldwrkr = *m;
						}
					}
					itau = ir + ldwrkr * *m;
					iwork = itau + *m;

					/* Compute A=L*Q (Workspace: need M*M+2*M, prefer
					   M*M+M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

					/* Copy L to WORK(IR) and zero out above it */

					NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
					i__2 = *m - 1;
					i__3 = *m - 1;
					NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[ir + ldwrkr], &ldwrkr);

					/* Generate Q in A (Workspace: need M*M+2*M, prefer
					   M*M+M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
					ie = itau;
					itauq = ie + *m;
					itaup = itauq + *m;
					iwork = itaup + *m;

					/* Bidiagonalize L in WORK(IR) (Workspace: need M*M+4*M,
					   prefer M*M+3*M+2*M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgebrd (m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__2, &ierr);

					/* Generate right vectors bidiagonalizing L (Workspace:
					   need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__2,
					                  &ierr);
					iwork = ie + *m;

					/* Perform bidiagonal QR iteration, computing right
					   singular vectors of L in WORK(IR) (Workspace: need
					   M*M+BDSPAC) */

					NUMlapack_dbdsqr ("U", m, m, &c__0, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, dum, &c__1,
					                  dum, &c__1, &work[iwork], info);
					iu = ie + *m;

					/* Multiply right singular vectors of L in WORK(IR) by Q
					   in A, storing result in WORK(IU) and copying to A
					   (Workspace: need M*M+2*M, prefer M*M+M*N+M) */

					i__2 = *n;
					i__3 = chunk;
					for (i__ = 1; i__3 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__3) {
						/* Computing MIN */
						i__4 = *n - i__ + 1;
						blk = MIN (i__4, chunk);
						NUMblas_dgemm ("N", "N", m, &blk, m, &c_b438, &work[ir], &ldwrkr, &a_ref (1, i__), lda,
						               &c_b416, &work[iu], &ldwrku);
						NUMlapack_dlacpy ("F", m, &blk, &work[iu], &ldwrku, &a_ref (1, i__), lda);
						/* L30: */
					}

				} else {

					/* Insufficient workspace for a fast algorithm */

					ie = 1;
					itauq = ie + *m;
					itaup = itauq + *m;
					iwork = itaup + *m;

					/* Bidiagonalize A (Workspace: need 3*M+N, prefer
					   3*M+(M+N)*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgebrd (m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__3, &ierr);

					/* Generate right vectors bidiagonalizing A (Workspace:
					   need 4*M, prefer 3*M+M*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", m, n, m, &a[a_offset], lda, &work[itaup], &work[iwork], &i__3,
					                  &ierr);
					iwork = ie + *m;

					/* Perform bidiagonal QR iteration, computing right
					   singular vectors of A in A (Workspace: need BDSPAC) */

					NUMlapack_dbdsqr ("L", m, n, &c__0, &c__0, &s[1], &work[ie], &a[a_offset], lda, dum, &c__1,
					                  dum, &c__1, &work[iwork], info);

				}

			} else if (wntvo && wntuas) {

				/* Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
				   M right singular vectors to be overwritten on A and M left
				   singular vectors to be computed in U

				   Computing MAX */
				i__3 = *m << 2;
				if (*lwork >= *m * *m + MAX (i__3, bdspac)) {

					/* Sufficient workspace for a fast algorithm */

					ir = 1;
					/* Computing MAX */
					i__3 = wrkbl, i__2 = *lda * *n + *m;
					if (*lwork >= MAX (i__3, i__2) + *lda * *m) {

						/* WORK(IU) is LDA by N and WORK(IR) is LDA by M */

						ldwrku = *lda;
						chunk = *n;
						ldwrkr = *lda;
					} else {	/* if(complicated condition) */
						/* Computing MAX */
						i__3 = wrkbl, i__2 = *lda * *n + *m;
						if (*lwork >= MAX (i__3, i__2) + *m * *m) {

							/* WORK(IU) is LDA by N and WORK(IR) is M by M */

							ldwrku = *lda;
							chunk = *n;
							ldwrkr = *m;
						} else {

							/* WORK(IU) is M by CHUNK and WORK(IR) is M by M */

							ldwrku = *m;
							chunk = (*lwork - *m * *m - *m) / *m;
							ldwrkr = *m;
						}
					}
					itau = ir + ldwrkr * *m;
					iwork = itau + *m;

					/* Compute A=L*Q (Workspace: need M*M+2*M, prefer
					   M*M+M+M*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__3, &ierr);

					/* Copy L to U, zeroing about above it */

					NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
					i__3 = *m - 1;
					i__2 = *m - 1;
					NUMlapack_dlaset ("U", &i__3, &i__2, &c_b416, &c_b416, &u_ref (1, 2), ldu);

					/* Generate Q in A (Workspace: need M*M+2*M, prefer
					   M*M+M+M*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__3, &ierr);
					ie = itau;
					itauq = ie + *m;
					itaup = itauq + *m;
					iwork = itaup + *m;

					/* Bidiagonalize L in U, copying result to WORK(IR)
					   (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dgebrd (m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__3, &ierr);
					NUMlapack_dlacpy ("U", m, m, &u[u_offset], ldu, &work[ir], &ldwrkr);

					/* Generate right vectors bidiagonalizing L in WORK(IR)
					   (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__3,
					                  &ierr);

					/* Generate left vectors bidiagonalizing L in U
					   (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */

					i__3 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__3,
					                  &ierr);
					iwork = ie + *m;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of L in U, and computing right
					   singular vectors of L in WORK(IR) (Workspace: need
					   M*M+BDSPAC) */

					NUMlapack_dbdsqr ("U", m, m, m, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, &u[u_offset],
					                  ldu, dum, &c__1, &work[iwork], info);
					iu = ie + *m;

					/* Multiply right singular vectors of L in WORK(IR) by Q
					   in A, storing result in WORK(IU) and copying to A
					   (Workspace: need M*M+2*M, prefer M*M+M*N+M)) */

					i__3 = *n;
					i__2 = chunk;
					for (i__ = 1; i__2 < 0 ? i__ >= i__3 : i__ <= i__3; i__ += i__2) {
						/* Computing MIN */
						i__4 = *n - i__ + 1;
						blk = MIN (i__4, chunk);
						NUMblas_dgemm ("N", "N", m, &blk, m, &c_b438, &work[ir], &ldwrkr, &a_ref (1, i__), lda,
						               &c_b416, &work[iu], &ldwrku);
						NUMlapack_dlacpy ("F", m, &blk, &work[iu], &ldwrku, &a_ref (1, i__), lda);
						/* L40: */
					}

				} else {

					/* Insufficient workspace for a fast algorithm */

					itau = 1;
					iwork = itau + *m;

					/* Compute A=L*Q (Workspace: need 2*M, prefer M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

					/* Copy L to U, zeroing out above it */

					NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
					i__2 = *m - 1;
					i__3 = *m - 1;
					NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &u_ref (1, 2), ldu);

					/* Generate Q in A (Workspace: need 2*M, prefer M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
					ie = itau;
					itauq = ie + *m;
					itaup = itauq + *m;
					iwork = itaup + *m;

					/* Bidiagonalize L in U (Workspace: need 4*M, prefer
					   3*M+2*M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dgebrd (m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[itauq], &work[itaup],
					                  &work[iwork], &i__2, &ierr);

					/* Multiply right vectors bidiagonalizing L by Q in A
					   (Workspace: need 3*M+N, prefer 3*M+N*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dormbr ("P", "L", "T", m, n, m, &u[u_offset], ldu, &work[itaup], &a[a_offset],
					                  lda, &work[iwork], &i__2, &ierr);

					/* Generate left vectors bidiagonalizing L in U
					   (Workspace: need 4*M, prefer 3*M+M*NB) */

					i__2 = *lwork - iwork + 1;
					NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2,
					                  &ierr);
					iwork = ie + *m;

					/* Perform bidiagonal QR iteration, computing left
					   singular vectors of A in U and computing right
					   singular vectors of A in A (Workspace: need BDSPAC) */

					NUMlapack_dbdsqr ("U", m, n, m, &c__0, &s[1], &work[ie], &a[a_offset], lda, &u[u_offset],
					                  ldu, dum, &c__1, &work[iwork], info);

				}

			} else if (wntvs) {

				if (wntun) {

					/* Path 4t(N much larger than M, JOBU='N', JOBVT='S') M
					   right singular vectors to be computed in VT and no
					   left singular vectors to be computed

					   Computing MAX */
					i__2 = *m << 2;
					if (*lwork >= *m * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						ir = 1;
						if (*lwork >= wrkbl + *lda * *m) {

							/* WORK(IR) is LDA by M */

							ldwrkr = *lda;
						} else {

							/* WORK(IR) is M by M */

							ldwrkr = *m;
						}
						itau = ir + ldwrkr * *m;
						iwork = itau + *m;

						/* Compute A=L*Q (Workspace: need M*M+2*M, prefer
						   M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy L to WORK(IR), zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[ir + ldwrkr], &ldwrkr);

						/* Generate Q in A (Workspace: need M*M+2*M, prefer
						   M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IR) (Workspace: need
						   M*M+4*M, prefer M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Generate right vectors bidiagonalizing L in
						   WORK(IR) (Workspace: need M*M+4*M, prefer
						   M*M+3*M+(M-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing right
						   singular vectors of L in WORK(IR) (Workspace: need
						   M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, &c__0, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, dum,
						                  &c__1, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IR)
						   by Q in A, storing result in VT (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[ir], &ldwrkr, &a[a_offset], lda, &c_b416,
						               &vt[vt_offset], ldvt);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q (Workspace: need 2*M, prefer M+M*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy result to VT */

						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need 2*M, prefer
						   M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Zero out above L in A */

						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &a_ref (1, 2), lda);

						/* Bidiagonalize L in A (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right vectors bidiagonalizing L by Q in
						   VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &a[a_offset], lda, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, &c__0, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt, dum,
						                  &c__1, dum, &c__1, &work[iwork], info);

					}

				} else if (wntuo) {

					/* Path 5t(N much larger than M, JOBU='O', JOBVT='S') M
					   right singular vectors to be computed in VT and M left
					   singular vectors to be overwritten on A

					   Computing MAX */
					i__2 = *m << 2;
					if (*lwork >= (*m << 1) * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + (*lda << 1) * *m) {

							/* WORK(IU) is LDA by M and WORK(IR) is LDA by M */

							ldwrku = *lda;
							ir = iu + ldwrku * *m;
							ldwrkr = *lda;
						} else if (*lwork >= wrkbl + (*lda + *m) * *m) {

							/* WORK(IU) is LDA by M and WORK(IR) is M by M */

							ldwrku = *lda;
							ir = iu + ldwrku * *m;
							ldwrkr = *m;
						} else {

							/* WORK(IU) is M by M and WORK(IR) is M by M */

							ldwrku = *m;
							ir = iu + ldwrku * *m;
							ldwrkr = *m;
						}
						itau = ir + ldwrkr * *m;
						iwork = itau + *m;

						/* Compute A=L*Q (Workspace: need 2*M*M+2*M, prefer
						   2*M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy L to WORK(IU), zeroing out below it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[iu + ldwrku], &ldwrku);

						/* Generate Q in A (Workspace: need 2*M*M+2*M, prefer
						   2*M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IU), copying result to
						   WORK(IR) (Workspace: need 2*M*M+4*M, prefer
						   2*M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, m, &work[iu], &ldwrku, &work[ir], &ldwrkr);

						/* Generate right bidiagonalizing vectors in WORK(IU)
						   (Workspace: need 2*M*M+4*M-1, prefer
						   2*M*M+3*M+(M-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[iu], &ldwrku, &work[itaup], &work[iwork], &i__2,
						                  &ierr);

						/* Generate left bidiagonalizing vectors in WORK(IR)
						   (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of L in WORK(IR) and computing
						   right singular vectors of L in WORK(IU)
						   (Workspace: need 2*M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, m, &c__0, &s[1], &work[ie], &work[iu], &ldwrku, &work[ir],
						                  &ldwrkr, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IU)
						   by Q in A, storing result in VT (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[iu], &ldwrku, &a[a_offset], lda, &c_b416,
						               &vt[vt_offset], ldvt);

						/* Copy left singular vectors of L to A (Workspace:
						   need M*M) */

						NUMlapack_dlacpy ("F", m, m, &work[ir], &ldwrkr, &a[a_offset], lda);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M, prefer M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need 2*M, prefer
						   M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Zero out above L in A */

						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &a_ref (1, 2), lda);

						/* Bidiagonalize L in A (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right vectors bidiagonalizing L by Q in
						   VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &a[a_offset], lda, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);

						/* Generate left bidiagonalizing vectors of L in A
						   (Workspace: need 4*M, prefer 3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, compute left
						   singular vectors of A in A and compute right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &a[a_offset], lda, dum, &c__1, &work[iwork], info);

					}

				} else if (wntuas) {

					/* Path 6t(N much larger than M, JOBU='S' or 'A',
					   JOBVT='S') M right singular vectors to be computed in
					   VT and M left singular vectors to be computed in U

					   Computing MAX */
					i__2 = *m << 2;
					if (*lwork >= *m * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + *lda * *m) {

							/* WORK(IU) is LDA by N */

							ldwrku = *lda;
						} else {

							/* WORK(IU) is LDA by M */

							ldwrku = *m;
						}
						itau = iu + ldwrku * *m;
						iwork = itau + *m;

						/* Compute A=L*Q (Workspace: need M*M+2*M, prefer
						   M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);

						/* Copy L to WORK(IU), zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[iu + ldwrku], &ldwrku);

						/* Generate Q in A (Workspace: need M*M+2*M, prefer
						   M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IU), copying result to U
						   (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, m, &work[iu], &ldwrku, &u[u_offset], ldu);

						/* Generate right bidiagonalizing vectors in WORK(IU)
						   (Workspace: need M*M+4*M-1, prefer
						   M*M+3*M+(M-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[iu], &ldwrku, &work[itaup], &work[iwork], &i__2,
						                  &ierr);

						/* Generate left bidiagonalizing vectors in U
						   (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of L in U and computing right
						   singular vectors of L in WORK(IU) (Workspace: need
						   M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, m, &c__0, &s[1], &work[ie], &work[iu], &ldwrku,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IU)
						   by Q in A, storing result in VT (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[iu], &ldwrku, &a[a_offset], lda, &c_b416,
						               &vt[vt_offset], ldvt);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M, prefer M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need 2*M, prefer
						   M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (m, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);

						/* Copy L to U, zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &u_ref (1, 2), ldu);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in U (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right bidiagonalizing vectors in U by Q
						   in VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &u[u_offset], ldu, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);

						/* Generate left bidiagonalizing vectors in U
						   (Workspace: need 4*M, prefer 3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				}

			} else if (wntva) {

				if (wntun) {

					/* Path 7t(N much larger than M, JOBU='N', JOBVT='A') N
					   right singular vectors to be computed in VT and no
					   left singular vectors to be computed

					   Computing MAX */
					i__2 = *n + *m, i__3 = *m << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= *m * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						ir = 1;
						if (*lwork >= wrkbl + *lda * *m) {

							/* WORK(IR) is LDA by M */

							ldwrkr = *lda;
						} else {

							/* WORK(IR) is M by M */

							ldwrkr = *m;
						}
						itau = ir + ldwrkr * *m;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need M*M+2*M, prefer M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Copy L to WORK(IR), zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[ir], &ldwrkr);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[ir + ldwrkr], &ldwrkr);

						/* Generate Q in VT (Workspace: need M*M+M+N, prefer
						   M*M+M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IR) (Workspace: need
						   M*M+4*M, prefer M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[ir], &ldwrkr, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);

						/* Generate right bidiagonalizing vectors in WORK(IR)
						   (Workspace: need M*M+4*M-1, prefer
						   M*M+3*M+(M-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[ir], &ldwrkr, &work[itaup], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing right
						   singular vectors of L in WORK(IR) (Workspace: need
						   M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, &c__0, &c__0, &s[1], &work[ie], &work[ir], &ldwrkr, dum,
						                  &c__1, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IR)
						   by Q in VT, storing result in A (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[ir], &ldwrkr, &vt[vt_offset], ldvt, &c_b416,
						               &a[a_offset], lda);

						/* Copy right singular vectors of A from A to VT */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M, prefer M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need M+N, prefer
						   M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Zero out above L in A */

						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &a_ref (1, 2), lda);

						/* Bidiagonalize L in A (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right bidiagonalizing vectors in A by Q
						   in VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &a[a_offset], lda, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, &c__0, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt, dum,
						                  &c__1, dum, &c__1, &work[iwork], info);

					}

				} else if (wntuo) {

					/* Path 8t(N much larger than M, JOBU='O', JOBVT='A') N
					   right singular vectors to be computed in VT and M left
					   singular vectors to be overwritten on A

					   Computing MAX */
					i__2 = *n + *m, i__3 = *m << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= (*m << 1) * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + (*lda << 1) * *m) {

							/* WORK(IU) is LDA by M and WORK(IR) is LDA by M */

							ldwrku = *lda;
							ir = iu + ldwrku * *m;
							ldwrkr = *lda;
						} else if (*lwork >= wrkbl + (*lda + *m) * *m) {

							/* WORK(IU) is LDA by M and WORK(IR) is M by M */

							ldwrku = *lda;
							ir = iu + ldwrku * *m;
							ldwrkr = *m;
						} else {

							/* WORK(IU) is M by M and WORK(IR) is M by M */

							ldwrku = *m;
							ir = iu + ldwrku * *m;
							ldwrkr = *m;
						}
						itau = ir + ldwrkr * *m;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M*M+2*M, prefer 2*M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need 2*M*M+M+N,
						   prefer 2*M*M+M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);

						/* Copy L to WORK(IU), zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[iu + ldwrku], &ldwrku);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IU), copying result to
						   WORK(IR) (Workspace: need 2*M*M+4*M, prefer
						   2*M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, m, &work[iu], &ldwrku, &work[ir], &ldwrkr);

						/* Generate right bidiagonalizing vectors in WORK(IU)
						   (Workspace: need 2*M*M+4*M-1, prefer
						   2*M*M+3*M+(M-1)*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[iu], &ldwrku, &work[itaup], &work[iwork], &i__2,
						                  &ierr);

						/* Generate left bidiagonalizing vectors in WORK(IR)
						   (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &work[ir], &ldwrkr, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of L in WORK(IR) and computing
						   right singular vectors of L in WORK(IU)
						   (Workspace: need 2*M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, m, &c__0, &s[1], &work[ie], &work[iu], &ldwrku, &work[ir],
						                  &ldwrkr, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IU)
						   by Q in VT, storing result in A (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[iu], &ldwrku, &vt[vt_offset], ldvt, &c_b416,
						               &a[a_offset], lda);

						/* Copy right singular vectors of A from A to VT */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Copy left singular vectors of A from WORK(IR) to A
						 */

						NUMlapack_dlacpy ("F", m, m, &work[ir], &ldwrkr, &a[a_offset], lda);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M, prefer M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need M+N, prefer
						   M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Zero out above L in A */

						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &a_ref (1, 2), lda);

						/* Bidiagonalize L in A (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right bidiagonalizing vectors in A by Q
						   in VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &a[a_offset], lda, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);

						/* Generate left bidiagonalizing vectors in A
						   (Workspace: need 4*M, prefer 3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &a[a_offset], lda, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in A and computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &a[a_offset], lda, dum, &c__1, &work[iwork], info);

					}

				} else if (wntuas) {

					/* Path 9t(N much larger than M, JOBU='S' or 'A',
					   JOBVT='A') N right singular vectors to be computed in
					   VT and M left singular vectors to be computed in U

					   Computing MAX */
					i__2 = *n + *m, i__3 = *m << 2, i__2 = MAX (i__2, i__3);
					if (*lwork >= *m * *m + MAX (i__2, bdspac)) {

						/* Sufficient workspace for a fast algorithm */

						iu = 1;
						if (*lwork >= wrkbl + *lda * *m) {

							/* WORK(IU) is LDA by M */

							ldwrku = *lda;
						} else {

							/* WORK(IU) is M by M */

							ldwrku = *m;
						}
						itau = iu + ldwrku * *m;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need M*M+2*M, prefer M*M+M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need M*M+M+N, prefer
						   M*M+M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);

						/* Copy L to WORK(IU), zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &work[iu], &ldwrku);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &work[iu + ldwrku], &ldwrku);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in WORK(IU), copying result to U
						   (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &work[iu], &ldwrku, &s[1], &work[ie], &work[itauq],
						                  &work[itaup], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("L", m, m, &work[iu], &ldwrku, &u[u_offset], ldu);

						/* Generate right bidiagonalizing vectors in WORK(IU)
						   (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
						 */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("P", m, m, m, &work[iu], &ldwrku, &work[itaup], &work[iwork], &i__2,
						                  &ierr);

						/* Generate left bidiagonalizing vectors in U
						   (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of L in U and computing right
						   singular vectors of L in WORK(IU) (Workspace: need
						   M*M+BDSPAC) */

						NUMlapack_dbdsqr ("U", m, m, m, &c__0, &s[1], &work[ie], &work[iu], &ldwrku,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

						/* Multiply right singular vectors of L in WORK(IU)
						   by Q in VT, storing result in A (Workspace: need
						   M*M) */

						NUMblas_dgemm ("N", "N", m, n, m, &c_b438, &work[iu], &ldwrku, &vt[vt_offset], ldvt, &c_b416,
						               &a[a_offset], lda);

						/* Copy right singular vectors of A from A to VT */

						NUMlapack_dlacpy ("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

					} else {

						/* Insufficient workspace for a fast algorithm */

						itau = 1;
						iwork = itau + *m;

						/* Compute A=L*Q, copying result to VT (Workspace:
						   need 2*M, prefer M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgelqf (m, n, &a[a_offset], lda, &work[itau], &work[iwork], &i__2, &ierr);
						NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);

						/* Generate Q in VT (Workspace: need M+N, prefer
						   M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorglq (n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[iwork], &i__2,
						                  &ierr);

						/* Copy L to U, zeroing out above it */

						NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
						i__2 = *m - 1;
						i__3 = *m - 1;
						NUMlapack_dlaset ("U", &i__2, &i__3, &c_b416, &c_b416, &u_ref (1, 2), ldu);
						ie = itau;
						itauq = ie + *m;
						itaup = itauq + *m;
						iwork = itaup + *m;

						/* Bidiagonalize L in U (Workspace: need 4*M, prefer
						   3*M+2*M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dgebrd (m, m, &u[u_offset], ldu, &s[1], &work[ie], &work[itauq], &work[itaup],
						                  &work[iwork], &i__2, &ierr);

						/* Multiply right bidiagonalizing vectors in U by Q
						   in VT (Workspace: need 3*M+N, prefer 3*M+N*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dormbr ("P", "L", "T", m, n, m, &u[u_offset], ldu, &work[itaup],
						                  &vt[vt_offset], ldvt, &work[iwork], &i__2, &ierr);

						/* Generate left bidiagonalizing vectors in U
						   (Workspace: need 4*M, prefer 3*M+M*NB) */

						i__2 = *lwork - iwork + 1;
						NUMlapack_dorgbr ("Q", m, m, m, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2,
						                  &ierr);
						iwork = ie + *m;

						/* Perform bidiagonal QR iteration, computing left
						   singular vectors of A in U and computing right
						   singular vectors of A in VT (Workspace: need
						   BDSPAC) */

						NUMlapack_dbdsqr ("U", m, n, m, &c__0, &s[1], &work[ie], &vt[vt_offset], ldvt,
						                  &u[u_offset], ldu, dum, &c__1, &work[iwork], info);

					}

				}

			}

		} else {

			/* N .LT. MNTHR

			   Path 10t(N greater than M, but not much larger) Reduce to
			   bidiagonal form without LQ decomposition */

			ie = 1;
			itauq = ie + *m;
			itaup = itauq + *m;
			iwork = itaup + *m;

			/* Bidiagonalize A (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */

			i__2 = *lwork - iwork + 1;
			NUMlapack_dgebrd (m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &work[itaup],
			                  &work[iwork], &i__2, &ierr);
			if (wntuas) {

				/* If left singular vectors desired in U, copy result to U
				   and generate left bidiagonalizing vectors in U (Workspace:
				   need 4*M-1, prefer 3*M+(M-1)*NB) */

				NUMlapack_dlacpy ("L", m, m, &a[a_offset], lda, &u[u_offset], ldu);
				i__2 = *lwork - iwork + 1;
				NUMlapack_dorgbr ("Q", m, m, n, &u[u_offset], ldu, &work[itauq], &work[iwork], &i__2, &ierr);
			}
			if (wntvas) {

				/* If right singular vectors desired in VT, copy result to VT
				   and generate right bidiagonalizing vectors in VT
				   (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB) */

				NUMlapack_dlacpy ("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
				if (wntva) {
					nrvt = *n;
				}
				if (wntvs)