266 SUBROUTINE zlalsa( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
267 $ ldu, vt, k, difl, difr, z, poles, givptr,
268 $ givcol, ldgcol, perm, givnum, c, s, rwork,
277 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
281 INTEGER GIVCOL( ldgcol, * ), GIVPTR( * ), IWORK( * ),
282 $ k( * ), perm( ldgcol, * )
283 DOUBLE PRECISION C( * ), DIFL( ldu, * ), DIFR( ldu, * ),
284 $ givnum( ldu, * ), poles( ldu, * ), rwork( * ),
285 $ s( * ), u( ldu, * ), vt( ldu, * ), z( ldu, * )
286 COMPLEX*16 B( ldb, * ), BX( ldbx, * )
292 DOUBLE PRECISION ZERO, ONE
293 parameter ( zero = 0.0d0, one = 1.0d0 )
296 INTEGER I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
297 $ jrow, lf, ll, lvl, lvl2, nd, ndb1, ndiml,
298 $ ndimr, nl, nlf, nlp1, nlvl, nr, nrf, nrp1, sqre
304 INTRINSIC dble, dcmplx, dimag
312 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) )
THEN
314 ELSE IF( smlsiz.LT.3 )
THEN
316 ELSE IF( n.LT.smlsiz )
THEN
318 ELSE IF( nrhs.LT.1 )
THEN
320 ELSE IF( ldb.LT.n )
THEN
322 ELSE IF( ldbx.LT.n )
THEN
324 ELSE IF( ldu.LT.n )
THEN
326 ELSE IF( ldgcol.LT.n )
THEN
330 CALL xerbla(
'ZLALSA', -info )
340 CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
341 $ iwork( ndimr ), smlsiz )
346 IF( icompq.EQ.1 )
THEN
365 ic = iwork( inode+i1 )
366 nl = iwork( ndiml+i1 )
367 nr = iwork( ndimr+i1 )
379 DO 10 jrow = nlf, nlf + nl - 1
381 rwork( j ) = dble( b( jrow, jcol ) )
384 CALL dgemm(
'T',
'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
385 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1 ), nl )
388 DO 30 jrow = nlf, nlf + nl - 1
390 rwork( j ) = dimag( b( jrow, jcol ) )
393 CALL dgemm(
'T',
'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
394 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1+nl*nrhs ),
399 DO 50 jrow = nlf, nlf + nl - 1
402 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
415 DO 70 jrow = nrf, nrf + nr - 1
417 rwork( j ) = dble( b( jrow, jcol ) )
420 CALL dgemm(
'T',
'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
421 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1 ), nr )
423 DO 100 jcol = 1, nrhs
424 DO 90 jrow = nrf, nrf + nr - 1
426 rwork( j ) = dimag( b( jrow, jcol ) )
429 CALL dgemm(
'T',
'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
430 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1+nr*nrhs ),
434 DO 120 jcol = 1, nrhs
435 DO 110 jrow = nrf, nrf + nr - 1
438 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
449 ic = iwork( inode+i-1 )
450 CALL zcopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
459 DO 160 lvl = nlvl, 1, -1
474 ic = iwork( inode+im1 )
475 nl = iwork( ndiml+im1 )
476 nr = iwork( ndimr+im1 )
480 CALL zlals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ), ldbx,
481 $ b( nlf, 1 ), ldb, perm( nlf, lvl ),
482 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
483 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
484 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
485 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
512 DO 180 i = ll, lf, -1
514 ic = iwork( inode+im1 )
515 nl = iwork( ndiml+im1 )
516 nr = iwork( ndimr+im1 )
525 CALL zlals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ), ldb,
526 $ bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
527 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
528 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
529 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
530 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
542 ic = iwork( inode+i1 )
543 nl = iwork( ndiml+i1 )
544 nr = iwork( ndimr+i1 )
561 DO 210 jcol = 1, nrhs
562 DO 200 jrow = nlf, nlf + nlp1 - 1
564 rwork( j ) = dble( b( jrow, jcol ) )
567 CALL dgemm(
'T',
'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
568 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero, rwork( 1 ),
571 DO 230 jcol = 1, nrhs
572 DO 220 jrow = nlf, nlf + nlp1 - 1
574 rwork( j ) = dimag( b( jrow, jcol ) )
577 CALL dgemm(
'T',
'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
578 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero,
579 $ rwork( 1+nlp1*nrhs ), nlp1 )
582 DO 250 jcol = 1, nrhs
583 DO 240 jrow = nlf, nlf + nlp1 - 1
586 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
598 DO 270 jcol = 1, nrhs
599 DO 260 jrow = nrf, nrf + nrp1 - 1
601 rwork( j ) = dble( b( jrow, jcol ) )
604 CALL dgemm(
'T',
'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
605 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero, rwork( 1 ),
608 DO 290 jcol = 1, nrhs
609 DO 280 jrow = nrf, nrf + nrp1 - 1
611 rwork( j ) = dimag( b( jrow, jcol ) )
614 CALL dgemm(
'T',
'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
615 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero,
616 $ rwork( 1+nrp1*nrhs ), nrp1 )
619 DO 310 jcol = 1, nrhs
620 DO 300 jrow = nrf, nrf + nrp1 - 1
623 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
subroutine zlals0(ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO)
ZLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dlasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlalsa(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)
ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.