SRC/slamswlq.f

   1 *
   2 *  Definition:
   3 *  ===========
   4 *
   5 *      SUBROUTINE SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T,
   6 *     $                LDT, C, LDC, WORK, LWORK, INFO )
   7 *
   8 *
   9 *     .. Scalar Arguments ..
  10 *      CHARACTER         SIDE, TRANS
  11 *      INTEGER           INFO, LDA, M, N, K, MB, NB, LDT, LWORK, LDC
  12 *     ..
  13 *     .. Array Arguments ..
  14 *      DOUBLE        A( LDA, * ), WORK( * ), C(LDC, * ),
  15 *     $                  T( LDT, * )
  16 *> \par Purpose:
  17 *  =============
  18 *>
  19 *> \verbatim
  20 *>
  21 *>    SLAMSWLQ overwrites the general real M-by-N matrix C with
  22 *>
  23 *>
  24 *>                    SIDE = 'L'     SIDE = 'R'
  25 *>    TRANS = 'N':      Q * C          C * Q
  26 *>    TRANS = 'T':      Q**T * C       C * Q**T
  27 *>    where Q is a real orthogonal matrix defined as the product of blocked
  28 *>    elementary reflectors computed by short wide LQ
  29 *>    factorization (SLASWLQ)
  30 *> \endverbatim
  31 *
  32 *  Arguments:
  33 *  ==========
  34 *
  35 *> \param[in] SIDE
  36 *> \verbatim
  37 *>          SIDE is CHARACTER*1
  38 *>          = 'L': apply Q or Q**T from the Left;
  39 *>          = 'R': apply Q or Q**T from the Right.
  40 *> \endverbatim
  41 *>
  42 *> \param[in] TRANS
  43 *> \verbatim
  44 *>          TRANS is CHARACTER*1
  45 *>          = 'N':  No transpose, apply Q;
  46 *>          = 'T':  Transpose, apply Q**T.
  47 *> \endverbatim
  48 *>
  49 *> \param[in] M
  50 *> \verbatim
  51 *>          M is INTEGER
  52 *>          The number of rows of the matrix C.  M >=0.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] N
  56 *> \verbatim
  57 *>          N is INTEGER
  58 *>          The number of columns of the matrix C. N >= M.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] K
  62 *> \verbatim
  63 *>          K is INTEGER
  64 *>          The number of elementary reflectors whose product defines
  65 *>          the matrix Q.
  66 *>          M >= K >= 0;
  67 *>
  68 *> \endverbatim
  69 *> \param[in] MB
  70 *> \verbatim
  71 *>          MB is INTEGER
  72 *>          The row block size to be used in the blocked QR.
  73 *>          M >= MB >= 1
  74 *> \endverbatim
  75 *>
  76 *> \param[in] NB
  77 *> \verbatim
  78 *>          NB is INTEGER
  79 *>          The column block size to be used in the blocked QR.
  80 *>          NB > M.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] NB
  84 *> \verbatim
  85 *>          NB is INTEGER
  86 *>          The block size to be used in the blocked QR.
  87 *>                MB > M.
  88 *>
  89 *> \endverbatim
  90 *>
  91 *> \param[in] A
  92 *> \verbatim
  93 *>          A is REAL array, dimension
  94 *>                               (LDA,M) if SIDE = 'L',
  95 *>                               (LDA,N) if SIDE = 'R'
  96 *>          The i-th row must contain the vector which defines the blocked
  97 *>          elementary reflector H(i), for i = 1,2,...,k, as returned by
  98 *>          SLASWLQ in the first k rows of its array argument A.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDA
 102 *> \verbatim
 103 *>          LDA is INTEGER
 104 *>          The leading dimension of the array A.
 105 *>          If SIDE = 'L', LDA >= max(1,M);
 106 *>          if SIDE = 'R', LDA >= max(1,N).
 107 *> \endverbatim
 108 *>
 109 *> \param[in] T
 110 *> \verbatim
 111 *>          T is REAL array, dimension
 112 *>          ( M * Number of blocks(CEIL(N-K/NB-K)),
 113 *>          The blocked upper triangular block reflectors stored in compact form
 114 *>          as a sequence of upper triangular blocks.  See below
 115 *>          for further details.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] LDT
 119 *> \verbatim
 120 *>          LDT is INTEGER
 121 *>          The leading dimension of the array T.  LDT >= MB.
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] C
 125 *> \verbatim
 126 *>          C is REAL array, dimension (LDC,N)
 127 *>          On entry, the M-by-N matrix C.
 128 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDC
 132 *> \verbatim
 133 *>          LDC is INTEGER
 134 *>          The leading dimension of the array C. LDC >= max(1,M).
 135 *> \endverbatim
 136 *>
 137 *> \param[out] WORK
 138 *> \verbatim
 139 *>         (workspace) REAL array, dimension (MAX(1,LWORK))
 140 *> \endverbatim
 141 *>
 142 *> \param[in] LWORK
 143 *> \verbatim
 144 *>          LWORK is INTEGER
 145 *>          The dimension of the array WORK.
 146 *>          If SIDE = 'L', LWORK >= max(1,NB) * MB;
 147 *>          if SIDE = 'R', LWORK >= max(1,M) * MB.
 148 *>          If LWORK = -1, then a workspace query is assumed; the routine
 149 *>          only calculates the optimal size of the WORK array, returns
 150 *>          this value as the first entry of the WORK array, and no error
 151 *>          message related to LWORK is issued by XERBLA.
 152 *> \endverbatim
 153 *>
 154 *> \param[out] INFO
 155 *> \verbatim
 156 *>          INFO is INTEGER
 157 *>          = 0:  successful exit
 158 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 159 *> \endverbatim
 160 *
 161 *  Authors:
 162 *  ========
 163 *
 164 *> \author Univ. of Tennessee
 165 *> \author Univ. of California Berkeley
 166 *> \author Univ. of Colorado Denver
 167 *> \author NAG Ltd.
 168 *
 169 *> \par Further Details:
 170 *  =====================
 171 *>
 172 *> \verbatim
 173 *> Short-Wide LQ (SWLQ) performs LQ by a sequence of orthogonal transformations,
 174 *> representing Q as a product of other orthogonal matrices
 175 *>   Q = Q(1) * Q(2) * . . . * Q(k)
 176 *> where each Q(i) zeros out upper diagonal entries of a block of NB rows of A:
 177 *>   Q(1) zeros out the upper diagonal entries of rows 1:NB of A
 178 *>   Q(2) zeros out the bottom MB-N rows of rows [1:M,NB+1:2*NB-M] of A
 179 *>   Q(3) zeros out the bottom MB-N rows of rows [1:M,2*NB-M+1:3*NB-2*M] of A
 180 *>   . . .
 181 *>
 182 *> Q(1) is computed by GELQT, which represents Q(1) by Householder vectors
 183 *> stored under the diagonal of rows 1:MB of A, and by upper triangular
 184 *> block reflectors, stored in array T(1:LDT,1:N).
 185 *> For more information see Further Details in GELQT.
 186 *>
 187 *> Q(i) for i>1 is computed by TPLQT, which represents Q(i) by Householder vectors
 188 *> stored in columns [(i-1)*(NB-M)+M+1:i*(NB-M)+M] of A, and by upper triangular
 189 *> block reflectors, stored in array T(1:LDT,(i-1)*M+1:i*M).
 190 *> The last Q(k) may use fewer rows.
 191 *> For more information see Further Details in TPQRT.
 192 *>
 193 *> For more details of the overall algorithm, see the description of
 194 *> Sequential TSQR in Section 2.2 of [1].
 195 *>
 196 *> [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,”
 197 *>     J. Demmel, L. Grigori, M. Hoemmen, J. Langou,
 198 *>     SIAM J. Sci. Comput, vol. 34, no. 1, 2012
 199 *> \endverbatim
 200 *>
 201 *  =====================================================================
 202       SUBROUTINE SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T,
 203      $    LDT, C, LDC, WORK, LWORK, INFO )
 204 *
 205 *  -- LAPACK computational routine (version 3.7.0) --
 206 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 207 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 208 *     December 2016
 209 *
 210 *     .. Scalar Arguments ..
 211       CHARACTER         SIDE, TRANS
 212       INTEGER           INFO, LDA, M, N, K, MB, NB, LDT, LWORK, LDC
 213 *     ..
 214 *     .. Array Arguments ..
 215       REAL              A( LDA, * ), WORK( * ), C(LDC, * ),
 216      $      T( LDT, * )
 217 *     ..
 218 *
 219 * =====================================================================
 220 *
 221 *     ..
 222 *     .. Local Scalars ..
 223       LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
 224       INTEGER    I, II, KK, LW, CTR
 225 *     ..
 226 *     .. External Functions ..
 227       LOGICAL            LSAME
 228       EXTERNAL           LSAME
 229 *     .. External Subroutines ..
 230       EXTERNAL           STPMLQT, SGEMLQT, XERBLA
 231 *     ..
 232 *     .. Executable Statements ..
 233 *
 234 *     Test the input arguments
 235 *
 236       LQUERY  = LWORK.LT.0
 237       NOTRAN  = LSAME( TRANS, 'N' )
 238       TRAN    = LSAME( TRANS, 'T' )
 239       LEFT    = LSAME( SIDE, 'L' )
 240       RIGHT   = LSAME( SIDE, 'R' )
 241       IF (LEFT) THEN
 242         LW = N * MB
 243       ELSE
 244         LW = M * MB
 245       END IF
 246 *
 247       INFO = 0
 248       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
 249          INFO = -1
 250       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
 251          INFO = -2
 252       ELSE IF( M.LT.0 ) THEN
 253         INFO = -3
 254       ELSE IF( N.LT.0 ) THEN
 255         INFO = -4
 256       ELSE IF( K.LT.0 ) THEN
 257         INFO = -5
 258       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
 259         INFO = -9
 260       ELSE IF( LDT.LT.MAX( 1, MB) ) THEN
 261         INFO = -11
 262       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
 263          INFO = -13
 264       ELSE IF(( LWORK.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
 265         INFO = -15
 266       END IF
 267 *
 268       IF( INFO.NE.0 ) THEN
 269         CALL XERBLA( 'SLAMSWLQ', -INFO )
 270         WORK(1) = LW
 271         RETURN
 272       ELSE IF (LQUERY) THEN
 273         WORK(1) = LW
 274         RETURN
 275       END IF
 276 *
 277 *     Quick return if possible
 278 *
 279       IF( MIN(M,N,K).EQ.0 ) THEN
 280         RETURN
 281       END IF
 282 *
 283       IF((NB.LE.K).OR.(NB.GE.MAX(M,N,K))) THEN
 284         CALL DGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
 285      $        T, LDT, C, LDC, WORK, INFO)
 286         RETURN
 287       END IF
 288 *
 289       IF(LEFT.AND.TRAN) THEN
 290 *
 291 *         Multiply Q to the last block of C
 292 *
 293           KK = MOD((M-K),(NB-K))
 294           CTR = (M-K)/(NB-K)
 295 *
 296           IF (KK.GT.0) THEN
 297             II=M-KK+1
 298             CALL STPMLQT('L','T',KK , N, K, 0, MB, A(1,II), LDA,
 299      $        T(1,CTR*K+1), LDT, C(1,1), LDC,
 300      $        C(II,1), LDC, WORK, INFO )
 301           ELSE
 302             II=M+1
 303           END IF
 304 *
 305           DO I=II-(NB-K),NB+1,-(NB-K)
 306 *
 307 *         Multiply Q to the current block of C (1:M,I:I+NB)
 308 *
 309             CTR = CTR - 1
 310             CALL STPMLQT('L','T',NB-K , N, K, 0,MB, A(1,I), LDA,
 311      $          T(1,CTR*K+1),LDT, C(1,1), LDC,
 312      $          C(I,1), LDC, WORK, INFO )
 313           END DO
 314 *
 315 *         Multiply Q to the first block of C (1:M,1:NB)
 316 *
 317           CALL SGEMLQT('L','T',NB , N, K, MB, A(1,1), LDA, T
 318      $              ,LDT ,C(1,1), LDC, WORK, INFO )
 319 *
 320       ELSE IF (LEFT.AND.NOTRAN) THEN
 321 *
 322 *         Multiply Q to the first block of C
 323 *
 324          KK = MOD((M-K),(NB-K))
 325          II=M-KK+1
 326          CTR = 1
 327          CALL SGEMLQT('L','N',NB , N, K, MB, A(1,1), LDA, T
 328      $              ,LDT ,C(1,1), LDC, WORK, INFO )
 329 *
 330          DO I=NB+1,II-NB+K,(NB-K)
 331 *
 332 *         Multiply Q to the current block of C (I:I+NB,1:N)
 333 *
 334           CALL STPMLQT('L','N',NB-K , N, K, 0,MB, A(1,I), LDA,
 335      $         T(1,CTR * K+1), LDT, C(1,1), LDC,
 336      $         C(I,1), LDC, WORK, INFO )
 337           CTR = CTR + 1
 338 *
 339          END DO
 340          IF(II.LE.M) THEN
 341 *
 342 *         Multiply Q to the last block of C
 343 *
 344           CALL STPMLQT('L','N',KK , N, K, 0, MB, A(1,II), LDA,
 345      $        T(1,CTR*K+1), LDT, C(1,1), LDC,
 346      $        C(II,1), LDC, WORK, INFO )
 347 *
 348          END IF
 349 *
 350       ELSE IF(RIGHT.AND.NOTRAN) THEN
 351 *
 352 *         Multiply Q to the last block of C
 353 *
 354           KK = MOD((N-K),(NB-K))
 355           CTR = (N-K)/(NB-K)
 356           IF (KK.GT.0) THEN
 357             II=N-KK+1
 358             CALL STPMLQT('R','N',M , KK, K, 0, MB, A(1, II), LDA,
 359      $        T(1,CTR*K+1), LDT, C(1,1), LDC,
 360      $        C(1,II), LDC, WORK, INFO )
 361           ELSE
 362             II=N+1
 363           END IF
 364 *
 365           DO I=II-(NB-K),NB+1,-(NB-K)
 366 *
 367 *         Multiply Q to the current block of C (1:M,I:I+MB)
 368 *
 369              CTR = CTR - 1
 370              CALL STPMLQT('R','N', M, NB-K, K, 0, MB, A(1, I), LDA,
 371      $            T(1,CTR*K+1), LDT, C(1,1), LDC,
 372      $            C(1,I), LDC, WORK, INFO )
 373
 374           END DO
 375 *
 376 *         Multiply Q to the first block of C (1:M,1:MB)
 377 *
 378           CALL SGEMLQT('R','N',M , NB, K, MB, A(1,1), LDA, T
 379      $            ,LDT ,C(1,1), LDC, WORK, INFO )
 380 *
 381       ELSE IF (RIGHT.AND.TRAN) THEN
 382 *
 383 *       Multiply Q to the first block of C
 384 *
 385          KK = MOD((N-K),(NB-K))
 386          II=N-KK+1
 387          CTR = 1
 388          CALL SGEMLQT('R','T',M , NB, K, MB, A(1,1), LDA, T
 389      $            ,LDT ,C(1,1), LDC, WORK, INFO )
 390 *
 391          DO I=NB+1,II-NB+K,(NB-K)
 392 *
 393 *         Multiply Q to the current block of C (1:M,I:I+MB)
 394 *
 395           CALL STPMLQT('R','T',M , NB-K, K, 0,MB, A(1,I), LDA,
 396      $       T(1, CTR*K+1), LDT, C(1,1), LDC,
 397      $       C(1,I), LDC, WORK, INFO )
 398           CTR = CTR + 1
 399 *
 400          END DO
 401          IF(II.LE.N) THEN
 402 *
 403 *       Multiply Q to the last block of C
 404 *
 405           CALL STPMLQT('R','T',M , KK, K, 0,MB, A(1,II), LDA,
 406      $      T(1,CTR*K+1),LDT, C(1,1), LDC,
 407      $      C(1,II), LDC, WORK, INFO )
 408 *
 409          END IF
 410 *
 411       END IF
 412 *
 413       WORK(1) = LW
 414       RETURN
 415 *
 416 *     End of SLAMSWLQ
 417 *
 418       END