LAPACK 3.3.0

sgelqf.f

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