176 SUBROUTINE clalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
177 $ RANK, WORK, RWORK, IWORK, INFO )
185 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
190 REAL D( * ), E( * ), RWORK( * )
191 COMPLEX B( LDB, * ), WORK( * )
198 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
200 parameter( czero = ( 0.0e0, 0.0e0 ) )
203 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
204 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
205 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
206 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
207 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
209 REAL CS, EPS, ORGNRM, R, RCND, SN, TOL
214 EXTERNAL isamax, slamch, slanst
223 INTRINSIC abs, aimag, cmplx, int, log, real, sign
233 ELSE IF( nrhs.LT.1 )
THEN
235 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
239 CALL xerbla(
'CLALSD', -info )
243 eps = slamch(
'Epsilon' )
247 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
259 ELSE IF( n.EQ.1 )
THEN
260 IF( d( 1 ).EQ.zero )
THEN
261 CALL claset(
'A', 1, nrhs, czero, czero, b, ldb )
264 CALL clascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb,
266 d( 1 ) = abs( d( 1 ) )
273 IF( uplo.EQ.
'L' )
THEN
275 CALL slartg( d( i ), e( i ), cs, sn, r )
278 d( i+1 ) = cs*d( i+1 )
280 CALL csrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
291 CALL csrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs,
301 orgnrm = slanst(
'M', n, d, e )
302 IF( orgnrm.EQ.zero )
THEN
303 CALL claset(
'A', n, nrhs, czero, czero, b, ldb )
307 CALL slascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
308 CALL slascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
313 IF( n.LE.smlsiz )
THEN
318 irwib = irwrb + n*nrhs
319 irwb = irwib + n*nrhs
320 CALL slaset(
'A', n, n, zero, one, rwork( irwu ), n )
321 CALL slaset(
'A', n, n, zero, one, rwork( irwvt ), n )
322 CALL slasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
323 $ rwork( irwu ), n, rwork( irwwrk ), 1,
324 $ rwork( irwwrk ), info )
337 rwork( j ) = real( b( jrow, jcol ) )
340 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
341 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
346 rwork( j ) = aimag( b( jrow, jcol ) )
349 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
350 $ rwork( irwb ), n, zero, rwork( irwib ), n )
357 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
361 tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
363 IF( d( i ).LE.tol )
THEN
364 CALL claset(
'A', 1, nrhs, czero, czero, b( i, 1 ),
367 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i,
382 DO 120 jcol = 1, nrhs
385 rwork( j ) = real( b( jrow, jcol ) )
388 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
389 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
391 DO 140 jcol = 1, nrhs
394 rwork( j ) = aimag( b( jrow, jcol ) )
397 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
398 $ rwork( irwb ), n, zero, rwork( irwib ), n )
401 DO 160 jcol = 1, nrhs
405 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
411 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
412 CALL slasrt(
'D', n, d, info )
413 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
420 nlvl = int( log( real( n ) / real( smlsiz+1 ) ) / log( two ) ) + 1
432 givnum = poles + 2*nlvl*n
433 nrwork = givnum + 2*nlvl*n
437 irwib = irwrb + smlsiz*nrhs
438 irwb = irwib + smlsiz*nrhs
444 givcol = perm + nlvl*n
445 iwk = givcol + nlvl*n*2
454 IF( abs( d( i ) ).LT.eps )
THEN
455 d( i ) = sign( eps, d( i ) )
460 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
472 iwork( sizei+nsub-1 ) = nsize
473 ELSE IF( abs( e( i ) ).GE.eps )
THEN
478 iwork( sizei+nsub-1 ) = nsize
486 iwork( sizei+nsub-1 ) = nsize
489 iwork( sizei+nsub-1 ) = 1
490 CALL ccopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
493 IF( nsize.EQ.1 )
THEN
498 CALL ccopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
499 ELSE IF( nsize.LE.smlsiz )
THEN
503 CALL slaset(
'A', nsize, nsize, zero, one,
504 $ rwork( vt+st1 ), n )
505 CALL slaset(
'A', nsize, nsize, zero, one,
506 $ rwork( u+st1 ), n )
507 CALL slasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
508 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
509 $ n, rwork( nrwork ), 1, rwork( nrwork ),
520 DO 190 jcol = 1, nrhs
521 DO 180 jrow = st, st + nsize - 1
523 rwork( j ) = real( b( jrow, jcol ) )
526 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
527 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
528 $ zero, rwork( irwrb ), nsize )
530 DO 210 jcol = 1, nrhs
531 DO 200 jrow = st, st + nsize - 1
533 rwork( j ) = aimag( b( jrow, jcol ) )
536 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
537 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
538 $ zero, rwork( irwib ), nsize )
541 DO 230 jcol = 1, nrhs
542 DO 220 jrow = st, st + nsize - 1
545 b( jrow, jcol ) = cmplx( rwork( jreal ),
550 CALL clacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
551 $ work( bx+st1 ), n )
556 CALL slasda( icmpq1, smlsiz, nsize, sqre, d( st ),
557 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
558 $ iwork( k+st1 ), rwork( difl+st1 ),
559 $ rwork( difr+st1 ), rwork( z+st1 ),
560 $ rwork( poles+st1 ), iwork( givptr+st1 ),
561 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
562 $ rwork( givnum+st1 ), rwork( c+st1 ),
563 $ rwork( s+st1 ), rwork( nrwork ),
564 $ iwork( iwk ), info )
569 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
570 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
571 $ rwork( vt+st1 ), iwork( k+st1 ),
572 $ rwork( difl+st1 ), rwork( difr+st1 ),
573 $ rwork( z+st1 ), rwork( poles+st1 ),
574 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
575 $ iwork( perm+st1 ), rwork( givnum+st1 ),
576 $ rwork( c+st1 ), rwork( s+st1 ),
577 $ rwork( nrwork ), iwork( iwk ), info )
588 tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
595 IF( abs( d( i ) ).LE.tol )
THEN
596 CALL claset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ),
600 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs,
601 $ work( bx+i-1 ), n, info )
603 d( i ) = abs( d( i ) )
612 nsize = iwork( sizei+i-1 )
614 IF( nsize.EQ.1 )
THEN
615 CALL ccopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
616 ELSE IF( nsize.LE.smlsiz )
THEN
627 DO 270 jcol = 1, nrhs
629 DO 260 jrow = 1, nsize
631 rwork( jreal ) = real( work( j+jrow ) )
634 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
635 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
636 $ rwork( irwrb ), nsize )
639 DO 290 jcol = 1, nrhs
641 DO 280 jrow = 1, nsize
643 rwork( jimag ) = aimag( work( j+jrow ) )
646 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
647 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
648 $ rwork( irwib ), nsize )
651 DO 310 jcol = 1, nrhs
652 DO 300 jrow = st, st + nsize - 1
655 b( jrow, jcol ) = cmplx( rwork( jreal ),
660 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ),
662 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
663 $ rwork( vt+st1 ), iwork( k+st1 ),
664 $ rwork( difl+st1 ), rwork( difr+st1 ),
665 $ rwork( z+st1 ), rwork( poles+st1 ),
666 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
667 $ iwork( perm+st1 ), rwork( givnum+st1 ),
668 $ rwork( c+st1 ), rwork( s+st1 ),
669 $ rwork( nrwork ), iwork( iwk ), info )
678 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
679 CALL slasrt(
'D', n, d, info )
680 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )