LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgesdd()

subroutine dgesdd ( character jobz,
integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) s,
double precision, dimension( ldu, * ) u,
integer ldu,
double precision, dimension( ldvt, * ) vt,
integer ldvt,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer info )

DGESDD

Download DGESDD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGESDD computes the singular value decomposition (SVD) of a real
!> M-by-N matrix A, optionally computing the left and right singular
!> vectors.  If singular vectors are desired, it uses a
!> divide-and-conquer algorithm.
!>
!> The SVD is written
!>
!>      A = U * SIGMA * transpose(V)
!>
!> where SIGMA is an M-by-N matrix which is zero except for its
!> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
!> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
!> are the singular values of A; they are real and non-negative, and
!> are returned in descending order.  The first min(m,n) columns of
!> U and V are the left and right singular vectors of A.
!>
!> Note that the routine returns VT = V**T, not V.
!>
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          Specifies options for computing all or part of the matrix U:
!>          = 'A':  all M columns of U and all N rows of V**T are
!>                  returned in the arrays U and VT;
!>          = 'S':  the first min(M,N) columns of U and the first
!>                  min(M,N) rows of V**T are returned in the arrays U
!>                  and VT;
!>          = 'O':  If M >= N, the first N columns of U are overwritten
!>                  on the array A and all rows of V**T are returned in
!>                  the array VT;
!>                  otherwise, all columns of U are returned in the
!>                  array U and the first M rows of V**T are overwritten
!>                  in the array A;
!>          = 'N':  no columns of U or rows of V**T are computed.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit,
!>          if JOBZ = 'O',  A is overwritten with the first N columns
!>                          of U (the left singular vectors, stored
!>                          columnwise) if M >= N;
!>                          A is overwritten with the first M rows
!>                          of V**T (the right singular vectors, stored
!>                          rowwise) otherwise.
!>          if JOBZ .ne. 'O', the contents of A are destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]S
!>          S is DOUBLE PRECISION array, dimension (min(M,N))
!>          The singular values of A, sorted so that S(i) >= S(i+1).
!>
[out]U
!>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
!>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
!>          UCOL = min(M,N) if JOBZ = 'S'.
!>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
!>          orthogonal matrix U;
!>          if JOBZ = 'S', U contains the first min(M,N) columns of U
!>          (the left singular vectors, stored columnwise);
!>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
!>
[in]LDU
!>          LDU is INTEGER
!>          The leading dimension of the array U.  LDU >= 1; if
!>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
!>
[out]VT
!>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
!>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
!>          N-by-N orthogonal matrix V**T;
!>          if JOBZ = 'S', VT contains the first min(M,N) rows of
!>          V**T (the right singular vectors, stored rowwise);
!>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
!>
[in]LDVT
!>          LDVT is INTEGER
!>          The leading dimension of the array VT.  LDVT >= 1;
!>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
!>          if JOBZ = 'S', LDVT >= min(M,N).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= 1.
!>          If LWORK = -1, a workspace query is assumed.  The optimal
!>          size for the WORK array is calculated and stored in WORK(1),
!>          and no other work except argument checking is performed.
!>
!>          Let mx = max(M,N) and mn = min(M,N).
!>          If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
!>          If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
!>          If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
!>          If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
!>          These are not tight minimums in all cases; see comments inside code.
!>          For good performance, LWORK should generally be larger;
!>          a query is recommended.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (8*min(M,N))
!>
[out]INFO
!>          INFO is INTEGER
!>          <  0:  if INFO = -i, the i-th argument had an illegal value.
!>          = -4:  if A had a NAN entry.
!>          >  0:  DBDSDC did not converge, updating process failed.
!>          =  0:  successful exit.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 209 of file dgesdd.f.

211 implicit none
212*
213* -- LAPACK driver routine --
214* -- LAPACK is a software package provided by Univ. of Tennessee, --
215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217* .. Scalar Arguments ..
218 CHARACTER JOBZ
219 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
220* ..
221* .. Array Arguments ..
222 INTEGER IWORK( * )
223 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
224 $ VT( LDVT, * ), WORK( * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 DOUBLE PRECISION ZERO, ONE
231 parameter( zero = 0.0d0, one = 1.0d0 )
232* ..
233* .. Local Scalars ..
234 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
235 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
236 $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
237 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
238 $ MNTHR, NWORK, WRKBL
239 INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
240 $ LWORK_DGEBRD_NN, LWORK_DGELQF_MN,
241 $ LWORK_DGEQRF_MN,
242 $ LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN,
243 $ LWORK_DORGLQ_MN, LWORK_DORGLQ_NN,
244 $ LWORK_DORGQR_MM, LWORK_DORGQR_MN,
245 $ LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM,
246 $ LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN,
247 $ LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN
248 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
249* ..
250* .. Local Arrays ..
251 INTEGER IDUM( 1 )
252 DOUBLE PRECISION DUM( 1 )
253* ..
254* .. External Subroutines ..
255 EXTERNAL dbdsdc, dgebrd, dgelqf, dgemm, dgeqrf,
256 $ dlacpy,
257 $ dlascl, dlaset, dorgbr, dorglq, dorgqr, dormbr,
258 $ xerbla
259* ..
260* .. External Functions ..
261 LOGICAL LSAME, DISNAN
262 DOUBLE PRECISION DLAMCH, DLANGE, DROUNDUP_LWORK
263 EXTERNAL dlamch, dlange, lsame, disnan,
264 $ droundup_lwork
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC int, max, min, sqrt
268* ..
269* .. Executable Statements ..
270*
271* Test the input arguments
272*
273 info = 0
274 minmn = min( m, n )
275 wntqa = lsame( jobz, 'A' )
276 wntqs = lsame( jobz, 'S' )
277 wntqas = wntqa .OR. wntqs
278 wntqo = lsame( jobz, 'O' )
279 wntqn = lsame( jobz, 'N' )
280 lquery = ( lwork.EQ.-1 )
281*
282 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) ) THEN
283 info = -1
284 ELSE IF( m.LT.0 ) THEN
285 info = -2
286 ELSE IF( n.LT.0 ) THEN
287 info = -3
288 ELSE IF( lda.LT.max( 1, m ) ) THEN
289 info = -5
290 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
291 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) ) THEN
292 info = -8
293 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
294 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
295 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) ) THEN
296 info = -10
297 END IF
298*
299* Compute workspace
300* Note: Comments in the code beginning "Workspace:" describe the
301* minimal amount of workspace allocated at that point in the code,
302* as well as the preferred amount for good performance.
303* NB refers to the optimal block size for the immediately
304* following subroutine, as returned by ILAENV.
305*
306 IF( info.EQ.0 ) THEN
307 minwrk = 1
308 maxwrk = 1
309 bdspac = 0
310 mnthr = int( minmn*11.0d0 / 6.0d0 )
311 IF( m.GE.n .AND. minmn.GT.0 ) THEN
312*
313* Compute space needed for DBDSDC
314*
315 IF( wntqn ) THEN
316* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
317* keep 7*N for backwards compatibility.
318 bdspac = 7*n
319 ELSE
320 bdspac = 3*n*n + 4*n
321 END IF
322*
323* Compute space preferred for each routine
324 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
325 $ dum(1), dum(1), -1, ierr )
326 lwork_dgebrd_mn = int( dum(1) )
327*
328 CALL dgebrd( n, n, dum(1), n, dum(1), dum(1), dum(1),
329 $ dum(1), dum(1), -1, ierr )
330 lwork_dgebrd_nn = int( dum(1) )
331*
332 CALL dgeqrf( m, n, dum(1), m, dum(1), dum(1), -1, ierr )
333 lwork_dgeqrf_mn = int( dum(1) )
334*
335 CALL dorgbr( 'Q', n, n, n, dum(1), n, dum(1), dum(1), -1,
336 $ ierr )
337 lwork_dorgbr_q_nn = int( dum(1) )
338*
339 CALL dorgqr( m, m, n, dum(1), m, dum(1), dum(1), -1,
340 $ ierr )
341 lwork_dorgqr_mm = int( dum(1) )
342*
343 CALL dorgqr( m, n, n, dum(1), m, dum(1), dum(1), -1,
344 $ ierr )
345 lwork_dorgqr_mn = int( dum(1) )
346*
347 CALL dormbr( 'P', 'R', 'T', n, n, n, dum(1), n,
348 $ dum(1), dum(1), n, dum(1), -1, ierr )
349 lwork_dormbr_prt_nn = int( dum(1) )
350*
351 CALL dormbr( 'Q', 'L', 'N', n, n, n, dum(1), n,
352 $ dum(1), dum(1), n, dum(1), -1, ierr )
353 lwork_dormbr_qln_nn = int( dum(1) )
354*
355 CALL dormbr( 'Q', 'L', 'N', m, n, n, dum(1), m,
356 $ dum(1), dum(1), m, dum(1), -1, ierr )
357 lwork_dormbr_qln_mn = int( dum(1) )
358*
359 CALL dormbr( 'Q', 'L', 'N', m, m, n, dum(1), m,
360 $ dum(1), dum(1), m, dum(1), -1, ierr )
361 lwork_dormbr_qln_mm = int( dum(1) )
362*
363 IF( m.GE.mnthr ) THEN
364 IF( wntqn ) THEN
365*
366* Path 1 (M >> N, JOBZ='N')
367*
368 wrkbl = n + lwork_dgeqrf_mn
369 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
370 maxwrk = max( wrkbl, bdspac + n )
371 minwrk = bdspac + n
372 ELSE IF( wntqo ) THEN
373*
374* Path 2 (M >> N, JOBZ='O')
375*
376 wrkbl = n + lwork_dgeqrf_mn
377 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
378 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
379 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
380 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
381 wrkbl = max( wrkbl, 3*n + bdspac )
382 maxwrk = wrkbl + 2*n*n
383 minwrk = bdspac + 2*n*n + 3*n
384 ELSE IF( wntqs ) THEN
385*
386* Path 3 (M >> N, JOBZ='S')
387*
388 wrkbl = n + lwork_dgeqrf_mn
389 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
390 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
391 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
392 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
393 wrkbl = max( wrkbl, 3*n + bdspac )
394 maxwrk = wrkbl + n*n
395 minwrk = bdspac + n*n + 3*n
396 ELSE IF( wntqa ) THEN
397*
398* Path 4 (M >> N, JOBZ='A')
399*
400 wrkbl = n + lwork_dgeqrf_mn
401 wrkbl = max( wrkbl, n + lwork_dorgqr_mm )
402 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
403 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
404 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
405 wrkbl = max( wrkbl, 3*n + bdspac )
406 maxwrk = wrkbl + n*n
407 minwrk = n*n + max( 3*n + bdspac, n + m )
408 END IF
409 ELSE
410*
411* Path 5 (M >= N, but not much larger)
412*
413 wrkbl = 3*n + lwork_dgebrd_mn
414 IF( wntqn ) THEN
415* Path 5n (M >= N, jobz='N')
416 maxwrk = max( wrkbl, 3*n + bdspac )
417 minwrk = 3*n + max( m, bdspac )
418 ELSE IF( wntqo ) THEN
419* Path 5o (M >= N, jobz='O')
420 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
421 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
422 wrkbl = max( wrkbl, 3*n + bdspac )
423 maxwrk = wrkbl + m*n
424 minwrk = 3*n + max( m, n*n + bdspac )
425 ELSE IF( wntqs ) THEN
426* Path 5s (M >= N, jobz='S')
427 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
428 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
429 maxwrk = max( wrkbl, 3*n + bdspac )
430 minwrk = 3*n + max( m, bdspac )
431 ELSE IF( wntqa ) THEN
432* Path 5a (M >= N, jobz='A')
433 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mm )
434 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
435 maxwrk = max( wrkbl, 3*n + bdspac )
436 minwrk = 3*n + max( m, bdspac )
437 END IF
438 END IF
439 ELSE IF( minmn.GT.0 ) THEN
440*
441* Compute space needed for DBDSDC
442*
443 IF( wntqn ) THEN
444* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
445* keep 7*N for backwards compatibility.
446 bdspac = 7*m
447 ELSE
448 bdspac = 3*m*m + 4*m
449 END IF
450*
451* Compute space preferred for each routine
452 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
453 $ dum(1), dum(1), -1, ierr )
454 lwork_dgebrd_mn = int( dum(1) )
455*
456 CALL dgebrd( m, m, a, m, s, dum(1), dum(1),
457 $ dum(1), dum(1), -1, ierr )
458 lwork_dgebrd_mm = int( dum(1) )
459*
460 CALL dgelqf( m, n, a, m, dum(1), dum(1), -1, ierr )
461 lwork_dgelqf_mn = int( dum(1) )
462*
463 CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1,
464 $ ierr )
465 lwork_dorglq_nn = int( dum(1) )
466*
467 CALL dorglq( m, n, m, a, m, dum(1), dum(1), -1, ierr )
468 lwork_dorglq_mn = int( dum(1) )
469*
470 CALL dorgbr( 'P', m, m, m, a, n, dum(1), dum(1), -1,
471 $ ierr )
472 lwork_dorgbr_p_mm = int( dum(1) )
473*
474 CALL dormbr( 'P', 'R', 'T', m, m, m, dum(1), m,
475 $ dum(1), dum(1), m, dum(1), -1, ierr )
476 lwork_dormbr_prt_mm = int( dum(1) )
477*
478 CALL dormbr( 'P', 'R', 'T', m, n, m, dum(1), m,
479 $ dum(1), dum(1), m, dum(1), -1, ierr )
480 lwork_dormbr_prt_mn = int( dum(1) )
481*
482 CALL dormbr( 'P', 'R', 'T', n, n, m, dum(1), n,
483 $ dum(1), dum(1), n, dum(1), -1, ierr )
484 lwork_dormbr_prt_nn = int( dum(1) )
485*
486 CALL dormbr( 'Q', 'L', 'N', m, m, m, dum(1), m,
487 $ dum(1), dum(1), m, dum(1), -1, ierr )
488 lwork_dormbr_qln_mm = int( dum(1) )
489*
490 IF( n.GE.mnthr ) THEN
491 IF( wntqn ) THEN
492*
493* Path 1t (N >> M, JOBZ='N')
494*
495 wrkbl = m + lwork_dgelqf_mn
496 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
497 maxwrk = max( wrkbl, bdspac + m )
498 minwrk = bdspac + m
499 ELSE IF( wntqo ) THEN
500*
501* Path 2t (N >> M, JOBZ='O')
502*
503 wrkbl = m + lwork_dgelqf_mn
504 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
505 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
506 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
507 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
508 wrkbl = max( wrkbl, 3*m + bdspac )
509 maxwrk = wrkbl + 2*m*m
510 minwrk = bdspac + 2*m*m + 3*m
511 ELSE IF( wntqs ) THEN
512*
513* Path 3t (N >> M, JOBZ='S')
514*
515 wrkbl = m + lwork_dgelqf_mn
516 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
517 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
518 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
519 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
520 wrkbl = max( wrkbl, 3*m + bdspac )
521 maxwrk = wrkbl + m*m
522 minwrk = bdspac + m*m + 3*m
523 ELSE IF( wntqa ) THEN
524*
525* Path 4t (N >> M, JOBZ='A')
526*
527 wrkbl = m + lwork_dgelqf_mn
528 wrkbl = max( wrkbl, m + lwork_dorglq_nn )
529 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
530 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
531 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
532 wrkbl = max( wrkbl, 3*m + bdspac )
533 maxwrk = wrkbl + m*m
534 minwrk = m*m + max( 3*m + bdspac, m + n )
535 END IF
536 ELSE
537*
538* Path 5t (N > M, but not much larger)
539*
540 wrkbl = 3*m + lwork_dgebrd_mn
541 IF( wntqn ) THEN
542* Path 5tn (N > M, jobz='N')
543 maxwrk = max( wrkbl, 3*m + bdspac )
544 minwrk = 3*m + max( n, bdspac )
545 ELSE IF( wntqo ) THEN
546* Path 5to (N > M, jobz='O')
547 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
548 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
549 wrkbl = max( wrkbl, 3*m + bdspac )
550 maxwrk = wrkbl + m*n
551 minwrk = 3*m + max( n, m*m + bdspac )
552 ELSE IF( wntqs ) THEN
553* Path 5ts (N > M, jobz='S')
554 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
555 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
556 maxwrk = max( wrkbl, 3*m + bdspac )
557 minwrk = 3*m + max( n, bdspac )
558 ELSE IF( wntqa ) THEN
559* Path 5ta (N > M, jobz='A')
560 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
561 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_nn )
562 maxwrk = max( wrkbl, 3*m + bdspac )
563 minwrk = 3*m + max( n, bdspac )
564 END IF
565 END IF
566 END IF
567
568 maxwrk = max( maxwrk, minwrk )
569 work( 1 ) = droundup_lwork( maxwrk )
570*
571 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
572 info = -12
573 END IF
574 END IF
575*
576 IF( info.NE.0 ) THEN
577 CALL xerbla( 'DGESDD', -info )
578 RETURN
579 ELSE IF( lquery ) THEN
580 RETURN
581 END IF
582*
583* Quick return if possible
584*
585 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
586 RETURN
587 END IF
588*
589* Get machine constants
590*
591 eps = dlamch( 'P' )
592 smlnum = sqrt( dlamch( 'S' ) ) / eps
593 bignum = one / smlnum
594*
595* Scale A if max element outside range [SMLNUM,BIGNUM]
596*
597 anrm = dlange( 'M', m, n, a, lda, dum )
598 IF( disnan( anrm ) ) THEN
599 info = -4
600 RETURN
601 END IF
602 iscl = 0
603 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
604 iscl = 1
605 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
606 ELSE IF( anrm.GT.bignum ) THEN
607 iscl = 1
608 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
609 END IF
610*
611 IF( m.GE.n ) THEN
612*
613* A has at least as many rows as columns. If A has sufficiently
614* more rows than columns, first reduce using the QR
615* decomposition (if sufficient workspace available)
616*
617 IF( m.GE.mnthr ) THEN
618*
619 IF( wntqn ) THEN
620*
621* Path 1 (M >> N, JOBZ='N')
622* No singular vectors to be computed
623*
624 itau = 1
625 nwork = itau + n
626*
627* Compute A=Q*R
628* Workspace: need N [tau] + N [work]
629* Workspace: prefer N [tau] + N*NB [work]
630*
631 CALL dgeqrf( m, n, a, lda, work( itau ),
632 $ work( nwork ),
633 $ lwork - nwork + 1, ierr )
634*
635* Zero out below R
636*
637 CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
638 $ lda )
639 ie = 1
640 itauq = ie + n
641 itaup = itauq + n
642 nwork = itaup + n
643*
644* Bidiagonalize R in A
645* Workspace: need 3*N [e, tauq, taup] + N [work]
646* Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
647*
648 CALL dgebrd( n, n, a, lda, s, work( ie ),
649 $ work( itauq ),
650 $ work( itaup ), work( nwork ), lwork-nwork+1,
651 $ ierr )
652 nwork = ie + n
653*
654* Perform bidiagonal SVD, computing singular values only
655* Workspace: need N [e] + BDSPAC
656*
657 CALL dbdsdc( 'U', 'N', n, s, work( ie ), dum, 1, dum,
658 $ 1,
659 $ dum, idum, work( nwork ), iwork, info )
660*
661 ELSE IF( wntqo ) THEN
662*
663* Path 2 (M >> N, JOBZ = 'O')
664* N left singular vectors to be overwritten on A and
665* N right singular vectors to be computed in VT
666*
667 ir = 1
668*
669* WORK(IR) is LDWRKR by N
670*
671 IF( lwork .GE. lda*n + n*n + 3*n + bdspac ) THEN
672 ldwrkr = lda
673 ELSE
674 ldwrkr = ( lwork - n*n - 3*n - bdspac ) / n
675 END IF
676 itau = ir + ldwrkr*n
677 nwork = itau + n
678*
679* Compute A=Q*R
680* Workspace: need N*N [R] + N [tau] + N [work]
681* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
682*
683 CALL dgeqrf( m, n, a, lda, work( itau ),
684 $ work( nwork ),
685 $ lwork - nwork + 1, ierr )
686*
687* Copy R to WORK(IR), zeroing out below it
688*
689 CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
690 CALL dlaset( 'L', n - 1, n - 1, zero, zero,
691 $ work(ir+1),
692 $ ldwrkr )
693*
694* Generate Q in A
695* Workspace: need N*N [R] + N [tau] + N [work]
696* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
697*
698 CALL dorgqr( m, n, n, a, lda, work( itau ),
699 $ work( nwork ), lwork - nwork + 1, ierr )
700 ie = itau
701 itauq = ie + n
702 itaup = itauq + n
703 nwork = itaup + n
704*
705* Bidiagonalize R in WORK(IR)
706* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
707* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
708*
709 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
710 $ work( itauq ), work( itaup ), work( nwork ),
711 $ lwork - nwork + 1, ierr )
712*
713* WORK(IU) is N by N
714*
715 iu = nwork
716 nwork = iu + n*n
717*
718* Perform bidiagonal SVD, computing left singular vectors
719* of bidiagonal matrix in WORK(IU) and computing right
720* singular vectors of bidiagonal matrix in VT
721* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
722*
723 CALL dbdsdc( 'U', 'I', n, s, work( ie ), work( iu ),
724 $ n,
725 $ vt, ldvt, dum, idum, work( nwork ), iwork,
726 $ info )
727*
728* Overwrite WORK(IU) by left singular vectors of R
729* and VT by right singular vectors of R
730* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N [work]
731* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
732*
733 CALL dormbr( 'Q', 'L', 'N', n, n, n, work( ir ),
734 $ ldwrkr,
735 $ work( itauq ), work( iu ), n, work( nwork ),
736 $ lwork - nwork + 1, ierr )
737 CALL dormbr( 'P', 'R', 'T', n, n, n, work( ir ),
738 $ ldwrkr,
739 $ work( itaup ), vt, ldvt, work( nwork ),
740 $ lwork - nwork + 1, ierr )
741*
742* Multiply Q in A by left singular vectors of R in
743* WORK(IU), storing result in WORK(IR) and copying to A
744* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U]
745* Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U]
746*
747 DO 10 i = 1, m, ldwrkr
748 chunk = min( m - i + 1, ldwrkr )
749 CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
750 $ lda, work( iu ), n, zero, work( ir ),
751 $ ldwrkr )
752 CALL dlacpy( 'F', chunk, n, work( ir ), ldwrkr,
753 $ a( i, 1 ), lda )
754 10 CONTINUE
755*
756 ELSE IF( wntqs ) THEN
757*
758* Path 3 (M >> N, JOBZ='S')
759* N left singular vectors to be computed in U and
760* N right singular vectors to be computed in VT
761*
762 ir = 1
763*
764* WORK(IR) is N by N
765*
766 ldwrkr = n
767 itau = ir + ldwrkr*n
768 nwork = itau + n
769*
770* Compute A=Q*R
771* Workspace: need N*N [R] + N [tau] + N [work]
772* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
773*
774 CALL dgeqrf( m, n, a, lda, work( itau ),
775 $ work( nwork ),
776 $ lwork - nwork + 1, ierr )
777*
778* Copy R to WORK(IR), zeroing out below it
779*
780 CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
781 CALL dlaset( 'L', n - 1, n - 1, zero, zero,
782 $ work(ir+1),
783 $ ldwrkr )
784*
785* Generate Q in A
786* Workspace: need N*N [R] + N [tau] + N [work]
787* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
788*
789 CALL dorgqr( m, n, n, a, lda, work( itau ),
790 $ work( nwork ), lwork - nwork + 1, ierr )
791 ie = itau
792 itauq = ie + n
793 itaup = itauq + n
794 nwork = itaup + n
795*
796* Bidiagonalize R in WORK(IR)
797* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
798* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
799*
800 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
801 $ work( itauq ), work( itaup ), work( nwork ),
802 $ lwork - nwork + 1, ierr )
803*
804* Perform bidiagonal SVD, computing left singular vectors
805* of bidiagoal matrix in U and computing right singular
806* vectors of bidiagonal matrix in VT
807* Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC
808*
809 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
810 $ ldvt, dum, idum, work( nwork ), iwork,
811 $ info )
812*
813* Overwrite U by left singular vectors of R and VT
814* by right singular vectors of R
815* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
816* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work]
817*
818 CALL dormbr( 'Q', 'L', 'N', n, n, n, work( ir ),
819 $ ldwrkr,
820 $ work( itauq ), u, ldu, work( nwork ),
821 $ lwork - nwork + 1, ierr )
822*
823 CALL dormbr( 'P', 'R', 'T', n, n, n, work( ir ),
824 $ ldwrkr,
825 $ work( itaup ), vt, ldvt, work( nwork ),
826 $ lwork - nwork + 1, ierr )
827*
828* Multiply Q in A by left singular vectors of R in
829* WORK(IR), storing result in U
830* Workspace: need N*N [R]
831*
832 CALL dlacpy( 'F', n, n, u, ldu, work( ir ), ldwrkr )
833 CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
834 $ work( ir ),
835 $ ldwrkr, zero, u, ldu )
836*
837 ELSE IF( wntqa ) THEN
838*
839* Path 4 (M >> N, JOBZ='A')
840* M left singular vectors to be computed in U and
841* N right singular vectors to be computed in VT
842*
843 iu = 1
844*
845* WORK(IU) is N by N
846*
847 ldwrku = n
848 itau = iu + ldwrku*n
849 nwork = itau + n
850*
851* Compute A=Q*R, copying result to U
852* Workspace: need N*N [U] + N [tau] + N [work]
853* Workspace: prefer N*N [U] + N [tau] + N*NB [work]
854*
855 CALL dgeqrf( m, n, a, lda, work( itau ),
856 $ work( nwork ),
857 $ lwork - nwork + 1, ierr )
858 CALL dlacpy( 'L', m, n, a, lda, u, ldu )
859*
860* Generate Q in U
861* Workspace: need N*N [U] + N [tau] + M [work]
862* Workspace: prefer N*N [U] + N [tau] + M*NB [work]
863 CALL dorgqr( m, m, n, u, ldu, work( itau ),
864 $ work( nwork ), lwork - nwork + 1, ierr )
865*
866* Produce R in A, zeroing out other entries
867*
868 CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
869 $ lda )
870 ie = itau
871 itauq = ie + n
872 itaup = itauq + n
873 nwork = itaup + n
874*
875* Bidiagonalize R in A
876* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
877* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work]
878*
879 CALL dgebrd( n, n, a, lda, s, work( ie ),
880 $ work( itauq ),
881 $ work( itaup ), work( nwork ), lwork-nwork+1,
882 $ ierr )
883*
884* Perform bidiagonal SVD, computing left singular vectors
885* of bidiagonal matrix in WORK(IU) and computing right
886* singular vectors of bidiagonal matrix in VT
887* Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC
888*
889 CALL dbdsdc( 'U', 'I', n, s, work( ie ), work( iu ),
890 $ n,
891 $ vt, ldvt, dum, idum, work( nwork ), iwork,
892 $ info )
893*
894* Overwrite WORK(IU) by left singular vectors of R and VT
895* by right singular vectors of R
896* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
897* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work]
898*
899 CALL dormbr( 'Q', 'L', 'N', n, n, n, a, lda,
900 $ work( itauq ), work( iu ), ldwrku,
901 $ work( nwork ), lwork - nwork + 1, ierr )
902 CALL dormbr( 'P', 'R', 'T', n, n, n, a, lda,
903 $ work( itaup ), vt, ldvt, work( nwork ),
904 $ lwork - nwork + 1, ierr )
905*
906* Multiply Q in U by left singular vectors of R in
907* WORK(IU), storing result in A
908* Workspace: need N*N [U]
909*
910 CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
911 $ work( iu ),
912 $ ldwrku, zero, a, lda )
913*
914* Copy left singular vectors of A from A to U
915*
916 CALL dlacpy( 'F', m, n, a, lda, u, ldu )
917*
918 END IF
919*
920 ELSE
921*
922* M .LT. MNTHR
923*
924* Path 5 (M >= N, but not much larger)
925* Reduce to bidiagonal form without QR decomposition
926*
927 ie = 1
928 itauq = ie + n
929 itaup = itauq + n
930 nwork = itaup + n
931*
932* Bidiagonalize A
933* Workspace: need 3*N [e, tauq, taup] + M [work]
934* Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
935*
936 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
937 $ work( itaup ), work( nwork ), lwork-nwork+1,
938 $ ierr )
939 IF( wntqn ) THEN
940*
941* Path 5n (M >= N, JOBZ='N')
942* Perform bidiagonal SVD, only computing singular values
943* Workspace: need 3*N [e, tauq, taup] + BDSPAC
944*
945 CALL dbdsdc( 'U', 'N', n, s, work( ie ), dum, 1, dum,
946 $ 1,
947 $ dum, idum, work( nwork ), iwork, info )
948 ELSE IF( wntqo ) THEN
949* Path 5o (M >= N, JOBZ='O')
950 iu = nwork
951 IF( lwork .GE. m*n + 3*n + bdspac ) THEN
952*
953* WORK( IU ) is M by N
954*
955 ldwrku = m
956 nwork = iu + ldwrku*n
957 CALL dlaset( 'F', m, n, zero, zero, work( iu ),
958 $ ldwrku )
959* IR is unused; silence compile warnings
960 ir = -1
961 ELSE
962*
963* WORK( IU ) is N by N
964*
965 ldwrku = n
966 nwork = iu + ldwrku*n
967*
968* WORK(IR) is LDWRKR by N
969*
970 ir = nwork
971 ldwrkr = ( lwork - n*n - 3*n ) / n
972 END IF
973 nwork = iu + ldwrku*n
974*
975* Perform bidiagonal SVD, computing left singular vectors
976* of bidiagonal matrix in WORK(IU) and computing right
977* singular vectors of bidiagonal matrix in VT
978* Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC
979*
980 CALL dbdsdc( 'U', 'I', n, s, work( ie ), work( iu ),
981 $ ldwrku, vt, ldvt, dum, idum, work( nwork ),
982 $ iwork, info )
983*
984* Overwrite VT by right singular vectors of A
985* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
986* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
987*
988 CALL dormbr( 'P', 'R', 'T', n, n, n, a, lda,
989 $ work( itaup ), vt, ldvt, work( nwork ),
990 $ lwork - nwork + 1, ierr )
991*
992 IF( lwork .GE. m*n + 3*n + bdspac ) THEN
993*
994* Path 5o-fast
995* Overwrite WORK(IU) by left singular vectors of A
996* Workspace: need 3*N [e, tauq, taup] + M*N [U] + N [work]
997* Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work]
998*
999 CALL dormbr( 'Q', 'L', 'N', m, n, n, a, lda,
1000 $ work( itauq ), work( iu ), ldwrku,
1001 $ work( nwork ), lwork - nwork + 1, ierr )
1002*
1003* Copy left singular vectors of A from WORK(IU) to A
1004*
1005 CALL dlacpy( 'F', m, n, work( iu ), ldwrku, a,
1006 $ lda )
1007 ELSE
1008*
1009* Path 5o-slow
1010* Generate Q in A
1011* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
1012* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
1013*
1014 CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
1015 $ work( nwork ), lwork - nwork + 1, ierr )
1016*
1017* Multiply Q in A by left singular vectors of
1018* bidiagonal matrix in WORK(IU), storing result in
1019* WORK(IR) and copying to A
1020* Workspace: need 3*N [e, tauq, taup] + N*N [U] + NB*N [R]
1021* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N [R]
1022*
1023 DO 20 i = 1, m, ldwrkr
1024 chunk = min( m - i + 1, ldwrkr )
1025 CALL dgemm( 'N', 'N', chunk, n, n, one, a( i,
1026 $ 1 ),
1027 $ lda, work( iu ), ldwrku, zero,
1028 $ work( ir ), ldwrkr )
1029 CALL dlacpy( 'F', chunk, n, work( ir ), ldwrkr,
1030 $ a( i, 1 ), lda )
1031 20 CONTINUE
1032 END IF
1033*
1034 ELSE IF( wntqs ) THEN
1035*
1036* Path 5s (M >= N, JOBZ='S')
1037* Perform bidiagonal SVD, computing left singular vectors
1038* of bidiagonal matrix in U and computing right singular
1039* vectors of bidiagonal matrix in VT
1040* Workspace: need 3*N [e, tauq, taup] + BDSPAC
1041*
1042 CALL dlaset( 'F', m, n, zero, zero, u, ldu )
1043 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
1044 $ ldvt, dum, idum, work( nwork ), iwork,
1045 $ info )
1046*
1047* Overwrite U by left singular vectors of A and VT
1048* by right singular vectors of A
1049* Workspace: need 3*N [e, tauq, taup] + N [work]
1050* Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1051*
1052 CALL dormbr( 'Q', 'L', 'N', m, n, n, a, lda,
1053 $ work( itauq ), u, ldu, work( nwork ),
1054 $ lwork - nwork + 1, ierr )
1055 CALL dormbr( 'P', 'R', 'T', n, n, n, a, lda,
1056 $ work( itaup ), vt, ldvt, work( nwork ),
1057 $ lwork - nwork + 1, ierr )
1058 ELSE IF( wntqa ) THEN
1059*
1060* Path 5a (M >= N, JOBZ='A')
1061* Perform bidiagonal SVD, computing left singular vectors
1062* of bidiagonal matrix in U and computing right singular
1063* vectors of bidiagonal matrix in VT
1064* Workspace: need 3*N [e, tauq, taup] + BDSPAC
1065*
1066 CALL dlaset( 'F', m, m, zero, zero, u, ldu )
1067 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
1068 $ ldvt, dum, idum, work( nwork ), iwork,
1069 $ info )
1070*
1071* Set the right corner of U to identity matrix
1072*
1073 IF( m.GT.n ) THEN
1074 CALL dlaset( 'F', m - n, m - n, zero, one, u(n+1,
1075 $ n+1),
1076 $ ldu )
1077 END IF
1078*
1079* Overwrite U by left singular vectors of A and VT
1080* by right singular vectors of A
1081* Workspace: need 3*N [e, tauq, taup] + M [work]
1082* Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1083*
1084 CALL dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1085 $ work( itauq ), u, ldu, work( nwork ),
1086 $ lwork - nwork + 1, ierr )
1087 CALL dormbr( 'P', 'R', 'T', n, n, m, a, lda,
1088 $ work( itaup ), vt, ldvt, work( nwork ),
1089 $ lwork - nwork + 1, ierr )
1090 END IF
1091*
1092 END IF
1093*
1094 ELSE
1095*
1096* A has more columns than rows. If A has sufficiently more
1097* columns than rows, first reduce using the LQ decomposition (if
1098* sufficient workspace available)
1099*
1100 IF( n.GE.mnthr ) THEN
1101*
1102 IF( wntqn ) THEN
1103*
1104* Path 1t (N >> M, JOBZ='N')
1105* No singular vectors to be computed
1106*
1107 itau = 1
1108 nwork = itau + m
1109*
1110* Compute A=L*Q
1111* Workspace: need M [tau] + M [work]
1112* Workspace: prefer M [tau] + M*NB [work]
1113*
1114 CALL dgelqf( m, n, a, lda, work( itau ),
1115 $ work( nwork ),
1116 $ lwork - nwork + 1, ierr )
1117*
1118* Zero out above L
1119*
1120 CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
1121 $ lda )
1122 ie = 1
1123 itauq = ie + m
1124 itaup = itauq + m
1125 nwork = itaup + m
1126*
1127* Bidiagonalize L in A
1128* Workspace: need 3*M [e, tauq, taup] + M [work]
1129* Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1130*
1131 CALL dgebrd( m, m, a, lda, s, work( ie ),
1132 $ work( itauq ),
1133 $ work( itaup ), work( nwork ), lwork-nwork+1,
1134 $ ierr )
1135 nwork = ie + m
1136*
1137* Perform bidiagonal SVD, computing singular values only
1138* Workspace: need M [e] + BDSPAC
1139*
1140 CALL dbdsdc( 'U', 'N', m, s, work( ie ), dum, 1, dum,
1141 $ 1,
1142 $ dum, idum, work( nwork ), iwork, info )
1143*
1144 ELSE IF( wntqo ) THEN
1145*
1146* Path 2t (N >> M, JOBZ='O')
1147* M right singular vectors to be overwritten on A and
1148* M left singular vectors to be computed in U
1149*
1150 ivt = 1
1151*
1152* WORK(IVT) is M by M
1153* WORK(IL) is M by M; it is later resized to M by chunk for gemm
1154*
1155 il = ivt + m*m
1156 IF( lwork .GE. m*n + m*m + 3*m + bdspac ) THEN
1157 ldwrkl = m
1158 chunk = n
1159 ELSE
1160 ldwrkl = m
1161 chunk = ( lwork - m*m ) / m
1162 END IF
1163 itau = il + ldwrkl*m
1164 nwork = itau + m
1165*
1166* Compute A=L*Q
1167* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1168* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1169*
1170 CALL dgelqf( m, n, a, lda, work( itau ),
1171 $ work( nwork ),
1172 $ lwork - nwork + 1, ierr )
1173*
1174* Copy L to WORK(IL), zeroing about above it
1175*
1176 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1177 CALL dlaset( 'U', m - 1, m - 1, zero, zero,
1178 $ work( il + ldwrkl ), ldwrkl )
1179*
1180* Generate Q in A
1181* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1182* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1183*
1184 CALL dorglq( m, n, m, a, lda, work( itau ),
1185 $ work( nwork ), lwork - nwork + 1, ierr )
1186 ie = itau
1187 itauq = ie + m
1188 itaup = itauq + m
1189 nwork = itaup + m
1190*
1191* Bidiagonalize L in WORK(IL)
1192* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1193* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1194*
1195 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1196 $ work( itauq ), work( itaup ), work( nwork ),
1197 $ lwork - nwork + 1, ierr )
1198*
1199* Perform bidiagonal SVD, computing left singular vectors
1200* of bidiagonal matrix in U, and computing right singular
1201* vectors of bidiagonal matrix in WORK(IVT)
1202* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1203*
1204 CALL dbdsdc( 'U', 'I', m, s, work( ie ), u, ldu,
1205 $ work( ivt ), m, dum, idum, work( nwork ),
1206 $ iwork, info )
1207*
1208* Overwrite U by left singular vectors of L and WORK(IVT)
1209* by right singular vectors of L
1210* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1211* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1212*
1213 CALL dormbr( 'Q', 'L', 'N', m, m, m, work( il ),
1214 $ ldwrkl,
1215 $ work( itauq ), u, ldu, work( nwork ),
1216 $ lwork - nwork + 1, ierr )
1217 CALL dormbr( 'P', 'R', 'T', m, m, m, work( il ),
1218 $ ldwrkl,
1219 $ work( itaup ), work( ivt ), m,
1220 $ work( nwork ), lwork - nwork + 1, ierr )
1221*
1222* Multiply right singular vectors of L in WORK(IVT) by Q
1223* in A, storing result in WORK(IL) and copying to A
1224* Workspace: need M*M [VT] + M*M [L]
1225* Workspace: prefer M*M [VT] + M*N [L]
1226* At this point, L is resized as M by chunk.
1227*
1228 DO 30 i = 1, n, chunk
1229 blk = min( n - i + 1, chunk )
1230 CALL dgemm( 'N', 'N', m, blk, m, one, work( ivt ),
1231 $ m,
1232 $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1233 CALL dlacpy( 'F', m, blk, work( il ), ldwrkl,
1234 $ a( 1, i ), lda )
1235 30 CONTINUE
1236*
1237 ELSE IF( wntqs ) THEN
1238*
1239* Path 3t (N >> M, JOBZ='S')
1240* M right singular vectors to be computed in VT and
1241* M left singular vectors to be computed in U
1242*
1243 il = 1
1244*
1245* WORK(IL) is M by M
1246*
1247 ldwrkl = m
1248 itau = il + ldwrkl*m
1249 nwork = itau + m
1250*
1251* Compute A=L*Q
1252* Workspace: need M*M [L] + M [tau] + M [work]
1253* Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1254*
1255 CALL dgelqf( m, n, a, lda, work( itau ),
1256 $ work( nwork ),
1257 $ lwork - nwork + 1, ierr )
1258*
1259* Copy L to WORK(IL), zeroing out above it
1260*
1261 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1262 CALL dlaset( 'U', m - 1, m - 1, zero, zero,
1263 $ work( il + ldwrkl ), ldwrkl )
1264*
1265* Generate Q in A
1266* Workspace: need M*M [L] + M [tau] + M [work]
1267* Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1268*
1269 CALL dorglq( m, n, m, a, lda, work( itau ),
1270 $ work( nwork ), lwork - nwork + 1, ierr )
1271 ie = itau
1272 itauq = ie + m
1273 itaup = itauq + m
1274 nwork = itaup + m
1275*
1276* Bidiagonalize L in WORK(IU).
1277* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1278* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1279*
1280 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1281 $ work( itauq ), work( itaup ), work( nwork ),
1282 $ lwork - nwork + 1, ierr )
1283*
1284* Perform bidiagonal SVD, computing left singular vectors
1285* of bidiagonal matrix in U and computing right singular
1286* vectors of bidiagonal matrix in VT
1287* Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1288*
1289 CALL dbdsdc( 'U', 'I', m, s, work( ie ), u, ldu, vt,
1290 $ ldvt, dum, idum, work( nwork ), iwork,
1291 $ info )
1292*
1293* Overwrite U by left singular vectors of L and VT
1294* by right singular vectors of L
1295* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1296* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1297*
1298 CALL dormbr( 'Q', 'L', 'N', m, m, m, work( il ),
1299 $ ldwrkl,
1300 $ work( itauq ), u, ldu, work( nwork ),
1301 $ lwork - nwork + 1, ierr )
1302 CALL dormbr( 'P', 'R', 'T', m, m, m, work( il ),
1303 $ ldwrkl,
1304 $ work( itaup ), vt, ldvt, work( nwork ),
1305 $ lwork - nwork + 1, ierr )
1306*
1307* Multiply right singular vectors of L in WORK(IL) by
1308* Q in A, storing result in VT
1309* Workspace: need M*M [L]
1310*
1311 CALL dlacpy( 'F', m, m, vt, ldvt, work( il ), ldwrkl )
1312 CALL dgemm( 'N', 'N', m, n, m, one, work( il ),
1313 $ ldwrkl,
1314 $ a, lda, zero, vt, ldvt )
1315*
1316 ELSE IF( wntqa ) THEN
1317*
1318* Path 4t (N >> M, JOBZ='A')
1319* N right singular vectors to be computed in VT and
1320* M left singular vectors to be computed in U
1321*
1322 ivt = 1
1323*
1324* WORK(IVT) is M by M
1325*
1326 ldwkvt = m
1327 itau = ivt + ldwkvt*m
1328 nwork = itau + m
1329*
1330* Compute A=L*Q, copying result to VT
1331* Workspace: need M*M [VT] + M [tau] + M [work]
1332* Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1333*
1334 CALL dgelqf( m, n, a, lda, work( itau ),
1335 $ work( nwork ),
1336 $ lwork - nwork + 1, ierr )
1337 CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
1338*
1339* Generate Q in VT
1340* Workspace: need M*M [VT] + M [tau] + N [work]
1341* Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1342*
1343 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
1344 $ work( nwork ), lwork - nwork + 1, ierr )
1345*
1346* Produce L in A, zeroing out other entries
1347*
1348 CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
1349 $ lda )
1350 ie = itau
1351 itauq = ie + m
1352 itaup = itauq + m
1353 nwork = itaup + m
1354*
1355* Bidiagonalize L in A
1356* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + M [work]
1357* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work]
1358*
1359 CALL dgebrd( m, m, a, lda, s, work( ie ),
1360 $ work( itauq ),
1361 $ work( itaup ), work( nwork ), lwork-nwork+1,
1362 $ ierr )
1363*
1364* Perform bidiagonal SVD, computing left singular vectors
1365* of bidiagonal matrix in U and computing right singular
1366* vectors of bidiagonal matrix in WORK(IVT)
1367* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1368*
1369 CALL dbdsdc( 'U', 'I', m, s, work( ie ), u, ldu,
1370 $ work( ivt ), ldwkvt, dum, idum,
1371 $ work( nwork ), iwork, info )
1372*
1373* Overwrite U by left singular vectors of L and WORK(IVT)
1374* by right singular vectors of L
1375* Workspace: need M*M [VT] + 3*M [e, tauq, taup]+ M [work]
1376* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work]
1377*
1378 CALL dormbr( 'Q', 'L', 'N', m, m, m, a, lda,
1379 $ work( itauq ), u, ldu, work( nwork ),
1380 $ lwork - nwork + 1, ierr )
1381 CALL dormbr( 'P', 'R', 'T', m, m, m, a, lda,
1382 $ work( itaup ), work( ivt ), ldwkvt,
1383 $ work( nwork ), lwork - nwork + 1, ierr )
1384*
1385* Multiply right singular vectors of L in WORK(IVT) by
1386* Q in VT, storing result in A
1387* Workspace: need M*M [VT]
1388*
1389 CALL dgemm( 'N', 'N', m, n, m, one, work( ivt ),
1390 $ ldwkvt,
1391 $ vt, ldvt, zero, a, lda )
1392*
1393* Copy right singular vectors of A from A to VT
1394*
1395 CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
1396*
1397 END IF
1398*
1399 ELSE
1400*
1401* N .LT. MNTHR
1402*
1403* Path 5t (N > M, but not much larger)
1404* Reduce to bidiagonal form without LQ decomposition
1405*
1406 ie = 1
1407 itauq = ie + m
1408 itaup = itauq + m
1409 nwork = itaup + m
1410*
1411* Bidiagonalize A
1412* Workspace: need 3*M [e, tauq, taup] + N [work]
1413* Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1414*
1415 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1416 $ work( itaup ), work( nwork ), lwork-nwork+1,
1417 $ ierr )
1418 IF( wntqn ) THEN
1419*
1420* Path 5tn (N > M, JOBZ='N')
1421* Perform bidiagonal SVD, only computing singular values
1422* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1423*
1424 CALL dbdsdc( 'L', 'N', m, s, work( ie ), dum, 1, dum,
1425 $ 1,
1426 $ dum, idum, work( nwork ), iwork, info )
1427 ELSE IF( wntqo ) THEN
1428* Path 5to (N > M, JOBZ='O')
1429 ldwkvt = m
1430 ivt = nwork
1431 IF( lwork .GE. m*n + 3*m + bdspac ) THEN
1432*
1433* WORK( IVT ) is M by N
1434*
1435 CALL dlaset( 'F', m, n, zero, zero, work( ivt ),
1436 $ ldwkvt )
1437 nwork = ivt + ldwkvt*n
1438* IL is unused; silence compile warnings
1439 il = -1
1440 ELSE
1441*
1442* WORK( IVT ) is M by M
1443*
1444 nwork = ivt + ldwkvt*m
1445 il = nwork
1446*
1447* WORK(IL) is M by CHUNK
1448*
1449 chunk = ( lwork - m*m - 3*m ) / m
1450 END IF
1451*
1452* Perform bidiagonal SVD, computing left singular vectors
1453* of bidiagonal matrix in U and computing right singular
1454* vectors of bidiagonal matrix in WORK(IVT)
1455* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1456*
1457 CALL dbdsdc( 'L', 'I', m, s, work( ie ), u, ldu,
1458 $ work( ivt ), ldwkvt, dum, idum,
1459 $ work( nwork ), iwork, info )
1460*
1461* Overwrite U by left singular vectors of A
1462* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1463* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1464*
1465 CALL dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1466 $ work( itauq ), u, ldu, work( nwork ),
1467 $ lwork - nwork + 1, ierr )
1468*
1469 IF( lwork .GE. m*n + 3*m + bdspac ) THEN
1470*
1471* Path 5to-fast
1472* Overwrite WORK(IVT) by left singular vectors of A
1473* Workspace: need 3*M [e, tauq, taup] + M*N [VT] + M [work]
1474* Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work]
1475*
1476 CALL dormbr( 'P', 'R', 'T', m, n, m, a, lda,
1477 $ work( itaup ), work( ivt ), ldwkvt,
1478 $ work( nwork ), lwork - nwork + 1, ierr )
1479*
1480* Copy right singular vectors of A from WORK(IVT) to A
1481*
1482 CALL dlacpy( 'F', m, n, work( ivt ), ldwkvt, a,
1483 $ lda )
1484 ELSE
1485*
1486* Path 5to-slow
1487* Generate P**T in A
1488* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1489* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1490*
1491 CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
1492 $ work( nwork ), lwork - nwork + 1, ierr )
1493*
1494* Multiply Q in A by right singular vectors of
1495* bidiagonal matrix in WORK(IVT), storing result in
1496* WORK(IL) and copying to A
1497* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M*NB [L]
1498* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N [L]
1499*
1500 DO 40 i = 1, n, chunk
1501 blk = min( n - i + 1, chunk )
1502 CALL dgemm( 'N', 'N', m, blk, m, one,
1503 $ work( ivt ),
1504 $ ldwkvt, a( 1, i ), lda, zero,
1505 $ work( il ), m )
1506 CALL dlacpy( 'F', m, blk, work( il ), m, a( 1,
1507 $ i ),
1508 $ lda )
1509 40 CONTINUE
1510 END IF
1511 ELSE IF( wntqs ) THEN
1512*
1513* Path 5ts (N > M, JOBZ='S')
1514* Perform bidiagonal SVD, computing left singular vectors
1515* of bidiagonal matrix in U and computing right singular
1516* vectors of bidiagonal matrix in VT
1517* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1518*
1519 CALL dlaset( 'F', m, n, zero, zero, vt, ldvt )
1520 CALL dbdsdc( 'L', 'I', m, s, work( ie ), u, ldu, vt,
1521 $ ldvt, dum, idum, work( nwork ), iwork,
1522 $ info )
1523*
1524* Overwrite U by left singular vectors of A and VT
1525* by right singular vectors of A
1526* Workspace: need 3*M [e, tauq, taup] + M [work]
1527* Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1528*
1529 CALL dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1530 $ work( itauq ), u, ldu, work( nwork ),
1531 $ lwork - nwork + 1, ierr )
1532 CALL dormbr( 'P', 'R', 'T', m, n, m, a, lda,
1533 $ work( itaup ), vt, ldvt, work( nwork ),
1534 $ lwork - nwork + 1, ierr )
1535 ELSE IF( wntqa ) THEN
1536*
1537* Path 5ta (N > M, JOBZ='A')
1538* Perform bidiagonal SVD, computing left singular vectors
1539* of bidiagonal matrix in U and computing right singular
1540* vectors of bidiagonal matrix in VT
1541* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1542*
1543 CALL dlaset( 'F', n, n, zero, zero, vt, ldvt )
1544 CALL dbdsdc( 'L', 'I', m, s, work( ie ), u, ldu, vt,
1545 $ ldvt, dum, idum, work( nwork ), iwork,
1546 $ info )
1547*
1548* Set the right corner of VT to identity matrix
1549*
1550 IF( n.GT.m ) THEN
1551 CALL dlaset( 'F', n-m, n-m, zero, one, vt(m+1,m+1),
1552 $ ldvt )
1553 END IF
1554*
1555* Overwrite U by left singular vectors of A and VT
1556* by right singular vectors of A
1557* Workspace: need 3*M [e, tauq, taup] + N [work]
1558* Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1559*
1560 CALL dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1561 $ work( itauq ), u, ldu, work( nwork ),
1562 $ lwork - nwork + 1, ierr )
1563 CALL dormbr( 'P', 'R', 'T', n, n, m, a, lda,
1564 $ work( itaup ), vt, ldvt, work( nwork ),
1565 $ lwork - nwork + 1, ierr )
1566 END IF
1567*
1568 END IF
1569*
1570 END IF
1571*
1572* Undo scaling if necessary
1573*
1574 IF( iscl.EQ.1 ) THEN
1575 IF( anrm.GT.bignum )
1576 $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1577 $ ierr )
1578 IF( anrm.LT.smlnum )
1579 $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
1580 $ ierr )
1581 END IF
1582*
1583* Return optimal workspace in WORK(1)
1584*
1585 work( 1 ) = droundup_lwork( maxwrk )
1586*
1587 RETURN
1588*
1589* End of DGESDD
1590*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dbdsdc
subroutine dbdsdc(uplo, compq, n, d, e, u, ldu, vt, ldvt, q, iq, work, iwork, info)
DBDSDC
Definition dbdsdc.f:197
dgebrd
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:204
dgelqf
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:142
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgeqrf
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:144
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:57
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:112
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:142
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
droundup_lwork
double precision function droundup_lwork(lwork)
DROUNDUP_LWORK
Definition droundup_lwork.f:59
dorgbr
subroutine dorgbr(vect, m, n, k, a, lda, tau, work, lwork, info)
DORGBR
Definition dorgbr.f:156
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:125
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:126
dormbr
subroutine dormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMBR
Definition dormbr.f:193
Here is the call graph for this function:
Here is the caller graph for this function: