SRC/ctpqrt2.f

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