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 )
325 lwork_cgeqrf = int( cdum(1) )
327 CALL cungqr( m, n, n, a, lda, cdum(1), cdum(1), -1, ierr )
328 lwork_cungqr_n = int( cdum(1) )
329 CALL cungqr( m, m, n, a, lda, cdum(1), cdum(1), -1, ierr )
330 lwork_cungqr_m = int( cdum(1) )
332 CALL cgebrd( n, n, a, lda, s, dum(1), cdum(1),
333 $ cdum(1), cdum(1), -1, ierr )
334 lwork_cgebrd = int( cdum(1) )
336 CALL cungbr(
'P', n, n, n, a, lda, cdum(1),
337 $ cdum(1), -1, ierr )
338 lwork_cungbr_p = int( cdum(1) )
339 CALL cungbr(
'Q', n, n, n, a, lda, cdum(1),
340 $ cdum(1), -1, ierr )
341 lwork_cungbr_q = int( 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 )
449 lwork_cgebrd = int( cdum(1) )
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 = int( 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 = int( 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 )
475 lwork_cgelqf = int( cdum(1) )
477 CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
479 lwork_cunglq_n = int( cdum(1) )
480 CALL cunglq( m, n, m, a, lda, cdum(1), cdum(1), -1, ierr )
481 lwork_cunglq_m = int( cdum(1) )
483 CALL cgebrd( m, m, a, lda, s, dum(1), cdum(1),
484 $ cdum(1), cdum(1), -1, ierr )
485 lwork_cgebrd = int( cdum(1) )
487 CALL cungbr(
'P', m, m, m, a, n, cdum(1),
488 $ cdum(1), -1, ierr )
489 lwork_cungbr_p = int( cdum(1) )
491 CALL cungbr(
'Q', m, m, m, a, n, cdum(1),
492 $ cdum(1), -1, ierr )
493 lwork_cungbr_q = int( 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 )
599 lwork_cgebrd = int( cdum(1) )
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 = int( 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 = int( 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 )
685 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
697 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
698 $ work( itaup ), work( iwork ), lwork-iwork+1,
701 IF( wntvo .OR. wntvas )
THEN
707 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
708 $ work( iwork ), lwork-iwork+1, ierr )
718 CALL cbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
719 $ cdum, 1, cdum, 1, rwork( irwork ), info )
724 $
CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
726 ELSE IF( wntuo .AND. wntvn )
THEN
732 IF( lwork.GE.n*n+3*n )
THEN
737 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
743 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
753 ldwrku = ( lwork-n*n ) / n
763 CALL cgeqrf( m, n, a, lda, work( itau ),
764 $ work( iwork ), lwork-iwork+1, ierr )
768 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
769 CALL claset(
'L', n-1, n-1, czero, czero,
770 $ work( ir+1 ), ldwrkr )
776 CALL cungqr( m, n, n, a, lda, work( itau ),
777 $ work( iwork ), lwork-iwork+1, ierr )
787 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
788 $ work( itauq ), work( itaup ),
789 $ work( iwork ), lwork-iwork+1, ierr )
795 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
796 $ work( itauq ), work( iwork ),
797 $ lwork-iwork+1, ierr )
805 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
806 $ work( ir ), ldwrkr, cdum, 1,
807 $ rwork( irwork ), info )
815 DO 10 i = 1, m, ldwrku
816 chunk = min( m-i+1, ldwrku )
817 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
818 $ lda, work( ir ), ldwrkr, czero,
819 $ work( iu ), ldwrku )
820 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
837 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
838 $ work( itauq ), work( itaup ),
839 $ work( iwork ), lwork-iwork+1, ierr )
845 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
846 $ work( iwork ), lwork-iwork+1, ierr )
854 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
855 $ a, lda, cdum, 1, rwork( irwork ), info )
859 ELSE IF( wntuo .AND. wntvas )
THEN
865 IF( lwork.GE.n*n+3*n )
THEN
870 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
876 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
886 ldwrku = ( lwork-n*n ) / n
896 CALL cgeqrf( m, n, a, lda, work( itau ),
897 $ work( iwork ), lwork-iwork+1, ierr )
901 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
903 $
CALL claset(
'L', n-1, n-1, czero, czero,
910 CALL cungqr( m, n, n, a, lda, work( itau ),
911 $ work( iwork ), lwork-iwork+1, ierr )
921 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
922 $ work( itauq ), work( itaup ),
923 $ work( iwork ), lwork-iwork+1, ierr )
924 CALL clacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
930 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
931 $ work( itauq ), work( iwork ),
932 $ lwork-iwork+1, ierr )
938 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
939 $ work( iwork ), lwork-iwork+1, ierr )
948 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
949 $ ldvt, work( ir ), ldwrkr, cdum, 1,
950 $ rwork( irwork ), info )
958 DO 20 i = 1, m, ldwrku
959 chunk = min( m-i+1, ldwrku )
960 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
961 $ lda, work( ir ), ldwrkr, czero,
962 $ work( iu ), ldwrku )
963 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
978 CALL cgeqrf( m, n, a, lda, work( itau ),
979 $ work( iwork ), lwork-iwork+1, ierr )
983 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
985 $
CALL claset(
'L', n-1, n-1, czero, czero,
992 CALL cungqr( m, n, n, a, lda, work( itau ),
993 $ work( iwork ), lwork-iwork+1, ierr )
1003 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1004 $ work( itauq ), work( itaup ),
1005 $ work( iwork ), lwork-iwork+1, ierr )
1011 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1012 $ work( itauq ), a, lda, work( iwork ),
1013 $ lwork-iwork+1, ierr )
1019 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1020 $ work( iwork ), lwork-iwork+1, ierr )
1029 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1030 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1035 ELSE IF( wntus )
THEN
1043 IF( lwork.GE.n*n+3*n )
THEN
1048 IF( lwork.GE.wrkbl+lda*n )
THEN
1059 itau = ir + ldwrkr*n
1066 CALL cgeqrf( m, n, a, lda, work( itau ),
1067 $ work( iwork ), lwork-iwork+1, ierr )
1071 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1073 CALL claset(
'L', n-1, n-1, czero, czero,
1074 $ work( ir+1 ), ldwrkr )
1080 CALL cungqr( m, n, n, a, lda, work( itau ),
1081 $ work( iwork ), lwork-iwork+1, ierr )
1091 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1092 $ rwork( ie ), work( itauq ),
1093 $ work( itaup ), work( iwork ),
1094 $ lwork-iwork+1, ierr )
1100 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1101 $ work( itauq ), work( iwork ),
1102 $ lwork-iwork+1, ierr )
1110 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1111 $ 1, work( ir ), ldwrkr, cdum, 1,
1112 $ rwork( irwork ), info )
1119 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1120 $ work( ir ), ldwrkr, czero, u, ldu )
1133 CALL cgeqrf( m, n, a, lda, work( itau ),
1134 $ work( iwork ), lwork-iwork+1, ierr )
1135 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1141 CALL cungqr( m, n, n, u, ldu, work( itau ),
1142 $ work( iwork ), lwork-iwork+1, ierr )
1151 CALL claset(
'L', n-1, n-1, czero, czero,
1159 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1160 $ work( itauq ), work( itaup ),
1161 $ work( iwork ), lwork-iwork+1, ierr )
1167 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1168 $ work( itauq ), u, ldu, work( iwork ),
1169 $ lwork-iwork+1, ierr )
1177 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1178 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1183 ELSE IF( wntvo )
THEN
1189 IF( lwork.GE.2*n*n+3*n )
THEN
1194 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1201 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1216 itau = ir + ldwrkr*n
1223 CALL cgeqrf( m, n, a, lda, work( itau ),
1224 $ work( iwork ), lwork-iwork+1, ierr )
1228 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1230 CALL claset(
'L', n-1, n-1, czero, czero,
1231 $ work( iu+1 ), ldwrku )
1237 CALL cungqr( m, n, n, a, lda, work( itau ),
1238 $ work( iwork ), lwork-iwork+1, ierr )
1250 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1251 $ rwork( ie ), work( itauq ),
1252 $ work( itaup ), work( iwork ),
1253 $ lwork-iwork+1, ierr )
1254 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1255 $ work( ir ), ldwrkr )
1261 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1262 $ work( itauq ), work( iwork ),
1263 $ lwork-iwork+1, ierr )
1270 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1271 $ work( itaup ), work( iwork ),
1272 $ lwork-iwork+1, ierr )
1281 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1282 $ work( ir ), ldwrkr, work( iu ),
1283 $ ldwrku, cdum, 1, rwork( irwork ),
1291 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1292 $ work( iu ), ldwrku, czero, u, ldu )
1298 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1312 CALL cgeqrf( m, n, a, lda, work( itau ),
1313 $ work( iwork ), lwork-iwork+1, ierr )
1314 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1320 CALL cungqr( m, n, n, u, ldu, work( itau ),
1321 $ work( iwork ), lwork-iwork+1, ierr )
1330 CALL claset(
'L', n-1, n-1, czero, czero,
1338 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1339 $ work( itauq ), work( itaup ),
1340 $ work( iwork ), lwork-iwork+1, ierr )
1346 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1347 $ work( itauq ), u, ldu, work( iwork ),
1348 $ lwork-iwork+1, ierr )
1354 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
1355 $ work( iwork ), lwork-iwork+1, ierr )
1364 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1365 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1370 ELSE IF( wntvas )
THEN
1377 IF( lwork.GE.n*n+3*n )
THEN
1382 IF( lwork.GE.wrkbl+lda*n )
THEN
1393 itau = iu + ldwrku*n
1400 CALL cgeqrf( m, n, a, lda, work( itau ),
1401 $ work( iwork ), lwork-iwork+1, ierr )
1405 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1407 CALL claset(
'L', n-1, n-1, czero, czero,
1408 $ work( iu+1 ), ldwrku )
1414 CALL cungqr( m, n, n, a, lda, work( itau ),
1415 $ work( iwork ), lwork-iwork+1, ierr )
1425 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1426 $ rwork( ie ), work( itauq ),
1427 $ work( itaup ), work( iwork ),
1428 $ lwork-iwork+1, ierr )
1429 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1436 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1437 $ work( itauq ), work( iwork ),
1438 $ lwork-iwork+1, ierr )
1445 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1446 $ work( iwork ), lwork-iwork+1, ierr )
1455 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1456 $ ldvt, work( iu ), ldwrku, cdum, 1,
1457 $ rwork( irwork ), info )
1464 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1465 $ work( iu ), ldwrku, czero, u, ldu )
1478 CALL cgeqrf( m, n, a, lda, work( itau ),
1479 $ work( iwork ), lwork-iwork+1, ierr )
1480 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1486 CALL cungqr( m, n, n, u, ldu, work( itau ),
1487 $ work( iwork ), lwork-iwork+1, ierr )
1491 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1493 $
CALL claset(
'L', n-1, n-1, czero, czero,
1494 $ vt( 2, 1 ), ldvt )
1504 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1505 $ work( itauq ), work( itaup ),
1506 $ work( iwork ), lwork-iwork+1, ierr )
1513 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1514 $ work( itauq ), u, ldu, work( iwork ),
1515 $ lwork-iwork+1, ierr )
1521 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1522 $ work( iwork ), lwork-iwork+1, ierr )
1531 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1532 $ ldvt, u, ldu, cdum, 1,
1533 $ rwork( irwork ), info )
1539 ELSE IF( wntua )
THEN
1547 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1552 IF( lwork.GE.wrkbl+lda*n )
THEN
1563 itau = ir + ldwrkr*n
1570 CALL cgeqrf( m, n, a, lda, work( itau ),
1571 $ work( iwork ), lwork-iwork+1, ierr )
1572 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1576 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1578 CALL claset(
'L', n-1, n-1, czero, czero,
1579 $ work( ir+1 ), ldwrkr )
1585 CALL cungqr( m, m, n, u, ldu, work( itau ),
1586 $ work( iwork ), lwork-iwork+1, ierr )
1596 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1597 $ rwork( ie ), work( itauq ),
1598 $ work( itaup ), work( iwork ),
1599 $ lwork-iwork+1, ierr )
1605 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1606 $ work( itauq ), work( iwork ),
1607 $ lwork-iwork+1, ierr )
1615 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1616 $ 1, work( ir ), ldwrkr, cdum, 1,
1617 $ rwork( irwork ), info )
1624 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1625 $ work( ir ), ldwrkr, czero, a, lda )
1629 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1642 CALL cgeqrf( m, n, a, lda, work( itau ),
1643 $ work( iwork ), lwork-iwork+1, ierr )
1644 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1650 CALL cungqr( m, m, n, u, ldu, work( itau ),
1651 $ work( iwork ), lwork-iwork+1, ierr )
1660 CALL claset(
'L', n-1, n-1, czero, czero,
1668 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1669 $ work( itauq ), work( itaup ),
1670 $ work( iwork ), lwork-iwork+1, ierr )
1677 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1678 $ work( itauq ), u, ldu, work( iwork ),
1679 $ lwork-iwork+1, ierr )
1687 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1688 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1693 ELSE IF( wntvo )
THEN
1699 IF( lwork.GE.2*n*n+max( n+m, 3*n ) )
THEN
1704 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1711 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1726 itau = ir + ldwrkr*n
1733 CALL cgeqrf( m, n, a, lda, work( itau ),
1734 $ work( iwork ), lwork-iwork+1, ierr )
1735 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1741 CALL cungqr( m, m, n, u, ldu, work( itau ),
1742 $ work( iwork ), lwork-iwork+1, ierr )
1746 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1748 CALL claset(
'L', n-1, n-1, czero, czero,
1749 $ work( iu+1 ), ldwrku )
1761 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1762 $ rwork( ie ), work( itauq ),
1763 $ work( itaup ), work( iwork ),
1764 $ lwork-iwork+1, ierr )
1765 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1766 $ work( ir ), ldwrkr )
1772 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1773 $ work( itauq ), work( iwork ),
1774 $ lwork-iwork+1, ierr )
1781 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1782 $ work( itaup ), work( iwork ),
1783 $ lwork-iwork+1, ierr )
1792 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1793 $ work( ir ), ldwrkr, work( iu ),
1794 $ ldwrku, cdum, 1, rwork( irwork ),
1802 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1803 $ work( iu ), ldwrku, czero, a, lda )
1807 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1811 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1825 CALL cgeqrf( m, n, a, lda, work( itau ),
1826 $ work( iwork ), lwork-iwork+1, ierr )
1827 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1833 CALL cungqr( m, m, n, u, ldu, work( itau ),
1834 $ work( iwork ), lwork-iwork+1, ierr )
1843 CALL claset(
'L', n-1, n-1, czero, czero,
1851 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1852 $ work( itauq ), work( itaup ),
1853 $ work( iwork ), lwork-iwork+1, ierr )
1860 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1861 $ work( itauq ), u, ldu, work( iwork ),
1862 $ lwork-iwork+1, ierr )
1868 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
1869 $ work( iwork ), lwork-iwork+1, ierr )
1878 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1879 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1884 ELSE IF( wntvas )
THEN
1891 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1896 IF( lwork.GE.wrkbl+lda*n )
THEN
1907 itau = iu + ldwrku*n
1914 CALL cgeqrf( m, n, a, lda, work( itau ),
1915 $ work( iwork ), lwork-iwork+1, ierr )
1916 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1922 CALL cungqr( m, m, n, u, ldu, work( itau ),
1923 $ work( iwork ), lwork-iwork+1, ierr )
1927 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1929 CALL claset(
'L', n-1, n-1, czero, czero,
1930 $ work( iu+1 ), ldwrku )
1940 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1941 $ rwork( ie ), work( itauq ),
1942 $ work( itaup ), work( iwork ),
1943 $ lwork-iwork+1, ierr )
1944 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1951 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1952 $ work( itauq ), work( iwork ),
1953 $ lwork-iwork+1, ierr )
1960 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1961 $ work( iwork ), lwork-iwork+1, ierr )
1970 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1971 $ ldvt, work( iu ), ldwrku, cdum, 1,
1972 $ rwork( irwork ), info )
1979 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1980 $ work( iu ), ldwrku, czero, a, lda )
1984 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1997 CALL cgeqrf( m, n, a, lda, work( itau ),
1998 $ work( iwork ), lwork-iwork+1, ierr )
1999 CALL clacpy(
'L', m, n, a, lda, u, ldu )
2005 CALL cungqr( m, m, n, u, ldu, work( itau ),
2006 $ work( iwork ), lwork-iwork+1, ierr )
2010 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2012 $
CALL claset(
'L', n-1, n-1, czero, czero,
2013 $ vt( 2, 1 ), ldvt )
2023 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2024 $ work( itauq ), work( itaup ),
2025 $ work( iwork ), lwork-iwork+1, ierr )
2032 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2033 $ work( itauq ), u, ldu, work( iwork ),
2034 $ lwork-iwork+1, ierr )
2040 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2041 $ work( iwork ), lwork-iwork+1, ierr )
2050 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
2051 $ ldvt, u, ldu, cdum, 1,
2052 $ rwork( irwork ), info )
2076 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2077 $ work( itaup ), work( iwork ), lwork-iwork+1,
2086 CALL clacpy(
'L', m, n, a, lda, u, ldu )
2091 CALL cungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2092 $ work( iwork ), lwork-iwork+1, ierr )
2101 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2102 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2103 $ work( iwork ), lwork-iwork+1, ierr )
2112 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
2113 $ work( iwork ), lwork-iwork+1, ierr )
2122 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
2123 $ work( iwork ), lwork-iwork+1, ierr )
2126 IF( wntuas .OR. wntuo )
2130 IF( wntvas .OR. wntvo )
2134 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2142 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2143 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2145 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2153 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2154 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2164 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2165 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2177 IF( n.GE.mnthr )
THEN
2191 CALL cgelqf( m, n, a, lda, work( itau ), work( iwork ),
2192 $ lwork-iwork+1, ierr )
2196 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2207 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2208 $ work( itaup ), work( iwork ), lwork-iwork+1,
2210 IF( wntuo .OR. wntuas )
THEN
2216 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2217 $ work( iwork ), lwork-iwork+1, ierr )
2221 IF( wntuo .OR. wntuas )
2229 CALL cbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2230 $ a, lda, cdum, 1, rwork( irwork ), info )
2235 $
CALL clacpy(
'F', m, m, a, lda, u, ldu )
2237 ELSE IF( wntvo .AND. wntun )
THEN
2243 IF( lwork.GE.m*m+3*m )
THEN
2248 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2255 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2267 chunk = ( lwork-m*m ) / m
2270 itau = ir + ldwrkr*m
2277 CALL cgelqf( m, n, a, lda, work( itau ),
2278 $ work( iwork ), lwork-iwork+1, ierr )
2282 CALL clacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2283 CALL claset(
'U', m-1, m-1, czero, czero,
2284 $ work( ir+ldwrkr ), ldwrkr )
2290 CALL cunglq( m, n, m, a, lda, work( itau ),
2291 $ work( iwork ), lwork-iwork+1, ierr )
2301 CALL cgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2302 $ work( itauq ), work( itaup ),
2303 $ work( iwork ), lwork-iwork+1, ierr )
2309 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2310 $ work( itaup ), work( iwork ),
2311 $ lwork-iwork+1, ierr )
2319 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2320 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2321 $ rwork( irwork ), info )
2329 DO 30 i = 1, n, chunk
2330 blk = min( n-i+1, chunk )
2331 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2332 $ ldwrkr, a( 1, i ), lda, czero,
2333 $ work( iu ), ldwrku )
2334 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2351 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
2352 $ work( itauq ), work( itaup ),
2353 $ work( iwork ), lwork-iwork+1, ierr )
2359 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
2360 $ work( iwork ), lwork-iwork+1, ierr )
2368 CALL cbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2369 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2373 ELSE IF( wntvo .AND. wntuas )
THEN
2379 IF( lwork.GE.m*m+3*m )
THEN
2384 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2391 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2403 chunk = ( lwork-m*m ) / m
2406 itau = ir + ldwrkr*m
2413 CALL cgelqf( m, n, a, lda, work( itau ),
2414 $ work( iwork ), lwork-iwork+1, ierr )
2418 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2419 CALL claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2426 CALL cunglq( m, n, m, a, lda, work( itau ),
2427 $ work( iwork ), lwork-iwork+1, ierr )
2437 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2438 $ work( itauq ), work( itaup ),
2439 $ work( iwork ), lwork-iwork+1, ierr )
2440 CALL clacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2446 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2447 $ work( itaup ), work( iwork ),
2448 $ lwork-iwork+1, ierr )
2454 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2455 $ work( iwork ), lwork-iwork+1, ierr )
2464 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2465 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2466 $ rwork( irwork ), info )
2474 DO 40 i = 1, n, chunk
2475 blk = min( n-i+1, chunk )
2476 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2477 $ ldwrkr, a( 1, i ), lda, czero,
2478 $ work( iu ), ldwrku )
2479 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2494 CALL cgelqf( m, n, a, lda, work( itau ),
2495 $ work( iwork ), lwork-iwork+1, ierr )
2499 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2500 CALL claset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2507 CALL cunglq( m, n, m, a, lda, work( itau ),
2508 $ work( iwork ), lwork-iwork+1, ierr )
2518 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2519 $ work( itauq ), work( itaup ),
2520 $ work( iwork ), lwork-iwork+1, ierr )
2526 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2527 $ work( itaup ), a, lda, work( iwork ),
2528 $ lwork-iwork+1, ierr )
2534 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2535 $ work( iwork ), lwork-iwork+1, ierr )
2544 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a, lda,
2545 $ u, ldu, cdum, 1, rwork( irwork ), info )
2549 ELSE IF( wntvs )
THEN
2557 IF( lwork.GE.m*m+3*m )
THEN
2562 IF( lwork.GE.wrkbl+lda*m )
THEN
2573 itau = ir + ldwrkr*m
2580 CALL cgelqf( m, n, a, lda, work( itau ),
2581 $ work( iwork ), lwork-iwork+1, ierr )
2585 CALL clacpy(
'L', m, m, a, lda, work( ir ),
2587 CALL claset(
'U', m-1, m-1, czero, czero,
2588 $ work( ir+ldwrkr ), ldwrkr )
2594 CALL cunglq( m, n, m, a, lda, work( itau ),
2595 $ work( iwork ), lwork-iwork+1, ierr )
2605 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2606 $ rwork( ie ), work( itauq ),
2607 $ work( itaup ), work( iwork ),
2608 $ lwork-iwork+1, ierr )
2615 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2616 $ work( itaup ), work( iwork ),
2617 $ lwork-iwork+1, ierr )
2625 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2626 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2627 $ rwork( irwork ), info )
2634 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
2635 $ ldwrkr, a, lda, czero, vt, ldvt )
2648 CALL cgelqf( m, n, a, lda, work( itau ),
2649 $ work( iwork ), lwork-iwork+1, ierr )
2653 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2659 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2660 $ work( iwork ), lwork-iwork+1, ierr )
2668 CALL claset(
'U', m-1, m-1, czero, czero,
2675 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2676 $ work( itauq ), work( itaup ),
2677 $ work( iwork ), lwork-iwork+1, ierr )
2683 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2684 $ work( itaup ), vt, ldvt,
2685 $ work( iwork ), lwork-iwork+1, ierr )
2693 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
2694 $ ldvt, cdum, 1, cdum, 1,
2695 $ rwork( irwork ), info )
2699 ELSE IF( wntuo )
THEN
2705 IF( lwork.GE.2*m*m+3*m )
THEN
2710 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2717 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2732 itau = ir + ldwrkr*m
2739 CALL cgelqf( m, n, a, lda, work( itau ),
2740 $ work( iwork ), lwork-iwork+1, ierr )
2744 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2746 CALL claset(
'U', m-1, m-1, czero, czero,
2747 $ work( iu+ldwrku ), ldwrku )
2753 CALL cunglq( m, n, m, a, lda, work( itau ),
2754 $ work( iwork ), lwork-iwork+1, ierr )
2766 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2767 $ rwork( ie ), work( itauq ),
2768 $ work( itaup ), work( iwork ),
2769 $ lwork-iwork+1, ierr )
2770 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
2771 $ work( ir ), ldwrkr )
2778 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2779 $ work( itaup ), work( iwork ),
2780 $ lwork-iwork+1, ierr )
2786 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2787 $ work( itauq ), work( iwork ),
2788 $ lwork-iwork+1, ierr )
2797 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2798 $ work( iu ), ldwrku, work( ir ),
2799 $ ldwrkr, cdum, 1, rwork( irwork ),
2807 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2808 $ ldwrku, a, lda, czero, vt, ldvt )
2814 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
2828 CALL cgelqf( m, n, a, lda, work( itau ),
2829 $ work( iwork ), lwork-iwork+1, ierr )
2830 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2836 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2837 $ work( iwork ), lwork-iwork+1, ierr )
2845 CALL claset(
'U', m-1, m-1, czero, czero,
2852 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2853 $ work( itauq ), work( itaup ),
2854 $ work( iwork ), lwork-iwork+1, ierr )
2860 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2861 $ work( itaup ), vt, ldvt,
2862 $ work( iwork ), lwork-iwork+1, ierr )
2868 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2869 $ work( iwork ), lwork-iwork+1, ierr )
2878 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
2879 $ ldvt, a, lda, cdum, 1,
2880 $ rwork( irwork ), info )
2884 ELSE IF( wntuas )
THEN
2891 IF( lwork.GE.m*m+3*m )
THEN
2896 IF( lwork.GE.wrkbl+lda*m )
THEN
2907 itau = iu + ldwrku*m
2914 CALL cgelqf( m, n, a, lda, work( itau ),
2915 $ work( iwork ), lwork-iwork+1, ierr )
2919 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2921 CALL claset(
'U', m-1, m-1, czero, czero,
2922 $ work( iu+ldwrku ), ldwrku )
2928 CALL cunglq( m, n, m, a, lda, work( itau ),
2929 $ work( iwork ), lwork-iwork+1, ierr )
2939 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2940 $ rwork( ie ), work( itauq ),
2941 $ work( itaup ), work( iwork ),
2942 $ lwork-iwork+1, ierr )
2943 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
2951 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2952 $ work( itaup ), work( iwork ),
2953 $ lwork-iwork+1, ierr )
2959 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2960 $ work( iwork ), lwork-iwork+1, ierr )
2969 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2970 $ work( iu ), ldwrku, u, ldu, cdum, 1,
2971 $ rwork( irwork ), info )
2978 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2979 $ ldwrku, a, lda, czero, vt, ldvt )
2992 CALL cgelqf( m, n, a, lda, work( itau ),
2993 $ work( iwork ), lwork-iwork+1, ierr )
2994 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3000 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
3001 $ work( iwork ), lwork-iwork+1, ierr )
3005 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3006 CALL claset(
'U', m-1, m-1, czero, czero,
3017 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3018 $ work( itauq ), work( itaup ),
3019 $ work( iwork ), lwork-iwork+1, ierr )
3026 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3027 $ work( itaup ), vt, ldvt,
3028 $ work( iwork ), lwork-iwork+1, ierr )
3034 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3035 $ work( iwork ), lwork-iwork+1, ierr )
3044 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3045 $ ldvt, u, ldu, cdum, 1,
3046 $ rwork( irwork ), info )
3052 ELSE IF( wntva )
THEN
3060 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3065 IF( lwork.GE.wrkbl+lda*m )
THEN
3076 itau = ir + ldwrkr*m
3083 CALL cgelqf( m, n, a, lda, work( itau ),
3084 $ work( iwork ), lwork-iwork+1, ierr )
3085 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3089 CALL clacpy(
'L', m, m, a, lda, work( ir ),
3091 CALL claset(
'U', m-1, m-1, czero, czero,
3092 $ work( ir+ldwrkr ), ldwrkr )
3098 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3099 $ work( iwork ), lwork-iwork+1, ierr )
3109 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
3110 $ rwork( ie ), work( itauq ),
3111 $ work( itaup ), work( iwork ),
3112 $ lwork-iwork+1, ierr )
3119 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
3120 $ work( itaup ), work( iwork ),
3121 $ lwork-iwork+1, ierr )
3129 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3130 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3131 $ rwork( irwork ), info )
3138 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
3139 $ ldwrkr, vt, ldvt, czero, a, lda )
3143 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3156 CALL cgelqf( m, n, a, lda, work( itau ),
3157 $ work( iwork ), lwork-iwork+1, ierr )
3158 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3164 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3165 $ work( iwork ), lwork-iwork+1, ierr )
3173 CALL claset(
'U', m-1, m-1, czero, czero,
3180 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3181 $ work( itauq ), work( itaup ),
3182 $ work( iwork ), lwork-iwork+1, ierr )
3189 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3190 $ work( itaup ), vt, ldvt,
3191 $ work( iwork ), lwork-iwork+1, ierr )
3199 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
3200 $ ldvt, cdum, 1, cdum, 1,
3201 $ rwork( irwork ), info )
3205 ELSE IF( wntuo )
THEN
3211 IF( lwork.GE.2*m*m+max( n+m, 3*m ) )
THEN
3216 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3223 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3238 itau = ir + ldwrkr*m
3245 CALL cgelqf( m, n, a, lda, work( itau ),
3246 $ work( iwork ), lwork-iwork+1, ierr )
3247 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3253 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3254 $ work( iwork ), lwork-iwork+1, ierr )
3258 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3260 CALL claset(
'U', m-1, m-1, czero, czero,
3261 $ work( iu+ldwrku ), ldwrku )
3273 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3274 $ rwork( ie ), work( itauq ),
3275 $ work( itaup ), work( iwork ),
3276 $ lwork-iwork+1, ierr )
3277 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
3278 $ work( ir ), ldwrkr )
3285 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3286 $ work( itaup ), work( iwork ),
3287 $ lwork-iwork+1, ierr )
3293 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3294 $ work( itauq ), work( iwork ),
3295 $ lwork-iwork+1, ierr )
3304 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3305 $ work( iu ), ldwrku, work( ir ),
3306 $ ldwrkr, cdum, 1, rwork( irwork ),
3314 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3315 $ ldwrku, vt, ldvt, czero, a, lda )
3319 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3323 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
3337 CALL cgelqf( m, n, a, lda, work( itau ),
3338 $ work( iwork ), lwork-iwork+1, ierr )
3339 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3345 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3346 $ work( iwork ), lwork-iwork+1, ierr )
3354 CALL claset(
'U', m-1, m-1, czero, czero,
3361 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3362 $ work( itauq ), work( itaup ),
3363 $ work( iwork ), lwork-iwork+1, ierr )
3370 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3371 $ work( itaup ), vt, ldvt,
3372 $ work( iwork ), lwork-iwork+1, ierr )
3378 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
3379 $ work( iwork ), lwork-iwork+1, ierr )
3388 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3389 $ ldvt, a, lda, cdum, 1,
3390 $ rwork( irwork ), info )
3394 ELSE IF( wntuas )
THEN
3401 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3406 IF( lwork.GE.wrkbl+lda*m )
THEN
3417 itau = iu + ldwrku*m
3424 CALL cgelqf( m, n, a, lda, work( itau ),
3425 $ work( iwork ), lwork-iwork+1, ierr )
3426 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3432 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3433 $ work( iwork ), lwork-iwork+1, ierr )
3437 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3439 CALL claset(
'U', m-1, m-1, czero, czero,
3440 $ work( iu+ldwrku ), ldwrku )
3450 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3451 $ rwork( ie ), work( itauq ),
3452 $ work( itaup ), work( iwork ),
3453 $ lwork-iwork+1, ierr )
3454 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
3461 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3462 $ work( itaup ), work( iwork ),
3463 $ lwork-iwork+1, ierr )
3469 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3470 $ work( iwork ), lwork-iwork+1, ierr )
3479 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3480 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3481 $ rwork( irwork ), info )
3488 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3489 $ ldwrku, vt, ldvt, czero, a, lda )
3493 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3506 CALL cgelqf( m, n, a, lda, work( itau ),
3507 $ work( iwork ), lwork-iwork+1, ierr )
3508 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3514 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3515 $ work( iwork ), lwork-iwork+1, ierr )
3519 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3520 CALL claset(
'U', m-1, m-1, czero, czero,
3531 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3532 $ work( itauq ), work( itaup ),
3533 $ work( iwork ), lwork-iwork+1, ierr )
3540 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3541 $ work( itaup ), vt, ldvt,
3542 $ work( iwork ), lwork-iwork+1, ierr )
3548 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3549 $ work( iwork ), lwork-iwork+1, ierr )
3558 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3559 $ ldvt, u, ldu, cdum, 1,
3560 $ rwork( irwork ), info )
3584 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3585 $ work( itaup ), work( iwork ), lwork-iwork+1,
3594 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3595 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3596 $ work( iwork ), lwork-iwork+1, ierr )
3605 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3610 CALL cungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3611 $ work( iwork ), lwork-iwork+1, ierr )
3620 CALL cungbr(
'Q', m, m, n, a, lda, work( itauq ),
3621 $ work( iwork ), lwork-iwork+1, ierr )
3630 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
3631 $ work( iwork ), lwork-iwork+1, ierr )
3634 IF( wntuas .OR. wntuo )
3638 IF( wntvas .OR. wntvo )
3642 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3650 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3651 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3653 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3661 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3662 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3672 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3673 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3683 IF( iscl.EQ.1 )
THEN
3684 IF( anrm.GT.bignum )
3685 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3687 IF( info.NE.0 .AND. anrm.GT.bignum )
3688 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3689 $ rwork( ie ), minmn, ierr )
3690 IF( anrm.LT.smlnum )
3691 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3693 IF( info.NE.0 .AND. anrm.LT.smlnum )
3694 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3695 $ rwork( ie ), minmn, ierr )
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
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 cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
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.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
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 cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
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.
logical function lsame(CA, CB)
LSAME
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
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