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