SRC/ztplqt2.f

   1 *> \brief \b ZTPLQT2 computes a LQ 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 ZTPLQT2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztplqt2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztplqt2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztplqt2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZTPLQT2( 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*16   A( LDA, * ), B( LDB, * ), T( LDT, * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> ZTPLQT2 computes a LQ a 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 lower 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*16 array, dimension (LDA,N)
  69 *>          On entry, the lower triangular M-by-M matrix A.
  70 *>          On exit, the elements on and below the diagonal of the array
  71 *>          contain the lower triangular matrix L.
  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*16 array, dimension (LDB,N)
  83 *>          On entry, the pentagonal M-by-N matrix B.  The first N-L columns
  84 *>          are rectangular, and the last L columns are lower 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*16 array, dimension (LDT,M)
  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,M)
 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 doubleOTHERcomputational
 125 *
 126 *> \par Further Details:
 127 *  =====================
 128 *>
 129 *> \verbatim
 130 *>
 131 *>  The input matrix C is a M-by-(M+N) matrix
 132 *>
 133 *>               C = [ A ][ B ]
 134 *>
 135 *>
 136 *>  where A is an lower triangular N-by-N matrix, and B is M-by-N pentagonal
 137 *>  matrix consisting of a M-by-(N-L) rectangular matrix B1 left of a M-by-L
 138 *>  upper trapezoidal matrix B2:
 139 *>
 140 *>               B = [ B1 ][ B2 ]
 141 *>                   [ B1 ]  <-     M-by-(N-L) rectangular
 142 *>                   [ B2 ]  <-     M-by-L lower trapezoidal.
 143 *>
 144 *>  The lower trapezoidal matrix B2 consists of the first L columns of a
 145 *>  N-by-N lower triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,
 146 *>  B is rectangular M-by-N; if M=L=N, B is lower triangular.
 147 *>
 148 *>  The matrix W stores the elementary reflectors H(i) in the i-th row
 149 *>  above the diagonal (of A) in the M-by-(M+N) input matrix C
 150 *>
 151 *>               C = [ A ][ B ]
 152 *>                   [ A ]  <- lower triangular N-by-N
 153 *>                   [ B ]  <- M-by-N pentagonal
 154 *>
 155 *>  so that W can be represented as
 156 *>
 157 *>               W = [ I ][ V ]
 158 *>                   [ I ]  <- identity, N-by-N
 159 *>                   [ V ]  <- M-by-N, same form as B.
 160 *>
 161 *>  Thus, all of information needed for W is contained on exit in B, which
 162 *>  we call V above.  Note that V has the same form as B; that is,
 163 *>
 164 *>               W = [ V1 ][ V2 ]
 165 *>                   [ V1 ] <-     M-by-(N-L) rectangular
 166 *>                   [ V2 ] <-     M-by-L lower trapezoidal.
 167 *>
 168 *>  The rows of V represent the vectors which define the H(i)'s.
 169 *>  The (M+N)-by-(M+N) block reflector H is then given by
 170 *>
 171 *>               H = I - W**T * T * W
 172 *>
 173 *>  where W^H is the conjugate transpose of W and T is the upper triangular
 174 *>  factor of the block reflector.
 175 *> \endverbatim
 176 *>
 177 *  =====================================================================
 178       SUBROUTINE ZTPLQT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
 179 *
 180 *  -- LAPACK computational routine (version 3.4.2) --
 181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 183 *     September 2012
 184 *
 185 *     .. Scalar Arguments ..
 186       INTEGER        INFO, LDA, LDB, LDT, N, M, L
 187 *     ..
 188 *     .. Array Arguments ..
 189       COMPLEX*16     A( LDA, * ), B( LDB, * ), T( LDT, * )
 190 *     ..
 191 *
 192 *  =====================================================================
 193 *
 194 *     .. Parameters ..
 195       COMPLEX*16  ONE, ZERO
 196       PARAMETER( ZERO = ( 0.0D+0, 0.0D+0 ),ONE  = ( 1.0D+0, 0.0D+0 ) )
 197 *     ..
 198 *     .. Local Scalars ..
 199       INTEGER   I, J, P, MP, NP
 200       COMPLEX*16   ALPHA
 201 *     ..
 202 *     .. External Subroutines ..
 203       EXTERNAL  ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA
 204 *     ..
 205 *     .. Intrinsic Functions ..
 206       INTRINSIC MAX, MIN
 207 *     ..
 208 *     .. Executable Statements ..
 209 *
 210 *     Test the input arguments
 211 *
 212       INFO = 0
 213       IF( M.LT.0 ) THEN
 214          INFO = -1
 215       ELSE IF( N.LT.0 ) THEN
 216          INFO = -2
 217       ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
 218          INFO = -3
 219       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 220          INFO = -5
 221       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
 222          INFO = -7
 223       ELSE IF( LDT.LT.MAX( 1, M ) ) THEN
 224          INFO = -9
 225       END IF
 226       IF( INFO.NE.0 ) THEN
 227          CALL XERBLA( 'ZTPLQT2', -INFO )
 228          RETURN
 229       END IF
 230 *
 231 *     Quick return if possible
 232 *
 233       IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
 234 *
 235       DO I = 1, M
 236 *
 237 *        Generate elementary reflector H(I) to annihilate B(I,:)
 238 *
 239          P = N-L+MIN( L, I )
 240          CALL ZLARFG( P+1, A( I, I ), B( I, 1 ), LDB, T( 1, I ) )
 241          T(1,I)=CONJG(T(1,I))
 242          IF( I.LT.M ) THEN
 243             DO J = 1, P
 244                B( I, J ) = CONJG(B(I,J))
 245             END DO
 246 *
 247 *           W(M-I:1) := C(I+1:M,I:N) * C(I,I:N) [use W = T(M,:)]
 248 *
 249             DO J = 1, M-I
 250                T( M, J ) = (A( I+J, I ))
 251             END DO
 252             CALL ZGEMV( 'N', M-I, P, ONE, B( I+1, 1 ), LDB,
 253      $                  B( I, 1 ), LDB, ONE, T( M, 1 ), LDT )
 254 *
 255 *           C(I+1:M,I:N) = C(I+1:M,I:N) + alpha * C(I,I:N)*W(M-1:1)^H
 256 *
 257             ALPHA = -(T( 1, I ))
 258             DO J = 1, M-I
 259                A( I+J, I ) = A( I+J, I ) + ALPHA*(T( M, J ))
 260             END DO
 261             CALL ZGERC( M-I, P, (ALPHA),  T( M, 1 ), LDT,
 262      $          B( I, 1 ), LDB, B( I+1, 1 ), LDB )
 263             DO J = 1, P
 264                B( I, J ) = CONJG(B(I,J))
 265             END DO
 266          END IF
 267       END DO
 268 *
 269       DO I = 2, M
 270 *
 271 *        T(I,1:I-1) := C(I:I-1,1:N)**H * (alpha * C(I,I:N))
 272 *
 273          ALPHA = -(T( 1, I ))
 274          DO J = 1, I-1
 275             T( I, J ) = ZERO
 276          END DO
 277          P = MIN( I-1, L )
 278          NP = MIN( N-L+1, N )
 279          MP = MIN( P+1, M )
 280          DO J = 1, N-L+P
 281            B(I,J)=CONJG(B(I,J))
 282          END DO
 283 *
 284 *        Triangular part of B2
 285 *
 286          DO J = 1, P
 287             T( I, J ) = (ALPHA*B( I, N-L+J ))
 288          END DO
 289          CALL ZTRMV( 'L', 'N', 'N', P, B( 1, NP ), LDB,
 290      $               T( I, 1 ), LDT )
 291 *
 292 *        Rectangular part of B2
 293 *
 294          CALL ZGEMV( 'N', I-1-P, L,  ALPHA, B( MP, NP ), LDB,
 295      $               B( I, NP ), LDB, ZERO, T( I,MP ), LDT )
 296 *
 297 *        B1
 298
 299 *
 300          CALL ZGEMV( 'N', I-1, N-L, ALPHA, B, LDB, B( I, 1 ), LDB,
 301      $               ONE, T( I, 1 ), LDT )
 302 *
 303
 304 *
 305 *        T(1:I-1,I) := T(1:I-1,1:I-1) * T(I,1:I-1)
 306 *
 307          DO J = 1, I-1
 308             T(I,J)=CONJG(T(I,J))
 309          END DO
 310          CALL ZTRMV( 'L', 'C', 'N', I-1, T, LDT, T( I, 1 ), LDT )
 311          DO J = 1, I-1
 312             T(I,J)=CONJG(T(I,J))
 313          END DO
 314          DO J = 1, N-L+P
 315             B(I,J)=CONJG(B(I,J))
 316          END DO
 317 *
 318 *        T(I,I) = tau(I)
 319 *
 320          T( I, I ) = T( 1, I )
 321          T( 1, I ) = ZERO
 322       END DO
 323       DO I=1,M
 324          DO J= I+1,M
 325             T(I,J)=(T(J,I))
 326             T(J,I)=ZERO
 327          END DO
 328       END DO
 329
 330 *
 331 *     End of ZTPLQT2
 332 *
 333       END