![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTPQRT2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTPQRT2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztpqrt2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztpqrt2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztpqrt2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, LDB, LDT, N, M, L 00025 * .. 00026 * .. Array Arguments .. 00027 * COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> ZTPQRT2 computes a QR factorization of a complex "triangular-pentagonal" 00037 *> matrix C, which is composed of a triangular block A and pentagonal block B, 00038 *> using the compact WY representation for Q. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] M 00045 *> \verbatim 00046 *> M is INTEGER 00047 *> The total number of rows of the matrix B. 00048 *> M >= 0. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] N 00052 *> \verbatim 00053 *> N is INTEGER 00054 *> The number of columns of the matrix B, and the order of 00055 *> the triangular matrix A. 00056 *> N >= 0. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] L 00060 *> \verbatim 00061 *> L is INTEGER 00062 *> The number of rows of the upper trapezoidal part of B. 00063 *> MIN(M,N) >= L >= 0. See Further Details. 00064 *> \endverbatim 00065 *> 00066 *> \param[in,out] A 00067 *> \verbatim 00068 *> A is COMPLEX*16 array, dimension (LDA,N) 00069 *> On entry, the upper triangular N-by-N matrix A. 00070 *> On exit, the elements on and above the diagonal of the array 00071 *> contain the upper triangular matrix R. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] LDA 00075 *> \verbatim 00076 *> LDA is INTEGER 00077 *> The leading dimension of the array A. LDA >= max(1,N). 00078 *> \endverbatim 00079 *> 00080 *> \param[in,out] B 00081 *> \verbatim 00082 *> B is COMPLEX*16 array, dimension (LDB,N) 00083 *> On entry, the pentagonal M-by-N matrix B. The first M-L rows 00084 *> are rectangular, and the last L rows are upper trapezoidal. 00085 *> On exit, B contains the pentagonal matrix V. See Further Details. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] LDB 00089 *> \verbatim 00090 *> LDB is INTEGER 00091 *> The leading dimension of the array B. LDB >= max(1,M). 00092 *> \endverbatim 00093 *> 00094 *> \param[out] T 00095 *> \verbatim 00096 *> T is COMPLEX*16 array, dimension (LDT,N) 00097 *> The N-by-N upper triangular factor T of the block reflector. 00098 *> See Further Details. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDT 00102 *> \verbatim 00103 *> LDT is INTEGER 00104 *> The leading dimension of the array T. LDT >= max(1,N) 00105 *> \endverbatim 00106 *> 00107 *> \param[out] INFO 00108 *> \verbatim 00109 *> INFO is INTEGER 00110 *> = 0: successful exit 00111 *> < 0: if INFO = -i, the i-th argument had an illegal value 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date April 2012 00123 * 00124 *> \ingroup complex16OTHERcomputational 00125 * 00126 *> \par Further Details: 00127 * ===================== 00128 *> 00129 *> \verbatim 00130 *> 00131 *> The input matrix C is a (N+M)-by-N matrix 00132 *> 00133 *> C = [ A ] 00134 *> [ B ] 00135 *> 00136 *> where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal 00137 *> matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N 00138 *> upper trapezoidal matrix B2: 00139 *> 00140 *> B = [ B1 ] <- (M-L)-by-N rectangular 00141 *> [ B2 ] <- L-by-N upper trapezoidal. 00142 *> 00143 *> The upper trapezoidal matrix B2 consists of the first L rows of a 00144 *> N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N). If L=0, 00145 *> B is rectangular M-by-N; if M=L=N, B is upper triangular. 00146 *> 00147 *> The matrix W stores the elementary reflectors H(i) in the i-th column 00148 *> below the diagonal (of A) in the (N+M)-by-N input matrix C 00149 *> 00150 *> C = [ A ] <- upper triangular N-by-N 00151 *> [ B ] <- M-by-N pentagonal 00152 *> 00153 *> so that W can be represented as 00154 *> 00155 *> W = [ I ] <- identity, N-by-N 00156 *> [ V ] <- M-by-N, same form as B. 00157 *> 00158 *> Thus, all of information needed for W is contained on exit in B, which 00159 *> we call V above. Note that V has the same form as B; that is, 00160 *> 00161 *> V = [ V1 ] <- (M-L)-by-N rectangular 00162 *> [ V2 ] <- L-by-N upper trapezoidal. 00163 *> 00164 *> The columns of V represent the vectors which define the H(i)'s. 00165 *> The (M+N)-by-(M+N) block reflector H is then given by 00166 *> 00167 *> H = I - W * T * W**H 00168 *> 00169 *> where W**H is the conjugate transpose of W and T is the upper triangular 00170 *> factor of the block reflector. 00171 *> \endverbatim 00172 *> 00173 * ===================================================================== 00174 SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) 00175 * 00176 * -- LAPACK computational routine (version 3.4.1) -- 00177 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00178 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00179 * April 2012 00180 * 00181 * .. Scalar Arguments .. 00182 INTEGER INFO, LDA, LDB, LDT, N, M, L 00183 * .. 00184 * .. Array Arguments .. 00185 COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ) 00186 * .. 00187 * 00188 * ===================================================================== 00189 * 00190 * .. Parameters .. 00191 COMPLEX*16 ONE, ZERO 00192 PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) ) 00193 * .. 00194 * .. Local Scalars .. 00195 INTEGER I, J, P, MP, NP 00196 COMPLEX*16 ALPHA 00197 * .. 00198 * .. External Subroutines .. 00199 EXTERNAL ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA 00200 * .. 00201 * .. Intrinsic Functions .. 00202 INTRINSIC MAX, MIN 00203 * .. 00204 * .. Executable Statements .. 00205 * 00206 * Test the input arguments 00207 * 00208 INFO = 0 00209 IF( M.LT.0 ) THEN 00210 INFO = -1 00211 ELSE IF( N.LT.0 ) THEN 00212 INFO = -2 00213 ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN 00214 INFO = -3 00215 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00216 INFO = -5 00217 ELSE IF( LDB.LT.MAX( 1, M ) ) THEN 00218 INFO = -7 00219 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00220 INFO = -9 00221 END IF 00222 IF( INFO.NE.0 ) THEN 00223 CALL XERBLA( 'ZTPQRT2', -INFO ) 00224 RETURN 00225 END IF 00226 * 00227 * Quick return if possible 00228 * 00229 IF( N.EQ.0 .OR. M.EQ.0 ) RETURN 00230 * 00231 DO I = 1, N 00232 * 00233 * Generate elementary reflector H(I) to annihilate B(:,I) 00234 * 00235 P = M-L+MIN( L, I ) 00236 CALL ZLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) ) 00237 IF( I.LT.N ) THEN 00238 * 00239 * W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)] 00240 * 00241 DO J = 1, N-I 00242 T( J, N ) = CONJG(A( I, I+J )) 00243 END DO 00244 CALL ZGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB, 00245 $ B( 1, I ), 1, ONE, T( 1, N ), 1 ) 00246 * 00247 * C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H 00248 * 00249 ALPHA = -CONJG(T( I, 1 )) 00250 DO J = 1, N-I 00251 A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N )) 00252 END DO 00253 CALL ZGERC( P, N-I, ALPHA, B( 1, I ), 1, 00254 $ T( 1, N ), 1, B( 1, I+1 ), LDB ) 00255 END IF 00256 END DO 00257 * 00258 DO I = 2, N 00259 * 00260 * T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I)) 00261 * 00262 ALPHA = -T( I, 1 ) 00263 00264 DO J = 1, I-1 00265 T( J, I ) = ZERO 00266 END DO 00267 P = MIN( I-1, L ) 00268 MP = MIN( M-L+1, M ) 00269 NP = MIN( P+1, N ) 00270 * 00271 * Triangular part of B2 00272 * 00273 DO J = 1, P 00274 T( J, I ) = ALPHA*B( M-L+J, I ) 00275 END DO 00276 CALL ZTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB, 00277 $ T( 1, I ), 1 ) 00278 * 00279 * Rectangular part of B2 00280 * 00281 CALL ZGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB, 00282 $ B( MP, I ), 1, ZERO, T( NP, I ), 1 ) 00283 * 00284 * B1 00285 * 00286 CALL ZGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 00287 $ ONE, T( 1, I ), 1 ) 00288 * 00289 * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I) 00290 * 00291 CALL ZTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 ) 00292 * 00293 * T(I,I) = tau(I) 00294 * 00295 T( I, I ) = T( I, 1 ) 00296 T( I, 1 ) = ZERO 00297 END DO 00298 00299 * 00300 * End of ZTPQRT2 00301 * 00302 END