Annotation of rpl/lapack/lapack/dorgtsqr_row.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DORGTSQR_ROW
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DORGTSQR_ROW + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorgtsqr_row.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorgtsqr_row.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorgtsqr_row.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DORGTSQR_ROW( M, N, MB, NB, A, LDA, T, LDT, WORK,
! 22: * $ LWORK, INFO )
! 23: * IMPLICIT NONE
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * INTEGER INFO, LDA, LDT, LWORK, M, N, MB, NB
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *> \par Purpose:
! 33: * =============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *> DORGTSQR_ROW generates an M-by-N real matrix Q_out with
! 38: *> orthonormal columns from the output of DLATSQR. These N orthonormal
! 39: *> columns are the first N columns of a product of complex unitary
! 40: *> matrices Q(k)_in of order M, which are returned by DLATSQR in
! 41: *> a special format.
! 42: *>
! 43: *> Q_out = first_N_columns_of( Q(1)_in * Q(2)_in * ... * Q(k)_in ).
! 44: *>
! 45: *> The input matrices Q(k)_in are stored in row and column blocks in A.
! 46: *> See the documentation of DLATSQR for more details on the format of
! 47: *> Q(k)_in, where each Q(k)_in is represented by block Householder
! 48: *> transformations. This routine calls an auxiliary routine DLARFB_GETT,
! 49: *> where the computation is performed on each individual block. The
! 50: *> algorithm first sweeps NB-sized column blocks from the right to left
! 51: *> starting in the bottom row block and continues to the top row block
! 52: *> (hence _ROW in the routine name). This sweep is in reverse order of
! 53: *> the order in which DLATSQR generates the output blocks.
! 54: *> \endverbatim
! 55: *
! 56: * Arguments:
! 57: * ==========
! 58: *
! 59: *> \param[in] M
! 60: *> \verbatim
! 61: *> M is INTEGER
! 62: *> The number of rows of the matrix A. M >= 0.
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in] N
! 66: *> \verbatim
! 67: *> N is INTEGER
! 68: *> The number of columns of the matrix A. M >= N >= 0.
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] MB
! 72: *> \verbatim
! 73: *> MB is INTEGER
! 74: *> The row block size used by DLATSQR to return
! 75: *> arrays A and T. MB > N.
! 76: *> (Note that if MB > M, then M is used instead of MB
! 77: *> as the row block size).
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] NB
! 81: *> \verbatim
! 82: *> NB is INTEGER
! 83: *> The column block size used by DLATSQR to return
! 84: *> arrays A and T. NB >= 1.
! 85: *> (Note that if NB > N, then N is used instead of NB
! 86: *> as the column block size).
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in,out] A
! 90: *> \verbatim
! 91: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 92: *>
! 93: *> On entry:
! 94: *>
! 95: *> The elements on and above the diagonal are not used as
! 96: *> input. The elements below the diagonal represent the unit
! 97: *> lower-trapezoidal blocked matrix V computed by DLATSQR
! 98: *> that defines the input matrices Q_in(k) (ones on the
! 99: *> diagonal are not stored). See DLATSQR for more details.
! 100: *>
! 101: *> On exit:
! 102: *>
! 103: *> The array A contains an M-by-N orthonormal matrix Q_out,
! 104: *> i.e the columns of A are orthogonal unit vectors.
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[in] LDA
! 108: *> \verbatim
! 109: *> LDA is INTEGER
! 110: *> The leading dimension of the array A. LDA >= max(1,M).
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] T
! 114: *> \verbatim
! 115: *> T is DOUBLE PRECISION array,
! 116: *> dimension (LDT, N * NIRB)
! 117: *> where NIRB = Number_of_input_row_blocks
! 118: *> = MAX( 1, CEIL((M-N)/(MB-N)) )
! 119: *> Let NICB = Number_of_input_col_blocks
! 120: *> = CEIL(N/NB)
! 121: *>
! 122: *> The upper-triangular block reflectors used to define the
! 123: *> input matrices Q_in(k), k=(1:NIRB*NICB). The block
! 124: *> reflectors are stored in compact form in NIRB block
! 125: *> reflector sequences. Each of the NIRB block reflector
! 126: *> sequences is stored in a larger NB-by-N column block of T
! 127: *> and consists of NICB smaller NB-by-NB upper-triangular
! 128: *> column blocks. See DLATSQR for more details on the format
! 129: *> of T.
! 130: *> \endverbatim
! 131: *>
! 132: *> \param[in] LDT
! 133: *> \verbatim
! 134: *> LDT is INTEGER
! 135: *> The leading dimension of the array T.
! 136: *> LDT >= max(1,min(NB,N)).
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[out] WORK
! 140: *> \verbatim
! 141: *> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 142: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[in] LWORK
! 146: *> \verbatim
! 147: *> The dimension of the array WORK.
! 148: *> LWORK >= NBLOCAL * MAX(NBLOCAL,(N-NBLOCAL)),
! 149: *> where NBLOCAL=MIN(NB,N).
! 150: *> If LWORK = -1, then a workspace query is assumed.
! 151: *> The routine only calculates the optimal size of the WORK
! 152: *> array, returns this value as the first entry of the WORK
! 153: *> array, and no error message related to LWORK is issued
! 154: *> by XERBLA.
! 155: *> \endverbatim
! 156: *>
! 157: *> \param[out] INFO
! 158: *> \verbatim
! 159: *> INFO is INTEGER
! 160: *> = 0: successful exit
! 161: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 162: *> \endverbatim
! 163: *>
! 164: * Authors:
! 165: * ========
! 166: *
! 167: *> \author Univ. of Tennessee
! 168: *> \author Univ. of California Berkeley
! 169: *> \author Univ. of Colorado Denver
! 170: *> \author NAG Ltd.
! 171: *
! 172: *> \ingroup doubleOTHERcomputational
! 173: *
! 174: *> \par Contributors:
! 175: * ==================
! 176: *>
! 177: *> \verbatim
! 178: *>
! 179: *> November 2020, Igor Kozachenko,
! 180: *> Computer Science Division,
! 181: *> University of California, Berkeley
! 182: *>
! 183: *> \endverbatim
! 184: *>
! 185: * =====================================================================
! 186: SUBROUTINE DORGTSQR_ROW( M, N, MB, NB, A, LDA, T, LDT, WORK,
! 187: $ LWORK, INFO )
! 188: IMPLICIT NONE
! 189: *
! 190: * -- LAPACK computational routine --
! 191: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 192: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 193: *
! 194: * .. Scalar Arguments ..
! 195: INTEGER INFO, LDA, LDT, LWORK, M, N, MB, NB
! 196: * ..
! 197: * .. Array Arguments ..
! 198: DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
! 199: * ..
! 200: *
! 201: * =====================================================================
! 202: *
! 203: * .. Parameters ..
! 204: DOUBLE PRECISION ONE, ZERO
! 205: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 206: * ..
! 207: * .. Local Scalars ..
! 208: LOGICAL LQUERY
! 209: INTEGER NBLOCAL, MB2, M_PLUS_ONE, ITMP, IB_BOTTOM,
! 210: $ LWORKOPT, NUM_ALL_ROW_BLOCKS, JB_T, IB, IMB,
! 211: $ KB, KB_LAST, KNB, MB1
! 212: * ..
! 213: * .. Local Arrays ..
! 214: DOUBLE PRECISION DUMMY( 1, 1 )
! 215: * ..
! 216: * .. External Subroutines ..
! 217: EXTERNAL DLARFB_GETT, DLASET, XERBLA
! 218: * ..
! 219: * .. Intrinsic Functions ..
! 220: INTRINSIC DBLE, MAX, MIN
! 221: * ..
! 222: * .. Executable Statements ..
! 223: *
! 224: * Test the input parameters
! 225: *
! 226: INFO = 0
! 227: LQUERY = LWORK.EQ.-1
! 228: IF( M.LT.0 ) THEN
! 229: INFO = -1
! 230: ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
! 231: INFO = -2
! 232: ELSE IF( MB.LE.N ) THEN
! 233: INFO = -3
! 234: ELSE IF( NB.LT.1 ) THEN
! 235: INFO = -4
! 236: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 237: INFO = -6
! 238: ELSE IF( LDT.LT.MAX( 1, MIN( NB, N ) ) ) THEN
! 239: INFO = -8
! 240: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
! 241: INFO = -10
! 242: END IF
! 243: *
! 244: NBLOCAL = MIN( NB, N )
! 245: *
! 246: * Determine the workspace size.
! 247: *
! 248: IF( INFO.EQ.0 ) THEN
! 249: LWORKOPT = NBLOCAL * MAX( NBLOCAL, ( N - NBLOCAL ) )
! 250: END IF
! 251: *
! 252: * Handle error in the input parameters and handle the workspace query.
! 253: *
! 254: IF( INFO.NE.0 ) THEN
! 255: CALL XERBLA( 'DORGTSQR_ROW', -INFO )
! 256: RETURN
! 257: ELSE IF ( LQUERY ) THEN
! 258: WORK( 1 ) = DBLE( LWORKOPT )
! 259: RETURN
! 260: END IF
! 261: *
! 262: * Quick return if possible
! 263: *
! 264: IF( MIN( M, N ).EQ.0 ) THEN
! 265: WORK( 1 ) = DBLE( LWORKOPT )
! 266: RETURN
! 267: END IF
! 268: *
! 269: * (0) Set the upper-triangular part of the matrix A to zero and
! 270: * its diagonal elements to one.
! 271: *
! 272: CALL DLASET('U', M, N, ZERO, ONE, A, LDA )
! 273: *
! 274: * KB_LAST is the column index of the last column block reflector
! 275: * in the matrices T and V.
! 276: *
! 277: KB_LAST = ( ( N-1 ) / NBLOCAL ) * NBLOCAL + 1
! 278: *
! 279: *
! 280: * (1) Bottom-up loop over row blocks of A, except the top row block.
! 281: * NOTE: If MB>=M, then the loop is never executed.
! 282: *
! 283: IF ( MB.LT.M ) THEN
! 284: *
! 285: * MB2 is the row blocking size for the row blocks before the
! 286: * first top row block in the matrix A. IB is the row index for
! 287: * the row blocks in the matrix A before the first top row block.
! 288: * IB_BOTTOM is the row index for the last bottom row block
! 289: * in the matrix A. JB_T is the column index of the corresponding
! 290: * column block in the matrix T.
! 291: *
! 292: * Initialize variables.
! 293: *
! 294: * NUM_ALL_ROW_BLOCKS is the number of row blocks in the matrix A
! 295: * including the first row block.
! 296: *
! 297: MB2 = MB - N
! 298: M_PLUS_ONE = M + 1
! 299: ITMP = ( M - MB - 1 ) / MB2
! 300: IB_BOTTOM = ITMP * MB2 + MB + 1
! 301: NUM_ALL_ROW_BLOCKS = ITMP + 2
! 302: JB_T = NUM_ALL_ROW_BLOCKS * N + 1
! 303: *
! 304: DO IB = IB_BOTTOM, MB+1, -MB2
! 305: *
! 306: * Determine the block size IMB for the current row block
! 307: * in the matrix A.
! 308: *
! 309: IMB = MIN( M_PLUS_ONE - IB, MB2 )
! 310: *
! 311: * Determine the column index JB_T for the current column block
! 312: * in the matrix T.
! 313: *
! 314: JB_T = JB_T - N
! 315: *
! 316: * Apply column blocks of H in the row block from right to left.
! 317: *
! 318: * KB is the column index of the current column block reflector
! 319: * in the matrices T and V.
! 320: *
! 321: DO KB = KB_LAST, 1, -NBLOCAL
! 322: *
! 323: * Determine the size of the current column block KNB in
! 324: * the matrices T and V.
! 325: *
! 326: KNB = MIN( NBLOCAL, N - KB + 1 )
! 327: *
! 328: CALL DLARFB_GETT( 'I', IMB, N-KB+1, KNB,
! 329: $ T( 1, JB_T+KB-1 ), LDT, A( KB, KB ), LDA,
! 330: $ A( IB, KB ), LDA, WORK, KNB )
! 331: *
! 332: END DO
! 333: *
! 334: END DO
! 335: *
! 336: END IF
! 337: *
! 338: * (2) Top row block of A.
! 339: * NOTE: If MB>=M, then we have only one row block of A of size M
! 340: * and we work on the entire matrix A.
! 341: *
! 342: MB1 = MIN( MB, M )
! 343: *
! 344: * Apply column blocks of H in the top row block from right to left.
! 345: *
! 346: * KB is the column index of the current block reflector in
! 347: * the matrices T and V.
! 348: *
! 349: DO KB = KB_LAST, 1, -NBLOCAL
! 350: *
! 351: * Determine the size of the current column block KNB in
! 352: * the matrices T and V.
! 353: *
! 354: KNB = MIN( NBLOCAL, N - KB + 1 )
! 355: *
! 356: IF( MB1-KB-KNB+1.EQ.0 ) THEN
! 357: *
! 358: * In SLARFB_GETT parameters, when M=0, then the matrix B
! 359: * does not exist, hence we need to pass a dummy array
! 360: * reference DUMMY(1,1) to B with LDDUMMY=1.
! 361: *
! 362: CALL DLARFB_GETT( 'N', 0, N-KB+1, KNB,
! 363: $ T( 1, KB ), LDT, A( KB, KB ), LDA,
! 364: $ DUMMY( 1, 1 ), 1, WORK, KNB )
! 365: ELSE
! 366: CALL DLARFB_GETT( 'N', MB1-KB-KNB+1, N-KB+1, KNB,
! 367: $ T( 1, KB ), LDT, A( KB, KB ), LDA,
! 368: $ A( KB+KNB, KB), LDA, WORK, KNB )
! 369:
! 370: END IF
! 371: *
! 372: END DO
! 373: *
! 374: WORK( 1 ) = DBLE( LWORKOPT )
! 375: RETURN
! 376: *
! 377: * End of DORGTSQR_ROW
! 378: *
! 379: END
CVSweb interface <joel.bertrand@systella.fr>