LAPACK 3.3.0

stzrzf.f

Go to the documentation of this file.
00001       SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, LDA, LWORK, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       REAL               A( LDA, * ), TAU( * ), WORK( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
00019 *  to upper triangular form by means of orthogonal transformations.
00020 *
00021 *  The upper trapezoidal matrix A is factored as
00022 *
00023 *     A = ( R  0 ) * Z,
00024 *
00025 *  where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
00026 *  triangular matrix.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  M       (input) INTEGER
00032 *          The number of rows of the matrix A.  M >= 0.
00033 *
00034 *  N       (input) INTEGER
00035 *          The number of columns of the matrix A.  N >= M.
00036 *
00037 *  A       (input/output) REAL array, dimension (LDA,N)
00038 *          On entry, the leading M-by-N upper trapezoidal part of the
00039 *          array A must contain the matrix to be factorized.
00040 *          On exit, the leading M-by-M upper triangular part of A
00041 *          contains the upper triangular matrix R, and elements M+1 to
00042 *          N of the first M rows of A, with the array TAU, represent the
00043 *          orthogonal matrix Z as a product of M elementary reflectors.
00044 *
00045 *  LDA     (input) INTEGER
00046 *          The leading dimension of the array A.  LDA >= max(1,M).
00047 *
00048 *  TAU     (output) REAL array, dimension (M)
00049 *          The scalar factors of the elementary reflectors.
00050 *
00051 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00052 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00053 *
00054 *  LWORK   (input) INTEGER
00055 *          The dimension of the array WORK.  LWORK >= max(1,M).
00056 *          For optimum performance LWORK >= M*NB, where NB is
00057 *          the optimal blocksize.
00058 *
00059 *          If LWORK = -1, then a workspace query is assumed; the routine
00060 *          only calculates the optimal size of the WORK array, returns
00061 *          this value as the first entry of the WORK array, and no error
00062 *          message related to LWORK is issued by XERBLA.
00063 *
00064 *  INFO    (output) INTEGER
00065 *          = 0:  successful exit
00066 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00067 *
00068 *  Further Details
00069 *  ===============
00070 *
00071 *  Based on contributions by
00072 *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00073 *
00074 *  The factorization is obtained by Householder's method.  The kth
00075 *  transformation matrix, Z( k ), which is used to introduce zeros into
00076 *  the ( m - k + 1 )th row of A, is given in the form
00077 *
00078 *     Z( k ) = ( I     0   ),
00079 *              ( 0  T( k ) )
00080 *
00081 *  where
00082 *
00083 *     T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
00084 *                                                 (   0    )
00085 *                                                 ( z( k ) )
00086 *
00087 *  tau is a scalar and z( k ) is an ( n - m ) element vector.
00088 *  tau and z( k ) are chosen to annihilate the elements of the kth row
00089 *  of X.
00090 *
00091 *  The scalar tau is returned in the kth element of TAU and the vector
00092 *  u( k ) in the kth row of A, such that the elements of z( k ) are
00093 *  in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
00094 *  the upper triangular part of A.
00095 *
00096 *  Z is given by
00097 *
00098 *     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
00099 *
00100 *  =====================================================================
00101 *
00102 *     .. Parameters ..
00103       REAL               ZERO
00104       PARAMETER          ( ZERO = 0.0E+0 )
00105 *     ..
00106 *     .. Local Scalars ..
00107       LOGICAL            LQUERY
00108       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKOPT, M1, MU, NB,
00109      $                   NBMIN, NX
00110 *     ..
00111 *     .. External Subroutines ..
00112       EXTERNAL           SLARZB, SLARZT, SLATRZ, XERBLA
00113 *     ..
00114 *     .. Intrinsic Functions ..
00115       INTRINSIC          MAX, MIN
00116 *     ..
00117 *     .. External Functions ..
00118       INTEGER            ILAENV
00119       EXTERNAL           ILAENV
00120 *     ..
00121 *     .. Executable Statements ..
00122 *
00123 *     Test the input arguments
00124 *
00125       INFO = 0
00126       LQUERY = ( LWORK.EQ.-1 )
00127       IF( M.LT.0 ) THEN
00128          INFO = -1
00129       ELSE IF( N.LT.M ) THEN
00130          INFO = -2
00131       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00132          INFO = -4
00133       END IF
00134 *
00135       IF( INFO.EQ.0 ) THEN
00136          IF( M.EQ.0 .OR. M.EQ.N ) THEN
00137             LWKOPT = 1
00138          ELSE
00139 *
00140 *           Determine the block size.
00141 *
00142             NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
00143             LWKOPT = M*NB
00144          END IF
00145          WORK( 1 ) = LWKOPT
00146 *
00147          IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
00148             INFO = -7
00149          END IF
00150       END IF
00151 *
00152       IF( INFO.NE.0 ) THEN
00153          CALL XERBLA( 'STZRZF', -INFO )
00154          RETURN
00155       ELSE IF( LQUERY ) THEN
00156          RETURN
00157       END IF
00158 *
00159 *     Quick return if possible
00160 *
00161       IF( M.EQ.0 ) THEN
00162          RETURN
00163       ELSE IF( M.EQ.N ) THEN
00164          DO 10 I = 1, N
00165             TAU( I ) = ZERO
00166    10    CONTINUE
00167          RETURN
00168       END IF
00169 *
00170       NBMIN = 2
00171       NX = 1
00172       IWS = M
00173       IF( NB.GT.1 .AND. NB.LT.M ) THEN
00174 *
00175 *        Determine when to cross over from blocked to unblocked code.
00176 *
00177          NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) )
00178          IF( NX.LT.M ) THEN
00179 *
00180 *           Determine if workspace is large enough for blocked code.
00181 *
00182             LDWORK = M
00183             IWS = LDWORK*NB
00184             IF( LWORK.LT.IWS ) THEN
00185 *
00186 *              Not enough workspace to use optimal NB:  reduce NB and
00187 *              determine the minimum value of NB.
00188 *
00189                NB = LWORK / LDWORK
00190                NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1,
00191      $                 -1 ) )
00192             END IF
00193          END IF
00194       END IF
00195 *
00196       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
00197 *
00198 *        Use blocked code initially.
00199 *        The last kk rows are handled by the block method.
00200 *
00201          M1 = MIN( M+1, N )
00202          KI = ( ( M-NX-1 ) / NB )*NB
00203          KK = MIN( M, KI+NB )
00204 *
00205          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
00206             IB = MIN( M-I+1, NB )
00207 *
00208 *           Compute the TZ factorization of the current block
00209 *           A(i:i+ib-1,i:n)
00210 *
00211             CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
00212      $                   WORK )
00213             IF( I.GT.1 ) THEN
00214 *
00215 *              Form the triangular factor of the block reflector
00216 *              H = H(i+ib-1) . . . H(i+1) H(i)
00217 *
00218                CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
00219      $                      LDA, TAU( I ), WORK, LDWORK )
00220 *
00221 *              Apply H to A(1:i-1,i:n) from the right
00222 *
00223                CALL SLARZB( 'Right', 'No transpose', 'Backward',
00224      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
00225      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
00226      $                      WORK( IB+1 ), LDWORK )
00227             END IF
00228    20    CONTINUE
00229          MU = I + NB - 1
00230       ELSE
00231          MU = M
00232       END IF
00233 *
00234 *     Use unblocked code to factor the last or only block
00235 *
00236       IF( MU.GT.0 )
00237      $   CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
00238 *
00239       WORK( 1 ) = LWKOPT
00240 *
00241       RETURN
00242 *
00243 *     End of STZRZF
00244 *
00245       END
 All Files Functions