263 SUBROUTINE clalsa( 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 REAL C( * ), DIFL( ldu, * ), DIFR( ldu, * ),
280 $ givnum( ldu, * ), poles( ldu, * ), rwork( * ),
281 $ s( * ), u( ldu, * ), vt( ldu, * ), z( ldu, * )
282 COMPLEX B( ldb, * ), BX( ldbx, * )
289 parameter( zero = 0.0e0, one = 1.0e0 )
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 aimag, cmplx, real
308 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) )
THEN 310 ELSE IF( smlsiz.LT.3 )
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(
'CLALSA', -info )
336 CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
337 $ iwork( ndimr ), smlsiz )
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 ) =
REAL( B( JROW, JCOL ) )
380 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
389 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
411 DO 70 jrow = nrf, nrf + nr - 1
413 rwork( j ) =
REAL( B( JROW, JCOL ) )
416 CALL sgemm(
'T',
'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
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 ) = aimag( b( jrow, jcol ) )
425 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
445 ic = iwork( inode+i-1 )
446 CALL ccopy( 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 clals0( 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 clals0( 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 ) =
REAL( B( JROW, JCOL ) )
563 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
573 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
594 DO 270 jcol = 1, nrhs
595 DO 260 jrow = nrf, nrf + nrp1 - 1
597 rwork( j ) =
REAL( B( JROW, JCOL ) )
600 CALL sgemm(
'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 ) = aimag( b( jrow, jcol ) )
610 CALL sgemm(
'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 ) = cmplx( rwork( jreal ),
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 xerbla(SRNAME, INFO)
XERBLA
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
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 slasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.