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, 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...
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 ...
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
logical function lsame(CA, CB)
LSAME
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
real function slamch(CMACH)
SLAMCH
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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