SRC/sorbdb4.f

   1 *> \brief \b SORBDB4
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SORBDB4 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorbdb4.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorbdb4.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorbdb4.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SORBDB4( 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 *       REAL               PHI(*), THETA(*)
  30 *       REAL               PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
  31 *      $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *> =============
  37 *>
  38 *>\verbatim
  39 *>
  40 *> SORBDB4 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 SORBDB1, SORBDB2, and SORBDB3 handle cases in
  51 *> which M-Q is not the minimum dimension.
  52 *>
  53 *> The orthogonal 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 *
 179 *  Authors:
 180 *  ========
 181 *
 182 *> \author Univ. of Tennessee
 183 *> \author Univ. of California Berkeley
 184 *> \author Univ. of Colorado Denver
 185 *> \author NAG Ltd.
 186 *
 187 *> \date July 2012
 188 *
 189 *> \ingroup realOTHERcomputational
 190 *
 191 *> \par Further Details:
 192 *  =====================
 193 *>
 194 *> \verbatim
 195 *>
 196 *>  The upper-bidiagonal blocks B11, B21 are represented implicitly by
 197 *>  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
 198 *>  in each bidiagonal band is a product of a sine or cosine of a THETA
 199 *>  with a sine or cosine of a PHI. See [1] or SORCSD for details.
 200 *>
 201 *>  P1, P2, and Q1 are represented as products of elementary reflectors.
 202 *>  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
 203 *>  and SORGLQ.
 204 *> \endverbatim
 205 *
 206 *> \par References:
 207 *  ================
 208 *>
 209 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
 210 *>      Algorithms, 50(1):33-65, 2009.
 211 *>
 212 *  =====================================================================
 213       SUBROUTINE SORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
 214      $                    TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
 215      $                    INFO )
 216 *
 217 *  -- LAPACK computational routine (version 3.6.1) --
 218 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 219 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 220 *     July 2012
 221 *
 222 *     .. Scalar Arguments ..
 223       INTEGER            INFO, LWORK, M, P, Q, LDX11, LDX21
 224 *     ..
 225 *     .. Array Arguments ..
 226       REAL               PHI(*), THETA(*)
 227       REAL               PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
 228      $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
 229 *     ..
 230 *
 231 *  ====================================================================
 232 *
 233 *     .. Parameters ..
 234       REAL               NEGONE, ONE, ZERO
 235       PARAMETER          ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 )
 236 *     ..
 237 *     .. Local Scalars ..
 238       REAL               C, S
 239       INTEGER            CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
 240      $                   LORBDB5, LWORKMIN, LWORKOPT
 241       LOGICAL            LQUERY
 242 *     ..
 243 *     .. External Subroutines ..
 244       EXTERNAL           SLARF, SLARFGP, SORBDB5, SROT, SSCAL, XERBLA
 245 *     ..
 246 *     .. External Functions ..
 247       REAL               SNRM2
 248       EXTERNAL           SNRM2
 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( 'SORBDB4', -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 SORBDB5( P, M-P, Q, PHANTOM(1), 1, PHANTOM(P+1), 1,
 303      $                    X11, LDX11, X21, LDX21, WORK(IORBDB5),
 304      $                    LORBDB5, CHILDINFO )
 305             CALL SSCAL( P, NEGONE, PHANTOM(1), 1 )
 306             CALL SLARFGP( P, PHANTOM(1), PHANTOM(2), 1, TAUP1(1) )
 307             CALL SLARFGP( M-P, PHANTOM(P+1), PHANTOM(P+2), 1, TAUP2(1) )
 308             THETA(I) = ATAN2( PHANTOM(1), PHANTOM(P+1) )
 309             C = COS( THETA(I) )
 310             S = SIN( THETA(I) )
 311             PHANTOM(1) = ONE
 312             PHANTOM(P+1) = ONE
 313             CALL SLARF( 'L', P, Q, PHANTOM(1), 1, TAUP1(1), X11, LDX11,
 314      $                  WORK(ILARF) )
 315             CALL SLARF( 'L', M-P, Q, PHANTOM(P+1), 1, TAUP2(1), X21,
 316      $                  LDX21, WORK(ILARF) )
 317          ELSE
 318             CALL SORBDB5( 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 SSCAL( P-I+1, NEGONE, X11(I,I-1), 1 )
 322             CALL SLARFGP( P-I+1, X11(I,I-1), X11(I+1,I-1), 1, TAUP1(I) )
 323             CALL SLARFGP( M-P-I+1, X21(I,I-1), X21(I+1,I-1), 1,
 324      $                    TAUP2(I) )
 325             THETA(I) = ATAN2( X11(I,I-1), 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 SLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, TAUP1(I),
 331      $                  X11(I,I), LDX11, WORK(ILARF) )
 332             CALL SLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, TAUP2(I),
 333      $                  X21(I,I), LDX21, WORK(ILARF) )
 334          END IF
 335 *
 336          CALL SROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C )
 337          CALL SLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) )
 338          C = X21(I,I)
 339          X21(I,I) = ONE
 340          CALL SLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
 341      $               X11(I+1,I), LDX11, WORK(ILARF) )
 342          CALL SLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
 343      $               X21(I+1,I), LDX21, WORK(ILARF) )
 344          IF( I .LT. M-Q ) THEN
 345             S = SQRT( SNRM2( P-I, X11(I+1,I), 1 )**2
 346      $              + SNRM2( M-P-I, X21(I+1,I), 1 )**2 )
 347             PHI(I) = ATAN2( S, C )
 348          END IF
 349 *
 350       END DO
 351 *
 352 *     Reduce the bottom-right portion of X11 to [ I 0 ]
 353 *
 354       DO I = M - Q + 1, P
 355          CALL SLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) )
 356          X11(I,I) = ONE
 357          CALL SLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
 358      $               X11(I+1,I), LDX11, WORK(ILARF) )
 359          CALL SLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
 360      $               X21(M-Q+1,I), LDX21, WORK(ILARF) )
 361       END DO
 362 *
 363 *     Reduce the bottom-right portion of X21 to [ 0 I ]
 364 *
 365       DO I = P + 1, Q
 366          CALL SLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), LDX21,
 367      $                 TAUQ1(I) )
 368          X21(M-Q+I-P,I) = ONE
 369          CALL SLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, TAUQ1(I),
 370      $               X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) )
 371       END DO
 372 *
 373       RETURN
 374 *
 375 *     End of SORBDB4
 376 *
 377       END
 378