211 SUBROUTINE sgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
212 $ work, lwork, info )
220 CHARACTER JOBU, JOBVT
221 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
224 REAL A( lda, * ), S( * ), U( ldu, * ),
225 $ vt( ldvt, * ), work( * )
232 parameter( zero = 0.0e0, one = 1.0e0 )
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_SGEQRF, LWORK_SORGQR_N, LWORK_SORGQR_M,
242 $ lwork_sgebrd, lwork_sorgbr_p, lwork_sorgbr_q,
243 $ lwork_sgelqf, lwork_sorglq_n, lwork_sorglq_m
244 REAL ANRM, BIGNUM, EPS, SMLNUM
258 EXTERNAL lsame, ilaenv, slamch, slange
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,
'SGESVD', jobu // jobvt, m, n, 0, 0 )
316 CALL sgeqrf( m, n, a, lda, dum(1), dum(1), -1, ierr )
319 CALL sorgqr( m, n, n, a, lda, dum(1), dum(1), -1, ierr )
320 lwork_sorgqr_n=dum(1)
321 CALL sorgqr( m, m, n, a, lda, dum(1), dum(1), -1, ierr )
322 lwork_sorgqr_m=dum(1)
324 CALL sgebrd( n, n, a, lda, s, dum(1), dum(1),
325 $ dum(1), dum(1), -1, ierr )
328 CALL sorgbr(
'P', n, n, n, a, lda, dum(1),
330 lwork_sorgbr_p=dum(1)
332 CALL sorgbr(
'Q', n, n, n, a, lda, dum(1),
334 lwork_sorgbr_q=dum(1)
336 IF( m.GE.mnthr )
THEN
341 maxwrk = n + lwork_sgeqrf
342 maxwrk = max( maxwrk, 3*n+lwork_sgebrd )
343 IF( wntvo .OR. wntvas )
344 $ maxwrk = max( maxwrk, 3*n+lwork_sorgbr_p )
345 maxwrk = max( maxwrk, bdspac )
346 minwrk = max( 4*n, bdspac )
347 ELSE IF( wntuo .AND. wntvn )
THEN
351 wrkbl = n + lwork_sgeqrf
352 wrkbl = max( wrkbl, n+lwork_sorgqr_n )
353 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
354 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_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_sgeqrf
364 wrkbl = max( wrkbl, n+lwork_sorgqr_n )
365 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
366 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
367 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_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_sgeqrf
376 wrkbl = max( wrkbl, n+lwork_sorgqr_n )
377 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
378 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
379 wrkbl = max( wrkbl, bdspac )
381 minwrk = max( 3*n+m, bdspac )
382 ELSE IF( wntus .AND. wntvo )
THEN
386 wrkbl = n + lwork_sgeqrf
387 wrkbl = max( wrkbl, n+lwork_sorgqr_n )
388 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
389 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
390 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_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_sgeqrf
400 wrkbl = max( wrkbl, n+lwork_sorgqr_n )
401 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
402 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
403 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_p )
404 wrkbl = max( wrkbl, bdspac )
406 minwrk = max( 3*n+m, bdspac )
407 ELSE IF( wntua .AND. wntvn )
THEN
411 wrkbl = n + lwork_sgeqrf
412 wrkbl = max( wrkbl, n+lwork_sorgqr_m )
413 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
414 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
415 wrkbl = max( wrkbl, bdspac )
417 minwrk = max( 3*n+m, bdspac )
418 ELSE IF( wntua .AND. wntvo )
THEN
422 wrkbl = n + lwork_sgeqrf
423 wrkbl = max( wrkbl, n+lwork_sorgqr_m )
424 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
425 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
426 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_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_sgeqrf
436 wrkbl = max( wrkbl, n+lwork_sorgqr_m )
437 wrkbl = max( wrkbl, 3*n+lwork_sgebrd )
438 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_q )
439 wrkbl = max( wrkbl, 3*n+lwork_sorgbr_p )
440 wrkbl = max( wrkbl, bdspac )
442 minwrk = max( 3*n+m, bdspac )
448 CALL sgebrd( m, n, a, lda, s, dum(1), dum(1),
449 $ dum(1), dum(1), -1, ierr )
451 maxwrk = 3*n + lwork_sgebrd
452 IF( wntus .OR. wntuo )
THEN
453 CALL sorgbr(
'Q', m, n, n, a, lda, dum(1),
455 lwork_sorgbr_q=dum(1)
456 maxwrk = max( maxwrk, 3*n+lwork_sorgbr_q )
459 CALL sorgbr(
'Q', m, m, n, a, lda, dum(1),
461 lwork_sorgbr_q=dum(1)
462 maxwrk = max( maxwrk, 3*n+lwork_sorgbr_q )
464 IF( .NOT.wntvn )
THEN
465 maxwrk = max( maxwrk, 3*n+lwork_sorgbr_p )
467 maxwrk = max( maxwrk, bdspac )
468 minwrk = max( 3*n+m, bdspac )
470 ELSE IF( minmn.GT.0 )
THEN
474 mnthr = ilaenv( 6,
'SGESVD', jobu // jobvt, m, n, 0, 0 )
477 CALL sgelqf( m, n, a, lda, dum(1), dum(1), -1, ierr )
480 CALL sorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
481 lwork_sorglq_n=dum(1)
482 CALL sorglq( m, n, m, a, lda, dum(1), dum(1), -1, ierr )
483 lwork_sorglq_m=dum(1)
485 CALL sgebrd( m, m, a, lda, s, dum(1), dum(1),
486 $ dum(1), dum(1), -1, ierr )
489 CALL sorgbr(
'P', m, m, m, a, n, dum(1),
491 lwork_sorgbr_p=dum(1)
493 CALL sorgbr(
'Q', m, m, m, a, n, dum(1),
495 lwork_sorgbr_q=dum(1)
496 IF( n.GE.mnthr )
THEN
501 maxwrk = m + lwork_sgelqf
502 maxwrk = max( maxwrk, 3*m+lwork_sgebrd )
503 IF( wntuo .OR. wntuas )
504 $ maxwrk = max( maxwrk, 3*m+lwork_sorgbr_q )
505 maxwrk = max( maxwrk, bdspac )
506 minwrk = max( 4*m, bdspac )
507 ELSE IF( wntvo .AND. wntun )
THEN
511 wrkbl = m + lwork_sgelqf
512 wrkbl = max( wrkbl, m+lwork_sorglq_m )
513 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
514 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_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_sgelqf
524 wrkbl = max( wrkbl, m+lwork_sorglq_m )
525 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
526 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
527 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_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_sgelqf
536 wrkbl = max( wrkbl, m+lwork_sorglq_m )
537 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
538 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
539 wrkbl = max( wrkbl, bdspac )
541 minwrk = max( 3*m+n, bdspac )
542 ELSE IF( wntvs .AND. wntuo )
THEN
546 wrkbl = m + lwork_sgelqf
547 wrkbl = max( wrkbl, m+lwork_sorglq_m )
548 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
549 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
550 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_q )
551 wrkbl = max( wrkbl, bdspac )
552 maxwrk = 2*m*m + wrkbl
553 minwrk = max( 3*m+n, bdspac )
554 maxwrk = max( maxwrk, minwrk )
555 ELSE IF( wntvs .AND. wntuas )
THEN
560 wrkbl = m + lwork_sgelqf
561 wrkbl = max( wrkbl, m+lwork_sorglq_m )
562 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
563 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
564 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_q )
565 wrkbl = max( wrkbl, bdspac )
567 minwrk = max( 3*m+n, bdspac )
568 ELSE IF( wntva .AND. wntun )
THEN
572 wrkbl = m + lwork_sgelqf
573 wrkbl = max( wrkbl, m+lwork_sorglq_n )
574 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
575 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
576 wrkbl = max( wrkbl, bdspac )
578 minwrk = max( 3*m+n, bdspac )
579 ELSE IF( wntva .AND. wntuo )
THEN
583 wrkbl = m + lwork_sgelqf
584 wrkbl = max( wrkbl, m+lwork_sorglq_n )
585 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
586 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
587 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_q )
588 wrkbl = max( wrkbl, bdspac )
589 maxwrk = 2*m*m + wrkbl
590 minwrk = max( 3*m+n, bdspac )
591 ELSE IF( wntva .AND. wntuas )
THEN
596 wrkbl = m + lwork_sgelqf
597 wrkbl = max( wrkbl, m+lwork_sorglq_n )
598 wrkbl = max( wrkbl, 3*m+lwork_sgebrd )
599 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_p )
600 wrkbl = max( wrkbl, 3*m+lwork_sorgbr_q )
601 wrkbl = max( wrkbl, bdspac )
603 minwrk = max( 3*m+n, bdspac )
609 CALL sgebrd( m, n, a, lda, s, dum(1), dum(1),
610 $ dum(1), dum(1), -1, ierr )
612 maxwrk = 3*m + lwork_sgebrd
613 IF( wntvs .OR. wntvo )
THEN
615 CALL sorgbr(
'P', m, n, m, a, n, dum(1),
617 lwork_sorgbr_p=dum(1)
618 maxwrk = max( maxwrk, 3*m+lwork_sorgbr_p )
621 CALL sorgbr(
'P', n, n, m, a, n, dum(1),
623 lwork_sorgbr_p=dum(1)
624 maxwrk = max( maxwrk, 3*m+lwork_sorgbr_p )
626 IF( .NOT.wntun )
THEN
627 maxwrk = max( maxwrk, 3*m+lwork_sorgbr_q )
629 maxwrk = max( maxwrk, bdspac )
630 minwrk = max( 3*m+n, bdspac )
633 maxwrk = max( maxwrk, minwrk )
636 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
642 CALL xerbla(
'SGESVD', -info )
644 ELSE IF( lquery )
THEN
650 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
657 smlnum = sqrt( slamch(
'S' ) ) / eps
658 bignum = one / smlnum
662 anrm = slange(
'M', m, n, a, lda, dum )
664 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
666 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
667 ELSE IF( anrm.GT.bignum )
THEN
669 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
678 IF( m.GE.mnthr )
THEN
691 CALL sgeqrf( m, n, a, lda, work( itau ), work( iwork ),
692 $ lwork-iwork+1, ierr )
696 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
705 CALL sgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
706 $ work( itaup ), work( iwork ), lwork-iwork+1,
709 IF( wntvo .OR. wntvas )
THEN
714 CALL sorgbr(
'P', n, n, n, a, lda, work( itaup ),
715 $ work( iwork ), lwork-iwork+1, ierr )
724 CALL sbdsqr(
'U', n, ncvt, 0, 0, s, work( ie ), a, lda,
725 $ dum, 1, dum, 1, work( iwork ), info )
730 $
CALL slacpy(
'F', n, n, a, lda, vt, ldvt )
732 ELSE IF( wntuo .AND. wntvn )
THEN
738 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
743 IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n )
THEN
749 ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n )
THEN
759 ldwrku = ( lwork-n*n-n ) / n
768 CALL sgeqrf( m, n, a, lda, work( itau ),
769 $ work( iwork ), lwork-iwork+1, ierr )
773 CALL slacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
774 CALL slaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
780 CALL sorgqr( m, n, n, a, lda, work( itau ),
781 $ work( iwork ), lwork-iwork+1, ierr )
790 CALL sgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
791 $ work( itauq ), work( itaup ),
792 $ work( iwork ), lwork-iwork+1, ierr )
797 CALL sorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
798 $ work( itauq ), work( iwork ),
799 $ lwork-iwork+1, ierr )
806 CALL sbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum, 1,
807 $ work( ir ), ldwrkr, dum, 1,
808 $ work( iwork ), info )
815 DO 10 i = 1, m, ldwrku
816 chunk = min( m-i+1, ldwrku )
817 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
818 $ lda, work( ir ), ldwrkr, zero,
819 $ work( iu ), ldwrku )
820 CALL slacpy(
'F', chunk, n, work( iu ), ldwrku,
836 CALL sgebrd( m, n, a, lda, s, work( ie ),
837 $ work( itauq ), work( itaup ),
838 $ work( iwork ), lwork-iwork+1, ierr )
843 CALL sorgbr(
'Q', m, n, n, a, lda, work( itauq ),
844 $ work( iwork ), lwork-iwork+1, ierr )
851 CALL sbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum, 1,
852 $ a, lda, dum, 1, work( iwork ), info )
856 ELSE IF( wntuo .AND. wntvas )
THEN
862 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
867 IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n )
THEN
873 ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n )
THEN
883 ldwrku = ( lwork-n*n-n ) / n
892 CALL sgeqrf( m, n, a, lda, work( itau ),
893 $ work( iwork ), lwork-iwork+1, ierr )
897 CALL slacpy(
'U', n, n, a, lda, vt, ldvt )
899 $
CALL slaset(
'L', n-1, n-1, zero, zero,
905 CALL sorgqr( m, n, n, a, lda, work( itau ),
906 $ work( iwork ), lwork-iwork+1, ierr )
915 CALL sgebrd( n, n, vt, ldvt, s, work( ie ),
916 $ work( itauq ), work( itaup ),
917 $ work( iwork ), lwork-iwork+1, ierr )
918 CALL slacpy(
'L', n, n, vt, ldvt, work( ir ), ldwrkr )
923 CALL sorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
924 $ work( itauq ), work( iwork ),
925 $ lwork-iwork+1, ierr )
930 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
931 $ work( iwork ), lwork-iwork+1, ierr )
939 CALL sbdsqr(
'U', n, n, n, 0, s, work( ie ), vt, ldvt,
940 $ work( ir ), ldwrkr, dum, 1,
941 $ work( iwork ), info )
948 DO 20 i = 1, m, ldwrku
949 chunk = min( m-i+1, ldwrku )
950 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
951 $ lda, work( ir ), ldwrkr, zero,
952 $ work( iu ), ldwrku )
953 CALL slacpy(
'F', chunk, n, work( iu ), ldwrku,
967 CALL sgeqrf( m, n, a, lda, work( itau ),
968 $ work( iwork ), lwork-iwork+1, ierr )
972 CALL slacpy(
'U', n, n, a, lda, vt, ldvt )
974 $
CALL slaset(
'L', n-1, n-1, zero, zero,
980 CALL sorgqr( m, n, n, a, lda, work( itau ),
981 $ work( iwork ), lwork-iwork+1, ierr )
990 CALL sgebrd( n, n, vt, ldvt, s, work( ie ),
991 $ work( itauq ), work( itaup ),
992 $ work( iwork ), lwork-iwork+1, ierr )
997 CALL sormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
998 $ work( itauq ), a, lda, work( iwork ),
999 $ lwork-iwork+1, ierr )
1004 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1005 $ work( iwork ), lwork-iwork+1, ierr )
1013 CALL sbdsqr(
'U', n, n, m, 0, s, work( ie ), vt, ldvt,
1014 $ a, lda, dum, 1, work( iwork ), info )
1018 ELSE IF( wntus )
THEN
1026 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1031 IF( lwork.GE.wrkbl+lda*n )
THEN
1042 itau = ir + ldwrkr*n
1048 CALL sgeqrf( m, n, a, lda, work( itau ),
1049 $ work( iwork ), lwork-iwork+1, ierr )
1053 CALL slacpy(
'U', n, n, a, lda, work( ir ),
1055 CALL slaset(
'L', n-1, n-1, zero, zero,
1056 $ work( ir+1 ), ldwrkr )
1061 CALL sorgqr( m, n, n, a, lda, work( itau ),
1062 $ work( iwork ), lwork-iwork+1, ierr )
1071 CALL sgebrd( n, n, work( ir ), ldwrkr, s,
1072 $ work( ie ), work( itauq ),
1073 $ work( itaup ), work( iwork ),
1074 $ lwork-iwork+1, ierr )
1079 CALL sorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1080 $ work( itauq ), work( iwork ),
1081 $ lwork-iwork+1, ierr )
1088 CALL sbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1089 $ 1, work( ir ), ldwrkr, dum, 1,
1090 $ work( iwork ), info )
1096 CALL sgemm(
'N',
'N', m, n, n, one, a, lda,
1097 $ work( ir ), ldwrkr, zero, u, ldu )
1109 CALL sgeqrf( m, n, a, lda, work( itau ),
1110 $ work( iwork ), lwork-iwork+1, ierr )
1111 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1116 CALL sorgqr( m, n, n, u, ldu, work( itau ),
1117 $ work( iwork ), lwork-iwork+1, ierr )
1125 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1131 CALL sgebrd( n, n, a, lda, s, work( ie ),
1132 $ work( itauq ), work( itaup ),
1133 $ work( iwork ), lwork-iwork+1, ierr )
1138 CALL sormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1139 $ work( itauq ), u, ldu, work( iwork ),
1140 $ lwork-iwork+1, ierr )
1147 CALL sbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1148 $ 1, u, ldu, dum, 1, work( iwork ),
1153 ELSE IF( wntvo )
THEN
1159 IF( lwork.GE.2*n*n+max( 4*n, bdspac ) )
THEN
1164 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1171 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1186 itau = ir + ldwrkr*n
1192 CALL sgeqrf( m, n, a, lda, work( itau ),
1193 $ work( iwork ), lwork-iwork+1, ierr )
1197 CALL slacpy(
'U', n, n, a, lda, work( iu ),
1199 CALL slaset(
'L', n-1, n-1, zero, zero,
1200 $ work( iu+1 ), ldwrku )
1205 CALL sorgqr( m, n, n, a, lda, work( itau ),
1206 $ work( iwork ), lwork-iwork+1, ierr )
1217 CALL sgebrd( n, n, work( iu ), ldwrku, s,
1218 $ work( ie ), work( itauq ),
1219 $ work( itaup ), work( iwork ),
1220 $ lwork-iwork+1, ierr )
1221 CALL slacpy(
'U', n, n, work( iu ), ldwrku,
1222 $ work( ir ), ldwrkr )
1227 CALL sorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1228 $ work( itauq ), work( iwork ),
1229 $ lwork-iwork+1, ierr )
1235 CALL sorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1236 $ work( itaup ), work( iwork ),
1237 $ lwork-iwork+1, ierr )
1245 CALL sbdsqr(
'U', n, n, n, 0, s, work( ie ),
1246 $ work( ir ), ldwrkr, work( iu ),
1247 $ ldwrku, dum, 1, work( iwork ), info )
1253 CALL sgemm(
'N',
'N', m, n, n, one, a, lda,
1254 $ work( iu ), ldwrku, zero, u, ldu )
1259 CALL slacpy(
'F', n, n, work( ir ), ldwrkr, a,
1272 CALL sgeqrf( m, n, a, lda, work( itau ),
1273 $ work( iwork ), lwork-iwork+1, ierr )
1274 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1279 CALL sorgqr( m, n, n, u, ldu, work( itau ),
1280 $ work( iwork ), lwork-iwork+1, ierr )
1288 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1294 CALL sgebrd( n, n, a, lda, s, work( ie ),
1295 $ work( itauq ), work( itaup ),
1296 $ work( iwork ), lwork-iwork+1, ierr )
1301 CALL sormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1302 $ work( itauq ), u, ldu, work( iwork ),
1303 $ lwork-iwork+1, ierr )
1308 CALL sorgbr(
'P', n, n, n, a, lda, work( itaup ),
1309 $ work( iwork ), lwork-iwork+1, ierr )
1317 CALL sbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1318 $ lda, u, ldu, dum, 1, work( iwork ),
1323 ELSE IF( wntvas )
THEN
1330 IF( lwork.GE.n*n+max( 4*n, bdspac ) )
THEN
1335 IF( lwork.GE.wrkbl+lda*n )
THEN
1346 itau = iu + ldwrku*n
1352 CALL sgeqrf( m, n, a, lda, work( itau ),
1353 $ work( iwork ), lwork-iwork+1, ierr )
1357 CALL slacpy(
'U', n, n, a, lda, work( iu ),
1359 CALL slaset(
'L', n-1, n-1, zero, zero,
1360 $ work( iu+1 ), ldwrku )
1365 CALL sorgqr( m, n, n, a, lda, work( itau ),
1366 $ work( iwork ), lwork-iwork+1, ierr )
1375 CALL sgebrd( n, n, work( iu ), ldwrku, s,
1376 $ work( ie ), work( itauq ),
1377 $ work( itaup ), work( iwork ),
1378 $ lwork-iwork+1, ierr )
1379 CALL slacpy(
'U', n, n, work( iu ), ldwrku, vt,
1385 CALL sorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1386 $ work( itauq ), work( iwork ),
1387 $ lwork-iwork+1, ierr )
1393 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1394 $ work( iwork ), lwork-iwork+1, ierr )
1402 CALL sbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1403 $ ldvt, work( iu ), ldwrku, dum, 1,
1404 $ work( iwork ), info )
1410 CALL sgemm(
'N',
'N', m, n, n, one, a, lda,
1411 $ work( iu ), ldwrku, zero, u, ldu )
1423 CALL sgeqrf( m, n, a, lda, work( itau ),
1424 $ work( iwork ), lwork-iwork+1, ierr )
1425 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1430 CALL sorgqr( m, n, n, u, ldu, work( itau ),
1431 $ work( iwork ), lwork-iwork+1, ierr )
1435 CALL slacpy(
'U', n, n, a, lda, vt, ldvt )
1437 $
CALL slaset(
'L', n-1, n-1, zero, zero,
1438 $ vt( 2, 1 ), ldvt )
1447 CALL sgebrd( n, n, vt, ldvt, s, work( ie ),
1448 $ work( itauq ), work( itaup ),
1449 $ work( iwork ), lwork-iwork+1, ierr )
1455 CALL sormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1456 $ work( itauq ), u, ldu, work( iwork ),
1457 $ lwork-iwork+1, ierr )
1462 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1463 $ work( iwork ), lwork-iwork+1, ierr )
1471 CALL sbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1472 $ ldvt, u, ldu, dum, 1, work( iwork ),
1479 ELSE IF( wntua )
THEN
1487 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1492 IF( lwork.GE.wrkbl+lda*n )
THEN
1503 itau = ir + ldwrkr*n
1509 CALL sgeqrf( m, n, a, lda, work( itau ),
1510 $ work( iwork ), lwork-iwork+1, ierr )
1511 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1515 CALL slacpy(
'U', n, n, a, lda, work( ir ),
1517 CALL slaset(
'L', n-1, n-1, zero, zero,
1518 $ work( ir+1 ), ldwrkr )
1523 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1524 $ work( iwork ), lwork-iwork+1, ierr )
1533 CALL sgebrd( n, n, work( ir ), ldwrkr, s,
1534 $ work( ie ), work( itauq ),
1535 $ work( itaup ), work( iwork ),
1536 $ lwork-iwork+1, ierr )
1541 CALL sorgbr(
'Q', n, n, n, work( ir ), ldwrkr,
1542 $ work( itauq ), work( iwork ),
1543 $ lwork-iwork+1, ierr )
1550 CALL sbdsqr(
'U', n, 0, n, 0, s, work( ie ), dum,
1551 $ 1, work( ir ), ldwrkr, dum, 1,
1552 $ work( iwork ), info )
1558 CALL sgemm(
'N',
'N', m, n, n, one, u, ldu,
1559 $ work( ir ), ldwrkr, zero, a, lda )
1563 CALL slacpy(
'F', m, n, a, lda, u, ldu )
1575 CALL sgeqrf( m, n, a, lda, work( itau ),
1576 $ work( iwork ), lwork-iwork+1, ierr )
1577 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1582 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1583 $ work( iwork ), lwork-iwork+1, ierr )
1591 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1597 CALL sgebrd( n, n, a, lda, s, work( ie ),
1598 $ work( itauq ), work( itaup ),
1599 $ work( iwork ), lwork-iwork+1, ierr )
1605 CALL sormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1606 $ work( itauq ), u, ldu, work( iwork ),
1607 $ lwork-iwork+1, ierr )
1614 CALL sbdsqr(
'U', n, 0, m, 0, s, work( ie ), dum,
1615 $ 1, u, ldu, dum, 1, work( iwork ),
1620 ELSE IF( wntvo )
THEN
1626 IF( lwork.GE.2*n*n+max( n+m, 4*n, bdspac ) )
THEN
1631 IF( lwork.GE.wrkbl+2*lda*n )
THEN
1638 ELSE IF( lwork.GE.wrkbl+( lda+n )*n )
THEN
1653 itau = ir + ldwrkr*n
1659 CALL sgeqrf( m, n, a, lda, work( itau ),
1660 $ work( iwork ), lwork-iwork+1, ierr )
1661 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1666 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1667 $ work( iwork ), lwork-iwork+1, ierr )
1671 CALL slacpy(
'U', n, n, a, lda, work( iu ),
1673 CALL slaset(
'L', n-1, n-1, zero, zero,
1674 $ work( iu+1 ), ldwrku )
1685 CALL sgebrd( n, n, work( iu ), ldwrku, s,
1686 $ work( ie ), work( itauq ),
1687 $ work( itaup ), work( iwork ),
1688 $ lwork-iwork+1, ierr )
1689 CALL slacpy(
'U', n, n, work( iu ), ldwrku,
1690 $ work( ir ), ldwrkr )
1695 CALL sorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1696 $ work( itauq ), work( iwork ),
1697 $ lwork-iwork+1, ierr )
1703 CALL sorgbr(
'P', n, n, n, work( ir ), ldwrkr,
1704 $ work( itaup ), work( iwork ),
1705 $ lwork-iwork+1, ierr )
1713 CALL sbdsqr(
'U', n, n, n, 0, s, work( ie ),
1714 $ work( ir ), ldwrkr, work( iu ),
1715 $ ldwrku, dum, 1, work( iwork ), info )
1721 CALL sgemm(
'N',
'N', m, n, n, one, u, ldu,
1722 $ work( iu ), ldwrku, zero, a, lda )
1726 CALL slacpy(
'F', m, n, a, lda, u, ldu )
1730 CALL slacpy(
'F', n, n, work( ir ), ldwrkr, a,
1743 CALL sgeqrf( m, n, a, lda, work( itau ),
1744 $ work( iwork ), lwork-iwork+1, ierr )
1745 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1750 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1751 $ work( iwork ), lwork-iwork+1, ierr )
1759 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
1765 CALL sgebrd( n, n, a, lda, s, work( ie ),
1766 $ work( itauq ), work( itaup ),
1767 $ work( iwork ), lwork-iwork+1, ierr )
1773 CALL sormbr(
'Q',
'R',
'N', m, n, n, a, lda,
1774 $ work( itauq ), u, ldu, work( iwork ),
1775 $ lwork-iwork+1, ierr )
1780 CALL sorgbr(
'P', n, n, n, a, lda, work( itaup ),
1781 $ work( iwork ), lwork-iwork+1, ierr )
1789 CALL sbdsqr(
'U', n, n, m, 0, s, work( ie ), a,
1790 $ lda, u, ldu, dum, 1, work( iwork ),
1795 ELSE IF( wntvas )
THEN
1802 IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) )
THEN
1807 IF( lwork.GE.wrkbl+lda*n )
THEN
1818 itau = iu + ldwrku*n
1824 CALL sgeqrf( m, n, a, lda, work( itau ),
1825 $ work( iwork ), lwork-iwork+1, ierr )
1826 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1831 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1832 $ work( iwork ), lwork-iwork+1, ierr )
1836 CALL slacpy(
'U', n, n, a, lda, work( iu ),
1838 CALL slaset(
'L', n-1, n-1, zero, zero,
1839 $ work( iu+1 ), ldwrku )
1848 CALL sgebrd( n, n, work( iu ), ldwrku, s,
1849 $ work( ie ), work( itauq ),
1850 $ work( itaup ), work( iwork ),
1851 $ lwork-iwork+1, ierr )
1852 CALL slacpy(
'U', n, n, work( iu ), ldwrku, vt,
1858 CALL sorgbr(
'Q', n, n, n, work( iu ), ldwrku,
1859 $ work( itauq ), work( iwork ),
1860 $ lwork-iwork+1, ierr )
1866 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1867 $ work( iwork ), lwork-iwork+1, ierr )
1875 CALL sbdsqr(
'U', n, n, n, 0, s, work( ie ), vt,
1876 $ ldvt, work( iu ), ldwrku, dum, 1,
1877 $ work( iwork ), info )
1883 CALL sgemm(
'N',
'N', m, n, n, one, u, ldu,
1884 $ work( iu ), ldwrku, zero, a, lda )
1888 CALL slacpy(
'F', m, n, a, lda, u, ldu )
1900 CALL sgeqrf( m, n, a, lda, work( itau ),
1901 $ work( iwork ), lwork-iwork+1, ierr )
1902 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1907 CALL sorgqr( m, m, n, u, ldu, work( itau ),
1908 $ work( iwork ), lwork-iwork+1, ierr )
1912 CALL slacpy(
'U', n, n, a, lda, vt, ldvt )
1914 $
CALL slaset(
'L', n-1, n-1, zero, zero,
1915 $ vt( 2, 1 ), ldvt )
1924 CALL sgebrd( n, n, vt, ldvt, s, work( ie ),
1925 $ work( itauq ), work( itaup ),
1926 $ work( iwork ), lwork-iwork+1, ierr )
1932 CALL sormbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1933 $ work( itauq ), u, ldu, work( iwork ),
1934 $ lwork-iwork+1, ierr )
1939 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1940 $ work( iwork ), lwork-iwork+1, ierr )
1948 CALL sbdsqr(
'U', n, n, m, 0, s, work( ie ), vt,
1949 $ ldvt, u, ldu, dum, 1, work( iwork ),
1973 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1974 $ work( itaup ), work( iwork ), lwork-iwork+1,
1982 CALL slacpy(
'L', m, n, a, lda, u, ldu )
1987 CALL sorgbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
1988 $ work( iwork ), lwork-iwork+1, ierr )
1996 CALL slacpy(
'U', n, n, a, lda, vt, ldvt )
1997 CALL sorgbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1998 $ work( iwork ), lwork-iwork+1, ierr )
2006 CALL sorgbr(
'Q', m, n, n, a, lda, work( itauq ),
2007 $ work( iwork ), lwork-iwork+1, ierr )
2015 CALL sorgbr(
'P', n, n, n, a, lda, work( itaup ),
2016 $ work( iwork ), lwork-iwork+1, ierr )
2019 IF( wntuas .OR. wntuo )
2023 IF( wntvas .OR. wntvo )
2027 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
2034 CALL sbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2035 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
2036 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
2043 CALL sbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), a, lda,
2044 $ u, ldu, dum, 1, work( iwork ), info )
2052 CALL sbdsqr(
'U', n, ncvt, nru, 0, s, work( ie ), vt,
2053 $ ldvt, a, lda, dum, 1, work( iwork ), info )
2064 IF( n.GE.mnthr )
THEN
2077 CALL sgelqf( m, n, a, lda, work( itau ), work( iwork ),
2078 $ lwork-iwork+1, ierr )
2082 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
2091 CALL sgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
2092 $ work( itaup ), work( iwork ), lwork-iwork+1,
2094 IF( wntuo .OR. wntuas )
THEN
2099 CALL sorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2100 $ work( iwork ), lwork-iwork+1, ierr )
2104 IF( wntuo .OR. wntuas )
2111 CALL sbdsqr(
'U', m, 0, nru, 0, s, work( ie ), dum, 1, a,
2112 $ lda, dum, 1, work( iwork ), info )
2117 $
CALL slacpy(
'F', m, m, a, lda, u, ldu )
2119 ELSE IF( wntvo .AND. wntun )
THEN
2125 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2130 IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m )
THEN
2137 ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m )
THEN
2149 chunk = ( lwork-m*m-m ) / m
2152 itau = ir + ldwrkr*m
2158 CALL sgelqf( m, n, a, lda, work( itau ),
2159 $ work( iwork ), lwork-iwork+1, ierr )
2163 CALL slacpy(
'L', m, m, a, lda, work( ir ), ldwrkr )
2164 CALL slaset(
'U', m-1, m-1, zero, zero,
2165 $ work( ir+ldwrkr ), ldwrkr )
2170 CALL sorglq( m, n, m, a, lda, work( itau ),
2171 $ work( iwork ), lwork-iwork+1, ierr )
2180 CALL sgebrd( m, m, work( ir ), ldwrkr, s, work( ie ),
2181 $ work( itauq ), work( itaup ),
2182 $ work( iwork ), lwork-iwork+1, ierr )
2187 CALL sorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2188 $ work( itaup ), work( iwork ),
2189 $ lwork-iwork+1, ierr )
2196 CALL sbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2197 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2198 $ work( iwork ), info )
2205 DO 30 i = 1, n, chunk
2206 blk = min( n-i+1, chunk )
2207 CALL sgemm(
'N',
'N', m, blk, m, one, work( ir ),
2208 $ ldwrkr, a( 1, i ), lda, zero,
2209 $ work( iu ), ldwrku )
2210 CALL slacpy(
'F', m, blk, work( iu ), ldwrku,
2226 CALL sgebrd( m, n, a, lda, s, work( ie ),
2227 $ work( itauq ), work( itaup ),
2228 $ work( iwork ), lwork-iwork+1, ierr )
2233 CALL sorgbr(
'P', m, n, m, a, lda, work( itaup ),
2234 $ work( iwork ), lwork-iwork+1, ierr )
2241 CALL sbdsqr(
'L', m, n, 0, 0, s, work( ie ), a, lda,
2242 $ dum, 1, dum, 1, work( iwork ), info )
2246 ELSE IF( wntvo .AND. wntuas )
THEN
2252 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2257 IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m )
THEN
2264 ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m )
THEN
2276 chunk = ( lwork-m*m-m ) / m
2279 itau = ir + ldwrkr*m
2285 CALL sgelqf( m, n, a, lda, work( itau ),
2286 $ work( iwork ), lwork-iwork+1, ierr )
2290 CALL slacpy(
'L', m, m, a, lda, u, ldu )
2291 CALL slaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2297 CALL sorglq( m, n, m, a, lda, work( itau ),
2298 $ work( iwork ), lwork-iwork+1, ierr )
2307 CALL sgebrd( m, m, u, ldu, s, work( ie ),
2308 $ work( itauq ), work( itaup ),
2309 $ work( iwork ), lwork-iwork+1, ierr )
2310 CALL slacpy(
'U', m, m, u, ldu, work( ir ), ldwrkr )
2315 CALL sorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2316 $ work( itaup ), work( iwork ),
2317 $ lwork-iwork+1, ierr )
2322 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2323 $ work( iwork ), lwork-iwork+1, ierr )
2331 CALL sbdsqr(
'U', m, m, m, 0, s, work( ie ),
2332 $ work( ir ), ldwrkr, u, ldu, dum, 1,
2333 $ work( iwork ), info )
2340 DO 40 i = 1, n, chunk
2341 blk = min( n-i+1, chunk )
2342 CALL sgemm(
'N',
'N', m, blk, m, one, work( ir ),
2343 $ ldwrkr, a( 1, i ), lda, zero,
2344 $ work( iu ), ldwrku )
2345 CALL slacpy(
'F', m, blk, work( iu ), ldwrku,
2359 CALL sgelqf( m, n, a, lda, work( itau ),
2360 $ work( iwork ), lwork-iwork+1, ierr )
2364 CALL slacpy(
'L', m, m, a, lda, u, ldu )
2365 CALL slaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2371 CALL sorglq( m, n, m, a, lda, work( itau ),
2372 $ work( iwork ), lwork-iwork+1, ierr )
2381 CALL sgebrd( m, m, u, ldu, s, work( ie ),
2382 $ work( itauq ), work( itaup ),
2383 $ work( iwork ), lwork-iwork+1, ierr )
2388 CALL sormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2389 $ work( itaup ), a, lda, work( iwork ),
2390 $ lwork-iwork+1, ierr )
2395 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2396 $ work( iwork ), lwork-iwork+1, ierr )
2404 CALL sbdsqr(
'U', m, n, m, 0, s, work( ie ), a, lda,
2405 $ u, ldu, dum, 1, work( iwork ), info )
2409 ELSE IF( wntvs )
THEN
2417 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2422 IF( lwork.GE.wrkbl+lda*m )
THEN
2433 itau = ir + ldwrkr*m
2439 CALL sgelqf( m, n, a, lda, work( itau ),
2440 $ work( iwork ), lwork-iwork+1, ierr )
2444 CALL slacpy(
'L', m, m, a, lda, work( ir ),
2446 CALL slaset(
'U', m-1, m-1, zero, zero,
2447 $ work( ir+ldwrkr ), ldwrkr )
2452 CALL sorglq( m, n, m, a, lda, work( itau ),
2453 $ work( iwork ), lwork-iwork+1, ierr )
2462 CALL sgebrd( m, m, work( ir ), ldwrkr, s,
2463 $ work( ie ), work( itauq ),
2464 $ work( itaup ), work( iwork ),
2465 $ lwork-iwork+1, ierr )
2471 CALL sorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2472 $ work( itaup ), work( iwork ),
2473 $ lwork-iwork+1, ierr )
2480 CALL sbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2481 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2482 $ work( iwork ), info )
2488 CALL sgemm(
'N',
'N', m, n, m, one, work( ir ),
2489 $ ldwrkr, a, lda, zero, vt, ldvt )
2501 CALL sgelqf( m, n, a, lda, work( itau ),
2502 $ work( iwork ), lwork-iwork+1, ierr )
2506 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
2511 CALL sorglq( m, n, m, vt, ldvt, work( itau ),
2512 $ work( iwork ), lwork-iwork+1, ierr )
2520 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2526 CALL sgebrd( m, m, a, lda, s, work( ie ),
2527 $ work( itauq ), work( itaup ),
2528 $ work( iwork ), lwork-iwork+1, ierr )
2533 CALL sormbr(
'P',
'L',
'T', m, n, m, a, lda,
2534 $ work( itaup ), vt, ldvt,
2535 $ work( iwork ), lwork-iwork+1, ierr )
2542 CALL sbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
2543 $ ldvt, dum, 1, dum, 1, work( iwork ),
2548 ELSE IF( wntuo )
THEN
2554 IF( lwork.GE.2*m*m+max( 4*m, bdspac ) )
THEN
2559 IF( lwork.GE.wrkbl+2*lda*m )
THEN
2566 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
2581 itau = ir + ldwrkr*m
2587 CALL sgelqf( m, n, a, lda, work( itau ),
2588 $ work( iwork ), lwork-iwork+1, ierr )
2592 CALL slacpy(
'L', m, m, a, lda, work( iu ),
2594 CALL slaset(
'U', m-1, m-1, zero, zero,
2595 $ work( iu+ldwrku ), ldwrku )
2600 CALL sorglq( m, n, m, a, lda, work( itau ),
2601 $ work( iwork ), lwork-iwork+1, ierr )
2612 CALL sgebrd( m, m, work( iu ), ldwrku, s,
2613 $ work( ie ), work( itauq ),
2614 $ work( itaup ), work( iwork ),
2615 $ lwork-iwork+1, ierr )
2616 CALL slacpy(
'L', m, m, work( iu ), ldwrku,
2617 $ work( ir ), ldwrkr )
2623 CALL sorgbr(
'P', m, m, m, work( iu ), ldwrku,
2624 $ work( itaup ), work( iwork ),
2625 $ lwork-iwork+1, ierr )
2630 CALL sorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
2631 $ work( itauq ), work( iwork ),
2632 $ lwork-iwork+1, ierr )
2640 CALL sbdsqr(
'U', m, m, m, 0, s, work( ie ),
2641 $ work( iu ), ldwrku, work( ir ),
2642 $ ldwrkr, dum, 1, work( iwork ), info )
2648 CALL sgemm(
'N',
'N', m, n, m, one, work( iu ),
2649 $ ldwrku, a, lda, zero, vt, ldvt )
2654 CALL slacpy(
'F', m, m, work( ir ), ldwrkr, a,
2667 CALL sgelqf( m, n, a, lda, work( itau ),
2668 $ work( iwork ), lwork-iwork+1, ierr )
2669 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
2674 CALL sorglq( m, n, m, vt, ldvt, work( itau ),
2675 $ work( iwork ), lwork-iwork+1, ierr )
2683 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2689 CALL sgebrd( m, m, a, lda, s, work( ie ),
2690 $ work( itauq ), work( itaup ),
2691 $ work( iwork ), lwork-iwork+1, ierr )
2696 CALL sormbr(
'P',
'L',
'T', m, n, m, a, lda,
2697 $ work( itaup ), vt, ldvt,
2698 $ work( iwork ), lwork-iwork+1, ierr )
2703 CALL sorgbr(
'Q', m, m, m, a, lda, work( itauq ),
2704 $ work( iwork ), lwork-iwork+1, ierr )
2712 CALL sbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2713 $ ldvt, a, lda, dum, 1, work( iwork ),
2718 ELSE IF( wntuas )
THEN
2725 IF( lwork.GE.m*m+max( 4*m, bdspac ) )
THEN
2730 IF( lwork.GE.wrkbl+lda*m )
THEN
2741 itau = iu + ldwrku*m
2747 CALL sgelqf( m, n, a, lda, work( itau ),
2748 $ work( iwork ), lwork-iwork+1, ierr )
2752 CALL slacpy(
'L', m, m, a, lda, work( iu ),
2754 CALL slaset(
'U', m-1, m-1, zero, zero,
2755 $ work( iu+ldwrku ), ldwrku )
2760 CALL sorglq( m, n, m, a, lda, work( itau ),
2761 $ work( iwork ), lwork-iwork+1, ierr )
2770 CALL sgebrd( m, m, work( iu ), ldwrku, s,
2771 $ work( ie ), work( itauq ),
2772 $ work( itaup ), work( iwork ),
2773 $ lwork-iwork+1, ierr )
2774 CALL slacpy(
'L', m, m, work( iu ), ldwrku, u,
2781 CALL sorgbr(
'P', m, m, m, work( iu ), ldwrku,
2782 $ work( itaup ), work( iwork ),
2783 $ lwork-iwork+1, ierr )
2788 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2789 $ work( iwork ), lwork-iwork+1, ierr )
2797 CALL sbdsqr(
'U', m, m, m, 0, s, work( ie ),
2798 $ work( iu ), ldwrku, u, ldu, dum, 1,
2799 $ work( iwork ), info )
2805 CALL sgemm(
'N',
'N', m, n, m, one, work( iu ),
2806 $ ldwrku, a, lda, zero, vt, ldvt )
2818 CALL sgelqf( m, n, a, lda, work( itau ),
2819 $ work( iwork ), lwork-iwork+1, ierr )
2820 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
2825 CALL sorglq( m, n, m, vt, ldvt, work( itau ),
2826 $ work( iwork ), lwork-iwork+1, ierr )
2830 CALL slacpy(
'L', m, m, a, lda, u, ldu )
2831 CALL slaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
2841 CALL sgebrd( m, m, u, ldu, s, work( ie ),
2842 $ work( itauq ), work( itaup ),
2843 $ work( iwork ), lwork-iwork+1, ierr )
2849 CALL sormbr(
'P',
'L',
'T', m, n, m, u, ldu,
2850 $ work( itaup ), vt, ldvt,
2851 $ work( iwork ), lwork-iwork+1, ierr )
2856 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
2857 $ work( iwork ), lwork-iwork+1, ierr )
2865 CALL sbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
2866 $ ldvt, u, ldu, dum, 1, work( iwork ),
2873 ELSE IF( wntva )
THEN
2881 IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) )
THEN
2886 IF( lwork.GE.wrkbl+lda*m )
THEN
2897 itau = ir + ldwrkr*m
2903 CALL sgelqf( m, n, a, lda, work( itau ),
2904 $ work( iwork ), lwork-iwork+1, ierr )
2905 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
2909 CALL slacpy(
'L', m, m, a, lda, work( ir ),
2911 CALL slaset(
'U', m-1, m-1, zero, zero,
2912 $ work( ir+ldwrkr ), ldwrkr )
2917 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
2918 $ work( iwork ), lwork-iwork+1, ierr )
2927 CALL sgebrd( m, m, work( ir ), ldwrkr, s,
2928 $ work( ie ), work( itauq ),
2929 $ work( itaup ), work( iwork ),
2930 $ lwork-iwork+1, ierr )
2936 CALL sorgbr(
'P', m, m, m, work( ir ), ldwrkr,
2937 $ work( itaup ), work( iwork ),
2938 $ lwork-iwork+1, ierr )
2945 CALL sbdsqr(
'U', m, m, 0, 0, s, work( ie ),
2946 $ work( ir ), ldwrkr, dum, 1, dum, 1,
2947 $ work( iwork ), info )
2953 CALL sgemm(
'N',
'N', m, n, m, one, work( ir ),
2954 $ ldwrkr, vt, ldvt, zero, a, lda )
2958 CALL slacpy(
'F', m, n, a, lda, vt, ldvt )
2970 CALL sgelqf( m, n, a, lda, work( itau ),
2971 $ work( iwork ), lwork-iwork+1, ierr )
2972 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
2977 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
2978 $ work( iwork ), lwork-iwork+1, ierr )
2986 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
2992 CALL sgebrd( m, m, a, lda, s, work( ie ),
2993 $ work( itauq ), work( itaup ),
2994 $ work( iwork ), lwork-iwork+1, ierr )
3000 CALL sormbr(
'P',
'L',
'T', m, n, m, a, lda,
3001 $ work( itaup ), vt, ldvt,
3002 $ work( iwork ), lwork-iwork+1, ierr )
3009 CALL sbdsqr(
'U', m, n, 0, 0, s, work( ie ), vt,
3010 $ ldvt, dum, 1, dum, 1, work( iwork ),
3015 ELSE IF( wntuo )
THEN
3021 IF( lwork.GE.2*m*m+max( n+m, 4*m, bdspac ) )
THEN
3026 IF( lwork.GE.wrkbl+2*lda*m )
THEN
3033 ELSE IF( lwork.GE.wrkbl+( lda+m )*m )
THEN
3048 itau = ir + ldwrkr*m
3054 CALL sgelqf( m, n, a, lda, work( itau ),
3055 $ work( iwork ), lwork-iwork+1, ierr )
3056 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
3061 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
3062 $ work( iwork ), lwork-iwork+1, ierr )
3066 CALL slacpy(
'L', m, m, a, lda, work( iu ),
3068 CALL slaset(
'U', m-1, m-1, zero, zero,
3069 $ work( iu+ldwrku ), ldwrku )
3080 CALL sgebrd( m, m, work( iu ), ldwrku, s,
3081 $ work( ie ), work( itauq ),
3082 $ work( itaup ), work( iwork ),
3083 $ lwork-iwork+1, ierr )
3084 CALL slacpy(
'L', m, m, work( iu ), ldwrku,
3085 $ work( ir ), ldwrkr )
3091 CALL sorgbr(
'P', m, m, m, work( iu ), ldwrku,
3092 $ work( itaup ), work( iwork ),
3093 $ lwork-iwork+1, ierr )
3098 CALL sorgbr(
'Q', m, m, m, work( ir ), ldwrkr,
3099 $ work( itauq ), work( iwork ),
3100 $ lwork-iwork+1, ierr )
3108 CALL sbdsqr(
'U', m, m, m, 0, s, work( ie ),
3109 $ work( iu ), ldwrku, work( ir ),
3110 $ ldwrkr, dum, 1, work( iwork ), info )
3116 CALL sgemm(
'N',
'N', m, n, m, one, work( iu ),
3117 $ ldwrku, vt, ldvt, zero, a, lda )
3121 CALL slacpy(
'F', m, n, a, lda, vt, ldvt )
3125 CALL slacpy(
'F', m, m, work( ir ), ldwrkr, a,
3138 CALL sgelqf( m, n, a, lda, work( itau ),
3139 $ work( iwork ), lwork-iwork+1, ierr )
3140 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
3145 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
3146 $ work( iwork ), lwork-iwork+1, ierr )
3154 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
3160 CALL sgebrd( m, m, a, lda, s, work( ie ),
3161 $ work( itauq ), work( itaup ),
3162 $ work( iwork ), lwork-iwork+1, ierr )
3168 CALL sormbr(
'P',
'L',
'T', m, n, m, a, lda,
3169 $ work( itaup ), vt, ldvt,
3170 $ work( iwork ), lwork-iwork+1, ierr )
3175 CALL sorgbr(
'Q', m, m, m, a, lda, work( itauq ),
3176 $ work( iwork ), lwork-iwork+1, ierr )
3184 CALL sbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3185 $ ldvt, a, lda, dum, 1, work( iwork ),
3190 ELSE IF( wntuas )
THEN
3197 IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) )
THEN
3202 IF( lwork.GE.wrkbl+lda*m )
THEN
3213 itau = iu + ldwrku*m
3219 CALL sgelqf( m, n, a, lda, work( itau ),
3220 $ work( iwork ), lwork-iwork+1, ierr )
3221 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
3226 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
3227 $ work( iwork ), lwork-iwork+1, ierr )
3231 CALL slacpy(
'L', m, m, a, lda, work( iu ),
3233 CALL slaset(
'U', m-1, m-1, zero, zero,
3234 $ work( iu+ldwrku ), ldwrku )
3243 CALL sgebrd( m, m, work( iu ), ldwrku, s,
3244 $ work( ie ), work( itauq ),
3245 $ work( itaup ), work( iwork ),
3246 $ lwork-iwork+1, ierr )
3247 CALL slacpy(
'L', m, m, work( iu ), ldwrku, u,
3253 CALL sorgbr(
'P', m, m, m, work( iu ), ldwrku,
3254 $ work( itaup ), work( iwork ),
3255 $ lwork-iwork+1, ierr )
3260 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3261 $ work( iwork ), lwork-iwork+1, ierr )
3269 CALL sbdsqr(
'U', m, m, m, 0, s, work( ie ),
3270 $ work( iu ), ldwrku, u, ldu, dum, 1,
3271 $ work( iwork ), info )
3277 CALL sgemm(
'N',
'N', m, n, m, one, work( iu ),
3278 $ ldwrku, vt, ldvt, zero, a, lda )
3282 CALL slacpy(
'F', m, n, a, lda, vt, ldvt )
3294 CALL sgelqf( m, n, a, lda, work( itau ),
3295 $ work( iwork ), lwork-iwork+1, ierr )
3296 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
3301 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
3302 $ work( iwork ), lwork-iwork+1, ierr )
3306 CALL slacpy(
'L', m, m, a, lda, u, ldu )
3307 CALL slaset(
'U', m-1, m-1, zero, zero, u( 1, 2 ),
3317 CALL sgebrd( m, m, u, ldu, s, work( ie ),
3318 $ work( itauq ), work( itaup ),
3319 $ work( iwork ), lwork-iwork+1, ierr )
3325 CALL sormbr(
'P',
'L',
'T', m, n, m, u, ldu,
3326 $ work( itaup ), vt, ldvt,
3327 $ work( iwork ), lwork-iwork+1, ierr )
3332 CALL sorgbr(
'Q', m, m, m, u, ldu, work( itauq ),
3333 $ work( iwork ), lwork-iwork+1, ierr )
3341 CALL sbdsqr(
'U', m, n, m, 0, s, work( ie ), vt,
3342 $ ldvt, u, ldu, dum, 1, work( iwork ),
3366 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
3367 $ work( itaup ), work( iwork ), lwork-iwork+1,
3375 CALL slacpy(
'L', m, m, a, lda, u, ldu )
3376 CALL sorgbr(
'Q', m, m, n, u, ldu, work( itauq ),
3377 $ work( iwork ), lwork-iwork+1, ierr )
3385 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
3390 CALL sorgbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3391 $ work( iwork ), lwork-iwork+1, ierr )
3399 CALL sorgbr(
'Q', m, m, n, a, lda, work( itauq ),
3400 $ work( iwork ), lwork-iwork+1, ierr )
3408 CALL sorgbr(
'P', m, n, m, a, lda, work( itaup ),
3409 $ work( iwork ), lwork-iwork+1, ierr )
3412 IF( wntuas .OR. wntuo )
3416 IF( wntvas .OR. wntvo )
3420 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) )
THEN
3427 CALL sbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3428 $ ldvt, u, ldu, dum, 1, work( iwork ), info )
3429 ELSE IF( ( .NOT.wntuo ) .AND. wntvo )
THEN
3436 CALL sbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), a, lda,
3437 $ u, ldu, dum, 1, work( iwork ), info )
3445 CALL sbdsqr(
'L', m, ncvt, nru, 0, s, work( ie ), vt,
3446 $ ldvt, a, lda, dum, 1, work( iwork ), info )
3456 IF( info.NE.0 )
THEN
3458 DO 50 i = 1, minmn - 1
3459 work( i+1 ) = work( i+ie-1 )
3463 DO 60 i = minmn - 1, 1, -1
3464 work( i+1 ) = work( i+ie-1 )
3471 IF( iscl.EQ.1 )
THEN
3472 IF( anrm.GT.bignum )
3473 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3475 IF( info.NE.0 .AND. anrm.GT.bignum )
3476 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1, work( 2 ),
3478 IF( anrm.LT.smlnum )
3479 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3481 IF( info.NE.0 .AND. anrm.LT.smlnum )
3482 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1, work( 2 ),
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
SGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ