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 xerbla(srname, info)
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
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 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 slasdt(n, lvl, nd, inode, ndiml, ndimr, msub)
SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.