214 SUBROUTINE cgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
215 $ work, lwork, rwork, info )
223 CHARACTER jobu, jobvt
224 INTEGER info, lda, ldu, ldvt, lwork, m, n
227 REAL rwork( * ), s( * )
228 COMPLEX a( lda, * ), u( ldu, * ), vt( ldvt, * ),
236 parameter( czero = ( 0.0e0, 0.0e0 ),
237 $ cone = ( 1.0e0, 0.0e0 ) )
239 parameter( zero = 0.0e0, one = 1.0e0 )
242 LOGICAL lquery, wntua, wntuas, wntun, wntuo, wntus,
243 $ wntva, wntvas, wntvn, wntvo, wntvs
244 INTEGER blk, chunk, i, ie, ierr, ir, irwork, iscl,
245 $ itau, itaup, itauq, iu, iwork, ldwrkr, ldwrku,
246 $ maxwrk, minmn, minwrk, mnthr, ncu, ncvt, nru,
248 INTEGER lwork_cgeqrf, lwork_cungqr_n, lwork_cungqr_m,
249 $ lwork_cgebrd, lwork_cungbr_p, lwork_cungbr_q,
250 $ lwork_cgelqf, lwork_cunglq_n, lwork_cunglq_m
251 REAL anrm, bignum, eps, smlnum
269 INTRINSIC max, min, sqrt
277 wntua =
lsame( jobu,
'A' )
278 wntus =
lsame( jobu,
'S' )
279 wntuas = wntua .OR. wntus
280 wntuo =
lsame( jobu,
'O' )
281 wntun =
lsame( jobu,
'N' )
282 wntva =
lsame( jobvt,
'A' )
283 wntvs =
lsame( jobvt,
'S' )
284 wntvas = wntva .OR. wntvs
285 wntvo =
lsame( jobvt,
'O' )
286 wntvn =
lsame( jobvt,
'N' )
287 lquery = ( lwork.EQ.-1 )
289 IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) )
THEN
291 ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
292 $ ( wntvo .AND. wntuo ) )
THEN
294 ELSE IF( m.LT.0 )
THEN
296 ELSE IF( n.LT.0 )
THEN
298 ELSE IF( lda.LT.max( 1, m ) )
THEN
300 ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) )
THEN
302 ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
303 $ ( wntvs .AND. ldvt.LT.minmn ) )
THEN
318 IF( m.GE.n .AND. minmn.GT.0 )
THEN
322 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
324 CALL
cgeqrf( m, n, a, lda, dum(1), dum(1), -1, ierr )
327 CALL
cungqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
328 lwork_cungqr_n=dum(1)
329 CALL
cungqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
330 lwork_cungqr_m=dum(1)
332 CALL
cgebrd( n, n, a, lda, s, dum(1), dum(1),
333 $ dum(1), dum(1), -1, ierr )
336 CALL
cungbr(
'P', n, n, n, a, lda, dum(1),
338 lwork_cungbr_p=dum(1)
339 CALL
cungbr(
'Q', n, n, n, a, lda, dum(1),
341 lwork_cungbr_q=dum(1)
343 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
344 IF( m.GE.mnthr )
THEN
349 maxwrk = n + lwork_cgeqrf
350 maxwrk = max( maxwrk, 2*n+lwork_cgebrd )
351 IF( wntvo .OR. wntvas )
352 $ maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
354 ELSE IF( wntuo .AND. wntvn )
THEN
358 wrkbl = n + lwork_cgeqrf
359 wrkbl = max( wrkbl, n+lwork_cungqr_n )
360 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
361 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
362 maxwrk = max( n*n+wrkbl, n*n+m*n )
364 ELSE IF( wntuo .AND. wntvas )
THEN
369 wrkbl = n + lwork_cgeqrf
370 wrkbl = max( wrkbl, n+lwork_cungqr_n )
371 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
372 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
373 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
374 maxwrk = max( n*n+wrkbl, n*n+m*n )
376 ELSE IF( wntus .AND. wntvn )
THEN
380 wrkbl = n + lwork_cgeqrf
381 wrkbl = max( wrkbl, n+lwork_cungqr_n )
382 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
383 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
386 ELSE IF( wntus .AND. wntvo )
THEN
390 wrkbl = n + lwork_cgeqrf
391 wrkbl = max( wrkbl, n+lwork_cungqr_n )
392 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
393 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
394 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
395 maxwrk = 2*n*n + wrkbl
397 ELSE IF( wntus .AND. wntvas )
THEN
402 wrkbl = n + lwork_cgeqrf
403 wrkbl = max( wrkbl, n+lwork_cungqr_n )
404 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
405 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
406 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
409 ELSE IF( wntua .AND. wntvn )
THEN
413 wrkbl = n + lwork_cgeqrf
414 wrkbl = max( wrkbl, n+lwork_cungqr_m )
415 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
416 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
419 ELSE IF( wntua .AND. wntvo )
THEN
423 wrkbl = n + lwork_cgeqrf
424 wrkbl = max( wrkbl, n+lwork_cungqr_m )
425 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
426 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
427 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
428 maxwrk = 2*n*n + wrkbl
430 ELSE IF( wntua .AND. wntvas )
THEN
435 wrkbl = n + lwork_cgeqrf
436 wrkbl = max( wrkbl, n+lwork_cungqr_m )
437 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
438 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
439 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
447 CALL
cgebrd( m, n, a, lda, s, dum(1), dum(1),
448 $ dum(1), dum(1), -1, ierr )
450 maxwrk = 2*n + lwork_cgebrd
451 IF( wntus .OR. wntuo )
THEN
452 CALL
cungbr(
'Q', m, n, n, a, lda, dum(1),
454 lwork_cungbr_q=dum(1)
455 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
458 CALL
cungbr(
'Q', m, m, n, a, lda, dum(1),
460 lwork_cungbr_q=dum(1)
461 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
463 IF( .NOT.wntvn )
THEN
464 maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
468 ELSE IF( minmn.GT.0 )
THEN
472 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
474 CALL
cgelqf( m, n, a, lda, dum(1), dum(1), -1, ierr )
477 CALL
cunglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
478 lwork_cunglq_n=dum(1)
479 CALL
cunglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
480 lwork_cunglq_m=dum(1)
482 CALL
cgebrd( m, m, a, lda, s, dum(1), dum(1),
483 $ dum(1), dum(1), -1, ierr )
486 CALL
cungbr(
'P', m, m, m, a, n, dum(1),
488 lwork_cungbr_p=dum(1)
490 CALL
cungbr(
'Q', m, m, m, a, n, dum(1),
492 lwork_cungbr_q=dum(1)
493 IF( n.GE.mnthr )
THEN
498 maxwrk = m + lwork_cgelqf
499 maxwrk = max( maxwrk, 2*m+lwork_cgebrd )
500 IF( wntuo .OR. wntuas )
501 $ maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
503 ELSE IF( wntvo .AND. wntun )
THEN
507 wrkbl = m + lwork_cgelqf
508 wrkbl = max( wrkbl, m+lwork_cunglq_m )
509 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
510 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
511 maxwrk = max( m*m+wrkbl, m*m+m*n )
513 ELSE IF( wntvo .AND. wntuas )
THEN
518 wrkbl = m + lwork_cgelqf
519 wrkbl = max( wrkbl, m+lwork_cunglq_m )
520 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
521 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
522 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
523 maxwrk = max( m*m+wrkbl, m*m+m*n )
525 ELSE IF( wntvs .AND. wntun )
THEN
529 wrkbl = m + lwork_cgelqf
530 wrkbl = max( wrkbl, m+lwork_cunglq_m )
531 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
532 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
535 ELSE IF( wntvs .AND. wntuo )
THEN
539 wrkbl = m + lwork_cgelqf
540 wrkbl = max( wrkbl, m+lwork_cunglq_m )
541 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
542 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
543 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
544 maxwrk = 2*m*m + wrkbl
546 ELSE IF( wntvs .AND. wntuas )
THEN
551 wrkbl = m + lwork_cgelqf
552 wrkbl = max( wrkbl, m+lwork_cunglq_m )
553 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
554 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
555 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
558 ELSE IF( wntva .AND. wntun )
THEN
562 wrkbl = m + lwork_cgelqf
563 wrkbl = max( wrkbl, m+lwork_cunglq_n )
564 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
565 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
568 ELSE IF( wntva .AND. wntuo )
THEN
572 wrkbl = m + lwork_cgelqf
573 wrkbl = max( wrkbl, m+lwork_cunglq_n )
574 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
575 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
576 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
577 maxwrk = 2*m*m + wrkbl
579 ELSE IF( wntva .AND. wntuas )
THEN
584 wrkbl = m + lwork_cgelqf
585 wrkbl = max( wrkbl, m+lwork_cunglq_n )
586 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
587 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
588 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
596 CALL
cgebrd( m, n, a, lda, s, dum(1), dum(1),
597 $ dum(1), dum(1), -1, ierr )
599 maxwrk = 2*m + lwork_cgebrd
600 IF( wntvs .OR. wntvo )
THEN
602 CALL
cungbr(
'P', m, n, m, a, n, dum(1),
604 lwork_cungbr_p=dum(1)
605 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
608 CALL
cungbr(
'P', n, n, m, a, n, dum(1),
610 lwork_cungbr_p=dum(1)
611 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
613 IF( .NOT.wntun )
THEN
614 maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
619 maxwrk = max( minwrk, maxwrk )
622 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
628 CALL
xerbla(
'CGESVD', -info )
630 ELSE IF( lquery )
THEN
636 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
643 smlnum = sqrt(
slamch(
'S' ) ) / eps
644 bignum = one / smlnum
648 anrm =
clange(
'M', m, n, a, lda, dum )
650 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
652 CALL
clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
653 ELSE IF( anrm.GT.bignum )
THEN
655 CALL
clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
664 IF( m.GE.mnthr )
THEN
678 CALL
cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
679 $ lwork-iwork+1, ierr )
683 CALL
claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
694 CALL
cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
695 $ work( itaup ), work( iwork ), lwork-iwork+1,
698 IF( wntvo .OR. wntvas )
THEN
704 CALL
cungbr(
'P', n, n, n, a, lda, work( itaup ),
705 $ work( iwork ), lwork-iwork+1, ierr )
715 CALL
cbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
716 $ cdum, 1, cdum, 1, rwork( irwork ), info )
721 $ CALL
clacpy(
'F', n, n, a, lda, vt, ldvt )
723 ELSE IF( wntuo .AND. wntvn )
THEN
729 IF( lwork.GE.n*n+3*n )
THEN
734 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
740 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
750 ldwrku = ( lwork-n*n ) / n
760 CALL
cgeqrf( m, n, a, lda, work( itau ),
761 $ work( iwork ), lwork-iwork+1, ierr )
765 CALL
clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
766 CALL
claset(
'L', n-1, n-1, czero, czero,
767 $ work( ir+1 ), ldwrkr )
773 CALL
cungqr( m, n, n, a, lda, work( itau ),
774 $ work( iwork ), lwork-iwork+1, ierr )
784 CALL
cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
785 $ work( itauq ), work( itaup ),
786 $ work( iwork ), lwork-iwork+1, ierr )
792 CALL
cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
793 $ work( itauq ), work( iwork ),
794 $ lwork-iwork+1, ierr )
802 CALL
cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
803 $ work( ir ), ldwrkr, cdum, 1,
804 $ rwork( irwork ), info )
812 DO 10 i = 1, m, ldwrku
813 chunk = min( m-i+1, ldwrku )
814 CALL
cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
815 $ lda, work( ir ), ldwrkr, czero,
816 $ work( iu ), ldwrku )
817 CALL
clacpy(
'F', chunk, n, work( iu ), ldwrku,
834 CALL
cgebrd( m, n, a, lda, s, rwork( ie ),
835 $ work( itauq ), work( itaup ),
836 $ work( iwork ), lwork-iwork+1, ierr )
842 CALL
cungbr(
'Q', m, n, n, a, lda, work( itauq ),
843 $ work( iwork ), lwork-iwork+1, ierr )
851 CALL
cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
852 $ a, lda, cdum, 1, rwork( irwork ), info )
856 ELSE IF( wntuo .AND. wntvas )
THEN
862 IF( lwork.GE.n*n+3*n )
THEN
867 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
873 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
883 ldwrku = ( lwork-n*n ) / n
893 CALL
cgeqrf( m, n, a, lda, work( itau ),
894 $ work( iwork ), lwork-iwork+1, ierr )
898 CALL
clacpy(
'U', n, n, a, lda, vt, ldvt )
900 $ CALL
claset(
'L', n-1, n-1, czero, czero,
907 CALL
cungqr( m, n, n, a, lda, work( itau ),
908 $ work( iwork ), lwork-iwork+1, ierr )
918 CALL
cgebrd( n, n, vt, ldvt, s, rwork( ie ),
919 $ work( itauq ), work( itaup ),
920 $ work( iwork ), lwork-iwork+1, ierr )
921 CALL
clacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
927 CALL
cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
928 $ work( itauq ), work( iwork ),
929 $ lwork-iwork+1, ierr )
935 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
936 $ work( iwork ), lwork-iwork+1, ierr )
945 CALL
cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
946 $ ldvt, work( ir ), ldwrkr, cdum, 1,
947 $ rwork( irwork ), info )
955 DO 20 i = 1, m, ldwrku
956 chunk = min( m-i+1, ldwrku )
957 CALL
cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
958 $ lda, work( ir ), ldwrkr, czero,
959 $ work( iu ), ldwrku )
960 CALL
clacpy(
'F', chunk, n, work( iu ), ldwrku,
975 CALL
cgeqrf( m, n, a, lda, work( itau ),
976 $ work( iwork ), lwork-iwork+1, ierr )
980 CALL
clacpy(
'U', n, n, a, lda, vt, ldvt )
982 $ CALL
claset(
'L', n-1, n-1, czero, czero,
989 CALL
cungqr( m, n, n, a, lda, work( itau ),
990 $ work( iwork ), lwork-iwork+1, ierr )
1000 CALL
cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1001 $ work( itauq ), work( itaup ),
1002 $ work( iwork ), lwork-iwork+1, ierr )
1008 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1009 $ work( itauq ), a, lda, work( iwork ),
1010 $ lwork-iwork+1, ierr )
1016 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1017 $ work( iwork ), lwork-iwork+1, ierr )
1026 CALL
cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1027 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1032 ELSE IF( wntus )
THEN
1040 IF( lwork.GE.n*n+3*n )
THEN
1045 IF( lwork.GE.wrkbl+lda*n )
THEN
1056 itau = ir + ldwrkr*n
1063 CALL
cgeqrf( m, n, a, lda, work( itau ),
1064 $ work( iwork ), lwork-iwork+1, ierr )
1068 CALL
clacpy(
'U', n, n, a, lda, work( ir ),
1070 CALL
claset(
'L', n-1, n-1, czero, czero,
1071 $ work( ir+1 ), ldwrkr )
1077 CALL
cungqr( m, n, n, a, lda, work( itau ),
1078 $ work( iwork ), lwork-iwork+1, ierr )
1088 CALL
cgebrd( n, n, work( ir ), ldwrkr, s,
1089 $ rwork( ie ), work( itauq ),
1090 $ work( itaup ), work( iwork ),
1091 $ lwork-iwork+1, ierr )
1097 CALL
cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1098 $ work( itauq ), work( iwork ),
1099 $ lwork-iwork+1, ierr )
1107 CALL
cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1108 $ 1, work( ir ), ldwrkr, cdum, 1,
1109 $ rwork( irwork ), info )
1116 CALL
cgemm(
'N',
'N', m, n, n, cone, a, lda,
1117 $ work( ir ), ldwrkr, czero, u, ldu )
1130 CALL
cgeqrf( m, n, a, lda, work( itau ),
1131 $ work( iwork ), lwork-iwork+1, ierr )
1132 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1138 CALL
cungqr( m, n, n, u, ldu, work( itau ),
1139 $ work( iwork ), lwork-iwork+1, ierr )
1147 CALL
claset(
'L', n-1, n-1, czero, czero,
1154 CALL
cgebrd( n, n, a, lda, s, rwork( ie ),
1155 $ work( itauq ), work( itaup ),
1156 $ work( iwork ), lwork-iwork+1, ierr )
1162 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1163 $ work( itauq ), u, ldu, work( iwork ),
1164 $ lwork-iwork+1, ierr )
1172 CALL
cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1173 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1178 ELSE IF( wntvo )
THEN
1184 IF( lwork.GE.2*n*n+3*n )
THEN
1189 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1196 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1211 itau = ir + ldwrkr*n
1218 CALL
cgeqrf( m, n, a, lda, work( itau ),
1219 $ work( iwork ), lwork-iwork+1, ierr )
1223 CALL
clacpy(
'U', n, n, a, lda, work( iu ),
1225 CALL
claset(
'L', n-1, n-1, czero, czero,
1226 $ work( iu+1 ), ldwrku )
1232 CALL
cungqr( m, n, n, a, lda, work( itau ),
1233 $ work( iwork ), lwork-iwork+1, ierr )
1245 CALL
cgebrd( n, n, work( iu ), ldwrku, s,
1246 $ rwork( ie ), work( itauq ),
1247 $ work( itaup ), work( iwork ),
1248 $ lwork-iwork+1, ierr )
1249 CALL
clacpy(
'U', n, n, work( iu ), ldwrku,
1250 $ work( ir ), ldwrkr )
1256 CALL
cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1257 $ work( itauq ), work( iwork ),
1258 $ lwork-iwork+1, ierr )
1265 CALL
cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1266 $ work( itaup ), work( iwork ),
1267 $ lwork-iwork+1, ierr )
1276 CALL
cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1277 $ work( ir ), ldwrkr, work( iu ),
1278 $ ldwrku, cdum, 1, rwork( irwork ),
1286 CALL
cgemm(
'N',
'N', m, n, n, cone, a, lda,
1287 $ work( iu ), ldwrku, czero, u, ldu )
1293 CALL
clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1307 CALL
cgeqrf( m, n, a, lda, work( itau ),
1308 $ work( iwork ), lwork-iwork+1, ierr )
1309 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1315 CALL
cungqr( m, n, n, u, ldu, work( itau ),
1316 $ work( iwork ), lwork-iwork+1, ierr )
1324 CALL
claset(
'L', n-1, n-1, czero, czero,
1331 CALL
cgebrd( n, n, a, lda, s, rwork( ie ),
1332 $ work( itauq ), work( itaup ),
1333 $ work( iwork ), lwork-iwork+1, ierr )
1339 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1340 $ work( itauq ), u, ldu, work( iwork ),
1341 $ lwork-iwork+1, ierr )
1347 CALL
cungbr(
'P', n, n, n, a, lda, work( itaup ),
1348 $ work( iwork ), lwork-iwork+1, ierr )
1357 CALL
cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1358 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1363 ELSE IF( wntvas )
THEN
1370 IF( lwork.GE.n*n+3*n )
THEN
1375 IF( lwork.GE.wrkbl+lda*n )
THEN
1386 itau = iu + ldwrku*n
1393 CALL
cgeqrf( m, n, a, lda, work( itau ),
1394 $ work( iwork ), lwork-iwork+1, ierr )
1398 CALL
clacpy(
'U', n, n, a, lda, work( iu ),
1400 CALL
claset(
'L', n-1, n-1, czero, czero,
1401 $ work( iu+1 ), ldwrku )
1407 CALL
cungqr( m, n, n, a, lda, work( itau ),
1408 $ work( iwork ), lwork-iwork+1, ierr )
1418 CALL
cgebrd( n, n, work( iu ), ldwrku, s,
1419 $ rwork( ie ), work( itauq ),
1420 $ work( itaup ), work( iwork ),
1421 $ lwork-iwork+1, ierr )
1422 CALL
clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1429 CALL
cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1430 $ work( itauq ), work( iwork ),
1431 $ lwork-iwork+1, ierr )
1438 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1439 $ work( iwork ), lwork-iwork+1, ierr )
1448 CALL
cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1449 $ ldvt, work( iu ), ldwrku, cdum, 1,
1450 $ rwork( irwork ), info )
1457 CALL
cgemm(
'N',
'N', m, n, n, cone, a, lda,
1458 $ work( iu ), ldwrku, czero, u, ldu )
1471 CALL
cgeqrf( m, n, a, lda, work( itau ),
1472 $ work( iwork ), lwork-iwork+1, ierr )
1473 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1479 CALL
cungqr( m, n, n, u, ldu, work( itau ),
1480 $ work( iwork ), lwork-iwork+1, ierr )
1484 CALL
clacpy(
'U', n, n, a, lda, vt, ldvt )
1486 $ CALL
claset(
'L', n-1, n-1, czero, czero,
1487 $ vt( 2, 1 ), ldvt )
1497 CALL
cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1498 $ work( itauq ), work( itaup ),
1499 $ work( iwork ), lwork-iwork+1, ierr )
1506 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1507 $ work( itauq ), u, ldu, work( iwork ),
1508 $ lwork-iwork+1, ierr )
1514 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1515 $ work( iwork ), lwork-iwork+1, ierr )
1524 CALL
cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1525 $ ldvt, u, ldu, cdum, 1,
1526 $ rwork( irwork ), info )
1532 ELSE IF( wntua )
THEN
1540 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1545 IF( lwork.GE.wrkbl+lda*n )
THEN
1556 itau = ir + ldwrkr*n
1563 CALL
cgeqrf( m, n, a, lda, work( itau ),
1564 $ work( iwork ), lwork-iwork+1, ierr )
1565 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1569 CALL
clacpy(
'U', n, n, a, lda, work( ir ),
1571 CALL
claset(
'L', n-1, n-1, czero, czero,
1572 $ work( ir+1 ), ldwrkr )
1578 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1579 $ work( iwork ), lwork-iwork+1, ierr )
1589 CALL
cgebrd( n, n, work( ir ), ldwrkr, s,
1590 $ rwork( ie ), work( itauq ),
1591 $ work( itaup ), work( iwork ),
1592 $ lwork-iwork+1, ierr )
1598 CALL
cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1599 $ work( itauq ), work( iwork ),
1600 $ lwork-iwork+1, ierr )
1608 CALL
cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1609 $ 1, work( ir ), ldwrkr, cdum, 1,
1610 $ rwork( irwork ), info )
1617 CALL
cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1618 $ work( ir ), ldwrkr, czero, a, lda )
1622 CALL
clacpy(
'F', m, n, a, lda, u, ldu )
1635 CALL
cgeqrf( m, n, a, lda, work( itau ),
1636 $ work( iwork ), lwork-iwork+1, ierr )
1637 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1643 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1644 $ work( iwork ), lwork-iwork+1, ierr )
1652 CALL
claset(
'L', n-1, n-1, czero, czero,
1659 CALL
cgebrd( n, n, a, lda, s, rwork( ie ),
1660 $ work( itauq ), work( itaup ),
1661 $ work( iwork ), lwork-iwork+1, ierr )
1668 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1669 $ work( itauq ), u, ldu, work( iwork ),
1670 $ lwork-iwork+1, ierr )
1678 CALL
cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1679 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1684 ELSE IF( wntvo )
THEN
1690 IF( lwork.GE.2*n*n+max( n+m, 3*n ) )
THEN
1695 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1702 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1717 itau = ir + ldwrkr*n
1724 CALL
cgeqrf( m, n, a, lda, work( itau ),
1725 $ work( iwork ), lwork-iwork+1, ierr )
1726 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1732 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1733 $ work( iwork ), lwork-iwork+1, ierr )
1737 CALL
clacpy(
'U', n, n, a, lda, work( iu ),
1739 CALL
claset(
'L', n-1, n-1, czero, czero,
1740 $ work( iu+1 ), ldwrku )
1752 CALL
cgebrd( n, n, work( iu ), ldwrku, s,
1753 $ rwork( ie ), work( itauq ),
1754 $ work( itaup ), work( iwork ),
1755 $ lwork-iwork+1, ierr )
1756 CALL
clacpy(
'U', n, n, work( iu ), ldwrku,
1757 $ work( ir ), ldwrkr )
1763 CALL
cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1764 $ work( itauq ), work( iwork ),
1765 $ lwork-iwork+1, ierr )
1772 CALL
cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1773 $ work( itaup ), work( iwork ),
1774 $ lwork-iwork+1, ierr )
1783 CALL
cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1784 $ work( ir ), ldwrkr, work( iu ),
1785 $ ldwrku, cdum, 1, rwork( irwork ),
1793 CALL
cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1794 $ work( iu ), ldwrku, czero, a, lda )
1798 CALL
clacpy(
'F', m, n, a, lda, u, ldu )
1802 CALL
clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1816 CALL
cgeqrf( m, n, a, lda, work( itau ),
1817 $ work( iwork ), lwork-iwork+1, ierr )
1818 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1824 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1825 $ work( iwork ), lwork-iwork+1, ierr )
1833 CALL
claset(
'L', n-1, n-1, czero, czero,
1840 CALL
cgebrd( n, n, a, lda, s, rwork( ie ),
1841 $ work( itauq ), work( itaup ),
1842 $ work( iwork ), lwork-iwork+1, ierr )
1849 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1850 $ work( itauq ), u, ldu, work( iwork ),
1851 $ lwork-iwork+1, ierr )
1857 CALL
cungbr(
'P', n, n, n, a, lda, work( itaup ),
1858 $ work( iwork ), lwork-iwork+1, ierr )
1867 CALL
cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1868 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1873 ELSE IF( wntvas )
THEN
1880 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1885 IF( lwork.GE.wrkbl+lda*n )
THEN
1896 itau = iu + ldwrku*n
1903 CALL
cgeqrf( m, n, a, lda, work( itau ),
1904 $ work( iwork ), lwork-iwork+1, ierr )
1905 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1911 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1912 $ work( iwork ), lwork-iwork+1, ierr )
1916 CALL
clacpy(
'U', n, n, a, lda, work( iu ),
1918 CALL
claset(
'L', n-1, n-1, czero, czero,
1919 $ work( iu+1 ), ldwrku )
1929 CALL
cgebrd( n, n, work( iu ), ldwrku, s,
1930 $ rwork( ie ), work( itauq ),
1931 $ work( itaup ), work( iwork ),
1932 $ lwork-iwork+1, ierr )
1933 CALL
clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1940 CALL
cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1941 $ work( itauq ), work( iwork ),
1942 $ lwork-iwork+1, ierr )
1949 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1950 $ work( iwork ), lwork-iwork+1, ierr )
1959 CALL
cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1960 $ ldvt, work( iu ), ldwrku, cdum, 1,
1961 $ rwork( irwork ), info )
1968 CALL
cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1969 $ work( iu ), ldwrku, czero, a, lda )
1973 CALL
clacpy(
'F', m, n, a, lda, u, ldu )
1986 CALL
cgeqrf( m, n, a, lda, work( itau ),
1987 $ work( iwork ), lwork-iwork+1, ierr )
1988 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
1994 CALL
cungqr( m, m, n, u, ldu, work( itau ),
1995 $ work( iwork ), lwork-iwork+1, ierr )
1999 CALL
clacpy(
'U', n, n, a, lda, vt, ldvt )
2001 $ CALL
claset(
'L', n-1, n-1, czero, czero,
2002 $ vt( 2, 1 ), ldvt )
2012 CALL
cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2013 $ work( itauq ), work( itaup ),
2014 $ work( iwork ), lwork-iwork+1, ierr )
2021 CALL
cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2022 $ work( itauq ), u, ldu, work( iwork ),
2023 $ lwork-iwork+1, ierr )
2029 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2030 $ work( iwork ), lwork-iwork+1, ierr )
2039 CALL
cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
2040 $ ldvt, u, ldu, cdum, 1,
2041 $ rwork( irwork ), info )
2065 CALL
cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2066 $ work( itaup ), work( iwork ), lwork-iwork+1,
2075 CALL
clacpy(
'L', m, n, a, lda, u, ldu )
2080 CALL
cungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2081 $ work( iwork ), lwork-iwork+1, ierr )
2090 CALL
clacpy(
'U', n, n, a, lda, vt, ldvt )
2091 CALL
cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2092 $ work( iwork ), lwork-iwork+1, ierr )
2101 CALL
cungbr(
'Q', m, n, n, a, lda, work( itauq ),
2102 $ work( iwork ), lwork-iwork+1, ierr )
2111 CALL
cungbr(
'P', n, n, n, a, lda, work( itaup ),
2112 $ work( iwork ), lwork-iwork+1, ierr )
2115 IF( wntuas .OR. wntuo )
2119 IF( wntvas .OR. wntvo )
2123 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2131 CALL
cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2132 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2134 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2142 CALL
cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2143 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2153 CALL
cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2154 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2166 IF( n.GE.mnthr )
THEN
2180 CALL
cgelqf( m, n, a, lda, work( itau ), work( iwork ),
2181 $ lwork-iwork+1, ierr )
2185 CALL
claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2196 CALL
cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2197 $ work( itaup ), work( iwork ), lwork-iwork+1,
2199 IF( wntuo .OR. wntuas )
THEN
2205 CALL
cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2206 $ work( iwork ), lwork-iwork+1, ierr )
2210 IF( wntuo .OR. wntuas )
2218 CALL
cbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2219 $ a, lda, cdum, 1, rwork( irwork ), info )
2224 $ CALL
clacpy(
'F', m, m, a, lda, u, ldu )
2226 ELSE IF( wntvo .AND. wntun )
THEN
2232 IF( lwork.GE.m*m+3*m )
THEN
2237 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2244 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2256 chunk = ( lwork-m*m ) / m
2259 itau = ir + ldwrkr*m
2266 CALL
cgelqf( m, n, a, lda, work( itau ),
2267 $ work( iwork ), lwork-iwork+1, ierr )
2271 CALL
clacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2272 CALL
claset(
'U', m-1, m-1, czero, czero,
2273 $ work( ir+ldwrkr ), ldwrkr )
2279 CALL
cunglq( m, n, m, a, lda, work( itau ),
2280 $ work( iwork ), lwork-iwork+1, ierr )
2290 CALL
cgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2291 $ work( itauq ), work( itaup ),
2292 $ work( iwork ), lwork-iwork+1, ierr )
2298 CALL
cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2299 $ work( itaup ), work( iwork ),
2300 $ lwork-iwork+1, ierr )
2308 CALL
cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2309 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2310 $ rwork( irwork ), info )
2318 DO 30 i = 1, n, chunk
2319 blk = min( n-i+1, chunk )
2320 CALL
cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2321 $ ldwrkr, a( 1, i ), lda, czero,
2322 $ work( iu ), ldwrku )
2323 CALL
clacpy(
'F', m, blk, work( iu ), ldwrku,
2340 CALL
cgebrd( m, n, a, lda, s, rwork( ie ),
2341 $ work( itauq ), work( itaup ),
2342 $ work( iwork ), lwork-iwork+1, ierr )
2348 CALL
cungbr(
'P', m, n, m, a, lda, work( itaup ),
2349 $ work( iwork ), lwork-iwork+1, ierr )
2357 CALL
cbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2358 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2362 ELSE IF( wntvo .AND. wntuas )
THEN
2368 IF( lwork.GE.m*m+3*m )
THEN
2373 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2380 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2392 chunk = ( lwork-m*m ) / m
2395 itau = ir + ldwrkr*m
2402 CALL
cgelqf( m, n, a, lda, work( itau ),
2403 $ work( iwork ), lwork-iwork+1, ierr )
2407 CALL
clacpy(
'L', m, m, a, lda, u, ldu )
2408 CALL
claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2415 CALL
cunglq( m, n, m, a, lda, work( itau ),
2416 $ work( iwork ), lwork-iwork+1, ierr )
2426 CALL
cgebrd( m, m, u, ldu, s, rwork( ie ),
2427 $ work( itauq ), work( itaup ),
2428 $ work( iwork ), lwork-iwork+1, ierr )
2429 CALL
clacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2435 CALL
cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2436 $ work( itaup ), work( iwork ),
2437 $ lwork-iwork+1, ierr )
2443 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2444 $ work( iwork ), lwork-iwork+1, ierr )
2453 CALL
cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2454 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2455 $ rwork( irwork ), info )
2463 DO 40 i = 1, n, chunk
2464 blk = min( n-i+1, chunk )
2465 CALL
cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2466 $ ldwrkr, a( 1, i ), lda, czero,
2467 $ work( iu ), ldwrku )
2468 CALL
clacpy(
'F', m, blk, work( iu ), ldwrku,
2483 CALL
cgelqf( m, n, a, lda, work( itau ),
2484 $ work( iwork ), lwork-iwork+1, ierr )
2488 CALL
clacpy(
'L', m, m, a, lda, u, ldu )
2489 CALL
claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2496 CALL
cunglq( m, n, m, a, lda, work( itau ),
2497 $ work( iwork ), lwork-iwork+1, ierr )
2507 CALL
cgebrd( m, m, u, ldu, s, rwork( ie ),
2508 $ work( itauq ), work( itaup ),
2509 $ work( iwork ), lwork-iwork+1, ierr )
2515 CALL
cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2516 $ work( itaup ), a, lda, work( iwork ),
2517 $ lwork-iwork+1, ierr )
2523 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2524 $ work( iwork ), lwork-iwork+1, ierr )
2533 CALL
cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a, lda,
2534 $ u, ldu, cdum, 1, rwork( irwork ), info )
2538 ELSE IF( wntvs )
THEN
2546 IF( lwork.GE.m*m+3*m )
THEN
2551 IF( lwork.GE.wrkbl+lda*m )
THEN
2562 itau = ir + ldwrkr*m
2569 CALL
cgelqf( m, n, a, lda, work( itau ),
2570 $ work( iwork ), lwork-iwork+1, ierr )
2574 CALL
clacpy(
'L', m, m, a, lda, work( ir ),
2576 CALL
claset(
'U', m-1, m-1, czero, czero,
2577 $ work( ir+ldwrkr ), ldwrkr )
2583 CALL
cunglq( m, n, m, a, lda, work( itau ),
2584 $ work( iwork ), lwork-iwork+1, ierr )
2594 CALL
cgebrd( m, m, work( ir ), ldwrkr, s,
2595 $ rwork( ie ), work( itauq ),
2596 $ work( itaup ), work( iwork ),
2597 $ lwork-iwork+1, ierr )
2604 CALL
cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2605 $ work( itaup ), work( iwork ),
2606 $ lwork-iwork+1, ierr )
2614 CALL
cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2615 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2616 $ rwork( irwork ), info )
2623 CALL
cgemm(
'N',
'N', m, n, m, cone, work( ir ),
2624 $ ldwrkr, a, lda, czero, vt, ldvt )
2637 CALL
cgelqf( m, n, a, lda, work( itau ),
2638 $ work( iwork ), lwork-iwork+1, ierr )
2642 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
2648 CALL
cunglq( m, n, m, vt, ldvt, work( itau ),
2649 $ work( iwork ), lwork-iwork+1, ierr )
2657 CALL
claset(
'U', m-1, m-1, czero, czero,
2664 CALL
cgebrd( m, m, a, lda, s, rwork( ie ),
2665 $ work( itauq ), work( itaup ),
2666 $ work( iwork ), lwork-iwork+1, ierr )
2672 CALL
cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2673 $ work( itaup ), vt, ldvt,
2674 $ work( iwork ), lwork-iwork+1, ierr )
2682 CALL
cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
2683 $ ldvt, cdum, 1, cdum, 1,
2684 $ rwork( irwork ), info )
2688 ELSE IF( wntuo )
THEN
2694 IF( lwork.GE.2*m*m+3*m )
THEN
2699 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2706 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2721 itau = ir + ldwrkr*m
2728 CALL
cgelqf( m, n, a, lda, work( itau ),
2729 $ work( iwork ), lwork-iwork+1, ierr )
2733 CALL
clacpy(
'L', m, m, a, lda, work( iu ),
2735 CALL
claset(
'U', m-1, m-1, czero, czero,
2736 $ work( iu+ldwrku ), ldwrku )
2742 CALL
cunglq( m, n, m, a, lda, work( itau ),
2743 $ work( iwork ), lwork-iwork+1, ierr )
2755 CALL
cgebrd( m, m, work( iu ), ldwrku, s,
2756 $ rwork( ie ), work( itauq ),
2757 $ work( itaup ), work( iwork ),
2758 $ lwork-iwork+1, ierr )
2759 CALL
clacpy(
'L', m, m, work( iu ), ldwrku,
2760 $ work( ir ), ldwrkr )
2767 CALL
cungbr(
'P', m, m, m, work( iu ), ldwrku,
2768 $ work( itaup ), work( iwork ),
2769 $ lwork-iwork+1, ierr )
2775 CALL
cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2776 $ work( itauq ), work( iwork ),
2777 $ lwork-iwork+1, ierr )
2786 CALL
cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2787 $ work( iu ), ldwrku, work( ir ),
2788 $ ldwrkr, cdum, 1, rwork( irwork ),
2796 CALL
cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2797 $ ldwrku, a, lda, czero, vt, ldvt )
2803 CALL
clacpy(
'F', m, m, work( ir ), ldwrkr, a,
2817 CALL
cgelqf( m, n, a, lda, work( itau ),
2818 $ work( iwork ), lwork-iwork+1, ierr )
2819 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
2825 CALL
cunglq( m, n, m, vt, ldvt, work( itau ),
2826 $ work( iwork ), lwork-iwork+1, ierr )
2834 CALL
claset(
'U', m-1, m-1, czero, czero,
2841 CALL
cgebrd( m, m, a, lda, s, rwork( ie ),
2842 $ work( itauq ), work( itaup ),
2843 $ work( iwork ), lwork-iwork+1, ierr )
2849 CALL
cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2850 $ work( itaup ), vt, ldvt,
2851 $ work( iwork ), lwork-iwork+1, ierr )
2857 CALL
cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2858 $ work( iwork ), lwork-iwork+1, ierr )
2867 CALL
cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
2868 $ ldvt, a, lda, cdum, 1,
2869 $ rwork( irwork ), info )
2873 ELSE IF( wntuas )
THEN
2880 IF( lwork.GE.m*m+3*m )
THEN
2885 IF( lwork.GE.wrkbl+lda*m )
THEN
2896 itau = iu + ldwrku*m
2903 CALL
cgelqf( m, n, a, lda, work( itau ),
2904 $ work( iwork ), lwork-iwork+1, ierr )
2908 CALL
clacpy(
'L', m, m, a, lda, work( iu ),
2910 CALL
claset(
'U', m-1, m-1, czero, czero,
2911 $ work( iu+ldwrku ), ldwrku )
2917 CALL
cunglq( m, n, m, a, lda, work( itau ),
2918 $ work( iwork ), lwork-iwork+1, ierr )
2928 CALL
cgebrd( m, m, work( iu ), ldwrku, s,
2929 $ rwork( ie ), work( itauq ),
2930 $ work( itaup ), work( iwork ),
2931 $ lwork-iwork+1, ierr )
2932 CALL
clacpy(
'L', m, m, work( iu ), ldwrku, u,
2940 CALL
cungbr(
'P', m, m, m, work( iu ), ldwrku,
2941 $ work( itaup ), work( iwork ),
2942 $ lwork-iwork+1, ierr )
2948 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2949 $ work( iwork ), lwork-iwork+1, ierr )
2958 CALL
cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2959 $ work( iu ), ldwrku, u, ldu, cdum, 1,
2960 $ rwork( irwork ), info )
2967 CALL
cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2968 $ ldwrku, a, lda, czero, vt, ldvt )
2981 CALL
cgelqf( m, n, a, lda, work( itau ),
2982 $ work( iwork ), lwork-iwork+1, ierr )
2983 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
2989 CALL
cunglq( m, n, m, vt, ldvt, work( itau ),
2990 $ work( iwork ), lwork-iwork+1, ierr )
2994 CALL
clacpy(
'L', m, m, a, lda, u, ldu )
2995 CALL
claset(
'U', m-1, m-1, czero, czero,
3006 CALL
cgebrd( m, m, u, ldu, s, rwork( ie ),
3007 $ work( itauq ), work( itaup ),
3008 $ work( iwork ), lwork-iwork+1, ierr )
3015 CALL
cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3016 $ work( itaup ), vt, ldvt,
3017 $ work( iwork ), lwork-iwork+1, ierr )
3023 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3024 $ work( iwork ), lwork-iwork+1, ierr )
3033 CALL
cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3034 $ ldvt, u, ldu, cdum, 1,
3035 $ rwork( irwork ), info )
3041 ELSE IF( wntva )
THEN
3049 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3054 IF( lwork.GE.wrkbl+lda*m )
THEN
3065 itau = ir + ldwrkr*m
3072 CALL
cgelqf( m, n, a, lda, work( itau ),
3073 $ work( iwork ), lwork-iwork+1, ierr )
3074 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3078 CALL
clacpy(
'L', m, m, a, lda, work( ir ),
3080 CALL
claset(
'U', m-1, m-1, czero, czero,
3081 $ work( ir+ldwrkr ), ldwrkr )
3087 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3088 $ work( iwork ), lwork-iwork+1, ierr )
3098 CALL
cgebrd( m, m, work( ir ), ldwrkr, s,
3099 $ rwork( ie ), work( itauq ),
3100 $ work( itaup ), work( iwork ),
3101 $ lwork-iwork+1, ierr )
3108 CALL
cungbr(
'P', m, m, m, work( ir ), ldwrkr,
3109 $ work( itaup ), work( iwork ),
3110 $ lwork-iwork+1, ierr )
3118 CALL
cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3119 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3120 $ rwork( irwork ), info )
3127 CALL
cgemm(
'N',
'N', m, n, m, cone, work( ir ),
3128 $ ldwrkr, vt, ldvt, czero, a, lda )
3132 CALL
clacpy(
'F', m, n, a, lda, vt, ldvt )
3145 CALL
cgelqf( m, n, a, lda, work( itau ),
3146 $ work( iwork ), lwork-iwork+1, ierr )
3147 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3153 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3154 $ work( iwork ), lwork-iwork+1, ierr )
3162 CALL
claset(
'U', m-1, m-1, czero, czero,
3169 CALL
cgebrd( m, m, a, lda, s, rwork( ie ),
3170 $ work( itauq ), work( itaup ),
3171 $ work( iwork ), lwork-iwork+1, ierr )
3178 CALL
cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3179 $ work( itaup ), vt, ldvt,
3180 $ work( iwork ), lwork-iwork+1, ierr )
3188 CALL
cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
3189 $ ldvt, cdum, 1, cdum, 1,
3190 $ rwork( irwork ), info )
3194 ELSE IF( wntuo )
THEN
3200 IF( lwork.GE.2*m*m+max( n+m, 3*m ) )
THEN
3205 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3212 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3227 itau = ir + ldwrkr*m
3234 CALL
cgelqf( m, n, a, lda, work( itau ),
3235 $ work( iwork ), lwork-iwork+1, ierr )
3236 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3242 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3243 $ work( iwork ), lwork-iwork+1, ierr )
3247 CALL
clacpy(
'L', m, m, a, lda, work( iu ),
3249 CALL
claset(
'U', m-1, m-1, czero, czero,
3250 $ work( iu+ldwrku ), ldwrku )
3262 CALL
cgebrd( m, m, work( iu ), ldwrku, s,
3263 $ rwork( ie ), work( itauq ),
3264 $ work( itaup ), work( iwork ),
3265 $ lwork-iwork+1, ierr )
3266 CALL
clacpy(
'L', m, m, work( iu ), ldwrku,
3267 $ work( ir ), ldwrkr )
3274 CALL
cungbr(
'P', m, m, m, work( iu ), ldwrku,
3275 $ work( itaup ), work( iwork ),
3276 $ lwork-iwork+1, ierr )
3282 CALL
cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3283 $ work( itauq ), work( iwork ),
3284 $ lwork-iwork+1, ierr )
3293 CALL
cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3294 $ work( iu ), ldwrku, work( ir ),
3295 $ ldwrkr, cdum, 1, rwork( irwork ),
3303 CALL
cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3304 $ ldwrku, vt, ldvt, czero, a, lda )
3308 CALL
clacpy(
'F', m, n, a, lda, vt, ldvt )
3312 CALL
clacpy(
'F', m, m, work( ir ), ldwrkr, a,
3326 CALL
cgelqf( m, n, a, lda, work( itau ),
3327 $ work( iwork ), lwork-iwork+1, ierr )
3328 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3334 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3335 $ work( iwork ), lwork-iwork+1, ierr )
3343 CALL
claset(
'U', m-1, m-1, czero, czero,
3350 CALL
cgebrd( m, m, a, lda, s, rwork( ie ),
3351 $ work( itauq ), work( itaup ),
3352 $ work( iwork ), lwork-iwork+1, ierr )
3359 CALL
cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3360 $ work( itaup ), vt, ldvt,
3361 $ work( iwork ), lwork-iwork+1, ierr )
3367 CALL
cungbr(
'Q', m, m, m, a, lda, work( itauq ),
3368 $ work( iwork ), lwork-iwork+1, ierr )
3377 CALL
cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3378 $ ldvt, a, lda, cdum, 1,
3379 $ rwork( irwork ), info )
3383 ELSE IF( wntuas )
THEN
3390 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3395 IF( lwork.GE.wrkbl+lda*m )
THEN
3406 itau = iu + ldwrku*m
3413 CALL
cgelqf( m, n, a, lda, work( itau ),
3414 $ work( iwork ), lwork-iwork+1, ierr )
3415 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3421 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3422 $ work( iwork ), lwork-iwork+1, ierr )
3426 CALL
clacpy(
'L', m, m, a, lda, work( iu ),
3428 CALL
claset(
'U', m-1, m-1, czero, czero,
3429 $ work( iu+ldwrku ), ldwrku )
3439 CALL
cgebrd( m, m, work( iu ), ldwrku, s,
3440 $ rwork( ie ), work( itauq ),
3441 $ work( itaup ), work( iwork ),
3442 $ lwork-iwork+1, ierr )
3443 CALL
clacpy(
'L', m, m, work( iu ), ldwrku, u,
3450 CALL
cungbr(
'P', m, m, m, work( iu ), ldwrku,
3451 $ work( itaup ), work( iwork ),
3452 $ lwork-iwork+1, ierr )
3458 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3459 $ work( iwork ), lwork-iwork+1, ierr )
3468 CALL
cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3469 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3470 $ rwork( irwork ), info )
3477 CALL
cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3478 $ ldwrku, vt, ldvt, czero, a, lda )
3482 CALL
clacpy(
'F', m, n, a, lda, vt, ldvt )
3495 CALL
cgelqf( m, n, a, lda, work( itau ),
3496 $ work( iwork ), lwork-iwork+1, ierr )
3497 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3503 CALL
cunglq( n, n, m, vt, ldvt, work( itau ),
3504 $ work( iwork ), lwork-iwork+1, ierr )
3508 CALL
clacpy(
'L', m, m, a, lda, u, ldu )
3509 CALL
claset(
'U', m-1, m-1, czero, czero,
3520 CALL
cgebrd( m, m, u, ldu, s, rwork( ie ),
3521 $ work( itauq ), work( itaup ),
3522 $ work( iwork ), lwork-iwork+1, ierr )
3529 CALL
cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3530 $ work( itaup ), vt, ldvt,
3531 $ work( iwork ), lwork-iwork+1, ierr )
3537 CALL
cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3538 $ work( iwork ), lwork-iwork+1, ierr )
3547 CALL
cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3548 $ ldvt, u, ldu, cdum, 1,
3549 $ rwork( irwork ), info )
3573 CALL
cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3574 $ work( itaup ), work( iwork ), lwork-iwork+1,
3583 CALL
clacpy(
'L', m, m, a, lda, u, ldu )
3584 CALL
cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3585 $ work( iwork ), lwork-iwork+1, ierr )
3594 CALL
clacpy(
'U', m, n, a, lda, vt, ldvt )
3599 CALL
cungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3600 $ work( iwork ), lwork-iwork+1, ierr )
3609 CALL
cungbr(
'Q', m, m, n, a, lda, work( itauq ),
3610 $ work( iwork ), lwork-iwork+1, ierr )
3619 CALL
cungbr(
'P', m, n, m, a, lda, work( itaup ),
3620 $ work( iwork ), lwork-iwork+1, ierr )
3623 IF( wntuas .OR. wntuo )
3627 IF( wntvas .OR. wntvo )
3631 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3639 CALL
cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3640 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3642 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3650 CALL
cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3651 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3661 CALL
cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3662 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3672 IF( iscl.EQ.1 )
THEN
3673 IF( anrm.GT.bignum )
3674 $ CALL
slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3676 IF( info.NE.0 .AND. anrm.GT.bignum )
3677 $ CALL
slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3678 $ rwork( ie ), minmn, ierr )
3679 IF( anrm.LT.smlnum )
3680 $ CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3682 IF( info.NE.0 .AND. anrm.LT.smlnum )
3683 $ CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3684 $ rwork( ie ), minmn, ierr )
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
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...
LOGICAL function lsame(CA, CB)
LSAME
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
REAL function slamch(CMACH)
SLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
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.
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 ...
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
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.
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR