summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorkonstantin.i.arturov <karturov@mugwort.jf.intel.com>2016-12-14 01:02:35 -0800
committerkonstantin.i.arturov <karturov@mugwort.jf.intel.com>2016-12-14 01:02:35 -0800
commitae90d35f0d85e236fdf3fa77d508d8859d568f5b (patch)
tree02e6a596d5510e2d524df51e69a4fba9d861faff
parentb80d7b3365caa35838a86e19615b13d913cb2fed (diff)
downloadlapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.tar.gz
lapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.tar.bz2
lapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.zip
Changes in TS QR API
-rw-r--r--SRC/cgelq.f228
-rw-r--r--SRC/cgemlq.f148
-rw-r--r--SRC/cgemqr.f135
-rw-r--r--SRC/cgeqr.f223
-rw-r--r--SRC/cgetsls.f163
-rw-r--r--SRC/dgelq.f232
-rw-r--r--SRC/dgemlq.f135
-rw-r--r--SRC/dgemqr.f133
-rw-r--r--SRC/dgeqr.f228
-rw-r--r--SRC/dgetsls.f114
-rw-r--r--SRC/sgelq.f228
-rw-r--r--SRC/sgemlq.f140
-rw-r--r--SRC/sgemqr.f132
-rw-r--r--SRC/sgeqr.f222
-rw-r--r--SRC/sgetsls.f121
-rw-r--r--SRC/zgelq.f227
-rw-r--r--SRC/zgemlq.f148
-rw-r--r--SRC/zgemqr.f135
-rw-r--r--SRC/zgeqr.f222
-rw-r--r--SRC/zgetsls.f160
20 files changed, 1924 insertions, 1550 deletions
diff --git a/SRC/cgelq.f b/SRC/cgelq.f
index c6c962d7..0abd2d72 100644
--- a/SRC/cgelq.f
+++ b/SRC/cgelq.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-* INFO)
+* SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* COMPLEX A( LDA, * ), WORK1( * ), WORK2( * )
+* COMPLEX A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> CGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using CLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise CGELQT:
-*> A = L * Q .
+*> CGELQ computes a LQ factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,8 +42,8 @@
*> On exit, the elements on and below the diagonal of the array
*> contain the M-by-min(M,N) lower trapezoidal matrix L
*> (L is lower triangular if M <= N);
-*> the elements above the diagonal are the rows of
-*> blocked V representing Q (see Further Details).
+*> the elements above the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -56,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is COMPLEX array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> CLASWLQ or CGELQT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): horizontal block size
-*> WORK1(5): vertical block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> CLASWLQ or CGELQT
+*> T is COMPLEX array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
-*>
+*> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -114,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLASWLQ or DGELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*> the LQ factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -142,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- COMPLEX A( LDA, * ), WORK1( * ), WORK2( * )
+ COMPLEX A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -173,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -184,33 +218,31 @@
MB = 1
NB = N
END IF
- IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
- IF( NB.GT.N.OR.NB.LE.M) NB = N
- MINLW1 = M + 5
- IF ((NB.GT.M).AND.(N.GT.M)) THEN
- IF(MOD(N-M, NB-M).EQ.0) THEN
- NBLCKS = (N-M)/(NB-M)
+ IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+ IF( NB.GT.N .OR. NB.LE.M ) NB = N
+ MINTSZ = M + 5
+ IF ( NB.GT.M .AND. N.GT.M ) THEN
+ IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+ NBLCKS = ( N - M ) / ( NB - M )
ELSE
- NBLCKS = (N-M)/(NB-M) + 1
+ NBLCKS = ( N - M ) / ( NB - M ) + 1
END IF
ELSE
NBLCKS = 1
END IF
-* Determine if the workspace size satisfies minimum size
+*
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
- $ .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+ $ .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+ $ .AND. ( .NOT.LQUERY) ) THEN
+ IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
MB = 1
- END IF
- IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
- LMINWS = .TRUE.
NB = N
END IF
- IF (LWORK2.LT.MB*M) THEN
+ IF ( LWORK.LT.MB*M ) THEN
LMINWS = .TRUE.
MB = 1
END IF
@@ -222,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
- $ .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS) ) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = MB*M*NBLCKS+5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = MB * M
- WORK2(2) = M
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = MB*M*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, MB*M )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGELQ', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The LQ Decomposition
*
- IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
- CALL CGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+ IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+ CALL CGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
ELSE
- CALL CLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
- $ LWORK2, INFO)
+ CALL CLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+ $ LWORK, INFO)
END IF
+ WORK(1) = MAX( 1, MB*M )
RETURN
*
* End of CGELQ
diff --git a/SRC/cgemlq.f b/SRC/cgemlq.f
index 1a551ca3..03dae76d 100644
--- a/SRC/cgemlq.f
+++ b/SRC/cgemlq.f
@@ -2,17 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* COMPLEX A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* COMPLEX A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -21,12 +20,12 @@
*> CGEMLQ overwrites the general real M-by-N matrix C with
*>
*>
-*> SIDE = 'L' SIDE = 'R'
-*> TRANS = 'N': Q * C C * Q
-*> TRANS = 'T': Q**T * C C * Q**T
-*> where Q is a complex orthogonal matrix defined as the product
-*> of blocked elementary reflectors computed by short wide LQ
-*> factorization (DGELQ)
+*> SIDE = 'L' SIDE = 'R'
+*> TRANS = 'N': Q * C C * Q
+*> TRANS = 'C': Q**H * C C * Q**H
+*> where Q is a complex unitary matrix defined as the product
+*> of blocked elementary reflectors computed by short wide
+*> LQ factorization (CGELQ)
*> \endverbatim
*
* Arguments:
@@ -62,12 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is COMPLEX array, dimension (LDA,K)
-*> The i-th row must contain the vector which defines the blocked
-*> elementary reflector H(i), for i = 1,2,...,k, as returned by
-*> DLASWLQ in the first k rows of its array argument A.
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -78,41 +75,42 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is COMPLEX array, dimension (MAX(1,LWORK1)) is
-*> returned by GEQR.
+*> T is COMPLEX array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
*> C is COMPLEX array, dimension (LDC,N)
*> On entry, the M-by-N matrix C.
*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
*> \param[in] LDC
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -128,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LASWQR or GELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*> the LQ factorization.
+*> This version of GEMLQ will use either LAMSWLQ or GEMLQT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LAMSWLQ or GEMLQT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -157,10 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- COMPLEX A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+ COMPLEX A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -174,7 +180,7 @@
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
- EXTERNAL CLAMSWLQ, CGEMLQT, XERBLA
+ EXTERNAL ZLAMSWLQ, ZGEMLQT, XERBLA
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, MOD
* ..
@@ -182,26 +188,27 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'C' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
- IF (LEFT) THEN
+ MB = INT(T(2))
+ NB = INT(T(3))
+ IF ( LEFT ) THEN
LW = N * MB
MN = M
ELSE
LW = M * MB
MN = N
END IF
- IF ((NB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, NB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(NB-K)
+*
+ IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( NB - K )
ELSE
- NBLCKS = (MN-K)/(NB-K) + 1
+ NBLCKS = ( MN - K ) / ( NB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -220,40 +227,43 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
- INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ INFO = -11
+ ELSE IF(( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = REAL( LW )
END IF
+*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGEMLQ', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
- $ (NB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
CALL CGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
- $ WORK1(6), MB, C, LDC, WORK2, INFO)
+ $ T(4), MB, C, LDC, WORK, INFO)
ELSE
- CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ MB, C, LDC, WORK2, LWORK2, INFO )
+ CALL CLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ MB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = REAL ( LW )
+*
RETURN
*
* End of CGEMLQ
diff --git a/SRC/cgemqr.f b/SRC/cgemqr.f
index 51d38b85..a5d6263f 100644
--- a/SRC/cgemqr.f
+++ b/SRC/cgemqr.f
@@ -2,17 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* COMPLEX A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* COMPLEX A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -23,8 +22,8 @@
*>
*> SIDE = 'L' SIDE = 'R'
*> TRANS = 'N': Q * C C * Q
-*> TRANS = 'T': Q**T * C C * Q**T
-*> where Q is a complex orthogonal matrix defined as the product
+*> TRANS = 'T': Q**H * C C * Q**H
+*> where Q is a complex unitary matrix defined as the product
*> of blocked elementary reflectors computed by tall skinny
*> QR factorization (CGEQR)
*> \endverbatim
@@ -62,13 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is COMPLEX array, dimension (LDA,K)
-*> The i-th column must contain the vector which defines the
-*> blockedelementary reflector H(i), for i = 1,2,...,k, as
-*> returned by DGETSQR in the first k columns of
-*> its array argument A.
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
*> \param[in] LDA
@@ -79,16 +75,16 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is COMPLEX array, dimension (MAX(1,LWORK1)) as
-*> it is returned by GEQR.
+*> T is COMPLEX array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
@@ -100,21 +96,21 @@
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -130,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LATSQR or GEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
+*> This version of GEMQR will use either LAMTSQR or GEMQRT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LATMSQR or GEMQRT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE CGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -159,11 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- COMPLEX A( LDA, * ), WORK1( * ), C(LDC, * ),
- $ WORK2( * )
+ COMPLEX A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -185,14 +188,14 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'C' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
+ MB = INT(T(2))
+ NB = INT(T(3))
IF(LEFT) THEN
LW = N * NB
MN = M
@@ -201,12 +204,12 @@
MN = N
END IF
*
- IF ((MB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, MB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(MB-K)
- ELSE
- NBLCKS = (MN-K)/(MB-K) + 1
- END IF
+ IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( MB - K )
+ ELSE
+ NBLCKS = ( MN - K ) / ( MB - K ) + 1
+ END IF
ELSE
NBLCKS = 1
END IF
@@ -224,24 +227,25 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
- ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+ ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
* Determine the block size if it is tall skinny or short and wide
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGEMQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
@@ -251,16 +255,17 @@
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
- $ (MB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
CALL CGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
- $ WORK1(6), NB, C, LDC, WORK2, INFO)
+ $ T(4), NB, C, LDC, WORK, INFO )
ELSE
- CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ NB, C, LDC, WORK2, LWORK2, INFO )
+ CALL CLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ NB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
+*
RETURN
*
* End of CGEMQR
diff --git a/SRC/cgeqr.f b/SRC/cgeqr.f
index 330fda5c..1743960f 100644
--- a/SRC/cgeqr.f
+++ b/SRC/cgeqr.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE CGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+* SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* COMPLEX A( LDA, * ), WORK1( * ), WORK2( * )
+* COMPLEX A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> CGEQR computes a QR factorization of an M-by-N matrix A,
-*> using CLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise CGEQRT:
-*> A = Q * R .
+*> CGEQR computes a QR factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,7 +42,8 @@
*> On exit, the elements on and above the diagonal of the array
*> contain the min(M,N)-by-N upper trapezoidal matrix R
*> (R is upper triangular if M >= N);
-*> the elements below the diagonal represent Q (see Further Details).
+*> the elements below the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -55,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is COMPLEX array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> CLATSQR or CGEQRT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): row block size
-*> WORK1(5): column block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> CLATSQR or CGEQRT
+*> T is COMPLEX array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -113,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLATSQR or DGEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE CGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- COMPLEX A( LDA, * ), WORK1( * ), WORK2( * )
+ COMPLEX A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -172,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -183,34 +218,31 @@
MB = M
NB = 1
END IF
- IF( MB.GT.M.OR.MB.LE.N) MB = M
- IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
- MINLW1 = N + 5
- IF ((MB.GT.N).AND.(M.GT.N)) THEN
- IF(MOD(M-N, MB-N).EQ.0) THEN
- NBLCKS = (M-N)/(MB-N)
+ IF( MB.GT.M .OR. MB.LE.N ) MB = M
+ IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+ MINTSZ = N + 5
+ IF ( MB.GT.N .AND. M.GT.N ) THEN
+ IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+ NBLCKS = ( M - N ) / ( MB - N )
ELSE
- NBLCKS = (M-N)/(MB-N) + 1
+ NBLCKS = ( M - N ) / ( MB - N ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
- $ .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+ $ .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
+ IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
NB = 1
- END IF
- IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
- LMINWS = .TRUE.
MB = M
END IF
- IF (LWORK2.LT.NB*N) THEN
+ IF ( LWORK.LT.NB*N ) THEN
LMINWS = .TRUE.
NB = 1
END IF
@@ -222,45 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
- $ .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = NB * N * NBLCKS + 5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = NB * N
- WORK2(2) = N
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = NB*N*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, NB*N )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGEQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
- RETURN
+ ELSE IF ( LQUERY ) THEN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The QR Decomposition
*
- IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
- CALL CGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
- RETURN
+ IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+ CALL CGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
ELSE
- CALL CLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
- $ LWORK2, INFO)
+ CALL CLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+ $ LWORK, INFO )
END IF
+ WORK(1) = MAX( 1, NB*N )
RETURN
*
* End of CGEQR
diff --git a/SRC/cgetsls.f b/SRC/cgetsls.f
index af5bd2cb..c2d5368f 100644
--- a/SRC/cgetsls.f
+++ b/SRC/cgetsls.f
@@ -1,16 +1,15 @@
* Definition:
* ===========
*
-* SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-* $ , WORK, LWORK, INFO )
-
+* SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+* $ WORK, LWORK, INFO )
*
* .. Scalar Arguments ..
* CHARACTER TRANS
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
-* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
+* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
* ..
*
*
@@ -19,10 +18,11 @@
*>
*> \verbatim
*>
-*> CGETSLS solves overdetermined or underdetermined real linear systems
-*> involving an M-by-N matrix A, or its transpose, using a tall skinny
-*> QR or short wide LQfactorization of A. It is assumed that A has
-*> full rank.
+*> CGETSLS solves overdetermined or underdetermined complex linear systems
+*> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
+*> factorization of A. It is assumed that A has full rank.
+*>
+*>
*>
*> The following options are provided:
*>
@@ -80,10 +80,8 @@
*> A is COMPLEX array, dimension (LDA,N)
*> On entry, the M-by-N matrix A.
*> On exit,
-*> if M >= N, A is overwritten by details of its QR
-*> factorization as returned by DGEQRF;
-*> if M < N, A is overwritten by details of its LQ
-*> factorization as returned by DGELQF.
+*> A is overwritten by details of its QR or LQ
+*> factorization as returned by CGEQR or CGEQL.
*> \endverbatim
*>
*> \param[in] LDA
@@ -97,21 +95,17 @@
*> B is COMPLEX array, dimension (LDB,NRHS)
*> On entry, the matrix B of right hand side vectors, stored
*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
-*> if TRANS = 'T'.
+*> if TRANS = 'C'.
*> On exit, if INFO = 0, B is overwritten by the solution
*> vectors, stored columnwise:
*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
-*> squares solution vectors; the residual sum of squares for the
-*> solution in each column is given by the sum of squares of
-*> elements N+1 to M in that column;
+*> squares solution vectors.
*> if TRANS = 'N' and m < n, rows 1 to N of B contain the
*> minimum norm solution vectors;
-*> if TRANS = 'T' and m >= n, rows 1 to M of B contain the
+*> if TRANS = 'C' and m >= n, rows 1 to M of B contain the
*> minimum norm solution vectors;
-*> if TRANS = 'T' and m < n, rows 1 to M of B contain the
-*> least squares solution vectors; the residual sum of squares
-*> for the solution in each column is given by the sum of
-*> squares of elements M+1 to N in that column.
+*> if TRANS = 'C' and m < n, rows 1 to M of B contain the
+*> least squares solution vectors.
*> \endverbatim
*>
*> \param[in] LDB
@@ -122,23 +116,21 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
-*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
-*> LWORK >= max( 1, MN + max( MN, NRHS ) ).
-*> For optimal performance,
-*> LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
-*> where MN = min(M,N) and NB is the optimum block size.
-*>
-*> 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.
+*> If LWORK = -1 or -2, then a workspace query is assumed.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -162,9 +154,11 @@
*
*> \date November 2011
*
+*> \ingroup complexGEsolve
+*
* =====================================================================
- SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
- $ , WORK, LWORK, INFO )
+ SUBROUTINE CGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+ $ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -173,7 +167,7 @@
*
* .. Scalar Arguments ..
CHARACTER TRANS
- INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, MB
+ INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
@@ -190,10 +184,11 @@
* ..
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
- INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
- $ SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2,
- $ INFO2, NB
+ INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+ $ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2,
+ $ INFO2
REAL ANRM, BIGNUM, BNRM, SMLNUM
+ COMPLEX WRK1(5), WRK2
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -206,7 +201,7 @@
$ CTRTRS, XERBLA, CGELQ, CGEMLQ
* ..
* .. Intrinsic Functions ..
- INTRINSIC REAL, MAX, MIN
+ INTRINSIC REAL, MAX, MIN, INT
* ..
* .. Executable Statements ..
*
@@ -218,7 +213,7 @@
MNK = MAX(MINMN,NRHS)
TRAN = LSAME( TRANS, 'C' )
*
- LQUERY = ( LWORK.EQ.-1 )
+ LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
$ LSAME( TRANS, 'C' ) ) ) THEN
INFO = -1
@@ -234,53 +229,58 @@
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
+ IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
- IF ( M.GE.N ) THEN
- CALL CGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
- ELSE
- CALL CGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
- END IF
+ IF ( M.GE.N ) THEN
+ CALL CGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1 , INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL CGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL CGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZM, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
+ ELSE
+ CALL CGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1, INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL CGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL CGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
+ END IF
*
- IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+ IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO=-10
END IF
END IF
*
IF( INFO.NE.0 ) THEN
- CALL XERBLA( 'DGETSLS', -INFO )
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ CALL XERBLA( 'CGETSLS', -INFO )
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
- ELSE IF (LQUERY) THEN
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ END IF
+ IF ( LQUERY ) THEN
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
END IF
- IF(LWORK.LT.WSIZEO) THEN
- LW1=INT(WORK(3))
- LW2=MAX(LW,INT(WORK(6)))
+ IF( LWORK.LT.( TSZ + LW ) ) THEN
+ LW1 = TSZM
+ LW2 = LWM
ELSE
- LW1=INT(WORK(2))
- LW2=MAX(LW,INT(WORK(6)))
+ LW1 = TSZ
+ LW2 = LW
END IF
*
* Quick return if possible
@@ -347,9 +347,9 @@
*
* compute QR factorization of A
*
- CALL CGEQR( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
- IF (.NOT.TRAN) THEN
+ CALL CGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
+ IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
*
@@ -401,8 +401,8 @@
*
* Compute LQ factorization of A
*
- CALL CGELQ( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
+ CALL CGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -482,10 +482,9 @@
END IF
*
50 CONTINUE
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
*
-* End of CGETSLS
+* End of ZGETSLS
*
END
diff --git a/SRC/dgelq.f b/SRC/dgelq.f
index d73f7454..59d9fa91 100644
--- a/SRC/dgelq.f
+++ b/SRC/dgelq.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE DGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-* INFO)
+* SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* DOUBLE PRECISION A( LDA, * ), WORK1( * ), WORK2( * )
+* DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> DGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using DLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise DGELQT:
-*> A = L * Q .
+*> DGELQ computes a LQ factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,8 +42,8 @@
*> On exit, the elements on and below the diagonal of the array
*> contain the M-by-min(M,N) lower trapezoidal matrix L
*> (L is lower triangular if M <= N);
-*> the elements above the diagonal are the rows of
-*> blocked V representing Q (see Further Details).
+*> the elements above the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -56,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> DLASWLQ or DGELQT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): horizontal block size
-*> WORK1(5): vertical block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> DLASWLQ or DGELQT
+*> T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
-*>
+*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -114,27 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLASWLQ or DGELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*> the LQ factorization.
*> \endverbatim
*>
-*>
* =====================================================================
- SUBROUTINE DGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -143,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- DOUBLE PRECISION A( LDA, * ), WORK1( * ), WORK2( * )
+ DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -174,45 +196,53 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
IF ( MIN(M,N).GT.0 ) THEN
- MB = ILAENV( 1, 'DGELQ ', ' ', M, N, 1, -1)
- NB = ILAENV( 1, 'DGELQ ', ' ', M, N, 2, -1)
+ MB = ILAENV( 1, 'DGELQ ', ' ', M, N, 1, -1 )
+ NB = ILAENV( 1, 'DGELQ ', ' ', M, N, 2, -1 )
ELSE
MB = 1
NB = N
END IF
- IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
- IF( NB.GT.N.OR.NB.LE.M) NB = N
- MINLW1 = M + 5
- IF ((NB.GT.M).AND.(N.GT.M)) THEN
- IF(MOD(N-M, NB-M).EQ.0) THEN
- NBLCKS = (N-M)/(NB-M)
+ IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+ IF( NB.GT.N .OR. NB.LE.M ) NB = N
+ MINTSZ = M + 5
+ IF ( NB.GT.M .AND. N.GT.M ) THEN
+ IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+ NBLCKS = ( N - M ) / ( NB - M )
ELSE
- NBLCKS = (N-M)/(NB-M) + 1
+ NBLCKS = ( N - M ) / ( NB - M ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
- $ .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+ $ .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+ $ .AND. ( .NOT.LQUERY) ) THEN
+ IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
MB = 1
- END IF
- IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
- LMINWS = .TRUE.
NB = N
END IF
- IF (LWORK2.LT.MB*M) THEN
+ IF ( LWORK.LT.MB*M ) THEN
LMINWS = .TRUE.
MB = 1
END IF
@@ -224,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
- $ .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS) ) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = MB*M*NBLCKS+5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = MB * M
- WORK2(2) = M
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = MB*M*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, MB*M )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGELQ', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The LQ Decomposition
*
- IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
- CALL DGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+ IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+ CALL DGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
ELSE
- CALL DLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
- $ LWORK2, INFO)
+ CALL DLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+ $ LWORK, INFO)
END IF
+ WORK(1) = MAX( 1, MB*M )
RETURN
*
* End of DGELQ
diff --git a/SRC/dgemlq.f b/SRC/dgemlq.f
index 7bdf97a1..17c4de5c 100644
--- a/SRC/dgemlq.f
+++ b/SRC/dgemlq.f
@@ -2,17 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* DOUBLE A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* DOUBLE PRECISION A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -62,12 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,K)
-*> The i-th row must contain the vector which defines the blocked
-*> elementary reflector H(i), for i = 1,2,...,k, as returned by
-*> DLASWLQ in the first k rows of its array argument A.
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -78,41 +75,42 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1)) is
-*> returned by GEQR.
+*> T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
*> C is DOUBLE PRECISION array, dimension (LDC,N)
*> On entry, the M-by-N matrix C.
*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
*> \param[in] LDC
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -128,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LASWQR or GELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*> the LQ factorization.
+*> This version of GEMLQ will use either LAMSWLQ or GEMLQT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LAMSWLQ or GEMLQT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -157,10 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- DOUBLE PRECISION A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+ DOUBLE PRECISION A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -174,7 +180,7 @@
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
- EXTERNAL DTPMLQT, DGEMLQT, XERBLA
+ EXTERNAL DLAMSWLQ, DGEMLQT, XERBLA
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, MOD
* ..
@@ -182,26 +188,27 @@
*
* Test the input arguments
*
- LQUERY = (LWORK2.LT.0)
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'T' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
- IF (LEFT) THEN
+ MB = INT(T(2))
+ NB = INT(T(3))
+ IF ( LEFT ) THEN
LW = N * MB
MN = M
ELSE
LW = M * MB
MN = N
END IF
- IF ((NB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, NB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(NB-K)
+*
+ IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( NB - K )
ELSE
- NBLCKS = (MN-K)/(NB-K) + 1
+ NBLCKS = ( MN - K ) / ( NB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -220,40 +227,42 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
- INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ INFO = -11
+ ELSE IF(( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
+*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEMLQ', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
- $ (NB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
CALL DGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
- $ WORK1(6), MB, C, LDC, WORK2, INFO)
+ $ T(4), MB, C, LDC, WORK, INFO)
ELSE
- CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ MB, C, LDC, WORK2, LWORK2, INFO )
+ CALL DLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ MB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
*
RETURN
*
diff --git a/SRC/dgemqr.f b/SRC/dgemqr.f
index f47e6a87..ee02d98e 100644
--- a/SRC/dgemqr.f
+++ b/SRC/dgemqr.f
@@ -2,20 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* DOUBLE PRECISION A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
-* ..
-*
-*
+* DOUBLE PRECISION A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -65,13 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,K)
-*> The i-th column must contain the vector which defines the
-*> blockedelementary reflector H(i), for i = 1,2,...,k, as
-*> returned by DGETSQR in the first k columns of
-*> its array argument A.
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
*> \param[in] LDA
@@ -82,16 +75,16 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1)) as
-*> it is returned by GEQR.
+*> T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
@@ -103,21 +96,21 @@
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -133,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLATSQR or DGEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
+*> This version of GEMQR will use either LAMTSQR or GEMQRT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LATMSQR or GEMQRT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -162,11 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- DOUBLE PRECISION A( LDA, * ), WORK1( * ), C(LDC, * ),
- $ WORK2( * )
+ DOUBLE PRECISION A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -180,7 +180,7 @@
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
- EXTERNAL DGEMQRT, DTPMQRT, XERBLA
+ EXTERNAL DGEMQRT, DLAMTSQR, XERBLA
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, MOD
* ..
@@ -188,14 +188,14 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'T' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
+ MB = INT(T(2))
+ NB = INT(T(3))
IF(LEFT) THEN
LW = N * NB
MN = M
@@ -204,11 +204,11 @@
MN = N
END IF
*
- IF ((MB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, MB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(MB-K)
+ IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( MB - K )
ELSE
- NBLCKS = (MN-K)/(MB-K) + 1
+ NBLCKS = ( MN - K ) / ( MB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -227,43 +227,44 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
- ELSE IF( LDC.LT.MAX( 1, M ).AND.MIN(M,N,K).NE.0 ) THEN
+ ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
* Determine the block size if it is tall skinny or short and wide
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEMQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
- $ (MB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
CALL DGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
- $ WORK1(6), NB, C, LDC, WORK2, INFO)
+ $ T(4), NB, C, LDC, WORK, INFO )
ELSE
- CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ NB, C, LDC, WORK2, LWORK2, INFO )
+ CALL DLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ NB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
*
RETURN
*
diff --git a/SRC/dgeqr.f b/SRC/dgeqr.f
index da0fc4ad..93f02bd8 100644
--- a/SRC/dgeqr.f
+++ b/SRC/dgeqr.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE DGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+* SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* DOUBLE PRECISION A( LDA, * ), WORK1( * ), WORK2( * )
+* DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> DGEQR computes a QR factorization of an M-by-N matrix A,
-*> using DLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise DGEQRT:
-*> A = Q * R .
+*> DGEQR computes a QR factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,7 +42,8 @@
*> On exit, the elements on and above the diagonal of the array
*> contain the min(M,N)-by-N upper trapezoidal matrix R
*> (R is upper triangular if M >= N);
-*> the elements below the diagonal represent Q (see Further Details).
+*> the elements below the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -55,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is DOUBLE PRECISION array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> DLATSQR or DGEQRT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): row block size
-*> WORK1(5): column block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> DLATSQR or DGEQRT
+*> T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK2))
+*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -113,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLATSQR or DGEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE DGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- DOUBLE PRECISION A( LDA, * ), WORK1( * ), WORK2( * )
+ DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -172,45 +196,53 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
IF ( MIN(M,N).GT.0 ) THEN
- MB = ILAENV( 1, 'DGEQR ', ' ', M, N, 1, -1)
- NB = ILAENV( 1, 'DGEQR ', ' ', M, N, 2, -1)
+ MB = ILAENV( 1, 'DGEQR ', ' ', M, N, 1, -1 )
+ NB = ILAENV( 1, 'DGEQR ', ' ', M, N, 2, -1 )
ELSE
MB = M
NB = 1
END IF
- IF( MB.GT.M.OR.MB.LE.N) MB = M
- IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
- MINLW1 = N + 5
- IF ((MB.GT.N).AND.(M.GT.N)) THEN
- IF(MOD(M-N, MB-N).EQ.0) THEN
- NBLCKS = (M-N)/(MB-N)
+ IF( MB.GT.M .OR. MB.LE.N ) MB = M
+ IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+ MINTSZ = N + 5
+ IF ( MB.GT.N .AND. M.GT.N ) THEN
+ IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+ NBLCKS = ( M - N ) / ( MB - N )
ELSE
- NBLCKS = (M-N)/(MB-N) + 1
+ NBLCKS = ( M - N ) / ( MB - N ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
- $ .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+ $ .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
+ IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
NB = 1
- END IF
- IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
- LMINWS = .TRUE.
MB = M
END IF
- IF (LWORK2.LT.NB*N) THEN
+ IF ( LWORK.LT.NB*N ) THEN
LMINWS = .TRUE.
NB = 1
END IF
@@ -222,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
- $ .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
-
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = NB * N * NBLCKS + 5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = NB * N
- WORK2(2) = N
+*
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = NB*N*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, NB*N )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
- RETURN
+ ELSE IF ( LQUERY ) THEN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The QR Decomposition
*
- IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
- CALL DGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+ IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+ CALL DGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
ELSE
- CALL DLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
- $ LWORK2, INFO)
+ CALL DLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+ $ LWORK, INFO )
END IF
+ WORK(1) = MAX( 1, NB*N )
RETURN
*
* End of DGEQR
diff --git a/SRC/dgetsls.f b/SRC/dgetsls.f
index f2a7b3a7..0a94694e 100644
--- a/SRC/dgetsls.f
+++ b/SRC/dgetsls.f
@@ -29,16 +29,21 @@
*> 1. If TRANS = 'N' and m >= n: find the least squares solution of
*> an overdetermined system, i.e., solve the least squares problem
*> minimize || B - A*X ||.
-
+*>
*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
*> an underdetermined system A * X = B.
-
+*>
*> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
*> an undetermined system A**T * X = B.
-
+*>
*> 4. If TRANS = 'T' and m < n: find the least squares solution of
*> an overdetermined system, i.e., solve the least squares problem
*> minimize || B - A**T * X ||.
+*>
+*> Several right hand side vectors b and solution vectors x can be
+*> handled in a single call; they are stored as the columns of the
+*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
+*> matrix X.
*> \endverbatim
*
* Arguments:
@@ -76,7 +81,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
-*> factorization as returned by DGETSQR.
+*> factorization as returned by DGEQR or DGEQL.
*> \endverbatim
*>
*> \param[in] LDA
@@ -111,18 +116,21 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is DOUBLE PRECISION array, dimension (MAX(12,LWORK))
-*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK,
-*> and WORK(2) returns the minimum LWORK.
+*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
-*> IF LWORK=-1, workspace query is assumed, and
-*> WORK(1) returns the optimal LWORK,
-*> and WORK(2) returns the minimum LWORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -174,10 +182,10 @@
* ..
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
- INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
- $ SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2,
+ INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+ $ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
$ INFO2
- DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
+ DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, WRK1(5), WRK2
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -190,7 +198,7 @@
$ DTRTRS, XERBLA, DGELQ, DGEMLQ
* ..
* .. Intrinsic Functions ..
- INTRINSIC DBLE, MAX, MIN
+ INTRINSIC DBLE, MAX, MIN, INT
* ..
* .. Executable Statements ..
*
@@ -199,10 +207,10 @@
INFO=0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
- MNK = MAX(MINMN,NRHS)
+ MNK = MAX( MINMN, NRHS )
TRAN = LSAME( TRANS, 'T' )
*
- LQUERY = ( LWORK.EQ.-1 )
+ LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
$ LSAME( TRANS, 'T' ) ) ) THEN
INFO = -1
@@ -218,49 +226,58 @@
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
+ IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
- CALL DGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- LW = INT(WORK(6))
- CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+ CALL DGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1 , INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL DGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL DGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZM, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
ELSE
- CALL DGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- LW = INT(WORK(6))
- CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+ CALL DGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1, INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL DGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL DGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
END IF
*
- IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+ IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
END IF
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGETSLS', -INFO )
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
- ELSE IF (LQUERY) THEN
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ END IF
+ IF ( LQUERY ) THEN
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
END IF
- IF(LWORK.LT.WSIZEO) THEN
- LW1=INT(WORK(3))
- LW2=MAX(LW,INT(WORK(6)))
+ IF( LWORK.LT.( TSZ + LW ) ) THEN
+ LW1 = TSZM
+ LW2 = LWM
ELSE
- LW1=INT(WORK(2))
- LW2=MAX(LW,INT(WORK(6)))
+ LW1 = TSZ
+ LW2 = LW
END IF
*
* Quick return if possible
@@ -323,13 +340,13 @@
IBSCL = 2
END IF
*
- IF ( M.GE.N) THEN
+ IF ( M.GE.N ) THEN
*
* compute QR factorization of A
*
CALL DGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
- IF (.NOT.TRAN) THEN
+ IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
*
@@ -381,7 +398,7 @@
*
* Compute LQ factorization of A
*
- CALL DGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+ CALL DGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
$ WORK(1), LW2, INFO )
*
* workspace at least M, optimally M*NB.
@@ -454,16 +471,15 @@
$ INFO )
END IF
IF( IBSCL.EQ.1 ) THEN
- CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
+ CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
ELSE IF( IBSCL.EQ.2 ) THEN
- CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
+ CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
$ INFO )
END IF
*
50 CONTINUE
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
*
* End of DGETSLS
diff --git a/SRC/sgelq.f b/SRC/sgelq.f
index 8a759834..adc606d9 100644
--- a/SRC/sgelq.f
+++ b/SRC/sgelq.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE SGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-* INFO)
+* SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* REAL A( LDA, * ), WORK1( * ), WORK2( * )
+* REAL A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> SGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using SLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise SGELQT:
-*> A = L * Q .
+*> SGELQ computes a LQ factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,8 +42,8 @@
*> On exit, the elements on and below the diagonal of the array
*> contain the M-by-min(M,N) lower trapezoidal matrix L
*> (L is lower triangular if M <= N);
-*> the elements above the diagonal are the rows of
-*> blocked V representing Q (see Further Details).
+*> the elements above the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -56,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is REAL array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> SLASWLQ or SGELQT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): horizontal block size
-*> WORK1(5): vertical block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> SLASWLQ or SGELQT
+*> T is REAL array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) REAL array, dimension (MAX(1,LWORK2))
-*>
+*> (workspace) REAL array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -114,27 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
+*>
+*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLASWLQ or DGELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*> the LQ factorization.
*> \endverbatim
*>
-*>
* =====================================================================
- SUBROUTINE SGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -143,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- REAL A( LDA, * ), WORK1( * ), WORK2( * )
+ REAL A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -174,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -185,34 +218,31 @@
MB = 1
NB = N
END IF
- IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
- IF( NB.GT.N.OR.NB.LE.M) NB = N
- MINLW1 = M + 5
- IF ((NB.GT.M).AND.(N.GT.M)) THEN
- IF(MOD(N-M, NB-M).EQ.0) THEN
- NBLCKS = (N-M)/(NB-M)
+ IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+ IF( NB.GT.N .OR. NB.LE.M ) NB = N
+ MINTSZ = M + 5
+ IF ( NB.GT.M .AND. N.GT.M ) THEN
+ IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+ NBLCKS = ( N - M ) / ( NB - M )
ELSE
- NBLCKS = (N-M)/(NB-M) + 1
+ NBLCKS = ( N - M ) / ( NB - M ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
- $ .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+ $ .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+ $ .AND. ( .NOT.LQUERY) ) THEN
+ IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
MB = 1
- END IF
- IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
- LMINWS = .TRUE.
NB = N
END IF
- IF (LWORK2.LT.MB*M) THEN
+ IF ( LWORK.LT.MB*M ) THEN
LMINWS = .TRUE.
MB = 1
END IF
@@ -224,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
- $ .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS) ) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = MB*M*NBLCKS+5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = MB * M
- WORK2(2) = M
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = MB*M*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, MB*M )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGELQ', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The LQ Decomposition
*
- IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
- CALL SGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+ IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+ CALL SGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
ELSE
- CALL SLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
- $ LWORK2, INFO)
+ CALL SLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+ $ LWORK, INFO)
END IF
+ WORK(1) = MAX( 1, MB*M )
RETURN
*
* End of SGELQ
diff --git a/SRC/sgemlq.f b/SRC/sgemlq.f
index 14a37a4d..a9cd54bd 100644
--- a/SRC/sgemlq.f
+++ b/SRC/sgemlq.f
@@ -2,23 +2,22 @@
* Definition:
* ===========
*
-* SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* REAL A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* REAL A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
*> \verbatim
*>
-*> DGEMLQ overwrites the general real M-by-N matrix C with
+*> SGEMLQ overwrites the general real M-by-N matrix C with
*>
*>
*> SIDE = 'L' SIDE = 'R'
@@ -26,7 +25,7 @@
*> TRANS = 'T': Q**T * C C * Q**T
*> where Q is a real orthogonal matrix defined as the product
*> of blocked elementary reflectors computed by short wide LQ
-*> factorization (DGELQ)
+*> factorization (SGELQ)
*> \endverbatim
*
* Arguments:
@@ -62,12 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is REAL array, dimension (LDA,K)
-*> The i-th row must contain the vector which defines the blocked
-*> elementary reflector H(i), for i = 1,2,...,k, as returned by
-*> DLASWLQ in the first k rows of its array argument A.
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -78,41 +75,42 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is REAL array, dimension (MAX(1,LWORK1)) is
-*> returned by GEQR.
+*> T is REAL array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
*> C is REAL array, dimension (LDC,N)
*> On entry, the M-by-N matrix C.
*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
*> \param[in] LDC
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) REAL array, dimension (MAX(1,LWORK2))
+*> (workspace) REAL array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -128,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LASWLQ or GELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*> the LQ factorization.
+*> This version of GEMLQ will use either LAMSWLQ or GEMLQT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LAMSWLQ or GEMLQT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -157,10 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- REAL A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+ REAL A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -174,7 +180,7 @@
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
- EXTERNAL STPMLQT, SGEMLQT, XERBLA
+ EXTERNAL SLAMSWLQ, SGEMLQT, XERBLA
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, MOD
* ..
@@ -182,26 +188,27 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'T' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
- IF (LEFT) THEN
+ MB = INT(T(2))
+ NB = INT(T(3))
+ IF ( LEFT ) THEN
LW = N * MB
MN = M
ELSE
LW = M * MB
MN = N
END IF
- IF ((NB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, NB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(NB-K)
+*
+ IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( NB - K )
ELSE
- NBLCKS = (MN-K)/(NB-K) + 1
+ NBLCKS = ( MN - K ) / ( NB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -220,40 +227,43 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
- INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ INFO = -11
+ ELSE IF(( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = REAL(LW)
END IF
+*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGEMLQ', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
- $ (NB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
CALL SGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
- $ WORK1(6), MB, C, LDC, WORK2, INFO)
+ $ T(4), MB, C, LDC, WORK, INFO)
ELSE
- CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ MB, C, LDC, WORK2, LWORK2, INFO )
+ CALL SLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ MB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = REAL(LW)
+*
RETURN
*
* End of SGEMLQ
diff --git a/SRC/sgemqr.f b/SRC/sgemqr.f
index cda7990c..78eb75c8 100644
--- a/SRC/sgemqr.f
+++ b/SRC/sgemqr.f
@@ -2,17 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* REAL A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* REAL A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -26,7 +25,7 @@
*> TRANS = 'T': Q**T * C C * Q**T
*> where Q is a real orthogonal matrix defined as the product
*> of blocked elementary reflectors computed by tall skinny
-*> QR factorization (DGEQR)
+*> QR factorization (SGEQR)
*> \endverbatim
*
* Arguments:
@@ -62,13 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is REAL array, dimension (LDA,K)
-*> The i-th column must contain the vector which defines the
-*> blockedelementary reflector H(i), for i = 1,2,...,k, as
-*> returned by DGETSQR in the first k columns of
-*> its array argument A.
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
*> \param[in] LDA
@@ -79,16 +75,16 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is REAL array, dimension (MAX(1,LWORK1)) as
-*> it is returned by GEQR.
+*> T is REAL array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
@@ -100,21 +96,21 @@
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) REAL array, dimension (MAX(1,LWORK2))
+*> (workspace) REAL array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -130,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LATSQR or GEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
+*> This version of GEMQR will use either LAMTSQR or GEMQRT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LATMSQR or GEMQRT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -159,11 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- REAL A( LDA, * ), WORK1( * ), C(LDC, * ),
- $ WORK2( * )
+ REAL A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -177,7 +180,7 @@
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
- EXTERNAL SGEMQRT, STPMQRT, XERBLA
+ EXTERNAL SGEMQRT, SLAMTSQR, XERBLA
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN, MOD
* ..
@@ -185,14 +188,14 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'T' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
+ MB = INT(T(2))
+ NB = INT(T(3))
IF(LEFT) THEN
LW = N * NB
MN = M
@@ -201,11 +204,11 @@
MN = N
END IF
*
- IF ((MB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, MB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(MB-K)
+ IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( MB - K )
ELSE
- NBLCKS = (MN-K)/(MB-K) + 1
+ NBLCKS = ( MN - K ) / ( MB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -224,43 +227,44 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
- ELSE IF( LDC.LT.MAX( 1, M ).AND.MIN(M,N,K).NE.0 ) THEN
+ ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
* Determine the block size if it is tall skinny or short and wide
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGEMQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
- $ (MB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
CALL SGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
- $ WORK1(6), NB, C, LDC, WORK2, INFO)
+ $ T(4), NB, C, LDC, WORK, INFO )
ELSE
- CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ NB, C, LDC, WORK2, LWORK2, INFO )
+ CALL SLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ NB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
*
RETURN
*
diff --git a/SRC/sgeqr.f b/SRC/sgeqr.f
index 41e04622..ad18e8a2 100644
--- a/SRC/sgeqr.f
+++ b/SRC/sgeqr.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE SGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+* SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* REAL A( LDA, * ), WORK1( * ), WORK2( * )
+* REAL A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> SGEQR computes a QR factorization of an M-by-N matrix A,
-*> using SLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise SGEQRT:
-*> A = Q * R .
+*> SGEQR computes a QR factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,7 +42,8 @@
*> On exit, the elements on and above the diagonal of the array
*> contain the min(M,N)-by-N upper trapezoidal matrix R
*> (R is upper triangular if M >= N);
-*> the elements below the diagonal represent Q (see Further Details).
+*> the elements below the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -55,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is REAL array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> DLATSQR or DGEQRT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): row block size
-*> WORK1(5): column block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> SLATSQR or SGEQRT
+*> T is REAL array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) REAL array, dimension (MAX(1,LWORK2))
+*> (workspace) REAL array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -113,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLATSQR or DGEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE SGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- REAL A( LDA, * ), WORK1( * ), WORK2( * )
+ REAL A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -172,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -183,34 +218,31 @@
MB = M
NB = 1
END IF
- IF( MB.GT.M.OR.MB.LE.N) MB = M
- IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
- MINLW1 = N + 5
- IF ((MB.GT.N).AND.(M.GT.N)) THEN
- IF(MOD(M-N, MB-N).EQ.0) THEN
- NBLCKS = (M-N)/(MB-N)
+ IF( MB.GT.M .OR. MB.LE.N ) MB = M
+ IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+ MINTSZ = N + 5
+ IF ( MB.GT.N .AND. M.GT.N ) THEN
+ IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+ NBLCKS = ( M - N ) / ( MB - N )
ELSE
- NBLCKS = (M-N)/(MB-N) + 1
+ NBLCKS = ( M - N ) / ( MB - N ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
- $ .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+ $ .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
+ IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
NB = 1
- END IF
- IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
- LMINWS = .TRUE.
MB = M
END IF
- IF (LWORK2.LT.NB*N) THEN
+ IF ( LWORK.LT.NB*N ) THEN
LMINWS = .TRUE.
NB = 1
END IF
@@ -222,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
- $ .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
-
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = NB * N * NBLCKS + 5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = NB * N
- WORK2(2) = N
+*
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = NB*N*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, NB*N )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGEQR', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The QR Decomposition
*
- IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
- CALL SGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+ IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+ CALL SGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
ELSE
- CALL SLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
- $ LWORK2, INFO)
+ CALL SLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+ $ LWORK, INFO )
END IF
+ WORK(1) = MAX( 1, NB*N )
RETURN
*
* End of SGEQR
diff --git a/SRC/sgetsls.f b/SRC/sgetsls.f
index b7bcd0f0..dcbaf36a 100644
--- a/SRC/sgetsls.f
+++ b/SRC/sgetsls.f
@@ -1,16 +1,15 @@
* Definition:
* ===========
*
-* SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-* $ , WORK, LWORK, INFO )
-
+* SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+* $ WORK, LWORK, INFO )
*
* .. Scalar Arguments ..
* CHARACTER TRANS
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
-* REAL A( LDA, * ), B( LDB, * ), WORK( * )
+* REAL A( LDA, * ), B( LDB, * ), WORK( * )
* ..
*
*
@@ -77,7 +76,7 @@
*> On entry, the M-by-N matrix A.
*> On exit,
*> A is overwritten by details of its QR or LQ
-*> factorization as returned by DGETSQR.
+*> factorization as returned by SGEQR or SGEQL.
*> \endverbatim
*>
*> \param[in] LDA
@@ -112,18 +111,21 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is REAL array, dimension (MAX(1,LWORK))
-*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK,
-*> and WORK(2) returns the minimum LWORK.
+*> (workspace) REAL array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
-*> IF LWORK=-1, workspace query is assumed, and
-*> WORK(1) returns the optimal LWORK,
-*> and WORK(2) returns the minimum LWORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -150,8 +152,8 @@
*> \ingroup doubleGEsolve
*
* =====================================================================
- SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
- $ , WORK, LWORK, INFO )
+ SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+ $ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -160,7 +162,7 @@
*
* .. Scalar Arguments ..
CHARACTER TRANS
- INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, MB
+ INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
REAL A( LDA, * ), B( LDB, * ), WORK( * )
@@ -175,10 +177,11 @@
* ..
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
- INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
- $ SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2, INFO2,
- $ NB
+ INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+ $ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
+ $ INFO2
REAL ANRM, BIGNUM, BNRM, SMLNUM
+ COMPLEX WRK1(5), WRK2
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -191,7 +194,7 @@
$ STRTRS, XERBLA, SGELQ, SGEMLQ
* ..
* .. Intrinsic Functions ..
- INTRINSIC REAL, MAX, MIN
+ INTRINSIC REAL, MAX, MIN, INT
* ..
* .. Executable Statements ..
*
@@ -200,10 +203,10 @@
INFO=0
MINMN = MIN( M, N )
MAXMN = MAX( M, N )
- MNK = MAX(MINMN,NRHS)
+ MNK = MAX( MINMN, NRHS )
TRAN = LSAME( TRANS, 'T' )
*
- LQUERY = ( LWORK.EQ.-1 )
+ LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
$ LSAME( TRANS, 'T' ) ) ) THEN
INFO = -1
@@ -219,53 +222,58 @@
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
+ IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
IF ( M.GE.N ) THEN
- CALL SGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+ CALL SGEQR( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1 , INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL SGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL SGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZM, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
ELSE
- CALL SGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
+ CALL SGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1, INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL SGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL SGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
END IF
*
- IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+ IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
END IF
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGETSLS', -INFO )
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
- ELSE IF (LQUERY) THEN
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ END IF
+ IF ( LQUERY ) THEN
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
END IF
- IF(LWORK.LT.WSIZEO) THEN
- LW1=INT(WORK(3))
- LW2=MAX(LW,INT(WORK(6)))
+ IF( LWORK.LT.( TSZ + LW ) ) THEN
+ LW1 = TSZM
+ LW2 = LWM
ELSE
- LW1=INT(WORK(2))
- LW2=MAX(LW,INT(WORK(6)))
+ LW1 = TSZ
+ LW2 = LW
END IF
*
* Quick return if possible
@@ -332,9 +340,9 @@
*
* compute QR factorization of A
*
- CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
- IF (.NOT.TRAN) THEN
+ CALL SGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
+ IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
*
@@ -386,8 +394,8 @@
*
* Compute LQ factorization of A
*
- CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
+ CALL SGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -467,8 +475,7 @@
END IF
*
50 CONTINUE
- WORK( 1 ) = REAL( WSIZEO )
- WORK( 2 ) = REAL( WSIZEM )
+ WORK( 1 ) = REAL( TSZ + LW )
RETURN
*
* End of SGETSLS
diff --git a/SRC/zgelq.f b/SRC/zgelq.f
index 33125b3d..5c51cf52 100644
--- a/SRC/zgelq.f
+++ b/SRC/zgelq.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE CGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
-* INFO)
+* SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
+* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), WORK1( * ), WORK2( * )
+* COMPLEX*16 A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> ZGELQ computes an LQ factorization of an M-by-N matrix A,
-*> using ZLASWLQ when A is short and wide
-*> (N sufficiently greater than M), and otherwise ZGELQT:
-*> A = L * Q .
+*> ZGELQ computes a LQ factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,8 +42,8 @@
*> On exit, the elements on and below the diagonal of the array
*> contain the M-by-min(M,N) lower trapezoidal matrix L
*> (L is lower triangular if M <= N);
-*> the elements above the diagonal are the rows of
-*> blocked V representing Q (see Further Details).
+*> the elements above the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -56,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> ZLASWLQ or ZGELQT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): horizontal block size
-*> WORK1(5): vertical block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> ZLASWLQ or ZGELQT
+*> T is COMPLEX*16 array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
-*>
+*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
-*> \param[in] LWORK2
+*>
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -114,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any LQ factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The upper part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLASWLQ or DGELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is short-and-wide) or GELQT to compute
+*> the LQ factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE ZGELQ( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE ZGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -142,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), WORK1( * ), WORK2( * )
+ COMPLEX*16 A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -173,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -184,34 +218,31 @@
MB = 1
NB = N
END IF
- IF( MB.GT.MIN(M,N).OR.MB.LT.1) MB = 1
- IF( NB.GT.N.OR.NB.LE.M) NB = N
- MINLW1 = M + 5
- IF ((NB.GT.M).AND.(N.GT.M)) THEN
- IF(MOD(N-M, NB-M).EQ.0) THEN
- NBLCKS = (N-M)/(NB-M)
+ IF( MB.GT.MIN( M, N ) .OR. MB.LT.1 ) MB = 1
+ IF( NB.GT.N .OR. NB.LE.M ) NB = N
+ MINTSZ = M + 5
+ IF ( NB.GT.M .AND. N.GT.M ) THEN
+ IF( MOD( N - M, NB - M ).EQ.0 ) THEN
+ NBLCKS = ( N - M ) / ( NB - M )
ELSE
- NBLCKS = (N-M)/(NB-M) + 1
+ NBLCKS = ( N - M ) / ( NB - M ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1,MB*M*NBLCKS+5)
- $ .OR.(LWORK2.LT.MB*M)).AND.(LWORK2.GE.M).AND.(LWORK1.GE.M+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1,MB*M*NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) .OR. LWORK.LT.MB*M )
+ $ .AND. ( LWORK.GE.M ) .AND. ( TSIZE.GE.M + 5 )
+ $ .AND. ( .NOT.LQUERY) ) THEN
+ IF ( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
MB = 1
- END IF
- IF (LWORK1.LT.MAX(1,M*NBLCKS+5)) THEN
- LMINWS = .TRUE.
NB = N
END IF
- IF (LWORK2.LT.MB*M) THEN
+ IF ( LWORK.LT.MB*M ) THEN
LMINWS = .TRUE.
MB = 1
END IF
@@ -223,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, MB*M*NBLCKS+5 )
- $ .AND.(.NOT.LQUERY).AND. (.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*M*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,M*MB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS) ) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, M*MB ) ) .AND .( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = MB*M*NBLCKS+5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = MB * M
- WORK2(2) = M
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = MB*M*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, MB*M )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGELQ', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The LQ Decomposition
*
- IF((N.LE.M).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
- CALL ZGELQT( M, N, MB, A, LDA, WORK1(6), MB, WORK2, INFO)
+ IF( ( N.LE.M ) .OR. ( NB.LE.M ) .OR. ( NB.GE.N ) ) THEN
+ CALL ZGELQT( M, N, MB, A, LDA, T(4), MB, WORK, INFO)
ELSE
- CALL ZLASWLQ( M, N, MB, NB, A, LDA, WORK1(6), MB, WORK2,
- $ LWORK2, INFO)
+ CALL ZLASWLQ( M, N, MB, NB, A, LDA, T(4), MB, WORK,
+ $ LWORK, INFO)
END IF
+ WORK(1) = MAX( 1, MB*M )
RETURN
*
* End of ZGELQ
diff --git a/SRC/zgemlq.f b/SRC/zgemlq.f
index 10d3a5e4..f02d7b1a 100644
--- a/SRC/zgemlq.f
+++ b/SRC/zgemlq.f
@@ -2,31 +2,30 @@
* Definition:
* ===========
*
-* SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, MB, NB, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* COMPLEX*16 A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
*> \verbatim
*>
-*> ZGEMLQ overwrites the general real M-by-N matrix C with
+*> ZGEMLQ overwrites the general real M-by-N matrix C with
*>
*>
-*> SIDE = 'L' SIDE = 'R'
-*> TRANS = 'N': Q * C C * Q
-*> TRANS = 'T': Q**T * C C * Q**T
-*> where Q is a complex orthogonal matrix defined as the product
-*> of blocked elementary reflectors computed by short wide LQ
-*> factorization (DGELQ)
+*> SIDE = 'L' SIDE = 'R'
+*> TRANS = 'N': Q * C C * Q
+*> TRANS = 'C': Q**H * C C * Q**H
+*> where Q is a complex unitary matrix defined as the product
+*> of blocked elementary reflectors computed by short wide
+*> LQ factorization (ZGELQ)
*> \endverbatim
*
* Arguments:
@@ -62,12 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,K)
-*> The i-th row must contain the vector which defines the blocked
-*> elementary reflector H(i), for i = 1,2,...,k, as returned by
-*> DLASWLQ in the first k rows of its array argument A.
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
*> \param[in] LDA
@@ -78,41 +75,42 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1)) is
-*> returned by GEQR.
+*> T is COMPLEX*16 array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by ZGELQ.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
*> C is COMPLEX*16 array, dimension (LDC,N)
*> On entry, the M-by-N matrix C.
*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+*>
*> \param[in] LDC
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -128,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LASWLQ or GELQT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GELQ will use either
-*> LASWLQ(if the matrix is short-and-wide) or GELQT to compute
-*> the LQ decomposition.
-*> The output of LASWLQ or GELQT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB, NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LASWLQ or GELQT was used is the same as used below in
-*> GELQ. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LASWLQ or GELQT.
+*> LASWLQ (if the matrix is wide-and-short) or GELQT to compute
+*> the LQ factorization.
+*> This version of GEMLQ will use either LAMSWLQ or GEMLQT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LAMSWLQ or GEMLQT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE ZGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -157,10 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), C( LDC, * ), WORK1( * ), WORK2( * )
+ COMPLEX*16 A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -182,26 +188,27 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'C' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
- IF (LEFT) THEN
+ MB = INT(T(2))
+ NB = INT(T(3))
+ IF ( LEFT ) THEN
LW = N * MB
MN = M
ELSE
LW = M * MB
MN = N
END IF
- IF ((NB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, NB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(NB-K)
+*
+ IF ( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( NB - K )
ELSE
- NBLCKS = (MN-K)/(NB-K) + 1
+ NBLCKS = ( MN - K ) / ( NB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -220,40 +227,43 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, MB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, MB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
- INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ INFO = -11
+ ELSE IF(( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
+*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGEMLQ', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(NB.LE.K).OR.
- $ (NB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( NB.LE.K ) .OR. ( NB.GE.MAX( M, N, K ) ) ) THEN
CALL ZGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
- $ WORK1(6), MB, C, LDC, WORK2, INFO)
+ $ T(4), MB, C, LDC, WORK, INFO)
ELSE
- CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ MB, C, LDC, WORK2, LWORK2, INFO )
+ CALL ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ MB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
+*
RETURN
*
* End of ZGEMLQ
diff --git a/SRC/zgemqr.f b/SRC/zgemqr.f
index 3141067f..71614d5e 100644
--- a/SRC/zgemqr.f
+++ b/SRC/zgemqr.f
@@ -2,17 +2,16 @@
* Definition:
* ===========
*
-* SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1,
-* $ LWORK1, C, LDC, WORK2, LWORK2, INFO )
+* SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, T,
+* $ TSIZE, C, LDC, WORK, LWORK, INFO )
*
*
* .. Scalar Arguments ..
* CHARACTER SIDE, TRANS
-* INTEGER INFO, LDA, M, N, K, LDT, LWORK1, LWORK2, LDC
+* INTEGER INFO, LDA, M, N, K, LDT, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), WORK1( * ), C(LDC, * ),
-* $ WORK2( * )
+* COMPLEX*16 A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
*> \par Purpose:
* =============
*>
@@ -23,8 +22,8 @@
*>
*> SIDE = 'L' SIDE = 'R'
*> TRANS = 'N': Q * C C * Q
-*> TRANS = 'T': Q**T * C C * Q**T
-*> where Q is a complex orthogonal matrix defined as the product
+*> TRANS = 'T': Q**H * C C * Q**H
+*> where Q is a complex unitary matrix defined as the product
*> of blocked elementary reflectors computed by tall skinny
*> QR factorization (ZGEQR)
*> \endverbatim
@@ -62,13 +61,10 @@
*>
*> \endverbatim
*>
-*> \param[in,out] A
+*> \param[in] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,K)
-*> The i-th column must contain the vector which defines the
-*> blockedelementary reflector H(i), for i = 1,2,...,k, as
-*> returned by DGETSQR in the first k columns of
-*> its array argument A.
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
*> \param[in] LDA
@@ -79,16 +75,16 @@
*> if SIDE = 'R', LDA >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK1
+*> \param[in] T
*> \verbatim
-*> WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1)) as
-*> it is returned by GEQR.
+*> T is COMPLEX*16 array, dimension (MAX(5,TSIZE)).
+*> Part of the data structure to represent Q as returned by SGEQR.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
+*> TSIZE is INTEGER
+*> The dimension of the array T. TSIZE >= 5.
*> \endverbatim
*>
*> \param[in,out] C
@@ -100,21 +96,21 @@
*> LDC is INTEGER
*> The leading dimension of the array C. LDC >= max(1,M).
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
*>
*> \endverbatim
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a workspace query is assumed; the routine
-*> only calculates the optimal size of the WORK2 array, returns
-*> this value as the third entry of the WORK2 array (WORK2(1)),
-*> and no error message related to LWORK2 is issued by XERBLA.
-*>
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1, then a workspace query is assumed. The routine
+*> only calculates the size of the WORK array, returns this
+*> value as WORK(1), and no error message related to WORK
+*> is issued by XERBLA.
*> \endverbatim
+*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@@ -130,27 +126,35 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> LATSQR or GEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
+*> This version of GEMQR will use either LAMTSQR or GEMQRT to
+*> multiply matrix Q by another matrix.
+*> Further Details in LATMSQR or GEMQRT.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, WORK1, LWORK1,
- $ C, LDC, WORK2, LWORK2, INFO )
+ SUBROUTINE ZGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
+ $ C, LDC, WORK, LWORK, INFO )
*
* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -159,11 +163,10 @@
*
* .. Scalar Arguments ..
CHARACTER SIDE, TRANS
- INTEGER INFO, LDA, M, N, K, LWORK1, LWORK2, LDC
+ INTEGER INFO, LDA, M, N, K, TSIZE, LWORK, LDC
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), WORK1( * ), C(LDC, * ),
- $ WORK2( * )
+ COMPLEX*16 A( LDA, * ), T( * ), C(LDC, * ), WORK( * )
* ..
*
* =====================================================================
@@ -185,14 +188,14 @@
*
* Test the input arguments
*
- LQUERY = LWORK2.LT.0
+ LQUERY = LWORK.LT.0
NOTRAN = LSAME( TRANS, 'N' )
TRAN = LSAME( TRANS, 'C' )
LEFT = LSAME( SIDE, 'L' )
RIGHT = LSAME( SIDE, 'R' )
*
- MB = INT(WORK1(4))
- NB = INT(WORK1(5))
+ MB = INT(T(2))
+ NB = INT(T(3))
IF(LEFT) THEN
LW = N * NB
MN = M
@@ -201,11 +204,11 @@
MN = N
END IF
*
- IF ((MB.GT.K).AND.(MN.GT.K)) THEN
- IF(MOD(MN-K, MB-K).EQ.0) THEN
- NBLCKS = (MN-K)/(MB-K)
+ IF ( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
+ IF( MOD( MN - K, MB - K ) .EQ. 0 ) THEN
+ NBLCKS = ( MN - K ) / ( MB - K )
ELSE
- NBLCKS = (MN-K)/(MB-K) + 1
+ NBLCKS = ( MN - K ) / ( MB - K ) + 1
END IF
ELSE
NBLCKS = 1
@@ -224,45 +227,47 @@
INFO = -5
ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
INFO = -7
- ELSE IF( LWORK1.LT.MAX( 1, NB*K*NBLCKS+5 )) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*K*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
INFO = -9
- ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
+ ELSE IF( LDC.LT.MAX( 1, M ) .AND. MIN( M, N, K ).NE.0 ) THEN
INFO = -11
- ELSE IF(( LWORK2.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
INFO = -13
END IF
*
* Determine the block size if it is tall skinny or short and wide
*
- IF( INFO.EQ.0) THEN
- WORK2(1) = LW
+ IF( INFO.EQ.0 ) THEN
+ WORK(1) = LW
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGEMQR', -INFO )
RETURN
- ELSE IF (LQUERY) THEN
+ ELSE IF ( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
- IF( MIN(M,N,K).EQ.0 ) THEN
+ IF( MIN( M, N, K ).EQ.0 ) THEN
RETURN
END IF
*
- IF((LEFT.AND.M.LE.K).OR.(RIGHT.AND.N.LE.K).OR.(MB.LE.K).OR.
- $ (MB.GE.MAX(M,N,K))) THEN
+ IF( ( LEFT .AND. M.LE.K ) .OR. ( RIGHT .AND. N.LE.K )
+ $ .OR. ( MB.LE.K ) .OR. ( MB.GE.MAX( M, N, K ) ) ) THEN
CALL ZGEMQRT( SIDE, TRANS, M, N, K, NB, A, LDA,
- $ WORK1(6), NB, C, LDC, WORK2, INFO)
+ $ T(4), NB, C, LDC, WORK, INFO )
ELSE
- CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, WORK1(6),
- $ NB, C, LDC, WORK2, LWORK2, INFO )
+ CALL ZLAMTSQR( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T(4),
+ $ NB, C, LDC, WORK, LWORK, INFO )
END IF
*
- WORK2(1) = LW
+ WORK(1) = LW
+*
RETURN
*
-* End of DGEMQR
+* End of ZGEMQR
*
END
diff --git a/SRC/zgeqr.f b/SRC/zgeqr.f
index 10fab97f..226db33e 100644
--- a/SRC/zgeqr.f
+++ b/SRC/zgeqr.f
@@ -2,14 +2,14 @@
* Definition:
* ===========
*
-* SUBROUTINE ZGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+* SUBROUTINE ZGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
* INFO)
*
* .. Scalar Arguments ..
-* INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+* INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), WORK1( * ), WORK2( * )
+* COMPLEX*16 A( LDA, * ), T( * ), WORK( * )
* ..
*
*
@@ -17,11 +17,7 @@
* =============
*>
*> \verbatim
-*>
-*> ZGEQR computes a QR factorization of an M-by-N matrix A,
-*> using ZLATSQR when A is tall and skinny
-*> (M sufficiently greater than N), and otherwise ZGEQRT:
-*> A = Q * R .
+*> ZGEQR computes a QR factorization of an M-by-N matrix A.
*> \endverbatim
*
* Arguments:
@@ -46,7 +42,8 @@
*> On exit, the elements on and above the diagonal of the array
*> contain the min(M,N)-by-N upper trapezoidal matrix R
*> (R is upper triangular if M >= N);
-*> the elements below the diagonal represent Q (see Further Details).
+*> the elements below the diagonal are used to store part of the
+*> data structure to represent Q.
*> \endverbatim
*>
*> \param[in] LDA
@@ -55,47 +52,50 @@
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
-*> \param[out] WORK1
+*> \param[out] T
*> \verbatim
-*> WORK1 is COMPLEX*16 array, dimension (MAX(1,LWORK1))
-*> WORK1 contains part of the data structure used to store Q.
-*> WORK1(1): algorithm type = 1, to indicate output from
-*> ZLATSQR or ZGEQRT
-*> WORK1(2): optimum size of WORK1
-*> WORK1(3): minimum size of WORK1
-*> WORK1(4): row block size
-*> WORK1(5): column block size
-*> WORK1(6:LWORK1): data structure needed for Q, computed by
-*> CLATSQR or CGEQRT
+*> T is COMPLEX*16 array, dimension (MAX(5,TSIZE))
+*> On exit, if INFO = 0, T(1) returns optimal (or either minimal
+*> or optimal, if query is assumed) TSIZE. See TSIZE for details.
+*> Remaining T contains part of the data structure used to represent Q.
+*> If one wants to apply or construct Q, then one needs to keep T
+*> (in addition to A) and pass it to further subroutines.
*> \endverbatim
*>
-*> \param[in] LWORK1
+*> \param[in] TSIZE
*> \verbatim
-*> LWORK1 is INTEGER
-*> The dimension of the array WORK1.
-*> If LWORK1 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK1 and
-*> returns this value in WORK1(2), and calculates the minimum
-*> size of WORK1 and returns this value in WORK1(3).
-*> No error message related to LWORK1 is issued by XERBLA when
-*> LWORK1 = -1.
+*> TSIZE is INTEGER
+*> If TSIZE >= 5, the dimension of the array T.
+*> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If TSIZE = -1, the routine calculates optimal size of T for the
+*> optimum performance and returns this value in T(1).
+*> If TSIZE = -2, the routine calculates minimal size of T and
+*> returns this value in T(1).
*> \endverbatim
*>
-*> \param[out] WORK2
+*> \param[out] WORK
*> \verbatim
-*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK2))
+*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
-*> \param[in] LWORK2
+*> \param[in] LWORK
*> \verbatim
-*> LWORK2 is INTEGER
-*> The dimension of the array WORK2.
-*> If LWORK2 = -1, then a query is assumed. In this case the
-*> routine calculates the optimal size of WORK2 and
-*> returns this value in WORK2(1), and calculates the minimum
-*> size of WORK2 and returns this value in WORK2(2).
-*> No error message related to LWORK2 is issued by XERBLA when
-*> LWORK2 = -1.
+*> LWORK is INTEGER
+*> The dimension of the array WORK.
+*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
+*> only calculates the sizes of the T and WORK arrays, returns these
+*> values as the first entries of the T and WORK arrays, and no error
+*> message related to T or WORK is issued by XERBLA.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -113,26 +113,50 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \par Further Details:
-* =====================
+*> \par Further Details
+* ====================
*>
*> \verbatim
+*>
+*> The goal of the interface is to give maximum freedom to the developers for
+*> creating any QR factorization algorithm they wish. The triangular
+*> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
+*> and the array T can be used to store any relevant information for applying or
+*> constructing the Q factor. The WORK array can safely be discarded after exit.
+*>
+*> Caution: One should not expect the sizes of T and WORK to be the same from one
+*> LAPACK implementation to the other, or even from one execution to the other.
+*> A workspace query (for T and WORK) is needed at each execution. However,
+*> for a given execution, the size of T and WORK are fixed and will not change
+*> from one query to the next.
+*>
+*> \endverbatim
+*>
+*> \par Further Details particular to this LAPACK implementation:
+* ==============================================================
+*>
+*> \verbatim
+*>
+*> These details are particular for this LAPACK implementation. Users should not
+*> take them for granted. These details may change in the future, and are unlikely not
+*> true for another LAPACK implementation. These details are relevant if one wants
+*> to try to understand the code. They are not part of the interface.
+*>
+*> In this version,
+*>
+*> T(2): row block size (MB)
+*> T(3): column block size (NB)
+*> T(4:TSIZE): data structure needed for Q, computed by
+*> DLATSQR or DGEQRT
+*>
*> Depending on the matrix dimensions M and N, and row and column
*> block sizes MB and NB returned by ILAENV, GEQR will use either
*> LATSQR (if the matrix is tall-and-skinny) or GEQRT to compute
-*> the QR decomposition.
-*> The output of LATSQR or GEQRT representing Q is stored in A and in
-*> array WORK1(6:LWORK1) for later use.
-*> WORK1(2:5) contains the matrix dimensions M,N and block sizes MB,NB
-*> which are needed to interpret A and WORK1(6:LWORK1) for later use.
-*> WORK1(1)=1 indicates that the code needed to take WORK1(2:5) and
-*> decide whether LATSQR or GEQRT was used is the same as used below in
-*> GEQR. For a detailed description of A and WORK1(6:LWORK1), see
-*> Further Details in LATSQR or GEQRT.
+*> the QR factorization.
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE ZGEQR( M, N, A, LDA, WORK1, LWORK1, WORK2, LWORK2,
+ SUBROUTINE ZGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO)
*
* -- LAPACK computational routine (version 3.5.0) --
@@ -141,18 +165,18 @@
* November 2013
*
* .. Scalar Arguments ..
- INTEGER INFO, LDA, M, N, LWORK1, LWORK2
+ INTEGER INFO, LDA, M, N, TSIZE, LWORK
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), WORK1( * ), WORK2( * )
+ COMPLEX*16 A( LDA, * ), T( * ), WORK( * )
* ..
*
* =====================================================================
*
* ..
* .. Local Scalars ..
- LOGICAL LQUERY, LMINWS
- INTEGER MB, NB, I, II, KK, MINLW1, NBLCKS
+ LOGICAL LQUERY, LMINWS, MINT, MINW
+ INTEGER MB, NB, I, II, KK, MINTSZ, NBLCKS
* ..
* .. EXTERNAL FUNCTIONS ..
LOGICAL LSAME
@@ -172,7 +196,18 @@
*
INFO = 0
*
- LQUERY = ( LWORK1.EQ.-1 .OR. LWORK2.EQ.-1 )
+ LQUERY = ( TSIZE.EQ.-1 .OR. TSIZE.EQ.-2 .OR.
+ $ LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
+*
+ MINT = .FALSE.
+ IF ( TSIZE.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINT = .TRUE.
+ ENDIF
+*
+ MINW = .FALSE.
+ IF ( LWORK.NE.-1 .AND. ( TSIZE.EQ.-2 .OR. LWORK.EQ.-2 ) ) THEN
+ MINW = .TRUE.
+ ENDIF
*
* Determine the block size
*
@@ -183,34 +218,31 @@
MB = M
NB = 1
END IF
- IF( MB.GT.M.OR.MB.LE.N) MB = M
- IF( NB.GT.MIN(M,N).OR.NB.LT.1) NB = 1
- MINLW1 = N + 5
- IF ((MB.GT.N).AND.(M.GT.N)) THEN
- IF(MOD(M-N, MB-N).EQ.0) THEN
- NBLCKS = (M-N)/(MB-N)
+ IF( MB.GT.M .OR. MB.LE.N ) MB = M
+ IF( NB.GT.MIN( M, N ) .OR. NB.LT.1 ) NB = 1
+ MINTSZ = N + 5
+ IF ( MB.GT.N .AND. M.GT.N ) THEN
+ IF( MOD( M - N, MB - N ).EQ.0 ) THEN
+ NBLCKS = ( M - N ) / ( MB - N )
ELSE
- NBLCKS = (M-N)/(MB-N) + 1
+ NBLCKS = ( M - N ) / ( MB - N ) + 1
END IF
ELSE
NBLCKS = 1
END IF
*
-* Determine if the workspace size satisfies minimum size
+* Determine if the workspace size satisfies minimal size
*
LMINWS = .FALSE.
- IF((LWORK1.LT.MAX(1, NB*N*NBLCKS+5)
- $ .OR.(LWORK2.LT.NB*N)).AND.(LWORK2.GE.N).AND.(LWORK1.GT.N+5)
- $ .AND.(.NOT.LQUERY)) THEN
- IF (LWORK1.LT.MAX(1, NB * N * NBLCKS+5)) THEN
+ IF( ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) .OR. LWORK.LT.NB*N )
+ $ .AND. ( LWORK.GE.N ) .AND. ( TSIZE.GE.N + 5 )
+ $ .AND. ( .NOT.LQUERY ) ) THEN
+ IF ( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 ) ) THEN
LMINWS = .TRUE.
NB = 1
- END IF
- IF (LWORK1.LT.MAX(1, N * NBLCKS+5)) THEN
- LMINWS = .TRUE.
MB = M
END IF
- IF (LWORK2.LT.NB*N) THEN
+ IF ( LWORK.LT.NB*N ) THEN
LMINWS = .TRUE.
NB = 1
END IF
@@ -222,44 +254,50 @@
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
- ELSE IF( LWORK1.LT.MAX( 1, NB * N * NBLCKS + 5 )
- $ .AND.(.NOT.LQUERY).AND.(.NOT.LMINWS)) THEN
+ ELSE IF( TSIZE.LT.MAX( 1, NB*N*NBLCKS + 5 )
+ $ .AND. ( .NOT.LQUERY ) .AND. ( .NOT.LMINWS ) ) THEN
INFO = -6
- ELSE IF( (LWORK2.LT.MAX(1,N*NB)).AND.(.NOT.LQUERY)
- $ .AND.(.NOT.LMINWS)) THEN
+ ELSE IF( ( LWORK.LT.MAX( 1, N*NB ) ) .AND. ( .NOT.LQUERY )
+ $ .AND. ( .NOT.LMINWS ) ) THEN
INFO = -8
END IF
-
- IF( INFO.EQ.0) THEN
- WORK1(1) = 1
- WORK1(2) = NB * N * NBLCKS + 5
- WORK1(3) = MINLW1
- WORK1(4) = MB
- WORK1(5) = NB
- WORK2(1) = NB * N
- WORK2(2) = N
+*
+ IF( INFO.EQ.0 ) THEN
+ IF ( MINT ) THEN
+ T(1) = MINTSZ
+ ELSE
+ T(1) = NB*N*NBLCKS + 5
+ ENDIF
+ T(2) = MB
+ T(3) = NB
+ IF ( MINW ) THEN
+ WORK(1) = MAX( 1, N )
+ ELSE
+ WORK(1) = MAX( 1, NB*N )
+ ENDIF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGEQR', -INFO )
RETURN
ELSE IF (LQUERY) THEN
- RETURN
+ RETURN
END IF
*
* Quick return if possible
*
IF( MIN(M,N).EQ.0 ) THEN
- RETURN
+ RETURN
END IF
*
* The QR Decomposition
*
- IF((M.LE.N).OR.(MB.LE.N).OR.(MB.GE.M)) THEN
- CALL ZGEQRT( M, N, NB, A, LDA, WORK1(6), NB, WORK2, INFO)
+ IF( ( M.LE.N ) .OR. ( MB.LE.N ) .OR. ( MB.GE.M ) ) THEN
+ CALL ZGEQRT( M, N, NB, A, LDA, T(4), NB, WORK, INFO )
ELSE
- CALL ZLATSQR( M, N, MB, NB, A, LDA, WORK1(6), NB, WORK2,
- $ LWORK2, INFO)
+ CALL ZLATSQR( M, N, MB, NB, A, LDA, T(4), NB, WORK,
+ $ LWORK, INFO )
END IF
+ WORK(1) = MAX( 1, NB*N )
RETURN
*
* End of ZGEQR
diff --git a/SRC/zgetsls.f b/SRC/zgetsls.f
index d61b88c3..2a46fc1b 100644
--- a/SRC/zgetsls.f
+++ b/SRC/zgetsls.f
@@ -1,16 +1,15 @@
* Definition:
* ===========
*
-* SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
-* $ , WORK, LWORK, INFO )
-
+* SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+* $ WORK, LWORK, INFO )
*
* .. Scalar Arguments ..
* CHARACTER TRANS
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
+* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
* ..
*
*
@@ -19,10 +18,11 @@
*>
*> \verbatim
*>
-*> ZGETSLS solves overdetermined or underdetermined real linear systems
-*> involving an M-by-N matrix A, or its transpose, using a tall skinny
-*> QR or short wide LQfactorization of A. It is assumed that A has
-*> full rank.
+*> ZGETSLS solves overdetermined or underdetermined complex linear systems
+*> involving an M-by-N matrix A, using a tall skinny QR or short wide LQ
+*> factorization of A. It is assumed that A has full rank.
+*>
+*>
*>
*> The following options are provided:
*>
@@ -80,10 +80,8 @@
*> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the M-by-N matrix A.
*> On exit,
-*> if M >= N, A is overwritten by details of its QR
-*> factorization as returned by DGEQRF;
-*> if M < N, A is overwritten by details of its LQ
-*> factorization as returned by DGELQF.
+*> A is overwritten by details of its QR or LQ
+*> factorization as returned by ZGEQR or ZGEQL.
*> \endverbatim
*>
*> \param[in] LDA
@@ -97,21 +95,17 @@
*> B is COMPLEX*16 array, dimension (LDB,NRHS)
*> On entry, the matrix B of right hand side vectors, stored
*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
-*> if TRANS = 'T'.
+*> if TRANS = 'C'.
*> On exit, if INFO = 0, B is overwritten by the solution
*> vectors, stored columnwise:
*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
-*> squares solution vectors; the residual sum of squares for the
-*> solution in each column is given by the sum of squares of
-*> elements N+1 to M in that column;
+*> squares solution vectors.
*> if TRANS = 'N' and m < n, rows 1 to N of B contain the
*> minimum norm solution vectors;
-*> if TRANS = 'T' and m >= n, rows 1 to M of B contain the
+*> if TRANS = 'C' and m >= n, rows 1 to M of B contain the
*> minimum norm solution vectors;
-*> if TRANS = 'T' and m < n, rows 1 to M of B contain the
-*> least squares solution vectors; the residual sum of squares
-*> for the solution in each column is given by the sum of
-*> squares of elements M+1 to N in that column.
+*> if TRANS = 'C' and m < n, rows 1 to M of B contain the
+*> least squares solution vectors.
*> \endverbatim
*>
*> \param[in] LDB
@@ -122,23 +116,21 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
-*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
+*> or optimal, if query was assumed) LWORK.
+*> See LWORK for details.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
-*> LWORK >= max( 1, MN + max( MN, NRHS ) ).
-*> For optimal performance,
-*> LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
-*> where MN = min(M,N) and NB is the optimum block size.
-*>
-*> 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.
+*> If LWORK = -1 or -2, then a workspace query is assumed.
+*> If LWORK = -1, the routine calculates optimal size of WORK for the
+*> optimal performance and returns this value in WORK(1).
+*> If LWORK = -2, the routine calculates minimal size of WORK and
+*> returns this value in WORK(1).
*> \endverbatim
*>
*> \param[out] INFO
@@ -162,9 +154,11 @@
*
*> \date November 2011
*
+*> \ingroup complex16GEsolve
+*
* =====================================================================
- SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB
- $ , WORK, LWORK, INFO )
+ SUBROUTINE ZGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
+ $ WORK, LWORK, INFO )
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -173,10 +167,10 @@
*
* .. Scalar Arguments ..
CHARACTER TRANS
- INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, MB, NB
+ INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
+ COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
*
* ..
*
@@ -190,9 +184,11 @@
* ..
* .. Local Scalars ..
LOGICAL LQUERY, TRAN
- INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW, LW,
- $ SCLLEN, MNK, WSIZEO, WSIZEM, LW1, LW2, INFO2
+ INTEGER I, IASCL, IBSCL, J, MINMN, MAXMN, BROW,
+ $ SCLLEN, MNK, TSZ, TSZM, LW, LWM, LW1, LW2
+ $ INFO2
DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
+ COMPLEX*16 WRK1(5), WRK2
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -205,7 +201,7 @@
$ ZTRTRS, XERBLA, ZGELQ, ZGEMLQ
* ..
* .. Intrinsic Functions ..
- INTRINSIC DBLE, MAX, MIN
+ INTRINSIC DBLE, MAX, MIN, INT
* ..
* .. Executable Statements ..
*
@@ -217,7 +213,7 @@
MNK = MAX(MINMN,NRHS)
TRAN = LSAME( TRANS, 'C' )
*
- LQUERY = ( LWORK.EQ.-1 )
+ LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
$ LSAME( TRANS, 'C' ) ) ) THEN
INFO = -1
@@ -233,53 +229,58 @@
INFO = -8
END IF
*
- IF( INFO.EQ.0) THEN
+ IF( INFO.EQ.0 ) THEN
*
* Determine the block size and minimum LWORK
*
- IF ( M.GE.N ) THEN
- CALL ZGEQR( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
- ELSE
- CALL ZGELQ( M, N, A, LDA, WORK(1), -1, WORK(6), -1,
- $ INFO2)
- MB = INT(WORK(4))
- NB = INT(WORK(5))
- LW = INT(WORK(6))
- CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WORK(1),
- $ INT(WORK(2)), B, LDB, WORK(6), -1 , INFO2 )
- WSIZEO = INT(WORK(2))+MAX(LW,INT(WORK(6)))
- WSIZEM = INT(WORK(3))+MAX(LW,INT(WORK(6)))
- END IF
+ IF ( M.GE.N ) THEN
+ CALL ZGEQR( M, N, A, LDA, WRK1, -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL ZGEQR( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1 , INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, WRK1(1),
+ $ TSZM, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
+ ELSE
+ CALL ZGELQ( M, N, A, LDA, WRK1(1), -1, WRK2, -1, INFO2 )
+ TSZ = INT( WRK1(1) )
+ LW = INT( WRK2 )
+ CALL ZGELQ( M, N, A, LDA, WRK1(1), -2, WRK2, -2, INFO2 )
+ TSZM = INT( WRK1(1) )
+ LWM = INT( WRK2 )
+ CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -1, INFO2 )
+ LW = MAX( LW, INT(WRK2) )
+ CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, WRK1(1),
+ $ TSZ, B, LDB, WRK2, -2, INFO2 )
+ LWM = MAX( LWM, INT(WRK2) )
+ END IF
*
- IF((LWORK.LT.WSIZEO).AND.(.NOT.LQUERY)) THEN
+ IF( ( LWORK.LT.( TSZM + LWM ) ) .AND. (.NOT.LQUERY)) THEN
INFO=-10
END IF
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZGETSLS', -INFO )
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
- ELSE IF (LQUERY) THEN
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ END IF
+ IF ( LQUERY ) THEN
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
END IF
- IF(LWORK.LT.WSIZEO) THEN
- LW1=INT(WORK(3))
- LW2=MAX(LW,INT(WORK(6)))
+ IF( LWORK.LT.( TSZ + LW ) ) THEN
+ LW1 = TSZM
+ LW2 = LWM
ELSE
- LW1=INT(WORK(2))
- LW2=MAX(LW,INT(WORK(6)))
+ LW1 = TSZ
+ LW2 = LW
END IF
*
* Quick return if possible
@@ -346,9 +347,9 @@
*
* compute QR factorization of A
*
- CALL ZGEQR( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
- IF (.NOT.TRAN) THEN
+ CALL ZGEQR( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
+ IF ( .NOT.TRAN ) THEN
*
* Least-Squares Problem min || A * X - B ||
*
@@ -400,8 +401,8 @@
*
* Compute LQ factorization of A
*
- CALL ZGELQ( M, N, A, LDA, WORK(LW2+1), LW1
- $ , WORK(1), LW2, INFO )
+ CALL ZGELQ( M, N, A, LDA, WORK(LW2+1), LW1,
+ $ WORK(1), LW2, INFO )
*
* workspace at least M, optimally M*NB.
*
@@ -481,8 +482,7 @@
END IF
*
50 CONTINUE
- WORK( 1 ) = DBLE( WSIZEO )
- WORK( 2 ) = DBLE( WSIZEM )
+ WORK( 1 ) = DBLE( TSZ + LW )
RETURN
*
* End of ZGETSLS