SRC/dlatsqr.f

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