![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGEQRT3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGEQRT3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeqrt3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeqrt3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeqrt3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * RECURSIVE SUBROUTINE CGEQRT3( M, N, A, LDA, T, LDT, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, M, N, LDT 00025 * .. 00026 * .. Array Arguments .. 00027 * COMPLEX A( LDA, * ), T( LDT, * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> CGEQRT3 recursively computes a QR factorization of a complex M-by-N matrix A, 00037 *> using the compact WY representation of Q. 00038 *> 00039 *> Based on the algorithm of Elmroth and Gustavson, 00040 *> IBM J. Res. Develop. Vol 44 No. 4 July 2000. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] M 00047 *> \verbatim 00048 *> M is INTEGER 00049 *> The number of rows of the matrix A. M >= N. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] N 00053 *> \verbatim 00054 *> N is INTEGER 00055 *> The number of columns of the matrix A. N >= 0. 00056 *> \endverbatim 00057 *> 00058 *> \param[in,out] A 00059 *> \verbatim 00060 *> A is COMPLEX array, dimension (LDA,N) 00061 *> On entry, the complex M-by-N matrix A. On exit, the elements on and 00062 *> above the diagonal contain the N-by-N upper triangular matrix R; the 00063 *> elements below the diagonal are the columns of V. See below for 00064 *> further details. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] LDA 00068 *> \verbatim 00069 *> LDA is INTEGER 00070 *> The leading dimension of the array A. LDA >= max(1,M). 00071 *> \endverbatim 00072 *> 00073 *> \param[out] T 00074 *> \verbatim 00075 *> T is COMPLEX array, dimension (LDT,N) 00076 *> The N-by-N upper triangular factor of the block reflector. 00077 *> The elements on and above the diagonal contain the block 00078 *> reflector T; the elements below the diagonal are not used. 00079 *> See below for further details. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDT 00083 *> \verbatim 00084 *> LDT is INTEGER 00085 *> The leading dimension of the array T. LDT >= max(1,N). 00086 *> \endverbatim 00087 *> 00088 *> \param[out] INFO 00089 *> \verbatim 00090 *> INFO is INTEGER 00091 *> = 0: successful exit 00092 *> < 0: if INFO = -i, the i-th argument had an illegal value 00093 *> \endverbatim 00094 * 00095 * Authors: 00096 * ======== 00097 * 00098 *> \author Univ. of Tennessee 00099 *> \author Univ. of California Berkeley 00100 *> \author Univ. of Colorado Denver 00101 *> \author NAG Ltd. 00102 * 00103 *> \date November 2011 00104 * 00105 *> \ingroup complexGEcomputational 00106 * 00107 *> \par Further Details: 00108 * ===================== 00109 *> 00110 *> \verbatim 00111 *> 00112 *> The matrix V stores the elementary reflectors H(i) in the i-th column 00113 *> below the diagonal. For example, if M=5 and N=3, the matrix V is 00114 *> 00115 *> V = ( 1 ) 00116 *> ( v1 1 ) 00117 *> ( v1 v2 1 ) 00118 *> ( v1 v2 v3 ) 00119 *> ( v1 v2 v3 ) 00120 *> 00121 *> where the vi's represent the vectors which define H(i), which are returned 00122 *> in the matrix A. The 1's along the diagonal of V are not stored in A. The 00123 *> block reflector H is then given by 00124 *> 00125 *> H = I - V * T * V**H 00126 *> 00127 *> where V**H is the conjugate transpose of V. 00128 *> 00129 *> For details of the algorithm, see Elmroth and Gustavson (cited above). 00130 *> \endverbatim 00131 *> 00132 * ===================================================================== 00133 RECURSIVE SUBROUTINE CGEQRT3( M, N, A, LDA, T, LDT, INFO ) 00134 * 00135 * -- LAPACK computational routine (version 3.4.0) -- 00136 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00138 * November 2011 00139 * 00140 * .. Scalar Arguments .. 00141 INTEGER INFO, LDA, M, N, LDT 00142 * .. 00143 * .. Array Arguments .. 00144 COMPLEX A( LDA, * ), T( LDT, * ) 00145 * .. 00146 * 00147 * ===================================================================== 00148 * 00149 * .. Parameters .. 00150 COMPLEX ONE 00151 PARAMETER ( ONE = (1.0,0.0) ) 00152 * .. 00153 * .. Local Scalars .. 00154 INTEGER I, I1, J, J1, N1, N2, IINFO 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL CLARFG, CTRMM, CGEMM, XERBLA 00158 * .. 00159 * .. Executable Statements .. 00160 * 00161 INFO = 0 00162 IF( N .LT. 0 ) THEN 00163 INFO = -2 00164 ELSE IF( M .LT. N ) THEN 00165 INFO = -1 00166 ELSE IF( LDA .LT. MAX( 1, M ) ) THEN 00167 INFO = -4 00168 ELSE IF( LDT .LT. MAX( 1, N ) ) THEN 00169 INFO = -6 00170 END IF 00171 IF( INFO.NE.0 ) THEN 00172 CALL XERBLA( 'CGEQRT3', -INFO ) 00173 RETURN 00174 END IF 00175 * 00176 IF( N.EQ.1 ) THEN 00177 * 00178 * Compute Householder transform when N=1 00179 * 00180 CALL CLARFG( M, A, A( MIN( 2, M ), 1 ), 1, T ) 00181 * 00182 ELSE 00183 * 00184 * Otherwise, split A into blocks... 00185 * 00186 N1 = N/2 00187 N2 = N-N1 00188 J1 = MIN( N1+1, N ) 00189 I1 = MIN( N+1, M ) 00190 * 00191 * Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1**H 00192 * 00193 CALL CGEQRT3( M, N1, A, LDA, T, LDT, IINFO ) 00194 * 00195 * Compute A(1:M,J1:N) = Q1**H A(1:M,J1:N) [workspace: T(1:N1,J1:N)] 00196 * 00197 DO J=1,N2 00198 DO I=1,N1 00199 T( I, J+N1 ) = A( I, J+N1 ) 00200 END DO 00201 END DO 00202 CALL CTRMM( 'L', 'L', 'C', 'U', N1, N2, ONE, 00203 & A, LDA, T( 1, J1 ), LDT ) 00204 * 00205 CALL CGEMM( 'C', 'N', N1, N2, M-N1, ONE, A( J1, 1 ), LDA, 00206 & A( J1, J1 ), LDA, ONE, T( 1, J1 ), LDT) 00207 * 00208 CALL CTRMM( 'L', 'U', 'C', 'N', N1, N2, ONE, 00209 & T, LDT, T( 1, J1 ), LDT ) 00210 * 00211 CALL CGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( J1, 1 ), LDA, 00212 & T( 1, J1 ), LDT, ONE, A( J1, J1 ), LDA ) 00213 * 00214 CALL CTRMM( 'L', 'L', 'N', 'U', N1, N2, ONE, 00215 & A, LDA, T( 1, J1 ), LDT ) 00216 * 00217 DO J=1,N2 00218 DO I=1,N1 00219 A( I, J+N1 ) = A( I, J+N1 ) - T( I, J+N1 ) 00220 END DO 00221 END DO 00222 * 00223 * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2**H 00224 * 00225 CALL CGEQRT3( M-N1, N2, A( J1, J1 ), LDA, 00226 & T( J1, J1 ), LDT, IINFO ) 00227 * 00228 * Compute T3 = T(1:N1,J1:N) = -T1 Y1**H Y2 T2 00229 * 00230 DO I=1,N1 00231 DO J=1,N2 00232 T( I, J+N1 ) = CONJG(A( J+N1, I )) 00233 END DO 00234 END DO 00235 * 00236 CALL CTRMM( 'R', 'L', 'N', 'U', N1, N2, ONE, 00237 & A( J1, J1 ), LDA, T( 1, J1 ), LDT ) 00238 * 00239 CALL CGEMM( 'C', 'N', N1, N2, M-N, ONE, A( I1, 1 ), LDA, 00240 & A( I1, J1 ), LDA, ONE, T( 1, J1 ), LDT ) 00241 * 00242 CALL CTRMM( 'L', 'U', 'N', 'N', N1, N2, -ONE, T, LDT, 00243 & T( 1, J1 ), LDT ) 00244 * 00245 CALL CTRMM( 'R', 'U', 'N', 'N', N1, N2, ONE, 00246 & T( J1, J1 ), LDT, T( 1, J1 ), LDT ) 00247 * 00248 * Y = (Y1,Y2); R = [ R1 A(1:N1,J1:N) ]; T = [T1 T3] 00249 * [ 0 R2 ] [ 0 T2] 00250 * 00251 END IF 00252 * 00253 RETURN 00254 * 00255 * End of CGEQRT3 00256 * 00257 END