File:  [local] / rpl / lapack / lapack / dlaswlq.f
Revision 1.2: download - view: text, annotated - select for diffs - revision graph
Sat Jun 17 11:06:27 2017 UTC (6 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_27, rpl-4_1_26, HEAD
Cohérence.

    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.
   58: *>          On exit, the elements on and bleow the diagonal
   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: *
  153: *  -- LAPACK computational routine (version 3.7.0) --
  154: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  155: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
  156: *     December 2016
  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>