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
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...
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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
logical function lsame(CA, CB)
LSAME
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
real function slamch(CMACH)
SLAMCH
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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