Annotation of rpl/lapack/lapack/zunbdb4.f, revision 1.1
1.1 ! bertrand 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.5.0) --
! 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, X11(I+1,I),
! 348: $ 1 )**2 + DZNRM2( M-P-I, X21(I+1,I), 1, X21(I+1,I),
! 349: $ 1 )**2 )
! 350: PHI(I) = ATAN2( S, C )
! 351: END IF
! 352: *
! 353: END DO
! 354: *
! 355: * Reduce the bottom-right portion of X11 to [ I 0 ]
! 356: *
! 357: DO I = M - Q + 1, P
! 358: CALL ZLACGV( Q-I+1, X11(I,I), LDX11 )
! 359: CALL ZLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) )
! 360: X11(I,I) = ONE
! 361: CALL ZLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
! 362: $ X11(I+1,I), LDX11, WORK(ILARF) )
! 363: CALL ZLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
! 364: $ X21(M-Q+1,I), LDX21, WORK(ILARF) )
! 365: CALL ZLACGV( Q-I+1, X11(I,I), LDX11 )
! 366: END DO
! 367: *
! 368: * Reduce the bottom-right portion of X21 to [ 0 I ]
! 369: *
! 370: DO I = P + 1, Q
! 371: CALL ZLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
! 372: CALL ZLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), LDX21,
! 373: $ TAUQ1(I) )
! 374: X21(M-Q+I-P,I) = ONE
! 375: CALL ZLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, TAUQ1(I),
! 376: $ X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) )
! 377: CALL ZLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
! 378: END DO
! 379: *
! 380: RETURN
! 381: *
! 382: * End of ZUNBDB4
! 383: *
! 384: END
! 385:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zunbdb4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack