zlalsa.f

Go to the documentation of this file.

*> \brief \b ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*> \htmlonly
*> Download ZLALSA + dependencies 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsa.f"> 
*> [TGZ]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsa.f"> 
*> [ZIP]</a> 
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsa.f"> 
*> [TXT]</a>
*> \endhtmlonly 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLALSA( 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 )
* 
*       .. 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, * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZLALSA is an itermediate 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 orthorgonal
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ICOMPQ
*> \verbatim
*>          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
*> \endverbatim
*>
*> \param[in] SMLSIZ
*> \verbatim
*>          SMLSIZ is INTEGER
*>         The maximum size of the subproblems at the bottom of the
*>         computation tree.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>         The row and column dimensions of the upper bidiagonal matrix.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>         The number of columns of B and BX. NRHS must be at least 1.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>         The leading dimension of B in the calling subprogram.
*>         LDB must be at least max(1,MAX( M, N ) ).
*> \endverbatim
*>
*> \param[out] BX
*> \verbatim
*>          BX is COMPLEX*16 array, dimension ( LDBX, NRHS )
*>         On exit, the result of applying the left or right singular
*>         vector matrix to B.
*> \endverbatim
*>
*> \param[in] LDBX
*> \verbatim
*>          LDBX is INTEGER
*>         The leading dimension of BX.
*> \endverbatim
*>
*> \param[in] U
*> \verbatim
*>          U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
*>         On entry, U contains the left singular vector matrices of all
*>         subproblems at the bottom level.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER, LDU = > N.
*>         The leading dimension of arrays U, VT, DIFL, DIFR,
*>         POLES, GIVNUM, and Z.
*> \endverbatim
*>
*> \param[in] VT
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER array, dimension ( N ).
*> \endverbatim
*>
*> \param[in] DIFL
*> \verbatim
*>          DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
*>         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
*> \endverbatim
*>
*> \param[in] DIFR
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] Z
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] POLES
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] GIVPTR
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] GIVCOL
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDGCOL
*> \verbatim
*>          LDGCOL is INTEGER, LDGCOL = > N.
*>         The leading dimension of arrays GIVCOL and PERM.
*> \endverbatim
*>
*> \param[in] PERM
*> \verbatim
*>          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
*>         On entry, PERM(*, I) records permutations done on the I-th
*>         level of the computation tree.
*> \endverbatim
*>
*> \param[in] GIVNUM
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] S
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension at least
*>         MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ).
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array.
*>         The dimension must be at least 3 * N
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date September 2012
*
*> \ingroup complex16OTHERcomputational
*
*> \par Contributors:
*  ==================
*>
*>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
*>       California at Berkeley, USA \n
*>     Osni Marques, LBNL/NERSC, USA \n
*
*  =====================================================================
      SUBROUTINE zlalsa( 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 )
*
*  -- LAPACK computational routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. 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
*
      END