[22] | 1 | *> \brief \b STZRZF |
---|
| 2 | * |
---|
| 3 | * =========== DOCUMENTATION =========== |
---|
| 4 | * |
---|
| 5 | * Online html documentation available at |
---|
| 6 | * http://www.netlib.org/lapack/explore-html/ |
---|
| 7 | * |
---|
| 8 | *> \htmlonly |
---|
| 9 | *> Download STZRZF + dependencies |
---|
| 10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
| 11 | *> [TGZ]</a> |
---|
| 12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
| 13 | *> [ZIP]</a> |
---|
| 14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
| 15 | *> [TXT]</a> |
---|
| 16 | *> \endhtmlonly |
---|
| 17 | * |
---|
| 18 | * Definition: |
---|
| 19 | * =========== |
---|
| 20 | * |
---|
| 21 | * SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) |
---|
| 22 | * |
---|
| 23 | * .. Scalar Arguments .. |
---|
| 24 | * INTEGER INFO, LDA, LWORK, M, N |
---|
| 25 | * .. |
---|
| 26 | * .. Array Arguments .. |
---|
| 27 | * REAL A( LDA, * ), TAU( * ), WORK( * ) |
---|
| 28 | * .. |
---|
| 29 | * |
---|
| 30 | * |
---|
| 31 | *> \par Purpose: |
---|
| 32 | * ============= |
---|
| 33 | *> |
---|
| 34 | *> \verbatim |
---|
| 35 | *> |
---|
| 36 | *> STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A |
---|
| 37 | *> to upper triangular form by means of orthogonal transformations. |
---|
| 38 | *> |
---|
| 39 | *> The upper trapezoidal matrix A is factored as |
---|
| 40 | *> |
---|
| 41 | *> A = ( R 0 ) * Z, |
---|
| 42 | *> |
---|
| 43 | *> where Z is an N-by-N orthogonal matrix and R is an M-by-M upper |
---|
| 44 | *> triangular matrix. |
---|
| 45 | *> \endverbatim |
---|
| 46 | * |
---|
| 47 | * Arguments: |
---|
| 48 | * ========== |
---|
| 49 | * |
---|
| 50 | *> \param[in] M |
---|
| 51 | *> \verbatim |
---|
| 52 | *> M is INTEGER |
---|
| 53 | *> The number of rows of the matrix A. M >= 0. |
---|
| 54 | *> \endverbatim |
---|
| 55 | *> |
---|
| 56 | *> \param[in] N |
---|
| 57 | *> \verbatim |
---|
| 58 | *> N is INTEGER |
---|
| 59 | *> The number of columns of the matrix A. N >= M. |
---|
| 60 | *> \endverbatim |
---|
| 61 | *> |
---|
| 62 | *> \param[in,out] A |
---|
| 63 | *> \verbatim |
---|
| 64 | *> A is REAL array, dimension (LDA,N) |
---|
| 65 | *> On entry, the leading M-by-N upper trapezoidal part of the |
---|
| 66 | *> array A must contain the matrix to be factorized. |
---|
| 67 | *> On exit, the leading M-by-M upper triangular part of A |
---|
| 68 | *> contains the upper triangular matrix R, and elements M+1 to |
---|
| 69 | *> N of the first M rows of A, with the array TAU, represent the |
---|
| 70 | *> orthogonal matrix Z as a product of M elementary reflectors. |
---|
| 71 | *> \endverbatim |
---|
| 72 | *> |
---|
| 73 | *> \param[in] LDA |
---|
| 74 | *> \verbatim |
---|
| 75 | *> LDA is INTEGER |
---|
| 76 | *> The leading dimension of the array A. LDA >= max(1,M). |
---|
| 77 | *> \endverbatim |
---|
| 78 | *> |
---|
| 79 | *> \param[out] TAU |
---|
| 80 | *> \verbatim |
---|
| 81 | *> TAU is REAL array, dimension (M) |
---|
| 82 | *> The scalar factors of the elementary reflectors. |
---|
| 83 | *> \endverbatim |
---|
| 84 | *> |
---|
| 85 | *> \param[out] WORK |
---|
| 86 | *> \verbatim |
---|
| 87 | *> WORK is REAL array, dimension (MAX(1,LWORK)) |
---|
| 88 | *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
---|
| 89 | *> \endverbatim |
---|
| 90 | *> |
---|
| 91 | *> \param[in] LWORK |
---|
| 92 | *> \verbatim |
---|
| 93 | *> LWORK is INTEGER |
---|
| 94 | *> The dimension of the array WORK. LWORK >= max(1,M). |
---|
| 95 | *> For optimum performance LWORK >= M*NB, where NB is |
---|
| 96 | *> the optimal blocksize. |
---|
| 97 | *> |
---|
| 98 | *> If LWORK = -1, then a workspace query is assumed; the routine |
---|
| 99 | *> only calculates the optimal size of the WORK array, returns |
---|
| 100 | *> this value as the first entry of the WORK array, and no error |
---|
| 101 | *> message related to LWORK is issued by XERBLA. |
---|
| 102 | *> \endverbatim |
---|
| 103 | *> |
---|
| 104 | *> \param[out] INFO |
---|
| 105 | *> \verbatim |
---|
| 106 | *> INFO is INTEGER |
---|
| 107 | *> = 0: successful exit |
---|
| 108 | *> < 0: if INFO = -i, the i-th argument had an illegal value |
---|
| 109 | *> \endverbatim |
---|
| 110 | * |
---|
| 111 | * Authors: |
---|
| 112 | * ======== |
---|
| 113 | * |
---|
| 114 | *> \author Univ. of Tennessee |
---|
| 115 | *> \author Univ. of California Berkeley |
---|
| 116 | *> \author Univ. of Colorado Denver |
---|
| 117 | *> \author NAG Ltd. |
---|
| 118 | * |
---|
| 119 | *> \date November 2011 |
---|
| 120 | * |
---|
| 121 | *> \ingroup realOTHERcomputational |
---|
| 122 | * |
---|
| 123 | *> \par Contributors: |
---|
| 124 | * ================== |
---|
| 125 | *> |
---|
| 126 | *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA |
---|
| 127 | * |
---|
| 128 | *> \par Further Details: |
---|
| 129 | * ===================== |
---|
| 130 | *> |
---|
| 131 | *> \verbatim |
---|
| 132 | *> |
---|
| 133 | *> The factorization is obtained by Householder's method. The kth |
---|
| 134 | *> transformation matrix, Z( k ), which is used to introduce zeros into |
---|
| 135 | *> the ( m - k + 1 )th row of A, is given in the form |
---|
| 136 | *> |
---|
| 137 | *> Z( k ) = ( I 0 ), |
---|
| 138 | *> ( 0 T( k ) ) |
---|
| 139 | *> |
---|
| 140 | *> where |
---|
| 141 | *> |
---|
| 142 | *> T( k ) = I - tau*u( k )*u( k )**T, u( k ) = ( 1 ), |
---|
| 143 | *> ( 0 ) |
---|
| 144 | *> ( z( k ) ) |
---|
| 145 | *> |
---|
| 146 | *> tau is a scalar and z( k ) is an ( n - m ) element vector. |
---|
| 147 | *> tau and z( k ) are chosen to annihilate the elements of the kth row |
---|
| 148 | *> of X. |
---|
| 149 | *> |
---|
| 150 | *> The scalar tau is returned in the kth element of TAU and the vector |
---|
| 151 | *> u( k ) in the kth row of A, such that the elements of z( k ) are |
---|
| 152 | *> in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in |
---|
| 153 | *> the upper triangular part of A. |
---|
| 154 | *> |
---|
| 155 | *> Z is given by |
---|
| 156 | *> |
---|
| 157 | *> Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). |
---|
| 158 | *> \endverbatim |
---|
| 159 | *> |
---|
| 160 | * ===================================================================== |
---|
| 161 | SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) |
---|
| 162 | * |
---|
| 163 | * -- LAPACK computational routine (version 3.4.0) -- |
---|
| 164 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
| 165 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
| 166 | * November 2011 |
---|
| 167 | * |
---|
| 168 | * .. Scalar Arguments .. |
---|
| 169 | INTEGER INFO, LDA, LWORK, M, N |
---|
| 170 | * .. |
---|
| 171 | * .. Array Arguments .. |
---|
| 172 | REAL A( LDA, * ), TAU( * ), WORK( * ) |
---|
| 173 | * .. |
---|
| 174 | * |
---|
| 175 | * ===================================================================== |
---|
| 176 | * |
---|
| 177 | * .. Parameters .. |
---|
| 178 | REAL ZERO |
---|
| 179 | PARAMETER ( ZERO = 0.0E+0 ) |
---|
| 180 | * .. |
---|
| 181 | * .. Local Scalars .. |
---|
| 182 | LOGICAL LQUERY |
---|
| 183 | INTEGER I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT, |
---|
| 184 | $ M1, MU, NB, NBMIN, NX |
---|
| 185 | * .. |
---|
| 186 | * .. External Subroutines .. |
---|
| 187 | EXTERNAL XERBLA, SLARZB, SLARZT, SLATRZ |
---|
| 188 | * .. |
---|
| 189 | * .. Intrinsic Functions .. |
---|
| 190 | INTRINSIC MAX, MIN |
---|
| 191 | * .. |
---|
| 192 | * .. External Functions .. |
---|
| 193 | INTEGER ILAENV |
---|
| 194 | EXTERNAL ILAENV |
---|
| 195 | * .. |
---|
| 196 | * .. Executable Statements .. |
---|
| 197 | * |
---|
| 198 | * Test the input arguments |
---|
| 199 | * |
---|
| 200 | INFO = 0 |
---|
| 201 | LQUERY = ( LWORK.EQ.-1 ) |
---|
| 202 | IF( M.LT.0 ) THEN |
---|
| 203 | INFO = -1 |
---|
| 204 | ELSE IF( N.LT.M ) THEN |
---|
| 205 | INFO = -2 |
---|
| 206 | ELSE IF( LDA.LT.MAX( 1, M ) ) THEN |
---|
| 207 | INFO = -4 |
---|
| 208 | END IF |
---|
| 209 | * |
---|
| 210 | IF( INFO.EQ.0 ) THEN |
---|
| 211 | IF( M.EQ.0 .OR. M.EQ.N ) THEN |
---|
| 212 | LWKOPT = 1 |
---|
| 213 | LWKMIN = 1 |
---|
| 214 | ELSE |
---|
| 215 | * |
---|
| 216 | * Determine the block size. |
---|
| 217 | * |
---|
| 218 | NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 ) |
---|
| 219 | LWKOPT = M*NB |
---|
| 220 | LWKMIN = MAX( 1, M ) |
---|
| 221 | END IF |
---|
| 222 | WORK( 1 ) = LWKOPT |
---|
| 223 | * |
---|
| 224 | IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN |
---|
| 225 | INFO = -7 |
---|
| 226 | END IF |
---|
| 227 | END IF |
---|
| 228 | * |
---|
| 229 | IF( INFO.NE.0 ) THEN |
---|
| 230 | CALL XERBLA( 'STZRZF', -INFO ) |
---|
| 231 | RETURN |
---|
| 232 | ELSE IF( LQUERY ) THEN |
---|
| 233 | RETURN |
---|
| 234 | END IF |
---|
| 235 | * |
---|
| 236 | * Quick return if possible |
---|
| 237 | * |
---|
| 238 | IF( M.EQ.0 ) THEN |
---|
| 239 | RETURN |
---|
| 240 | ELSE IF( M.EQ.N ) THEN |
---|
| 241 | DO 10 I = 1, N |
---|
| 242 | TAU( I ) = ZERO |
---|
| 243 | 10 CONTINUE |
---|
| 244 | RETURN |
---|
| 245 | END IF |
---|
| 246 | * |
---|
| 247 | NBMIN = 2 |
---|
| 248 | NX = 1 |
---|
| 249 | IWS = M |
---|
| 250 | IF( NB.GT.1 .AND. NB.LT.M ) THEN |
---|
| 251 | * |
---|
| 252 | * Determine when to cross over from blocked to unblocked code. |
---|
| 253 | * |
---|
| 254 | NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) ) |
---|
| 255 | IF( NX.LT.M ) THEN |
---|
| 256 | * |
---|
| 257 | * Determine if workspace is large enough for blocked code. |
---|
| 258 | * |
---|
| 259 | LDWORK = M |
---|
| 260 | IWS = LDWORK*NB |
---|
| 261 | IF( LWORK.LT.IWS ) THEN |
---|
| 262 | * |
---|
| 263 | * Not enough workspace to use optimal NB: reduce NB and |
---|
| 264 | * determine the minimum value of NB. |
---|
| 265 | * |
---|
| 266 | NB = LWORK / LDWORK |
---|
| 267 | NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1, |
---|
| 268 | $ -1 ) ) |
---|
| 269 | END IF |
---|
| 270 | END IF |
---|
| 271 | END IF |
---|
| 272 | * |
---|
| 273 | IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN |
---|
| 274 | * |
---|
| 275 | * Use blocked code initially. |
---|
| 276 | * The last kk rows are handled by the block method. |
---|
| 277 | * |
---|
| 278 | M1 = MIN( M+1, N ) |
---|
| 279 | KI = ( ( M-NX-1 ) / NB )*NB |
---|
| 280 | KK = MIN( M, KI+NB ) |
---|
| 281 | * |
---|
| 282 | DO 20 I = M - KK + KI + 1, M - KK + 1, -NB |
---|
| 283 | IB = MIN( M-I+1, NB ) |
---|
| 284 | * |
---|
| 285 | * Compute the TZ factorization of the current block |
---|
| 286 | * A(i:i+ib-1,i:n) |
---|
| 287 | * |
---|
| 288 | CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ), |
---|
| 289 | $ WORK ) |
---|
| 290 | IF( I.GT.1 ) THEN |
---|
| 291 | * |
---|
| 292 | * Form the triangular factor of the block reflector |
---|
| 293 | * H = H(i+ib-1) . . . H(i+1) H(i) |
---|
| 294 | * |
---|
| 295 | CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ), |
---|
| 296 | $ LDA, TAU( I ), WORK, LDWORK ) |
---|
| 297 | * |
---|
| 298 | * Apply H to A(1:i-1,i:n) from the right |
---|
| 299 | * |
---|
| 300 | CALL SLARZB( 'Right', 'No transpose', 'Backward', |
---|
| 301 | $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ), |
---|
| 302 | $ LDA, WORK, LDWORK, A( 1, I ), LDA, |
---|
| 303 | $ WORK( IB+1 ), LDWORK ) |
---|
| 304 | END IF |
---|
| 305 | 20 CONTINUE |
---|
| 306 | MU = I + NB - 1 |
---|
| 307 | ELSE |
---|
| 308 | MU = M |
---|
| 309 | END IF |
---|
| 310 | * |
---|
| 311 | * Use unblocked code to factor the last or only block |
---|
| 312 | * |
---|
| 313 | IF( MU.GT.0 ) |
---|
| 314 | $ CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK ) |
---|
| 315 | * |
---|
| 316 | WORK( 1 ) = LWKOPT |
---|
| 317 | * |
---|
| 318 | RETURN |
---|
| 319 | * |
---|
| 320 | * End of STZRZF |
---|
| 321 | * |
---|
| 322 | END |
---|