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

◆ cgesdd()

subroutine cgesdd ( character  JOBZ,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  S,
complex, dimension( ldu, * )  U,
integer  LDU,
complex, dimension( ldvt, * )  VT,
integer  LDVT,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

CGESDD

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

Purpose:
 CGESDD 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.

 The divide and conquer algorithm makes very mild assumptions about
 floating point arithmetic. It will work on machines with a guard
 digit in add/subtract, or on those binary machines without guard
 digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
 Cray-2. It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.
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 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 REAL array, dimension (min(M,N))
          The singular values of A, sorted so that S(i) >= S(i+1).
[out]U
          U is COMPLEX 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 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 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 REAL 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 SBDSDC 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 225 of file cgesdd.f.

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