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