SRC/zunbdb4.f

   1 *> \brief \b ZUNBDB4
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZUNBDB4 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb4.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb4.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb4.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
  22 *                           TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
  23 *                           INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       INTEGER            INFO, LWORK, M, P, Q, LDX11, LDX21
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   PHI(*), THETA(*)
  30 *       COMPLEX*16         PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
  31 *      $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *> =============
  37 *>
  38 *>\verbatim
  39 *>
  40 *> ZUNBDB4 simultaneously bidiagonalizes the blocks of a tall and skinny
  41 *> matrix X with orthonomal columns:
  42 *>
  43 *>                            [ B11 ]
  44 *>      [ X11 ]   [ P1 |    ] [  0  ]
  45 *>      [-----] = [---------] [-----] Q1**T .
  46 *>      [ X21 ]   [    | P2 ] [ B21 ]
  47 *>                            [  0  ]
  48 *>
  49 *> X11 is P-by-Q, and X21 is (M-P)-by-Q. M-Q must be no larger than P,
  50 *> M-P, or Q. Routines ZUNBDB1, ZUNBDB2, and ZUNBDB3 handle cases in
  51 *> which M-Q is not the minimum dimension.
  52 *>
  53 *> The unitary matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
  54 *> and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
  55 *> Householder vectors.
  56 *>
  57 *> B11 and B12 are (M-Q)-by-(M-Q) bidiagonal matrices represented
  58 *> implicitly by angles THETA, PHI.
  59 *>
  60 *>\endverbatim
  61 *
  62 *  Arguments:
  63 *  ==========
  64 *
  65 *> \param[in] M
  66 *> \verbatim
  67 *>          M is INTEGER
  68 *>           The number of rows X11 plus the number of rows in X21.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] P
  72 *> \verbatim
  73 *>          P is INTEGER
  74 *>           The number of rows in X11. 0 <= P <= M.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] Q
  78 *> \verbatim
  79 *>          Q is INTEGER
  80 *>           The number of columns in X11 and X21. 0 <= Q <= M and
  81 *>           M-Q <= min(P,M-P,Q).
  82 *> \endverbatim
  83 *>
  84 *> \param[in,out] X11
  85 *> \verbatim
  86 *>          X11 is COMPLEX*16 array, dimension (LDX11,Q)
  87 *>           On entry, the top block of the matrix X to be reduced. On
  88 *>           exit, the columns of tril(X11) specify reflectors for P1 and
  89 *>           the rows of triu(X11,1) specify reflectors for Q1.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] LDX11
  93 *> \verbatim
  94 *>          LDX11 is INTEGER
  95 *>           The leading dimension of X11. LDX11 >= P.
  96 *> \endverbatim
  97 *>
  98 *> \param[in,out] X21
  99 *> \verbatim
 100 *>          X21 is COMPLEX*16 array, dimension (LDX21,Q)
 101 *>           On entry, the bottom block of the matrix X to be reduced. On
 102 *>           exit, the columns of tril(X21) specify reflectors for P2.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] LDX21
 106 *> \verbatim
 107 *>          LDX21 is INTEGER
 108 *>           The leading dimension of X21. LDX21 >= M-P.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] THETA
 112 *> \verbatim
 113 *>          THETA is DOUBLE PRECISION array, dimension (Q)
 114 *>           The entries of the bidiagonal blocks B11, B21 are defined by
 115 *>           THETA and PHI. See Further Details.
 116 *> \endverbatim
 117 *>
 118 *> \param[out] PHI
 119 *> \verbatim
 120 *>          PHI is DOUBLE PRECISION array, dimension (Q-1)
 121 *>           The entries of the bidiagonal blocks B11, B21 are defined by
 122 *>           THETA and PHI. See Further Details.
 123 *> \endverbatim
 124 *>
 125 *> \param[out] TAUP1
 126 *> \verbatim
 127 *>          TAUP1 is COMPLEX*16 array, dimension (P)
 128 *>           The scalar factors of the elementary reflectors that define
 129 *>           P1.
 130 *> \endverbatim
 131 *>
 132 *> \param[out] TAUP2
 133 *> \verbatim
 134 *>          TAUP2 is COMPLEX*16 array, dimension (M-P)
 135 *>           The scalar factors of the elementary reflectors that define
 136 *>           P2.
 137 *> \endverbatim
 138 *>
 139 *> \param[out] TAUQ1
 140 *> \verbatim
 141 *>          TAUQ1 is COMPLEX*16 array, dimension (Q)
 142 *>           The scalar factors of the elementary reflectors that define
 143 *>           Q1.
 144 *> \endverbatim
 145 *>
 146 *> \param[out] PHANTOM
 147 *> \verbatim
 148 *>          PHANTOM is COMPLEX*16 array, dimension (M)
 149 *>           The routine computes an M-by-1 column vector Y that is
 150 *>           orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
 151 *>           PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
 152 *>           Y(P+1:M), respectively.
 153 *> \endverbatim
 154 *>
 155 *> \param[out] WORK
 156 *> \verbatim
 157 *>          WORK is COMPLEX*16 array, dimension (LWORK)
 158 *> \endverbatim
 159 *>
 160 *> \param[in] LWORK
 161 *> \verbatim
 162 *>          LWORK is INTEGER
 163 *>           The dimension of the array WORK. LWORK >= M-Q.
 164 *>
 165 *>           If LWORK = -1, then a workspace query is assumed; the routine
 166 *>           only calculates the optimal size of the WORK array, returns
 167 *>           this value as the first entry of the WORK array, and no error
 168 *>           message related to LWORK is issued by XERBLA.
 169 *> \endverbatim
 170 *>
 171 *> \param[out] INFO
 172 *> \verbatim
 173 *>          INFO is INTEGER
 174 *>           = 0:  successful exit.
 175 *>           < 0:  if INFO = -i, the i-th argument had an illegal value.
 176 *> \endverbatim
 177 *
 178 *  Authors:
 179 *  ========
 180 *
 181 *> \author Univ. of Tennessee
 182 *> \author Univ. of California Berkeley
 183 *> \author Univ. of Colorado Denver
 184 *> \author NAG Ltd.
 185 *
 186 *> \date July 2012
 187 *
 188 *> \ingroup complex16OTHERcomputational
 189 *
 190 *> \par Further Details:
 191 *  =====================
 192 *>
 193 *> \verbatim
 194 *>
 195 *>  The upper-bidiagonal blocks B11, B21 are represented implicitly by
 196 *>  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
 197 *>  in each bidiagonal band is a product of a sine or cosine of a THETA
 198 *>  with a sine or cosine of a PHI. See [1] or ZUNCSD for details.
 199 *>
 200 *>  P1, P2, and Q1 are represented as products of elementary reflectors.
 201 *>  See ZUNCSD2BY1 for details on generating P1, P2, and Q1 using ZUNGQR
 202 *>  and ZUNGLQ.
 203 *> \endverbatim
 204 *
 205 *> \par References:
 206 *  ================
 207 *>
 208 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
 209 *>      Algorithms, 50(1):33-65, 2009.
 210 *>
 211 *  =====================================================================
 212       SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
 213      $                    TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
 214      $                    INFO )
 215 *
 216 *  -- LAPACK computational routine (version 3.6.1) --
 217 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 218 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 219 *     July 2012
 220 *
 221 *     .. Scalar Arguments ..
 222       INTEGER            INFO, LWORK, M, P, Q, LDX11, LDX21
 223 *     ..
 224 *     .. Array Arguments ..
 225       DOUBLE PRECISION   PHI(*), THETA(*)
 226       COMPLEX*16         PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
 227      $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
 228 *     ..
 229 *
 230 *  ====================================================================
 231 *
 232 *     .. Parameters ..
 233       COMPLEX*16         NEGONE, ONE, ZERO
 234       PARAMETER          ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0),
 235      $                     ZERO = (0.0D0,0.0D0) )
 236 *     ..
 237 *     .. Local Scalars ..
 238       DOUBLE PRECISION   C, S
 239       INTEGER            CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
 240      $                   LORBDB5, LWORKMIN, LWORKOPT
 241       LOGICAL            LQUERY
 242 *     ..
 243 *     .. External Subroutines ..
 244       EXTERNAL           ZLARF, ZLARFGP, ZUNBDB5, ZDROT, ZSCAL, XERBLA
 245 *     ..
 246 *     .. External Functions ..
 247       DOUBLE PRECISION   DZNRM2
 248       EXTERNAL           DZNRM2
 249 *     ..
 250 *     .. Intrinsic Function ..
 251       INTRINSIC          ATAN2, COS, MAX, SIN, SQRT
 252 *     ..
 253 *     .. Executable Statements ..
 254 *
 255 *     Test input arguments
 256 *
 257       INFO = 0
 258       LQUERY = LWORK .EQ. -1
 259 *
 260       IF( M .LT. 0 ) THEN
 261          INFO = -1
 262       ELSE IF( P .LT. M-Q .OR. M-P .LT. M-Q ) THEN
 263          INFO = -2
 264       ELSE IF( Q .LT. M-Q .OR. Q .GT. M ) THEN
 265          INFO = -3
 266       ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
 267          INFO = -5
 268       ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
 269          INFO = -7
 270       END IF
 271 *
 272 *     Compute workspace
 273 *
 274       IF( INFO .EQ. 0 ) THEN
 275          ILARF = 2
 276          LLARF = MAX( Q-1, P-1, M-P-1 )
 277          IORBDB5 = 2
 278          LORBDB5 = Q
 279          LWORKOPT = ILARF + LLARF - 1
 280          LWORKOPT = MAX( LWORKOPT, IORBDB5 + LORBDB5 - 1 )
 281          LWORKMIN = LWORKOPT
 282          WORK(1) = LWORKOPT
 283          IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
 284            INFO = -14
 285          END IF
 286       END IF
 287       IF( INFO .NE. 0 ) THEN
 288          CALL XERBLA( 'ZUNBDB4', -INFO )
 289          RETURN
 290       ELSE IF( LQUERY ) THEN
 291          RETURN
 292       END IF
 293 *
 294 *     Reduce columns 1, ..., M-Q of X11 and X21
 295 *
 296       DO I = 1, M-Q
 297 *
 298          IF( I .EQ. 1 ) THEN
 299             DO J = 1, M
 300                PHANTOM(J) = ZERO
 301             END DO
 302             CALL ZUNBDB5( P, M-P, Q, PHANTOM(1), 1, PHANTOM(P+1), 1,
 303      $                    X11, LDX11, X21, LDX21, WORK(IORBDB5),
 304      $                    LORBDB5, CHILDINFO )
 305             CALL ZSCAL( P, NEGONE, PHANTOM(1), 1 )
 306             CALL ZLARFGP( P, PHANTOM(1), PHANTOM(2), 1, TAUP1(1) )
 307             CALL ZLARFGP( M-P, PHANTOM(P+1), PHANTOM(P+2), 1, TAUP2(1) )
 308             THETA(I) = ATAN2( DBLE( PHANTOM(1) ), DBLE( PHANTOM(P+1) ) )
 309             C = COS( THETA(I) )
 310             S = SIN( THETA(I) )
 311             PHANTOM(1) = ONE
 312             PHANTOM(P+1) = ONE
 313             CALL ZLARF( 'L', P, Q, PHANTOM(1), 1, DCONJG(TAUP1(1)), X11,
 314      $                  LDX11, WORK(ILARF) )
 315             CALL ZLARF( 'L', M-P, Q, PHANTOM(P+1), 1, DCONJG(TAUP2(1)),
 316      $                  X21, LDX21, WORK(ILARF) )
 317          ELSE
 318             CALL ZUNBDB5( P-I+1, M-P-I+1, Q-I+1, X11(I,I-1), 1,
 319      $                    X21(I,I-1), 1, X11(I,I), LDX11, X21(I,I),
 320      $                    LDX21, WORK(IORBDB5), LORBDB5, CHILDINFO )
 321             CALL ZSCAL( P-I+1, NEGONE, X11(I,I-1), 1 )
 322             CALL ZLARFGP( P-I+1, X11(I,I-1), X11(I+1,I-1), 1, TAUP1(I) )
 323             CALL ZLARFGP( M-P-I+1, X21(I,I-1), X21(I+1,I-1), 1,
 324      $                    TAUP2(I) )
 325             THETA(I) = ATAN2( DBLE( X11(I,I-1) ), DBLE( X21(I,I-1) ) )
 326             C = COS( THETA(I) )
 327             S = SIN( THETA(I) )
 328             X11(I,I-1) = ONE
 329             X21(I,I-1) = ONE
 330             CALL ZLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1,
 331      $                  DCONJG(TAUP1(I)), X11(I,I), LDX11, WORK(ILARF) )
 332             CALL ZLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1,
 333      $                  DCONJG(TAUP2(I)), X21(I,I), LDX21, WORK(ILARF) )
 334          END IF
 335 *
 336          CALL ZDROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C )
 337          CALL ZLACGV( Q-I+1, X21(I,I), LDX21 )
 338          CALL ZLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) )
 339          C = DBLE( X21(I,I) )
 340          X21(I,I) = ONE
 341          CALL ZLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
 342      $               X11(I+1,I), LDX11, WORK(ILARF) )
 343          CALL ZLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
 344      $               X21(I+1,I), LDX21, WORK(ILARF) )
 345          CALL ZLACGV( Q-I+1, X21(I,I), LDX21 )
 346          IF( I .LT. M-Q ) THEN
 347             S = SQRT( DZNRM2( P-I, X11(I+1,I), 1 )**2
 348      $              + DZNRM2( M-P-I, X21(I+1,I), 1 )**2 )
 349             PHI(I) = ATAN2( S, C )
 350          END IF
 351 *
 352       END DO
 353 *
 354 *     Reduce the bottom-right portion of X11 to [ I 0 ]
 355 *
 356       DO I = M - Q + 1, P
 357          CALL ZLACGV( Q-I+1, X11(I,I), LDX11 )
 358          CALL ZLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) )
 359          X11(I,I) = ONE
 360          CALL ZLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
 361      $               X11(I+1,I), LDX11, WORK(ILARF) )
 362          CALL ZLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
 363      $               X21(M-Q+1,I), LDX21, WORK(ILARF) )
 364          CALL ZLACGV( Q-I+1, X11(I,I), LDX11 )
 365       END DO
 366 *
 367 *     Reduce the bottom-right portion of X21 to [ 0 I ]
 368 *
 369       DO I = P + 1, Q
 370          CALL ZLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
 371          CALL ZLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), LDX21,
 372      $                 TAUQ1(I) )
 373          X21(M-Q+I-P,I) = ONE
 374          CALL ZLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, TAUQ1(I),
 375      $               X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) )
 376          CALL ZLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
 377       END DO
 378 *
 379       RETURN
 380 *
 381 *     End of ZUNBDB4
 382 *
 383       END
 384