Changes in TS QR API

1 files changed, 130 insertions, 93 deletions

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