LAPACK 3.12.0
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 211 of file dgesdd.f.

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