184 SUBROUTINE clalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
185 $ RANK, WORK, RWORK, IWORK, INFO )
193 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
198 REAL D( * ), E( * ), RWORK( * )
199 COMPLEX B( LDB, * ), WORK( * )
206 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
208 parameter( czero = ( 0.0e0, 0.0e0 ) )
211 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
212 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
213 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
214 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
215 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
217 REAL CS, EPS, ORGNRM, R, RCND, SN, TOL
222 EXTERNAL isamax, slamch, slanst
230 INTRINSIC abs, aimag, cmplx, int, log, real, sign
240 ELSE IF( nrhs.LT.1 )
THEN
242 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
246 CALL xerbla(
'CLALSD', -info )
250 eps = slamch(
'Epsilon' )
254 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
266 ELSE IF( n.EQ.1 )
THEN
267 IF( d( 1 ).EQ.zero )
THEN
268 CALL claset(
'A', 1, nrhs, czero, czero, b, ldb )
271 CALL clascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
272 d( 1 ) = abs( d( 1 ) )
279 IF( uplo.EQ.
'L' )
THEN
281 CALL slartg( d( i ), e( i ), cs, sn, r )
284 d( i+1 ) = cs*d( i+1 )
286 CALL csrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
297 CALL csrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
306 orgnrm = slanst(
'M', n, d, e )
307 IF( orgnrm.EQ.zero )
THEN
308 CALL claset(
'A', n, nrhs, czero, czero, b, ldb )
312 CALL slascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
313 CALL slascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
318 IF( n.LE.smlsiz )
THEN
323 irwib = irwrb + n*nrhs
324 irwb = irwib + n*nrhs
325 CALL slaset(
'A', n, n, zero, one, rwork( irwu ), n )
326 CALL slaset(
'A', n, n, zero, one, rwork( irwvt ), n )
327 CALL slasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
328 $ rwork( irwu ), n, rwork( irwwrk ), 1,
329 $ rwork( irwwrk ), info )
342 rwork( j ) = real( b( jrow, jcol ) )
345 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
346 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
351 rwork( j ) = aimag( b( jrow, jcol ) )
354 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
355 $ rwork( irwb ), n, zero, rwork( irwib ), n )
362 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
366 tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
368 IF( d( i ).LE.tol )
THEN
369 CALL claset(
'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
371 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
385 DO 120 jcol = 1, nrhs
388 rwork( j ) = real( b( jrow, jcol ) )
391 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
392 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
394 DO 140 jcol = 1, nrhs
397 rwork( j ) = aimag( b( jrow, jcol ) )
400 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
401 $ rwork( irwb ), n, zero, rwork( irwib ), n )
404 DO 160 jcol = 1, nrhs
408 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
414 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
415 CALL slasrt(
'D', n, d, info )
416 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
423 nlvl = int( log( real( n ) / real( smlsiz+1 ) ) / log( two ) ) + 1
435 givnum = poles + 2*nlvl*n
436 nrwork = givnum + 2*nlvl*n
440 irwib = irwrb + smlsiz*nrhs
441 irwb = irwib + smlsiz*nrhs
447 givcol = perm + nlvl*n
448 iwk = givcol + nlvl*n*2
457 IF( abs( d( i ) ).LT.eps )
THEN
458 d( i ) = sign( eps, d( i ) )
463 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
475 iwork( sizei+nsub-1 ) = nsize
476 ELSE IF( abs( e( i ) ).GE.eps )
THEN
481 iwork( sizei+nsub-1 ) = nsize
489 iwork( sizei+nsub-1 ) = nsize
492 iwork( sizei+nsub-1 ) = 1
493 CALL ccopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
496 IF( nsize.EQ.1 )
THEN
501 CALL ccopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
502 ELSE IF( nsize.LE.smlsiz )
THEN
506 CALL slaset(
'A', nsize, nsize, zero, one,
507 $ rwork( vt+st1 ), n )
508 CALL slaset(
'A', nsize, nsize, zero, one,
509 $ rwork( u+st1 ), n )
510 CALL slasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
511 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
512 $ n, rwork( nrwork ), 1, rwork( nrwork ),
523 DO 190 jcol = 1, nrhs
524 DO 180 jrow = st, st + nsize - 1
526 rwork( j ) = real( b( jrow, jcol ) )
529 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
530 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
531 $ zero, rwork( irwrb ), nsize )
533 DO 210 jcol = 1, nrhs
534 DO 200 jrow = st, st + nsize - 1
536 rwork( j ) = aimag( b( jrow, jcol ) )
539 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
540 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
541 $ zero, rwork( irwib ), nsize )
544 DO 230 jcol = 1, nrhs
545 DO 220 jrow = st, st + nsize - 1
548 b( jrow, jcol ) = cmplx( rwork( jreal ),
553 CALL clacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
554 $ work( bx+st1 ), n )
559 CALL slasda( icmpq1, smlsiz, nsize, sqre, d( st ),
560 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
561 $ iwork( k+st1 ), rwork( difl+st1 ),
562 $ rwork( difr+st1 ), rwork( z+st1 ),
563 $ rwork( poles+st1 ), iwork( givptr+st1 ),
564 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
565 $ rwork( givnum+st1 ), rwork( c+st1 ),
566 $ rwork( s+st1 ), rwork( nrwork ),
567 $ iwork( iwk ), info )
572 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
573 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
574 $ rwork( vt+st1 ), iwork( k+st1 ),
575 $ rwork( difl+st1 ), rwork( difr+st1 ),
576 $ rwork( z+st1 ), rwork( poles+st1 ),
577 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
578 $ iwork( perm+st1 ), rwork( givnum+st1 ),
579 $ rwork( c+st1 ), rwork( s+st1 ),
580 $ rwork( nrwork ), iwork( iwk ), info )
591 tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
598 IF( abs( d( i ) ).LE.tol )
THEN
599 CALL claset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
602 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs,
603 $ work( bx+i-1 ), n, info )
605 d( i ) = abs( d( i ) )
614 nsize = iwork( sizei+i-1 )
616 IF( nsize.EQ.1 )
THEN
617 CALL ccopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
618 ELSE IF( nsize.LE.smlsiz )
THEN
629 DO 270 jcol = 1, nrhs
631 DO 260 jrow = 1, nsize
633 rwork( jreal ) = real( work( j+jrow ) )
636 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
637 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
638 $ rwork( irwrb ), nsize )
641 DO 290 jcol = 1, nrhs
643 DO 280 jrow = 1, nsize
645 rwork( jimag ) = aimag( work( j+jrow ) )
648 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
649 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
650 $ rwork( irwib ), nsize )
653 DO 310 jcol = 1, nrhs
654 DO 300 jrow = st, st + nsize - 1
657 b( jrow, jcol ) = cmplx( rwork( jreal ),
662 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
663 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
664 $ rwork( vt+st1 ), iwork( k+st1 ),
665 $ rwork( difl+st1 ), rwork( difr+st1 ),
666 $ rwork( z+st1 ), rwork( poles+st1 ),
667 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
668 $ iwork( perm+st1 ), rwork( givnum+st1 ),
669 $ rwork( c+st1 ), rwork( s+st1 ),
670 $ rwork( nrwork ), iwork( iwk ), info )
679 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
680 CALL slasrt(
'D', n, d, info )
681 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
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 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 slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine clalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
CLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine clalsa(ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK, IWORK, INFO)
CLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM