diff options
author | konstantin.i.arturov <karturov@mugwort.jf.intel.com> | 2016-12-14 01:02:35 -0800 |
---|---|---|
committer | konstantin.i.arturov <karturov@mugwort.jf.intel.com> | 2016-12-14 01:02:35 -0800 |
commit | ae90d35f0d85e236fdf3fa77d508d8859d568f5b (patch) | |
tree | 02e6a596d5510e2d524df51e69a4fba9d861faff | |
parent | b80d7b3365caa35838a86e19615b13d913cb2fed (diff) | |
download | lapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.tar.gz lapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.tar.bz2 lapack-ae90d35f0d85e236fdf3fa77d508d8859d568f5b.zip |
Changes in TS QR API
-rw-r--r-- | SRC/cgelq.f | 228 | ||||
-rw-r--r-- | SRC/cgemlq.f | 148 | ||||
-rw-r--r-- | SRC/cgemqr.f | 135 | ||||
-rw-r--r-- | SRC/cgeqr.f | 223 | ||||
-rw-r--r-- | SRC/cgetsls.f | 163 | ||||
-rw-r--r-- | SRC/dgelq.f | 232 | ||||
-rw-r--r-- | SRC/dgemlq.f | 135 | ||||
-rw-r--r-- | SRC/dgemqr.f | 133 | ||||
-rw-r--r-- | SRC/dgeqr.f | 228 | ||||
-rw-r--r-- | SRC/dgetsls.f | 114 | ||||
-rw-r--r-- | SRC/sgelq.f | 228 | ||||
-rw-r--r-- | SRC/sgemlq.f | 140 | ||||
-rw-r--r-- | SRC/sgemqr.f | 132 | ||||
-rw-r--r-- | SRC/sgeqr.f | 222 | ||||
-rw-r--r-- | SRC/sgetsls.f | 121 | ||||
-rw-r--r-- | SRC/zgelq.f | 227 | ||||
-rw-r--r-- | SRC/zgemlq.f | 148 | ||||
-rw-r--r-- | SRC/zgemqr.f | 135 | ||||
-rw-r--r-- | SRC/zgeqr.f | 222 | ||||
-rw-r--r-- | SRC/zgetsls.f | 160 |
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 |