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