211 SUBROUTINE dgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
212 $ vt, ldvt, work, lwork, info )
220 CHARACTER JOBU, JOBVT
221 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
224 DOUBLE PRECISION A( lda, * ), S( * ), U( ldu, * ),
225 $ vt( ldvt, * ), work( * )
231 DOUBLE PRECISION ZERO, ONE
232 parameter( zero = 0.0d0, one = 1.0d0 )
235 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
236 $ wntva, wntvas, wntvn, wntvo, wntvs
237 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IR, ISCL,
238 $ itau, itaup, itauq, iu, iwork, ldwrkr, ldwrku,
239 $ maxwrk, minmn, minwrk, mnthr, ncu, ncvt, nru,
241 INTEGER LWORK_DGEQRF, LWORK_DORGQR_N, LWORK_DORGQR_M,
242 $ lwork_dgebrd, lwork_dorgbr_p, lwork_dorgbr_q,
243 $ lwork_dgelqf, lwork_dorglq_n, lwork_dorglq_m
244 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
247 DOUBLE PRECISION DUM( 1 )
257 DOUBLE PRECISION DLAMCH, DLANGE
258 EXTERNAL lsame, ilaenv, dlamch, dlange
261 INTRINSIC max, min, sqrt
269 wntua = lsame( jobu,
'A' )
270 wntus = lsame( jobu,
'S' )
271 wntuas = wntua .OR. wntus
272 wntuo = lsame( jobu,
'O' )
273 wntun = lsame( jobu,
'N' )
274 wntva = lsame( jobvt,
'A' )
275 wntvs = lsame( jobvt,
'S' )
276 wntvas = wntva .OR. wntvs
277 wntvo = lsame( jobvt,
'O' )
278 wntvn = lsame( jobvt,
'N' )
279 lquery = ( lwork.EQ.-1 )
281 IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) )
THEN
283 ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
284 $ ( wntvo .AND. wntuo ) )
THEN
286 ELSE IF( m.LT.0 )
THEN
288 ELSE IF( n.LT.0 )
THEN
290 ELSE IF( lda.LT.max( 1, m ) )
THEN
292 ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) )
THEN
294 ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
295 $ ( wntvs .AND. ldvt.LT.minmn ) )
THEN
309 IF( m.GE.n .AND. minmn.GT.0 )
THEN
313 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0, 0 )
316 CALL dgeqrf( m, n, a, lda, dum(1), dum(1), -1, ierr )
319 CALL dorgqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
320 lwork_dorgqr_n=dum(1)
321 CALL dorgqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
322 lwork_dorgqr_m=dum(1)
324 CALL dgebrd( n, n, a, lda, s, dum(1), dum(1),
325 $ dum(1), dum(1), -1, ierr )
328 CALL dorgbr(
'P', n, n, n, a, lda, dum(1),
330 lwork_dorgbr_p=dum(1)
332 CALL dorgbr(
'Q', n, n, n, a, lda, dum(1),
334 lwork_dorgbr_q=dum(1)
336 IF( m.GE.mnthr )
THEN
341 maxwrk = n + lwork_dgeqrf
342 maxwrk = max( maxwrk, 3*n+lwork_dgebrd )
343 IF( wntvo .OR. wntvas )
344 $ maxwrk = max( maxwrk, 3*n+lwork_dorgbr_p )
345 maxwrk = max( maxwrk, bdspac )
346 minwrk = max( 4*n, bdspac )
347 ELSE IF( wntuo .AND. wntvn )
THEN
351 wrkbl = n + lwork_dgeqrf
352 wrkbl = max( wrkbl, n+lwork_dorgqr_n )
353 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
354 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
355 wrkbl = max( wrkbl, bdspac )
356 maxwrk = max( n*n+wrkbl, n*n+m*n+n )
357 minwrk = max( 3*n+m, bdspac )
358 ELSE IF( wntuo .AND. wntvas )
THEN
363 wrkbl = n + lwork_dgeqrf
364 wrkbl = max( wrkbl, n+lwork_dorgqr_n )
365 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
366 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
367 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
368 wrkbl = max( wrkbl, bdspac )
369 maxwrk = max( n*n+wrkbl, n*n+m*n+n )
370 minwrk = max( 3*n+m, bdspac )
371 ELSE IF( wntus .AND. wntvn )
THEN
375 wrkbl = n + lwork_dgeqrf
376 wrkbl = max( wrkbl, n+lwork_dorgqr_n )
377 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
378 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
379 wrkbl = max( wrkbl, bdspac )
381 minwrk = max( 3*n+m, bdspac )
382 ELSE IF( wntus .AND. wntvo )
THEN
386 wrkbl = n + lwork_dgeqrf
387 wrkbl = max( wrkbl, n+lwork_dorgqr_n )
388 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
389 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
390 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
391 wrkbl = max( wrkbl, bdspac )
392 maxwrk = 2*n*n + wrkbl
393 minwrk = max( 3*n+m, bdspac )
394 ELSE IF( wntus .AND. wntvas )
THEN
399 wrkbl = n + lwork_dgeqrf
400 wrkbl = max( wrkbl, n+lwork_dorgqr_n )
401 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
402 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
403 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
404 wrkbl = max( wrkbl, bdspac )
406 minwrk = max( 3*n+m, bdspac )
407 ELSE IF( wntua .AND. wntvn )
THEN
411 wrkbl = n + lwork_dgeqrf
412 wrkbl = max( wrkbl, n+lwork_dorgqr_m )
413 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
414 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
415 wrkbl = max( wrkbl, bdspac )
417 minwrk = max( 3*n+m, bdspac )
418 ELSE IF( wntua .AND. wntvo )
THEN
422 wrkbl = n + lwork_dgeqrf
423 wrkbl = max( wrkbl, n+lwork_dorgqr_m )
424 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
425 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
426 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
427 wrkbl = max( wrkbl, bdspac )
428 maxwrk = 2*n*n + wrkbl
429 minwrk = max( 3*n+m, bdspac )
430 ELSE IF( wntua .AND. wntvas )
THEN
435 wrkbl = n + lwork_dgeqrf
436 wrkbl = max( wrkbl, n+lwork_dorgqr_m )
437 wrkbl = max( wrkbl, 3*n+lwork_dgebrd )
438 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_q )
439 wrkbl = max( wrkbl, 3*n+lwork_dorgbr_p )
440 wrkbl = max( wrkbl, bdspac )
442 minwrk = max( 3*n+m, bdspac )
448 CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
449 $ dum(1), dum(1), -1, ierr )
451 maxwrk = 3*n + lwork_dgebrd
452 IF( wntus .OR. wntuo )
THEN
453 CALL dorgbr(
'Q', m, n, n, a, lda, dum(1),
455 lwork_dorgbr_q=dum(1)
456 maxwrk = max( maxwrk, 3*n+lwork_dorgbr_q )
459 CALL dorgbr(
'Q', m, m, n, a, lda, dum(1),
461 lwork_dorgbr_q=dum(1)
462 maxwrk = max( maxwrk, 3*n+lwork_dorgbr_q )
464 IF( .NOT.wntvn )
THEN
465 maxwrk = max( maxwrk, 3*n+lwork_dorgbr_p )
467 maxwrk = max( maxwrk, bdspac )
468 minwrk = max( 3*n+m, bdspac )
470 ELSE IF( minmn.GT.0 )
THEN
474 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0, 0 )
477 CALL dgelqf( m, n, a, lda, dum(1), dum(1), -1, ierr )
480 CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
481 lwork_dorglq_n=dum(1)
482 CALL dorglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
483 lwork_dorglq_m=dum(1)
485 CALL dgebrd( m, m, a, lda, s, dum(1), dum(1),
486 $ dum(1), dum(1), -1, ierr )
489 CALL dorgbr(
'P', m, m, m, a, n, dum(1),
491 lwork_dorgbr_p=dum(1)
493 CALL dorgbr(
'Q', m, m, m, a, n, dum(1),
495 lwork_dorgbr_q=dum(1)
496 IF( n.GE.mnthr )
THEN
501 maxwrk = m + lwork_dgelqf
502 maxwrk = max( maxwrk, 3*m+lwork_dgebrd )
503 IF( wntuo .OR. wntuas )
504 $ maxwrk = max( maxwrk, 3*m+lwork_dorgbr_q )
505 maxwrk = max( maxwrk, bdspac )
506 minwrk = max( 4*m, bdspac )
507 ELSE IF( wntvo .AND. wntun )
THEN
511 wrkbl = m + lwork_dgelqf
512 wrkbl = max( wrkbl, m+lwork_dorglq_m )
513 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
514 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
515 wrkbl = max( wrkbl, bdspac )
516 maxwrk = max( m*m+wrkbl, m*m+m*n+m )
517 minwrk = max( 3*m+n, bdspac )
518 ELSE IF( wntvo .AND. wntuas )
THEN
523 wrkbl = m + lwork_dgelqf
524 wrkbl = max( wrkbl, m+lwork_dorglq_m )
525 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
526 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
527 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
528 wrkbl = max( wrkbl, bdspac )
529 maxwrk = max( m*m+wrkbl, m*m+m*n+m )
530 minwrk = max( 3*m+n, bdspac )
531 ELSE IF( wntvs .AND. wntun )
THEN
535 wrkbl = m + lwork_dgelqf
536 wrkbl = max( wrkbl, m+lwork_dorglq_m )
537 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
538 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
539 wrkbl = max( wrkbl, bdspac )
541 minwrk = max( 3*m+n, bdspac )
542 ELSE IF( wntvs .AND. wntuo )
THEN
546 wrkbl = m + lwork_dgelqf
547 wrkbl = max( wrkbl, m+lwork_dorglq_m )
548 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
549 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
550 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
551 wrkbl = max( wrkbl, bdspac )
552 maxwrk = 2*m*m + wrkbl
553 minwrk = max( 3*m+n, bdspac )
554 ELSE IF( wntvs .AND. wntuas )
THEN
559 wrkbl = m + lwork_dgelqf
560 wrkbl = max( wrkbl, m+lwork_dorglq_m )
561 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
562 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
563 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
564 wrkbl = max( wrkbl, bdspac )
566 minwrk = max( 3*m+n, bdspac )
567 ELSE IF( wntva .AND. wntun )
THEN
571 wrkbl = m + lwork_dgelqf
572 wrkbl = max( wrkbl, m+lwork_dorglq_n )
573 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
574 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
575 wrkbl = max( wrkbl, bdspac )
577 minwrk = max( 3*m+n, bdspac )
578 ELSE IF( wntva .AND. wntuo )
THEN
582 wrkbl = m + lwork_dgelqf
583 wrkbl = max( wrkbl, m+lwork_dorglq_n )
584 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
585 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
586 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
587 wrkbl = max( wrkbl, bdspac )
588 maxwrk = 2*m*m + wrkbl
589 minwrk = max( 3*m+n, bdspac )
590 ELSE IF( wntva .AND. wntuas )
THEN
595 wrkbl = m + lwork_dgelqf
596 wrkbl = max( wrkbl, m+lwork_dorglq_n )
597 wrkbl = max( wrkbl, 3*m+lwork_dgebrd )
598 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_p )
599 wrkbl = max( wrkbl, 3*m+lwork_dorgbr_q )
600 wrkbl = max( wrkbl, bdspac )
602 minwrk = max( 3*m+n, bdspac )
608 CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
609 $ dum(1), dum(1), -1, ierr )
611 maxwrk = 3*m + lwork_dgebrd
612 IF( wntvs .OR. wntvo )
THEN
614 CALL dorgbr(
'P', m, n, m, a, n, dum(1),
616 lwork_dorgbr_p=dum(1)
617 maxwrk = max( maxwrk, 3*m+lwork_dorgbr_p )
620 CALL dorgbr(
'P', n, n, m, a, n, dum(1),
622 lwork_dorgbr_p=dum(1)
623 maxwrk = max( maxwrk, 3*m+lwork_dorgbr_p )
625 IF( .NOT.wntun )
THEN
626 maxwrk = max( maxwrk, 3*m+lwork_dorgbr_q )
628 maxwrk = max( maxwrk, bdspac )
629 minwrk = max( 3*m+n, bdspac )
632 maxwrk = max( maxwrk, minwrk )
635 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
641 CALL xerbla(
'DGESVD', -info )
643 ELSE IF( lquery )
THEN
649 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
656 smlnum = sqrt( dlamch(
'S' ) ) / eps
657 bignum = one / smlnum
661 anrm = dlange(
'M', m, n, a, lda, dum )
663 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
665 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
666 ELSE IF( anrm.GT.bignum )
THEN
668 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
677 IF( m.GE.mnthr )
THEN
690 CALL dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
691 $ lwork-iwork+1, ierr )
695 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
704 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
705 $ work( itaup ), work( iwork ), lwork-iwork+1,
708 IF( wntvo .OR. wntvas )
THEN
713 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
714 $ work( iwork ), lwork-iwork+1, ierr )
723 CALL dbdsqr(
'U', n, ncvt, 0, 0, s, work( ie ), a, lda,
724 $ dum, 1, dum, 1, work( iwork ), info )
729 $
CALL dlacpy(
'F', n, n, a, lda, vt, ldvt )
731 ELSE IF( wntuo .AND. wntvn )
THEN
737 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
742 IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n )
THEN
748 ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n )
THEN
758 ldwrku = ( lwork-n*n-n ) / n
767 CALL dgeqrf( m, n, a, lda, work( itau ),
768 $ work( iwork ), lwork-iwork+1, ierr )
772 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
773 CALL dlaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
779 CALL dorgqr( m, n, n, a, lda, work( itau ),
780 $ work( iwork ), lwork-iwork+1, ierr )
789 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
790 $ work( itauq ), work( itaup ),
791 $ work( iwork ), lwork-iwork+1, ierr )
796 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
797 $ work( itauq ), work( iwork ),
798 $ lwork-iwork+1, ierr )
805 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum, 1,
806 $ work( ir ), ldwrkr, dum, 1,
807 $ work( iwork ), info )
814 DO 10 i = 1, m, ldwrku
815 chunk = min( m-i+1, ldwrku )
816 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
817 $ lda, work( ir ), ldwrkr, zero,
818 $ work( iu ), ldwrku )
819 CALL dlacpy(
'F', chunk, n, work( iu ), ldwrku,
835 CALL dgebrd( m, n, a, lda, s, work( ie ),
836 $ work( itauq ), work( itaup ),
837 $ work( iwork ), lwork-iwork+1, ierr )
842 CALL dorgbr(
'Q', m, n, n, a, lda, work( itauq ),
843 $ work( iwork ), lwork-iwork+1, ierr )
850 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum, 1,
851 $ a, lda, dum, 1, work( iwork ), info )
855 ELSE IF( wntuo .AND. wntvas )
THEN
861 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
866 IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n )
THEN
872 ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n )
THEN
882 ldwrku = ( lwork-n*n-n ) / n
891 CALL dgeqrf( m, n, a, lda, work( itau ),
892 $ work( iwork ), lwork-iwork+1, ierr )
896 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
898 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
904 CALL dorgqr( m, n, n, a, lda, work( itau ),
905 $ work( iwork ), lwork-iwork+1, ierr )
914 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
915 $ work( itauq ), work( itaup ),
916 $ work( iwork ), lwork-iwork+1, ierr )
917 CALL dlacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
922 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
923 $ work( itauq ), work( iwork ),
924 $ lwork-iwork+1, ierr )
929 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
930 $ work( iwork ), lwork-iwork+1, ierr )
938 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt, ldvt,
939 $ work( ir ), ldwrkr, dum, 1,
940 $ work( iwork ), info )
947 DO 20 i = 1, m, ldwrku
948 chunk = min( m-i+1, ldwrku )
949 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
950 $ lda, work( ir ), ldwrkr, zero,
951 $ work( iu ), ldwrku )
952 CALL dlacpy(
'F', chunk, n, work( iu ), ldwrku,
966 CALL dgeqrf( m, n, a, lda, work( itau ),
967 $ work( iwork ), lwork-iwork+1, ierr )
971 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
973 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
979 CALL dorgqr( m, n, n, a, lda, work( itau ),
980 $ work( iwork ), lwork-iwork+1, ierr )
989 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
990 $ work( itauq ), work( itaup ),
991 $ work( iwork ), lwork-iwork+1, ierr )
996 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
997 $ work( itauq ), a, lda, work( iwork ),
998 $ lwork-iwork+1, ierr )
1003 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1004 $ work( iwork ), lwork-iwork+1, ierr )
1012 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt, ldvt,
1013 $ a, lda, dum, 1, work( iwork ), info )
1017 ELSE IF( wntus )
THEN
1025 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1030 IF( lwork.GE.wrkbl+lda*n )
THEN
1041 itau = ir + ldwrkr*n
1047 CALL dgeqrf( m, n, a, lda, work( itau ),
1048 $ work( iwork ), lwork-iwork+1, ierr )
1052 CALL dlacpy(
'U', n, n, a, lda, work( ir ),
1054 CALL dlaset(
'L', n-1, n-1, zero, zero,
1055 $ work( ir+1 ), ldwrkr )
1060 CALL dorgqr( m, n, n, a, lda, work( itau ),
1061 $ work( iwork ), lwork-iwork+1, ierr )
1070 CALL dgebrd( n, n, work( ir ), ldwrkr, s,
1071 $ work( ie ), work( itauq ),
1072 $ work( itaup ), work( iwork ),
1073 $ lwork-iwork+1, ierr )
1078 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1079 $ work( itauq ), work( iwork ),
1080 $ lwork-iwork+1, ierr )
1087 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1088 $ 1, work( ir ), ldwrkr, dum, 1,
1089 $ work( iwork ), info )
1095 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1096 $ work( ir ), ldwrkr, zero, u, ldu )
1108 CALL dgeqrf( m, n, a, lda, work( itau ),
1109 $ work( iwork ), lwork-iwork+1, ierr )
1110 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1115 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1116 $ work( iwork ), lwork-iwork+1, ierr )
1124 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1130 CALL dgebrd( n, n, a, lda, s, work( ie ),
1131 $ work( itauq ), work( itaup ),
1132 $ work( iwork ), lwork-iwork+1, ierr )
1137 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1138 $ work( itauq ), u, ldu, work( iwork ),
1139 $ lwork-iwork+1, ierr )
1146 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1147 $ 1, u, ldu, dum, 1, work( iwork ),
1152 ELSE IF( wntvo )
THEN
1158 IF( lwork.GE.2*n*n+max( 4*n, bdspac ) )
THEN
1163 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1170 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1185 itau = ir + ldwrkr*n
1191 CALL dgeqrf( m, n, a, lda, work( itau ),
1192 $ work( iwork ), lwork-iwork+1, ierr )
1196 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1198 CALL dlaset(
'L', n-1, n-1, zero, zero,
1199 $ work( iu+1 ), ldwrku )
1204 CALL dorgqr( m, n, n, a, lda, work( itau ),
1205 $ work( iwork ), lwork-iwork+1, ierr )
1216 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1217 $ work( ie ), work( itauq ),
1218 $ work( itaup ), work( iwork ),
1219 $ lwork-iwork+1, ierr )
1220 CALL dlacpy(
'U', n, n, work( iu ), ldwrku,
1221 $ work( ir ), ldwrkr )
1226 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1227 $ work( itauq ), work( iwork ),
1228 $ lwork-iwork+1, ierr )
1234 CALL dorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1235 $ work( itaup ), work( iwork ),
1236 $ lwork-iwork+1, ierr )
1244 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ),
1245 $ work( ir ), ldwrkr, work( iu ),
1246 $ ldwrku, dum, 1, work( iwork ), info )
1252 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1253 $ work( iu ), ldwrku, zero, u, ldu )
1258 CALL dlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1271 CALL dgeqrf( m, n, a, lda, work( itau ),
1272 $ work( iwork ), lwork-iwork+1, ierr )
1273 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1278 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1279 $ work( iwork ), lwork-iwork+1, ierr )
1287 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1293 CALL dgebrd( n, n, a, lda, s, work( ie ),
1294 $ work( itauq ), work( itaup ),
1295 $ work( iwork ), lwork-iwork+1, ierr )
1300 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1301 $ work( itauq ), u, ldu, work( iwork ),
1302 $ lwork-iwork+1, ierr )
1307 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
1308 $ work( iwork ), lwork-iwork+1, ierr )
1316 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1317 $ lda, u, ldu, dum, 1, work( iwork ),
1322 ELSE IF( wntvas )
THEN
1329 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1334 IF( lwork.GE.wrkbl+lda*n )
THEN
1345 itau = iu + ldwrku*n
1351 CALL dgeqrf( m, n, a, lda, work( itau ),
1352 $ work( iwork ), lwork-iwork+1, ierr )
1356 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1358 CALL dlaset(
'L', n-1, n-1, zero, zero,
1359 $ work( iu+1 ), ldwrku )
1364 CALL dorgqr( m, n, n, a, lda, work( itau ),
1365 $ work( iwork ), lwork-iwork+1, ierr )
1374 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1375 $ work( ie ), work( itauq ),
1376 $ work( itaup ), work( iwork ),
1377 $ lwork-iwork+1, ierr )
1378 CALL dlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1384 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1385 $ work( itauq ), work( iwork ),
1386 $ lwork-iwork+1, ierr )
1392 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1393 $ work( iwork ), lwork-iwork+1, ierr )
1401 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1402 $ ldvt, work( iu ), ldwrku, dum, 1,
1403 $ work( iwork ), info )
1409 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1410 $ work( iu ), ldwrku, zero, u, ldu )
1422 CALL dgeqrf( m, n, a, lda, work( itau ),
1423 $ work( iwork ), lwork-iwork+1, ierr )
1424 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1429 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1430 $ work( iwork ), lwork-iwork+1, ierr )
1434 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
1436 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
1437 $ vt( 2, 1 ), ldvt )
1446 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1447 $ work( itauq ), work( itaup ),
1448 $ work( iwork ), lwork-iwork+1, ierr )
1454 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1455 $ work( itauq ), u, ldu, work( iwork ),
1456 $ lwork-iwork+1, ierr )
1461 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1462 $ work( iwork ), lwork-iwork+1, ierr )
1470 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1471 $ ldvt, u, ldu, dum, 1, work( iwork ),
1478 ELSE IF( wntua )
THEN
1486 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1491 IF( lwork.GE.wrkbl+lda*n )
THEN
1502 itau = ir + ldwrkr*n
1508 CALL dgeqrf( m, n, a, lda, work( itau ),
1509 $ work( iwork ), lwork-iwork+1, ierr )
1510 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1514 CALL dlacpy(
'U', n, n, a, lda, work( ir ),
1516 CALL dlaset(
'L', n-1, n-1, zero, zero,
1517 $ work( ir+1 ), ldwrkr )
1522 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1523 $ work( iwork ), lwork-iwork+1, ierr )
1532 CALL dgebrd( n, n, work( ir ), ldwrkr, s,
1533 $ work( ie ), work( itauq ),
1534 $ work( itaup ), work( iwork ),
1535 $ lwork-iwork+1, ierr )
1540 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1541 $ work( itauq ), work( iwork ),
1542 $ lwork-iwork+1, ierr )
1549 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1550 $ 1, work( ir ), ldwrkr, dum, 1,
1551 $ work( iwork ), info )
1557 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1558 $ work( ir ), ldwrkr, zero, a, lda )
1562 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1574 CALL dgeqrf( m, n, a, lda, work( itau ),
1575 $ work( iwork ), lwork-iwork+1, ierr )
1576 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1581 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1582 $ work( iwork ), lwork-iwork+1, ierr )
1590 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1596 CALL dgebrd( n, n, a, lda, s, work( ie ),
1597 $ work( itauq ), work( itaup ),
1598 $ work( iwork ), lwork-iwork+1, ierr )
1604 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1605 $ work( itauq ), u, ldu, work( iwork ),
1606 $ lwork-iwork+1, ierr )
1613 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1614 $ 1, u, ldu, dum, 1, work( iwork ),
1619 ELSE IF( wntvo )
THEN
1625 IF( lwork.GE.2*n*n+max( n+m, 4*n, bdspac ) )
THEN
1630 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1637 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1652 itau = ir + ldwrkr*n
1658 CALL dgeqrf( m, n, a, lda, work( itau ),
1659 $ work( iwork ), lwork-iwork+1, ierr )
1660 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1665 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1666 $ work( iwork ), lwork-iwork+1, ierr )
1670 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1672 CALL dlaset(
'L', n-1, n-1, zero, zero,
1673 $ work( iu+1 ), ldwrku )
1684 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1685 $ work( ie ), work( itauq ),
1686 $ work( itaup ), work( iwork ),
1687 $ lwork-iwork+1, ierr )
1688 CALL dlacpy(
'U', n, n, work( iu ), ldwrku,
1689 $ work( ir ), ldwrkr )
1694 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1695 $ work( itauq ), work( iwork ),
1696 $ lwork-iwork+1, ierr )
1702 CALL dorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1703 $ work( itaup ), work( iwork ),
1704 $ lwork-iwork+1, ierr )
1712 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ),
1713 $ work( ir ), ldwrkr, work( iu ),
1714 $ ldwrku, dum, 1, work( iwork ), info )
1720 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1721 $ work( iu ), ldwrku, zero, a, lda )
1725 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1729 CALL dlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1742 CALL dgeqrf( m, n, a, lda, work( itau ),
1743 $ work( iwork ), lwork-iwork+1, ierr )
1744 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1749 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1750 $ work( iwork ), lwork-iwork+1, ierr )
1758 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1764 CALL dgebrd( n, n, a, lda, s, work( ie ),
1765 $ work( itauq ), work( itaup ),
1766 $ work( iwork ), lwork-iwork+1, ierr )
1772 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1773 $ work( itauq ), u, ldu, work( iwork ),
1774 $ lwork-iwork+1, ierr )
1779 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
1780 $ work( iwork ), lwork-iwork+1, ierr )
1788 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1789 $ lda, u, ldu, dum, 1, work( iwork ),
1794 ELSE IF( wntvas )
THEN
1801 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1806 IF( lwork.GE.wrkbl+lda*n )
THEN
1817 itau = iu + ldwrku*n
1823 CALL dgeqrf( m, n, a, lda, work( itau ),
1824 $ work( iwork ), lwork-iwork+1, ierr )
1825 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1830 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1831 $ work( iwork ), lwork-iwork+1, ierr )
1835 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1837 CALL dlaset(
'L', n-1, n-1, zero, zero,
1838 $ work( iu+1 ), ldwrku )
1847 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1848 $ work( ie ), work( itauq ),
1849 $ work( itaup ), work( iwork ),
1850 $ lwork-iwork+1, ierr )
1851 CALL dlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1857 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1858 $ work( itauq ), work( iwork ),
1859 $ lwork-iwork+1, ierr )
1865 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1866 $ work( iwork ), lwork-iwork+1, ierr )
1874 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1875 $ ldvt, work( iu ), ldwrku, dum, 1,
1876 $ work( iwork ), info )
1882 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1883 $ work( iu ), ldwrku, zero, a, lda )
1887 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1899 CALL dgeqrf( m, n, a, lda, work( itau ),
1900 $ work( iwork ), lwork-iwork+1, ierr )
1901 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1906 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1907 $ work( iwork ), lwork-iwork+1, ierr )
1911 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
1913 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
1914 $ vt( 2, 1 ), ldvt )
1923 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1924 $ work( itauq ), work( itaup ),
1925 $ work( iwork ), lwork-iwork+1, ierr )
1931 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1932 $ work( itauq ), u, ldu, work( iwork ),
1933 $ lwork-iwork+1, ierr )
1938 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1939 $ work( iwork ), lwork-iwork+1, ierr )
1947 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1948 $ ldvt, u, ldu, dum, 1, work( iwork ),
1972 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1973 $ work( itaup ), work( iwork ), lwork-iwork+1,
1981 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1986 CALL dorgbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
1987 $ work( iwork ), lwork-iwork+1, ierr )
1995 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
1996 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1997 $ work( iwork ), lwork-iwork+1, ierr )
2005 CALL dorgbr(
'Q', m, n, n, a, lda, work( itauq ),
2006 $ work( iwork ), lwork-iwork+1, ierr )
2014 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
2015 $ work( iwork ), lwork-iwork+1, ierr )
2018 IF( wntuas .OR. wntuo )
2022 IF( wntvas .OR. wntvo )
2026 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2033 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2034 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
2035 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2042 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), a, lda,
2043 $ u, ldu, dum, 1, work( iwork ), info )
2051 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2052 $ ldvt, a, lda, dum, 1, work( iwork ), info )
2063 IF( n.GE.mnthr )
THEN
2076 CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
2077 $ lwork-iwork+1, ierr )
2081 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
2090 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
2091 $ work( itaup ), work( iwork ), lwork-iwork+1,
2093 IF( wntuo .OR. wntuas )
THEN
2098 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2099 $ work( iwork ), lwork-iwork+1, ierr )
2103 IF( wntuo .OR. wntuas )
2110 CALL dbdsqr(
'U', m, 0, nru, 0, s, work( ie ), dum, 1, a,
2111 $ lda, dum, 1, work( iwork ), info )
2116 $
CALL dlacpy(
'F', m, m, a, lda, u, ldu )
2118 ELSE IF( wntvo .AND. wntun )
THEN
2124 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2129 IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m )
THEN
2136 ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m )
THEN
2148 chunk = ( lwork-m*m-m ) / m
2151 itau = ir + ldwrkr*m
2157 CALL dgelqf( m, n, a, lda, work( itau ),
2158 $ work( iwork ), lwork-iwork+1, ierr )
2162 CALL dlacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2163 CALL dlaset(
'U', m-1, m-1, zero, zero,
2164 $ work( ir+ldwrkr ), ldwrkr )
2169 CALL dorglq( m, n, m, a, lda, work( itau ),
2170 $ work( iwork ), lwork-iwork+1, ierr )
2179 CALL dgebrd( m, m, work( ir ), ldwrkr, s, work( ie ),
2180 $ work( itauq ), work( itaup ),
2181 $ work( iwork ), lwork-iwork+1, ierr )
2186 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2187 $ work( itaup ), work( iwork ),
2188 $ lwork-iwork+1, ierr )
2195 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2196 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2197 $ work( iwork ), info )
2204 DO 30 i = 1, n, chunk
2205 blk = min( n-i+1, chunk )
2206 CALL dgemm(
'N',
'N', m, blk, m, one, work( ir ),
2207 $ ldwrkr, a( 1, i ), lda, zero,
2208 $ work( iu ), ldwrku )
2209 CALL dlacpy(
'F', m, blk, work( iu ), ldwrku,
2225 CALL dgebrd( m, n, a, lda, s, work( ie ),
2226 $ work( itauq ), work( itaup ),
2227 $ work( iwork ), lwork-iwork+1, ierr )
2232 CALL dorgbr(
'P', m, n, m, a, lda, work( itaup ),
2233 $ work( iwork ), lwork-iwork+1, ierr )
2240 CALL dbdsqr(
'L', m, n, 0, 0, s, work( ie ), a, lda,
2241 $ dum, 1, dum, 1, work( iwork ), info )
2245 ELSE IF( wntvo .AND. wntuas )
THEN
2251 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2256 IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m )
THEN
2263 ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m )
THEN
2275 chunk = ( lwork-m*m-m ) / m
2278 itau = ir + ldwrkr*m
2284 CALL dgelqf( m, n, a, lda, work( itau ),
2285 $ work( iwork ), lwork-iwork+1, ierr )
2289 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2290 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2296 CALL dorglq( m, n, m, a, lda, work( itau ),
2297 $ work( iwork ), lwork-iwork+1, ierr )
2306 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2307 $ work( itauq ), work( itaup ),
2308 $ work( iwork ), lwork-iwork+1, ierr )
2309 CALL dlacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2314 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2315 $ work( itaup ), work( iwork ),
2316 $ lwork-iwork+1, ierr )
2321 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2322 $ work( iwork ), lwork-iwork+1, ierr )
2330 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2331 $ work( ir ), ldwrkr, u, ldu, dum, 1,
2332 $ work( iwork ), info )
2339 DO 40 i = 1, n, chunk
2340 blk = min( n-i+1, chunk )
2341 CALL dgemm(
'N',
'N', m, blk, m, one, work( ir ),
2342 $ ldwrkr, a( 1, i ), lda, zero,
2343 $ work( iu ), ldwrku )
2344 CALL dlacpy(
'F', m, blk, work( iu ), ldwrku,
2358 CALL dgelqf( m, n, a, lda, work( itau ),
2359 $ work( iwork ), lwork-iwork+1, ierr )
2363 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2364 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2370 CALL dorglq( m, n, m, a, lda, work( itau ),
2371 $ work( iwork ), lwork-iwork+1, ierr )
2380 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2381 $ work( itauq ), work( itaup ),
2382 $ work( iwork ), lwork-iwork+1, ierr )
2387 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2388 $ work( itaup ), a, lda, work( iwork ),
2389 $ lwork-iwork+1, ierr )
2394 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2395 $ work( iwork ), lwork-iwork+1, ierr )
2403 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), a, lda,
2404 $ u, ldu, dum, 1, work( iwork ), info )
2408 ELSE IF( wntvs )
THEN
2416 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2421 IF( lwork.GE.wrkbl+lda*m )
THEN
2432 itau = ir + ldwrkr*m
2438 CALL dgelqf( m, n, a, lda, work( itau ),
2439 $ work( iwork ), lwork-iwork+1, ierr )
2443 CALL dlacpy(
'L', m, m, a, lda, work( ir ),
2445 CALL dlaset(
'U', m-1, m-1, zero, zero,
2446 $ work( ir+ldwrkr ), ldwrkr )
2451 CALL dorglq( m, n, m, a, lda, work( itau ),
2452 $ work( iwork ), lwork-iwork+1, ierr )
2461 CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2462 $ work( ie ), work( itauq ),
2463 $ work( itaup ), work( iwork ),
2464 $ lwork-iwork+1, ierr )
2470 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2471 $ work( itaup ), work( iwork ),
2472 $ lwork-iwork+1, ierr )
2479 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2480 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2481 $ work( iwork ), info )
2487 CALL dgemm(
'N',
'N', m, n, m, one, work( ir ),
2488 $ ldwrkr, a, lda, zero, vt, ldvt )
2500 CALL dgelqf( m, n, a, lda, work( itau ),
2501 $ work( iwork ), lwork-iwork+1, ierr )
2505 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2510 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2511 $ work( iwork ), lwork-iwork+1, ierr )
2519 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2525 CALL dgebrd( m, m, a, lda, s, work( ie ),
2526 $ work( itauq ), work( itaup ),
2527 $ work( iwork ), lwork-iwork+1, ierr )
2532 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
2533 $ work( itaup ), vt, ldvt,
2534 $ work( iwork ), lwork-iwork+1, ierr )
2541 CALL dbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
2542 $ ldvt, dum, 1, dum, 1, work( iwork ),
2547 ELSE IF( wntuo )
THEN
2553 IF( lwork.GE.2*m*m+max( 4*m, bdspac ) )
THEN
2558 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2565 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2580 itau = ir + ldwrkr*m
2586 CALL dgelqf( m, n, a, lda, work( itau ),
2587 $ work( iwork ), lwork-iwork+1, ierr )
2591 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
2593 CALL dlaset(
'U', m-1, m-1, zero, zero,
2594 $ work( iu+ldwrku ), ldwrku )
2599 CALL dorglq( m, n, m, a, lda, work( itau ),
2600 $ work( iwork ), lwork-iwork+1, ierr )
2611 CALL dgebrd( m, m, work( iu ), ldwrku, s,
2612 $ work( ie ), work( itauq ),
2613 $ work( itaup ), work( iwork ),
2614 $ lwork-iwork+1, ierr )
2615 CALL dlacpy(
'L', m, m, work( iu ), ldwrku,
2616 $ work( ir ), ldwrkr )
2622 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
2623 $ work( itaup ), work( iwork ),
2624 $ lwork-iwork+1, ierr )
2629 CALL dorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
2630 $ work( itauq ), work( iwork ),
2631 $ lwork-iwork+1, ierr )
2639 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2640 $ work( iu ), ldwrku, work( ir ),
2641 $ ldwrkr, dum, 1, work( iwork ), info )
2647 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
2648 $ ldwrku, a, lda, zero, vt, ldvt )
2653 CALL dlacpy(
'F', m, m, work( ir ), ldwrkr, a,
2666 CALL dgelqf( m, n, a, lda, work( itau ),
2667 $ work( iwork ), lwork-iwork+1, ierr )
2668 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2673 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2674 $ work( iwork ), lwork-iwork+1, ierr )
2682 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2688 CALL dgebrd( m, m, a, lda, s, work( ie ),
2689 $ work( itauq ), work( itaup ),
2690 $ work( iwork ), lwork-iwork+1, ierr )
2695 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
2696 $ work( itaup ), vt, ldvt,
2697 $ work( iwork ), lwork-iwork+1, ierr )
2702 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2703 $ work( iwork ), lwork-iwork+1, ierr )
2711 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2712 $ ldvt, a, lda, dum, 1, work( iwork ),
2717 ELSE IF( wntuas )
THEN
2724 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2729 IF( lwork.GE.wrkbl+lda*m )
THEN
2740 itau = iu + ldwrku*m
2746 CALL dgelqf( m, n, a, lda, work( itau ),
2747 $ work( iwork ), lwork-iwork+1, ierr )
2751 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
2753 CALL dlaset(
'U', m-1, m-1, zero, zero,
2754 $ work( iu+ldwrku ), ldwrku )
2759 CALL dorglq( m, n, m, a, lda, work( itau ),
2760 $ work( iwork ), lwork-iwork+1, ierr )
2769 CALL dgebrd( m, m, work( iu ), ldwrku, s,
2770 $ work( ie ), work( itauq ),
2771 $ work( itaup ), work( iwork ),
2772 $ lwork-iwork+1, ierr )
2773 CALL dlacpy(
'L', m, m, work( iu ), ldwrku, u,
2780 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
2781 $ work( itaup ), work( iwork ),
2782 $ lwork-iwork+1, ierr )
2787 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2788 $ work( iwork ), lwork-iwork+1, ierr )
2796 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2797 $ work( iu ), ldwrku, u, ldu, dum, 1,
2798 $ work( iwork ), info )
2804 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
2805 $ ldwrku, a, lda, zero, vt, ldvt )
2817 CALL dgelqf( m, n, a, lda, work( itau ),
2818 $ work( iwork ), lwork-iwork+1, ierr )
2819 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2824 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2825 $ work( iwork ), lwork-iwork+1, ierr )
2829 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2830 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2840 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2841 $ work( itauq ), work( itaup ),
2842 $ work( iwork ), lwork-iwork+1, ierr )
2848 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2849 $ work( itaup ), vt, ldvt,
2850 $ work( iwork ), lwork-iwork+1, ierr )
2855 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2856 $ work( iwork ), lwork-iwork+1, ierr )
2864 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2865 $ ldvt, u, ldu, dum, 1, work( iwork ),
2872 ELSE IF( wntva )
THEN
2880 IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) )
THEN
2885 IF( lwork.GE.wrkbl+lda*m )
THEN
2896 itau = ir + ldwrkr*m
2902 CALL dgelqf( m, n, a, lda, work( itau ),
2903 $ work( iwork ), lwork-iwork+1, ierr )
2904 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2908 CALL dlacpy(
'L', m, m, a, lda, work( ir ),
2910 CALL dlaset(
'U', m-1, m-1, zero, zero,
2911 $ work( ir+ldwrkr ), ldwrkr )
2916 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2917 $ work( iwork ), lwork-iwork+1, ierr )
2926 CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2927 $ work( ie ), work( itauq ),
2928 $ work( itaup ), work( iwork ),
2929 $ lwork-iwork+1, ierr )
2935 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2936 $ work( itaup ), work( iwork ),
2937 $ lwork-iwork+1, ierr )
2944 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2945 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2946 $ work( iwork ), info )
2952 CALL dgemm(
'N',
'N', m, n, m, one, work( ir ),
2953 $ ldwrkr, vt, ldvt, zero, a, lda )
2957 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
2969 CALL dgelqf( m, n, a, lda, work( itau ),
2970 $ work( iwork ), lwork-iwork+1, ierr )
2971 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2976 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2977 $ work( iwork ), lwork-iwork+1, ierr )
2985 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2991 CALL dgebrd( m, m, a, lda, s, work( ie ),
2992 $ work( itauq ), work( itaup ),
2993 $ work( iwork ), lwork-iwork+1, ierr )
2999 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
3000 $ work( itaup ), vt, ldvt,
3001 $ work( iwork ), lwork-iwork+1, ierr )
3008 CALL dbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
3009 $ ldvt, dum, 1, dum, 1, work( iwork ),
3014 ELSE IF( wntuo )
THEN
3020 IF( lwork.GE.2*m*m+max( n+m, 4*m, bdspac ) )
THEN
3025 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3032 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3047 itau = ir + ldwrkr*m
3053 CALL dgelqf( m, n, a, lda, work( itau ),
3054 $ work( iwork ), lwork-iwork+1, ierr )
3055 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3060 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3061 $ work( iwork ), lwork-iwork+1, ierr )
3065 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
3067 CALL dlaset(
'U', m-1, m-1, zero, zero,
3068 $ work( iu+ldwrku ), ldwrku )
3079 CALL dgebrd( m, m, work( iu ), ldwrku, s,
3080 $ work( ie ), work( itauq ),
3081 $ work( itaup ), work( iwork ),
3082 $ lwork-iwork+1, ierr )
3083 CALL dlacpy(
'L', m, m, work( iu ), ldwrku,
3084 $ work( ir ), ldwrkr )
3090 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
3091 $ work( itaup ), work( iwork ),
3092 $ lwork-iwork+1, ierr )
3097 CALL dorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
3098 $ work( itauq ), work( iwork ),
3099 $ lwork-iwork+1, ierr )
3107 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
3108 $ work( iu ), ldwrku, work( ir ),
3109 $ ldwrkr, dum, 1, work( iwork ), info )
3115 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
3116 $ ldwrku, vt, ldvt, zero, a, lda )
3120 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
3124 CALL dlacpy(
'F', m, m, work( ir ), ldwrkr, a,
3137 CALL dgelqf( m, n, a, lda, work( itau ),
3138 $ work( iwork ), lwork-iwork+1, ierr )
3139 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3144 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3145 $ work( iwork ), lwork-iwork+1, ierr )
3153 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
3159 CALL dgebrd( m, m, a, lda, s, work( ie ),
3160 $ work( itauq ), work( itaup ),
3161 $ work( iwork ), lwork-iwork+1, ierr )
3167 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
3168 $ work( itaup ), vt, ldvt,
3169 $ work( iwork ), lwork-iwork+1, ierr )
3174 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
3175 $ work( iwork ), lwork-iwork+1, ierr )
3183 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3184 $ ldvt, a, lda, dum, 1, work( iwork ),
3189 ELSE IF( wntuas )
THEN
3196 IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) )
THEN
3201 IF( lwork.GE.wrkbl+lda*m )
THEN
3212 itau = iu + ldwrku*m
3218 CALL dgelqf( m, n, a, lda, work( itau ),
3219 $ work( iwork ), lwork-iwork+1, ierr )
3220 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3225 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3226 $ work( iwork ), lwork-iwork+1, ierr )
3230 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
3232 CALL dlaset(
'U', m-1, m-1, zero, zero,
3233 $ work( iu+ldwrku ), ldwrku )
3242 CALL dgebrd( m, m, work( iu ), ldwrku, s,
3243 $ work( ie ), work( itauq ),
3244 $ work( itaup ), work( iwork ),
3245 $ lwork-iwork+1, ierr )
3246 CALL dlacpy(
'L', m, m, work( iu ), ldwrku, u,
3252 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
3253 $ work( itaup ), work( iwork ),
3254 $ lwork-iwork+1, ierr )
3259 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3260 $ work( iwork ), lwork-iwork+1, ierr )
3268 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
3269 $ work( iu ), ldwrku, u, ldu, dum, 1,
3270 $ work( iwork ), info )
3276 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
3277 $ ldwrku, vt, ldvt, zero, a, lda )
3281 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
3293 CALL dgelqf( m, n, a, lda, work( itau ),
3294 $ work( iwork ), lwork-iwork+1, ierr )
3295 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3300 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3301 $ work( iwork ), lwork-iwork+1, ierr )
3305 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
3306 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
3316 CALL dgebrd( m, m, u, ldu, s, work( ie ),
3317 $ work( itauq ), work( itaup ),
3318 $ work( iwork ), lwork-iwork+1, ierr )
3324 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
3325 $ work( itaup ), vt, ldvt,
3326 $ work( iwork ), lwork-iwork+1, ierr )
3331 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3332 $ work( iwork ), lwork-iwork+1, ierr )
3340 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3341 $ ldvt, u, ldu, dum, 1, work( iwork ),
3365 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
3366 $ work( itaup ), work( iwork ), lwork-iwork+1,
3374 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
3375 CALL dorgbr(
'Q', m, m, n, u, ldu, work( itauq ),
3376 $ work( iwork ), lwork-iwork+1, ierr )
3384 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3389 CALL dorgbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3390 $ work( iwork ), lwork-iwork+1, ierr )
3398 CALL dorgbr(
'Q', m, m, n, a, lda, work( itauq ),
3399 $ work( iwork ), lwork-iwork+1, ierr )
3407 CALL dorgbr(
'P', m, n, m, a, lda, work( itaup ),
3408 $ work( iwork ), lwork-iwork+1, ierr )
3411 IF( wntuas .OR. wntuo )
3415 IF( wntvas .OR. wntvo )
3419 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3426 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3427 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
3428 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3435 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), a, lda,
3436 $ u, ldu, dum, 1, work( iwork ), info )
3444 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3445 $ ldvt, a, lda, dum, 1, work( iwork ), info )
3455 IF( info.NE.0 )
THEN
3457 DO 50 i = 1, minmn - 1
3458 work( i+1 ) = work( i+ie-1 )
3462 DO 60 i = minmn - 1, 1, -1
3463 work( i+1 ) = work( i+ie-1 )
3470 IF( iscl.EQ.1 )
THEN
3471 IF( anrm.GT.bignum )
3472 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3474 IF( info.NE.0 .AND. anrm.GT.bignum )
3475 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn-1, 1, work( 2 ),
3477 IF( anrm.LT.smlnum )
3478 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3480 IF( info.NE.0 .AND. anrm.LT.smlnum )
3481 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1, work( 2 ),
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
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 dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR