Changes in TS QR API

1 files changed, 133 insertions, 95 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