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.
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