263 SUBROUTINE zlalsa( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
264 $ LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
265 $ GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK,
273 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
277 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
278 $ K( * ), PERM( LDGCOL, * )
279 DOUBLE PRECISION ( * ), DIFL( LDU, * ), DIFR( LDU, * ),
280 $ givnum( ldu, * ), poles( ldu, * ), rwork( * ),
281 $ s( * ), u( ldu, * ), vt( ldu, * ), z( ldu, * )
282 COMPLEX*16 B( LDB, * ), BX( LDBX, * )
288 DOUBLE PRECISION ZERO, ONE
289 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
292 INTEGER I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
293 $ JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
294 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
300 INTRINSIC dble, dcmplx, dimag
308 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) )
THEN
312 ELSE IF( n.LT.smlsiz )
THEN
314 ELSE IF( nrhs.LT.1 )
THEN
316 ELSE IF( ldb.LT.n )
THEN
318 ELSE IF( ldbx.LT.n )
THEN
320 ELSE IF( ldu.LT.n )
THEN
322 ELSE IF( ldgcol.LT.n )
THEN
326 CALL xerbla(
'ZLALSA', -info )
336 CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
342 IF( icompq.EQ.1 )
THEN
361 ic = iwork( inode+i1 )
362 nl = iwork( ndiml+i1 )
363 nr = iwork( ndimr+i1 )
375 DO 10 jrow = nlf, nlf + nl - 1
377 rwork( j ) = dble( b( jrow, jcol ) )
380 CALL dgemm(
'T',
'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
381 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1 ), nl )
384 DO 30 jrow = nlf, nlf + nl - 1
386 rwork( j ) = dimag( b( jrow, jcol ) )
389 CALL dgemm(
'T',
'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
390 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1+nl*nrhs ),
395 DO 50 jrow = nlf, nlf + nl - 1
398 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
411 DO 70 jrow = nrf, nrf + nr - 1
413 rwork( j ) = dble( b( jrow, jcol ) )
417 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1 ), nr )
419 DO 100 jcol = 1, nrhs
420 DO 90 jrow = nrf, nrf + nr - 1
422 rwork( j ) = dimag( b( jrow, jcol ) )
425 CALL dgemm(
'T',
'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
426 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1+nr*nrhs ),
430 DO 120 jcol = 1, nrhs
431 DO 110 jrow = nrf, nrf + nr - 1
434 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
445 ic = iwork( inode+i-1 )
446 CALL zcopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
455 DO 160 lvl = nlvl, 1, -1
470 ic = iwork( inode+im1 )
471 nl = iwork( ndiml+im1 )
472 nr = iwork( ndimr+im1 )
476 CALL zlals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ), ldbx,
477 $ b( nlf, 1 ), ldb, perm( nlf, lvl ),
478 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
479 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
480 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
481 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
508 DO 180 i = ll, lf, -1
510 ic = iwork( inode+im1 )
511 nl = iwork( ndiml+im1 )
512 nr = iwork( ndimr+im1 )
521 CALL zlals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ), ldb,
522 $ bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
523 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
524 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
525 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
526 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
538 ic = iwork( inode+i1 )
539 nl = iwork( ndiml+i1 )
540 nr = iwork( ndimr+i1 )
557 DO 210 jcol = 1, nrhs
558 DO 200 jrow = nlf, nlf + nlp1 - 1
560 rwork( j ) = dble( b( jrow, jcol ) )
563 CALL dgemm(
'T',
'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
564 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero, rwork( 1 ),
567 DO 230 jcol = 1, nrhs
568 DO 220 jrow = nlf, nlf + nlp1 - 1
570 rwork( j ) = dimag( b( jrow, jcol ) )
573 CALL dgemm(
'T',
'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
574 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero,
575 $ rwork( 1+nlp1*nrhs ), nlp1 )
578 DO 250 jcol = 1, nrhs
579 DO 240 jrow = nlf, nlf + nlp1 - 1
582 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
594 DO 270 jcol = 1, nrhs
595 DO 260 jrow = nrf, nrf + nrp1 - 1
597 rwork( j ) = dble( b( jrow, jcol ) )
600 CALL dgemm(
'T',
'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
601 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero, rwork( 1 ),
604 DO 290 jcol = 1, nrhs
605 DO 280 jrow = nrf, nrf + nrp1 - 1
607 rwork( j ) = dimag( b( jrow, jcol ) )
610 CALL dgemm(
'T',
'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
611 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero,
612 $ rwork( 1+nrp1*nrhs ), nrp1 )
615 DO 310 jcol = 1, nrhs
616 DO 300 jrow = nrf, nrf + nrp1 - 1
619 bx( jrow, jcol ) = dcmplx( rwork( jreal ),