SRC/dlaswlq.f

   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.5.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 *     November 2013
 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