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 )
265 DOUBLE PRECISION DLAMCH, ZLANGE
266 EXTERNAL lsame, ilaenv, dlamch, zlange
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 )
325 lwork_zgeqrf = int( cdum(1) )
327 CALL zungqr( m, n, n, a, lda, cdum(1), cdum(1), -1, ierr )
328 lwork_zungqr_n = int( cdum(1) )
329 CALL zungqr( m, m, n, a, lda, cdum(1), cdum(1), -1, ierr )
330 lwork_zungqr_m = int( cdum(1) )
332 CALL zgebrd( n, n, a, lda, s, dum(1), cdum(1),
333 $ cdum(1), cdum(1), -1, ierr )
334 lwork_zgebrd = int( cdum(1) )
336 CALL zungbr(
'P', n, n, n, a, lda, cdum(1),
337 $ cdum(1), -1, ierr )
338 lwork_zungbr_p = int( cdum(1) )
339 CALL zungbr(
'Q', n, n, n, a, lda, cdum(1),
340 $ cdum(1), -1, ierr )
341 lwork_zungbr_q = int( 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 )
448 lwork_zgebrd = int( cdum(1) )
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 = int( 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 = int( 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 )
474 lwork_zgelqf = int( cdum(1) )
476 CALL zunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
478 lwork_zunglq_n = int( cdum(1) )
479 CALL zunglq( m, n, m, a, lda, cdum(1), cdum(1), -1, ierr )
480 lwork_zunglq_m = int( cdum(1) )
482 CALL zgebrd( m, m, a, lda, s, dum(1), cdum(1),
483 $ cdum(1), cdum(1), -1, ierr )
484 lwork_zgebrd = int( cdum(1) )
486 CALL zungbr(
'P', m, m, m, a, n, cdum(1),
487 $ cdum(1), -1, ierr )
488 lwork_zungbr_p = int( cdum(1) )
490 CALL zungbr(
'Q', m, m, m, a, n, cdum(1),
491 $ cdum(1), -1, ierr )
492 lwork_zungbr_q = int( 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 )
598 lwork_zgebrd = int( cdum(1) )
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 = int( 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 = int( 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 )
684 CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
696 CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
697 $ work( itaup ), work( iwork ), lwork-iwork+1,
700 IF( wntvo .OR. wntvas )
THEN
706 CALL zungbr(
'P', n, n, n, a, lda, work( itaup ),
707 $ work( iwork ), lwork-iwork+1, ierr )
717 CALL zbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
718 $ cdum, 1, cdum, 1, rwork( irwork ), info )
723 $
CALL zlacpy(
'F', n, n, a, lda, vt, ldvt )
725 ELSE IF( wntuo .AND. wntvn )
THEN
731 IF( lwork.GE.n*n+3*n )
THEN
736 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
742 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
752 ldwrku = ( lwork-n*n ) / n
762 CALL zgeqrf( m, n, a, lda, work( itau ),
763 $ work( iwork ), lwork-iwork+1, ierr )
767 CALL zlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
768 CALL zlaset(
'L', n-1, n-1, czero, czero,
769 $ work( ir+1 ), ldwrkr )
775 CALL zungqr( m, n, n, a, lda, work( itau ),
776 $ work( iwork ), lwork-iwork+1, ierr )
786 CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
787 $ work( itauq ), work( itaup ),
788 $ work( iwork ), lwork-iwork+1, ierr )
794 CALL zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
795 $ work( itauq ), work( iwork ),
796 $ lwork-iwork+1, ierr )
804 CALL zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
805 $ work( ir ), ldwrkr, cdum, 1,
806 $ rwork( irwork ), info )
814 DO 10 i = 1, m, ldwrku
815 chunk = min( m-i+1, ldwrku )
816 CALL zgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
817 $ lda, work( ir ), ldwrkr, czero,
818 $ work( iu ), ldwrku )
819 CALL zlacpy(
'F', chunk, n, work( iu ), ldwrku,
836 CALL zgebrd( m, n, a, lda, s, rwork( ie ),
837 $ work( itauq ), work( itaup ),
838 $ work( iwork ), lwork-iwork+1, ierr )
844 CALL zungbr(
'Q', m, n, n, a, lda, work( itauq ),
845 $ work( iwork ), lwork-iwork+1, ierr )
853 CALL zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
854 $ a, lda, cdum, 1, rwork( irwork ), info )
858 ELSE IF( wntuo .AND. wntvas )
THEN
864 IF( lwork.GE.n*n+3*n )
THEN
869 IF( lwork.GE.max( wrkbl, lda*n )+lda*n )
THEN
875 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n )
THEN
885 ldwrku = ( lwork-n*n ) / n
895 CALL zgeqrf( m, n, a, lda, work( itau ),
896 $ work( iwork ), lwork-iwork+1, ierr )
900 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
902 $
CALL zlaset(
'L', n-1, n-1, czero, czero,
909 CALL zungqr( m, n, n, a, lda, work( itau ),
910 $ work( iwork ), lwork-iwork+1, ierr )
920 CALL zgebrd( n, n, vt, ldvt, s, rwork( ie ),
921 $ work( itauq ), work( itaup ),
922 $ work( iwork ), lwork-iwork+1, ierr )
923 CALL zlacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
929 CALL zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
930 $ work( itauq ), work( iwork ),
931 $ lwork-iwork+1, ierr )
937 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
938 $ work( iwork ), lwork-iwork+1, ierr )
947 CALL zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
948 $ ldvt, work( ir ), ldwrkr, cdum, 1,
949 $ rwork( irwork ), info )
957 DO 20 i = 1, m, ldwrku
958 chunk = min( m-i+1, ldwrku )
959 CALL zgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
960 $ lda, work( ir ), ldwrkr, czero,
961 $ work( iu ), ldwrku )
962 CALL zlacpy(
'F', chunk, n, work( iu ), ldwrku,
977 CALL zgeqrf( m, n, a, lda, work( itau ),
978 $ work( iwork ), lwork-iwork+1, ierr )
982 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
984 $
CALL zlaset(
'L', n-1, n-1, czero, czero,
991 CALL zungqr( m, n, n, a, lda, work( itau ),
992 $ work( iwork ), lwork-iwork+1, ierr )
1002 CALL zgebrd( n, n, vt, ldvt, s, rwork( ie ),
1003 $ work( itauq ), work( itaup ),
1004 $ work( iwork ), lwork-iwork+1, ierr )
1010 CALL zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1011 $ work( itauq ), a, lda, work( iwork ),
1012 $ lwork-iwork+1, ierr )
1018 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1019 $ work( iwork ), lwork-iwork+1, ierr )
1028 CALL zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1029 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1034 ELSE IF( wntus )
THEN
1042 IF( lwork.GE.n*n+3*n )
THEN
1047 IF( lwork.GE.wrkbl+lda*n )
THEN
1058 itau = ir + ldwrkr*n
1065 CALL zgeqrf( m, n, a, lda, work( itau ),
1066 $ work( iwork ), lwork-iwork+1, ierr )
1070 CALL zlacpy(
'U', n, n, a, lda, work( ir ),
1072 CALL zlaset(
'L', n-1, n-1, czero, czero,
1073 $ work( ir+1 ), ldwrkr )
1079 CALL zungqr( m, n, n, a, lda, work( itau ),
1080 $ work( iwork ), lwork-iwork+1, ierr )
1090 CALL zgebrd( n, n, work( ir ), ldwrkr, s,
1091 $ rwork( ie ), work( itauq ),
1092 $ work( itaup ), work( iwork ),
1093 $ lwork-iwork+1, ierr )
1099 CALL zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1100 $ work( itauq ), work( iwork ),
1101 $ lwork-iwork+1, ierr )
1109 CALL zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1110 $ 1, work( ir ), ldwrkr, cdum, 1,
1111 $ rwork( irwork ), info )
1118 CALL zgemm(
'N',
'N', m, n, n, cone, a, lda,
1119 $ work( ir ), ldwrkr, czero, u, ldu )
1132 CALL zgeqrf( m, n, a, lda, work( itau ),
1133 $ work( iwork ), lwork-iwork+1, ierr )
1134 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1140 CALL zungqr( m, n, n, u, ldu, work( itau ),
1141 $ work( iwork ), lwork-iwork+1, ierr )
1150 CALL zlaset(
'L', n-1, n-1, czero, czero,
1158 CALL zgebrd( n, n, a, lda, s, rwork( ie ),
1159 $ work( itauq ), work( itaup ),
1160 $ work( iwork ), lwork-iwork+1, ierr )
1166 CALL zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1167 $ work( itauq ), u, ldu, work( iwork ),
1168 $ lwork-iwork+1, ierr )
1176 CALL zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1177 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1182 ELSE IF( wntvo )
THEN
1188 IF( lwork.GE.2*n*n+3*n )
THEN
1193 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1200 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1215 itau = ir + ldwrkr*n
1222 CALL zgeqrf( m, n, a, lda, work( itau ),
1223 $ work( iwork ), lwork-iwork+1, ierr )
1227 CALL zlacpy(
'U', n, n, a, lda, work( iu ),
1229 CALL zlaset(
'L', n-1, n-1, czero, czero,
1230 $ work( iu+1 ), ldwrku )
1236 CALL zungqr( m, n, n, a, lda, work( itau ),
1237 $ work( iwork ), lwork-iwork+1, ierr )
1249 CALL zgebrd( n, n, work( iu ), ldwrku, s,
1250 $ rwork( ie ), work( itauq ),
1251 $ work( itaup ), work( iwork ),
1252 $ lwork-iwork+1, ierr )
1253 CALL zlacpy(
'U', n, n, work( iu ), ldwrku,
1254 $ work( ir ), ldwrkr )
1260 CALL zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1261 $ work( itauq ), work( iwork ),
1262 $ lwork-iwork+1, ierr )
1269 CALL zungbr(
'P', n, n, n, work( ir ), ldwrkr,
1270 $ work( itaup ), work( iwork ),
1271 $ lwork-iwork+1, ierr )
1280 CALL zbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1281 $ work( ir ), ldwrkr, work( iu ),
1282 $ ldwrku, cdum, 1, rwork( irwork ),
1290 CALL zgemm(
'N',
'N', m, n, n, cone, a, lda,
1291 $ work( iu ), ldwrku, czero, u, ldu )
1297 CALL zlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1311 CALL zgeqrf( m, n, a, lda, work( itau ),
1312 $ work( iwork ), lwork-iwork+1, ierr )
1313 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1319 CALL zungqr( m, n, n, u, ldu, work( itau ),
1320 $ work( iwork ), lwork-iwork+1, ierr )
1329 CALL zlaset(
'L', n-1, n-1, czero, czero,
1337 CALL zgebrd( n, n, a, lda, s, rwork( ie ),
1338 $ work( itauq ), work( itaup ),
1339 $ work( iwork ), lwork-iwork+1, ierr )
1345 CALL zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1346 $ work( itauq ), u, ldu, work( iwork ),
1347 $ lwork-iwork+1, ierr )
1353 CALL zungbr(
'P', n, n, n, a, lda, work( itaup ),
1354 $ work( iwork ), lwork-iwork+1, ierr )
1363 CALL zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1364 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1369 ELSE IF( wntvas )
THEN
1376 IF( lwork.GE.n*n+3*n )
THEN
1381 IF( lwork.GE.wrkbl+lda*n )
THEN
1392 itau = iu + ldwrku*n
1399 CALL zgeqrf( m, n, a, lda, work( itau ),
1400 $ work( iwork ), lwork-iwork+1, ierr )
1404 CALL zlacpy(
'U', n, n, a, lda, work( iu ),
1406 CALL zlaset(
'L', n-1, n-1, czero, czero,
1407 $ work( iu+1 ), ldwrku )
1413 CALL zungqr( m, n, n, a, lda, work( itau ),
1414 $ work( iwork ), lwork-iwork+1, ierr )
1424 CALL zgebrd( n, n, work( iu ), ldwrku, s,
1425 $ rwork( ie ), work( itauq ),
1426 $ work( itaup ), work( iwork ),
1427 $ lwork-iwork+1, ierr )
1428 CALL zlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1435 CALL zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1436 $ work( itauq ), work( iwork ),
1437 $ lwork-iwork+1, ierr )
1444 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1445 $ work( iwork ), lwork-iwork+1, ierr )
1454 CALL zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1455 $ ldvt, work( iu ), ldwrku, cdum, 1,
1456 $ rwork( irwork ), info )
1463 CALL zgemm(
'N',
'N', m, n, n, cone, a, lda,
1464 $ work( iu ), ldwrku, czero, u, ldu )
1477 CALL zgeqrf( m, n, a, lda, work( itau ),
1478 $ work( iwork ), lwork-iwork+1, ierr )
1479 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1485 CALL zungqr( m, n, n, u, ldu, work( itau ),
1486 $ work( iwork ), lwork-iwork+1, ierr )
1490 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
1492 $
CALL zlaset(
'L', n-1, n-1, czero, czero,
1493 $ vt( 2, 1 ), ldvt )
1503 CALL zgebrd( n, n, vt, ldvt, s, rwork( ie ),
1504 $ work( itauq ), work( itaup ),
1505 $ work( iwork ), lwork-iwork+1, ierr )
1512 CALL zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1513 $ work( itauq ), u, ldu, work( iwork ),
1514 $ lwork-iwork+1, ierr )
1520 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1521 $ work( iwork ), lwork-iwork+1, ierr )
1530 CALL zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1531 $ ldvt, u, ldu, cdum, 1,
1532 $ rwork( irwork ), info )
1538 ELSE IF( wntua )
THEN
1546 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1551 IF( lwork.GE.wrkbl+lda*n )
THEN
1562 itau = ir + ldwrkr*n
1569 CALL zgeqrf( m, n, a, lda, work( itau ),
1570 $ work( iwork ), lwork-iwork+1, ierr )
1571 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1575 CALL zlacpy(
'U', n, n, a, lda, work( ir ),
1577 CALL zlaset(
'L', n-1, n-1, czero, czero,
1578 $ work( ir+1 ), ldwrkr )
1584 CALL zungqr( m, m, n, u, ldu, work( itau ),
1585 $ work( iwork ), lwork-iwork+1, ierr )
1595 CALL zgebrd( n, n, work( ir ), ldwrkr, s,
1596 $ rwork( ie ), work( itauq ),
1597 $ work( itaup ), work( iwork ),
1598 $ lwork-iwork+1, ierr )
1604 CALL zungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1605 $ work( itauq ), work( iwork ),
1606 $ lwork-iwork+1, ierr )
1614 CALL zbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
1615 $ 1, work( ir ), ldwrkr, cdum, 1,
1616 $ rwork( irwork ), info )
1623 CALL zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1624 $ work( ir ), ldwrkr, czero, a, lda )
1628 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
1641 CALL zgeqrf( m, n, a, lda, work( itau ),
1642 $ work( iwork ), lwork-iwork+1, ierr )
1643 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1649 CALL zungqr( m, m, n, u, ldu, work( itau ),
1650 $ work( iwork ), lwork-iwork+1, ierr )
1659 CALL zlaset(
'L', n-1, n-1, czero, czero,
1667 CALL zgebrd( n, n, a, lda, s, rwork( ie ),
1668 $ work( itauq ), work( itaup ),
1669 $ work( iwork ), lwork-iwork+1, ierr )
1676 CALL zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1677 $ work( itauq ), u, ldu, work( iwork ),
1678 $ lwork-iwork+1, ierr )
1686 CALL zbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
1687 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1692 ELSE IF( wntvo )
THEN
1698 IF( lwork.GE.2*n*n+max( n+m, 3*n ) )
THEN
1703 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1710 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1725 itau = ir + ldwrkr*n
1732 CALL zgeqrf( m, n, a, lda, work( itau ),
1733 $ work( iwork ), lwork-iwork+1, ierr )
1734 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1740 CALL zungqr( m, m, n, u, ldu, work( itau ),
1741 $ work( iwork ), lwork-iwork+1, ierr )
1745 CALL zlacpy(
'U', n, n, a, lda, work( iu ),
1747 CALL zlaset(
'L', n-1, n-1, czero, czero,
1748 $ work( iu+1 ), ldwrku )
1760 CALL zgebrd( n, n, work( iu ), ldwrku, s,
1761 $ rwork( ie ), work( itauq ),
1762 $ work( itaup ), work( iwork ),
1763 $ lwork-iwork+1, ierr )
1764 CALL zlacpy(
'U', n, n, work( iu ), ldwrku,
1765 $ work( ir ), ldwrkr )
1771 CALL zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1772 $ work( itauq ), work( iwork ),
1773 $ lwork-iwork+1, ierr )
1780 CALL zungbr(
'P', n, n, n, work( ir ), ldwrkr,
1781 $ work( itaup ), work( iwork ),
1782 $ lwork-iwork+1, ierr )
1791 CALL zbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1792 $ work( ir ), ldwrkr, work( iu ),
1793 $ ldwrku, cdum, 1, rwork( irwork ),
1801 CALL zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1802 $ work( iu ), ldwrku, czero, a, lda )
1806 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
1810 CALL zlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1824 CALL zgeqrf( m, n, a, lda, work( itau ),
1825 $ work( iwork ), lwork-iwork+1, ierr )
1826 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1832 CALL zungqr( m, m, n, u, ldu, work( itau ),
1833 $ work( iwork ), lwork-iwork+1, ierr )
1842 CALL zlaset(
'L', n-1, n-1, czero, czero,
1850 CALL zgebrd( n, n, a, lda, s, rwork( ie ),
1851 $ work( itauq ), work( itaup ),
1852 $ work( iwork ), lwork-iwork+1, ierr )
1859 CALL zunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1860 $ work( itauq ), u, ldu, work( iwork ),
1861 $ lwork-iwork+1, ierr )
1867 CALL zungbr(
'P', n, n, n, a, lda, work( itaup ),
1868 $ work( iwork ), lwork-iwork+1, ierr )
1877 CALL zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1878 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1883 ELSE IF( wntvas )
THEN
1890 IF( lwork.GE.n*n+max( n+m, 3*n ) )
THEN
1895 IF( lwork.GE.wrkbl+lda*n )
THEN
1906 itau = iu + ldwrku*n
1913 CALL zgeqrf( m, n, a, lda, work( itau ),
1914 $ work( iwork ), lwork-iwork+1, ierr )
1915 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1921 CALL zungqr( m, m, n, u, ldu, work( itau ),
1922 $ work( iwork ), lwork-iwork+1, ierr )
1926 CALL zlacpy(
'U', n, n, a, lda, work( iu ),
1928 CALL zlaset(
'L', n-1, n-1, czero, czero,
1929 $ work( iu+1 ), ldwrku )
1939 CALL zgebrd( n, n, work( iu ), ldwrku, s,
1940 $ rwork( ie ), work( itauq ),
1941 $ work( itaup ), work( iwork ),
1942 $ lwork-iwork+1, ierr )
1943 CALL zlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1950 CALL zungbr(
'Q', n, n, n, work( iu ), ldwrku,
1951 $ work( itauq ), work( iwork ),
1952 $ lwork-iwork+1, ierr )
1959 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1960 $ work( iwork ), lwork-iwork+1, ierr )
1969 CALL zbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
1970 $ ldvt, work( iu ), ldwrku, cdum, 1,
1971 $ rwork( irwork ), info )
1978 CALL zgemm(
'N',
'N', m, n, n, cone, u, ldu,
1979 $ work( iu ), ldwrku, czero, a, lda )
1983 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
1996 CALL zgeqrf( m, n, a, lda, work( itau ),
1997 $ work( iwork ), lwork-iwork+1, ierr )
1998 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
2004 CALL zungqr( m, m, n, u, ldu, work( itau ),
2005 $ work( iwork ), lwork-iwork+1, ierr )
2009 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
2011 $
CALL zlaset(
'L', n-1, n-1, czero, czero,
2012 $ vt( 2, 1 ), ldvt )
2022 CALL zgebrd( n, n, vt, ldvt, s, rwork( ie ),
2023 $ work( itauq ), work( itaup ),
2024 $ work( iwork ), lwork-iwork+1, ierr )
2031 CALL zunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2032 $ work( itauq ), u, ldu, work( iwork ),
2033 $ lwork-iwork+1, ierr )
2039 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2040 $ work( iwork ), lwork-iwork+1, ierr )
2049 CALL zbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
2050 $ ldvt, u, ldu, cdum, 1,
2051 $ rwork( irwork ), info )
2075 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2076 $ work( itaup ), work( iwork ), lwork-iwork+1,
2085 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
2090 CALL zungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2091 $ work( iwork ), lwork-iwork+1, ierr )
2100 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
2101 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2102 $ work( iwork ), lwork-iwork+1, ierr )
2111 CALL zungbr(
'Q', m, n, n, a, lda, work( itauq ),
2112 $ work( iwork ), lwork-iwork+1, ierr )
2121 CALL zungbr(
'P', n, n, n, a, lda, work( itaup ),
2122 $ work( iwork ), lwork-iwork+1, ierr )
2125 IF( wntuas .OR. wntuo )
2129 IF( wntvas .OR. wntvo )
2133 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2141 CALL zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2142 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2144 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2152 CALL zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2153 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2163 CALL zbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2164 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2176 IF( n.GE.mnthr )
THEN
2190 CALL zgelqf( m, n, a, lda, work( itau ), work( iwork ),
2191 $ lwork-iwork+1, ierr )
2195 CALL zlaset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2206 CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2207 $ work( itaup ), work( iwork ), lwork-iwork+1,
2209 IF( wntuo .OR. wntuas )
THEN
2215 CALL zungbr(
'Q', m, m, m, a, lda, work( itauq ),
2216 $ work( iwork ), lwork-iwork+1, ierr )
2220 IF( wntuo .OR. wntuas )
2228 CALL zbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2229 $ a, lda, cdum, 1, rwork( irwork ), info )
2234 $
CALL zlacpy(
'F', m, m, a, lda, u, ldu )
2236 ELSE IF( wntvo .AND. wntun )
THEN
2242 IF( lwork.GE.m*m+3*m )
THEN
2247 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2254 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2266 chunk = ( lwork-m*m ) / m
2269 itau = ir + ldwrkr*m
2276 CALL zgelqf( m, n, a, lda, work( itau ),
2277 $ work( iwork ), lwork-iwork+1, ierr )
2281 CALL zlacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2282 CALL zlaset(
'U', m-1, m-1, czero, czero,
2283 $ work( ir+ldwrkr ), ldwrkr )
2289 CALL zunglq( m, n, m, a, lda, work( itau ),
2290 $ work( iwork ), lwork-iwork+1, ierr )
2300 CALL zgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2301 $ work( itauq ), work( itaup ),
2302 $ work( iwork ), lwork-iwork+1, ierr )
2308 CALL zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2309 $ work( itaup ), work( iwork ),
2310 $ lwork-iwork+1, ierr )
2318 CALL zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2319 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2320 $ rwork( irwork ), info )
2328 DO 30 i = 1, n, chunk
2329 blk = min( n-i+1, chunk )
2330 CALL zgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2331 $ ldwrkr, a( 1, i ), lda, czero,
2332 $ work( iu ), ldwrku )
2333 CALL zlacpy(
'F', m, blk, work( iu ), ldwrku,
2350 CALL zgebrd( m, n, a, lda, s, rwork( ie ),
2351 $ work( itauq ), work( itaup ),
2352 $ work( iwork ), lwork-iwork+1, ierr )
2358 CALL zungbr(
'P', m, n, m, a, lda, work( itaup ),
2359 $ work( iwork ), lwork-iwork+1, ierr )
2367 CALL zbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2368 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2372 ELSE IF( wntvo .AND. wntuas )
THEN
2378 IF( lwork.GE.m*m+3*m )
THEN
2383 IF( lwork.GE.max( wrkbl, lda*n )+lda*m )
THEN
2390 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m )
THEN
2402 chunk = ( lwork-m*m ) / m
2405 itau = ir + ldwrkr*m
2412 CALL zgelqf( m, n, a, lda, work( itau ),
2413 $ work( iwork ), lwork-iwork+1, ierr )
2417 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
2418 CALL zlaset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2425 CALL zunglq( m, n, m, a, lda, work( itau ),
2426 $ work( iwork ), lwork-iwork+1, ierr )
2436 CALL zgebrd( m, m, u, ldu, s, rwork( ie ),
2437 $ work( itauq ), work( itaup ),
2438 $ work( iwork ), lwork-iwork+1, ierr )
2439 CALL zlacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2445 CALL zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2446 $ work( itaup ), work( iwork ),
2447 $ lwork-iwork+1, ierr )
2453 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2454 $ work( iwork ), lwork-iwork+1, ierr )
2463 CALL zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2464 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2465 $ rwork( irwork ), info )
2473 DO 40 i = 1, n, chunk
2474 blk = min( n-i+1, chunk )
2475 CALL zgemm(
'N',
'N', m, blk, m, cone, work( ir ),
2476 $ ldwrkr, a( 1, i ), lda, czero,
2477 $ work( iu ), ldwrku )
2478 CALL zlacpy(
'F', m, blk, work( iu ), ldwrku,
2493 CALL zgelqf( m, n, a, lda, work( itau ),
2494 $ work( iwork ), lwork-iwork+1, ierr )
2498 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
2499 CALL zlaset(
'U', m-1, m-1, czero, czero, u( 1, 2 ),
2506 CALL zunglq( m, n, m, a, lda, work( itau ),
2507 $ work( iwork ), lwork-iwork+1, ierr )
2517 CALL zgebrd( m, m, u, ldu, s, rwork( ie ),
2518 $ work( itauq ), work( itaup ),
2519 $ work( iwork ), lwork-iwork+1, ierr )
2525 CALL zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2526 $ work( itaup ), a, lda, work( iwork ),
2527 $ lwork-iwork+1, ierr )
2533 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2534 $ work( iwork ), lwork-iwork+1, ierr )
2543 CALL zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a, lda,
2544 $ u, ldu, cdum, 1, rwork( irwork ), info )
2548 ELSE IF( wntvs )
THEN
2556 IF( lwork.GE.m*m+3*m )
THEN
2561 IF( lwork.GE.wrkbl+lda*m )
THEN
2572 itau = ir + ldwrkr*m
2579 CALL zgelqf( m, n, a, lda, work( itau ),
2580 $ work( iwork ), lwork-iwork+1, ierr )
2584 CALL zlacpy(
'L', m, m, a, lda, work( ir ),
2586 CALL zlaset(
'U', m-1, m-1, czero, czero,
2587 $ work( ir+ldwrkr ), ldwrkr )
2593 CALL zunglq( m, n, m, a, lda, work( itau ),
2594 $ work( iwork ), lwork-iwork+1, ierr )
2604 CALL zgebrd( m, m, work( ir ), ldwrkr, s,
2605 $ rwork( ie ), work( itauq ),
2606 $ work( itaup ), work( iwork ),
2607 $ lwork-iwork+1, ierr )
2614 CALL zungbr(
'P', m, m, m, work( ir ), ldwrkr,
2615 $ work( itaup ), work( iwork ),
2616 $ lwork-iwork+1, ierr )
2624 CALL zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2625 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2626 $ rwork( irwork ), info )
2633 CALL zgemm(
'N',
'N', m, n, m, cone, work( ir ),
2634 $ ldwrkr, a, lda, czero, vt, ldvt )
2647 CALL zgelqf( m, n, a, lda, work( itau ),
2648 $ work( iwork ), lwork-iwork+1, ierr )
2652 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
2658 CALL zunglq( m, n, m, vt, ldvt, work( itau ),
2659 $ work( iwork ), lwork-iwork+1, ierr )
2667 CALL zlaset(
'U', m-1, m-1, czero, czero,
2674 CALL zgebrd( m, m, a, lda, s, rwork( ie ),
2675 $ work( itauq ), work( itaup ),
2676 $ work( iwork ), lwork-iwork+1, ierr )
2682 CALL zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2683 $ work( itaup ), vt, ldvt,
2684 $ work( iwork ), lwork-iwork+1, ierr )
2692 CALL zbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
2693 $ ldvt, cdum, 1, cdum, 1,
2694 $ rwork( irwork ), info )
2698 ELSE IF( wntuo )
THEN
2704 IF( lwork.GE.2*m*m+3*m )
THEN
2709 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2716 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2731 itau = ir + ldwrkr*m
2738 CALL zgelqf( m, n, a, lda, work( itau ),
2739 $ work( iwork ), lwork-iwork+1, ierr )
2743 CALL zlacpy(
'L', m, m, a, lda, work( iu ),
2745 CALL zlaset(
'U', m-1, m-1, czero, czero,
2746 $ work( iu+ldwrku ), ldwrku )
2752 CALL zunglq( m, n, m, a, lda, work( itau ),
2753 $ work( iwork ), lwork-iwork+1, ierr )
2765 CALL zgebrd( m, m, work( iu ), ldwrku, s,
2766 $ rwork( ie ), work( itauq ),
2767 $ work( itaup ), work( iwork ),
2768 $ lwork-iwork+1, ierr )
2769 CALL zlacpy(
'L', m, m, work( iu ), ldwrku,
2770 $ work( ir ), ldwrkr )
2777 CALL zungbr(
'P', m, m, m, work( iu ), ldwrku,
2778 $ work( itaup ), work( iwork ),
2779 $ lwork-iwork+1, ierr )
2785 CALL zungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2786 $ work( itauq ), work( iwork ),
2787 $ lwork-iwork+1, ierr )
2796 CALL zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2797 $ work( iu ), ldwrku, work( ir ),
2798 $ ldwrkr, cdum, 1, rwork( irwork ),
2806 CALL zgemm(
'N',
'N', m, n, m, cone, work( iu ),
2807 $ ldwrku, a, lda, czero, vt, ldvt )
2813 CALL zlacpy(
'F', m, m, work( ir ), ldwrkr, a,
2827 CALL zgelqf( m, n, a, lda, work( itau ),
2828 $ work( iwork ), lwork-iwork+1, ierr )
2829 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
2835 CALL zunglq( m, n, m, vt, ldvt, work( itau ),
2836 $ work( iwork ), lwork-iwork+1, ierr )
2844 CALL zlaset(
'U', m-1, m-1, czero, czero,
2851 CALL zgebrd( m, m, a, lda, s, rwork( ie ),
2852 $ work( itauq ), work( itaup ),
2853 $ work( iwork ), lwork-iwork+1, ierr )
2859 CALL zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2860 $ work( itaup ), vt, ldvt,
2861 $ work( iwork ), lwork-iwork+1, ierr )
2867 CALL zungbr(
'Q', m, m, m, a, lda, work( itauq ),
2868 $ work( iwork ), lwork-iwork+1, ierr )
2877 CALL zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
2878 $ ldvt, a, lda, cdum, 1,
2879 $ rwork( irwork ), info )
2883 ELSE IF( wntuas )
THEN
2890 IF( lwork.GE.m*m+3*m )
THEN
2895 IF( lwork.GE.wrkbl+lda*m )
THEN
2906 itau = iu + ldwrku*m
2913 CALL zgelqf( m, n, a, lda, work( itau ),
2914 $ work( iwork ), lwork-iwork+1, ierr )
2918 CALL zlacpy(
'L', m, m, a, lda, work( iu ),
2920 CALL zlaset(
'U', m-1, m-1, czero, czero,
2921 $ work( iu+ldwrku ), ldwrku )
2927 CALL zunglq( m, n, m, a, lda, work( itau ),
2928 $ work( iwork ), lwork-iwork+1, ierr )
2938 CALL zgebrd( m, m, work( iu ), ldwrku, s,
2939 $ rwork( ie ), work( itauq ),
2940 $ work( itaup ), work( iwork ),
2941 $ lwork-iwork+1, ierr )
2942 CALL zlacpy(
'L', m, m, work( iu ), ldwrku, u,
2950 CALL zungbr(
'P', m, m, m, work( iu ), ldwrku,
2951 $ work( itaup ), work( iwork ),
2952 $ lwork-iwork+1, ierr )
2958 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2959 $ work( iwork ), lwork-iwork+1, ierr )
2968 CALL zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2969 $ work( iu ), ldwrku, u, ldu, cdum, 1,
2970 $ rwork( irwork ), info )
2977 CALL zgemm(
'N',
'N', m, n, m, cone, work( iu ),
2978 $ ldwrku, a, lda, czero, vt, ldvt )
2991 CALL zgelqf( m, n, a, lda, work( itau ),
2992 $ work( iwork ), lwork-iwork+1, ierr )
2993 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
2999 CALL zunglq( m, n, m, vt, ldvt, work( itau ),
3000 $ work( iwork ), lwork-iwork+1, ierr )
3004 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
3005 CALL zlaset(
'U', m-1, m-1, czero, czero,
3016 CALL zgebrd( m, m, u, ldu, s, rwork( ie ),
3017 $ work( itauq ), work( itaup ),
3018 $ work( iwork ), lwork-iwork+1, ierr )
3025 CALL zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3026 $ work( itaup ), vt, ldvt,
3027 $ work( iwork ), lwork-iwork+1, ierr )
3033 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3034 $ work( iwork ), lwork-iwork+1, ierr )
3043 CALL zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3044 $ ldvt, u, ldu, cdum, 1,
3045 $ rwork( irwork ), info )
3051 ELSE IF( wntva )
THEN
3059 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3064 IF( lwork.GE.wrkbl+lda*m )
THEN
3075 itau = ir + ldwrkr*m
3082 CALL zgelqf( m, n, a, lda, work( itau ),
3083 $ work( iwork ), lwork-iwork+1, ierr )
3084 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3088 CALL zlacpy(
'L', m, m, a, lda, work( ir ),
3090 CALL zlaset(
'U', m-1, m-1, czero, czero,
3091 $ work( ir+ldwrkr ), ldwrkr )
3097 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3098 $ work( iwork ), lwork-iwork+1, ierr )
3108 CALL zgebrd( m, m, work( ir ), ldwrkr, s,
3109 $ rwork( ie ), work( itauq ),
3110 $ work( itaup ), work( iwork ),
3111 $ lwork-iwork+1, ierr )
3118 CALL zungbr(
'P', m, m, m, work( ir ), ldwrkr,
3119 $ work( itaup ), work( iwork ),
3120 $ lwork-iwork+1, ierr )
3128 CALL zbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3129 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3130 $ rwork( irwork ), info )
3137 CALL zgemm(
'N',
'N', m, n, m, cone, work( ir ),
3138 $ ldwrkr, vt, ldvt, czero, a, lda )
3142 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
3155 CALL zgelqf( m, n, a, lda, work( itau ),
3156 $ work( iwork ), lwork-iwork+1, ierr )
3157 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3163 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3164 $ work( iwork ), lwork-iwork+1, ierr )
3172 CALL zlaset(
'U', m-1, m-1, czero, czero,
3179 CALL zgebrd( m, m, a, lda, s, rwork( ie ),
3180 $ work( itauq ), work( itaup ),
3181 $ work( iwork ), lwork-iwork+1, ierr )
3188 CALL zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3189 $ work( itaup ), vt, ldvt,
3190 $ work( iwork ), lwork-iwork+1, ierr )
3198 CALL zbdsqr(
'U', m, n, 0, 0, s, rwork( ie ), vt,
3199 $ ldvt, cdum, 1, cdum, 1,
3200 $ rwork( irwork ), info )
3204 ELSE IF( wntuo )
THEN
3210 IF( lwork.GE.2*m*m+max( n+m, 3*m ) )
THEN
3215 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3222 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3237 itau = ir + ldwrkr*m
3244 CALL zgelqf( m, n, a, lda, work( itau ),
3245 $ work( iwork ), lwork-iwork+1, ierr )
3246 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3252 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3253 $ work( iwork ), lwork-iwork+1, ierr )
3257 CALL zlacpy(
'L', m, m, a, lda, work( iu ),
3259 CALL zlaset(
'U', m-1, m-1, czero, czero,
3260 $ work( iu+ldwrku ), ldwrku )
3272 CALL zgebrd( m, m, work( iu ), ldwrku, s,
3273 $ rwork( ie ), work( itauq ),
3274 $ work( itaup ), work( iwork ),
3275 $ lwork-iwork+1, ierr )
3276 CALL zlacpy(
'L', m, m, work( iu ), ldwrku,
3277 $ work( ir ), ldwrkr )
3284 CALL zungbr(
'P', m, m, m, work( iu ), ldwrku,
3285 $ work( itaup ), work( iwork ),
3286 $ lwork-iwork+1, ierr )
3292 CALL zungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3293 $ work( itauq ), work( iwork ),
3294 $ lwork-iwork+1, ierr )
3303 CALL zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3304 $ work( iu ), ldwrku, work( ir ),
3305 $ ldwrkr, cdum, 1, rwork( irwork ),
3313 CALL zgemm(
'N',
'N', m, n, m, cone, work( iu ),
3314 $ ldwrku, vt, ldvt, czero, a, lda )
3318 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
3322 CALL zlacpy(
'F', m, m, work( ir ), ldwrkr, a,
3336 CALL zgelqf( m, n, a, lda, work( itau ),
3337 $ work( iwork ), lwork-iwork+1, ierr )
3338 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3344 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3345 $ work( iwork ), lwork-iwork+1, ierr )
3353 CALL zlaset(
'U', m-1, m-1, czero, czero,
3360 CALL zgebrd( m, m, a, lda, s, rwork( ie ),
3361 $ work( itauq ), work( itaup ),
3362 $ work( iwork ), lwork-iwork+1, ierr )
3369 CALL zunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3370 $ work( itaup ), vt, ldvt,
3371 $ work( iwork ), lwork-iwork+1, ierr )
3377 CALL zungbr(
'Q', m, m, m, a, lda, work( itauq ),
3378 $ work( iwork ), lwork-iwork+1, ierr )
3387 CALL zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3388 $ ldvt, a, lda, cdum, 1,
3389 $ rwork( irwork ), info )
3393 ELSE IF( wntuas )
THEN
3400 IF( lwork.GE.m*m+max( n+m, 3*m ) )
THEN
3405 IF( lwork.GE.wrkbl+lda*m )
THEN
3416 itau = iu + ldwrku*m
3423 CALL zgelqf( m, n, a, lda, work( itau ),
3424 $ work( iwork ), lwork-iwork+1, ierr )
3425 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3431 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3432 $ work( iwork ), lwork-iwork+1, ierr )
3436 CALL zlacpy(
'L', m, m, a, lda, work( iu ),
3438 CALL zlaset(
'U', m-1, m-1, czero, czero,
3439 $ work( iu+ldwrku ), ldwrku )
3449 CALL zgebrd( m, m, work( iu ), ldwrku, s,
3450 $ rwork( ie ), work( itauq ),
3451 $ work( itaup ), work( iwork ),
3452 $ lwork-iwork+1, ierr )
3453 CALL zlacpy(
'L', m, m, work( iu ), ldwrku, u,
3460 CALL zungbr(
'P', m, m, m, work( iu ), ldwrku,
3461 $ work( itaup ), work( iwork ),
3462 $ lwork-iwork+1, ierr )
3468 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3469 $ work( iwork ), lwork-iwork+1, ierr )
3478 CALL zbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3479 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3480 $ rwork( irwork ), info )
3487 CALL zgemm(
'N',
'N', m, n, m, cone, work( iu ),
3488 $ ldwrku, vt, ldvt, czero, a, lda )
3492 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
3505 CALL zgelqf( m, n, a, lda, work( itau ),
3506 $ work( iwork ), lwork-iwork+1, ierr )
3507 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3513 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
3514 $ work( iwork ), lwork-iwork+1, ierr )
3518 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
3519 CALL zlaset(
'U', m-1, m-1, czero, czero,
3530 CALL zgebrd( m, m, u, ldu, s, rwork( ie ),
3531 $ work( itauq ), work( itaup ),
3532 $ work( iwork ), lwork-iwork+1, ierr )
3539 CALL zunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3540 $ work( itaup ), vt, ldvt,
3541 $ work( iwork ), lwork-iwork+1, ierr )
3547 CALL zungbr(
'Q', m, m, m, u, ldu, work( itauq ),
3548 $ work( iwork ), lwork-iwork+1, ierr )
3557 CALL zbdsqr(
'U', m, n, m, 0, s, rwork( ie ), vt,
3558 $ ldvt, u, ldu, cdum, 1,
3559 $ rwork( irwork ), info )
3583 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3584 $ work( itaup ), work( iwork ), lwork-iwork+1,
3593 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
3594 CALL zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3595 $ work( iwork ), lwork-iwork+1, ierr )
3604 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
3609 CALL zungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3610 $ work( iwork ), lwork-iwork+1, ierr )
3619 CALL zungbr(
'Q', m, m, n, a, lda, work( itauq ),
3620 $ work( iwork ), lwork-iwork+1, ierr )
3629 CALL zungbr(
'P', m, n, m, a, lda, work( itaup ),
3630 $ work( iwork ), lwork-iwork+1, ierr )
3633 IF( wntuas .OR. wntuo )
3637 IF( wntvas .OR. wntvo )
3641 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3649 CALL zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3650 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3652 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3660 CALL zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3661 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3671 CALL zbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3672 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3682 IF( iscl.EQ.1 )
THEN
3683 IF( anrm.GT.bignum )
3684 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3686 IF( info.NE.0 .AND. anrm.GT.bignum )
3687 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3688 $ rwork( ie ), minmn, ierr )
3689 IF( anrm.LT.smlnum )
3690 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3692 IF( info.NE.0 .AND. anrm.LT.smlnum )
3693 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3694 $ rwork( ie ), minmn, ierr )
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
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
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
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 zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
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 xerbla(SRNAME, INFO)
XERBLA
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
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.