Annotation of rpl/lapack/lapack/dorbdb4.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DORBDB4
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DORBDB4 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorbdb4.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorbdb4.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorbdb4.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DORBDB4( 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: * DOUBLE PRECISION PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
! 31: * $ WORK(*), X11(LDX11,*), X21(LDX21,*)
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: *> =============
! 37: *>
! 38: *>\verbatim
! 39: *>
! 40: *> DORBDB4 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 DORBDB1, DORBDB2, and DORBDB3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational
! 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 DORCSD for details.
! 199: *>
! 200: *> P1, P2, and Q1 are represented as products of elementary reflectors.
! 201: *> See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR
! 202: *> and DORGLQ.
! 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 DORBDB4( 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: DOUBLE PRECISION PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
! 227: $ WORK(*), X11(LDX11,*), X21(LDX21,*)
! 228: * ..
! 229: *
! 230: * ====================================================================
! 231: *
! 232: * .. Parameters ..
! 233: DOUBLE PRECISION NEGONE, ONE, ZERO
! 234: PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
! 235: * ..
! 236: * .. Local Scalars ..
! 237: DOUBLE PRECISION C, S
! 238: INTEGER CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
! 239: $ LORBDB5, LWORKMIN, LWORKOPT
! 240: LOGICAL LQUERY
! 241: * ..
! 242: * .. External Subroutines ..
! 243: EXTERNAL DLARF, DLARFGP, DORBDB5, DROT, DSCAL, XERBLA
! 244: * ..
! 245: * .. External Functions ..
! 246: DOUBLE PRECISION DNRM2
! 247: EXTERNAL DNRM2
! 248: * ..
! 249: * .. Intrinsic Function ..
! 250: INTRINSIC ATAN2, COS, MAX, SIN, SQRT
! 251: * ..
! 252: * .. Executable Statements ..
! 253: *
! 254: * Test input arguments
! 255: *
! 256: INFO = 0
! 257: LQUERY = LWORK .EQ. -1
! 258: *
! 259: IF( M .LT. 0 ) THEN
! 260: INFO = -1
! 261: ELSE IF( P .LT. M-Q .OR. M-P .LT. M-Q ) THEN
! 262: INFO = -2
! 263: ELSE IF( Q .LT. M-Q .OR. Q .GT. M ) THEN
! 264: INFO = -3
! 265: ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
! 266: INFO = -5
! 267: ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
! 268: INFO = -7
! 269: END IF
! 270: *
! 271: * Compute workspace
! 272: *
! 273: IF( INFO .EQ. 0 ) THEN
! 274: ILARF = 2
! 275: LLARF = MAX( Q-1, P-1, M-P-1 )
! 276: IORBDB5 = 2
! 277: LORBDB5 = Q
! 278: LWORKOPT = ILARF + LLARF - 1
! 279: LWORKOPT = MAX( LWORKOPT, IORBDB5 + LORBDB5 - 1 )
! 280: LWORKMIN = LWORKOPT
! 281: WORK(1) = LWORKOPT
! 282: IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
! 283: INFO = -14
! 284: END IF
! 285: END IF
! 286: IF( INFO .NE. 0 ) THEN
! 287: CALL XERBLA( 'DORBDB4', -INFO )
! 288: RETURN
! 289: ELSE IF( LQUERY ) THEN
! 290: RETURN
! 291: END IF
! 292: *
! 293: * Reduce columns 1, ..., M-Q of X11 and X21
! 294: *
! 295: DO I = 1, M-Q
! 296: *
! 297: IF( I .EQ. 1 ) THEN
! 298: DO J = 1, M
! 299: PHANTOM(J) = ZERO
! 300: END DO
! 301: CALL DORBDB5( P, M-P, Q, PHANTOM(1), 1, PHANTOM(P+1), 1,
! 302: $ X11, LDX11, X21, LDX21, WORK(IORBDB5),
! 303: $ LORBDB5, CHILDINFO )
! 304: CALL DSCAL( P, NEGONE, PHANTOM(1), 1 )
! 305: CALL DLARFGP( P, PHANTOM(1), PHANTOM(2), 1, TAUP1(1) )
! 306: CALL DLARFGP( M-P, PHANTOM(P+1), PHANTOM(P+2), 1, TAUP2(1) )
! 307: THETA(I) = ATAN2( PHANTOM(1), PHANTOM(P+1) )
! 308: C = COS( THETA(I) )
! 309: S = SIN( THETA(I) )
! 310: PHANTOM(1) = ONE
! 311: PHANTOM(P+1) = ONE
! 312: CALL DLARF( 'L', P, Q, PHANTOM(1), 1, TAUP1(1), X11, LDX11,
! 313: $ WORK(ILARF) )
! 314: CALL DLARF( 'L', M-P, Q, PHANTOM(P+1), 1, TAUP2(1), X21,
! 315: $ LDX21, WORK(ILARF) )
! 316: ELSE
! 317: CALL DORBDB5( P-I+1, M-P-I+1, Q-I+1, X11(I,I-1), 1,
! 318: $ X21(I,I-1), 1, X11(I,I), LDX11, X21(I,I),
! 319: $ LDX21, WORK(IORBDB5), LORBDB5, CHILDINFO )
! 320: CALL DSCAL( P-I+1, NEGONE, X11(I,I-1), 1 )
! 321: CALL DLARFGP( P-I+1, X11(I,I-1), X11(I+1,I-1), 1, TAUP1(I) )
! 322: CALL DLARFGP( M-P-I+1, X21(I,I-1), X21(I+1,I-1), 1,
! 323: $ TAUP2(I) )
! 324: THETA(I) = ATAN2( X11(I,I-1), X21(I,I-1) )
! 325: C = COS( THETA(I) )
! 326: S = SIN( THETA(I) )
! 327: X11(I,I-1) = ONE
! 328: X21(I,I-1) = ONE
! 329: CALL DLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, TAUP1(I),
! 330: $ X11(I,I), LDX11, WORK(ILARF) )
! 331: CALL DLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, TAUP2(I),
! 332: $ X21(I,I), LDX21, WORK(ILARF) )
! 333: END IF
! 334: *
! 335: CALL DROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C )
! 336: CALL DLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) )
! 337: C = X21(I,I)
! 338: X21(I,I) = ONE
! 339: CALL DLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
! 340: $ X11(I+1,I), LDX11, WORK(ILARF) )
! 341: CALL DLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
! 342: $ X21(I+1,I), LDX21, WORK(ILARF) )
! 343: IF( I .LT. M-Q ) THEN
! 344: S = SQRT( DNRM2( P-I, X11(I+1,I), 1, X11(I+1,I),
! 345: $ 1 )**2 + DNRM2( M-P-I, X21(I+1,I), 1, X21(I+1,I),
! 346: $ 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 DLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) )
! 356: X11(I,I) = ONE
! 357: CALL DLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
! 358: $ X11(I+1,I), LDX11, WORK(ILARF) )
! 359: CALL DLARF( '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 DLARFGP( 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 DLARF( '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 DORBDB4
! 376: *
! 377: END
! 378:
CVSweb interface <joel.bertrand@systella.fr>