slalsa()

subroutine slalsa	(	integer	icompq,
		integer	smlsiz,
		integer	n,
		integer	nrhs,
		real, dimension( ldb, * )	b,
		integer	ldb,
		real, dimension( ldbx, * )	bx,
		integer	ldbx,
		real, dimension( ldu, * )	u,
		integer	ldu,
		real, dimension( ldu, * )	vt,
		integer, dimension( * )	k,
		real, dimension( ldu, * )	difl,
		real, dimension( ldu, * )	difr,
		real, dimension( ldu, * )	z,
		real, dimension( ldu, * )	poles,
		integer, dimension( * )	givptr,
		integer, dimension( ldgcol, * )	givcol,
		integer	ldgcol,
		integer, dimension( ldgcol, * )	perm,
		real, dimension( ldu, * )	givnum,
		real, dimension( * )	c,
		real, dimension( * )	s,
		real, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

SLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.

Download SLALSA + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> SLALSA is an intermediate step in solving the least squares problem
!> by computing the SVD of the coefficient matrix in compact form (The
!> singular vectors are computed as products of simple orthogonal
!> matrices.).
!>
!> If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector
!> matrix of an upper bidiagonal matrix to the right hand side; and if
!> ICOMPQ = 1, SLALSA applies the right singular vector matrix to the
!> right hand side. The singular vector matrices were generated in
!> compact form by SLALSA.
!> 

Parameters

[in]	ICOMPQ	!> ICOMPQ is INTEGER !> Specifies whether the left or the right singular vector !> matrix is involved. !> = 0: Left singular vector matrix !> = 1: Right singular vector matrix !>
[in]	SMLSIZ	!> SMLSIZ is INTEGER !> The maximum size of the subproblems at the bottom of the !> computation tree. !>
[in]	N	!> N is INTEGER !> The row and column dimensions of the upper bidiagonal matrix. !>
[in]	NRHS	!> NRHS is INTEGER !> The number of columns of B and BX. NRHS must be at least 1. !>
[in,out]	B	!> B is REAL array, dimension ( LDB, NRHS ) !> On input, B contains the right hand sides of the least !> squares problem in rows 1 through M. !> On output, B contains the solution X in rows 1 through N. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of B in the calling subprogram. !> LDB must be at least max(1,MAX( M, N ) ). !>
[out]	BX	!> BX is REAL array, dimension ( LDBX, NRHS ) !> On exit, the result of applying the left or right singular !> vector matrix to B. !>
[in]	LDBX	!> LDBX is INTEGER !> The leading dimension of BX. !>
[in]	U	!> U is REAL array, dimension ( LDU, SMLSIZ ). !> On entry, U contains the left singular vector matrices of all !> subproblems at the bottom level. !>
[in]	LDU	!> LDU is INTEGER, LDU = > N. !> The leading dimension of arrays U, VT, DIFL, DIFR, !> POLES, GIVNUM, and Z. !>
[in]	VT	!> VT is REAL array, dimension ( LDU, SMLSIZ+1 ). !> On entry, VT**T contains the right singular vector matrices of !> all subproblems at the bottom level. !>
[in]	K	!> K is INTEGER array, dimension ( N ). !>
[in]	DIFL	!> DIFL is REAL array, dimension ( LDU, NLVL ). !> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. !>
[in]	DIFR	!> DIFR is REAL array, dimension ( LDU, 2 * NLVL ). !> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record !> distances between singular values on the I-th level and !> singular values on the (I -1)-th level, and DIFR(*, 2 * I) !> record the normalizing factors of the right singular vectors !> matrices of subproblems on I-th level. !>
[in]	Z	!> Z is REAL array, dimension ( LDU, NLVL ). !> On entry, Z(1, I) contains the components of the deflation- !> adjusted updating row vector for subproblems on the I-th !> level. !>
[in]	POLES	!> POLES is REAL array, dimension ( LDU, 2 * NLVL ). !> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old !> singular values involved in the secular equations on the I-th !> level. !>
[in]	GIVPTR	!> GIVPTR is INTEGER array, dimension ( N ). !> On entry, GIVPTR( I ) records the number of Givens !> rotations performed on the I-th problem on the computation !> tree. !>
[in]	GIVCOL	!> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ). !> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the !> locations of Givens rotations performed on the I-th level on !> the computation tree. !>
[in]	LDGCOL	!> LDGCOL is INTEGER, LDGCOL = > N. !> The leading dimension of arrays GIVCOL and PERM. !>
[in]	PERM	!> PERM is INTEGER array, dimension ( LDGCOL, NLVL ). !> On entry, PERM(*, I) records permutations done on the I-th !> level of the computation tree. !>
[in]	GIVNUM	!> GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ). !> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- !> values of Givens rotations performed on the I-th level on the !> computation tree. !>
[in]	C	!> C is REAL array, dimension ( N ). !> On entry, if the I-th subproblem is not square, !> C( I ) contains the C-value of a Givens rotation related to !> the right null space of the I-th subproblem. !>
[in]	S	!> S is REAL array, dimension ( N ). !> On entry, if the I-th subproblem is not square, !> S( I ) contains the S-value of a Givens rotation related to !> the right null space of the I-th subproblem. !>
[out]	WORK	!> WORK is REAL array, dimension (N) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (3*N) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 261 of file slalsa.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
     $                   SMLSIZ
*     ..
*     .. Array Arguments ..
      INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
     $                   K( * ), PERM( LDGCOL, * )
      REAL               B( LDB, * ), BX( LDBX, * ), C( * ),
     $                   DIFL( LDU, * ), DIFR( LDU, * ),
     $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
     $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
     $                   Z( LDU, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
     $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
     $                   NR, NRF, NRP1, SQRE
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgemm, slals0, slasdt,
     $                   xerbla
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
         info = -1
      ELSE IF( smlsiz.LT.3 ) THEN
         info = -2
      ELSE IF( n.LT.smlsiz ) THEN
         info = -3
      ELSE IF( nrhs.LT.1 ) THEN
         info = -4
      ELSE IF( ldb.LT.n ) THEN
         info = -6
      ELSE IF( ldbx.LT.n ) THEN
         info = -8
      ELSE IF( ldu.LT.n ) THEN
         info = -10
      ELSE IF( ldgcol.LT.n ) THEN
         info = -19
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SLALSA', -info )
         RETURN
      END IF
*
*     Book-keeping and  setting up the computation tree.
*
      inode = 1
      ndiml = inode + n
      ndimr = ndiml + n
*
      CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
     $             iwork( ndimr ), smlsiz )
*
*     The following code applies back the left singular vector factors.
*     For applying back the right singular vector factors, go to 50.
*
      IF( icompq.EQ.1 ) THEN
         GO TO 50
      END IF
*
*     The nodes on the bottom level of the tree were solved
*     by SLASDQ. The corresponding left and right singular vector
*     matrices are in explicit form. First apply back the left
*     singular vector matrices.
*
      ndb1 = ( nd+1 ) / 2
      DO 10 i = ndb1, nd
*
*        IC : center row of each node
*        NL : number of rows of left  subproblem
*        NR : number of rows of right subproblem
*        NLF: starting row of the left   subproblem
*        NRF: starting row of the right  subproblem
*
         i1 = i - 1
         ic = iwork( inode+i1 )
         nl = iwork( ndiml+i1 )
         nr = iwork( ndimr+i1 )
         nlf = ic - nl
         nrf = ic + 1
         CALL sgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
     $               b( nlf, 1 ), ldb, zero, bx( nlf, 1 ), ldbx )
         CALL sgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
     $               b( nrf, 1 ), ldb, zero, bx( nrf, 1 ), ldbx )
   10 CONTINUE
*
*     Next copy the rows of B that correspond to unchanged rows
*     in the bidiagonal matrix to BX.
*
      DO 20 i = 1, nd
         ic = iwork( inode+i-1 )
         CALL scopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
   20 CONTINUE
*
*     Finally go through the left singular vector matrices of all
*     the other subproblems bottom-up on the tree.
*
      j = 2**nlvl
      sqre = 0
*
      DO 40 lvl = nlvl, 1, -1
         lvl2 = 2*lvl - 1
*
*        find the first node LF and last node LL on
*        the current level LVL
*
         IF( lvl.EQ.1 ) THEN
            lf = 1
            ll = 1
         ELSE
            lf = 2**( lvl-1 )
            ll = 2*lf - 1
         END IF
         DO 30 i = lf, ll
            im1 = i - 1
            ic = iwork( inode+im1 )
            nl = iwork( ndiml+im1 )
            nr = iwork( ndimr+im1 )
            nlf = ic - nl
            nrf = ic + 1
            j = j - 1
            CALL slals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ),
     $                   ldbx,
     $                   b( nlf, 1 ), ldb, perm( nlf, lvl ),
     $                   givptr( j ), givcol( nlf, lvl2 ), ldgcol,
     $                   givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
     $                   difl( nlf, lvl ), difr( nlf, lvl2 ),
     $                   z( nlf, lvl ), k( j ), c( j ), s( j ), work,
     $                   info )
   30    CONTINUE
   40 CONTINUE
      GO TO 90
*
*     ICOMPQ = 1: applying back the right singular vector factors.
*
   50 CONTINUE
*
*     First now go through the right singular vector matrices of all
*     the tree nodes top-down.
*
      j = 0
      DO 70 lvl = 1, nlvl
         lvl2 = 2*lvl - 1
*
*        Find the first node LF and last node LL on
*        the current level LVL.
*
         IF( lvl.EQ.1 ) THEN
            lf = 1
            ll = 1
         ELSE
            lf = 2**( lvl-1 )
            ll = 2*lf - 1
         END IF
         DO 60 i = ll, lf, -1
            im1 = i - 1
            ic = iwork( inode+im1 )
            nl = iwork( ndiml+im1 )
            nr = iwork( ndimr+im1 )
            nlf = ic - nl
            nrf = ic + 1
            IF( i.EQ.ll ) THEN
               sqre = 0
            ELSE
               sqre = 1
            END IF
            j = j + 1
            CALL slals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ),
     $                   ldb,
     $                   bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
     $                   givptr( j ), givcol( nlf, lvl2 ), ldgcol,
     $                   givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
     $                   difl( nlf, lvl ), difr( nlf, lvl2 ),
     $                   z( nlf, lvl ), k( j ), c( j ), s( j ), work,
     $                   info )
   60    CONTINUE
   70 CONTINUE
*
*     The nodes on the bottom level of the tree were solved
*     by SLASDQ. The corresponding right singular vector
*     matrices are in explicit form. Apply them back.
*
      ndb1 = ( nd+1 ) / 2
      DO 80 i = ndb1, nd
         i1 = i - 1
         ic = iwork( inode+i1 )
         nl = iwork( ndiml+i1 )
         nr = iwork( ndimr+i1 )
         nlp1 = nl + 1
         IF( i.EQ.nd ) THEN
            nrp1 = nr
         ELSE
            nrp1 = nr + 1
         END IF
         nlf = ic - nl
         nrf = ic + 1
         CALL sgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ),
     $               ldu,
     $               b( nlf, 1 ), ldb, zero, bx( nlf, 1 ), ldbx )
         CALL sgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ),
     $               ldu,
     $               b( nrf, 1 ), ldb, zero, bx( nrf, 1 ), ldbx )
   80 CONTINUE
*
   90 CONTINUE
*
      RETURN
*
*     End of SLALSA
*

Here is the call graph for this function:

Here is the caller graph for this function: