LAPACK 3.3.0

dgerqf.f

Go to the documentation of this file.
00001       SUBROUTINE DGERQF( 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       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DGERQF 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (min(M,N))
00045 *          The scalar factors of the elementary reflectors (see Further
00046 *          Details).
00047 *
00048 *  WORK    (workspace/output) DOUBLE PRECISION 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           DGERQ2, DLARFB, DLARFT, 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       END IF
00110 *
00111       IF( INFO.EQ.0 ) THEN
00112          K = MIN( M, N )
00113          IF( K.EQ.0 ) THEN
00114             LWKOPT = 1
00115          ELSE
00116             NB = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
00117             LWKOPT = M*NB
00118          END IF
00119          WORK( 1 ) = LWKOPT
00120 *
00121          IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
00122             INFO = -7
00123          END IF
00124       END IF
00125 *
00126       IF( INFO.NE.0 ) THEN
00127          CALL XERBLA( 'DGERQF', -INFO )
00128          RETURN
00129       ELSE IF( LQUERY ) THEN
00130          RETURN
00131       END IF
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( K.EQ.0 ) THEN
00136          RETURN
00137       END IF
00138 *
00139       NBMIN = 2
00140       NX = 1
00141       IWS = M
00142       IF( NB.GT.1 .AND. NB.LT.K ) THEN
00143 *
00144 *        Determine when to cross over from blocked to unblocked code.
00145 *
00146          NX = MAX( 0, ILAENV( 3, 'DGERQF', ' ', M, N, -1, -1 ) )
00147          IF( NX.LT.K ) THEN
00148 *
00149 *           Determine if workspace is large enough for blocked code.
00150 *
00151             LDWORK = M
00152             IWS = LDWORK*NB
00153             IF( LWORK.LT.IWS ) THEN
00154 *
00155 *              Not enough workspace to use optimal NB:  reduce NB and
00156 *              determine the minimum value of NB.
00157 *
00158                NB = LWORK / LDWORK
00159                NBMIN = MAX( 2, ILAENV( 2, 'DGERQF', ' ', M, N, -1,
00160      $                 -1 ) )
00161             END IF
00162          END IF
00163       END IF
00164 *
00165       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
00166 *
00167 *        Use blocked code initially.
00168 *        The last kk rows are handled by the block method.
00169 *
00170          KI = ( ( K-NX-1 ) / NB )*NB
00171          KK = MIN( K, KI+NB )
00172 *
00173          DO 10 I = K - KK + KI + 1, K - KK + 1, -NB
00174             IB = MIN( K-I+1, NB )
00175 *
00176 *           Compute the RQ factorization of the current block
00177 *           A(m-k+i:m-k+i+ib-1,1:n-k+i+ib-1)
00178 *
00179             CALL DGERQ2( IB, N-K+I+IB-1, A( M-K+I, 1 ), LDA, TAU( I ),
00180      $                   WORK, IINFO )
00181             IF( M-K+I.GT.1 ) THEN
00182 *
00183 *              Form the triangular factor of the block reflector
00184 *              H = H(i+ib-1) . . . H(i+1) H(i)
00185 *
00186                CALL DLARFT( 'Backward', 'Rowwise', N-K+I+IB-1, IB,
00187      $                      A( M-K+I, 1 ), LDA, TAU( I ), WORK, LDWORK )
00188 *
00189 *              Apply H to A(1:m-k+i-1,1:n-k+i+ib-1) from the right
00190 *
00191                CALL DLARFB( 'Right', 'No transpose', 'Backward',
00192      $                      'Rowwise', M-K+I-1, N-K+I+IB-1, IB,
00193      $                      A( M-K+I, 1 ), LDA, WORK, LDWORK, A, LDA,
00194      $                      WORK( IB+1 ), LDWORK )
00195             END IF
00196    10    CONTINUE
00197          MU = M - K + I + NB - 1
00198          NU = N - K + I + NB - 1
00199       ELSE
00200          MU = M
00201          NU = N
00202       END IF
00203 *
00204 *     Use unblocked code to factor the last or only block
00205 *
00206       IF( MU.GT.0 .AND. NU.GT.0 )
00207      $   CALL DGERQ2( MU, NU, A, LDA, TAU, WORK, IINFO )
00208 *
00209       WORK( 1 ) = IWS
00210       RETURN
00211 *
00212 *     End of DGERQF
00213 *
00214       END
 All Files Functions