266 SUBROUTINE clalsa( 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 REAL C( * ), DIFL( ldu, * ), DIFR( ldu, * ),
284 $ givnum( ldu, * ), poles( ldu, * ), rwork( * ),
285 $ s( * ), u( ldu, * ), vt( ldu, * ), z( ldu, * )
286 COMPLEX B( ldb, * ), BX( ldbx, * )
293 parameter ( zero = 0.0e0, one = 1.0e0 )
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 aimag, cmplx, real
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(
'CLALSA', -info )
340 CALL slasdt( 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 ) =
REAL( B( JROW, JCOL ) )
384 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
393 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
415 DO 70 jrow = nrf, nrf + nr - 1
417 rwork( j ) =
REAL( B( JROW, JCOL ) )
420 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
429 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
449 ic = iwork( inode+i-1 )
450 CALL ccopy( 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 clals0( 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 clals0( 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 ) =
REAL( B( JROW, JCOL ) )
567 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
577 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
598 DO 270 jcol = 1, nrhs
599 DO 260 jrow = nrf, nrf + nrp1 - 1
601 rwork( j ) =
REAL( B( JROW, JCOL ) )
604 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
614 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
subroutine slasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clals0(ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO)
CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY