214 SUBROUTINE zgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
215 $ vt, ldvt, work, lwork, rwork, info )
223 CHARACTER jobu, jobvt
224 INTEGER info, lda, ldu, ldvt, lwork, m, n
227 DOUBLE PRECISION rwork( * ), s( * )
228 COMPLEX*16 a( lda, * ), u( ldu, * ), vt( ldvt, * ),
235 COMPLEX*16 czero, cone
236 parameter( czero = ( 0.0d0, 0.0d0 ),
237 $ cone = ( 1.0d0, 0.0d0 ) )
238 DOUBLE PRECISION zero, one
239 parameter( zero = 0.0d0, one = 1.0d0 )
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_zgeqrf, lwork_zungqr_n, lwork_zungqr_m,
249 $ lwork_zgebrd, lwork_zungbr_p, lwork_zungbr_q,
250 $ lwork_zgelqf, lwork_zunglq_n, lwork_zunglq_m
251 DOUBLE PRECISION anrm, bignum, eps, smlnum
254 DOUBLE PRECISION dum( 1 )
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,
'ZGESVD', jobu // jobvt, m, n, 0, 0 )
324 CALL
zgeqrf( m, n, a, lda, dum(1), dum(1), -1, ierr )
327 CALL
zungqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
328 lwork_zungqr_n=dum(1)
329 CALL
zungqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
330 lwork_zungqr_m=dum(1)
332 CALL
zgebrd( n, n, a, lda, s, dum(1), dum(1),
333 $ dum(1), dum(1), -1, ierr )
336 CALL
zungbr(
'P', n, n, n, a, lda, dum(1),
338 lwork_zungbr_p=dum(1)
339 CALL
zungbr(
'Q', n, n, n, a, lda, dum(1),
341 lwork_zungbr_q=dum(1)
343 IF( m.GE.mnthr )
THEN
348 maxwrk = n + lwork_zgeqrf
349 maxwrk = max( maxwrk, 2*n+lwork_zgebrd )
350 IF( wntvo .OR. wntvas )
351 $ maxwrk = max( maxwrk, 2*n+lwork_zungbr_p )
353 ELSE IF( wntuo .AND. wntvn )
THEN
357 wrkbl = n + lwork_zgeqrf
358 wrkbl = max( wrkbl, n+lwork_zungqr_n )
359 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
360 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
361 maxwrk = max( n*n+wrkbl, n*n+m*n )
363 ELSE IF( wntuo .AND. wntvas )
THEN
368 wrkbl = n + lwork_zgeqrf
369 wrkbl = max( wrkbl, n+lwork_zungqr_n )
370 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
371 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
372 wrkbl = max( wrkbl, 2*n+lwork_zungbr_p )
373 maxwrk = max( n*n+wrkbl, n*n+m*n )
375 ELSE IF( wntus .AND. wntvn )
THEN
379 wrkbl = n + lwork_zgeqrf
380 wrkbl = max( wrkbl, n+lwork_zungqr_n )
381 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
382 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
385 ELSE IF( wntus .AND. wntvo )
THEN
389 wrkbl = n + lwork_zgeqrf
390 wrkbl = max( wrkbl, n+lwork_zungqr_n )
391 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
392 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
393 wrkbl = max( wrkbl, 2*n+lwork_zungbr_p )
394 maxwrk = 2*n*n + wrkbl
396 ELSE IF( wntus .AND. wntvas )
THEN
401 wrkbl = n + lwork_zgeqrf
402 wrkbl = max( wrkbl, n+lwork_zungqr_n )
403 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
404 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
405 wrkbl = max( wrkbl, 2*n+lwork_zungbr_p )
408 ELSE IF( wntua .AND. wntvn )
THEN
412 wrkbl = n + lwork_zgeqrf
413 wrkbl = max( wrkbl, n+lwork_zungqr_m )
414 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
415 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
418 ELSE IF( wntua .AND. wntvo )
THEN
422 wrkbl = n + lwork_zgeqrf
423 wrkbl = max( wrkbl, n+lwork_zungqr_m )
424 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
425 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
426 wrkbl = max( wrkbl, 2*n+lwork_zungbr_p )
427 maxwrk = 2*n*n + wrkbl
429 ELSE IF( wntua .AND. wntvas )
THEN
434 wrkbl = n + lwork_zgeqrf
435 wrkbl = max( wrkbl, n+lwork_zungqr_m )
436 wrkbl = max( wrkbl, 2*n+lwork_zgebrd )
437 wrkbl = max( wrkbl, 2*n+lwork_zungbr_q )
438 wrkbl = max( wrkbl, 2*n+lwork_zungbr_p )
446 CALL
zgebrd( m, n, a, lda, s, dum(1), dum(1),
447 $ dum(1), dum(1), -1, ierr )
449 maxwrk = 2*n + lwork_zgebrd
450 IF( wntus .OR. wntuo )
THEN
451 CALL
zungbr(
'Q', m, n, n, a, lda, dum(1),
453 lwork_zungbr_q=dum(1)
454 maxwrk = max( maxwrk, 2*n+lwork_zungbr_q )
457 CALL
zungbr(
'Q', m, m, n, a, lda, dum(1),
459 lwork_zungbr_q=dum(1)
460 maxwrk = max( maxwrk, 2*n+lwork_zungbr_q )
462 IF( .NOT.wntvn )
THEN
463 maxwrk = max( maxwrk, 2*n+lwork_zungbr_p )
467 ELSE IF( minmn.GT.0 )
THEN
471 mnthr =
ilaenv( 6,
'ZGESVD', jobu // jobvt, m, n, 0, 0 )
473 CALL
zgelqf( m, n, a, lda, dum(1), dum(1), -1, ierr )
476 CALL
zunglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
477 lwork_zunglq_n=dum(1)
478 CALL
zunglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
479 lwork_zunglq_m=dum(1)
481 CALL
zgebrd( m, m, a, lda, s, dum(1), dum(1),
482 $ dum(1), dum(1), -1, ierr )
485 CALL
zungbr(
'P', m, m, m, a, n, dum(1),
487 lwork_zungbr_p=dum(1)
489 CALL
zungbr(
'Q', m, m, m, a, n, dum(1),
491 lwork_zungbr_q=dum(1)
492 IF( n.GE.mnthr )
THEN
497 maxwrk = m + lwork_zgelqf
498 maxwrk = max( maxwrk, 2*m+lwork_zgebrd )
499 IF( wntuo .OR. wntuas )
500 $ maxwrk = max( maxwrk, 2*m+lwork_zungbr_q )
502 ELSE IF( wntvo .AND. wntun )
THEN
506 wrkbl = m + lwork_zgelqf
507 wrkbl = max( wrkbl, m+lwork_zunglq_m )
508 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
509 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
510 maxwrk = max( m*m+wrkbl, m*m+m*n )
512 ELSE IF( wntvo .AND. wntuas )
THEN
517 wrkbl = m + lwork_zgelqf
518 wrkbl = max( wrkbl, m+lwork_zunglq_m )
519 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
520 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
521 wrkbl = max( wrkbl, 2*m+lwork_zungbr_q )
522 maxwrk = max( m*m+wrkbl, m*m+m*n )
524 ELSE IF( wntvs .AND. wntun )
THEN
528 wrkbl = m + lwork_zgelqf
529 wrkbl = max( wrkbl, m+lwork_zunglq_m )
530 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
531 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
534 ELSE IF( wntvs .AND. wntuo )
THEN
538 wrkbl = m + lwork_zgelqf
539 wrkbl = max( wrkbl, m+lwork_zunglq_m )
540 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
541 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
542 wrkbl = max( wrkbl, 2*m+lwork_zungbr_q )
543 maxwrk = 2*m*m + wrkbl
545 ELSE IF( wntvs .AND. wntuas )
THEN
550 wrkbl = m + lwork_zgelqf
551 wrkbl = max( wrkbl, m+lwork_zunglq_m )
552 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
553 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
554 wrkbl = max( wrkbl, 2*m+lwork_zungbr_q )
557 ELSE IF( wntva .AND. wntun )
THEN
561 wrkbl = m + lwork_zgelqf
562 wrkbl = max( wrkbl, m+lwork_zunglq_n )
563 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
564 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
567 ELSE IF( wntva .AND. wntuo )
THEN
571 wrkbl = m + lwork_zgelqf
572 wrkbl = max( wrkbl, m+lwork_zunglq_n )
573 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
574 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
575 wrkbl = max( wrkbl, 2*m+lwork_zungbr_q )
576 maxwrk = 2*m*m + wrkbl
578 ELSE IF( wntva .AND. wntuas )
THEN
583 wrkbl = m + lwork_zgelqf
584 wrkbl = max( wrkbl, m+lwork_zunglq_n )
585 wrkbl = max( wrkbl, 2*m+lwork_zgebrd )
586 wrkbl = max( wrkbl, 2*m+lwork_zungbr_p )
587 wrkbl = max( wrkbl, 2*m+lwork_zungbr_q )
595 CALL
zgebrd( m, n, a, lda, s, dum(1), dum(1),
596 $ dum(1), dum(1), -1, ierr )
598 maxwrk = 2*m + lwork_zgebrd
599 IF( wntvs .OR. wntvo )
THEN
601 CALL
zungbr(
'P', m, n, m, a, n, dum(1),
603 lwork_zungbr_p=dum(1)
604 maxwrk = max( maxwrk, 2*m+lwork_zungbr_p )
607 CALL
zungbr(
'P', n, n, m, a, n, dum(1),
609 lwork_zungbr_p=dum(1)
610 maxwrk = max( maxwrk, 2*m+lwork_zungbr_p )
612 IF( .NOT.wntun )
THEN
613 maxwrk = max( maxwrk, 2*m+lwork_zungbr_q )
618 maxwrk = max( maxwrk, minwrk )
621 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
627 CALL
xerbla(
'ZGESVD', -info )
629 ELSE IF( lquery )
THEN
635 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
642 smlnum = sqrt(
dlamch(
'S' ) ) / eps
643 bignum = one / smlnum
647 anrm =
zlange(
'M', m, n, a, lda, dum )
649 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
651 CALL
zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
652 ELSE IF( anrm.GT.bignum )
THEN
654 CALL
zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
663 IF( m.GE.mnthr )
THEN
677 CALL
zgeqrf( m, n, a, lda, work( itau ), work( iwork ),
678 $ lwork-iwork+1, ierr )
682 CALL
zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
693 CALL
zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
694 $ work( itaup ), work( iwork ), lwork-iwork+1,
697 IF( wntvo .OR. wntvas )
THEN
703 CALL
zungbr(
'P', n, n, n, a, lda, work( itaup ),
704 $ work( iwork ), lwork-iwork+1, ierr )
714 CALL
zbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
715 $ cdum, 1, cdum, 1, rwork( irwork ), info )
720 $ CALL
zlacpy(
'F', n, n, a, lda, vt, ldvt )
722 ELSE IF( wntuo .AND. wntvn )
THEN
728 IF( lwork.GE.n*n+3*n )
THEN
733 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
739 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
749 ldwrku = ( lwork-n*n ) / n
759 CALL
zgeqrf( m, n, a, lda, work( itau ),
760 $ work( iwork ), lwork-iwork+1, ierr )
764 CALL
zlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
765 CALL
zlaset(
'L', n-1, n-1, czero, czero,
766 $ work( ir+1 ), ldwrkr )
772 CALL
zungqr( m, n, n, a, lda, work( itau ),
773 $ work( iwork ), lwork-iwork+1, ierr )
783 CALL
zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
784 $ work( itauq ), work( itaup ),
785 $ work( iwork ), lwork-iwork+1, ierr )
791 CALL
zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
792 $ work( itauq ), work( iwork ),
793 $ lwork-iwork+1, ierr )
801 CALL
zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
802 $ work( ir ), ldwrkr, cdum, 1,
803 $ rwork( irwork ), info )
811 DO 10 i = 1, m, ldwrku
812 chunk = min( m-i+1, ldwrku )
813 CALL
zgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
814 $ lda, work( ir ), ldwrkr, czero,
815 $ work( iu ), ldwrku )
816 CALL
zlacpy(
'F', chunk, n, work( iu ), ldwrku,
833 CALL
zgebrd( m, n, a, lda, s, rwork( ie ),
834 $ work( itauq ), work( itaup ),
835 $ work( iwork ), lwork-iwork+1, ierr )
841 CALL
zungbr(
'Q', m, n, n, a, lda, work( itauq ),
842 $ work( iwork ), lwork-iwork+1, ierr )
850 CALL
zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
851 $ a, lda, cdum, 1, rwork( irwork ), info )
855 ELSE IF( wntuo .AND. wntvas )
THEN
861 IF( lwork.GE.n*n+3*n )
THEN
866 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
872 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
882 ldwrku = ( lwork-n*n ) / n
892 CALL
zgeqrf( m, n, a, lda, work( itau ),
893 $ work( iwork ), lwork-iwork+1, ierr )
897 CALL
zlacpy(
'U', n, n, a, lda, vt, ldvt )
899 $ CALL
zlaset(
'L', n-1, n-1, czero, czero,
906 CALL
zungqr( m, n, n, a, lda, work( itau ),
907 $ work( iwork ), lwork-iwork+1, ierr )
917 CALL
zgebrd( n, n, vt, ldvt, s, rwork( ie ),
918 $ work( itauq ), work( itaup ),
919 $ work( iwork ), lwork-iwork+1, ierr )
920 CALL
zlacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
926 CALL
zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
927 $ work( itauq ), work( iwork ),
928 $ lwork-iwork+1, ierr )
934 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
935 $ work( iwork ), lwork-iwork+1, ierr )
944 CALL
zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
945 $ ldvt, work( ir ), ldwrkr, cdum, 1,
946 $ rwork( irwork ), info )
954 DO 20 i = 1, m, ldwrku
955 chunk = min( m-i+1, ldwrku )
956 CALL
zgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
957 $ lda, work( ir ), ldwrkr, czero,
958 $ work( iu ), ldwrku )
959 CALL
zlacpy(
'F', chunk, n, work( iu ), ldwrku,
974 CALL
zgeqrf( m, n, a, lda, work( itau ),
975 $ work( iwork ), lwork-iwork+1, ierr )
979 CALL
zlacpy(
'U', n, n, a, lda, vt, ldvt )
981 $ CALL
zlaset(
'L', n-1, n-1, czero, czero,
988 CALL
zungqr( m, n, n, a, lda, work( itau ),
989 $ work( iwork ), lwork-iwork+1, ierr )
999 CALL
zgebrd( n, n, vt, ldvt, s, rwork( ie ),
1000 $ work( itauq ), work( itaup ),
1001 $ work( iwork ), lwork-iwork+1, ierr )
1007 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1008 $ work( itauq ), a, lda, work( iwork ),
1009 $ lwork-iwork+1, ierr )
1015 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1016 $ work( iwork ), lwork-iwork+1, ierr )
1025 CALL
zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1026 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1031 ELSE IF( wntus )
THEN
1039 IF( lwork.GE.n*n+3*n )
THEN
1044 IF( lwork.GE.wrkbl+lda*n )
THEN
1055 itau = ir + ldwrkr*n
1062 CALL
zgeqrf( m, n, a, lda, work( itau ),
1063 $ work( iwork ), lwork-iwork+1, ierr )
1067 CALL
zlacpy(
'U', n, n, a, lda, work( ir ),
1069 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1070 $ work( ir+1 ), ldwrkr )
1076 CALL
zungqr( m, n, n, a, lda, work( itau ),
1077 $ work( iwork ), lwork-iwork+1, ierr )
1087 CALL
zgebrd( n, n, work( ir ), ldwrkr, s,
1088 $ rwork( ie ), work( itauq ),
1089 $ work( itaup ), work( iwork ),
1090 $ lwork-iwork+1, ierr )
1096 CALL
zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1097 $ work( itauq ), work( iwork ),
1098 $ lwork-iwork+1, ierr )
1106 CALL
zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1107 $ 1, work( ir ), ldwrkr, cdum, 1,
1108 $ rwork( irwork ), info )
1115 CALL
zgemm(
'N',
'N', m, n, n, cone, a, lda,
1116 $ work( ir ), ldwrkr, czero, u, ldu )
1129 CALL
zgeqrf( m, n, a, lda, work( itau ),
1130 $ work( iwork ), lwork-iwork+1, ierr )
1131 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1137 CALL
zungqr( m, n, n, u, ldu, work( itau ),
1138 $ work( iwork ), lwork-iwork+1, ierr )
1146 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1153 CALL
zgebrd( n, n, a, lda, s, rwork( ie ),
1154 $ work( itauq ), work( itaup ),
1155 $ work( iwork ), lwork-iwork+1, ierr )
1161 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1162 $ work( itauq ), u, ldu, work( iwork ),
1163 $ lwork-iwork+1, ierr )
1171 CALL
zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1172 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1177 ELSE IF( wntvo )
THEN
1183 IF( lwork.GE.2*n*n+3*n )
THEN
1188 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1195 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1210 itau = ir + ldwrkr*n
1217 CALL
zgeqrf( m, n, a, lda, work( itau ),
1218 $ work( iwork ), lwork-iwork+1, ierr )
1222 CALL
zlacpy(
'U', n, n, a, lda, work( iu ),
1224 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1225 $ work( iu+1 ), ldwrku )
1231 CALL
zungqr( m, n, n, a, lda, work( itau ),
1232 $ work( iwork ), lwork-iwork+1, ierr )
1244 CALL
zgebrd( n, n, work( iu ), ldwrku, s,
1245 $ rwork( ie ), work( itauq ),
1246 $ work( itaup ), work( iwork ),
1247 $ lwork-iwork+1, ierr )
1248 CALL
zlacpy(
'U', n, n, work( iu ), ldwrku,
1249 $ work( ir ), ldwrkr )
1255 CALL
zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1256 $ work( itauq ), work( iwork ),
1257 $ lwork-iwork+1, ierr )
1264 CALL
zungbr(
'P', n, n, n, work( ir ), ldwrkr,
1265 $ work( itaup ), work( iwork ),
1266 $ lwork-iwork+1, ierr )
1275 CALL
zbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1276 $ work( ir ), ldwrkr, work( iu ),
1277 $ ldwrku, cdum, 1, rwork( irwork ),
1285 CALL
zgemm(
'N',
'N', m, n, n, cone, a, lda,
1286 $ work( iu ), ldwrku, czero, u, ldu )
1292 CALL
zlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1306 CALL
zgeqrf( m, n, a, lda, work( itau ),
1307 $ work( iwork ), lwork-iwork+1, ierr )
1308 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1314 CALL
zungqr( m, n, n, u, ldu, work( itau ),
1315 $ work( iwork ), lwork-iwork+1, ierr )
1323 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1330 CALL
zgebrd( n, n, a, lda, s, rwork( ie ),
1331 $ work( itauq ), work( itaup ),
1332 $ work( iwork ), lwork-iwork+1, ierr )
1338 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1339 $ work( itauq ), u, ldu, work( iwork ),
1340 $ lwork-iwork+1, ierr )
1346 CALL
zungbr(
'P', n, n, n, a, lda, work( itaup ),
1347 $ work( iwork ), lwork-iwork+1, ierr )
1356 CALL
zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1357 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1362 ELSE IF( wntvas )
THEN
1369 IF( lwork.GE.n*n+3*n )
THEN
1374 IF( lwork.GE.wrkbl+lda*n )
THEN
1385 itau = iu + ldwrku*n
1392 CALL
zgeqrf( m, n, a, lda, work( itau ),
1393 $ work( iwork ), lwork-iwork+1, ierr )
1397 CALL
zlacpy(
'U', n, n, a, lda, work( iu ),
1399 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1400 $ work( iu+1 ), ldwrku )
1406 CALL
zungqr( m, n, n, a, lda, work( itau ),
1407 $ work( iwork ), lwork-iwork+1, ierr )
1417 CALL
zgebrd( n, n, work( iu ), ldwrku, s,
1418 $ rwork( ie ), work( itauq ),
1419 $ work( itaup ), work( iwork ),
1420 $ lwork-iwork+1, ierr )
1421 CALL
zlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1428 CALL
zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1429 $ work( itauq ), work( iwork ),
1430 $ lwork-iwork+1, ierr )
1437 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1438 $ work( iwork ), lwork-iwork+1, ierr )
1447 CALL
zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1448 $ ldvt, work( iu ), ldwrku, cdum, 1,
1449 $ rwork( irwork ), info )
1456 CALL
zgemm(
'N',
'N', m, n, n, cone, a, lda,
1457 $ work( iu ), ldwrku, czero, u, ldu )
1470 CALL
zgeqrf( m, n, a, lda, work( itau ),
1471 $ work( iwork ), lwork-iwork+1, ierr )
1472 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1478 CALL
zungqr( m, n, n, u, ldu, work( itau ),
1479 $ work( iwork ), lwork-iwork+1, ierr )
1483 CALL
zlacpy(
'U', n, n, a, lda, vt, ldvt )
1485 $ CALL
zlaset(
'L', n-1, n-1, czero, czero,
1486 $ vt( 2, 1 ), ldvt )
1496 CALL
zgebrd( n, n, vt, ldvt, s, rwork( ie ),
1497 $ work( itauq ), work( itaup ),
1498 $ work( iwork ), lwork-iwork+1, ierr )
1505 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1506 $ work( itauq ), u, ldu, work( iwork ),
1507 $ lwork-iwork+1, ierr )
1513 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1514 $ work( iwork ), lwork-iwork+1, ierr )
1523 CALL
zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1524 $ ldvt, u, ldu, cdum, 1,
1525 $ rwork( irwork ), info )
1531 ELSE IF( wntua )
THEN
1539 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1544 IF( lwork.GE.wrkbl+lda*n )
THEN
1555 itau = ir + ldwrkr*n
1562 CALL
zgeqrf( m, n, a, lda, work( itau ),
1563 $ work( iwork ), lwork-iwork+1, ierr )
1564 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1568 CALL
zlacpy(
'U', n, n, a, lda, work( ir ),
1570 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1571 $ work( ir+1 ), ldwrkr )
1577 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1578 $ work( iwork ), lwork-iwork+1, ierr )
1588 CALL
zgebrd( n, n, work( ir ), ldwrkr, s,
1589 $ rwork( ie ), work( itauq ),
1590 $ work( itaup ), work( iwork ),
1591 $ lwork-iwork+1, ierr )
1597 CALL
zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1598 $ work( itauq ), work( iwork ),
1599 $ lwork-iwork+1, ierr )
1607 CALL
zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1608 $ 1, work( ir ), ldwrkr, cdum, 1,
1609 $ rwork( irwork ), info )
1616 CALL
zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1617 $ work( ir ), ldwrkr, czero, a, lda )
1621 CALL
zlacpy(
'F', m, n, a, lda, u, ldu )
1634 CALL
zgeqrf( m, n, a, lda, work( itau ),
1635 $ work( iwork ), lwork-iwork+1, ierr )
1636 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1642 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1643 $ work( iwork ), lwork-iwork+1, ierr )
1651 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1658 CALL
zgebrd( n, n, a, lda, s, rwork( ie ),
1659 $ work( itauq ), work( itaup ),
1660 $ work( iwork ), lwork-iwork+1, ierr )
1667 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1668 $ work( itauq ), u, ldu, work( iwork ),
1669 $ lwork-iwork+1, ierr )
1677 CALL
zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1678 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1683 ELSE IF( wntvo )
THEN
1689 IF( lwork.GE.2*n*n+max( n+m, 3*n ) )
THEN
1694 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1701 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1716 itau = ir + ldwrkr*n
1723 CALL
zgeqrf( m, n, a, lda, work( itau ),
1724 $ work( iwork ), lwork-iwork+1, ierr )
1725 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1731 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1732 $ work( iwork ), lwork-iwork+1, ierr )
1736 CALL
zlacpy(
'U', n, n, a, lda, work( iu ),
1738 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1739 $ work( iu+1 ), ldwrku )
1751 CALL
zgebrd( n, n, work( iu ), ldwrku, s,
1752 $ rwork( ie ), work( itauq ),
1753 $ work( itaup ), work( iwork ),
1754 $ lwork-iwork+1, ierr )
1755 CALL
zlacpy(
'U', n, n, work( iu ), ldwrku,
1756 $ work( ir ), ldwrkr )
1762 CALL
zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1763 $ work( itauq ), work( iwork ),
1764 $ lwork-iwork+1, ierr )
1771 CALL
zungbr(
'P', n, n, n, work( ir ), ldwrkr,
1772 $ work( itaup ), work( iwork ),
1773 $ lwork-iwork+1, ierr )
1782 CALL
zbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1783 $ work( ir ), ldwrkr, work( iu ),
1784 $ ldwrku, cdum, 1, rwork( irwork ),
1792 CALL
zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1793 $ work( iu ), ldwrku, czero, a, lda )
1797 CALL
zlacpy(
'F', m, n, a, lda, u, ldu )
1801 CALL
zlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1815 CALL
zgeqrf( m, n, a, lda, work( itau ),
1816 $ work( iwork ), lwork-iwork+1, ierr )
1817 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1823 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1824 $ work( iwork ), lwork-iwork+1, ierr )
1832 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1839 CALL
zgebrd( n, n, a, lda, s, rwork( ie ),
1840 $ work( itauq ), work( itaup ),
1841 $ work( iwork ), lwork-iwork+1, ierr )
1848 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1849 $ work( itauq ), u, ldu, work( iwork ),
1850 $ lwork-iwork+1, ierr )
1856 CALL
zungbr(
'P', n, n, n, a, lda, work( itaup ),
1857 $ work( iwork ), lwork-iwork+1, ierr )
1866 CALL
zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1867 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1872 ELSE IF( wntvas )
THEN
1879 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1884 IF( lwork.GE.wrkbl+lda*n )
THEN
1895 itau = iu + ldwrku*n
1902 CALL
zgeqrf( m, n, a, lda, work( itau ),
1903 $ work( iwork ), lwork-iwork+1, ierr )
1904 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1910 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1911 $ work( iwork ), lwork-iwork+1, ierr )
1915 CALL
zlacpy(
'U', n, n, a, lda, work( iu ),
1917 CALL
zlaset(
'L', n-1, n-1, czero, czero,
1918 $ work( iu+1 ), ldwrku )
1928 CALL
zgebrd( n, n, work( iu ), ldwrku, s,
1929 $ rwork( ie ), work( itauq ),
1930 $ work( itaup ), work( iwork ),
1931 $ lwork-iwork+1, ierr )
1932 CALL
zlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1939 CALL
zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1940 $ work( itauq ), work( iwork ),
1941 $ lwork-iwork+1, ierr )
1948 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1949 $ work( iwork ), lwork-iwork+1, ierr )
1958 CALL
zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1959 $ ldvt, work( iu ), ldwrku, cdum, 1,
1960 $ rwork( irwork ), info )
1967 CALL
zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1968 $ work( iu ), ldwrku, czero, a, lda )
1972 CALL
zlacpy(
'F', m, n, a, lda, u, ldu )
1985 CALL
zgeqrf( m, n, a, lda, work( itau ),
1986 $ work( iwork ), lwork-iwork+1, ierr )
1987 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
1993 CALL
zungqr( m, m, n, u, ldu, work( itau ),
1994 $ work( iwork ), lwork-iwork+1, ierr )
1998 CALL
zlacpy(
'U', n, n, a, lda, vt, ldvt )
2000 $ CALL
zlaset(
'L', n-1, n-1, czero, czero,
2001 $ vt( 2, 1 ), ldvt )
2011 CALL
zgebrd( n, n, vt, ldvt, s, rwork( ie ),
2012 $ work( itauq ), work( itaup ),
2013 $ work( iwork ), lwork-iwork+1, ierr )
2020 CALL
zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2021 $ work( itauq ), u, ldu, work( iwork ),
2022 $ lwork-iwork+1, ierr )
2028 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2029 $ work( iwork ), lwork-iwork+1, ierr )
2038 CALL
zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
2039 $ ldvt, u, ldu, cdum, 1,
2040 $ rwork( irwork ), info )
2064 CALL
zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2065 $ work( itaup ), work( iwork ), lwork-iwork+1,
2074 CALL
zlacpy(
'L', m, n, a, lda, u, ldu )
2079 CALL
zungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2080 $ work( iwork ), lwork-iwork+1, ierr )
2089 CALL
zlacpy(
'U', n, n, a, lda, vt, ldvt )
2090 CALL
zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2091 $ work( iwork ), lwork-iwork+1, ierr )
2100 CALL
zungbr(
'Q', m, n, n, a, lda, work( itauq ),
2101 $ work( iwork ), lwork-iwork+1, ierr )
2110 CALL
zungbr(
'P', n, n, n, a, lda, work( itaup ),
2111 $ work( iwork ), lwork-iwork+1, ierr )
2114 IF( wntuas .OR. wntuo )
2118 IF( wntvas .OR. wntvo )
2122 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2130 CALL
zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2131 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2133 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2141 CALL
zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2142 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2152 CALL
zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2153 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2165 IF( n.GE.mnthr )
THEN
2179 CALL
zgelqf( m, n, a, lda, work( itau ), work( iwork ),
2180 $ lwork-iwork+1, ierr )
2184 CALL
zlaset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2195 CALL
zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2196 $ work( itaup ), work( iwork ), lwork-iwork+1,
2198 IF( wntuo .OR. wntuas )
THEN
2204 CALL
zungbr(
'Q', m, m, m, a, lda, work( itauq ),
2205 $ work( iwork ), lwork-iwork+1, ierr )
2209 IF( wntuo .OR. wntuas )
2217 CALL
zbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2218 $ a, lda, cdum, 1, rwork( irwork ), info )
2223 $ CALL
zlacpy(
'F', m, m, a, lda, u, ldu )
2225 ELSE IF( wntvo .AND. wntun )
THEN
2231 IF( lwork.GE.m*m+3*m )
THEN
2236 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2243 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2255 chunk = ( lwork-m*m ) / m
2258 itau = ir + ldwrkr*m
2265 CALL
zgelqf( m, n, a, lda, work( itau ),
2266 $ work( iwork ), lwork-iwork+1, ierr )
2270 CALL
zlacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2271 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2272 $ work( ir+ldwrkr ), ldwrkr )
2278 CALL
zunglq( m, n, m, a, lda, work( itau ),
2279 $ work( iwork ), lwork-iwork+1, ierr )
2289 CALL
zgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2290 $ work( itauq ), work( itaup ),
2291 $ work( iwork ), lwork-iwork+1, ierr )
2297 CALL
zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2298 $ work( itaup ), work( iwork ),
2299 $ lwork-iwork+1, ierr )
2307 CALL
zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2308 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2309 $ rwork( irwork ), info )
2317 DO 30 i = 1, n, chunk
2318 blk = min( n-i+1, chunk )
2319 CALL
zgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2320 $ ldwrkr, a( 1, i ), lda, czero,
2321 $ work( iu ), ldwrku )
2322 CALL
zlacpy(
'F', m, blk, work( iu ), ldwrku,
2339 CALL
zgebrd( m, n, a, lda, s, rwork( ie ),
2340 $ work( itauq ), work( itaup ),
2341 $ work( iwork ), lwork-iwork+1, ierr )
2347 CALL
zungbr(
'P', m, n, m, a, lda, work( itaup ),
2348 $ work( iwork ), lwork-iwork+1, ierr )
2356 CALL
zbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2357 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2361 ELSE IF( wntvo .AND. wntuas )
THEN
2367 IF( lwork.GE.m*m+3*m )
THEN
2372 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2379 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2391 chunk = ( lwork-m*m ) / m
2394 itau = ir + ldwrkr*m
2401 CALL
zgelqf( m, n, a, lda, work( itau ),
2402 $ work( iwork ), lwork-iwork+1, ierr )
2406 CALL
zlacpy(
'L', m, m, a, lda, u, ldu )
2407 CALL
zlaset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2414 CALL
zunglq( m, n, m, a, lda, work( itau ),
2415 $ work( iwork ), lwork-iwork+1, ierr )
2425 CALL
zgebrd( m, m, u, ldu, s, rwork( ie ),
2426 $ work( itauq ), work( itaup ),
2427 $ work( iwork ), lwork-iwork+1, ierr )
2428 CALL
zlacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2434 CALL
zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2435 $ work( itaup ), work( iwork ),
2436 $ lwork-iwork+1, ierr )
2442 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2443 $ work( iwork ), lwork-iwork+1, ierr )
2452 CALL
zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2453 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2454 $ rwork( irwork ), info )
2462 DO 40 i = 1, n, chunk
2463 blk = min( n-i+1, chunk )
2464 CALL
zgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2465 $ ldwrkr, a( 1, i ), lda, czero,
2466 $ work( iu ), ldwrku )
2467 CALL
zlacpy(
'F', m, blk, work( iu ), ldwrku,
2482 CALL
zgelqf( m, n, a, lda, work( itau ),
2483 $ work( iwork ), lwork-iwork+1, ierr )
2487 CALL
zlacpy(
'L', m, m, a, lda, u, ldu )
2488 CALL
zlaset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2495 CALL
zunglq( m, n, m, a, lda, work( itau ),
2496 $ work( iwork ), lwork-iwork+1, ierr )
2506 CALL
zgebrd( m, m, u, ldu, s, rwork( ie ),
2507 $ work( itauq ), work( itaup ),
2508 $ work( iwork ), lwork-iwork+1, ierr )
2514 CALL
zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2515 $ work( itaup ), a, lda, work( iwork ),
2516 $ lwork-iwork+1, ierr )
2522 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2523 $ work( iwork ), lwork-iwork+1, ierr )
2532 CALL
zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a, lda,
2533 $ u, ldu, cdum, 1, rwork( irwork ), info )
2537 ELSE IF( wntvs )
THEN
2545 IF( lwork.GE.m*m+3*m )
THEN
2550 IF( lwork.GE.wrkbl+lda*m )
THEN
2561 itau = ir + ldwrkr*m
2568 CALL
zgelqf( m, n, a, lda, work( itau ),
2569 $ work( iwork ), lwork-iwork+1, ierr )
2573 CALL
zlacpy(
'L', m, m, a, lda, work( ir ),
2575 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2576 $ work( ir+ldwrkr ), ldwrkr )
2582 CALL
zunglq( m, n, m, a, lda, work( itau ),
2583 $ work( iwork ), lwork-iwork+1, ierr )
2593 CALL
zgebrd( m, m, work( ir ), ldwrkr, s,
2594 $ rwork( ie ), work( itauq ),
2595 $ work( itaup ), work( iwork ),
2596 $ lwork-iwork+1, ierr )
2603 CALL
zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2604 $ work( itaup ), work( iwork ),
2605 $ lwork-iwork+1, ierr )
2613 CALL
zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2614 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2615 $ rwork( irwork ), info )
2622 CALL
zgemm(
'N',
'N', m, n, m, cone, work( ir ),
2623 $ ldwrkr, a, lda, czero, vt, ldvt )
2636 CALL
zgelqf( m, n, a, lda, work( itau ),
2637 $ work( iwork ), lwork-iwork+1, ierr )
2641 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
2647 CALL
zunglq( m, n, m, vt, ldvt, work( itau ),
2648 $ work( iwork ), lwork-iwork+1, ierr )
2656 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2663 CALL
zgebrd( m, m, a, lda, s, rwork( ie ),
2664 $ work( itauq ), work( itaup ),
2665 $ work( iwork ), lwork-iwork+1, ierr )
2671 CALL
zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2672 $ work( itaup ), vt, ldvt,
2673 $ work( iwork ), lwork-iwork+1, ierr )
2681 CALL
zbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
2682 $ ldvt, cdum, 1, cdum, 1,
2683 $ rwork( irwork ), info )
2687 ELSE IF( wntuo )
THEN
2693 IF( lwork.GE.2*m*m+3*m )
THEN
2698 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2705 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2720 itau = ir + ldwrkr*m
2727 CALL
zgelqf( m, n, a, lda, work( itau ),
2728 $ work( iwork ), lwork-iwork+1, ierr )
2732 CALL
zlacpy(
'L', m, m, a, lda, work( iu ),
2734 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2735 $ work( iu+ldwrku ), ldwrku )
2741 CALL
zunglq( m, n, m, a, lda, work( itau ),
2742 $ work( iwork ), lwork-iwork+1, ierr )
2754 CALL
zgebrd( m, m, work( iu ), ldwrku, s,
2755 $ rwork( ie ), work( itauq ),
2756 $ work( itaup ), work( iwork ),
2757 $ lwork-iwork+1, ierr )
2758 CALL
zlacpy(
'L', m, m, work( iu ), ldwrku,
2759 $ work( ir ), ldwrkr )
2766 CALL
zungbr(
'P', m, m, m, work( iu ), ldwrku,
2767 $ work( itaup ), work( iwork ),
2768 $ lwork-iwork+1, ierr )
2774 CALL
zungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2775 $ work( itauq ), work( iwork ),
2776 $ lwork-iwork+1, ierr )
2785 CALL
zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2786 $ work( iu ), ldwrku, work( ir ),
2787 $ ldwrkr, cdum, 1, rwork( irwork ),
2795 CALL
zgemm(
'N',
'N', m, n, m, cone, work( iu ),
2796 $ ldwrku, a, lda, czero, vt, ldvt )
2802 CALL
zlacpy(
'F', m, m, work( ir ), ldwrkr, a,
2816 CALL
zgelqf( m, n, a, lda, work( itau ),
2817 $ work( iwork ), lwork-iwork+1, ierr )
2818 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
2824 CALL
zunglq( m, n, m, vt, ldvt, work( itau ),
2825 $ work( iwork ), lwork-iwork+1, ierr )
2833 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2840 CALL
zgebrd( m, m, a, lda, s, rwork( ie ),
2841 $ work( itauq ), work( itaup ),
2842 $ work( iwork ), lwork-iwork+1, ierr )
2848 CALL
zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2849 $ work( itaup ), vt, ldvt,
2850 $ work( iwork ), lwork-iwork+1, ierr )
2856 CALL
zungbr(
'Q', m, m, m, a, lda, work( itauq ),
2857 $ work( iwork ), lwork-iwork+1, ierr )
2866 CALL
zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
2867 $ ldvt, a, lda, cdum, 1,
2868 $ rwork( irwork ), info )
2872 ELSE IF( wntuas )
THEN
2879 IF( lwork.GE.m*m+3*m )
THEN
2884 IF( lwork.GE.wrkbl+lda*m )
THEN
2895 itau = iu + ldwrku*m
2902 CALL
zgelqf( m, n, a, lda, work( itau ),
2903 $ work( iwork ), lwork-iwork+1, ierr )
2907 CALL
zlacpy(
'L', m, m, a, lda, work( iu ),
2909 CALL
zlaset(
'U', m-1, m-1, czero, czero,
2910 $ work( iu+ldwrku ), ldwrku )
2916 CALL
zunglq( m, n, m, a, lda, work( itau ),
2917 $ work( iwork ), lwork-iwork+1, ierr )
2927 CALL
zgebrd( m, m, work( iu ), ldwrku, s,
2928 $ rwork( ie ), work( itauq ),
2929 $ work( itaup ), work( iwork ),
2930 $ lwork-iwork+1, ierr )
2931 CALL
zlacpy(
'L', m, m, work( iu ), ldwrku, u,
2939 CALL
zungbr(
'P', m, m, m, work( iu ), ldwrku,
2940 $ work( itaup ), work( iwork ),
2941 $ lwork-iwork+1, ierr )
2947 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2948 $ work( iwork ), lwork-iwork+1, ierr )
2957 CALL
zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2958 $ work( iu ), ldwrku, u, ldu, cdum, 1,
2959 $ rwork( irwork ), info )
2966 CALL
zgemm(
'N',
'N', m, n, m, cone, work( iu ),
2967 $ ldwrku, a, lda, czero, vt, ldvt )
2980 CALL
zgelqf( m, n, a, lda, work( itau ),
2981 $ work( iwork ), lwork-iwork+1, ierr )
2982 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
2988 CALL
zunglq( m, n, m, vt, ldvt, work( itau ),
2989 $ work( iwork ), lwork-iwork+1, ierr )
2993 CALL
zlacpy(
'L', m, m, a, lda, u, ldu )
2994 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3005 CALL
zgebrd( m, m, u, ldu, s, rwork( ie ),
3006 $ work( itauq ), work( itaup ),
3007 $ work( iwork ), lwork-iwork+1, ierr )
3014 CALL
zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3015 $ work( itaup ), vt, ldvt,
3016 $ work( iwork ), lwork-iwork+1, ierr )
3022 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3023 $ work( iwork ), lwork-iwork+1, ierr )
3032 CALL
zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3033 $ ldvt, u, ldu, cdum, 1,
3034 $ rwork( irwork ), info )
3040 ELSE IF( wntva )
THEN
3048 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3053 IF( lwork.GE.wrkbl+lda*m )
THEN
3064 itau = ir + ldwrkr*m
3071 CALL
zgelqf( m, n, a, lda, work( itau ),
3072 $ work( iwork ), lwork-iwork+1, ierr )
3073 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3077 CALL
zlacpy(
'L', m, m, a, lda, work( ir ),
3079 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3080 $ work( ir+ldwrkr ), ldwrkr )
3086 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3087 $ work( iwork ), lwork-iwork+1, ierr )
3097 CALL
zgebrd( m, m, work( ir ), ldwrkr, s,
3098 $ rwork( ie ), work( itauq ),
3099 $ work( itaup ), work( iwork ),
3100 $ lwork-iwork+1, ierr )
3107 CALL
zungbr(
'P', m, m, m, work( ir ), ldwrkr,
3108 $ work( itaup ), work( iwork ),
3109 $ lwork-iwork+1, ierr )
3117 CALL
zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3118 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3119 $ rwork( irwork ), info )
3126 CALL
zgemm(
'N',
'N', m, n, m, cone, work( ir ),
3127 $ ldwrkr, vt, ldvt, czero, a, lda )
3131 CALL
zlacpy(
'F', m, n, a, lda, vt, ldvt )
3144 CALL
zgelqf( m, n, a, lda, work( itau ),
3145 $ work( iwork ), lwork-iwork+1, ierr )
3146 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3152 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3153 $ work( iwork ), lwork-iwork+1, ierr )
3161 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3168 CALL
zgebrd( m, m, a, lda, s, rwork( ie ),
3169 $ work( itauq ), work( itaup ),
3170 $ work( iwork ), lwork-iwork+1, ierr )
3177 CALL
zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3178 $ work( itaup ), vt, ldvt,
3179 $ work( iwork ), lwork-iwork+1, ierr )
3187 CALL
zbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
3188 $ ldvt, cdum, 1, cdum, 1,
3189 $ rwork( irwork ), info )
3193 ELSE IF( wntuo )
THEN
3199 IF( lwork.GE.2*m*m+max( n+m, 3*m ) )
THEN
3204 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3211 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3226 itau = ir + ldwrkr*m
3233 CALL
zgelqf( m, n, a, lda, work( itau ),
3234 $ work( iwork ), lwork-iwork+1, ierr )
3235 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3241 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3242 $ work( iwork ), lwork-iwork+1, ierr )
3246 CALL
zlacpy(
'L', m, m, a, lda, work( iu ),
3248 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3249 $ work( iu+ldwrku ), ldwrku )
3261 CALL
zgebrd( m, m, work( iu ), ldwrku, s,
3262 $ rwork( ie ), work( itauq ),
3263 $ work( itaup ), work( iwork ),
3264 $ lwork-iwork+1, ierr )
3265 CALL
zlacpy(
'L', m, m, work( iu ), ldwrku,
3266 $ work( ir ), ldwrkr )
3273 CALL
zungbr(
'P', m, m, m, work( iu ), ldwrku,
3274 $ work( itaup ), work( iwork ),
3275 $ lwork-iwork+1, ierr )
3281 CALL
zungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3282 $ work( itauq ), work( iwork ),
3283 $ lwork-iwork+1, ierr )
3292 CALL
zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3293 $ work( iu ), ldwrku, work( ir ),
3294 $ ldwrkr, cdum, 1, rwork( irwork ),
3302 CALL
zgemm(
'N',
'N', m, n, m, cone, work( iu ),
3303 $ ldwrku, vt, ldvt, czero, a, lda )
3307 CALL
zlacpy(
'F', m, n, a, lda, vt, ldvt )
3311 CALL
zlacpy(
'F', m, m, work( ir ), ldwrkr, a,
3325 CALL
zgelqf( m, n, a, lda, work( itau ),
3326 $ work( iwork ), lwork-iwork+1, ierr )
3327 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3333 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3334 $ work( iwork ), lwork-iwork+1, ierr )
3342 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3349 CALL
zgebrd( m, m, a, lda, s, rwork( ie ),
3350 $ work( itauq ), work( itaup ),
3351 $ work( iwork ), lwork-iwork+1, ierr )
3358 CALL
zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3359 $ work( itaup ), vt, ldvt,
3360 $ work( iwork ), lwork-iwork+1, ierr )
3366 CALL
zungbr(
'Q', m, m, m, a, lda, work( itauq ),
3367 $ work( iwork ), lwork-iwork+1, ierr )
3376 CALL
zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3377 $ ldvt, a, lda, cdum, 1,
3378 $ rwork( irwork ), info )
3382 ELSE IF( wntuas )
THEN
3389 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3394 IF( lwork.GE.wrkbl+lda*m )
THEN
3405 itau = iu + ldwrku*m
3412 CALL
zgelqf( m, n, a, lda, work( itau ),
3413 $ work( iwork ), lwork-iwork+1, ierr )
3414 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3420 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3421 $ work( iwork ), lwork-iwork+1, ierr )
3425 CALL
zlacpy(
'L', m, m, a, lda, work( iu ),
3427 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3428 $ work( iu+ldwrku ), ldwrku )
3438 CALL
zgebrd( m, m, work( iu ), ldwrku, s,
3439 $ rwork( ie ), work( itauq ),
3440 $ work( itaup ), work( iwork ),
3441 $ lwork-iwork+1, ierr )
3442 CALL
zlacpy(
'L', m, m, work( iu ), ldwrku, u,
3449 CALL
zungbr(
'P', m, m, m, work( iu ), ldwrku,
3450 $ work( itaup ), work( iwork ),
3451 $ lwork-iwork+1, ierr )
3457 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3458 $ work( iwork ), lwork-iwork+1, ierr )
3467 CALL
zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3468 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3469 $ rwork( irwork ), info )
3476 CALL
zgemm(
'N',
'N', m, n, m, cone, work( iu ),
3477 $ ldwrku, vt, ldvt, czero, a, lda )
3481 CALL
zlacpy(
'F', m, n, a, lda, vt, ldvt )
3494 CALL
zgelqf( m, n, a, lda, work( itau ),
3495 $ work( iwork ), lwork-iwork+1, ierr )
3496 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3502 CALL
zunglq( n, n, m, vt, ldvt, work( itau ),
3503 $ work( iwork ), lwork-iwork+1, ierr )
3507 CALL
zlacpy(
'L', m, m, a, lda, u, ldu )
3508 CALL
zlaset(
'U', m-1, m-1, czero, czero,
3519 CALL
zgebrd( m, m, u, ldu, s, rwork( ie ),
3520 $ work( itauq ), work( itaup ),
3521 $ work( iwork ), lwork-iwork+1, ierr )
3528 CALL
zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3529 $ work( itaup ), vt, ldvt,
3530 $ work( iwork ), lwork-iwork+1, ierr )
3536 CALL
zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3537 $ work( iwork ), lwork-iwork+1, ierr )
3546 CALL
zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3547 $ ldvt, u, ldu, cdum, 1,
3548 $ rwork( irwork ), info )
3572 CALL
zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3573 $ work( itaup ), work( iwork ), lwork-iwork+1,
3582 CALL
zlacpy(
'L', m, m, a, lda, u, ldu )
3583 CALL
zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3584 $ work( iwork ), lwork-iwork+1, ierr )
3593 CALL
zlacpy(
'U', m, n, a, lda, vt, ldvt )
3598 CALL
zungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3599 $ work( iwork ), lwork-iwork+1, ierr )
3608 CALL
zungbr(
'Q', m, m, n, a, lda, work( itauq ),
3609 $ work( iwork ), lwork-iwork+1, ierr )
3618 CALL
zungbr(
'P', m, n, m, a, lda, work( itaup ),
3619 $ work( iwork ), lwork-iwork+1, ierr )
3622 IF( wntuas .OR. wntuo )
3626 IF( wntvas .OR. wntvo )
3630 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3638 CALL
zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3639 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3641 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3649 CALL
zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3650 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3660 CALL
zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3661 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3671 IF( iscl.EQ.1 )
THEN
3672 IF( anrm.GT.bignum )
3673 $ CALL
dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3675 IF( info.NE.0 .AND. anrm.GT.bignum )
3676 $ CALL
dlascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3677 $ rwork( ie ), minmn, ierr )
3678 IF( anrm.LT.smlnum )
3679 $ CALL
dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3681 IF( info.NE.0 .AND. anrm.LT.smlnum )
3682 $ CALL
dlascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3683 $ rwork( ie ), minmn, ierr )
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
LOGICAL function lsame(CA, CB)
LSAME
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
subroutine zgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD