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
266 EXTERNAL lsame, ilaenv, clange, slamch
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, cdum(1), cdum(1), -1, ierr )
327 CALL cungqr( m, n, n, a, lda, cdum(1), cdum(1), -1, ierr )
328 lwork_cungqr_n=cdum(1)
329 CALL cungqr( m, m, n, a, lda, cdum(1), cdum(1), -1, ierr )
330 lwork_cungqr_m=cdum(1)
332 CALL cgebrd( n, n, a, lda, s, dum(1), cdum(1),
333 $ cdum(1), cdum(1), -1, ierr )
336 CALL cungbr(
'P', n, n, n, a, lda, cdum(1),
337 $ cdum(1), -1, ierr )
338 lwork_cungbr_p=cdum(1)
339 CALL cungbr(
'Q', n, n, n, a, lda, cdum(1),
340 $ cdum(1), -1, ierr )
341 lwork_cungbr_q=cdum(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), cdum(1),
448 $ cdum(1), cdum(1), -1, ierr )
450 maxwrk = 2*n + lwork_cgebrd
451 IF( wntus .OR. wntuo )
THEN
452 CALL cungbr(
'Q', m, n, n, a, lda, cdum(1),
453 $ cdum(1), -1, ierr )
454 lwork_cungbr_q=cdum(1)
455 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
458 CALL cungbr(
'Q', m, m, n, a, lda, cdum(1),
459 $ cdum(1), -1, ierr )
460 lwork_cungbr_q=cdum(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, cdum(1), cdum(1), -1, ierr )
477 CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
479 lwork_cunglq_n=cdum(1)
480 CALL cunglq( m, n, m, a, lda, cdum(1), cdum(1), -1, ierr )
481 lwork_cunglq_m=cdum(1)
483 CALL cgebrd( m, m, a, lda, s, dum(1), cdum(1),
484 $ cdum(1), cdum(1), -1, ierr )
487 CALL cungbr(
'P', m, m, m, a, n, cdum(1),
488 $ cdum(1), -1, ierr )
489 lwork_cungbr_p=cdum(1)
491 CALL cungbr(
'Q', m, m, m, a, n, cdum(1),
492 $ cdum(1), -1, ierr )
493 lwork_cungbr_q=cdum(1)
494 IF( n.GE.mnthr )
THEN
499 maxwrk = m + lwork_cgelqf
500 maxwrk = max( maxwrk, 2*m+lwork_cgebrd )
501 IF( wntuo .OR. wntuas )
502 $ maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
504 ELSE IF( wntvo .AND. wntun )
THEN
508 wrkbl = m + lwork_cgelqf
509 wrkbl = max( wrkbl, m+lwork_cunglq_m )
510 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
511 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
512 maxwrk = max( m*m+wrkbl, m*m+m*n )
514 ELSE IF( wntvo .AND. wntuas )
THEN
519 wrkbl = m + lwork_cgelqf
520 wrkbl = max( wrkbl, m+lwork_cunglq_m )
521 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
522 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
523 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
524 maxwrk = max( m*m+wrkbl, m*m+m*n )
526 ELSE IF( wntvs .AND. wntun )
THEN
530 wrkbl = m + lwork_cgelqf
531 wrkbl = max( wrkbl, m+lwork_cunglq_m )
532 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
533 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
536 ELSE IF( wntvs .AND. wntuo )
THEN
540 wrkbl = m + lwork_cgelqf
541 wrkbl = max( wrkbl, m+lwork_cunglq_m )
542 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
543 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
544 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
545 maxwrk = 2*m*m + wrkbl
547 ELSE IF( wntvs .AND. wntuas )
THEN
552 wrkbl = m + lwork_cgelqf
553 wrkbl = max( wrkbl, m+lwork_cunglq_m )
554 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
555 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
556 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
559 ELSE IF( wntva .AND. wntun )
THEN
563 wrkbl = m + lwork_cgelqf
564 wrkbl = max( wrkbl, m+lwork_cunglq_n )
565 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
566 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
569 ELSE IF( wntva .AND. wntuo )
THEN
573 wrkbl = m + lwork_cgelqf
574 wrkbl = max( wrkbl, m+lwork_cunglq_n )
575 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
576 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
577 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
578 maxwrk = 2*m*m + wrkbl
580 ELSE IF( wntva .AND. wntuas )
THEN
585 wrkbl = m + lwork_cgelqf
586 wrkbl = max( wrkbl, m+lwork_cunglq_n )
587 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
588 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
589 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
597 CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
598 $ cdum(1), cdum(1), -1, ierr )
600 maxwrk = 2*m + lwork_cgebrd
601 IF( wntvs .OR. wntvo )
THEN
603 CALL cungbr(
'P', m, n, m, a, n, cdum(1),
604 $ cdum(1), -1, ierr )
605 lwork_cungbr_p=cdum(1)
606 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
609 CALL cungbr(
'P', n, n, m, a, n, cdum(1),
610 $ cdum(1), -1, ierr )
611 lwork_cungbr_p=cdum(1)
612 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
614 IF( .NOT.wntun )
THEN
615 maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
620 maxwrk = max( minwrk, maxwrk )
623 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
629 CALL xerbla(
'CGESVD', -info )
631 ELSE IF( lquery )
THEN
637 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
644 smlnum = sqrt( slamch(
'S' ) ) / eps
645 bignum = one / smlnum
649 anrm = clange(
'M', m, n, a, lda, dum )
651 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
653 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
654 ELSE IF( anrm.GT.bignum )
THEN
656 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
665 IF( m.GE.mnthr )
THEN
679 CALL cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
680 $ lwork-iwork+1, ierr )
684 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
695 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
696 $ work( itaup ), work( iwork ), lwork-iwork+1,
699 IF( wntvo .OR. wntvas )
THEN
705 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
706 $ work( iwork ), lwork-iwork+1, ierr )
716 CALL cbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
717 $ cdum, 1, cdum, 1, rwork( irwork ), info )
722 $
CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
724 ELSE IF( wntuo .AND. wntvn )
THEN
730 IF( lwork.GE.n*n+3*n )
THEN
735 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
741 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
751 ldwrku = ( lwork-n*n ) / n
761 CALL cgeqrf( m, n, a, lda, work( itau ),
762 $ work( iwork ), lwork-iwork+1, ierr )
766 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
767 CALL claset(
'L', n-1, n-1, czero, czero,
768 $ work( ir+1 ), ldwrkr )
774 CALL cungqr( m, n, n, a, lda, work( itau ),
775 $ work( iwork ), lwork-iwork+1, ierr )
785 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
786 $ work( itauq ), work( itaup ),
787 $ work( iwork ), lwork-iwork+1, ierr )
793 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
794 $ work( itauq ), work( iwork ),
795 $ lwork-iwork+1, ierr )
803 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
804 $ work( ir ), ldwrkr, cdum, 1,
805 $ rwork( irwork ), info )
813 DO 10 i = 1, m, ldwrku
814 chunk = min( m-i+1, ldwrku )
815 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
816 $ lda, work( ir ), ldwrkr, czero,
817 $ work( iu ), ldwrku )
818 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
835 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
836 $ work( itauq ), work( itaup ),
837 $ work( iwork ), lwork-iwork+1, ierr )
843 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
844 $ work( iwork ), lwork-iwork+1, ierr )
852 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
853 $ a, lda, cdum, 1, rwork( irwork ), info )
857 ELSE IF( wntuo .AND. wntvas )
THEN
863 IF( lwork.GE.n*n+3*n )
THEN
868 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
874 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
884 ldwrku = ( lwork-n*n ) / n
894 CALL cgeqrf( m, n, a, lda, work( itau ),
895 $ work( iwork ), lwork-iwork+1, ierr )
899 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
901 $
CALL claset(
'L', n-1, n-1, czero, czero,
908 CALL cungqr( m, n, n, a, lda, work( itau ),
909 $ work( iwork ), lwork-iwork+1, ierr )
919 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
920 $ work( itauq ), work( itaup ),
921 $ work( iwork ), lwork-iwork+1, ierr )
922 CALL clacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
928 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
929 $ work( itauq ), work( iwork ),
930 $ lwork-iwork+1, ierr )
936 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
937 $ work( iwork ), lwork-iwork+1, ierr )
946 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
947 $ ldvt, work( ir ), ldwrkr, cdum, 1,
948 $ rwork( irwork ), info )
956 DO 20 i = 1, m, ldwrku
957 chunk = min( m-i+1, ldwrku )
958 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
959 $ lda, work( ir ), ldwrkr, czero,
960 $ work( iu ), ldwrku )
961 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
976 CALL cgeqrf( m, n, a, lda, work( itau ),
977 $ work( iwork ), lwork-iwork+1, ierr )
981 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
983 $
CALL claset(
'L', n-1, n-1, czero, czero,
990 CALL cungqr( m, n, n, a, lda, work( itau ),
991 $ work( iwork ), lwork-iwork+1, ierr )
1001 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1002 $ work( itauq ), work( itaup ),
1003 $ work( iwork ), lwork-iwork+1, ierr )
1009 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1010 $ work( itauq ), a, lda, work( iwork ),
1011 $ lwork-iwork+1, ierr )
1017 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1018 $ work( iwork ), lwork-iwork+1, ierr )
1027 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1028 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1033 ELSE IF( wntus )
THEN
1041 IF( lwork.GE.n*n+3*n )
THEN
1046 IF( lwork.GE.wrkbl+lda*n )
THEN
1057 itau = ir + ldwrkr*n
1064 CALL cgeqrf( m, n, a, lda, work( itau ),
1065 $ work( iwork ), lwork-iwork+1, ierr )
1069 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1071 CALL claset(
'L', n-1, n-1, czero, czero,
1072 $ work( ir+1 ), ldwrkr )
1078 CALL cungqr( m, n, n, a, lda, work( itau ),
1079 $ work( iwork ), lwork-iwork+1, ierr )
1089 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1090 $ rwork( ie ), work( itauq ),
1091 $ work( itaup ), work( iwork ),
1092 $ lwork-iwork+1, ierr )
1098 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1099 $ work( itauq ), work( iwork ),
1100 $ lwork-iwork+1, ierr )
1108 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1109 $ 1, work( ir ), ldwrkr, cdum, 1,
1110 $ rwork( irwork ), info )
1117 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1118 $ work( ir ), ldwrkr, czero, u, ldu )
1131 CALL cgeqrf( m, n, a, lda, work( itau ),
1132 $ work( iwork ), lwork-iwork+1, ierr )
1133 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1139 CALL cungqr( m, n, n, u, ldu, work( itau ),
1140 $ work( iwork ), lwork-iwork+1, ierr )
1148 CALL claset(
'L', n-1, n-1, czero, czero,
1155 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1156 $ work( itauq ), work( itaup ),
1157 $ work( iwork ), lwork-iwork+1, ierr )
1163 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1164 $ work( itauq ), u, ldu, work( iwork ),
1165 $ lwork-iwork+1, ierr )
1173 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1174 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1179 ELSE IF( wntvo )
THEN
1185 IF( lwork.GE.2*n*n+3*n )
THEN
1190 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1197 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1212 itau = ir + ldwrkr*n
1219 CALL cgeqrf( m, n, a, lda, work( itau ),
1220 $ work( iwork ), lwork-iwork+1, ierr )
1224 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1226 CALL claset(
'L', n-1, n-1, czero, czero,
1227 $ work( iu+1 ), ldwrku )
1233 CALL cungqr( m, n, n, a, lda, work( itau ),
1234 $ work( iwork ), lwork-iwork+1, ierr )
1246 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1247 $ rwork( ie ), work( itauq ),
1248 $ work( itaup ), work( iwork ),
1249 $ lwork-iwork+1, ierr )
1250 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1251 $ work( ir ), ldwrkr )
1257 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1258 $ work( itauq ), work( iwork ),
1259 $ lwork-iwork+1, ierr )
1266 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1267 $ work( itaup ), work( iwork ),
1268 $ lwork-iwork+1, ierr )
1277 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1278 $ work( ir ), ldwrkr, work( iu ),
1279 $ ldwrku, cdum, 1, rwork( irwork ),
1287 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1288 $ work( iu ), ldwrku, czero, u, ldu )
1294 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1308 CALL cgeqrf( m, n, a, lda, work( itau ),
1309 $ work( iwork ), lwork-iwork+1, ierr )
1310 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1316 CALL cungqr( m, n, n, u, ldu, work( itau ),
1317 $ work( iwork ), lwork-iwork+1, ierr )
1325 CALL claset(
'L', n-1, n-1, czero, czero,
1332 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1333 $ work( itauq ), work( itaup ),
1334 $ work( iwork ), lwork-iwork+1, ierr )
1340 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1341 $ work( itauq ), u, ldu, work( iwork ),
1342 $ lwork-iwork+1, ierr )
1348 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
1349 $ work( iwork ), lwork-iwork+1, ierr )
1358 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1359 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1364 ELSE IF( wntvas )
THEN
1371 IF( lwork.GE.n*n+3*n )
THEN
1376 IF( lwork.GE.wrkbl+lda*n )
THEN
1387 itau = iu + ldwrku*n
1394 CALL cgeqrf( m, n, a, lda, work( itau ),
1395 $ work( iwork ), lwork-iwork+1, ierr )
1399 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1401 CALL claset(
'L', n-1, n-1, czero, czero,
1402 $ work( iu+1 ), ldwrku )
1408 CALL cungqr( m, n, n, a, lda, work( itau ),
1409 $ work( iwork ), lwork-iwork+1, ierr )
1419 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1420 $ rwork( ie ), work( itauq ),
1421 $ work( itaup ), work( iwork ),
1422 $ lwork-iwork+1, ierr )
1423 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1430 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1431 $ work( itauq ), work( iwork ),
1432 $ lwork-iwork+1, ierr )
1439 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1440 $ work( iwork ), lwork-iwork+1, ierr )
1449 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1450 $ ldvt, work( iu ), ldwrku, cdum, 1,
1451 $ rwork( irwork ), info )
1458 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1459 $ work( iu ), ldwrku, czero, u, ldu )
1472 CALL cgeqrf( m, n, a, lda, work( itau ),
1473 $ work( iwork ), lwork-iwork+1, ierr )
1474 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1480 CALL cungqr( m, n, n, u, ldu, work( itau ),
1481 $ work( iwork ), lwork-iwork+1, ierr )
1485 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1487 $
CALL claset(
'L', n-1, n-1, czero, czero,
1488 $ vt( 2, 1 ), ldvt )
1498 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1499 $ work( itauq ), work( itaup ),
1500 $ work( iwork ), lwork-iwork+1, ierr )
1507 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1508 $ work( itauq ), u, ldu, work( iwork ),
1509 $ lwork-iwork+1, ierr )
1515 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1516 $ work( iwork ), lwork-iwork+1, ierr )
1525 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1526 $ ldvt, u, ldu, cdum, 1,
1527 $ rwork( irwork ), info )
1533 ELSE IF( wntua )
THEN
1541 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1546 IF( lwork.GE.wrkbl+lda*n )
THEN
1557 itau = ir + ldwrkr*n
1564 CALL cgeqrf( m, n, a, lda, work( itau ),
1565 $ work( iwork ), lwork-iwork+1, ierr )
1566 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1570 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1572 CALL claset(
'L', n-1, n-1, czero, czero,
1573 $ work( ir+1 ), ldwrkr )
1579 CALL cungqr( m, m, n, u, ldu, work( itau ),
1580 $ work( iwork ), lwork-iwork+1, ierr )
1590 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1591 $ rwork( ie ), work( itauq ),
1592 $ work( itaup ), work( iwork ),
1593 $ lwork-iwork+1, ierr )
1599 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1600 $ work( itauq ), work( iwork ),
1601 $ lwork-iwork+1, ierr )
1609 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1610 $ 1, work( ir ), ldwrkr, cdum, 1,
1611 $ rwork( irwork ), info )
1618 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1619 $ work( ir ), ldwrkr, czero, a, lda )
1623 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1636 CALL cgeqrf( m, n, a, lda, work( itau ),
1637 $ work( iwork ), lwork-iwork+1, ierr )
1638 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1644 CALL cungqr( m, m, n, u, ldu, work( itau ),
1645 $ work( iwork ), lwork-iwork+1, ierr )
1653 CALL claset(
'L', n-1, n-1, czero, czero,
1660 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1661 $ work( itauq ), work( itaup ),
1662 $ work( iwork ), lwork-iwork+1, ierr )
1669 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1670 $ work( itauq ), u, ldu, work( iwork ),
1671 $ lwork-iwork+1, ierr )
1679 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1680 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1685 ELSE IF( wntvo )
THEN
1691 IF( lwork.GE.2*n*n+max( n+m, 3*n ) )
THEN
1696 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1703 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1718 itau = ir + ldwrkr*n
1725 CALL cgeqrf( m, n, a, lda, work( itau ),
1726 $ work( iwork ), lwork-iwork+1, ierr )
1727 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1733 CALL cungqr( m, m, n, u, ldu, work( itau ),
1734 $ work( iwork ), lwork-iwork+1, ierr )
1738 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1740 CALL claset(
'L', n-1, n-1, czero, czero,
1741 $ work( iu+1 ), ldwrku )
1753 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1754 $ rwork( ie ), work( itauq ),
1755 $ work( itaup ), work( iwork ),
1756 $ lwork-iwork+1, ierr )
1757 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1758 $ work( ir ), ldwrkr )
1764 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1765 $ work( itauq ), work( iwork ),
1766 $ lwork-iwork+1, ierr )
1773 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1774 $ work( itaup ), work( iwork ),
1775 $ lwork-iwork+1, ierr )
1784 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1785 $ work( ir ), ldwrkr, work( iu ),
1786 $ ldwrku, cdum, 1, rwork( irwork ),
1794 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1795 $ work( iu ), ldwrku, czero, a, lda )
1799 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1803 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1817 CALL cgeqrf( m, n, a, lda, work( itau ),
1818 $ work( iwork ), lwork-iwork+1, ierr )
1819 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1825 CALL cungqr( m, m, n, u, ldu, work( itau ),
1826 $ work( iwork ), lwork-iwork+1, ierr )
1834 CALL claset(
'L', n-1, n-1, czero, czero,
1841 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1842 $ work( itauq ), work( itaup ),
1843 $ work( iwork ), lwork-iwork+1, ierr )
1850 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1851 $ work( itauq ), u, ldu, work( iwork ),
1852 $ lwork-iwork+1, ierr )
1858 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
1859 $ work( iwork ), lwork-iwork+1, ierr )
1868 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1869 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1874 ELSE IF( wntvas )
THEN
1881 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1886 IF( lwork.GE.wrkbl+lda*n )
THEN
1897 itau = iu + ldwrku*n
1904 CALL cgeqrf( m, n, a, lda, work( itau ),
1905 $ work( iwork ), lwork-iwork+1, ierr )
1906 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1912 CALL cungqr( m, m, n, u, ldu, work( itau ),
1913 $ work( iwork ), lwork-iwork+1, ierr )
1917 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1919 CALL claset(
'L', n-1, n-1, czero, czero,
1920 $ work( iu+1 ), ldwrku )
1930 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1931 $ rwork( ie ), work( itauq ),
1932 $ work( itaup ), work( iwork ),
1933 $ lwork-iwork+1, ierr )
1934 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1941 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1942 $ work( itauq ), work( iwork ),
1943 $ lwork-iwork+1, ierr )
1950 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1951 $ work( iwork ), lwork-iwork+1, ierr )
1960 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1961 $ ldvt, work( iu ), ldwrku, cdum, 1,
1962 $ rwork( irwork ), info )
1969 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1970 $ work( iu ), ldwrku, czero, a, lda )
1974 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1987 CALL cgeqrf( m, n, a, lda, work( itau ),
1988 $ work( iwork ), lwork-iwork+1, ierr )
1989 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1995 CALL cungqr( m, m, n, u, ldu, work( itau ),
1996 $ work( iwork ), lwork-iwork+1, ierr )
2000 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2002 $
CALL claset(
'L', n-1, n-1, czero, czero,
2003 $ vt( 2, 1 ), ldvt )
2013 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2014 $ work( itauq ), work( itaup ),
2015 $ work( iwork ), lwork-iwork+1, ierr )
2022 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2023 $ work( itauq ), u, ldu, work( iwork ),
2024 $ lwork-iwork+1, ierr )
2030 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2031 $ work( iwork ), lwork-iwork+1, ierr )
2040 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
2041 $ ldvt, u, ldu, cdum, 1,
2042 $ rwork( irwork ), info )
2066 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2067 $ work( itaup ), work( iwork ), lwork-iwork+1,
2076 CALL clacpy(
'L', m, n, a, lda, u, ldu )
2081 CALL cungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2082 $ work( iwork ), lwork-iwork+1, ierr )
2091 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2092 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2093 $ work( iwork ), lwork-iwork+1, ierr )
2102 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
2103 $ work( iwork ), lwork-iwork+1, ierr )
2112 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
2113 $ work( iwork ), lwork-iwork+1, ierr )
2116 IF( wntuas .OR. wntuo )
2120 IF( wntvas .OR. wntvo )
2124 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2132 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2133 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2135 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2143 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2144 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2154 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2155 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2167 IF( n.GE.mnthr )
THEN
2181 CALL cgelqf( m, n, a, lda, work( itau ), work( iwork ),
2182 $ lwork-iwork+1, ierr )
2186 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2197 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2198 $ work( itaup ), work( iwork ), lwork-iwork+1,
2200 IF( wntuo .OR. wntuas )
THEN
2206 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2207 $ work( iwork ), lwork-iwork+1, ierr )
2211 IF( wntuo .OR. wntuas )
2219 CALL cbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2220 $ a, lda, cdum, 1, rwork( irwork ), info )
2225 $
CALL clacpy(
'F', m, m, a, lda, u, ldu )
2227 ELSE IF( wntvo .AND. wntun )
THEN
2233 IF( lwork.GE.m*m+3*m )
THEN
2238 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2245 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2257 chunk = ( lwork-m*m ) / m
2260 itau = ir + ldwrkr*m
2267 CALL cgelqf( m, n, a, lda, work( itau ),
2268 $ work( iwork ), lwork-iwork+1, ierr )
2272 CALL clacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2273 CALL claset(
'U', m-1, m-1, czero, czero,
2274 $ work( ir+ldwrkr ), ldwrkr )
2280 CALL cunglq( m, n, m, a, lda, work( itau ),
2281 $ work( iwork ), lwork-iwork+1, ierr )
2291 CALL cgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2292 $ work( itauq ), work( itaup ),
2293 $ work( iwork ), lwork-iwork+1, ierr )
2299 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2300 $ work( itaup ), work( iwork ),
2301 $ lwork-iwork+1, ierr )
2309 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2310 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2311 $ rwork( irwork ), info )
2319 DO 30 i = 1, n, chunk
2320 blk = min( n-i+1, chunk )
2321 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2322 $ ldwrkr, a( 1, i ), lda, czero,
2323 $ work( iu ), ldwrku )
2324 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2341 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
2342 $ work( itauq ), work( itaup ),
2343 $ work( iwork ), lwork-iwork+1, ierr )
2349 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
2350 $ work( iwork ), lwork-iwork+1, ierr )
2358 CALL cbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2359 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2363 ELSE IF( wntvo .AND. wntuas )
THEN
2369 IF( lwork.GE.m*m+3*m )
THEN
2374 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2381 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2393 chunk = ( lwork-m*m ) / m
2396 itau = ir + ldwrkr*m
2403 CALL cgelqf( m, n, a, lda, work( itau ),
2404 $ work( iwork ), lwork-iwork+1, ierr )
2408 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2409 CALL claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2416 CALL cunglq( m, n, m, a, lda, work( itau ),
2417 $ work( iwork ), lwork-iwork+1, ierr )
2427 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2428 $ work( itauq ), work( itaup ),
2429 $ work( iwork ), lwork-iwork+1, ierr )
2430 CALL clacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2436 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2437 $ work( itaup ), work( iwork ),
2438 $ lwork-iwork+1, ierr )
2444 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2445 $ work( iwork ), lwork-iwork+1, ierr )
2454 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2455 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2456 $ rwork( irwork ), info )
2464 DO 40 i = 1, n, chunk
2465 blk = min( n-i+1, chunk )
2466 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2467 $ ldwrkr, a( 1, i ), lda, czero,
2468 $ work( iu ), ldwrku )
2469 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2484 CALL cgelqf( m, n, a, lda, work( itau ),
2485 $ work( iwork ), lwork-iwork+1, ierr )
2489 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2490 CALL claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2497 CALL cunglq( m, n, m, a, lda, work( itau ),
2498 $ work( iwork ), lwork-iwork+1, ierr )
2508 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2509 $ work( itauq ), work( itaup ),
2510 $ work( iwork ), lwork-iwork+1, ierr )
2516 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2517 $ work( itaup ), a, lda, work( iwork ),
2518 $ lwork-iwork+1, ierr )
2524 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2525 $ work( iwork ), lwork-iwork+1, ierr )
2534 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a, lda,
2535 $ u, ldu, cdum, 1, rwork( irwork ), info )
2539 ELSE IF( wntvs )
THEN
2547 IF( lwork.GE.m*m+3*m )
THEN
2552 IF( lwork.GE.wrkbl+lda*m )
THEN
2563 itau = ir + ldwrkr*m
2570 CALL cgelqf( m, n, a, lda, work( itau ),
2571 $ work( iwork ), lwork-iwork+1, ierr )
2575 CALL clacpy(
'L', m, m, a, lda, work( ir ),
2577 CALL claset(
'U', m-1, m-1, czero, czero,
2578 $ work( ir+ldwrkr ), ldwrkr )
2584 CALL cunglq( m, n, m, a, lda, work( itau ),
2585 $ work( iwork ), lwork-iwork+1, ierr )
2595 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2596 $ rwork( ie ), work( itauq ),
2597 $ work( itaup ), work( iwork ),
2598 $ lwork-iwork+1, ierr )
2605 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2606 $ work( itaup ), work( iwork ),
2607 $ lwork-iwork+1, ierr )
2615 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2616 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2617 $ rwork( irwork ), info )
2624 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
2625 $ ldwrkr, a, lda, czero, vt, ldvt )
2638 CALL cgelqf( m, n, a, lda, work( itau ),
2639 $ work( iwork ), lwork-iwork+1, ierr )
2643 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2649 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2650 $ work( iwork ), lwork-iwork+1, ierr )
2658 CALL claset(
'U', m-1, m-1, czero, czero,
2665 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2666 $ work( itauq ), work( itaup ),
2667 $ work( iwork ), lwork-iwork+1, ierr )
2673 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2674 $ work( itaup ), vt, ldvt,
2675 $ work( iwork ), lwork-iwork+1, ierr )
2683 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
2684 $ ldvt, cdum, 1, cdum, 1,
2685 $ rwork( irwork ), info )
2689 ELSE IF( wntuo )
THEN
2695 IF( lwork.GE.2*m*m+3*m )
THEN
2700 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2707 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2722 itau = ir + ldwrkr*m
2729 CALL cgelqf( m, n, a, lda, work( itau ),
2730 $ work( iwork ), lwork-iwork+1, ierr )
2734 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2736 CALL claset(
'U', m-1, m-1, czero, czero,
2737 $ work( iu+ldwrku ), ldwrku )
2743 CALL cunglq( m, n, m, a, lda, work( itau ),
2744 $ work( iwork ), lwork-iwork+1, ierr )
2756 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2757 $ rwork( ie ), work( itauq ),
2758 $ work( itaup ), work( iwork ),
2759 $ lwork-iwork+1, ierr )
2760 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
2761 $ work( ir ), ldwrkr )
2768 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2769 $ work( itaup ), work( iwork ),
2770 $ lwork-iwork+1, ierr )
2776 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2777 $ work( itauq ), work( iwork ),
2778 $ lwork-iwork+1, ierr )
2787 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2788 $ work( iu ), ldwrku, work( ir ),
2789 $ ldwrkr, cdum, 1, rwork( irwork ),
2797 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2798 $ ldwrku, a, lda, czero, vt, ldvt )
2804 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
2818 CALL cgelqf( m, n, a, lda, work( itau ),
2819 $ work( iwork ), lwork-iwork+1, ierr )
2820 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2826 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2827 $ work( iwork ), lwork-iwork+1, ierr )
2835 CALL claset(
'U', m-1, m-1, czero, czero,
2842 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2843 $ work( itauq ), work( itaup ),
2844 $ work( iwork ), lwork-iwork+1, ierr )
2850 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2851 $ work( itaup ), vt, ldvt,
2852 $ work( iwork ), lwork-iwork+1, ierr )
2858 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2859 $ work( iwork ), lwork-iwork+1, ierr )
2868 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
2869 $ ldvt, a, lda, cdum, 1,
2870 $ rwork( irwork ), info )
2874 ELSE IF( wntuas )
THEN
2881 IF( lwork.GE.m*m+3*m )
THEN
2886 IF( lwork.GE.wrkbl+lda*m )
THEN
2897 itau = iu + ldwrku*m
2904 CALL cgelqf( m, n, a, lda, work( itau ),
2905 $ work( iwork ), lwork-iwork+1, ierr )
2909 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2911 CALL claset(
'U', m-1, m-1, czero, czero,
2912 $ work( iu+ldwrku ), ldwrku )
2918 CALL cunglq( m, n, m, a, lda, work( itau ),
2919 $ work( iwork ), lwork-iwork+1, ierr )
2929 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2930 $ rwork( ie ), work( itauq ),
2931 $ work( itaup ), work( iwork ),
2932 $ lwork-iwork+1, ierr )
2933 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
2941 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2942 $ work( itaup ), work( iwork ),
2943 $ lwork-iwork+1, ierr )
2949 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2950 $ work( iwork ), lwork-iwork+1, ierr )
2959 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2960 $ work( iu ), ldwrku, u, ldu, cdum, 1,
2961 $ rwork( irwork ), info )
2968 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2969 $ ldwrku, a, lda, czero, vt, ldvt )
2982 CALL cgelqf( m, n, a, lda, work( itau ),
2983 $ work( iwork ), lwork-iwork+1, ierr )
2984 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2990 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2991 $ work( iwork ), lwork-iwork+1, ierr )
2995 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2996 CALL claset(
'U', m-1, m-1, czero, czero,
3007 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3008 $ work( itauq ), work( itaup ),
3009 $ work( iwork ), lwork-iwork+1, ierr )
3016 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3017 $ work( itaup ), vt, ldvt,
3018 $ work( iwork ), lwork-iwork+1, ierr )
3024 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3025 $ work( iwork ), lwork-iwork+1, ierr )
3034 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3035 $ ldvt, u, ldu, cdum, 1,
3036 $ rwork( irwork ), info )
3042 ELSE IF( wntva )
THEN
3050 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3055 IF( lwork.GE.wrkbl+lda*m )
THEN
3066 itau = ir + ldwrkr*m
3073 CALL cgelqf( m, n, a, lda, work( itau ),
3074 $ work( iwork ), lwork-iwork+1, ierr )
3075 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3079 CALL clacpy(
'L', m, m, a, lda, work( ir ),
3081 CALL claset(
'U', m-1, m-1, czero, czero,
3082 $ work( ir+ldwrkr ), ldwrkr )
3088 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3089 $ work( iwork ), lwork-iwork+1, ierr )
3099 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
3100 $ rwork( ie ), work( itauq ),
3101 $ work( itaup ), work( iwork ),
3102 $ lwork-iwork+1, ierr )
3109 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
3110 $ work( itaup ), work( iwork ),
3111 $ lwork-iwork+1, ierr )
3119 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3120 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3121 $ rwork( irwork ), info )
3128 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
3129 $ ldwrkr, vt, ldvt, czero, a, lda )
3133 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3146 CALL cgelqf( m, n, a, lda, work( itau ),
3147 $ work( iwork ), lwork-iwork+1, ierr )
3148 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3154 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3155 $ work( iwork ), lwork-iwork+1, ierr )
3163 CALL claset(
'U', m-1, m-1, czero, czero,
3170 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3171 $ work( itauq ), work( itaup ),
3172 $ work( iwork ), lwork-iwork+1, ierr )
3179 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3180 $ work( itaup ), vt, ldvt,
3181 $ work( iwork ), lwork-iwork+1, ierr )
3189 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
3190 $ ldvt, cdum, 1, cdum, 1,
3191 $ rwork( irwork ), info )
3195 ELSE IF( wntuo )
THEN
3201 IF( lwork.GE.2*m*m+max( n+m, 3*m ) )
THEN
3206 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3213 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3228 itau = ir + ldwrkr*m
3235 CALL cgelqf( m, n, a, lda, work( itau ),
3236 $ work( iwork ), lwork-iwork+1, ierr )
3237 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3243 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3244 $ work( iwork ), lwork-iwork+1, ierr )
3248 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3250 CALL claset(
'U', m-1, m-1, czero, czero,
3251 $ work( iu+ldwrku ), ldwrku )
3263 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3264 $ rwork( ie ), work( itauq ),
3265 $ work( itaup ), work( iwork ),
3266 $ lwork-iwork+1, ierr )
3267 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
3268 $ work( ir ), ldwrkr )
3275 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3276 $ work( itaup ), work( iwork ),
3277 $ lwork-iwork+1, ierr )
3283 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3284 $ work( itauq ), work( iwork ),
3285 $ lwork-iwork+1, ierr )
3294 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3295 $ work( iu ), ldwrku, work( ir ),
3296 $ ldwrkr, cdum, 1, rwork( irwork ),
3304 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3305 $ ldwrku, vt, ldvt, czero, a, lda )
3309 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3313 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
3327 CALL cgelqf( m, n, a, lda, work( itau ),
3328 $ work( iwork ), lwork-iwork+1, ierr )
3329 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3335 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3336 $ work( iwork ), lwork-iwork+1, ierr )
3344 CALL claset(
'U', m-1, m-1, czero, czero,
3351 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3352 $ work( itauq ), work( itaup ),
3353 $ work( iwork ), lwork-iwork+1, ierr )
3360 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3361 $ work( itaup ), vt, ldvt,
3362 $ work( iwork ), lwork-iwork+1, ierr )
3368 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
3369 $ work( iwork ), lwork-iwork+1, ierr )
3378 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3379 $ ldvt, a, lda, cdum, 1,
3380 $ rwork( irwork ), info )
3384 ELSE IF( wntuas )
THEN
3391 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3396 IF( lwork.GE.wrkbl+lda*m )
THEN
3407 itau = iu + ldwrku*m
3414 CALL cgelqf( m, n, a, lda, work( itau ),
3415 $ work( iwork ), lwork-iwork+1, ierr )
3416 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3422 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3423 $ work( iwork ), lwork-iwork+1, ierr )
3427 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3429 CALL claset(
'U', m-1, m-1, czero, czero,
3430 $ work( iu+ldwrku ), ldwrku )
3440 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3441 $ rwork( ie ), work( itauq ),
3442 $ work( itaup ), work( iwork ),
3443 $ lwork-iwork+1, ierr )
3444 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
3451 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3452 $ work( itaup ), work( iwork ),
3453 $ lwork-iwork+1, ierr )
3459 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3460 $ work( iwork ), lwork-iwork+1, ierr )
3469 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3470 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3471 $ rwork( irwork ), info )
3478 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3479 $ ldwrku, vt, ldvt, czero, a, lda )
3483 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3496 CALL cgelqf( m, n, a, lda, work( itau ),
3497 $ work( iwork ), lwork-iwork+1, ierr )
3498 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3504 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3505 $ work( iwork ), lwork-iwork+1, ierr )
3509 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3510 CALL claset(
'U', m-1, m-1, czero, czero,
3521 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3522 $ work( itauq ), work( itaup ),
3523 $ work( iwork ), lwork-iwork+1, ierr )
3530 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3531 $ work( itaup ), vt, ldvt,
3532 $ work( iwork ), lwork-iwork+1, ierr )
3538 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3539 $ work( iwork ), lwork-iwork+1, ierr )
3548 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3549 $ ldvt, u, ldu, cdum, 1,
3550 $ rwork( irwork ), info )
3574 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3575 $ work( itaup ), work( iwork ), lwork-iwork+1,
3584 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3585 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3586 $ work( iwork ), lwork-iwork+1, ierr )
3595 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3600 CALL cungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3601 $ work( iwork ), lwork-iwork+1, ierr )
3610 CALL cungbr(
'Q', m, m, n, a, lda, work( itauq ),
3611 $ work( iwork ), lwork-iwork+1, ierr )
3620 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
3621 $ work( iwork ), lwork-iwork+1, ierr )
3624 IF( wntuas .OR. wntuo )
3628 IF( wntvas .OR. wntvo )
3632 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3640 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3641 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3643 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3651 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3652 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3662 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3663 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3673 IF( iscl.EQ.1 )
THEN
3674 IF( anrm.GT.bignum )
3675 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3677 IF( info.NE.0 .AND. anrm.GT.bignum )
3678 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3679 $ rwork( ie ), minmn, ierr )
3680 IF( anrm.LT.smlnum )
3681 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3683 IF( info.NE.0 .AND. anrm.LT.smlnum )
3684 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3685 $ 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...
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
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.
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