summaryrefslogtreecommitdiff
path: root/src/lib/lapack/dgeqlf.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/lapack/dgeqlf.f')
-rw-r--r--src/lib/lapack/dgeqlf.f213
1 files changed, 0 insertions, 213 deletions
diff --git a/src/lib/lapack/dgeqlf.f b/src/lib/lapack/dgeqlf.f
deleted file mode 100644
index ec293574..00000000
--- a/src/lib/lapack/dgeqlf.f
+++ /dev/null
@@ -1,213 +0,0 @@
- SUBROUTINE DGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
-*
-* -- LAPACK routine (version 3.1) --
-* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-* November 2006
-*
-* .. Scalar Arguments ..
- INTEGER INFO, LDA, LWORK, M, N
-* ..
-* .. Array Arguments ..
- DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
-* ..
-*
-* Purpose
-* =======
-*
-* DGEQLF computes a QL factorization of a real M-by-N matrix A:
-* A = Q * L.
-*
-* Arguments
-* =========
-*
-* M (input) INTEGER
-* The number of rows of the matrix A. M >= 0.
-*
-* N (input) INTEGER
-* The number of columns of the matrix A. N >= 0.
-*
-* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-* On entry, the M-by-N matrix A.
-* On exit,
-* if m >= n, the lower triangle of the subarray
-* A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L;
-* if m <= n, the elements on and below the (n-m)-th
-* superdiagonal contain the M-by-N lower trapezoidal matrix L;
-* the remaining elements, with the array TAU, represent the
-* orthogonal matrix Q as a product of elementary reflectors
-* (see Further Details).
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,M).
-*
-* TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
-* The scalar factors of the elementary reflectors (see Further
-* Details).
-*
-* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
-* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
-*
-* LWORK (input) INTEGER
-* The dimension of the array WORK. LWORK >= max(1,N).
-* For optimum performance LWORK >= N*NB, where NB is the
-* optimal blocksize.
-*
-* If LWORK = -1, then a workspace query is assumed; the routine
-* only calculates the optimal size of the WORK array, returns
-* this value as the first entry of the WORK array, and no error
-* message related to LWORK is issued by XERBLA.
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -i, the i-th argument had an illegal value
-*
-* Further Details
-* ===============
-*
-* The matrix Q is represented as a product of elementary reflectors
-*
-* Q = H(k) . . . H(2) H(1), where k = min(m,n).
-*
-* Each H(i) has the form
-*
-* H(i) = I - tau * v * v'
-*
-* where tau is a real scalar, and v is a real vector with
-* v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in
-* A(1:m-k+i-1,n-k+i), and tau in TAU(i).
-*
-* =====================================================================
-*
-* .. Local Scalars ..
- LOGICAL LQUERY
- INTEGER I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT,
- $ MU, NB, NBMIN, NU, NX
-* ..
-* .. External Subroutines ..
- EXTERNAL DGEQL2, DLARFB, DLARFT, XERBLA
-* ..
-* .. Intrinsic Functions ..
- INTRINSIC MAX, MIN
-* ..
-* .. External Functions ..
- INTEGER ILAENV
- EXTERNAL ILAENV
-* ..
-* .. Executable Statements ..
-*
-* Test the input arguments
-*
- INFO = 0
- LQUERY = ( LWORK.EQ.-1 )
- IF( M.LT.0 ) THEN
- INFO = -1
- ELSE IF( N.LT.0 ) THEN
- INFO = -2
- ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
- INFO = -4
- END IF
-*
- IF( INFO.EQ.0 ) THEN
- K = MIN( M, N )
- IF( K.EQ.0 ) THEN
- LWKOPT = 1
- ELSE
- NB = ILAENV( 1, 'DGEQLF', ' ', M, N, -1, -1 )
- LWKOPT = N*NB
- END IF
- WORK( 1 ) = LWKOPT
-*
- IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
- INFO = -7
- END IF
- END IF
-*
- IF( INFO.NE.0 ) THEN
- CALL XERBLA( 'DGEQLF', -INFO )
- RETURN
- ELSE IF( LQUERY ) THEN
- RETURN
- END IF
-*
-* Quick return if possible
-*
- IF( K.EQ.0 ) THEN
- RETURN
- END IF
-*
- NBMIN = 2
- NX = 1
- IWS = N
- IF( NB.GT.1 .AND. NB.LT.K ) THEN
-*
-* Determine when to cross over from blocked to unblocked code.
-*
- NX = MAX( 0, ILAENV( 3, 'DGEQLF', ' ', M, N, -1, -1 ) )
- IF( NX.LT.K ) THEN
-*
-* Determine if workspace is large enough for blocked code.
-*
- LDWORK = N
- IWS = LDWORK*NB
- IF( LWORK.LT.IWS ) THEN
-*
-* Not enough workspace to use optimal NB: reduce NB and
-* determine the minimum value of NB.
-*
- NB = LWORK / LDWORK
- NBMIN = MAX( 2, ILAENV( 2, 'DGEQLF', ' ', M, N, -1,
- $ -1 ) )
- END IF
- END IF
- END IF
-*
- IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
-*
-* Use blocked code initially.
-* The last kk columns are handled by the block method.
-*
- KI = ( ( K-NX-1 ) / NB )*NB
- KK = MIN( K, KI+NB )
-*
- DO 10 I = K - KK + KI + 1, K - KK + 1, -NB
- IB = MIN( K-I+1, NB )
-*
-* Compute the QL factorization of the current block
-* A(1:m-k+i+ib-1,n-k+i:n-k+i+ib-1)
-*
- CALL DGEQL2( M-K+I+IB-1, IB, A( 1, N-K+I ), LDA, TAU( I ),
- $ WORK, IINFO )
- IF( N-K+I.GT.1 ) THEN
-*
-* Form the triangular factor of the block reflector
-* H = H(i+ib-1) . . . H(i+1) H(i)
-*
- CALL DLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB,
- $ A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK )
-*
-* Apply H' to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
-*
- CALL DLARFB( 'Left', 'Transpose', 'Backward',
- $ 'Columnwise', M-K+I+IB-1, N-K+I-1, IB,
- $ A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA,
- $ WORK( IB+1 ), LDWORK )
- END IF
- 10 CONTINUE
- MU = M - K + I + NB - 1
- NU = N - K + I + NB - 1
- ELSE
- MU = M
- NU = N
- END IF
-*
-* Use unblocked code to factor the last or only block
-*
- IF( MU.GT.0 .AND. NU.GT.0 )
- $ CALL DGEQL2( MU, NU, A, LDA, TAU, WORK, IINFO )
-*
- WORK( 1 ) = IWS
- RETURN
-*
-* End of DGEQLF
-*
- END