LAPACK: SRC/sgerqf.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGERQF( 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 *  SGERQF computes an RQ factorization of a real M-by-N matrix A:
00019 *  A = R * Q.
00020 *
00021 *  Arguments
00022 *  =========
00023 *
00024 *  M       (input) INTEGER
00025 *          The number of rows of the matrix A.  M >= 0.
00026 *
00027 *  N       (input) INTEGER
00028 *          The number of columns of the matrix A.  N >= 0.
00029 *
00030 *  A       (input/output) REAL array, dimension (LDA,N)
00031 *          On entry, the M-by-N matrix A.
00032 *          On exit,
00033 *          if m <= n, the upper triangle of the subarray
00034 *          A(1:m,n-m+1:n) contains the M-by-M upper triangular matrix R;
00035 *          if m >= n, the elements on and above the (m-n)-th subdiagonal
00036 *          contain the M-by-N upper trapezoidal matrix R;
00037 *          the remaining elements, with the array TAU, represent the
00038 *          orthogonal matrix Q as a product of min(m,n) elementary
00039 *          reflectors (see Further Details).
00040 *
00041 *  LDA     (input) INTEGER
00042 *          The leading dimension of the array A.  LDA >= max(1,M).
00043 *
00044 *  TAU     (output) REAL array, dimension (min(M,N))
00045 *          The scalar factors of the elementary reflectors (see Further
00046 *          Details).
00047 *
00048 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00049 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00050 *
00051 *  LWORK   (input) INTEGER
00052 *          The dimension of the array WORK.  LWORK >= max(1,M).
00053 *          For optimum performance LWORK >= M*NB, where NB is
00054 *          the optimal blocksize.
00055 *
00056 *          If LWORK = -1, then a workspace query is assumed; the routine
00057 *          only calculates the optimal size of the WORK array, returns
00058 *          this value as the first entry of the WORK array, and no error
00059 *          message related to LWORK is issued by XERBLA.
00060 *
00061 *  INFO    (output) INTEGER
00062 *          = 0:  successful exit
00063 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00064 *
00065 *  Further Details
00066 *  ===============
00067 *
00068 *  The matrix Q is represented as a product of elementary reflectors
00069 *
00070 *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
00071 *
00072 *  Each H(i) has the form
00073 *
00074 *     H(i) = I - tau * v * v'
00075 *
00076 *  where tau is a real scalar, and v is a real vector with
00077 *  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
00078 *  A(m-k+i,1:n-k+i-1), and tau in TAU(i).
00079 *
00080 *  =====================================================================
00081 *
00082 *     .. Local Scalars ..
00083       LOGICAL            LQUERY
00084       INTEGER            I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT,
00085      $                   MU, NB, NBMIN, NU, NX
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           SGERQ2, SLARFB, SLARFT, XERBLA
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          MAX, MIN
00092 *     ..
00093 *     .. External Functions ..
00094       INTEGER            ILAENV
00095       EXTERNAL           ILAENV
00096 *     ..
00097 *     .. Executable Statements ..
00098 *
00099 *     Test the input arguments
00100 *
00101       INFO = 0
00102       LQUERY = ( LWORK.EQ.-1 )
00103       IF( M.LT.0 ) THEN
00104          INFO = -1
00105       ELSE IF( N.LT.0 ) THEN
00106          INFO = -2
00107       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00108          INFO = -4
00109       ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
00110          INFO = -7
00111       END IF
00112 *
00113       IF( INFO.EQ.0 ) THEN
00114          K = MIN( M, N )
00115          IF( K.EQ.0 ) THEN
00116             LWKOPT = 1
00117          ELSE
00118             NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
00119             LWKOPT = M*NB
00120             WORK( 1 ) = LWKOPT
00121          END IF
00122          WORK( 1 ) = LWKOPT
00123 *
00124          IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
00125             INFO = -7
00126          END IF
00127       END IF
00128 *
00129       IF( INFO.NE.0 ) THEN
00130          CALL XERBLA( 'SGERQF', -INFO )
00131          RETURN
00132       ELSE IF( LQUERY ) THEN
00133          RETURN
00134       END IF
00135 *
00136 *     Quick return if possible
00137 *
00138       IF( K.EQ.0 ) THEN
00139          RETURN
00140       END IF
00141 *
00142       NBMIN = 2
00143       NX = 1
00144       IWS = M
00145       IF( NB.GT.1 .AND. NB.LT.K ) THEN
00146 *
00147 *        Determine when to cross over from blocked to unblocked code.
00148 *
00149          NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) )
00150          IF( NX.LT.K ) THEN
00151 *
00152 *           Determine if workspace is large enough for blocked code.
00153 *
00154             LDWORK = M
00155             IWS = LDWORK*NB
00156             IF( LWORK.LT.IWS ) THEN
00157 *
00158 *              Not enough workspace to use optimal NB:  reduce NB and
00159 *              determine the minimum value of NB.
00160 *
00161                NB = LWORK / LDWORK
00162                NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1,
00163      $                 -1 ) )
00164             END IF
00165          END IF
00166       END IF
00167 *
00168       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
00169 *
00170 *        Use blocked code initially.
00171 *        The last kk rows are handled by the block method.
00172 *
00173          KI = ( ( K-NX-1 ) / NB )*NB
00174          KK = MIN( K, KI+NB )
00175 *
00176          DO 10 I = K - KK + KI + 1, K - KK + 1, -NB
00177             IB = MIN( K-I+1, NB )
00178 *
00179 *           Compute the RQ factorization of the current block
00180 *           A(m-k+i:m-k+i+ib-1,1:n-k+i+ib-1)
00181 *
00182             CALL SGERQ2( IB, N-K+I+IB-1, A( M-K+I, 1 ), LDA, TAU( I ),
00183      $                   WORK, IINFO )
00184             IF( M-K+I.GT.1 ) THEN
00185 *
00186 *              Form the triangular factor of the block reflector
00187 *              H = H(i+ib-1) . . . H(i+1) H(i)
00188 *
00189                CALL SLARFT( 'Backward', 'Rowwise', N-K+I+IB-1, IB,
00190      $                      A( M-K+I, 1 ), LDA, TAU( I ), WORK, LDWORK )
00191 *
00192 *              Apply H to A(1:m-k+i-1,1:n-k+i+ib-1) from the right
00193 *
00194                CALL SLARFB( 'Right', 'No transpose', 'Backward',
00195      $                      'Rowwise', M-K+I-1, N-K+I+IB-1, IB,
00196      $                      A( M-K+I, 1 ), LDA, WORK, LDWORK, A, LDA,
00197      $                      WORK( IB+1 ), LDWORK )
00198             END IF
00199    10    CONTINUE
00200          MU = M - K + I + NB - 1
00201          NU = N - K + I + NB - 1
00202       ELSE
00203          MU = M
00204          NU = N
00205       END IF
00206 *
00207 *     Use unblocked code to factor the last or only block
00208 *
00209       IF( MU.GT.0 .AND. NU.GT.0 )
00210      $   CALL SGERQ2( MU, NU, A, LDA, TAU, WORK, IINFO )
00211 *
00212       WORK( 1 ) = IWS
00213       RETURN
00214 *
00215 *     End of SGERQF
00216 *
00217       END