186 SUBROUTINE clalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
187 $ rank, work, rwork, iwork, info )
196 INTEGER info, ldb, n, nrhs, rank, smlsiz
201 REAL d( * ), e( * ), rwork( * )
202 COMPLEX b( ldb, * ), work( * )
209 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
211 parameter( czero = ( 0.0e0, 0.0e0 ) )
214 INTEGER bx, bxst, c, difl, difr, givcol, givnum,
215 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
216 $ irwu, irwvt, irwwrk, iwk,
j, jcol, jimag,
217 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
218 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
220 REAL cs, eps, orgnrm, r, rcnd, sn, tol
233 INTRINSIC abs, aimag, cmplx, int, log,
REAL, sign
243 ELSE IF( nrhs.LT.1 )
THEN
245 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
249 CALL
xerbla(
'CLALSD', -info )
257 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
269 ELSE IF( n.EQ.1 )
THEN
270 IF( d( 1 ).EQ.zero )
THEN
271 CALL
claset(
'A', 1, nrhs, czero, czero,
b, ldb )
274 CALL
clascl(
'G', 0, 0, d( 1 ), one, 1, nrhs,
b, ldb, info )
275 d( 1 ) = abs( d( 1 ) )
282 IF( uplo.EQ.
'L' )
THEN
284 CALL
slartg( d( i ), e( i ), cs, sn, r )
287 d( i+1 ) = cs*d( i+1 )
289 CALL
csrot( 1,
b( i, 1 ), 1,
b( i+1, 1 ), 1, cs, sn )
300 CALL
csrot( 1,
b(
j, i ), 1,
b(
j+1, i ), 1, cs, sn )
309 orgnrm =
slanst(
'M', n, d, e )
310 IF( orgnrm.EQ.zero )
THEN
311 CALL
claset(
'A', n, nrhs, czero, czero,
b, ldb )
315 CALL
slascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
316 CALL
slascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
321 IF( n.LE.smlsiz )
THEN
326 irwib = irwrb + n*nrhs
327 irwb = irwib + n*nrhs
328 CALL
slaset(
'A', n, n, zero, one, rwork( irwu ), n )
329 CALL
slaset(
'A', n, n, zero, one, rwork( irwvt ), n )
330 CALL
slasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
331 $ rwork( irwu ), n, rwork( irwwrk ), 1,
332 $ rwork( irwwrk ), info )
345 rwork(
j ) =
REAL( B( JROW, JCOL ) )
348 CALL
sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
349 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
354 rwork(
j ) = aimag(
b( jrow, jcol ) )
357 CALL
sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
358 $ rwork( irwb ), n, zero, rwork( irwib ), n )
365 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
369 tol = rcnd*abs( d(
isamax( n, d, 1 ) ) )
371 IF( d( i ).LE.tol )
THEN
372 CALL
claset(
'A', 1, nrhs, czero, czero,
b( i, 1 ), ldb )
374 CALL
clascl(
'G', 0, 0, d( i ), one, 1, nrhs,
b( i, 1 ),
388 DO 120 jcol = 1, nrhs
391 rwork(
j ) =
REAL( B( JROW, JCOL ) )
394 CALL
sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
395 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
397 DO 140 jcol = 1, nrhs
400 rwork(
j ) = aimag(
b( jrow, jcol ) )
403 CALL
sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
404 $ rwork( irwb ), n, zero, rwork( irwib ), n )
407 DO 160 jcol = 1, nrhs
411 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
417 CALL
slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
418 CALL
slasrt(
'D', n, d, info )
419 CALL
clascl(
'G', 0, 0, orgnrm, one, n, nrhs,
b, ldb, info )
426 nlvl = int( log(
REAL( N ) /
REAL( SMLSIZ+1 ) ) / log( two ) ) + 1
438 givnum = poles + 2*nlvl*n
439 nrwork = givnum + 2*nlvl*n
443 irwib = irwrb + smlsiz*nrhs
444 irwb = irwib + smlsiz*nrhs
450 givcol = perm + nlvl*n
451 iwk = givcol + nlvl*n*2
460 IF( abs( d( i ) ).LT.eps )
THEN
461 d( i ) = sign( eps, d( i ) )
466 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
478 iwork( sizei+nsub-1 ) = nsize
479 ELSE IF( abs( e( i ) ).GE.eps )
THEN
484 iwork( sizei+nsub-1 ) = nsize
492 iwork( sizei+nsub-1 ) = nsize
495 iwork( sizei+nsub-1 ) = 1
496 CALL
ccopy( nrhs,
b( n, 1 ), ldb, work( bx+nm1 ), n )
499 IF( nsize.EQ.1 )
THEN
504 CALL
ccopy( nrhs,
b( st, 1 ), ldb, work( bx+st1 ), n )
505 ELSE IF( nsize.LE.smlsiz )
THEN
509 CALL
slaset(
'A', nsize, nsize, zero, one,
510 $ rwork( vt+st1 ), n )
511 CALL
slaset(
'A', nsize, nsize, zero, one,
512 $ rwork( u+st1 ), n )
513 CALL
slasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
514 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
515 $ n, rwork( nrwork ), 1, rwork( nrwork ),
526 DO 190 jcol = 1, nrhs
527 DO 180 jrow = st, st + nsize - 1
529 rwork(
j ) =
REAL( B( JROW, JCOL ) )
532 CALL
sgemm(
'T',
'N', nsize, nrhs, nsize, one,
533 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
534 $ zero, rwork( irwrb ), nsize )
536 DO 210 jcol = 1, nrhs
537 DO 200 jrow = st, st + nsize - 1
539 rwork(
j ) = aimag(
b( jrow, jcol ) )
542 CALL
sgemm(
'T',
'N', nsize, nrhs, nsize, one,
543 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
544 $ zero, rwork( irwib ), nsize )
547 DO 230 jcol = 1, nrhs
548 DO 220 jrow = st, st + nsize - 1
551 b( jrow, jcol ) = cmplx( rwork( jreal ),
556 CALL
clacpy(
'A', nsize, nrhs,
b( st, 1 ), ldb,
557 $ work( bx+st1 ), n )
562 CALL
slasda( icmpq1, smlsiz, nsize, sqre, d( st ),
563 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
564 $ iwork( k+st1 ), rwork( difl+st1 ),
565 $ rwork( difr+st1 ), rwork( z+st1 ),
566 $ rwork( poles+st1 ), iwork( givptr+st1 ),
567 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
568 $ rwork( givnum+st1 ), rwork( c+st1 ),
569 $ rwork( s+st1 ), rwork( nrwork ),
570 $ iwork( iwk ), info )
575 CALL
clalsa( icmpq2, smlsiz, nsize, nrhs,
b( st, 1 ),
576 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
577 $ rwork( vt+st1 ), iwork( k+st1 ),
578 $ rwork( difl+st1 ), rwork( difr+st1 ),
579 $ rwork( z+st1 ), rwork( poles+st1 ),
580 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
581 $ iwork( perm+st1 ), rwork( givnum+st1 ),
582 $ rwork( c+st1 ), rwork( s+st1 ),
583 $ rwork( nrwork ), iwork( iwk ), info )
594 tol = rcnd*abs( d(
isamax( n, d, 1 ) ) )
601 IF( abs( d( i ) ).LE.tol )
THEN
602 CALL
claset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
605 CALL
clascl(
'G', 0, 0, d( i ), one, 1, nrhs,
606 $ work( bx+i-1 ), n, info )
608 d( i ) = abs( d( i ) )
617 nsize = iwork( sizei+i-1 )
619 IF( nsize.EQ.1 )
THEN
620 CALL
ccopy( nrhs, work( bxst ), n,
b( st, 1 ), ldb )
621 ELSE IF( nsize.LE.smlsiz )
THEN
632 DO 270 jcol = 1, nrhs
634 DO 260 jrow = 1, nsize
636 rwork( jreal ) =
REAL( WORK( J+JROW ) )
639 CALL
sgemm(
'T',
'N', nsize, nrhs, nsize, one,
640 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
641 $ rwork( irwrb ), nsize )
644 DO 290 jcol = 1, nrhs
646 DO 280 jrow = 1, nsize
648 rwork( jimag ) = aimag( work(
j+jrow ) )
651 CALL
sgemm(
'T',
'N', nsize, nrhs, nsize, one,
652 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
653 $ rwork( irwib ), nsize )
656 DO 310 jcol = 1, nrhs
657 DO 300 jrow = st, st + nsize - 1
660 b( jrow, jcol ) = cmplx( rwork( jreal ),
665 CALL
clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
666 $
b( st, 1 ), ldb, rwork( u+st1 ), n,
667 $ rwork( vt+st1 ), iwork( k+st1 ),
668 $ rwork( difl+st1 ), rwork( difr+st1 ),
669 $ rwork( z+st1 ), rwork( poles+st1 ),
670 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
671 $ iwork( perm+st1 ), rwork( givnum+st1 ),
672 $ rwork( c+st1 ), rwork( s+st1 ),
673 $ rwork( nrwork ), iwork( iwk ), info )
682 CALL
slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
683 CALL
slasrt(
'D', n, d, info )
684 CALL
clascl(
'G', 0, 0, orgnrm, one, n, nrhs,
b, ldb, info )
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 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 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.
INTEGER function isamax(N, SX, INCX)
ISAMAX
REAL function slamch(CMACH)
SLAMCH
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
REAL function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
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...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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 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 clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY 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 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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.