Go to the documentation of this file.
00001 *> \brief \b STZRZF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STZRZF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stzrzf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stzrzf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stzrzf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, LWORK, M, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       REAL               A( LDA, * ), TAU( * ), WORK( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
00037 *> to upper triangular form by means of orthogonal transformations.
00038 *>
00039 *> The upper trapezoidal matrix A is factored as
00040 *>
00041 *>    A = ( R  0 ) * Z,
00042 *>
00043 *> where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
00044 *> triangular matrix.
00045 *> \endverbatim
00046 *
00047 *  Arguments:
00048 *  ==========
00049 *
00050 *> \param[in] M
00051 *> \verbatim
00052 *>          M is INTEGER
00053 *>          The number of rows of the matrix A.  M >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] N
00057 *> \verbatim
00058 *>          N is INTEGER
00059 *>          The number of columns of the matrix A.  N >= M.
00060 *> \endverbatim
00061 *>
00062 *> \param[in,out] A
00063 *> \verbatim
00064 *>          A is REAL array, dimension (LDA,N)
00065 *>          On entry, the leading M-by-N upper trapezoidal part of the
00066 *>          array A must contain the matrix to be factorized.
00067 *>          On exit, the leading M-by-M upper triangular part of A
00068 *>          contains the upper triangular matrix R, and elements M+1 to
00069 *>          N of the first M rows of A, with the array TAU, represent the
00070 *>          orthogonal matrix Z as a product of M elementary reflectors.
00071 *> \endverbatim
00072 *>
00073 *> \param[in] LDA
00074 *> \verbatim
00075 *>          LDA is INTEGER
00076 *>          The leading dimension of the array A.  LDA >= max(1,M).
00077 *> \endverbatim
00078 *>
00079 *> \param[out] TAU
00080 *> \verbatim
00081 *>          TAU is REAL array, dimension (M)
00082 *>          The scalar factors of the elementary reflectors.
00083 *> \endverbatim
00084 *>
00085 *> \param[out] WORK
00086 *> \verbatim
00087 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00088 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LWORK
00092 *> \verbatim
00093 *>          LWORK is INTEGER
00094 *>          The dimension of the array WORK.  LWORK >= max(1,M).
00095 *>          For optimum performance LWORK >= M*NB, where NB is
00096 *>          the optimal blocksize.
00097 *>
00098 *>          If LWORK = -1, then a workspace query is assumed; the routine
00099 *>          only calculates the optimal size of the WORK array, returns
00100 *>          this value as the first entry of the WORK array, and no error
00101 *>          message related to LWORK is issued by XERBLA.
00102 *> \endverbatim
00103 *>
00104 *> \param[out] INFO
00105 *> \verbatim
00106 *>          INFO is INTEGER
00107 *>          = 0:  successful exit
00108 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00109 *> \endverbatim
00110 *
00111 *  Authors:
00112 *  ========
00113 *
00114 *> \author Univ. of Tennessee 
00115 *> \author Univ. of California Berkeley 
00116 *> \author Univ. of Colorado Denver 
00117 *> \author NAG Ltd. 
00118 *
00119 *> \date November 2011
00120 *
00121 *> \ingroup realOTHERcomputational
00122 *
00123 *> \par Contributors:
00124 *  ==================
00125 *>
00126 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00127 *
00128 *> \par Further Details:
00129 *  =====================
00130 *>
00131 *> \verbatim
00132 *>
00133 *>  The factorization is obtained by Householder's method.  The kth
00134 *>  transformation matrix, Z( k ), which is used to introduce zeros into
00135 *>  the ( m - k + 1 )th row of A, is given in the form
00136 *>
00137 *>     Z( k ) = ( I     0   ),
00138 *>              ( 0  T( k ) )
00139 *>
00140 *>  where
00141 *>
00142 *>     T( k ) = I - tau*u( k )*u( k )**T,   u( k ) = (   1    ),
00143 *>                                                 (   0    )
00144 *>                                                 ( z( k ) )
00145 *>
00146 *>  tau is a scalar and z( k ) is an ( n - m ) element vector.
00147 *>  tau and z( k ) are chosen to annihilate the elements of the kth row
00148 *>  of X.
00149 *>
00150 *>  The scalar tau is returned in the kth element of TAU and the vector
00151 *>  u( k ) in the kth row of A, such that the elements of z( k ) are
00152 *>  in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
00153 *>  the upper triangular part of A.
00154 *>
00155 *>  Z is given by
00156 *>
00157 *>     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
00158 *> \endverbatim
00159 *>
00160 *  =====================================================================
00161       SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
00162 *
00163 *  -- LAPACK computational routine (version 3.4.0) --
00164 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00165 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00166 *     November 2011
00167 *
00168 *     .. Scalar Arguments ..
00169       INTEGER            INFO, LDA, LWORK, M, N
00170 *     ..
00171 *     .. Array Arguments ..
00172       REAL               A( LDA, * ), TAU( * ), WORK( * )
00173 *     ..
00174 *
00175 *  =====================================================================
00176 *
00177 *     .. Parameters ..
00178       REAL               ZERO
00179       PARAMETER          ( ZERO = 0.0E+0 )
00180 *     ..
00181 *     .. Local Scalars ..
00182       LOGICAL            LQUERY
00183       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
00184      $                   M1, MU, NB, NBMIN, NX
00185 *     ..
00186 *     .. External Subroutines ..
00187       EXTERNAL           XERBLA, SLARZB, SLARZT, SLATRZ
00188 *     ..
00189 *     .. Intrinsic Functions ..
00190       INTRINSIC          MAX, MIN
00191 *     ..
00192 *     .. External Functions ..
00193       INTEGER            ILAENV
00194       EXTERNAL           ILAENV
00195 *     ..
00196 *     .. Executable Statements ..
00197 *
00198 *     Test the input arguments
00199 *
00200       INFO = 0
00201       LQUERY = ( LWORK.EQ.-1 )
00202       IF( M.LT.0 ) THEN
00203          INFO = -1
00204       ELSE IF( N.LT.M ) THEN
00205          INFO = -2
00206       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00207          INFO = -4
00208       END IF
00209 *
00210       IF( INFO.EQ.0 ) THEN
00211          IF( M.EQ.0 .OR. M.EQ.N ) THEN
00212             LWKOPT = 1
00213             LWKMIN = 1
00214          ELSE
00215 *
00216 *           Determine the block size.
00217 *
00218             NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
00219             LWKOPT = M*NB
00220             LWKMIN = MAX( 1, M )
00221          END IF
00222          WORK( 1 ) = LWKOPT
00223 *
00224          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00225             INFO = -7
00226          END IF
00227       END IF
00228 *
00229       IF( INFO.NE.0 ) THEN
00230          CALL XERBLA( 'STZRZF', -INFO )
00231          RETURN
00232       ELSE IF( LQUERY ) THEN
00233          RETURN
00234       END IF
00235 *
00236 *     Quick return if possible
00237 *
00238       IF( M.EQ.0 ) THEN
00239          RETURN
00240       ELSE IF( M.EQ.N ) THEN
00241          DO 10 I = 1, N
00242             TAU( I ) = ZERO
00243    10    CONTINUE
00244          RETURN
00245       END IF
00246 *
00247       NBMIN = 2
00248       NX = 1
00249       IWS = M
00250       IF( NB.GT.1 .AND. NB.LT.M ) THEN
00251 *
00252 *        Determine when to cross over from blocked to unblocked code.
00253 *
00254          NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) )
00255          IF( NX.LT.M ) THEN
00256 *
00257 *           Determine if workspace is large enough for blocked code.
00258 *
00259             LDWORK = M
00260             IWS = LDWORK*NB
00261             IF( LWORK.LT.IWS ) THEN
00262 *
00263 *              Not enough workspace to use optimal NB:  reduce NB and
00264 *              determine the minimum value of NB.
00265 *
00266                NB = LWORK / LDWORK
00267                NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1,
00268      $                 -1 ) )
00269             END IF
00270          END IF
00271       END IF
00272 *
00273       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
00274 *
00275 *        Use blocked code initially.
00276 *        The last kk rows are handled by the block method.
00277 *
00278          M1 = MIN( M+1, N )
00279          KI = ( ( M-NX-1 ) / NB )*NB
00280          KK = MIN( M, KI+NB )
00281 *
00282          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
00283             IB = MIN( M-I+1, NB )
00284 *
00285 *           Compute the TZ factorization of the current block
00286 *           A(i:i+ib-1,i:n)
00287 *
00288             CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
00289      $                   WORK )
00290             IF( I.GT.1 ) THEN
00291 *
00292 *              Form the triangular factor of the block reflector
00293 *              H = H(i+ib-1) . . . H(i+1) H(i)
00294 *
00295                CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
00296      $                      LDA, TAU( I ), WORK, LDWORK )
00297 *
00298 *              Apply H to A(1:i-1,i:n) from the right
00299 *
00300                CALL SLARZB( 'Right', 'No transpose', 'Backward',
00301      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
00302      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
00303      $                      WORK( IB+1 ), LDWORK )
00304             END IF
00305    20    CONTINUE
00306          MU = I + NB - 1
00307       ELSE
00308          MU = M
00309       END IF
00310 *
00311 *     Use unblocked code to factor the last or only block
00312 *
00313       IF( MU.GT.0 )
00314      $   CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
00315 *
00316       WORK( 1 ) = LWKOPT
00317 *
00318       RETURN
00319 *
00320 *     End of STZRZF
00321 *
00322       END