Annotation of rpl/lapack/lapack/dlaswlq.f, revision 1.3

1.1       bertrand    1: *
                      2: *  Definition:
                      3: *  ===========
                      4: *
                      5: *       SUBROUTINE DLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK,
                      6: *                            LWORK, INFO)
                      7: *
                      8: *       .. Scalar Arguments ..
                      9: *       INTEGER           INFO, LDA, M, N, MB, NB, LDT, LWORK
                     10: *       ..
                     11: *       .. Array Arguments ..
                     12: *       DOUBLE PRECISION  A( LDA, * ), T( LDT, * ), WORK( * )
                     13: *       ..
                     14: *
                     15: *
                     16: *> \par Purpose:
                     17: *  =============
                     18: *>
                     19: *> \verbatim
                     20: *>
                     21: *>          DLASWLQ computes a blocked Short-Wide LQ factorization of a
                     22: *>          M-by-N matrix A, where N >= M:
                     23: *>          A = L * Q
                     24: *> \endverbatim
                     25: *
                     26: *  Arguments:
                     27: *  ==========
                     28: *
                     29: *> \param[in] M
                     30: *> \verbatim
                     31: *>          M is INTEGER
                     32: *>          The number of rows of the matrix A.  M >= 0.
                     33: *> \endverbatim
                     34: *>
                     35: *> \param[in] N
                     36: *> \verbatim
                     37: *>          N is INTEGER
                     38: *>          The number of columns of the matrix A.  N >= M >= 0.
                     39: *> \endverbatim
                     40: *>
                     41: *> \param[in] MB
                     42: *> \verbatim
                     43: *>          MB is INTEGER
                     44: *>          The row block size to be used in the blocked QR.
                     45: *>          M >= MB >= 1
                     46: *> \endverbatim
                     47: *> \param[in] NB
                     48: *> \verbatim
                     49: *>          NB is INTEGER
                     50: *>          The column block size to be used in the blocked QR.
                     51: *>          NB > M.
                     52: *> \endverbatim
                     53: *>
                     54: *> \param[in,out] A
                     55: *> \verbatim
                     56: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
                     57: *>          On entry, the M-by-N matrix A.
1.3     ! bertrand   58: *>          On exit, the elements on and below the diagonal
1.1       bertrand   59: *>          of the array contain the N-by-N lower triangular matrix L;
                     60: *>          the elements above the diagonal represent Q by the rows
                     61: *>          of blocked V (see Further Details).
                     62: *>
                     63: *> \endverbatim
                     64: *>
                     65: *> \param[in] LDA
                     66: *> \verbatim
                     67: *>          LDA is INTEGER
                     68: *>          The leading dimension of the array A.  LDA >= max(1,M).
                     69: *> \endverbatim
                     70: *>
                     71: *> \param[out] T
                     72: *> \verbatim
                     73: *>          T is DOUBLE PRECISION array,
                     74: *>          dimension (LDT, N * Number_of_row_blocks)
                     75: *>          where Number_of_row_blocks = CEIL((N-M)/(NB-M))
                     76: *>          The blocked upper triangular block reflectors stored in compact form
                     77: *>          as a sequence of upper triangular blocks.
                     78: *>          See Further Details below.
                     79: *> \endverbatim
                     80: *>
                     81: *> \param[in] LDT
                     82: *> \verbatim
                     83: *>          LDT is INTEGER
                     84: *>          The leading dimension of the array T.  LDT >= MB.
                     85: *> \endverbatim
                     86: *>
                     87: *>
                     88: *> \param[out] WORK
                     89: *> \verbatim
                     90: *>         (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                     91: *>
                     92: *> \endverbatim
                     93: *> \param[in] LWORK
                     94: *> \verbatim
                     95: *>          The dimension of the array WORK.  LWORK >= MB*M.
                     96: *>          If LWORK = -1, then a workspace query is assumed; the routine
                     97: *>          only calculates the optimal size of the WORK array, returns
                     98: *>          this value as the first entry of the WORK array, and no error
                     99: *>          message related to LWORK is issued by XERBLA.
                    100: *>
                    101: *> \endverbatim
                    102: *> \param[out] INFO
                    103: *> \verbatim
                    104: *>          INFO is INTEGER
                    105: *>          = 0:  successful exit
                    106: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                    107: *> \endverbatim
                    108: *
                    109: *  Authors:
                    110: *  ========
                    111: *
                    112: *> \author Univ. of Tennessee
                    113: *> \author Univ. of California Berkeley
                    114: *> \author Univ. of Colorado Denver
                    115: *> \author NAG Ltd.
                    116: *
                    117: *> \par Further Details:
                    118: *  =====================
                    119: *>
                    120: *> \verbatim
                    121: *> Short-Wide LQ (SWLQ) performs LQ by a sequence of orthogonal transformations,
                    122: *> representing Q as a product of other orthogonal matrices
                    123: *>   Q = Q(1) * Q(2) * . . . * Q(k)
                    124: *> where each Q(i) zeros out upper diagonal entries of a block of NB rows of A:
                    125: *>   Q(1) zeros out the upper diagonal entries of rows 1:NB of A
                    126: *>   Q(2) zeros out the bottom MB-N rows of rows [1:M,NB+1:2*NB-M] of A
                    127: *>   Q(3) zeros out the bottom MB-N rows of rows [1:M,2*NB-M+1:3*NB-2*M] of A
                    128: *>   . . .
                    129: *>
                    130: *> Q(1) is computed by GELQT, which represents Q(1) by Householder vectors
                    131: *> stored under the diagonal of rows 1:MB of A, and by upper triangular
                    132: *> block reflectors, stored in array T(1:LDT,1:N).
                    133: *> For more information see Further Details in GELQT.
                    134: *>
                    135: *> Q(i) for i>1 is computed by TPLQT, which represents Q(i) by Householder vectors
                    136: *> stored in columns [(i-1)*(NB-M)+M+1:i*(NB-M)+M] of A, and by upper triangular
                    137: *> block reflectors, stored in array T(1:LDT,(i-1)*M+1:i*M).
                    138: *> The last Q(k) may use fewer rows.
                    139: *> For more information see Further Details in TPQRT.
                    140: *>
                    141: *> For more details of the overall algorithm, see the description of
                    142: *> Sequential TSQR in Section 2.2 of [1].
                    143: *>
                    144: *> [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,”
                    145: *>     J. Demmel, L. Grigori, M. Hoemmen, J. Langou,
                    146: *>     SIAM J. Sci. Comput, vol. 34, no. 1, 2012
                    147: *> \endverbatim
                    148: *>
                    149: *  =====================================================================
                    150:       SUBROUTINE DLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
                    151:      $                  INFO)
                    152: *
1.3     ! bertrand  153: *  -- LAPACK computational routine (version 3.7.1) --
1.1       bertrand  154: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    155: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
1.3     ! bertrand  156: *     June 2017
1.1       bertrand  157: *
                    158: *     .. Scalar Arguments ..
                    159:       INTEGER           INFO, LDA, M, N, MB, NB, LWORK, LDT
                    160: *     ..
                    161: *     .. Array Arguments ..
                    162:       DOUBLE PRECISION  A( LDA, * ), WORK( * ), T( LDT, *)
                    163: *     ..
                    164: *
                    165: *  =====================================================================
                    166: *
                    167: *     ..
                    168: *     .. Local Scalars ..
                    169:       LOGICAL    LQUERY
                    170:       INTEGER    I, II, KK, CTR
                    171: *     ..
                    172: *     .. EXTERNAL FUNCTIONS ..
                    173:       LOGICAL            LSAME
                    174:       EXTERNAL           LSAME
                    175: *     .. EXTERNAL SUBROUTINES ..
                    176:       EXTERNAL           DGELQT, DTPLQT, XERBLA
                    177: *     .. INTRINSIC FUNCTIONS ..
                    178:       INTRINSIC          MAX, MIN, MOD
                    179: *     ..
                    180: *     .. EXECUTABLE STATEMENTS ..
                    181: *
                    182: *     TEST THE INPUT ARGUMENTS
                    183: *
                    184:       INFO = 0
                    185: *
                    186:       LQUERY = ( LWORK.EQ.-1 )
                    187: *
                    188:       IF( M.LT.0 ) THEN
                    189:         INFO = -1
                    190:       ELSE IF( N.LT.0 .OR. N.LT.M ) THEN
                    191:         INFO = -2
                    192:       ELSE IF( MB.LT.1 .OR. ( MB.GT.M .AND. M.GT.0 )) THEN
                    193:         INFO = -3
                    194:       ELSE IF( NB.LE.M ) THEN
                    195:         INFO = -4
                    196:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    197:         INFO = -5
                    198:       ELSE IF( LDT.LT.MB ) THEN
                    199:         INFO = -8
                    200:       ELSE IF( ( LWORK.LT.M*MB) .AND. (.NOT.LQUERY) ) THEN
                    201:         INFO = -10
                    202:       END IF
                    203:       IF( INFO.EQ.0)  THEN
                    204:       WORK(1) = MB*M
                    205:       END IF
                    206: *
                    207:       IF( INFO.NE.0 ) THEN
                    208:         CALL XERBLA( 'DLASWLQ', -INFO )
                    209:         RETURN
                    210:       ELSE IF (LQUERY) THEN
                    211:        RETURN
                    212:       END IF
                    213: *
                    214: *     Quick return if possible
                    215: *
                    216:       IF( MIN(M,N).EQ.0 ) THEN
                    217:           RETURN
                    218:       END IF
                    219: *
                    220: *     The LQ Decomposition
                    221: *
                    222:        IF((M.GE.N).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
                    223:         CALL DGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO)
                    224:         RETURN
                    225:        END IF
                    226: *
                    227:        KK = MOD((N-M),(NB-M))
                    228:        II=N-KK+1
                    229: *
                    230: *      Compute the LQ factorization of the first block A(1:M,1:NB)
                    231: *
                    232:        CALL DGELQT( M, NB, MB, A(1,1), LDA, T, LDT, WORK, INFO)
                    233:        CTR = 1
                    234: *
                    235:        DO I = NB+1, II-NB+M , (NB-M)
                    236: *
                    237: *      Compute the QR factorization of the current block A(1:M,I:I+NB-M)
                    238: *
                    239:          CALL DTPLQT( M, NB-M, 0, MB, A(1,1), LDA, A( 1, I ),
                    240:      $                  LDA, T(1, CTR * M + 1),
                    241:      $                  LDT, WORK, INFO )
                    242:          CTR = CTR + 1
                    243:        END DO
                    244: *
                    245: *     Compute the QR factorization of the last block A(1:M,II:N)
                    246: *
                    247:        IF (II.LE.N) THEN
                    248:         CALL DTPLQT( M, KK, 0, MB, A(1,1), LDA, A( 1, II ),
                    249:      $                  LDA, T(1, CTR * M + 1), LDT,
                    250:      $                  WORK, INFO )
                    251:        END IF
                    252: *
                    253:       WORK( 1 ) = M * MB
                    254:       RETURN
                    255: *
                    256: *     End of DLASWLQ
                    257: *
                    258:       END

CVSweb interface <joel.bertrand@systella.fr>