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 )
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 )
317 lwork_dgeqrf = int( dum(1) )
319 CALL dorgqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
320 lwork_dorgqr_n = int( dum(1) )
321 CALL dorgqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
322 lwork_dorgqr_m = int( dum(1) )
324 CALL dgebrd( n, n, a, lda, s, dum(1), dum(1),
325 $ dum(1), dum(1), -1, ierr )
326 lwork_dgebrd = int( dum(1) )
328 CALL dorgbr(
'P', n, n, n, a, lda, dum(1),
330 lwork_dorgbr_p = int( dum(1) )
332 CALL dorgbr(
'Q', n, n, n, a, lda, dum(1),
334 lwork_dorgbr_q = int( 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 )
450 lwork_dgebrd = int( dum(1) )
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 = int( 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 = int( 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 )
478 lwork_dgelqf = int( dum(1) )
480 CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
481 lwork_dorglq_n = int( dum(1) )
482 CALL dorglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
483 lwork_dorglq_m = int( dum(1) )
485 CALL dgebrd( m, m, a, lda, s, dum(1), dum(1),
486 $ dum(1), dum(1), -1, ierr )
487 lwork_dgebrd = int( dum(1) )
489 CALL dorgbr(
'P', m, m, m, a, n, dum(1),
491 lwork_dorgbr_p = int( dum(1) )
493 CALL dorgbr(
'Q', m, m, m, a, n, dum(1),
495 lwork_dorgbr_q = int( 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 )
610 lwork_dgebrd = int( dum(1) )
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 = int( 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 = int( 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 )
696 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
707 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
708 $ work( itaup ), work( iwork ), lwork-iwork+1,
711 IF( wntvo .OR. wntvas )
THEN
716 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
717 $ work( iwork ), lwork-iwork+1, ierr )
726 CALL dbdsqr(
'U', n, ncvt, 0, 0, s, work( ie ), a, lda,
727 $ dum, 1, dum, 1, work( iwork ), info )
732 $
CALL dlacpy(
'F', n, n, a, lda, vt, ldvt )
734 ELSE IF( wntuo .AND. wntvn )
THEN
740 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
745 IF( lwork.GE.max( wrkbl, lda*n + n ) + lda*n )
THEN
751 ELSE IF( lwork.GE.max( wrkbl, lda*n + n ) + n*n )
THEN
761 ldwrku = ( lwork-n*n-n ) / n
770 CALL dgeqrf( m, n, a, lda, work( itau ),
771 $ work( iwork ), lwork-iwork+1, ierr )
775 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
776 CALL dlaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
782 CALL dorgqr( m, n, n, a, lda, work( itau ),
783 $ work( iwork ), lwork-iwork+1, ierr )
792 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
793 $ work( itauq ), work( itaup ),
794 $ work( iwork ), lwork-iwork+1, ierr )
799 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
800 $ work( itauq ), work( iwork ),
801 $ lwork-iwork+1, ierr )
808 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum, 1,
809 $ work( ir ), ldwrkr, dum, 1,
810 $ work( iwork ), info )
817 DO 10 i = 1, m, ldwrku
818 chunk = min( m-i+1, ldwrku )
819 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
820 $ lda, work( ir ), ldwrkr, zero,
821 $ work( iu ), ldwrku )
822 CALL dlacpy(
'F', chunk, n, work( iu ), ldwrku,
838 CALL dgebrd( m, n, a, lda, s, work( ie ),
839 $ work( itauq ), work( itaup ),
840 $ work( iwork ), lwork-iwork+1, ierr )
845 CALL dorgbr(
'Q', m, n, n, a, lda, work( itauq ),
846 $ work( iwork ), lwork-iwork+1, ierr )
853 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum, 1,
854 $ a, lda, dum, 1, work( iwork ), info )
858 ELSE IF( wntuo .AND. wntvas )
THEN
864 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
869 IF( lwork.GE.max( wrkbl, lda*n + n ) + lda*n )
THEN
875 ELSE IF( lwork.GE.max( wrkbl, lda*n + n ) + n*n )
THEN
885 ldwrku = ( lwork-n*n-n ) / n
894 CALL dgeqrf( m, n, a, lda, work( itau ),
895 $ work( iwork ), lwork-iwork+1, ierr )
899 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
901 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
907 CALL dorgqr( m, n, n, a, lda, work( itau ),
908 $ work( iwork ), lwork-iwork+1, ierr )
917 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
918 $ work( itauq ), work( itaup ),
919 $ work( iwork ), lwork-iwork+1, ierr )
920 CALL dlacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
925 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
926 $ work( itauq ), work( iwork ),
927 $ lwork-iwork+1, ierr )
932 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
933 $ work( iwork ), lwork-iwork+1, ierr )
941 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt, ldvt,
942 $ work( ir ), ldwrkr, dum, 1,
943 $ work( iwork ), info )
950 DO 20 i = 1, m, ldwrku
951 chunk = min( m-i+1, ldwrku )
952 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
953 $ lda, work( ir ), ldwrkr, zero,
954 $ work( iu ), ldwrku )
955 CALL dlacpy(
'F', chunk, n, work( iu ), ldwrku,
969 CALL dgeqrf( m, n, a, lda, work( itau ),
970 $ work( iwork ), lwork-iwork+1, ierr )
974 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
976 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
982 CALL dorgqr( m, n, n, a, lda, work( itau ),
983 $ work( iwork ), lwork-iwork+1, ierr )
992 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
993 $ work( itauq ), work( itaup ),
994 $ work( iwork ), lwork-iwork+1, ierr )
999 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1000 $ work( itauq ), a, lda, work( iwork ),
1001 $ lwork-iwork+1, ierr )
1006 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1007 $ work( iwork ), lwork-iwork+1, ierr )
1015 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt, ldvt,
1016 $ a, lda, dum, 1, work( iwork ), info )
1020 ELSE IF( wntus )
THEN
1028 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1033 IF( lwork.GE.wrkbl+lda*n )
THEN
1044 itau = ir + ldwrkr*n
1050 CALL dgeqrf( m, n, a, lda, work( itau ),
1051 $ work( iwork ), lwork-iwork+1, ierr )
1055 CALL dlacpy(
'U', n, n, a, lda, work( ir ),
1057 CALL dlaset(
'L', n-1, n-1, zero, zero,
1058 $ work( ir+1 ), ldwrkr )
1063 CALL dorgqr( m, n, n, a, lda, work( itau ),
1064 $ work( iwork ), lwork-iwork+1, ierr )
1073 CALL dgebrd( n, n, work( ir ), ldwrkr, s,
1074 $ work( ie ), work( itauq ),
1075 $ work( itaup ), work( iwork ),
1076 $ lwork-iwork+1, ierr )
1081 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1082 $ work( itauq ), work( iwork ),
1083 $ lwork-iwork+1, ierr )
1090 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1091 $ 1, work( ir ), ldwrkr, dum, 1,
1092 $ work( iwork ), info )
1098 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1099 $ work( ir ), ldwrkr, zero, u, ldu )
1111 CALL dgeqrf( m, n, a, lda, work( itau ),
1112 $ work( iwork ), lwork-iwork+1, ierr )
1113 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1118 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1119 $ work( iwork ), lwork-iwork+1, ierr )
1128 CALL dlaset(
'L', n-1, n-1, zero, zero,
1135 CALL dgebrd( n, n, a, lda, s, work( ie ),
1136 $ work( itauq ), work( itaup ),
1137 $ work( iwork ), lwork-iwork+1, ierr )
1142 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1143 $ work( itauq ), u, ldu, work( iwork ),
1144 $ lwork-iwork+1, ierr )
1151 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1152 $ 1, u, ldu, dum, 1, work( iwork ),
1157 ELSE IF( wntvo )
THEN
1163 IF( lwork.GE.2*n*n+max( 4*n, bdspac ) )
THEN
1168 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1175 ELSE IF( lwork.GE.wrkbl+( lda + n )*n )
THEN
1190 itau = ir + ldwrkr*n
1196 CALL dgeqrf( m, n, a, lda, work( itau ),
1197 $ work( iwork ), lwork-iwork+1, ierr )
1201 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1203 CALL dlaset(
'L', n-1, n-1, zero, zero,
1204 $ work( iu+1 ), ldwrku )
1209 CALL dorgqr( m, n, n, a, lda, work( itau ),
1210 $ work( iwork ), lwork-iwork+1, ierr )
1221 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1222 $ work( ie ), work( itauq ),
1223 $ work( itaup ), work( iwork ),
1224 $ lwork-iwork+1, ierr )
1225 CALL dlacpy(
'U', n, n, work( iu ), ldwrku,
1226 $ work( ir ), ldwrkr )
1231 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1232 $ work( itauq ), work( iwork ),
1233 $ lwork-iwork+1, ierr )
1239 CALL dorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1240 $ work( itaup ), work( iwork ),
1241 $ lwork-iwork+1, ierr )
1249 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ),
1250 $ work( ir ), ldwrkr, work( iu ),
1251 $ ldwrku, dum, 1, work( iwork ), info )
1257 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1258 $ work( iu ), ldwrku, zero, u, ldu )
1263 CALL dlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1276 CALL dgeqrf( m, n, a, lda, work( itau ),
1277 $ work( iwork ), lwork-iwork+1, ierr )
1278 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1283 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1284 $ work( iwork ), lwork-iwork+1, ierr )
1293 CALL dlaset(
'L', n-1, n-1, zero, zero,
1300 CALL dgebrd( n, n, a, lda, s, work( ie ),
1301 $ work( itauq ), work( itaup ),
1302 $ work( iwork ), lwork-iwork+1, ierr )
1307 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1308 $ work( itauq ), u, ldu, work( iwork ),
1309 $ lwork-iwork+1, ierr )
1314 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
1315 $ work( iwork ), lwork-iwork+1, ierr )
1323 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1324 $ lda, u, ldu, dum, 1, work( iwork ),
1329 ELSE IF( wntvas )
THEN
1336 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1341 IF( lwork.GE.wrkbl+lda*n )
THEN
1352 itau = iu + ldwrku*n
1358 CALL dgeqrf( m, n, a, lda, work( itau ),
1359 $ work( iwork ), lwork-iwork+1, ierr )
1363 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1365 CALL dlaset(
'L', n-1, n-1, zero, zero,
1366 $ work( iu+1 ), ldwrku )
1371 CALL dorgqr( m, n, n, a, lda, work( itau ),
1372 $ work( iwork ), lwork-iwork+1, ierr )
1381 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1382 $ work( ie ), work( itauq ),
1383 $ work( itaup ), work( iwork ),
1384 $ lwork-iwork+1, ierr )
1385 CALL dlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1391 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1392 $ work( itauq ), work( iwork ),
1393 $ lwork-iwork+1, ierr )
1399 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1400 $ work( iwork ), lwork-iwork+1, ierr )
1408 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1409 $ ldvt, work( iu ), ldwrku, dum, 1,
1410 $ work( iwork ), info )
1416 CALL dgemm(
'N',
'N', m, n, n, one, a, lda,
1417 $ work( iu ), ldwrku, zero, u, ldu )
1429 CALL dgeqrf( m, n, a, lda, work( itau ),
1430 $ work( iwork ), lwork-iwork+1, ierr )
1431 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1436 CALL dorgqr( m, n, n, u, ldu, work( itau ),
1437 $ work( iwork ), lwork-iwork+1, ierr )
1441 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
1443 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
1444 $ vt( 2, 1 ), ldvt )
1453 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1454 $ work( itauq ), work( itaup ),
1455 $ work( iwork ), lwork-iwork+1, ierr )
1461 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1462 $ work( itauq ), u, ldu, work( iwork ),
1463 $ lwork-iwork+1, ierr )
1468 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1469 $ work( iwork ), lwork-iwork+1, ierr )
1477 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1478 $ ldvt, u, ldu, dum, 1, work( iwork ),
1485 ELSE IF( wntua )
THEN
1493 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1498 IF( lwork.GE.wrkbl+lda*n )
THEN
1509 itau = ir + ldwrkr*n
1515 CALL dgeqrf( m, n, a, lda, work( itau ),
1516 $ work( iwork ), lwork-iwork+1, ierr )
1517 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1521 CALL dlacpy(
'U', n, n, a, lda, work( ir ),
1523 CALL dlaset(
'L', n-1, n-1, zero, zero,
1524 $ work( ir+1 ), ldwrkr )
1529 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1530 $ work( iwork ), lwork-iwork+1, ierr )
1539 CALL dgebrd( n, n, work( ir ), ldwrkr, s,
1540 $ work( ie ), work( itauq ),
1541 $ work( itaup ), work( iwork ),
1542 $ lwork-iwork+1, ierr )
1547 CALL dorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1548 $ work( itauq ), work( iwork ),
1549 $ lwork-iwork+1, ierr )
1556 CALL dbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1557 $ 1, work( ir ), ldwrkr, dum, 1,
1558 $ work( iwork ), info )
1564 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1565 $ work( ir ), ldwrkr, zero, a, lda )
1569 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1581 CALL dgeqrf( m, n, a, lda, work( itau ),
1582 $ work( iwork ), lwork-iwork+1, ierr )
1583 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1588 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1589 $ work( iwork ), lwork-iwork+1, ierr )
1598 CALL dlaset(
'L', n-1, n-1, zero, zero,
1605 CALL dgebrd( n, n, a, lda, s, work( ie ),
1606 $ work( itauq ), work( itaup ),
1607 $ work( iwork ), lwork-iwork+1, ierr )
1613 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1614 $ work( itauq ), u, ldu, work( iwork ),
1615 $ lwork-iwork+1, ierr )
1622 CALL dbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1623 $ 1, u, ldu, dum, 1, work( iwork ),
1628 ELSE IF( wntvo )
THEN
1634 IF( lwork.GE.2*n*n+max( n+m, 4*n, bdspac ) )
THEN
1639 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1646 ELSE IF( lwork.GE.wrkbl+( lda + n )*n )
THEN
1661 itau = ir + ldwrkr*n
1667 CALL dgeqrf( m, n, a, lda, work( itau ),
1668 $ work( iwork ), lwork-iwork+1, ierr )
1669 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1674 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1675 $ work( iwork ), lwork-iwork+1, ierr )
1679 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1681 CALL dlaset(
'L', n-1, n-1, zero, zero,
1682 $ work( iu+1 ), ldwrku )
1693 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1694 $ work( ie ), work( itauq ),
1695 $ work( itaup ), work( iwork ),
1696 $ lwork-iwork+1, ierr )
1697 CALL dlacpy(
'U', n, n, work( iu ), ldwrku,
1698 $ work( ir ), ldwrkr )
1703 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1704 $ work( itauq ), work( iwork ),
1705 $ lwork-iwork+1, ierr )
1711 CALL dorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1712 $ work( itaup ), work( iwork ),
1713 $ lwork-iwork+1, ierr )
1721 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ),
1722 $ work( ir ), ldwrkr, work( iu ),
1723 $ ldwrku, dum, 1, work( iwork ), info )
1729 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1730 $ work( iu ), ldwrku, zero, a, lda )
1734 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1738 CALL dlacpy(
'F', n, n, work( ir ), ldwrkr, a,
1751 CALL dgeqrf( m, n, a, lda, work( itau ),
1752 $ work( iwork ), lwork-iwork+1, ierr )
1753 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1758 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1759 $ work( iwork ), lwork-iwork+1, ierr )
1768 CALL dlaset(
'L', n-1, n-1, zero, zero,
1775 CALL dgebrd( n, n, a, lda, s, work( ie ),
1776 $ work( itauq ), work( itaup ),
1777 $ work( iwork ), lwork-iwork+1, ierr )
1783 CALL dormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1784 $ work( itauq ), u, ldu, work( iwork ),
1785 $ lwork-iwork+1, ierr )
1790 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
1791 $ work( iwork ), lwork-iwork+1, ierr )
1799 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1800 $ lda, u, ldu, dum, 1, work( iwork ),
1805 ELSE IF( wntvas )
THEN
1812 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1817 IF( lwork.GE.wrkbl+lda*n )
THEN
1828 itau = iu + ldwrku*n
1834 CALL dgeqrf( m, n, a, lda, work( itau ),
1835 $ work( iwork ), lwork-iwork+1, ierr )
1836 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1841 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1842 $ work( iwork ), lwork-iwork+1, ierr )
1846 CALL dlacpy(
'U', n, n, a, lda, work( iu ),
1848 CALL dlaset(
'L', n-1, n-1, zero, zero,
1849 $ work( iu+1 ), ldwrku )
1858 CALL dgebrd( n, n, work( iu ), ldwrku, s,
1859 $ work( ie ), work( itauq ),
1860 $ work( itaup ), work( iwork ),
1861 $ lwork-iwork+1, ierr )
1862 CALL dlacpy(
'U', n, n, work( iu ), ldwrku, vt,
1868 CALL dorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1869 $ work( itauq ), work( iwork ),
1870 $ lwork-iwork+1, ierr )
1876 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1877 $ work( iwork ), lwork-iwork+1, ierr )
1885 CALL dbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1886 $ ldvt, work( iu ), ldwrku, dum, 1,
1887 $ work( iwork ), info )
1893 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu,
1894 $ work( iu ), ldwrku, zero, a, lda )
1898 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
1910 CALL dgeqrf( m, n, a, lda, work( itau ),
1911 $ work( iwork ), lwork-iwork+1, ierr )
1912 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1917 CALL dorgqr( m, m, n, u, ldu, work( itau ),
1918 $ work( iwork ), lwork-iwork+1, ierr )
1922 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
1924 $
CALL dlaset(
'L', n-1, n-1, zero, zero,
1925 $ vt( 2, 1 ), ldvt )
1934 CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1935 $ work( itauq ), work( itaup ),
1936 $ work( iwork ), lwork-iwork+1, ierr )
1942 CALL dormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1943 $ work( itauq ), u, ldu, work( iwork ),
1944 $ lwork-iwork+1, ierr )
1949 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1950 $ work( iwork ), lwork-iwork+1, ierr )
1958 CALL dbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1959 $ ldvt, u, ldu, dum, 1, work( iwork ),
1983 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1984 $ work( itaup ), work( iwork ), lwork-iwork+1,
1992 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
1997 CALL dorgbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
1998 $ work( iwork ), lwork-iwork+1, ierr )
2006 CALL dlacpy(
'U', n, n, a, lda, vt, ldvt )
2007 CALL dorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2008 $ work( iwork ), lwork-iwork+1, ierr )
2016 CALL dorgbr(
'Q', m, n, n, a, lda, work( itauq ),
2017 $ work( iwork ), lwork-iwork+1, ierr )
2025 CALL dorgbr(
'P', n, n, n, a, lda, work( itaup ),
2026 $ work( iwork ), lwork-iwork+1, ierr )
2029 IF( wntuas .OR. wntuo )
2033 IF( wntvas .OR. wntvo )
2037 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2044 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2045 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
2046 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2053 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), a, lda,
2054 $ u, ldu, dum, 1, work( iwork ), info )
2062 CALL dbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2063 $ ldvt, a, lda, dum, 1, work( iwork ), info )
2074 IF( n.GE.mnthr )
THEN
2087 CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
2088 $ lwork-iwork+1, ierr )
2092 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
2101 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
2102 $ work( itaup ), work( iwork ), lwork-iwork+1,
2104 IF( wntuo .OR. wntuas )
THEN
2109 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2110 $ work( iwork ), lwork-iwork+1, ierr )
2114 IF( wntuo .OR. wntuas )
2121 CALL dbdsqr(
'U', m, 0, nru, 0, s, work( ie ), dum, 1, a,
2122 $ lda, dum, 1, work( iwork ), info )
2127 $
CALL dlacpy(
'F', m, m, a, lda, u, ldu )
2129 ELSE IF( wntvo .AND. wntun )
THEN
2135 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2140 IF( lwork.GE.max( wrkbl, lda*n + m ) + lda*m )
THEN
2147 ELSE IF( lwork.GE.max( wrkbl, lda*n + m ) + m*m )
THEN
2159 chunk = ( lwork-m*m-m ) / m
2162 itau = ir + ldwrkr*m
2168 CALL dgelqf( m, n, a, lda, work( itau ),
2169 $ work( iwork ), lwork-iwork+1, ierr )
2173 CALL dlacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2174 CALL dlaset(
'U', m-1, m-1, zero, zero,
2175 $ work( ir+ldwrkr ), ldwrkr )
2180 CALL dorglq( m, n, m, a, lda, work( itau ),
2181 $ work( iwork ), lwork-iwork+1, ierr )
2190 CALL dgebrd( m, m, work( ir ), ldwrkr, s, work( ie ),
2191 $ work( itauq ), work( itaup ),
2192 $ work( iwork ), lwork-iwork+1, ierr )
2197 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2198 $ work( itaup ), work( iwork ),
2199 $ lwork-iwork+1, ierr )
2206 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2207 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2208 $ work( iwork ), info )
2215 DO 30 i = 1, n, chunk
2216 blk = min( n-i+1, chunk )
2217 CALL dgemm(
'N',
'N', m, blk, m, one, work( ir ),
2218 $ ldwrkr, a( 1, i ), lda, zero,
2219 $ work( iu ), ldwrku )
2220 CALL dlacpy(
'F', m, blk, work( iu ), ldwrku,
2236 CALL dgebrd( m, n, a, lda, s, work( ie ),
2237 $ work( itauq ), work( itaup ),
2238 $ work( iwork ), lwork-iwork+1, ierr )
2243 CALL dorgbr(
'P', m, n, m, a, lda, work( itaup ),
2244 $ work( iwork ), lwork-iwork+1, ierr )
2251 CALL dbdsqr(
'L', m, n, 0, 0, s, work( ie ), a, lda,
2252 $ dum, 1, dum, 1, work( iwork ), info )
2256 ELSE IF( wntvo .AND. wntuas )
THEN
2262 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2267 IF( lwork.GE.max( wrkbl, lda*n + m ) + lda*m )
THEN
2274 ELSE IF( lwork.GE.max( wrkbl, lda*n + m ) + m*m )
THEN
2286 chunk = ( lwork-m*m-m ) / m
2289 itau = ir + ldwrkr*m
2295 CALL dgelqf( m, n, a, lda, work( itau ),
2296 $ work( iwork ), lwork-iwork+1, ierr )
2300 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2301 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2307 CALL dorglq( m, n, m, a, lda, work( itau ),
2308 $ work( iwork ), lwork-iwork+1, ierr )
2317 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2318 $ work( itauq ), work( itaup ),
2319 $ work( iwork ), lwork-iwork+1, ierr )
2320 CALL dlacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2325 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2326 $ work( itaup ), work( iwork ),
2327 $ lwork-iwork+1, ierr )
2332 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2333 $ work( iwork ), lwork-iwork+1, ierr )
2341 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2342 $ work( ir ), ldwrkr, u, ldu, dum, 1,
2343 $ work( iwork ), info )
2350 DO 40 i = 1, n, chunk
2351 blk = min( n-i+1, chunk )
2352 CALL dgemm(
'N',
'N', m, blk, m, one, work( ir ),
2353 $ ldwrkr, a( 1, i ), lda, zero,
2354 $ work( iu ), ldwrku )
2355 CALL dlacpy(
'F', m, blk, work( iu ), ldwrku,
2369 CALL dgelqf( m, n, a, lda, work( itau ),
2370 $ work( iwork ), lwork-iwork+1, ierr )
2374 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2375 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2381 CALL dorglq( m, n, m, a, lda, work( itau ),
2382 $ work( iwork ), lwork-iwork+1, ierr )
2391 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2392 $ work( itauq ), work( itaup ),
2393 $ work( iwork ), lwork-iwork+1, ierr )
2398 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2399 $ work( itaup ), a, lda, work( iwork ),
2400 $ lwork-iwork+1, ierr )
2405 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2406 $ work( iwork ), lwork-iwork+1, ierr )
2414 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), a, lda,
2415 $ u, ldu, dum, 1, work( iwork ), info )
2419 ELSE IF( wntvs )
THEN
2427 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2432 IF( lwork.GE.wrkbl+lda*m )
THEN
2443 itau = ir + ldwrkr*m
2449 CALL dgelqf( m, n, a, lda, work( itau ),
2450 $ work( iwork ), lwork-iwork+1, ierr )
2454 CALL dlacpy(
'L', m, m, a, lda, work( ir ),
2456 CALL dlaset(
'U', m-1, m-1, zero, zero,
2457 $ work( ir+ldwrkr ), ldwrkr )
2462 CALL dorglq( m, n, m, a, lda, work( itau ),
2463 $ work( iwork ), lwork-iwork+1, ierr )
2472 CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2473 $ work( ie ), work( itauq ),
2474 $ work( itaup ), work( iwork ),
2475 $ lwork-iwork+1, ierr )
2481 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2482 $ work( itaup ), work( iwork ),
2483 $ lwork-iwork+1, ierr )
2490 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2491 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2492 $ work( iwork ), info )
2498 CALL dgemm(
'N',
'N', m, n, m, one, work( ir ),
2499 $ ldwrkr, a, lda, zero, vt, ldvt )
2511 CALL dgelqf( m, n, a, lda, work( itau ),
2512 $ work( iwork ), lwork-iwork+1, ierr )
2516 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2521 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2522 $ work( iwork ), lwork-iwork+1, ierr )
2530 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2536 CALL dgebrd( m, m, a, lda, s, work( ie ),
2537 $ work( itauq ), work( itaup ),
2538 $ work( iwork ), lwork-iwork+1, ierr )
2543 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
2544 $ work( itaup ), vt, ldvt,
2545 $ work( iwork ), lwork-iwork+1, ierr )
2552 CALL dbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
2553 $ ldvt, dum, 1, dum, 1, work( iwork ),
2558 ELSE IF( wntuo )
THEN
2564 IF( lwork.GE.2*m*m+max( 4*m, bdspac ) )
THEN
2569 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2576 ELSE IF( lwork.GE.wrkbl+( lda + m )*m )
THEN
2591 itau = ir + ldwrkr*m
2597 CALL dgelqf( m, n, a, lda, work( itau ),
2598 $ work( iwork ), lwork-iwork+1, ierr )
2602 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
2604 CALL dlaset(
'U', m-1, m-1, zero, zero,
2605 $ work( iu+ldwrku ), ldwrku )
2610 CALL dorglq( m, n, m, a, lda, work( itau ),
2611 $ work( iwork ), lwork-iwork+1, ierr )
2622 CALL dgebrd( m, m, work( iu ), ldwrku, s,
2623 $ work( ie ), work( itauq ),
2624 $ work( itaup ), work( iwork ),
2625 $ lwork-iwork+1, ierr )
2626 CALL dlacpy(
'L', m, m, work( iu ), ldwrku,
2627 $ work( ir ), ldwrkr )
2633 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
2634 $ work( itaup ), work( iwork ),
2635 $ lwork-iwork+1, ierr )
2640 CALL dorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
2641 $ work( itauq ), work( iwork ),
2642 $ lwork-iwork+1, ierr )
2650 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2651 $ work( iu ), ldwrku, work( ir ),
2652 $ ldwrkr, dum, 1, work( iwork ), info )
2658 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
2659 $ ldwrku, a, lda, zero, vt, ldvt )
2664 CALL dlacpy(
'F', m, m, work( ir ), ldwrkr, a,
2677 CALL dgelqf( m, n, a, lda, work( itau ),
2678 $ work( iwork ), lwork-iwork+1, ierr )
2679 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2684 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2685 $ work( iwork ), lwork-iwork+1, ierr )
2693 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2699 CALL dgebrd( m, m, a, lda, s, work( ie ),
2700 $ work( itauq ), work( itaup ),
2701 $ work( iwork ), lwork-iwork+1, ierr )
2706 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
2707 $ work( itaup ), vt, ldvt,
2708 $ work( iwork ), lwork-iwork+1, ierr )
2713 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2714 $ work( iwork ), lwork-iwork+1, ierr )
2722 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2723 $ ldvt, a, lda, dum, 1, work( iwork ),
2728 ELSE IF( wntuas )
THEN
2735 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2740 IF( lwork.GE.wrkbl+lda*m )
THEN
2751 itau = iu + ldwrku*m
2757 CALL dgelqf( m, n, a, lda, work( itau ),
2758 $ work( iwork ), lwork-iwork+1, ierr )
2762 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
2764 CALL dlaset(
'U', m-1, m-1, zero, zero,
2765 $ work( iu+ldwrku ), ldwrku )
2770 CALL dorglq( m, n, m, a, lda, work( itau ),
2771 $ work( iwork ), lwork-iwork+1, ierr )
2780 CALL dgebrd( m, m, work( iu ), ldwrku, s,
2781 $ work( ie ), work( itauq ),
2782 $ work( itaup ), work( iwork ),
2783 $ lwork-iwork+1, ierr )
2784 CALL dlacpy(
'L', m, m, work( iu ), ldwrku, u,
2791 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
2792 $ work( itaup ), work( iwork ),
2793 $ lwork-iwork+1, ierr )
2798 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2799 $ work( iwork ), lwork-iwork+1, ierr )
2807 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
2808 $ work( iu ), ldwrku, u, ldu, dum, 1,
2809 $ work( iwork ), info )
2815 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
2816 $ ldwrku, a, lda, zero, vt, ldvt )
2828 CALL dgelqf( m, n, a, lda, work( itau ),
2829 $ work( iwork ), lwork-iwork+1, ierr )
2830 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2835 CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2836 $ work( iwork ), lwork-iwork+1, ierr )
2840 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
2841 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2851 CALL dgebrd( m, m, u, ldu, s, work( ie ),
2852 $ work( itauq ), work( itaup ),
2853 $ work( iwork ), lwork-iwork+1, ierr )
2859 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2860 $ work( itaup ), vt, ldvt,
2861 $ work( iwork ), lwork-iwork+1, ierr )
2866 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2867 $ work( iwork ), lwork-iwork+1, ierr )
2875 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2876 $ ldvt, u, ldu, dum, 1, work( iwork ),
2883 ELSE IF( wntva )
THEN
2891 IF( lwork.GE.m*m+max( n + m, 4*m, bdspac ) )
THEN
2896 IF( lwork.GE.wrkbl+lda*m )
THEN
2907 itau = ir + ldwrkr*m
2913 CALL dgelqf( m, n, a, lda, work( itau ),
2914 $ work( iwork ), lwork-iwork+1, ierr )
2915 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2919 CALL dlacpy(
'L', m, m, a, lda, work( ir ),
2921 CALL dlaset(
'U', m-1, m-1, zero, zero,
2922 $ work( ir+ldwrkr ), ldwrkr )
2927 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2928 $ work( iwork ), lwork-iwork+1, ierr )
2937 CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2938 $ work( ie ), work( itauq ),
2939 $ work( itaup ), work( iwork ),
2940 $ lwork-iwork+1, ierr )
2946 CALL dorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2947 $ work( itaup ), work( iwork ),
2948 $ lwork-iwork+1, ierr )
2955 CALL dbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2956 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2957 $ work( iwork ), info )
2963 CALL dgemm(
'N',
'N', m, n, m, one, work( ir ),
2964 $ ldwrkr, vt, ldvt, zero, a, lda )
2968 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
2980 CALL dgelqf( m, n, a, lda, work( itau ),
2981 $ work( iwork ), lwork-iwork+1, ierr )
2982 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
2987 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2988 $ work( iwork ), lwork-iwork+1, ierr )
2996 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
3002 CALL dgebrd( m, m, a, lda, s, work( ie ),
3003 $ work( itauq ), work( itaup ),
3004 $ work( iwork ), lwork-iwork+1, ierr )
3010 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
3011 $ work( itaup ), vt, ldvt,
3012 $ work( iwork ), lwork-iwork+1, ierr )
3019 CALL dbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
3020 $ ldvt, dum, 1, dum, 1, work( iwork ),
3025 ELSE IF( wntuo )
THEN
3031 IF( lwork.GE.2*m*m+max( n + m, 4*m, bdspac ) )
THEN
3036 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3043 ELSE IF( lwork.GE.wrkbl+( lda + m )*m )
THEN
3058 itau = ir + ldwrkr*m
3064 CALL dgelqf( m, n, a, lda, work( itau ),
3065 $ work( iwork ), lwork-iwork+1, ierr )
3066 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3071 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3072 $ work( iwork ), lwork-iwork+1, ierr )
3076 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
3078 CALL dlaset(
'U', m-1, m-1, zero, zero,
3079 $ work( iu+ldwrku ), ldwrku )
3090 CALL dgebrd( m, m, work( iu ), ldwrku, s,
3091 $ work( ie ), work( itauq ),
3092 $ work( itaup ), work( iwork ),
3093 $ lwork-iwork+1, ierr )
3094 CALL dlacpy(
'L', m, m, work( iu ), ldwrku,
3095 $ work( ir ), ldwrkr )
3101 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
3102 $ work( itaup ), work( iwork ),
3103 $ lwork-iwork+1, ierr )
3108 CALL dorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
3109 $ work( itauq ), work( iwork ),
3110 $ lwork-iwork+1, ierr )
3118 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
3119 $ work( iu ), ldwrku, work( ir ),
3120 $ ldwrkr, dum, 1, work( iwork ), info )
3126 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
3127 $ ldwrku, vt, ldvt, zero, a, lda )
3131 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
3135 CALL dlacpy(
'F', m, m, work( ir ), ldwrkr, a,
3148 CALL dgelqf( m, n, a, lda, work( itau ),
3149 $ work( iwork ), lwork-iwork+1, ierr )
3150 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3155 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3156 $ work( iwork ), lwork-iwork+1, ierr )
3164 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
3170 CALL dgebrd( m, m, a, lda, s, work( ie ),
3171 $ work( itauq ), work( itaup ),
3172 $ work( iwork ), lwork-iwork+1, ierr )
3178 CALL dormbr(
'P',
'L',
'T', m, n, m, a, lda,
3179 $ work( itaup ), vt, ldvt,
3180 $ work( iwork ), lwork-iwork+1, ierr )
3185 CALL dorgbr(
'Q', m, m, m, a, lda, work( itauq ),
3186 $ work( iwork ), lwork-iwork+1, ierr )
3194 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3195 $ ldvt, a, lda, dum, 1, work( iwork ),
3200 ELSE IF( wntuas )
THEN
3207 IF( lwork.GE.m*m+max( n + m, 4*m, bdspac ) )
THEN
3212 IF( lwork.GE.wrkbl+lda*m )
THEN
3223 itau = iu + ldwrku*m
3229 CALL dgelqf( m, n, a, lda, work( itau ),
3230 $ work( iwork ), lwork-iwork+1, ierr )
3231 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3236 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3237 $ work( iwork ), lwork-iwork+1, ierr )
3241 CALL dlacpy(
'L', m, m, a, lda, work( iu ),
3243 CALL dlaset(
'U', m-1, m-1, zero, zero,
3244 $ work( iu+ldwrku ), ldwrku )
3253 CALL dgebrd( m, m, work( iu ), ldwrku, s,
3254 $ work( ie ), work( itauq ),
3255 $ work( itaup ), work( iwork ),
3256 $ lwork-iwork+1, ierr )
3257 CALL dlacpy(
'L', m, m, work( iu ), ldwrku, u,
3263 CALL dorgbr(
'P', m, m, m, work( iu ), ldwrku,
3264 $ work( itaup ), work( iwork ),
3265 $ lwork-iwork+1, ierr )
3270 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3271 $ work( iwork ), lwork-iwork+1, ierr )
3279 CALL dbdsqr(
'U', m, m, m, 0, s, work( ie ),
3280 $ work( iu ), ldwrku, u, ldu, dum, 1,
3281 $ work( iwork ), info )
3287 CALL dgemm(
'N',
'N', m, n, m, one, work( iu ),
3288 $ ldwrku, vt, ldvt, zero, a, lda )
3292 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
3304 CALL dgelqf( m, n, a, lda, work( itau ),
3305 $ work( iwork ), lwork-iwork+1, ierr )
3306 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3311 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3312 $ work( iwork ), lwork-iwork+1, ierr )
3316 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
3317 CALL dlaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
3327 CALL dgebrd( m, m, u, ldu, s, work( ie ),
3328 $ work( itauq ), work( itaup ),
3329 $ work( iwork ), lwork-iwork+1, ierr )
3335 CALL dormbr(
'P',
'L',
'T', m, n, m, u, ldu,
3336 $ work( itaup ), vt, ldvt,
3337 $ work( iwork ), lwork-iwork+1, ierr )
3342 CALL dorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3343 $ work( iwork ), lwork-iwork+1, ierr )
3351 CALL dbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3352 $ ldvt, u, ldu, dum, 1, work( iwork ),
3376 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
3377 $ work( itaup ), work( iwork ), lwork-iwork+1,
3385 CALL dlacpy(
'L', m, m, a, lda, u, ldu )
3386 CALL dorgbr(
'Q', m, m, n, u, ldu, work( itauq ),
3387 $ work( iwork ), lwork-iwork+1, ierr )
3395 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
3400 CALL dorgbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3401 $ work( iwork ), lwork-iwork+1, ierr )
3409 CALL dorgbr(
'Q', m, m, n, a, lda, work( itauq ),
3410 $ work( iwork ), lwork-iwork+1, ierr )
3418 CALL dorgbr(
'P', m, n, m, a, lda, work( itaup ),
3419 $ work( iwork ), lwork-iwork+1, ierr )
3422 IF( wntuas .OR. wntuo )
3426 IF( wntvas .OR. wntvo )
3430 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3437 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3438 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
3439 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3446 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), a, lda,
3447 $ u, ldu, dum, 1, work( iwork ), info )
3455 CALL dbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3456 $ ldvt, a, lda, dum, 1, work( iwork ), info )
3466 IF( info.NE.0 )
THEN
3468 DO 50 i = 1, minmn - 1
3469 work( i+1 ) = work( i+ie-1 )
3473 DO 60 i = minmn - 1, 1, -1
3474 work( i+1 ) = work( i+ie-1 )
3481 IF( iscl.EQ.1 )
THEN
3482 IF( anrm.GT.bignum )
3483 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3485 IF( info.NE.0 .AND. anrm.GT.bignum )
3486 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn-1, 1, work( 2 ),
3488 IF( anrm.LT.smlnum )
3489 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3491 IF( info.NE.0 .AND. anrm.LT.smlnum )
3492 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1, work( 2 ),
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
double precision function dlamch(CMACH)
DLAMCH
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
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 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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR