SRC/ctplqt.f

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