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

◆ zgesdd()

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

ZGESDD

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

Purpose:
 ZGESDD computes the singular value decomposition (SVD) of a complex
 M-by-N matrix A, optionally computing the left and/or right singular
 vectors, by using divide-and-conquer method. The SVD is written

      A = U * SIGMA * conjugate-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 unitary matrix, and
 V is an N-by-N unitary 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**H, 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**H 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**H are returned in the arrays U
                  and VT;
          = 'O':  If M >= N, the first N columns of U are overwritten
                  in the array A and all rows of V**H are returned in
                  the array VT;
                  otherwise, all columns of U are returned in the
                  array U and the first M rows of V**H are overwritten
                  in the array A;
          = 'N':  no columns of U or rows of V**H 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 COMPLEX*16 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**H (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 COMPLEX*16 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
          unitary 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 COMPLEX*16 array, dimension (LDVT,N)
          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
          N-by-N unitary matrix V**H;
          if JOBZ = 'S', VT contains the first min(M,N) rows of
          V**H (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 COMPLEX*16 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 >= 2*mn + mx.
          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn.
          If JOBZ = 'A', LWORK >=   mn*mn + 2*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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
          Let mx = max(M,N) and mn = min(M,N).
          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
          else              LRWORK >= max( 5*mn*mn + 5*mn,
                                           2*mx*mn + 2*mn*mn + mn ).
[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:  The updating process of DBDSDC did not converge.
          =  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 219 of file zgesdd.f.

221 implicit none
222*
223* -- LAPACK driver routine --
224* -- LAPACK is a software package provided by Univ. of Tennessee, --
225* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226*
227* .. Scalar Arguments ..
228 CHARACTER JOBZ
229 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
230* ..
231* .. Array Arguments ..
232 INTEGER IWORK( * )
233 DOUBLE PRECISION RWORK( * ), S( * )
234 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
235 $ WORK( * )
236* ..
237*
238* =====================================================================
239*
240* .. Parameters ..
241 COMPLEX*16 CZERO, CONE
242 parameter( czero = ( 0.0d+0, 0.0d+0 ),
243 $ cone = ( 1.0d+0, 0.0d+0 ) )
244 DOUBLE PRECISION ZERO, ONE
245 parameter( zero = 0.0d+0, one = 1.0d+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
249 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
250 $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
251 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
252 $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
253 INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
254 $ LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN,
255 $ LWORK_ZGEQRF_MN,
256 $ LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN,
257 $ LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM,
258 $ LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN,
259 $ LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN,
260 $ LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM,
261 $ LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN,
262 $ LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN
263 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
264* ..
265* .. Local Arrays ..
266 INTEGER IDUM( 1 )
267 DOUBLE PRECISION DUM( 1 )
268 COMPLEX*16 CDUM( 1 )
269* ..
270* .. External Subroutines ..
271 EXTERNAL dbdsdc, dlascl, xerbla, zgebrd, zgelqf, zgemm,
272 $ zgeqrf, zlacp2, zlacpy, zlacrm, zlarcm, zlascl,
273 $ zlaset, zungbr, zunglq, zungqr, zunmbr
274* ..
275* .. External Functions ..
276 LOGICAL LSAME, DISNAN
277 DOUBLE PRECISION DLAMCH, ZLANGE, DROUNDUP_LWORK
278 EXTERNAL lsame, dlamch, zlange, disnan,
279 $ droundup_lwork
280* ..
281* .. Intrinsic Functions ..
282 INTRINSIC int, max, min, sqrt
283* ..
284* .. Executable Statements ..
285*
286* Test the input arguments
287*
288 info = 0
289 minmn = min( m, n )
290 mnthr1 = int( minmn*17.0d0 / 9.0d0 )
291 mnthr2 = int( minmn*5.0d0 / 3.0d0 )
292 wntqa = lsame( jobz, 'A' )
293 wntqs = lsame( jobz, 'S' )
294 wntqas = wntqa .OR. wntqs
295 wntqo = lsame( jobz, 'O' )
296 wntqn = lsame( jobz, 'N' )
297 lquery = ( lwork.EQ.-1 )
298 minwrk = 1
299 maxwrk = 1
300*
301 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) ) THEN
302 info = -1
303 ELSE IF( m.LT.0 ) THEN
304 info = -2
305 ELSE IF( n.LT.0 ) THEN
306 info = -3
307 ELSE IF( lda.LT.max( 1, m ) ) THEN
308 info = -5
309 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
310 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) ) THEN
311 info = -8
312 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
313 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
314 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) ) THEN
315 info = -10
316 END IF
317*
318* Compute workspace
319* Note: Comments in the code beginning "Workspace:" describe the
320* minimal amount of workspace allocated at that point in the code,
321* as well as the preferred amount for good performance.
322* CWorkspace refers to complex workspace, and RWorkspace to
323* real workspace. NB refers to the optimal block size for the
324* immediately following subroutine, as returned by ILAENV.)
325*
326 IF( info.EQ.0 ) THEN
327 minwrk = 1
328 maxwrk = 1
329 IF( m.GE.n .AND. minmn.GT.0 ) THEN
330*
331* There is no complex work space needed for bidiagonal SVD
332* The real work space needed for bidiagonal SVD (dbdsdc) is
333* BDSPAC = 3*N*N + 4*N for singular values and vectors;
334* BDSPAC = 4*N for singular values only;
335* not including e, RU, and RVT matrices.
336*
337* Compute space preferred for each routine
338 CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
339 $ cdum(1), cdum(1), -1, ierr )
340 lwork_zgebrd_mn = int( cdum(1) )
341*
342 CALL zgebrd( n, n, cdum(1), n, dum(1), dum(1), cdum(1),
343 $ cdum(1), cdum(1), -1, ierr )
344 lwork_zgebrd_nn = int( cdum(1) )
345*
346 CALL zgeqrf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
347 lwork_zgeqrf_mn = int( cdum(1) )
348*
349 CALL zungbr( 'P', n, n, n, cdum(1), n, cdum(1), cdum(1),
350 $ -1, ierr )
351 lwork_zungbr_p_nn = int( cdum(1) )
352*
353 CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
354 $ -1, ierr )
355 lwork_zungbr_q_mm = int( cdum(1) )
356*
357 CALL zungbr( 'Q', m, n, n, cdum(1), m, cdum(1), cdum(1),
358 $ -1, ierr )
359 lwork_zungbr_q_mn = int( cdum(1) )
360*
361 CALL zungqr( m, m, n, cdum(1), m, cdum(1), cdum(1),
362 $ -1, ierr )
363 lwork_zungqr_mm = int( cdum(1) )
364*
365 CALL zungqr( m, n, n, cdum(1), m, cdum(1), cdum(1),
366 $ -1, ierr )
367 lwork_zungqr_mn = int( cdum(1) )
368*
369 CALL zunmbr( 'P', 'R', 'C', n, n, n, cdum(1), n, cdum(1),
370 $ cdum(1), n, cdum(1), -1, ierr )
371 lwork_zunmbr_prc_nn = int( cdum(1) )
372*
373 CALL zunmbr( 'Q', 'L', 'N', m, m, n, cdum(1), m, cdum(1),
374 $ cdum(1), m, cdum(1), -1, ierr )
375 lwork_zunmbr_qln_mm = int( cdum(1) )
376*
377 CALL zunmbr( 'Q', 'L', 'N', m, n, n, cdum(1), m, cdum(1),
378 $ cdum(1), m, cdum(1), -1, ierr )
379 lwork_zunmbr_qln_mn = int( cdum(1) )
380*
381 CALL zunmbr( 'Q', 'L', 'N', n, n, n, cdum(1), n, cdum(1),
382 $ cdum(1), n, cdum(1), -1, ierr )
383 lwork_zunmbr_qln_nn = int( cdum(1) )
384*
385 IF( m.GE.mnthr1 ) THEN
386 IF( wntqn ) THEN
387*
388* Path 1 (M >> N, JOBZ='N')
389*
390 maxwrk = n + lwork_zgeqrf_mn
391 maxwrk = max( maxwrk, 2*n + lwork_zgebrd_nn )
392 minwrk = 3*n
393 ELSE IF( wntqo ) THEN
394*
395* Path 2 (M >> N, JOBZ='O')
396*
397 wrkbl = n + lwork_zgeqrf_mn
398 wrkbl = max( wrkbl, n + lwork_zungqr_mn )
399 wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
400 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
401 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
402 maxwrk = m*n + n*n + wrkbl
403 minwrk = 2*n*n + 3*n
404 ELSE IF( wntqs ) THEN
405*
406* Path 3 (M >> N, JOBZ='S')
407*
408 wrkbl = n + lwork_zgeqrf_mn
409 wrkbl = max( wrkbl, n + lwork_zungqr_mn )
410 wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
411 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
412 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
413 maxwrk = n*n + wrkbl
414 minwrk = n*n + 3*n
415 ELSE IF( wntqa ) THEN
416*
417* Path 4 (M >> N, JOBZ='A')
418*
419 wrkbl = n + lwork_zgeqrf_mn
420 wrkbl = max( wrkbl, n + lwork_zungqr_mm )
421 wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
422 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
423 wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
424 maxwrk = n*n + wrkbl
425 minwrk = n*n + max( 3*n, n + m )
426 END IF
427 ELSE IF( m.GE.mnthr2 ) THEN
428*
429* Path 5 (M >> N, but not as much as MNTHR1)
430*
431 maxwrk = 2*n + lwork_zgebrd_mn
432 minwrk = 2*n + m
433 IF( wntqo ) THEN
434* Path 5o (M >> N, JOBZ='O')
435 maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
436 maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
437 maxwrk = maxwrk + m*n
438 minwrk = minwrk + n*n
439 ELSE IF( wntqs ) THEN
440* Path 5s (M >> N, JOBZ='S')
441 maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
442 maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
443 ELSE IF( wntqa ) THEN
444* Path 5a (M >> N, JOBZ='A')
445 maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
446 maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mm )
447 END IF
448 ELSE
449*
450* Path 6 (M >= N, but not much larger)
451*
452 maxwrk = 2*n + lwork_zgebrd_mn
453 minwrk = 2*n + m
454 IF( wntqo ) THEN
455* Path 6o (M >= N, JOBZ='O')
456 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
457 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
458 maxwrk = maxwrk + m*n
459 minwrk = minwrk + n*n
460 ELSE IF( wntqs ) THEN
461* Path 6s (M >= N, JOBZ='S')
462 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
463 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
464 ELSE IF( wntqa ) THEN
465* Path 6a (M >= N, JOBZ='A')
466 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mm )
467 maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
468 END IF
469 END IF
470 ELSE IF( minmn.GT.0 ) THEN
471*
472* There is no complex work space needed for bidiagonal SVD
473* The real work space needed for bidiagonal SVD (dbdsdc) is
474* BDSPAC = 3*M*M + 4*M for singular values and vectors;
475* BDSPAC = 4*M for singular values only;
476* not including e, RU, and RVT matrices.
477*
478* Compute space preferred for each routine
479 CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
480 $ cdum(1), cdum(1), -1, ierr )
481 lwork_zgebrd_mn = int( cdum(1) )
482*
483 CALL zgebrd( m, m, cdum(1), m, dum(1), dum(1), cdum(1),
484 $ cdum(1), cdum(1), -1, ierr )
485 lwork_zgebrd_mm = int( cdum(1) )
486*
487 CALL zgelqf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
488 lwork_zgelqf_mn = int( cdum(1) )
489*
490 CALL zungbr( 'P', m, n, m, cdum(1), m, cdum(1), cdum(1),
491 $ -1, ierr )
492 lwork_zungbr_p_mn = int( cdum(1) )
493*
494 CALL zungbr( 'P', n, n, m, cdum(1), n, cdum(1), cdum(1),
495 $ -1, ierr )
496 lwork_zungbr_p_nn = int( cdum(1) )
497*
498 CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
499 $ -1, ierr )
500 lwork_zungbr_q_mm = int( cdum(1) )
501*
502 CALL zunglq( m, n, m, cdum(1), m, cdum(1), cdum(1),
503 $ -1, ierr )
504 lwork_zunglq_mn = int( cdum(1) )
505*
506 CALL zunglq( n, n, m, cdum(1), n, cdum(1), cdum(1),
507 $ -1, ierr )
508 lwork_zunglq_nn = int( cdum(1) )
509*
510 CALL zunmbr( 'P', 'R', 'C', m, m, m, cdum(1), m, cdum(1),
511 $ cdum(1), m, cdum(1), -1, ierr )
512 lwork_zunmbr_prc_mm = int( cdum(1) )
513*
514 CALL zunmbr( 'P', 'R', 'C', m, n, m, cdum(1), m, cdum(1),
515 $ cdum(1), m, cdum(1), -1, ierr )
516 lwork_zunmbr_prc_mn = int( cdum(1) )
517*
518 CALL zunmbr( 'P', 'R', 'C', n, n, m, cdum(1), n, cdum(1),
519 $ cdum(1), n, cdum(1), -1, ierr )
520 lwork_zunmbr_prc_nn = int( cdum(1) )
521*
522 CALL zunmbr( 'Q', 'L', 'N', m, m, m, cdum(1), m, cdum(1),
523 $ cdum(1), m, cdum(1), -1, ierr )
524 lwork_zunmbr_qln_mm = int( cdum(1) )
525*
526 IF( n.GE.mnthr1 ) THEN
527 IF( wntqn ) THEN
528*
529* Path 1t (N >> M, JOBZ='N')
530*
531 maxwrk = m + lwork_zgelqf_mn
532 maxwrk = max( maxwrk, 2*m + lwork_zgebrd_mm )
533 minwrk = 3*m
534 ELSE IF( wntqo ) THEN
535*
536* Path 2t (N >> M, JOBZ='O')
537*
538 wrkbl = m + lwork_zgelqf_mn
539 wrkbl = max( wrkbl, m + lwork_zunglq_mn )
540 wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
541 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
542 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
543 maxwrk = m*n + m*m + wrkbl
544 minwrk = 2*m*m + 3*m
545 ELSE IF( wntqs ) THEN
546*
547* Path 3t (N >> M, JOBZ='S')
548*
549 wrkbl = m + lwork_zgelqf_mn
550 wrkbl = max( wrkbl, m + lwork_zunglq_mn )
551 wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
552 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
553 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
554 maxwrk = m*m + wrkbl
555 minwrk = m*m + 3*m
556 ELSE IF( wntqa ) THEN
557*
558* Path 4t (N >> M, JOBZ='A')
559*
560 wrkbl = m + lwork_zgelqf_mn
561 wrkbl = max( wrkbl, m + lwork_zunglq_nn )
562 wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
563 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
564 wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
565 maxwrk = m*m + wrkbl
566 minwrk = m*m + max( 3*m, m + n )
567 END IF
568 ELSE IF( n.GE.mnthr2 ) THEN
569*
570* Path 5t (N >> M, but not as much as MNTHR1)
571*
572 maxwrk = 2*m + lwork_zgebrd_mn
573 minwrk = 2*m + n
574 IF( wntqo ) THEN
575* Path 5to (N >> M, JOBZ='O')
576 maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
577 maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
578 maxwrk = maxwrk + m*n
579 minwrk = minwrk + m*m
580 ELSE IF( wntqs ) THEN
581* Path 5ts (N >> M, JOBZ='S')
582 maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
583 maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
584 ELSE IF( wntqa ) THEN
585* Path 5ta (N >> M, JOBZ='A')
586 maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
587 maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_nn )
588 END IF
589 ELSE
590*
591* Path 6t (N > M, but not much larger)
592*
593 maxwrk = 2*m + lwork_zgebrd_mn
594 minwrk = 2*m + n
595 IF( wntqo ) THEN
596* Path 6to (N > M, JOBZ='O')
597 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
598 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
599 maxwrk = maxwrk + m*n
600 minwrk = minwrk + m*m
601 ELSE IF( wntqs ) THEN
602* Path 6ts (N > M, JOBZ='S')
603 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
604 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
605 ELSE IF( wntqa ) THEN
606* Path 6ta (N > M, JOBZ='A')
607 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
608 maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_nn )
609 END IF
610 END IF
611 END IF
612 maxwrk = max( maxwrk, minwrk )
613 END IF
614 IF( info.EQ.0 ) THEN
615 work( 1 ) = droundup_lwork( maxwrk )
616 IF( lwork.LT.minwrk .AND. .NOT. lquery ) THEN
617 info = -12
618 END IF
619 END IF
620*
621 IF( info.NE.0 ) THEN
622 CALL xerbla( 'ZGESDD', -info )
623 RETURN
624 ELSE IF( lquery ) THEN
625 RETURN
626 END IF
627*
628* Quick return if possible
629*
630 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
631 RETURN
632 END IF
633*
634* Get machine constants
635*
636 eps = dlamch( 'P' )
637 smlnum = sqrt( dlamch( 'S' ) ) / eps
638 bignum = one / smlnum
639*
640* Scale A if max element outside range [SMLNUM,BIGNUM]
641*
642 anrm = zlange( 'M', m, n, a, lda, dum )
643 IF( disnan( anrm ) ) THEN
644 info = -4
645 RETURN
646 END IF
647 iscl = 0
648 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
649 iscl = 1
650 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
651 ELSE IF( anrm.GT.bignum ) THEN
652 iscl = 1
653 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
654 END IF
655*
656 IF( m.GE.n ) THEN
657*
658* A has at least as many rows as columns. If A has sufficiently
659* more rows than columns, first reduce using the QR
660* decomposition (if sufficient workspace available)
661*
662 IF( m.GE.mnthr1 ) THEN
663*
664 IF( wntqn ) THEN
665*
666* Path 1 (M >> N, JOBZ='N')
667* No singular vectors to be computed
668*
669 itau = 1
670 nwork = itau + n
671*
672* Compute A=Q*R
673* CWorkspace: need N [tau] + N [work]
674* CWorkspace: prefer N [tau] + N*NB [work]
675* RWorkspace: need 0
676*
677 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
678 $ lwork-nwork+1, ierr )
679*
680* Zero out below R
681*
682 CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
683 $ lda )
684 ie = 1
685 itauq = 1
686 itaup = itauq + n
687 nwork = itaup + n
688*
689* Bidiagonalize R in A
690* CWorkspace: need 2*N [tauq, taup] + N [work]
691* CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
692* RWorkspace: need N [e]
693*
694 CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
695 $ work( itaup ), work( nwork ), lwork-nwork+1,
696 $ ierr )
697 nrwork = ie + n
698*
699* Perform bidiagonal SVD, compute singular values only
700* CWorkspace: need 0
701* RWorkspace: need N [e] + BDSPAC
702*
703 CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
704 $ dum, idum, rwork( nrwork ), iwork, info )
705*
706 ELSE IF( wntqo ) THEN
707*
708* Path 2 (M >> N, JOBZ='O')
709* N left singular vectors to be overwritten on A and
710* N right singular vectors to be computed in VT
711*
712 iu = 1
713*
714* WORK(IU) is N by N
715*
716 ldwrku = n
717 ir = iu + ldwrku*n
718 IF( lwork .GE. m*n + n*n + 3*n ) THEN
719*
720* WORK(IR) is M by N
721*
722 ldwrkr = m
723 ELSE
724 ldwrkr = ( lwork - n*n - 3*n ) / n
725 END IF
726 itau = ir + ldwrkr*n
727 nwork = itau + n
728*
729* Compute A=Q*R
730* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
731* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
732* RWorkspace: need 0
733*
734 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
735 $ lwork-nwork+1, ierr )
736*
737* Copy R to WORK( IR ), zeroing out below it
738*
739 CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
740 CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
741 $ ldwrkr )
742*
743* Generate Q in A
744* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
745* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
746* RWorkspace: need 0
747*
748 CALL zungqr( m, n, n, a, lda, work( itau ),
749 $ work( nwork ), lwork-nwork+1, ierr )
750 ie = 1
751 itauq = itau
752 itaup = itauq + n
753 nwork = itaup + n
754*
755* Bidiagonalize R in WORK(IR)
756* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
757* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
758* RWorkspace: need N [e]
759*
760 CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
761 $ work( itauq ), work( itaup ), work( nwork ),
762 $ lwork-nwork+1, ierr )
763*
764* Perform bidiagonal SVD, computing left singular vectors
765* of R in WORK(IRU) and computing right singular vectors
766* of R in WORK(IRVT)
767* CWorkspace: need 0
768* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
769*
770 iru = ie + n
771 irvt = iru + n*n
772 nrwork = irvt + n*n
773 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
774 $ n, rwork( irvt ), n, dum, idum,
775 $ rwork( nrwork ), iwork, info )
776*
777* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
778* Overwrite WORK(IU) by the left singular vectors of R
779* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
780* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
781* RWorkspace: need 0
782*
783 CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
784 $ ldwrku )
785 CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
786 $ work( itauq ), work( iu ), ldwrku,
787 $ work( nwork ), lwork-nwork+1, ierr )
788*
789* Copy real matrix RWORK(IRVT) to complex matrix VT
790* Overwrite VT by the right singular vectors of R
791* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
792* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
793* RWorkspace: need 0
794*
795 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
796 CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
797 $ work( itaup ), vt, ldvt, work( nwork ),
798 $ lwork-nwork+1, ierr )
799*
800* Multiply Q in A by left singular vectors of R in
801* WORK(IU), storing result in WORK(IR) and copying to A
802* CWorkspace: need N*N [U] + N*N [R]
803* CWorkspace: prefer N*N [U] + M*N [R]
804* RWorkspace: need 0
805*
806 DO 10 i = 1, m, ldwrkr
807 chunk = min( m-i+1, ldwrkr )
808 CALL zgemm( 'N', 'N', chunk, n, n, cone, a( i, 1 ),
809 $ lda, work( iu ), ldwrku, czero,
810 $ work( ir ), ldwrkr )
811 CALL zlacpy( 'F', chunk, n, work( ir ), ldwrkr,
812 $ a( i, 1 ), lda )
813 10 CONTINUE
814*
815 ELSE IF( wntqs ) THEN
816*
817* Path 3 (M >> N, JOBZ='S')
818* N left singular vectors to be computed in U and
819* N right singular vectors to be computed in VT
820*
821 ir = 1
822*
823* WORK(IR) is N by N
824*
825 ldwrkr = n
826 itau = ir + ldwrkr*n
827 nwork = itau + n
828*
829* Compute A=Q*R
830* CWorkspace: need N*N [R] + N [tau] + N [work]
831* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
832* RWorkspace: need 0
833*
834 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
835 $ lwork-nwork+1, ierr )
836*
837* Copy R to WORK(IR), zeroing out below it
838*
839 CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
840 CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
841 $ ldwrkr )
842*
843* Generate Q in A
844* CWorkspace: need N*N [R] + N [tau] + N [work]
845* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
846* RWorkspace: need 0
847*
848 CALL zungqr( m, n, n, a, lda, work( itau ),
849 $ work( nwork ), lwork-nwork+1, ierr )
850 ie = 1
851 itauq = itau
852 itaup = itauq + n
853 nwork = itaup + n
854*
855* Bidiagonalize R in WORK(IR)
856* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
857* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
858* RWorkspace: need N [e]
859*
860 CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
861 $ work( itauq ), work( itaup ), work( nwork ),
862 $ lwork-nwork+1, ierr )
863*
864* Perform bidiagonal SVD, computing left singular vectors
865* of bidiagonal matrix in RWORK(IRU) and computing right
866* singular vectors of bidiagonal matrix in RWORK(IRVT)
867* CWorkspace: need 0
868* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
869*
870 iru = ie + n
871 irvt = iru + n*n
872 nrwork = irvt + n*n
873 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
874 $ n, rwork( irvt ), n, dum, idum,
875 $ rwork( nrwork ), iwork, info )
876*
877* Copy real matrix RWORK(IRU) to complex matrix U
878* Overwrite U by left singular vectors of R
879* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
880* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
881* RWorkspace: need 0
882*
883 CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
884 CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
885 $ work( itauq ), u, ldu, work( nwork ),
886 $ lwork-nwork+1, ierr )
887*
888* Copy real matrix RWORK(IRVT) to complex matrix VT
889* Overwrite VT by right singular vectors of R
890* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
891* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
892* RWorkspace: need 0
893*
894 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
895 CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
896 $ work( itaup ), vt, ldvt, work( nwork ),
897 $ lwork-nwork+1, ierr )
898*
899* Multiply Q in A by left singular vectors of R in
900* WORK(IR), storing result in U
901* CWorkspace: need N*N [R]
902* RWorkspace: need 0
903*
904 CALL zlacpy( 'F', n, n, u, ldu, work( ir ), ldwrkr )
905 CALL zgemm( 'N', 'N', m, n, n, cone, a, lda, work( ir ),
906 $ ldwrkr, czero, u, ldu )
907*
908 ELSE IF( wntqa ) THEN
909*
910* Path 4 (M >> N, JOBZ='A')
911* M left singular vectors to be computed in U and
912* N right singular vectors to be computed in VT
913*
914 iu = 1
915*
916* WORK(IU) is N by N
917*
918 ldwrku = n
919 itau = iu + ldwrku*n
920 nwork = itau + n
921*
922* Compute A=Q*R, copying result to U
923* CWorkspace: need N*N [U] + N [tau] + N [work]
924* CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
925* RWorkspace: need 0
926*
927 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
928 $ lwork-nwork+1, ierr )
929 CALL zlacpy( 'L', m, n, a, lda, u, ldu )
930*
931* Generate Q in U
932* CWorkspace: need N*N [U] + N [tau] + M [work]
933* CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
934* RWorkspace: need 0
935*
936 CALL zungqr( m, m, n, u, ldu, work( itau ),
937 $ work( nwork ), lwork-nwork+1, ierr )
938*
939* Produce R in A, zeroing out below it
940*
941 CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
942 $ lda )
943 ie = 1
944 itauq = itau
945 itaup = itauq + n
946 nwork = itaup + n
947*
948* Bidiagonalize R in A
949* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
950* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
951* RWorkspace: need N [e]
952*
953 CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
954 $ work( itaup ), work( nwork ), lwork-nwork+1,
955 $ ierr )
956 iru = ie + n
957 irvt = iru + n*n
958 nrwork = irvt + n*n
959*
960* Perform bidiagonal SVD, computing left singular vectors
961* of bidiagonal matrix in RWORK(IRU) and computing right
962* singular vectors of bidiagonal matrix in RWORK(IRVT)
963* CWorkspace: need 0
964* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
965*
966 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
967 $ n, rwork( irvt ), n, dum, idum,
968 $ rwork( nrwork ), iwork, info )
969*
970* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
971* Overwrite WORK(IU) by left singular vectors of R
972* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
973* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
974* RWorkspace: need 0
975*
976 CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
977 $ ldwrku )
978 CALL zunmbr( 'Q', 'L', 'N', n, n, n, a, lda,
979 $ work( itauq ), work( iu ), ldwrku,
980 $ work( nwork ), lwork-nwork+1, ierr )
981*
982* Copy real matrix RWORK(IRVT) to complex matrix VT
983* Overwrite VT by right singular vectors of R
984* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
985* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
986* RWorkspace: need 0
987*
988 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
989 CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
990 $ work( itaup ), vt, ldvt, work( nwork ),
991 $ lwork-nwork+1, ierr )
992*
993* Multiply Q in U by left singular vectors of R in
994* WORK(IU), storing result in A
995* CWorkspace: need N*N [U]
996* RWorkspace: need 0
997*
998 CALL zgemm( 'N', 'N', m, n, n, cone, u, ldu, work( iu ),
999 $ ldwrku, czero, a, lda )
1000*
1001* Copy left singular vectors of A from A to U
1002*
1003 CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1004*
1005 END IF
1006*
1007 ELSE IF( m.GE.mnthr2 ) THEN
1008*
1009* MNTHR2 <= M < MNTHR1
1010*
1011* Path 5 (M >> N, but not as much as MNTHR1)
1012* Reduce to bidiagonal form without QR decomposition, use
1013* ZUNGBR and matrix multiplication to compute singular vectors
1014*
1015 ie = 1
1016 nrwork = ie + n
1017 itauq = 1
1018 itaup = itauq + n
1019 nwork = itaup + n
1020*
1021* Bidiagonalize A
1022* CWorkspace: need 2*N [tauq, taup] + M [work]
1023* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1024* RWorkspace: need N [e]
1025*
1026 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1027 $ work( itaup ), work( nwork ), lwork-nwork+1,
1028 $ ierr )
1029 IF( wntqn ) THEN
1030*
1031* Path 5n (M >> N, JOBZ='N')
1032* Compute singular values only
1033* CWorkspace: need 0
1034* RWorkspace: need N [e] + BDSPAC
1035*
1036 CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum, 1,dum,1,
1037 $ dum, idum, rwork( nrwork ), iwork, info )
1038 ELSE IF( wntqo ) THEN
1039 iu = nwork
1040 iru = nrwork
1041 irvt = iru + n*n
1042 nrwork = irvt + n*n
1043*
1044* Path 5o (M >> N, JOBZ='O')
1045* Copy A to VT, generate P**H
1046* CWorkspace: need 2*N [tauq, taup] + N [work]
1047* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1048* RWorkspace: need 0
1049*
1050 CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1051 CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1052 $ work( nwork ), lwork-nwork+1, ierr )
1053*
1054* Generate Q in A
1055* CWorkspace: need 2*N [tauq, taup] + N [work]
1056* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1057* RWorkspace: need 0
1058*
1059 CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1060 $ work( nwork ), lwork-nwork+1, ierr )
1061*
1062 IF( lwork .GE. m*n + 3*n ) THEN
1063*
1064* WORK( IU ) is M by N
1065*
1066 ldwrku = m
1067 ELSE
1068*
1069* WORK(IU) is LDWRKU by N
1070*
1071 ldwrku = ( lwork - 3*n ) / n
1072 END IF
1073 nwork = iu + ldwrku*n
1074*
1075* Perform bidiagonal SVD, computing left singular vectors
1076* of bidiagonal matrix in RWORK(IRU) and computing right
1077* singular vectors of bidiagonal matrix in RWORK(IRVT)
1078* CWorkspace: need 0
1079* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1080*
1081 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1082 $ n, rwork( irvt ), n, dum, idum,
1083 $ rwork( nrwork ), iwork, info )
1084*
1085* Multiply real matrix RWORK(IRVT) by P**H in VT,
1086* storing the result in WORK(IU), copying to VT
1087* CWorkspace: need 2*N [tauq, taup] + N*N [U]
1088* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1089*
1090 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt,
1091 $ work( iu ), ldwrku, rwork( nrwork ) )
1092 CALL zlacpy( 'F', n, n, work( iu ), ldwrku, vt, ldvt )
1093*
1094* Multiply Q in A by real matrix RWORK(IRU), storing the
1095* result in WORK(IU), copying to A
1096* CWorkspace: need 2*N [tauq, taup] + N*N [U]
1097* CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1098* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1099* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1100*
1101 nrwork = irvt
1102 DO 20 i = 1, m, ldwrku
1103 chunk = min( m-i+1, ldwrku )
1104 CALL zlacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1105 $ n, work( iu ), ldwrku, rwork( nrwork ) )
1106 CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1107 $ a( i, 1 ), lda )
1108 20 CONTINUE
1109*
1110 ELSE IF( wntqs ) THEN
1111*
1112* Path 5s (M >> N, JOBZ='S')
1113* Copy A to VT, generate P**H
1114* CWorkspace: need 2*N [tauq, taup] + N [work]
1115* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1116* RWorkspace: need 0
1117*
1118 CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1119 CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1120 $ work( nwork ), lwork-nwork+1, ierr )
1121*
1122* Copy A to U, generate Q
1123* CWorkspace: need 2*N [tauq, taup] + N [work]
1124* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1125* RWorkspace: need 0
1126*
1127 CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1128 CALL zungbr( 'Q', m, n, n, u, ldu, work( itauq ),
1129 $ work( nwork ), lwork-nwork+1, ierr )
1130*
1131* Perform bidiagonal SVD, computing left singular vectors
1132* of bidiagonal matrix in RWORK(IRU) and computing right
1133* singular vectors of bidiagonal matrix in RWORK(IRVT)
1134* CWorkspace: need 0
1135* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1136*
1137 iru = nrwork
1138 irvt = iru + n*n
1139 nrwork = irvt + n*n
1140 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1141 $ n, rwork( irvt ), n, dum, idum,
1142 $ rwork( nrwork ), iwork, info )
1143*
1144* Multiply real matrix RWORK(IRVT) by P**H in VT,
1145* storing the result in A, copying to VT
1146* CWorkspace: need 0
1147* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1148*
1149 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1150 $ rwork( nrwork ) )
1151 CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1152*
1153* Multiply Q in U by real matrix RWORK(IRU), storing the
1154* result in A, copying to U
1155* CWorkspace: need 0
1156* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1157*
1158 nrwork = irvt
1159 CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1160 $ rwork( nrwork ) )
1161 CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1162 ELSE
1163*
1164* Path 5a (M >> N, JOBZ='A')
1165* Copy A to VT, generate P**H
1166* CWorkspace: need 2*N [tauq, taup] + N [work]
1167* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1168* RWorkspace: need 0
1169*
1170 CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1171 CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1172 $ work( nwork ), lwork-nwork+1, ierr )
1173*
1174* Copy A to U, generate Q
1175* CWorkspace: need 2*N [tauq, taup] + M [work]
1176* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1177* RWorkspace: need 0
1178*
1179 CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1180 CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1181 $ work( nwork ), lwork-nwork+1, ierr )
1182*
1183* Perform bidiagonal SVD, computing left singular vectors
1184* of bidiagonal matrix in RWORK(IRU) and computing right
1185* singular vectors of bidiagonal matrix in RWORK(IRVT)
1186* CWorkspace: need 0
1187* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1188*
1189 iru = nrwork
1190 irvt = iru + n*n
1191 nrwork = irvt + n*n
1192 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1193 $ n, rwork( irvt ), n, dum, idum,
1194 $ rwork( nrwork ), iwork, info )
1195*
1196* Multiply real matrix RWORK(IRVT) by P**H in VT,
1197* storing the result in A, copying to VT
1198* CWorkspace: need 0
1199* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1200*
1201 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1202 $ rwork( nrwork ) )
1203 CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1204*
1205* Multiply Q in U by real matrix RWORK(IRU), storing the
1206* result in A, copying to U
1207* CWorkspace: need 0
1208* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1209*
1210 nrwork = irvt
1211 CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1212 $ rwork( nrwork ) )
1213 CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1214 END IF
1215*
1216 ELSE
1217*
1218* M .LT. MNTHR2
1219*
1220* Path 6 (M >= N, but not much larger)
1221* Reduce to bidiagonal form without QR decomposition
1222* Use ZUNMBR to compute singular vectors
1223*
1224 ie = 1
1225 nrwork = ie + n
1226 itauq = 1
1227 itaup = itauq + n
1228 nwork = itaup + n
1229*
1230* Bidiagonalize A
1231* CWorkspace: need 2*N [tauq, taup] + M [work]
1232* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1233* RWorkspace: need N [e]
1234*
1235 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1236 $ work( itaup ), work( nwork ), lwork-nwork+1,
1237 $ ierr )
1238 IF( wntqn ) THEN
1239*
1240* Path 6n (M >= N, JOBZ='N')
1241* Compute singular values only
1242* CWorkspace: need 0
1243* RWorkspace: need N [e] + BDSPAC
1244*
1245 CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
1246 $ dum, idum, rwork( nrwork ), iwork, info )
1247 ELSE IF( wntqo ) THEN
1248 iu = nwork
1249 iru = nrwork
1250 irvt = iru + n*n
1251 nrwork = irvt + n*n
1252 IF( lwork .GE. m*n + 3*n ) THEN
1253*
1254* WORK( IU ) is M by N
1255*
1256 ldwrku = m
1257 ELSE
1258*
1259* WORK( IU ) is LDWRKU by N
1260*
1261 ldwrku = ( lwork - 3*n ) / n
1262 END IF
1263 nwork = iu + ldwrku*n
1264*
1265* Path 6o (M >= N, JOBZ='O')
1266* Perform bidiagonal SVD, computing left singular vectors
1267* of bidiagonal matrix in RWORK(IRU) and computing right
1268* singular vectors of bidiagonal matrix in RWORK(IRVT)
1269* CWorkspace: need 0
1270* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1271*
1272 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1273 $ n, rwork( irvt ), n, dum, idum,
1274 $ rwork( nrwork ), iwork, info )
1275*
1276* Copy real matrix RWORK(IRVT) to complex matrix VT
1277* Overwrite VT by right singular vectors of A
1278* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1279* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1280* RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1281*
1282 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1283 CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1284 $ work( itaup ), vt, ldvt, work( nwork ),
1285 $ lwork-nwork+1, ierr )
1286*
1287 IF( lwork .GE. m*n + 3*n ) THEN
1288*
1289* Path 6o-fast
1290* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1291* Overwrite WORK(IU) by left singular vectors of A, copying
1292* to A
1293* CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work]
1294* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
1295* RWorkspace: need N [e] + N*N [RU]
1296*
1297 CALL zlaset( 'F', m, n, czero, czero, work( iu ),
1298 $ ldwrku )
1299 CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
1300 $ ldwrku )
1301 CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1302 $ work( itauq ), work( iu ), ldwrku,
1303 $ work( nwork ), lwork-nwork+1, ierr )
1304 CALL zlacpy( 'F', m, n, work( iu ), ldwrku, a, lda )
1305 ELSE
1306*
1307* Path 6o-slow
1308* Generate Q in A
1309* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1310* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1311* RWorkspace: need 0
1312*
1313 CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1314 $ work( nwork ), lwork-nwork+1, ierr )
1315*
1316* Multiply Q in A by real matrix RWORK(IRU), storing the
1317* result in WORK(IU), copying to A
1318* CWorkspace: need 2*N [tauq, taup] + N*N [U]
1319* CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1320* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1321* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1322*
1323 nrwork = irvt
1324 DO 30 i = 1, m, ldwrku
1325 chunk = min( m-i+1, ldwrku )
1326 CALL zlacrm( chunk, n, a( i, 1 ), lda,
1327 $ rwork( iru ), n, work( iu ), ldwrku,
1328 $ rwork( nrwork ) )
1329 CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1330 $ a( i, 1 ), lda )
1331 30 CONTINUE
1332 END IF
1333*
1334 ELSE IF( wntqs ) THEN
1335*
1336* Path 6s (M >= N, JOBZ='S')
1337* Perform bidiagonal SVD, computing left singular vectors
1338* of bidiagonal matrix in RWORK(IRU) and computing right
1339* singular vectors of bidiagonal matrix in RWORK(IRVT)
1340* CWorkspace: need 0
1341* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1342*
1343 iru = nrwork
1344 irvt = iru + n*n
1345 nrwork = irvt + n*n
1346 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1347 $ n, rwork( irvt ), n, dum, idum,
1348 $ rwork( nrwork ), iwork, info )
1349*
1350* Copy real matrix RWORK(IRU) to complex matrix U
1351* Overwrite U by left singular vectors of A
1352* CWorkspace: need 2*N [tauq, taup] + N [work]
1353* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1354* RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1355*
1356 CALL zlaset( 'F', m, n, czero, czero, u, ldu )
1357 CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1358 CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1359 $ work( itauq ), u, ldu, work( nwork ),
1360 $ lwork-nwork+1, ierr )
1361*
1362* Copy real matrix RWORK(IRVT) to complex matrix VT
1363* Overwrite VT by right singular vectors of A
1364* CWorkspace: need 2*N [tauq, taup] + N [work]
1365* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1366* RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1367*
1368 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1369 CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1370 $ work( itaup ), vt, ldvt, work( nwork ),
1371 $ lwork-nwork+1, ierr )
1372 ELSE
1373*
1374* Path 6a (M >= N, JOBZ='A')
1375* Perform bidiagonal SVD, computing left singular vectors
1376* of bidiagonal matrix in RWORK(IRU) and computing right
1377* singular vectors of bidiagonal matrix in RWORK(IRVT)
1378* CWorkspace: need 0
1379* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1380*
1381 iru = nrwork
1382 irvt = iru + n*n
1383 nrwork = irvt + n*n
1384 CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1385 $ n, rwork( irvt ), n, dum, idum,
1386 $ rwork( nrwork ), iwork, info )
1387*
1388* Set the right corner of U to identity matrix
1389*
1390 CALL zlaset( 'F', m, m, czero, czero, u, ldu )
1391 IF( m.GT.n ) THEN
1392 CALL zlaset( 'F', m-n, m-n, czero, cone,
1393 $ u( n+1, n+1 ), ldu )
1394 END IF
1395*
1396* Copy real matrix RWORK(IRU) to complex matrix U
1397* Overwrite U by left singular vectors of A
1398* CWorkspace: need 2*N [tauq, taup] + M [work]
1399* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1400* RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1401*
1402 CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1403 CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
1404 $ work( itauq ), u, ldu, work( nwork ),
1405 $ lwork-nwork+1, ierr )
1406*
1407* Copy real matrix RWORK(IRVT) to complex matrix VT
1408* Overwrite VT by right singular vectors of A
1409* CWorkspace: need 2*N [tauq, taup] + N [work]
1410* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1411* RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1412*
1413 CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1414 CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1415 $ work( itaup ), vt, ldvt, work( nwork ),
1416 $ lwork-nwork+1, ierr )
1417 END IF
1418*
1419 END IF
1420*
1421 ELSE
1422*
1423* A has more columns than rows. If A has sufficiently more
1424* columns than rows, first reduce using the LQ decomposition (if
1425* sufficient workspace available)
1426*
1427 IF( n.GE.mnthr1 ) THEN
1428*
1429 IF( wntqn ) THEN
1430*
1431* Path 1t (N >> M, JOBZ='N')
1432* No singular vectors to be computed
1433*
1434 itau = 1
1435 nwork = itau + m
1436*
1437* Compute A=L*Q
1438* CWorkspace: need M [tau] + M [work]
1439* CWorkspace: prefer M [tau] + M*NB [work]
1440* RWorkspace: need 0
1441*
1442 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1443 $ lwork-nwork+1, ierr )
1444*
1445* Zero out above L
1446*
1447 CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1448 $ lda )
1449 ie = 1
1450 itauq = 1
1451 itaup = itauq + m
1452 nwork = itaup + m
1453*
1454* Bidiagonalize L in A
1455* CWorkspace: need 2*M [tauq, taup] + M [work]
1456* CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
1457* RWorkspace: need M [e]
1458*
1459 CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1460 $ work( itaup ), work( nwork ), lwork-nwork+1,
1461 $ ierr )
1462 nrwork = ie + m
1463*
1464* Perform bidiagonal SVD, compute singular values only
1465* CWorkspace: need 0
1466* RWorkspace: need M [e] + BDSPAC
1467*
1468 CALL dbdsdc( 'U', 'N', m, s, rwork( ie ), dum,1,dum,1,
1469 $ dum, idum, rwork( nrwork ), iwork, info )
1470*
1471 ELSE IF( wntqo ) THEN
1472*
1473* Path 2t (N >> M, JOBZ='O')
1474* M right singular vectors to be overwritten on A and
1475* M left singular vectors to be computed in U
1476*
1477 ivt = 1
1478 ldwkvt = m
1479*
1480* WORK(IVT) is M by M
1481*
1482 il = ivt + ldwkvt*m
1483 IF( lwork .GE. m*n + m*m + 3*m ) THEN
1484*
1485* WORK(IL) M by N
1486*
1487 ldwrkl = m
1488 chunk = n
1489 ELSE
1490*
1491* WORK(IL) is M by CHUNK
1492*
1493 ldwrkl = m
1494 chunk = ( lwork - m*m - 3*m ) / m
1495 END IF
1496 itau = il + ldwrkl*chunk
1497 nwork = itau + m
1498*
1499* Compute A=L*Q
1500* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1501* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1502* RWorkspace: need 0
1503*
1504 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1505 $ lwork-nwork+1, ierr )
1506*
1507* Copy L to WORK(IL), zeroing about above it
1508*
1509 CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1510 CALL zlaset( 'U', m-1, m-1, czero, czero,
1511 $ work( il+ldwrkl ), ldwrkl )
1512*
1513* Generate Q in A
1514* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1515* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1516* RWorkspace: need 0
1517*
1518 CALL zunglq( m, n, m, a, lda, work( itau ),
1519 $ work( nwork ), lwork-nwork+1, ierr )
1520 ie = 1
1521 itauq = itau
1522 itaup = itauq + m
1523 nwork = itaup + m
1524*
1525* Bidiagonalize L in WORK(IL)
1526* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1527* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1528* RWorkspace: need M [e]
1529*
1530 CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1531 $ work( itauq ), work( itaup ), work( nwork ),
1532 $ lwork-nwork+1, ierr )
1533*
1534* Perform bidiagonal SVD, computing left singular vectors
1535* of bidiagonal matrix in RWORK(IRU) and computing right
1536* singular vectors of bidiagonal matrix in RWORK(IRVT)
1537* CWorkspace: need 0
1538* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1539*
1540 iru = ie + m
1541 irvt = iru + m*m
1542 nrwork = irvt + m*m
1543 CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1544 $ m, rwork( irvt ), m, dum, idum,
1545 $ rwork( nrwork ), iwork, info )
1546*
1547* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1548* Overwrite WORK(IU) by the left singular vectors of L
1549* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1550* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1551* RWorkspace: need 0
1552*
1553 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1554 CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1555 $ work( itauq ), u, ldu, work( nwork ),
1556 $ lwork-nwork+1, ierr )
1557*
1558* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1559* Overwrite WORK(IVT) by the right singular vectors of L
1560* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1561* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1562* RWorkspace: need 0
1563*
1564 CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1565 $ ldwkvt )
1566 CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1567 $ work( itaup ), work( ivt ), ldwkvt,
1568 $ work( nwork ), lwork-nwork+1, ierr )
1569*
1570* Multiply right singular vectors of L in WORK(IL) by Q
1571* in A, storing result in WORK(IL) and copying to A
1572* CWorkspace: need M*M [VT] + M*M [L]
1573* CWorkspace: prefer M*M [VT] + M*N [L]
1574* RWorkspace: need 0
1575*
1576 DO 40 i = 1, n, chunk
1577 blk = min( n-i+1, chunk )
1578 CALL zgemm( 'N', 'N', m, blk, m, cone, work( ivt ), m,
1579 $ a( 1, i ), lda, czero, work( il ),
1580 $ ldwrkl )
1581 CALL zlacpy( 'F', m, blk, work( il ), ldwrkl,
1582 $ a( 1, i ), lda )
1583 40 CONTINUE
1584*
1585 ELSE IF( wntqs ) THEN
1586*
1587* Path 3t (N >> M, JOBZ='S')
1588* M right singular vectors to be computed in VT and
1589* M left singular vectors to be computed in U
1590*
1591 il = 1
1592*
1593* WORK(IL) is M by M
1594*
1595 ldwrkl = m
1596 itau = il + ldwrkl*m
1597 nwork = itau + m
1598*
1599* Compute A=L*Q
1600* CWorkspace: need M*M [L] + M [tau] + M [work]
1601* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1602* RWorkspace: need 0
1603*
1604 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1605 $ lwork-nwork+1, ierr )
1606*
1607* Copy L to WORK(IL), zeroing out above it
1608*
1609 CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1610 CALL zlaset( 'U', m-1, m-1, czero, czero,
1611 $ work( il+ldwrkl ), ldwrkl )
1612*
1613* Generate Q in A
1614* CWorkspace: need M*M [L] + M [tau] + M [work]
1615* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1616* RWorkspace: need 0
1617*
1618 CALL zunglq( m, n, m, a, lda, work( itau ),
1619 $ work( nwork ), lwork-nwork+1, ierr )
1620 ie = 1
1621 itauq = itau
1622 itaup = itauq + m
1623 nwork = itaup + m
1624*
1625* Bidiagonalize L in WORK(IL)
1626* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1627* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1628* RWorkspace: need M [e]
1629*
1630 CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1631 $ work( itauq ), work( itaup ), work( nwork ),
1632 $ lwork-nwork+1, ierr )
1633*
1634* Perform bidiagonal SVD, computing left singular vectors
1635* of bidiagonal matrix in RWORK(IRU) and computing right
1636* singular vectors of bidiagonal matrix in RWORK(IRVT)
1637* CWorkspace: need 0
1638* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1639*
1640 iru = ie + m
1641 irvt = iru + m*m
1642 nrwork = irvt + m*m
1643 CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1644 $ m, rwork( irvt ), m, dum, idum,
1645 $ rwork( nrwork ), iwork, info )
1646*
1647* Copy real matrix RWORK(IRU) to complex matrix U
1648* Overwrite U by left singular vectors of L
1649* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1650* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1651* RWorkspace: need 0
1652*
1653 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1654 CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1655 $ work( itauq ), u, ldu, work( nwork ),
1656 $ lwork-nwork+1, ierr )
1657*
1658* Copy real matrix RWORK(IRVT) to complex matrix VT
1659* Overwrite VT by left singular vectors of L
1660* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1661* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1662* RWorkspace: need 0
1663*
1664 CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
1665 CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1666 $ work( itaup ), vt, ldvt, work( nwork ),
1667 $ lwork-nwork+1, ierr )
1668*
1669* Copy VT to WORK(IL), multiply right singular vectors of L
1670* in WORK(IL) by Q in A, storing result in VT
1671* CWorkspace: need M*M [L]
1672* RWorkspace: need 0
1673*
1674 CALL zlacpy( 'F', m, m, vt, ldvt, work( il ), ldwrkl )
1675 CALL zgemm( 'N', 'N', m, n, m, cone, work( il ), ldwrkl,
1676 $ a, lda, czero, vt, ldvt )
1677*
1678 ELSE IF( wntqa ) THEN
1679*
1680* Path 4t (N >> M, JOBZ='A')
1681* N right singular vectors to be computed in VT and
1682* M left singular vectors to be computed in U
1683*
1684 ivt = 1
1685*
1686* WORK(IVT) is M by M
1687*
1688 ldwkvt = m
1689 itau = ivt + ldwkvt*m
1690 nwork = itau + m
1691*
1692* Compute A=L*Q, copying result to VT
1693* CWorkspace: need M*M [VT] + M [tau] + M [work]
1694* CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
1695* RWorkspace: need 0
1696*
1697 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1698 $ lwork-nwork+1, ierr )
1699 CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1700*
1701* Generate Q in VT
1702* CWorkspace: need M*M [VT] + M [tau] + N [work]
1703* CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
1704* RWorkspace: need 0
1705*
1706 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
1707 $ work( nwork ), lwork-nwork+1, ierr )
1708*
1709* Produce L in A, zeroing out above it
1710*
1711 CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1712 $ lda )
1713 ie = 1
1714 itauq = itau
1715 itaup = itauq + m
1716 nwork = itaup + m
1717*
1718* Bidiagonalize L in A
1719* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1720* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
1721* RWorkspace: need M [e]
1722*
1723 CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1724 $ work( itaup ), work( nwork ), lwork-nwork+1,
1725 $ ierr )
1726*
1727* Perform bidiagonal SVD, computing left singular vectors
1728* of bidiagonal matrix in RWORK(IRU) and computing right
1729* singular vectors of bidiagonal matrix in RWORK(IRVT)
1730* CWorkspace: need 0
1731* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1732*
1733 iru = ie + m
1734 irvt = iru + m*m
1735 nrwork = irvt + m*m
1736 CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1737 $ m, rwork( irvt ), m, dum, idum,
1738 $ rwork( nrwork ), iwork, info )
1739*
1740* Copy real matrix RWORK(IRU) to complex matrix U
1741* Overwrite U by left singular vectors of L
1742* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1743* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1744* RWorkspace: need 0
1745*
1746 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1747 CALL zunmbr( 'Q', 'L', 'N', m, m, m, a, lda,
1748 $ work( itauq ), u, ldu, work( nwork ),
1749 $ lwork-nwork+1, ierr )
1750*
1751* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1752* Overwrite WORK(IVT) by right singular vectors of L
1753* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1754* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1755* RWorkspace: need 0
1756*
1757 CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1758 $ ldwkvt )
1759 CALL zunmbr( 'P', 'R', 'C', m, m, m, a, lda,
1760 $ work( itaup ), work( ivt ), ldwkvt,
1761 $ work( nwork ), lwork-nwork+1, ierr )
1762*
1763* Multiply right singular vectors of L in WORK(IVT) by
1764* Q in VT, storing result in A
1765* CWorkspace: need M*M [VT]
1766* RWorkspace: need 0
1767*
1768 CALL zgemm( 'N', 'N', m, n, m, cone, work( ivt ), ldwkvt,
1769 $ vt, ldvt, czero, a, lda )
1770*
1771* Copy right singular vectors of A from A to VT
1772*
1773 CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1774*
1775 END IF
1776*
1777 ELSE IF( n.GE.mnthr2 ) THEN
1778*
1779* MNTHR2 <= N < MNTHR1
1780*
1781* Path 5t (N >> M, but not as much as MNTHR1)
1782* Reduce to bidiagonal form without QR decomposition, use
1783* ZUNGBR and matrix multiplication to compute singular vectors
1784*
1785 ie = 1
1786 nrwork = ie + m
1787 itauq = 1
1788 itaup = itauq + m
1789 nwork = itaup + m
1790*
1791* Bidiagonalize A
1792* CWorkspace: need 2*M [tauq, taup] + N [work]
1793* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
1794* RWorkspace: need M [e]
1795*
1796 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1797 $ work( itaup ), work( nwork ), lwork-nwork+1,
1798 $ ierr )
1799*
1800 IF( wntqn ) THEN
1801*
1802* Path 5tn (N >> M, JOBZ='N')
1803* Compute singular values only
1804* CWorkspace: need 0
1805* RWorkspace: need M [e] + BDSPAC
1806*
1807 CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
1808 $ dum, idum, rwork( nrwork ), iwork, info )
1809 ELSE IF( wntqo ) THEN
1810 irvt = nrwork
1811 iru = irvt + m*m
1812 nrwork = iru + m*m
1813 ivt = nwork
1814*
1815* Path 5to (N >> M, JOBZ='O')
1816* Copy A to U, generate Q
1817* CWorkspace: need 2*M [tauq, taup] + M [work]
1818* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1819* RWorkspace: need 0
1820*
1821 CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1822 CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1823 $ work( nwork ), lwork-nwork+1, ierr )
1824*
1825* Generate P**H in A
1826* CWorkspace: need 2*M [tauq, taup] + M [work]
1827* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1828* RWorkspace: need 0
1829*
1830 CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
1831 $ work( nwork ), lwork-nwork+1, ierr )
1832*
1833 ldwkvt = m
1834 IF( lwork .GE. m*n + 3*m ) THEN
1835*
1836* WORK( IVT ) is M by N
1837*
1838 nwork = ivt + ldwkvt*n
1839 chunk = n
1840 ELSE
1841*
1842* WORK( IVT ) is M by CHUNK
1843*
1844 chunk = ( lwork - 3*m ) / m
1845 nwork = ivt + ldwkvt*chunk
1846 END IF
1847*
1848* Perform bidiagonal SVD, computing left singular vectors
1849* of bidiagonal matrix in RWORK(IRU) and computing right
1850* singular vectors of bidiagonal matrix in RWORK(IRVT)
1851* CWorkspace: need 0
1852* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1853*
1854 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1855 $ m, rwork( irvt ), m, dum, idum,
1856 $ rwork( nrwork ), iwork, info )
1857*
1858* Multiply Q in U by real matrix RWORK(IRVT)
1859* storing the result in WORK(IVT), copying to U
1860* CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1861* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1862*
1863 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1864 $ ldwkvt, rwork( nrwork ) )
1865 CALL zlacpy( 'F', m, m, work( ivt ), ldwkvt, u, ldu )
1866*
1867* Multiply RWORK(IRVT) by P**H in A, storing the
1868* result in WORK(IVT), copying to A
1869* CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1870* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
1871* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
1872* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1873*
1874 nrwork = iru
1875 DO 50 i = 1, n, chunk
1876 blk = min( n-i+1, chunk )
1877 CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1878 $ work( ivt ), ldwkvt, rwork( nrwork ) )
1879 CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
1880 $ a( 1, i ), lda )
1881 50 CONTINUE
1882 ELSE IF( wntqs ) THEN
1883*
1884* Path 5ts (N >> M, JOBZ='S')
1885* Copy A to U, generate Q
1886* CWorkspace: need 2*M [tauq, taup] + M [work]
1887* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1888* RWorkspace: need 0
1889*
1890 CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1891 CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1892 $ work( nwork ), lwork-nwork+1, ierr )
1893*
1894* Copy A to VT, generate P**H
1895* CWorkspace: need 2*M [tauq, taup] + M [work]
1896* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1897* RWorkspace: need 0
1898*
1899 CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1900 CALL zungbr( 'P', m, n, m, vt, ldvt, work( itaup ),
1901 $ work( nwork ), lwork-nwork+1, ierr )
1902*
1903* Perform bidiagonal SVD, computing left singular vectors
1904* of bidiagonal matrix in RWORK(IRU) and computing right
1905* singular vectors of bidiagonal matrix in RWORK(IRVT)
1906* CWorkspace: need 0
1907* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1908*
1909 irvt = nrwork
1910 iru = irvt + m*m
1911 nrwork = iru + m*m
1912 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1913 $ m, rwork( irvt ), m, dum, idum,
1914 $ rwork( nrwork ), iwork, info )
1915*
1916* Multiply Q in U by real matrix RWORK(IRU), storing the
1917* result in A, copying to U
1918* CWorkspace: need 0
1919* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1920*
1921 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1922 $ rwork( nrwork ) )
1923 CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1924*
1925* Multiply real matrix RWORK(IRVT) by P**H in VT,
1926* storing the result in A, copying to VT
1927* CWorkspace: need 0
1928* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1929*
1930 nrwork = iru
1931 CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1932 $ rwork( nrwork ) )
1933 CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1934 ELSE
1935*
1936* Path 5ta (N >> M, JOBZ='A')
1937* Copy A to U, generate Q
1938* CWorkspace: need 2*M [tauq, taup] + M [work]
1939* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1940* RWorkspace: need 0
1941*
1942 CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1943 CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1944 $ work( nwork ), lwork-nwork+1, ierr )
1945*
1946* Copy A to VT, generate P**H
1947* CWorkspace: need 2*M [tauq, taup] + N [work]
1948* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
1949* RWorkspace: need 0
1950*
1951 CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1952 CALL zungbr( 'P', n, n, m, vt, ldvt, work( itaup ),
1953 $ work( nwork ), lwork-nwork+1, ierr )
1954*
1955* Perform bidiagonal SVD, computing left singular vectors
1956* of bidiagonal matrix in RWORK(IRU) and computing right
1957* singular vectors of bidiagonal matrix in RWORK(IRVT)
1958* CWorkspace: need 0
1959* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1960*
1961 irvt = nrwork
1962 iru = irvt + m*m
1963 nrwork = iru + m*m
1964 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1965 $ m, rwork( irvt ), m, dum, idum,
1966 $ rwork( nrwork ), iwork, info )
1967*
1968* Multiply Q in U by real matrix RWORK(IRU), storing the
1969* result in A, copying to U
1970* CWorkspace: need 0
1971* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1972*
1973 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1974 $ rwork( nrwork ) )
1975 CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1976*
1977* Multiply real matrix RWORK(IRVT) by P**H in VT,
1978* storing the result in A, copying to VT
1979* CWorkspace: need 0
1980* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1981*
1982 nrwork = iru
1983 CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1984 $ rwork( nrwork ) )
1985 CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1986 END IF
1987*
1988 ELSE
1989*
1990* N .LT. MNTHR2
1991*
1992* Path 6t (N > M, but not much larger)
1993* Reduce to bidiagonal form without LQ decomposition
1994* Use ZUNMBR to compute singular vectors
1995*
1996 ie = 1
1997 nrwork = ie + m
1998 itauq = 1
1999 itaup = itauq + m
2000 nwork = itaup + m
2001*
2002* Bidiagonalize A
2003* CWorkspace: need 2*M [tauq, taup] + N [work]
2004* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
2005* RWorkspace: need M [e]
2006*
2007 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2008 $ work( itaup ), work( nwork ), lwork-nwork+1,
2009 $ ierr )
2010 IF( wntqn ) THEN
2011*
2012* Path 6tn (N > M, JOBZ='N')
2013* Compute singular values only
2014* CWorkspace: need 0
2015* RWorkspace: need M [e] + BDSPAC
2016*
2017 CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
2018 $ dum, idum, rwork( nrwork ), iwork, info )
2019 ELSE IF( wntqo ) THEN
2020* Path 6to (N > M, JOBZ='O')
2021 ldwkvt = m
2022 ivt = nwork
2023 IF( lwork .GE. m*n + 3*m ) THEN
2024*
2025* WORK( IVT ) is M by N
2026*
2027 CALL zlaset( 'F', m, n, czero, czero, work( ivt ),
2028 $ ldwkvt )
2029 nwork = ivt + ldwkvt*n
2030 ELSE
2031*
2032* WORK( IVT ) is M by CHUNK
2033*
2034 chunk = ( lwork - 3*m ) / m
2035 nwork = ivt + ldwkvt*chunk
2036 END IF
2037*
2038* Perform bidiagonal SVD, computing left singular vectors
2039* of bidiagonal matrix in RWORK(IRU) and computing right
2040* singular vectors of bidiagonal matrix in RWORK(IRVT)
2041* CWorkspace: need 0
2042* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2043*
2044 irvt = nrwork
2045 iru = irvt + m*m
2046 nrwork = iru + m*m
2047 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2048 $ m, rwork( irvt ), m, dum, idum,
2049 $ rwork( nrwork ), iwork, info )
2050*
2051* Copy real matrix RWORK(IRU) to complex matrix U
2052* Overwrite U by left singular vectors of A
2053* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2054* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2055* RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2056*
2057 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2058 CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2059 $ work( itauq ), u, ldu, work( nwork ),
2060 $ lwork-nwork+1, ierr )
2061*
2062 IF( lwork .GE. m*n + 3*m ) THEN
2063*
2064* Path 6to-fast
2065* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
2066* Overwrite WORK(IVT) by right singular vectors of A,
2067* copying to A
2068* CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work]
2069* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
2070* RWorkspace: need M [e] + M*M [RVT]
2071*
2072 CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
2073 $ ldwkvt )
2074 CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2075 $ work( itaup ), work( ivt ), ldwkvt,
2076 $ work( nwork ), lwork-nwork+1, ierr )
2077 CALL zlacpy( 'F', m, n, work( ivt ), ldwkvt, a, lda )
2078 ELSE
2079*
2080* Path 6to-slow
2081* Generate P**H in A
2082* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2083* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2084* RWorkspace: need 0
2085*
2086 CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
2087 $ work( nwork ), lwork-nwork+1, ierr )
2088*
2089* Multiply Q in A by real matrix RWORK(IRU), storing the
2090* result in WORK(IU), copying to A
2091* CWorkspace: need 2*M [tauq, taup] + M*M [VT]
2092* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
2093* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
2094* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
2095*
2096 nrwork = iru
2097 DO 60 i = 1, n, chunk
2098 blk = min( n-i+1, chunk )
2099 CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ),
2100 $ lda, work( ivt ), ldwkvt,
2101 $ rwork( nrwork ) )
2102 CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
2103 $ a( 1, i ), lda )
2104 60 CONTINUE
2105 END IF
2106 ELSE IF( wntqs ) THEN
2107*
2108* Path 6ts (N > M, JOBZ='S')
2109* Perform bidiagonal SVD, computing left singular vectors
2110* of bidiagonal matrix in RWORK(IRU) and computing right
2111* singular vectors of bidiagonal matrix in RWORK(IRVT)
2112* CWorkspace: need 0
2113* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2114*
2115 irvt = nrwork
2116 iru = irvt + m*m
2117 nrwork = iru + m*m
2118 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2119 $ m, rwork( irvt ), m, dum, idum,
2120 $ rwork( nrwork ), iwork, info )
2121*
2122* Copy real matrix RWORK(IRU) to complex matrix U
2123* Overwrite U by left singular vectors of A
2124* CWorkspace: need 2*M [tauq, taup] + M [work]
2125* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2126* RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2127*
2128 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2129 CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2130 $ work( itauq ), u, ldu, work( nwork ),
2131 $ lwork-nwork+1, ierr )
2132*
2133* Copy real matrix RWORK(IRVT) to complex matrix VT
2134* Overwrite VT by right singular vectors of A
2135* CWorkspace: need 2*M [tauq, taup] + M [work]
2136* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2137* RWorkspace: need M [e] + M*M [RVT]
2138*
2139 CALL zlaset( 'F', m, n, czero, czero, vt, ldvt )
2140 CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2141 CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2142 $ work( itaup ), vt, ldvt, work( nwork ),
2143 $ lwork-nwork+1, ierr )
2144 ELSE
2145*
2146* Path 6ta (N > M, JOBZ='A')
2147* Perform bidiagonal SVD, computing left singular vectors
2148* of bidiagonal matrix in RWORK(IRU) and computing right
2149* singular vectors of bidiagonal matrix in RWORK(IRVT)
2150* CWorkspace: need 0
2151* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2152*
2153 irvt = nrwork
2154 iru = irvt + m*m
2155 nrwork = iru + m*m
2156*
2157 CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2158 $ m, rwork( irvt ), m, dum, idum,
2159 $ rwork( nrwork ), iwork, info )
2160*
2161* Copy real matrix RWORK(IRU) to complex matrix U
2162* Overwrite U by left singular vectors of A
2163* CWorkspace: need 2*M [tauq, taup] + M [work]
2164* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2165* RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2166*
2167 CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2168 CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2169 $ work( itauq ), u, ldu, work( nwork ),
2170 $ lwork-nwork+1, ierr )
2171*
2172* Set all of VT to identity matrix
2173*
2174 CALL zlaset( 'F', n, n, czero, cone, vt, ldvt )
2175*
2176* Copy real matrix RWORK(IRVT) to complex matrix VT
2177* Overwrite VT by right singular vectors of A
2178* CWorkspace: need 2*M [tauq, taup] + N [work]
2179* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
2180* RWorkspace: need M [e] + M*M [RVT]
2181*
2182 CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2183 CALL zunmbr( 'P', 'R', 'C', n, n, m, a, lda,
2184 $ work( itaup ), vt, ldvt, work( nwork ),
2185 $ lwork-nwork+1, ierr )
2186 END IF
2187*
2188 END IF
2189*
2190 END IF
2191*
2192* Undo scaling if necessary
2193*
2194 IF( iscl.EQ.1 ) THEN
2195 IF( anrm.GT.bignum )
2196 $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2197 $ ierr )
2198 IF( info.NE.0 .AND. anrm.GT.bignum )
2199 $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn-1, 1,
2200 $ rwork( ie ), minmn, ierr )
2201 IF( anrm.LT.smlnum )
2202 $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2203 $ ierr )
2204 IF( info.NE.0 .AND. anrm.LT.smlnum )
2205 $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1,
2206 $ rwork( ie ), minmn, ierr )
2207 END IF
2208*
2209* Return optimal workspace in WORK(1)
2210*
2211 work( 1 ) = droundup_lwork( maxwrk )
2212*
2213 RETURN
2214*
2215* End of ZGESDD
2216*
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
zgebrd
subroutine zgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
ZGEBRD
Definition zgebrd.f:205
zgelqf
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
Definition zgelqf.f:143
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:146
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:59
zlacp2
subroutine zlacp2(uplo, m, n, a, lda, b, ldb)
ZLACP2 copies all or part of a real two-dimensional array to a complex array.
Definition zlacp2.f:104
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zlacrm
subroutine zlacrm(m, n, a, lda, b, ldb, c, ldc, rwork)
ZLACRM multiplies a complex matrix by a square real matrix.
Definition zlacrm.f:114
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:115
zlarcm
subroutine zlarcm(m, n, a, lda, b, ldb, c, ldc, rwork)
ZLARCM copies all or part of a real two-dimensional array to a complex array.
Definition zlarcm.f:114
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
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
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
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
zungbr
subroutine zungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
ZUNGBR
Definition zungbr.f:157
zunglq
subroutine zunglq(m, n, k, a, lda, tau, work, lwork, info)
ZUNGLQ
Definition zunglq.f:127
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:128
zunmbr
subroutine zunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMBR
Definition zunmbr.f:196
Here is the call graph for this function:
Here is the caller graph for this function: