zlalsa()

subroutine zlalsa	(	integer	icompq,
		integer	smlsiz,
		integer	n,
		integer	nrhs,
		complex*16, dimension( ldb, * )	b,
		integer	ldb,
		complex*16, dimension( ldbx, * )	bx,
		integer	ldbx,
		double precision, dimension( ldu, * )	u,
		integer	ldu,
		double precision, dimension( ldu, * )	vt,
		integer, dimension( * )	k,
		double precision, dimension( ldu, * )	difl,
		double precision, dimension( ldu, * )	difr,
		double precision, dimension( ldu, * )	z,
		double precision, dimension( ldu, * )	poles,
		integer, dimension( * )	givptr,
		integer, dimension( ldgcol, * )	givcol,
		integer	ldgcol,
		integer, dimension( ldgcol, * )	perm,
		double precision, dimension( ldu, * )	givnum,
		double precision, dimension( * )	c,
		double precision, dimension( * )	s,
		double precision, dimension( * )	rwork,
		integer, dimension( * )	iwork,
		integer	info )

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

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

Purpose:

!>
!> ZLALSA 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, ZLALSA applies the inverse of the left singular vector
!> matrix of an upper bidiagonal matrix to the right hand side; and if
!> ICOMPQ = 1, ZLALSA applies the right singular vector matrix to the
!> right hand side. The singular vector matrices were generated in
!> compact form by ZLALSA.
!> 

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 COMPLEX*16 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 COMPLEX*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ). !> On entry, VT**H 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 DOUBLE PRECISION array, dimension ( LDU, NLVL ). !> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. !>
[in]	DIFR	!> DIFR is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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]	RWORK	!> RWORK is DOUBLE PRECISION array, dimension at least !> MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ). !>
[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 zlalsa.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, * )
      DOUBLE PRECISION   C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
     $                   GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
     $                   S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
      COMPLEX*16         B( LDB, * ), BX( LDBX, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
     $                   JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
     $                   NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemm, dlasdt, xerbla, zcopy,
     $                   zlals0
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, dcmplx, dimag
*     ..
*     .. 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( 'ZLALSA', -info )
         RETURN
      END IF
*
*     Book-keeping and  setting up the computation tree.
*
      inode = 1
      ndiml = inode + n
      ndimr = ndiml + n
*
      CALL dlasdt( 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 170.
*
      IF( icompq.EQ.1 ) THEN
         GO TO 170
      END IF
*
*     The nodes on the bottom level of the tree were solved
*     by DLASDQ. 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 130 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
*
*        Since B and BX are complex, the following call to DGEMM
*        is performed in two steps (real and imaginary parts).
*
*        CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
*     $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
*
         j = nl*nrhs*2
         DO 20 jcol = 1, nrhs
            DO 10 jrow = nlf, nlf + nl - 1
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
   10       CONTINUE
   20    CONTINUE
         CALL dgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
     $               rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1 ), nl )
         j = nl*nrhs*2
         DO 40 jcol = 1, nrhs
            DO 30 jrow = nlf, nlf + nl - 1
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
   30       CONTINUE
   40    CONTINUE
         CALL dgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
     $               rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1+nl*nrhs ),
     $               nl )
         jreal = 0
         jimag = nl*nrhs
         DO 60 jcol = 1, nrhs
            DO 50 jrow = nlf, nlf + nl - 1
               jreal = jreal + 1
               jimag = jimag + 1
               bx( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                            rwork( jimag ) )
   50       CONTINUE
   60    CONTINUE
*
*        Since B and BX are complex, the following call to DGEMM
*        is performed in two steps (real and imaginary parts).
*
*        CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
*    $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
*
         j = nr*nrhs*2
         DO 80 jcol = 1, nrhs
            DO 70 jrow = nrf, nrf + nr - 1
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
   70       CONTINUE
   80    CONTINUE
         CALL dgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
     $               rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1 ), nr )
         j = nr*nrhs*2
         DO 100 jcol = 1, nrhs
            DO 90 jrow = nrf, nrf + nr - 1
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
   90       CONTINUE
  100    CONTINUE
         CALL dgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
     $               rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1+nr*nrhs ),
     $               nr )
         jreal = 0
         jimag = nr*nrhs
         DO 120 jcol = 1, nrhs
            DO 110 jrow = nrf, nrf + nr - 1
               jreal = jreal + 1
               jimag = jimag + 1
               bx( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                            rwork( jimag ) )
  110       CONTINUE
  120    CONTINUE
*
  130 CONTINUE
*
*     Next copy the rows of B that correspond to unchanged rows
*     in the bidiagonal matrix to BX.
*
      DO 140 i = 1, nd
         ic = iwork( inode+i-1 )
         CALL zcopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
  140 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 160 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 150 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 zlals0( 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 ), rwork,
     $                   info )
  150    CONTINUE
  160 CONTINUE
      GO TO 330
*
*     ICOMPQ = 1: applying back the right singular vector factors.
*
  170 CONTINUE
*
*     First now go through the right singular vector matrices of all
*     the tree nodes top-down.
*
      j = 0
      DO 190 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 180 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 zlals0( 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 ), rwork,
     $                   info )
  180    CONTINUE
  190 CONTINUE
*
*     The nodes on the bottom level of the tree were solved
*     by DLASDQ. The corresponding right singular vector
*     matrices are in explicit form. Apply them back.
*
      ndb1 = ( nd+1 ) / 2
      DO 320 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
*
*        Since B and BX are complex, the following call to DGEMM is
*        performed in two steps (real and imaginary parts).
*
*        CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
*    $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
*
         j = nlp1*nrhs*2
         DO 210 jcol = 1, nrhs
            DO 200 jrow = nlf, nlf + nlp1 - 1
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
  200       CONTINUE
  210    CONTINUE
         CALL dgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ),
     $               ldu,
     $               rwork( 1+nlp1*nrhs*2 ), nlp1, zero, rwork( 1 ),
     $               nlp1 )
         j = nlp1*nrhs*2
         DO 230 jcol = 1, nrhs
            DO 220 jrow = nlf, nlf + nlp1 - 1
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
  220       CONTINUE
  230    CONTINUE
         CALL dgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ),
     $               ldu,
     $               rwork( 1+nlp1*nrhs*2 ), nlp1, zero,
     $               rwork( 1+nlp1*nrhs ), nlp1 )
         jreal = 0
         jimag = nlp1*nrhs
         DO 250 jcol = 1, nrhs
            DO 240 jrow = nlf, nlf + nlp1 - 1
               jreal = jreal + 1
               jimag = jimag + 1
               bx( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                            rwork( jimag ) )
  240       CONTINUE
  250    CONTINUE
*
*        Since B and BX are complex, the following call to DGEMM is
*        performed in two steps (real and imaginary parts).
*
*        CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
*    $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
*
         j = nrp1*nrhs*2
         DO 270 jcol = 1, nrhs
            DO 260 jrow = nrf, nrf + nrp1 - 1
               j = j + 1
               rwork( j ) = dble( b( jrow, jcol ) )
  260       CONTINUE
  270    CONTINUE
         CALL dgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ),
     $               ldu,
     $               rwork( 1+nrp1*nrhs*2 ), nrp1, zero, rwork( 1 ),
     $               nrp1 )
         j = nrp1*nrhs*2
         DO 290 jcol = 1, nrhs
            DO 280 jrow = nrf, nrf + nrp1 - 1
               j = j + 1
               rwork( j ) = dimag( b( jrow, jcol ) )
  280       CONTINUE
  290    CONTINUE
         CALL dgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ),
     $               ldu,
     $               rwork( 1+nrp1*nrhs*2 ), nrp1, zero,
     $               rwork( 1+nrp1*nrhs ), nrp1 )
         jreal = 0
         jimag = nrp1*nrhs
         DO 310 jcol = 1, nrhs
            DO 300 jrow = nrf, nrf + nrp1 - 1
               jreal = jreal + 1
               jimag = jimag + 1
               bx( jrow, jcol ) = dcmplx( rwork( jreal ),
     $                            rwork( jimag ) )
  300       CONTINUE
  310    CONTINUE
*
  320 CONTINUE
*
  330 CONTINUE
*
      RETURN
*
*     End of ZLALSA
*

Here is the call graph for this function:

Here is the caller graph for this function: