3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DGESDD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
21 * SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22 * WORK, LWORK, IWORK, INFO )
24 * .. Scalar Arguments ..
26 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
28 * .. Array Arguments ..
30 * DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
31 * $ VT( LDVT, * ), WORK( * )
40 *> DGESDD computes the singular value decomposition (SVD) of a real
41 *> M-by-N matrix A, optionally computing the left and right singular
42 *> vectors. If singular vectors are desired, it uses a
43 *> divide-and-conquer algorithm.
47 *> A = U * SIGMA * transpose(V)
49 *> where SIGMA is an M-by-N matrix which is zero except for its
50 *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
51 *> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
52 *> are the singular values of A; they are real and non-negative, and
53 *> are returned in descending order. The first min(m,n) columns of
54 *> U and V are the left and right singular vectors of A.
56 *> Note that the routine returns VT = V**T, not V.
58 *> The divide and conquer algorithm makes very mild assumptions about
59 *> floating point arithmetic. It will work on machines with a guard
60 *> digit in add/subtract, or on those binary machines without guard
61 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
62 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
63 *> without guard digits, but we know of none.
71 *> JOBZ is CHARACTER*1
72 *> Specifies options for computing all or part of the matrix U:
73 *> = 'A': all M columns of U and all N rows of V**T are
74 *> returned in the arrays U and VT;
75 *> = 'S': the first min(M,N) columns of U and the first
76 *> min(M,N) rows of V**T are returned in the arrays U
78 *> = 'O': If M >= N, the first N columns of U are overwritten
79 *> on the array A and all rows of V**T are returned in
81 *> otherwise, all columns of U are returned in the
82 *> array U and the first M rows of V**T are overwritten
84 *> = 'N': no columns of U or rows of V**T are computed.
90 *> The number of rows of the input matrix A. M >= 0.
96 *> The number of columns of the input matrix A. N >= 0.
101 *> A is DOUBLE PRECISION array, dimension (LDA,N)
102 *> On entry, the M-by-N matrix A.
104 *> if JOBZ = 'O', A is overwritten with the first N columns
105 *> of U (the left singular vectors, stored
106 *> columnwise) if M >= N;
107 *> A is overwritten with the first M rows
108 *> of V**T (the right singular vectors, stored
109 *> rowwise) otherwise.
110 *> if JOBZ .ne. 'O', the contents of A are destroyed.
116 *> The leading dimension of the array A. LDA >= max(1,M).
121 *> S is DOUBLE PRECISION array, dimension (min(M,N))
122 *> The singular values of A, sorted so that S(i) >= S(i+1).
127 *> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
128 *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
129 *> UCOL = min(M,N) if JOBZ = 'S'.
130 *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
131 *> orthogonal matrix U;
132 *> if JOBZ = 'S', U contains the first min(M,N) columns of U
133 *> (the left singular vectors, stored columnwise);
134 *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
140 *> The leading dimension of the array U. LDU >= 1; if
141 *> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
146 *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
147 *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
148 *> N-by-N orthogonal matrix V**T;
149 *> if JOBZ = 'S', VT contains the first min(M,N) rows of
150 *> V**T (the right singular vectors, stored rowwise);
151 *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
157 *> The leading dimension of the array VT. LDVT >= 1;
158 *> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
159 *> if JOBZ = 'S', LDVT >= min(M,N).
164 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
171 *> The dimension of the array WORK. LWORK >= 1.
172 *> If LWORK = -1, a workspace query is assumed. The optimal
173 *> size for the WORK array is calculated and stored in WORK(1),
174 *> and no other work except argument checking is performed.
176 *> Let mx = max(M,N) and mn = min(M,N).
177 *> If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
178 *> If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
179 *> If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
180 *> If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
181 *> These are not tight minimums in all cases; see comments inside code.
182 *> For good performance, LWORK should generally be larger;
183 *> a query is recommended.
188 *> IWORK is INTEGER array, dimension (8*min(M,N))
194 *> = 0: successful exit.
195 *> < 0: if INFO = -i, the i-th argument had an illegal value.
196 *> > 0: DBDSDC did not converge, updating process failed.
202 *> \author Univ. of Tennessee
203 *> \author Univ. of California Berkeley
204 *> \author Univ. of Colorado Denver
209 *> \ingroup doubleGEsing
211 *> \par Contributors:
214 *> Ming Gu and Huan Ren, Computer Science Division, University of
215 *> California at Berkeley, USA
217 *> @precisions fortran d -> s
218 * =====================================================================
219 SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
220 $ WORK, LWORK, IWORK, INFO )
223 * -- LAPACK driver routine (version 3.6.1) --
224 * -- LAPACK is a software package provided by Univ. of Tennessee, --
225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
228 * .. Scalar Arguments ..
230 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
232 * .. Array Arguments ..
234 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
235 $ VT( LDVT, * ), WORK( * )
238 * =====================================================================
241 DOUBLE PRECISION ZERO, ONE
242 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
244 * .. Local Scalars ..
245 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
246 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
247 $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
248 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
249 $ MNTHR, NWORK, WRKBL
250 INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
251 $ LWORK_DGEBRD_NN, LWORK_DGELQF_MN,
253 $ LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN,
254 $ LWORK_DORGLQ_MN, LWORK_DORGLQ_NN,
255 $ LWORK_DORGQR_MM, LWORK_DORGQR_MN,
256 $ LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM,
257 $ LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN,
258 $ LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN
259 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
263 DOUBLE PRECISION DUM( 1 )
265 * .. External Subroutines ..
266 EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
267 $ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
270 * .. External Functions ..
272 DOUBLE PRECISION DLAMCH, DLANGE
273 EXTERNAL DLAMCH, DLANGE, LSAME
275 * .. Intrinsic Functions ..
276 INTRINSIC INT, MAX, MIN, SQRT
278 * .. Executable Statements ..
280 * Test the input arguments
284 WNTQA = LSAME( JOBZ, 'A' )
285 WNTQS = LSAME( JOBZ, 'S' )
286 WNTQAS = WNTQA .OR. WNTQS
287 WNTQO = LSAME( JOBZ, 'O' )
288 WNTQN = LSAME( JOBZ, 'N' )
289 LQUERY = ( LWORK.EQ.-1 )
291 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
293 ELSE IF( M.LT.0 ) THEN
295 ELSE IF( N.LT.0 ) THEN
297 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
299 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
300 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
302 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
303 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
304 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
309 * Note: Comments in the code beginning "Workspace:" describe the
310 * minimal amount of workspace allocated at that point in the code,
311 * as well as the preferred amount for good performance.
312 * NB refers to the optimal block size for the immediately
313 * following subroutine, as returned by ILAENV.
319 MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
320 IF( M.GE.N .AND. MINMN.GT.0 ) THEN
322 * Compute space needed for DBDSDC
325 * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
326 * keep 7*N for backwards compatability.
332 * Compute space preferred for each routine
333 CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
334 $ DUM(1), DUM(1), -1, IERR )
335 LWORK_DGEBRD_MN = INT( DUM(1) )
337 CALL DGEBRD( N, N, DUM(1), N, DUM(1), DUM(1), DUM(1),
338 $ DUM(1), DUM(1), -1, IERR )
339 LWORK_DGEBRD_NN = INT( DUM(1) )
341 CALL DGEQRF( M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
342 LWORK_DGEQRF_MN = INT( DUM(1) )
344 CALL DORGBR( 'Q', N, N, N, DUM(1), N, DUM(1), DUM(1), -1,
346 LWORK_DORGBR_Q_NN = INT( DUM(1) )
348 CALL DORGQR( M, M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
349 LWORK_DORGQR_MM = INT( DUM(1) )
351 CALL DORGQR( M, N, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
352 LWORK_DORGQR_MN = INT( DUM(1) )
354 CALL DORMBR( 'P', 'R', 'T', N, N, N, DUM(1), N,
355 $ DUM(1), DUM(1), N, DUM(1), -1, IERR )
356 LWORK_DORMBR_PRT_NN = INT( DUM(1) )
358 CALL DORMBR( 'Q', 'L', 'N', N, N, N, DUM(1), N,
359 $ DUM(1), DUM(1), N, DUM(1), -1, IERR )
360 LWORK_DORMBR_QLN_NN = INT( DUM(1) )
362 CALL DORMBR( 'Q', 'L', 'N', M, N, N, DUM(1), M,
363 $ DUM(1), DUM(1), M, DUM(1), -1, IERR )
364 LWORK_DORMBR_QLN_MN = INT( DUM(1) )
366 CALL DORMBR( 'Q', 'L', 'N', M, M, N, DUM(1), M,
367 $ DUM(1), DUM(1), M, DUM(1), -1, IERR )
368 LWORK_DORMBR_QLN_MM = INT( DUM(1) )
370 IF( M.GE.MNTHR ) THEN
373 * Path 1 (M >> N, JOBZ='N')
375 WRKBL = N + LWORK_DGEQRF_MN
376 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
377 MAXWRK = MAX( WRKBL, BDSPAC + N )
379 ELSE IF( WNTQO ) THEN
381 * Path 2 (M >> N, JOBZ='O')
383 WRKBL = N + LWORK_DGEQRF_MN
384 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN )
385 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
386 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
387 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
388 WRKBL = MAX( WRKBL, 3*N + BDSPAC )
389 MAXWRK = WRKBL + 2*N*N
390 MINWRK = BDSPAC + 2*N*N + 3*N
391 ELSE IF( WNTQS ) THEN
393 * Path 3 (M >> N, JOBZ='S')
395 WRKBL = N + LWORK_DGEQRF_MN
396 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN )
397 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
398 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
399 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
400 WRKBL = MAX( WRKBL, 3*N + BDSPAC )
402 MINWRK = BDSPAC + N*N + 3*N
403 ELSE IF( WNTQA ) THEN
405 * Path 4 (M >> N, JOBZ='A')
407 WRKBL = N + LWORK_DGEQRF_MN
408 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MM )
409 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
410 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
411 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
412 WRKBL = MAX( WRKBL, 3*N + BDSPAC )
414 MINWRK = N*N + MAX( 3*N + BDSPAC, N + M )
418 * Path 5 (M >= N, but not much larger)
420 WRKBL = 3*N + LWORK_DGEBRD_MN
422 * Path 5n (M >= N, jobz='N')
423 MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
424 MINWRK = 3*N + MAX( M, BDSPAC )
425 ELSE IF( WNTQO ) THEN
426 * Path 5o (M >= N, jobz='O')
427 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
428 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
429 WRKBL = MAX( WRKBL, 3*N + BDSPAC )
431 MINWRK = 3*N + MAX( M, N*N + BDSPAC )
432 ELSE IF( WNTQS ) THEN
433 * Path 5s (M >= N, jobz='S')
434 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
435 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
436 MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
437 MINWRK = 3*N + MAX( M, BDSPAC )
438 ELSE IF( WNTQA ) THEN
439 * Path 5a (M >= N, jobz='A')
440 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MM )
441 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
442 MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
443 MINWRK = 3*N + MAX( M, BDSPAC )
446 ELSE IF( MINMN.GT.0 ) THEN
448 * Compute space needed for DBDSDC
451 * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
452 * keep 7*N for backwards compatability.
458 * Compute space preferred for each routine
459 CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
460 $ DUM(1), DUM(1), -1, IERR )
461 LWORK_DGEBRD_MN = INT( DUM(1) )
463 CALL DGEBRD( M, M, A, M, S, DUM(1), DUM(1),
464 $ DUM(1), DUM(1), -1, IERR )
465 LWORK_DGEBRD_MM = INT( DUM(1) )
467 CALL DGELQF( M, N, A, M, DUM(1), DUM(1), -1, IERR )
468 LWORK_DGELQF_MN = INT( DUM(1) )
470 CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
471 LWORK_DORGLQ_NN = INT( DUM(1) )
473 CALL DORGLQ( M, N, M, A, M, DUM(1), DUM(1), -1, IERR )
474 LWORK_DORGLQ_MN = INT( DUM(1) )
476 CALL DORGBR( 'P', M, M, M, A, N, DUM(1), DUM(1), -1, IERR )
477 LWORK_DORGBR_P_MM = INT( DUM(1) )
479 CALL DORMBR( 'P', 'R', 'T', M, M, M, DUM(1), M,
480 $ DUM(1), DUM(1), M, DUM(1), -1, IERR )
481 LWORK_DORMBR_PRT_MM = INT( DUM(1) )
483 CALL DORMBR( 'P', 'R', 'T', M, N, M, DUM(1), M,
484 $ DUM(1), DUM(1), M, DUM(1), -1, IERR )
485 LWORK_DORMBR_PRT_MN = INT( DUM(1) )
487 CALL DORMBR( 'P', 'R', 'T', N, N, M, DUM(1), N,
488 $ DUM(1), DUM(1), N, DUM(1), -1, IERR )
489 LWORK_DORMBR_PRT_NN = INT( DUM(1) )
491 CALL DORMBR( 'Q', 'L', 'N', M, M, M, DUM(1), M,
492 $ DUM(1), DUM(1), M, DUM(1), -1, IERR )
493 LWORK_DORMBR_QLN_MM = INT( DUM(1) )
495 IF( N.GE.MNTHR ) THEN
498 * Path 1t (N >> M, JOBZ='N')
500 WRKBL = M + LWORK_DGELQF_MN
501 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
502 MAXWRK = MAX( WRKBL, BDSPAC + M )
504 ELSE IF( WNTQO ) THEN
506 * Path 2t (N >> M, JOBZ='O')
508 WRKBL = M + LWORK_DGELQF_MN
509 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN )
510 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
511 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
512 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
513 WRKBL = MAX( WRKBL, 3*M + BDSPAC )
514 MAXWRK = WRKBL + 2*M*M
515 MINWRK = BDSPAC + 2*M*M + 3*M
516 ELSE IF( WNTQS ) THEN
518 * Path 3t (N >> M, JOBZ='S')
520 WRKBL = M + LWORK_DGELQF_MN
521 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN )
522 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
523 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
524 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
525 WRKBL = MAX( WRKBL, 3*M + BDSPAC )
527 MINWRK = BDSPAC + M*M + 3*M
528 ELSE IF( WNTQA ) THEN
530 * Path 4t (N >> M, JOBZ='A')
532 WRKBL = M + LWORK_DGELQF_MN
533 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_NN )
534 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
535 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
536 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
537 WRKBL = MAX( WRKBL, 3*M + BDSPAC )
539 MINWRK = M*M + MAX( 3*M + BDSPAC, M + N )
543 * Path 5t (N > M, but not much larger)
545 WRKBL = 3*M + LWORK_DGEBRD_MN
547 * Path 5tn (N > M, jobz='N')
548 MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
549 MINWRK = 3*M + MAX( N, BDSPAC )
550 ELSE IF( WNTQO ) THEN
551 * Path 5to (N > M, jobz='O')
552 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
553 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
554 WRKBL = MAX( WRKBL, 3*M + BDSPAC )
556 MINWRK = 3*M + MAX( N, M*M + BDSPAC )
557 ELSE IF( WNTQS ) THEN
558 * Path 5ts (N > M, jobz='S')
559 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
560 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
561 MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
562 MINWRK = 3*M + MAX( N, BDSPAC )
563 ELSE IF( WNTQA ) THEN
564 * Path 5ta (N > M, jobz='A')
565 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
566 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_NN )
567 MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
568 MINWRK = 3*M + MAX( N, BDSPAC )
573 MAXWRK = MAX( MAXWRK, MINWRK )
576 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
582 CALL XERBLA( 'DGESDD', -INFO )
584 ELSE IF( LQUERY ) THEN
588 * Quick return if possible
590 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
594 * Get machine constants
597 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
598 BIGNUM = ONE / SMLNUM
600 * Scale A if max element outside range [SMLNUM,BIGNUM]
602 ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
604 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
606 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
607 ELSE IF( ANRM.GT.BIGNUM ) THEN
609 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
614 * A has at least as many rows as columns. If A has sufficiently
615 * more rows than columns, first reduce using the QR
616 * decomposition (if sufficient workspace available)
618 IF( M.GE.MNTHR ) THEN
622 * Path 1 (M >> N, JOBZ='N')
623 * No singular vectors to be computed
629 * Workspace: need N [tau] + N [work]
630 * Workspace: prefer N [tau] + N*NB [work]
632 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
633 $ LWORK - NWORK + 1, IERR )
637 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
643 * Bidiagonalize R in A
644 * Workspace: need 3*N [e, tauq, taup] + N [work]
645 * Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
647 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
648 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
652 * Perform bidiagonal SVD, computing singular values only
653 * Workspace: need N [e] + BDSPAC
655 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
656 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
658 ELSE IF( WNTQO ) THEN
660 * Path 2 (M >> N, JOBZ = 'O')
661 * N left singular vectors to be overwritten on A and
662 * N right singular vectors to be computed in VT
666 * WORK(IR) is LDWRKR by N
668 IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN
671 LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N
677 * Workspace: need N*N [R] + N [tau] + N [work]
678 * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
680 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
681 $ LWORK - NWORK + 1, IERR )
683 * Copy R to WORK(IR), zeroing out below it
685 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
686 CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
690 * Workspace: need N*N [R] + N [tau] + N [work]
691 * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
693 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
694 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
700 * Bidiagonalize R in WORK(IR)
701 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
702 * Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
704 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
705 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
706 $ LWORK - NWORK + 1, IERR )
713 * Perform bidiagonal SVD, computing left singular vectors
714 * of bidiagonal matrix in WORK(IU) and computing right
715 * singular vectors of bidiagonal matrix in VT
716 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
718 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
719 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
722 * Overwrite WORK(IU) by left singular vectors of R
723 * and VT by right singular vectors of R
724 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N [work]
725 * Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
727 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
728 $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
729 $ LWORK - NWORK + 1, IERR )
730 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
731 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
732 $ LWORK - NWORK + 1, IERR )
734 * Multiply Q in A by left singular vectors of R in
735 * WORK(IU), storing result in WORK(IR) and copying to A
736 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U]
737 * Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U]
739 DO 10 I = 1, M, LDWRKR
740 CHUNK = MIN( M - I + 1, LDWRKR )
741 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
742 $ LDA, WORK( IU ), N, ZERO, WORK( IR ),
744 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
748 ELSE IF( WNTQS ) THEN
750 * Path 3 (M >> N, JOBZ='S')
751 * N left singular vectors to be computed in U and
752 * N right singular vectors to be computed in VT
763 * Workspace: need N*N [R] + N [tau] + N [work]
764 * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
766 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
767 $ LWORK - NWORK + 1, IERR )
769 * Copy R to WORK(IR), zeroing out below it
771 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
772 CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
776 * Workspace: need N*N [R] + N [tau] + N [work]
777 * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
779 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
780 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
786 * Bidiagonalize R in WORK(IR)
787 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
788 * Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
790 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
791 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
792 $ LWORK - NWORK + 1, IERR )
794 * Perform bidiagonal SVD, computing left singular vectors
795 * of bidiagoal matrix in U and computing right singular
796 * vectors of bidiagonal matrix in VT
797 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC
799 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
800 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
803 * Overwrite U by left singular vectors of R and VT
804 * by right singular vectors of R
805 * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
806 * Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work]
808 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
809 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
810 $ LWORK - NWORK + 1, IERR )
812 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
813 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
814 $ LWORK - NWORK + 1, IERR )
816 * Multiply Q in A by left singular vectors of R in
817 * WORK(IR), storing result in U
818 * Workspace: need N*N [R]
820 CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
821 CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
822 $ LDWRKR, ZERO, U, LDU )
824 ELSE IF( WNTQA ) THEN
826 * Path 4 (M >> N, JOBZ='A')
827 * M left singular vectors to be computed in U and
828 * N right singular vectors to be computed in VT
838 * Compute A=Q*R, copying result to U
839 * Workspace: need N*N [U] + N [tau] + N [work]
840 * Workspace: prefer N*N [U] + N [tau] + N*NB [work]
842 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
843 $ LWORK - NWORK + 1, IERR )
844 CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
847 * Workspace: need N*N [U] + N [tau] + M [work]
848 * Workspace: prefer N*N [U] + N [tau] + M*NB [work]
849 CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
850 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
852 * Produce R in A, zeroing out other entries
854 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
860 * Bidiagonalize R in A
861 * Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
862 * Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work]
864 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
865 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
868 * Perform bidiagonal SVD, computing left singular vectors
869 * of bidiagonal matrix in WORK(IU) and computing right
870 * singular vectors of bidiagonal matrix in VT
871 * Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC
873 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
874 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
877 * Overwrite WORK(IU) by left singular vectors of R and VT
878 * by right singular vectors of R
879 * Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
880 * Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work]
882 CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
883 $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
884 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
885 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
886 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
887 $ LWORK - NWORK + 1, IERR )
889 * Multiply Q in U by left singular vectors of R in
890 * WORK(IU), storing result in A
891 * Workspace: need N*N [U]
893 CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
894 $ LDWRKU, ZERO, A, LDA )
896 * Copy left singular vectors of A from A to U
898 CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
906 * Path 5 (M >= N, but not much larger)
907 * Reduce to bidiagonal form without QR decomposition
915 * Workspace: need 3*N [e, tauq, taup] + M [work]
916 * Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
918 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
919 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
923 * Path 5n (M >= N, JOBZ='N')
924 * Perform bidiagonal SVD, only computing singular values
925 * Workspace: need 3*N [e, tauq, taup] + BDSPAC
927 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
928 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
929 ELSE IF( WNTQO ) THEN
930 * Path 5o (M >= N, JOBZ='O')
932 IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
934 * WORK( IU ) is M by N
937 NWORK = IU + LDWRKU*N
938 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
940 * IR is unused; silence compile warnings
944 * WORK( IU ) is N by N
947 NWORK = IU + LDWRKU*N
949 * WORK(IR) is LDWRKR by N
952 LDWRKR = ( LWORK - N*N - 3*N ) / N
954 NWORK = IU + LDWRKU*N
956 * Perform bidiagonal SVD, computing left singular vectors
957 * of bidiagonal matrix in WORK(IU) and computing right
958 * singular vectors of bidiagonal matrix in VT
959 * Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC
961 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
962 $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
965 * Overwrite VT by right singular vectors of A
966 * Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
967 * Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
969 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
970 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
971 $ LWORK - NWORK + 1, IERR )
973 IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
976 * Overwrite WORK(IU) by left singular vectors of A
977 * Workspace: need 3*N [e, tauq, taup] + M*N [U] + N [work]
978 * Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work]
980 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
981 $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
982 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
984 * Copy left singular vectors of A from WORK(IU) to A
986 CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
991 * Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
992 * Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
994 CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
995 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
997 * Multiply Q in A by left singular vectors of
998 * bidiagonal matrix in WORK(IU), storing result in
999 * WORK(IR) and copying to A
1000 * Workspace: need 3*N [e, tauq, taup] + N*N [U] + NB*N [R]
1001 * Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N [R]
1003 DO 20 I = 1, M, LDWRKR
1004 CHUNK = MIN( M - I + 1, LDWRKR )
1005 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
1006 $ LDA, WORK( IU ), LDWRKU, ZERO,
1007 $ WORK( IR ), LDWRKR )
1008 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
1013 ELSE IF( WNTQS ) THEN
1015 * Path 5s (M >= N, JOBZ='S')
1016 * Perform bidiagonal SVD, computing left singular vectors
1017 * of bidiagonal matrix in U and computing right singular
1018 * vectors of bidiagonal matrix in VT
1019 * Workspace: need 3*N [e, tauq, taup] + BDSPAC
1021 CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
1022 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1023 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1026 * Overwrite U by left singular vectors of A and VT
1027 * by right singular vectors of A
1028 * Workspace: need 3*N [e, tauq, taup] + N [work]
1029 * Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1031 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1032 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1033 $ LWORK - NWORK + 1, IERR )
1034 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
1035 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1036 $ LWORK - NWORK + 1, IERR )
1037 ELSE IF( WNTQA ) THEN
1039 * Path 5a (M >= N, JOBZ='A')
1040 * Perform bidiagonal SVD, computing left singular vectors
1041 * of bidiagonal matrix in U and computing right singular
1042 * vectors of bidiagonal matrix in VT
1043 * Workspace: need 3*N [e, tauq, taup] + BDSPAC
1045 CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
1046 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1047 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1050 * Set the right corner of U to identity matrix
1053 CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
1057 * Overwrite U by left singular vectors of A and VT
1058 * by right singular vectors of A
1059 * Workspace: need 3*N [e, tauq, taup] + M [work]
1060 * Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1062 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1063 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1064 $ LWORK - NWORK + 1, IERR )
1065 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1066 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1067 $ LWORK - NWORK + 1, IERR )
1074 * A has more columns than rows. If A has sufficiently more
1075 * columns than rows, first reduce using the LQ decomposition (if
1076 * sufficient workspace available)
1078 IF( N.GE.MNTHR ) THEN
1082 * Path 1t (N >> M, JOBZ='N')
1083 * No singular vectors to be computed
1089 * Workspace: need M [tau] + M [work]
1090 * Workspace: prefer M [tau] + M*NB [work]
1092 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1093 $ LWORK - NWORK + 1, IERR )
1097 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1103 * Bidiagonalize L in A
1104 * Workspace: need 3*M [e, tauq, taup] + M [work]
1105 * Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1107 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1108 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1112 * Perform bidiagonal SVD, computing singular values only
1113 * Workspace: need M [e] + BDSPAC
1115 CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1116 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1118 ELSE IF( WNTQO ) THEN
1120 * Path 2t (N >> M, JOBZ='O')
1121 * M right singular vectors to be overwritten on A and
1122 * M left singular vectors to be computed in U
1126 * WORK(IVT) is M by M
1127 * WORK(IL) is M by M; it is later resized to M by chunk for gemm
1130 IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN
1135 CHUNK = ( LWORK - M*M ) / M
1137 ITAU = IL + LDWRKL*M
1141 * Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1142 * Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1144 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1145 $ LWORK - NWORK + 1, IERR )
1147 * Copy L to WORK(IL), zeroing about above it
1149 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1150 CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1151 $ WORK( IL + LDWRKL ), LDWRKL )
1154 * Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1155 * Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1157 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1158 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1164 * Bidiagonalize L in WORK(IL)
1165 * Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1166 * Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1168 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1169 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1170 $ LWORK - NWORK + 1, IERR )
1172 * Perform bidiagonal SVD, computing left singular vectors
1173 * of bidiagonal matrix in U, and computing right singular
1174 * vectors of bidiagonal matrix in WORK(IVT)
1175 * Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1177 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1178 $ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
1181 * Overwrite U by left singular vectors of L and WORK(IVT)
1182 * by right singular vectors of L
1183 * Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1184 * Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1186 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1187 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1188 $ LWORK - NWORK + 1, IERR )
1189 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1190 $ WORK( ITAUP ), WORK( IVT ), M,
1191 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1193 * Multiply right singular vectors of L in WORK(IVT) by Q
1194 * in A, storing result in WORK(IL) and copying to A
1195 * Workspace: need M*M [VT] + M*M [L]
1196 * Workspace: prefer M*M [VT] + M*N [L]
1197 * At this point, L is resized as M by chunk.
1199 DO 30 I = 1, N, CHUNK
1200 BLK = MIN( N - I + 1, CHUNK )
1201 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1202 $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1203 CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1207 ELSE IF( WNTQS ) THEN
1209 * Path 3t (N >> M, JOBZ='S')
1210 * M right singular vectors to be computed in VT and
1211 * M left singular vectors to be computed in U
1215 * WORK(IL) is M by M
1218 ITAU = IL + LDWRKL*M
1222 * Workspace: need M*M [L] + M [tau] + M [work]
1223 * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1225 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1226 $ LWORK - NWORK + 1, IERR )
1228 * Copy L to WORK(IL), zeroing out above it
1230 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1231 CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1232 $ WORK( IL + LDWRKL ), LDWRKL )
1235 * Workspace: need M*M [L] + M [tau] + M [work]
1236 * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1238 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1239 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1245 * Bidiagonalize L in WORK(IU).
1246 * Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1247 * Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1249 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1250 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1251 $ LWORK - NWORK + 1, IERR )
1253 * Perform bidiagonal SVD, computing left singular vectors
1254 * of bidiagonal matrix in U and computing right singular
1255 * vectors of bidiagonal matrix in VT
1256 * Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1258 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1259 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1262 * Overwrite U by left singular vectors of L and VT
1263 * by right singular vectors of L
1264 * Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1265 * Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1267 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1268 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1269 $ LWORK - NWORK + 1, IERR )
1270 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1271 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1272 $ LWORK - NWORK + 1, IERR )
1274 * Multiply right singular vectors of L in WORK(IL) by
1275 * Q in A, storing result in VT
1276 * Workspace: need M*M [L]
1278 CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1279 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1280 $ A, LDA, ZERO, VT, LDVT )
1282 ELSE IF( WNTQA ) THEN
1284 * Path 4t (N >> M, JOBZ='A')
1285 * N right singular vectors to be computed in VT and
1286 * M left singular vectors to be computed in U
1290 * WORK(IVT) is M by M
1293 ITAU = IVT + LDWKVT*M
1296 * Compute A=L*Q, copying result to VT
1297 * Workspace: need M*M [VT] + M [tau] + M [work]
1298 * Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1300 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1301 $ LWORK - NWORK + 1, IERR )
1302 CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1305 * Workspace: need M*M [VT] + M [tau] + N [work]
1306 * Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1308 CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1309 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1311 * Produce L in A, zeroing out other entries
1313 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1319 * Bidiagonalize L in A
1320 * Workspace: need M*M [VT] + 3*M [e, tauq, taup] + M [work]
1321 * Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work]
1323 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1324 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1327 * Perform bidiagonal SVD, computing left singular vectors
1328 * of bidiagonal matrix in U and computing right singular
1329 * vectors of bidiagonal matrix in WORK(IVT)
1330 * Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1332 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1333 $ WORK( IVT ), LDWKVT, DUM, IDUM,
1334 $ WORK( NWORK ), IWORK, INFO )
1336 * Overwrite U by left singular vectors of L and WORK(IVT)
1337 * by right singular vectors of L
1338 * Workspace: need M*M [VT] + 3*M [e, tauq, taup]+ M [work]
1339 * Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work]
1341 CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1342 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1343 $ LWORK - NWORK + 1, IERR )
1344 CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1345 $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1346 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1348 * Multiply right singular vectors of L in WORK(IVT) by
1349 * Q in VT, storing result in A
1350 * Workspace: need M*M [VT]
1352 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1353 $ VT, LDVT, ZERO, A, LDA )
1355 * Copy right singular vectors of A from A to VT
1357 CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
1365 * Path 5t (N > M, but not much larger)
1366 * Reduce to bidiagonal form without LQ decomposition
1374 * Workspace: need 3*M [e, tauq, taup] + N [work]
1375 * Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1377 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1378 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1382 * Path 5tn (N > M, JOBZ='N')
1383 * Perform bidiagonal SVD, only computing singular values
1384 * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1386 CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1387 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1388 ELSE IF( WNTQO ) THEN
1389 * Path 5to (N > M, JOBZ='O')
1392 IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1394 * WORK( IVT ) is M by N
1396 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1398 NWORK = IVT + LDWKVT*N
1399 * IL is unused; silence compile warnings
1403 * WORK( IVT ) is M by M
1405 NWORK = IVT + LDWKVT*M
1408 * WORK(IL) is M by CHUNK
1410 CHUNK = ( LWORK - M*M - 3*M ) / M
1413 * Perform bidiagonal SVD, computing left singular vectors
1414 * of bidiagonal matrix in U and computing right singular
1415 * vectors of bidiagonal matrix in WORK(IVT)
1416 * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1418 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1419 $ WORK( IVT ), LDWKVT, DUM, IDUM,
1420 $ WORK( NWORK ), IWORK, INFO )
1422 * Overwrite U by left singular vectors of A
1423 * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1424 * Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1426 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1427 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1428 $ LWORK - NWORK + 1, IERR )
1430 IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1433 * Overwrite WORK(IVT) by left singular vectors of A
1434 * Workspace: need 3*M [e, tauq, taup] + M*N [VT] + M [work]
1435 * Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work]
1437 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1438 $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1439 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1441 * Copy right singular vectors of A from WORK(IVT) to A
1443 CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1447 * Generate P**T in A
1448 * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1449 * Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1451 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1452 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1454 * Multiply Q in A by right singular vectors of
1455 * bidiagonal matrix in WORK(IVT), storing result in
1456 * WORK(IL) and copying to A
1457 * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M*NB [L]
1458 * Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N [L]
1460 DO 40 I = 1, N, CHUNK
1461 BLK = MIN( N - I + 1, CHUNK )
1462 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1463 $ LDWKVT, A( 1, I ), LDA, ZERO,
1465 CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1469 ELSE IF( WNTQS ) THEN
1471 * Path 5ts (N > M, JOBZ='S')
1472 * Perform bidiagonal SVD, computing left singular vectors
1473 * of bidiagonal matrix in U and computing right singular
1474 * vectors of bidiagonal matrix in VT
1475 * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1477 CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1478 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1479 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1482 * Overwrite U by left singular vectors of A and VT
1483 * by right singular vectors of A
1484 * Workspace: need 3*M [e, tauq, taup] + M [work]
1485 * Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1487 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1488 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1489 $ LWORK - NWORK + 1, IERR )
1490 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1491 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1492 $ LWORK - NWORK + 1, IERR )
1493 ELSE IF( WNTQA ) THEN
1495 * Path 5ta (N > M, JOBZ='A')
1496 * Perform bidiagonal SVD, computing left singular vectors
1497 * of bidiagonal matrix in U and computing right singular
1498 * vectors of bidiagonal matrix in VT
1499 * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1501 CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1502 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1503 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1506 * Set the right corner of VT to identity matrix
1509 CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1513 * Overwrite U by left singular vectors of A and VT
1514 * by right singular vectors of A
1515 * Workspace: need 3*M [e, tauq, taup] + N [work]
1516 * Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1518 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1519 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1520 $ LWORK - NWORK + 1, IERR )
1521 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1522 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1523 $ LWORK - NWORK + 1, IERR )
1530 * Undo scaling if necessary
1532 IF( ISCL.EQ.1 ) THEN
1533 IF( ANRM.GT.BIGNUM )
1534 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1536 IF( ANRM.LT.SMLNUM )
1537 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1541 * Return optimal workspace in WORK(1)