Annotation of rpl/lapack/lapack/zlatsqr.f, revision 1.5

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

CVSweb interface <joel.bertrand@systella.fr>