SRC/dtpqrt.f

   1 *> \brief \b DTPQRT
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTPQRT + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtpqrt.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtpqrt.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtpqrt.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DTPQRT computes a blocked QR factorization of a real
  38 *> "triangular-pentagonal" matrix C, which is composed of a
  39 *> triangular block A and pentagonal block B, using the compact
  40 *> WY representation for Q.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] M
  47 *> \verbatim
  48 *>          M is INTEGER
  49 *>          The number of rows of the matrix B.
  50 *>          M >= 0.
  51 *> \endverbatim
  52 *>
  53 *> \param[in] N
  54 *> \verbatim
  55 *>          N is INTEGER
  56 *>          The number of columns of the matrix B, and the order of the
  57 *>          triangular matrix A.
  58 *>          N >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] L
  62 *> \verbatim
  63 *>          L is INTEGER
  64 *>          The number of rows of the upper trapezoidal part of B.
  65 *>          MIN(M,N) >= L >= 0.  See Further Details.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] NB
  69 *> \verbatim
  70 *>          NB is INTEGER
  71 *>          The block size to be used in the blocked QR.  N >= NB >= 1.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] A
  75 *> \verbatim
  76 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  77 *>          On entry, the upper triangular N-by-N matrix A.
  78 *>          On exit, the elements on and above the diagonal of the array
  79 *>          contain the upper triangular matrix R.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LDA
  83 *> \verbatim
  84 *>          LDA is INTEGER
  85 *>          The leading dimension of the array A.  LDA >= max(1,N).
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] B
  89 *> \verbatim
  90 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
  91 *>          On entry, the pentagonal M-by-N matrix B.  The first M-L rows
  92 *>          are rectangular, and the last L rows are upper trapezoidal.
  93 *>          On exit, B contains the pentagonal matrix V.  See Further Details.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] LDB
  97 *> \verbatim
  98 *>          LDB is INTEGER
  99 *>          The leading dimension of the array B.  LDB >= max(1,M).
 100 *> \endverbatim
 101 *>
 102 *> \param[out] T
 103 *> \verbatim
 104 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
 105 *>          The upper triangular block reflectors stored in compact form
 106 *>          as a sequence of upper triangular blocks.  See Further Details.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LDT
 110 *> \verbatim
 111 *>          LDT is INTEGER
 112 *>          The leading dimension of the array T.  LDT >= NB.
 113 *> \endverbatim
 114 *>
 115 *> \param[out] WORK
 116 *> \verbatim
 117 *>          WORK is DOUBLE PRECISION array, dimension (NB*N)
 118 *> \endverbatim
 119 *>
 120 *> \param[out] INFO
 121 *> \verbatim
 122 *>          INFO is INTEGER
 123 *>          = 0:  successful exit
 124 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 125 *> \endverbatim
 126 *
 127 *  Authors:
 128 *  ========
 129 *
 130 *> \author Univ. of Tennessee
 131 *> \author Univ. of California Berkeley
 132 *> \author Univ. of Colorado Denver
 133 *> \author NAG Ltd.
 134 *
 135 *> \date November 2013
 136 *
 137 *> \ingroup doubleOTHERcomputational
 138 *
 139 *> \par Further Details:
 140 *  =====================
 141 *>
 142 *> \verbatim
 143 *>
 144 *>  The input matrix C is a (N+M)-by-N matrix
 145 *>
 146 *>               C = [ A ]
 147 *>                   [ B ]
 148 *>
 149 *>  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
 150 *>  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
 151 *>  upper trapezoidal matrix B2:
 152 *>
 153 *>               B = [ B1 ]  <- (M-L)-by-N rectangular
 154 *>                   [ B2 ]  <-     L-by-N upper trapezoidal.
 155 *>
 156 *>  The upper trapezoidal matrix B2 consists of the first L rows of a
 157 *>  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,
 158 *>  B is rectangular M-by-N; if M=L=N, B is upper triangular.
 159 *>
 160 *>  The matrix W stores the elementary reflectors H(i) in the i-th column
 161 *>  below the diagonal (of A) in the (N+M)-by-N input matrix C
 162 *>
 163 *>               C = [ A ]  <- upper triangular N-by-N
 164 *>                   [ B ]  <- M-by-N pentagonal
 165 *>
 166 *>  so that W can be represented as
 167 *>
 168 *>               W = [ I ]  <- identity, N-by-N
 169 *>                   [ V ]  <- M-by-N, same form as B.
 170 *>
 171 *>  Thus, all of information needed for W is contained on exit in B, which
 172 *>  we call V above.  Note that V has the same form as B; that is,
 173 *>
 174 *>               V = [ V1 ] <- (M-L)-by-N rectangular
 175 *>                   [ V2 ] <-     L-by-N upper trapezoidal.
 176 *>
 177 *>  The columns of V represent the vectors which define the H(i)'s.
 178 *>
 179 *>  The number of blocks is B = ceiling(N/NB), where each
 180 *>  block is of order NB except for the last block, which is of order
 181 *>  IB = N - (B-1)*NB.  For each of the B blocks, a upper triangular block
 182 *>  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB
 183 *>  for the last block) T's are stored in the NB-by-N matrix T as
 184 *>
 185 *>               T = [T1 T2 ... TB].
 186 *> \endverbatim
 187 *>
 188 *  =====================================================================
 189       SUBROUTINE DTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK,
 190      $                   INFO )
 191 *
 192 *  -- LAPACK computational routine (version 3.5.0) --
 193 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 194 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 195 *     November 2013
 196 *
 197 *     .. Scalar Arguments ..
 198       INTEGER INFO, LDA, LDB, LDT, N, M, L, NB
 199 *     ..
 200 *     .. Array Arguments ..
 201       DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * )
 202 *     ..
 203 *
 204 * =====================================================================
 205 *
 206 *     ..
 207 *     .. Local Scalars ..
 208       INTEGER    I, IB, LB, MB, IINFO
 209 *     ..
 210 *     .. External Subroutines ..
 211       EXTERNAL   DTPQRT2, DTPRFB, XERBLA
 212 *     ..
 213 *     .. Executable Statements ..
 214 *
 215 *     Test the input arguments
 216 *
 217       INFO = 0
 218       IF( M.LT.0 ) THEN
 219          INFO = -1
 220       ELSE IF( N.LT.0 ) THEN
 221          INFO = -2
 222       ELSE IF( L.LT.0 .OR. (L.GT.MIN(M,N) .AND. MIN(M,N).GE.0)) THEN
 223          INFO = -3
 224       ELSE IF( NB.LT.1 .OR. (NB.GT.N .AND. N.GT.0)) THEN
 225          INFO = -4
 226       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 227          INFO = -6
 228       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
 229          INFO = -8
 230       ELSE IF( LDT.LT.NB ) THEN
 231          INFO = -10
 232       END IF
 233       IF( INFO.NE.0 ) THEN
 234          CALL XERBLA( 'DTPQRT', -INFO )
 235          RETURN
 236       END IF
 237 *
 238 *     Quick return if possible
 239 *
 240       IF( M.EQ.0 .OR. N.EQ.0 ) RETURN
 241 *
 242       DO I = 1, N, NB
 243 *
 244 *     Compute the QR factorization of the current block
 245 *
 246          IB = MIN( N-I+1, NB )
 247          MB = MIN( M-L+I+IB-1, M )
 248          IF( I.GE.L ) THEN
 249             LB = 0
 250          ELSE
 251             LB = MB-M+L-I+1
 252          END IF
 253 *
 254          CALL DTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB,
 255      $                 T(1, I ), LDT, IINFO )
 256 *
 257 *     Update by applying H**T to B(:,I+IB:N) from the left
 258 *
 259          IF( I+IB.LE.N ) THEN
 260             CALL DTPRFB( 'L', 'T', 'F', 'C', MB, N-I-IB+1, IB, LB,
 261      $                    B( 1, I ), LDB, T( 1, I ), LDT,
 262      $                    A( I, I+IB ), LDA, B( 1, I+IB ), LDB,
 263      $                    WORK, IB )
 264          END IF
 265       END DO
 266       RETURN
 267 *
 268 *     End of DTPQRT
 269 *
 270       END