APPLYING INTEL PATCHES sent to Julie on Feb 19th 2016 by Dima from INTEL (dmitry.g.baksheev@intel.com)

Subject: [PATCH 01/42] Fix ?GEJSV: handle M=N=0 case correctly, amend
 documentation. ISSUE: need LQUERY

- Complex kind functions operate on complex matrices
- JOBA cannot be 'N', and we have no parameter named JOBE
- Typos like CGEJSV uses CUNMQR, not SUNMQR
- Case M=0 N=0 shall zero IWORK(1:3) and WORK(1:7) anyway (F90 syntax!)
- Comments fixed to not corrupt doxygen output
- Missing fill of IWORK(3) = 0
- Math is written from capital letter :-)
- Typos like use ZGELQF, not ZGELQ

4 files changed, 64 insertions, 47 deletions

diff --git a/SRC/cgejsv.f b/SRC/cgejsv.f
index 271f3ace..60ef786b 100644
--- a/SRC/cgejsv.f
+++ b/SRC/cgejsv.f
@@ -39,14 +39,14 @@
 *>
 *> \verbatim
 *>
-*> CGEJSV computes the singular value decomposition (SVD) of a real M-by-N
+*> CGEJSV computes the singular value decomposition (SVD) of a complex M-by-N
 *> matrix [A], where M >= N. The SVD of [A] is written as
 *>
 *>              [A] = [U] * [SIGMA] * [V]^*,
 *>
 *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
-*> diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
-*> [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
+*> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and
+*> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are
 *> the singular values of [A]. The columns of [U] and [V] are the left and
 *> the right singular vectors of [A], respectively. The matrices [U] and [V]
 *> are computed and stored in the arrays U and V, respectively. The diagonal
@@ -221,7 +221,7 @@
 *>
 *> \param[out] U
 *> \verbatim
-*>          U is COMPLEX array, dimension ( LDU, N )
+*>          U is COMPLEX array, dimension ( LDU, N ) or ( LDU, M )
 *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
 *>                         the left singular vectors.
 *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
@@ -278,7 +278,7 @@
 *>          LWORK depends on the job:
 *>
 *>          1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*>            1.1 .. no scaled condition estimate required (JOBA.EQ.'N'):
+*>            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
 *>               LWORK >= 2*N+1. This is the minimal requirement.
 *>               ->> For optimal performance (blocked code) the optimal value
 *>               is LWORK >= N + (N+1)*NB. Here NB is the optimal
@@ -318,7 +318,7 @@
 *>            4.2. if JOBV.EQ.'J' the minimal requirement is 
 *>               LWORK >= 4*N+N*N.
 *>            In both cases, the allocated CWORK can accommodate blocked runs
-*>            of CGEQP3, CGEQRF, CGELQF, SUNMQR, CUNMLQ.
+*>            of CGEQP3, CGEQRF, CGELQF, CUNMQR, CUNMLQ.
 *> \endverbatim
 *>
 *> \param[out] RWORK
@@ -498,7 +498,7 @@
 *>     LAPACK Working note 170.
 *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
 *>     factorization software - a case study.
-*>     ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
+*>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
 *>     LAPACK Working note 176.
 *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 *>     QSVD, (H,K)-SVD computations.
@@ -636,7 +636,11 @@
 *
 *     Quick return for void matrix (Y3K safe)
 * #:)
-      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+         IWORK(1:3) = 0
+         RWORK(1:7) = 0
+         RETURN
+      ENDIF
 *
 *     Determine whether the matrix U should be M x N or M x M
 *
diff --git a/SRC/dgejsv.f b/SRC/dgejsv.f
index 66d79bfc..0387ec7d 100644
--- a/SRC/dgejsv.f
+++ b/SRC/dgejsv.f
@@ -52,7 +52,8 @@
 *> are computed and stored in the arrays U and V, respectively. The diagonal
 *> of [SIGMA] is computed and stored in the array SVA.
 *> DGEJSV can sometimes compute tiny singular values and their singular vectors much
-*> more accurately than other SVD routines, see below under Further Details.*> \endverbatim
+*> more accurately than other SVD routines, see below under Further Details.
+*> \endverbatim
 *
 *  Arguments:
 *  ==========
@@ -332,10 +333,10 @@
 *>          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
 *>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
-*>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ,
+*>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
 *>               DORMLQ. In general, the optimal length LWORK is computed as
 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), 
-*>                       N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
+*>                       N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
 *>
 *>          If SIGMA and the left singular vectors are needed
 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
@@ -589,7 +590,11 @@
 *
 *     Quick return for void matrix (Y3K safe)
 * #:)
-      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+         IWORK(1:3) = 0
+         WORK(1:7) = 0
+         RETURN
+      ENDIF
 *
 *     Determine whether the matrix U should be M x N or M x M
 *
@@ -715,6 +720,7 @@
             IWORK(1) = 0
             IWORK(2) = 0
          END IF
+         IWORK(3) = 0
          IF ( ERREST ) WORK(3) = ONE
          IF ( LSVEC .AND. RSVEC ) THEN
             WORK(4) = ONE
diff --git a/SRC/sgejsv.f b/SRC/sgejsv.f
index b849c9ed..8d4065e4 100644
--- a/SRC/sgejsv.f
+++ b/SRC/sgejsv.f
@@ -53,7 +53,6 @@
 *> of [SIGMA] is computed and stored in the array SVA.
 *> SGEJSV can sometimes compute tiny singular values and their singular vectors much
 *> more accurately than other SVD routines, see below under Further Details.
-
 *> \endverbatim
 *
 *  Arguments:
@@ -459,7 +458,7 @@
 *>     LAPACK Working note 170.
 *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
 *>     factorization software - a case study.
-*>     ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
+*>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
 *>     LAPACK Working note 176.
 *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 *>     QSVD, (H,K)-SVD computations.
@@ -591,7 +590,11 @@
 *
 *     Quick return for void matrix (Y3K safe)
 * #:)
-      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+         IWORK(1:3) = 0
+         WORK(1:7) = 0
+         RETURN
+      ENDIF
 *
 *     Determine whether the matrix U should be M x N or M x M
 *
@@ -717,6 +720,7 @@
             IWORK(1) = 0
             IWORK(2) = 0
          END IF
+         IWORK(3) = 0
          IF ( ERREST ) WORK(3) = ONE
          IF ( LSVEC .AND. RSVEC ) THEN
             WORK(4) = ONE
diff --git a/SRC/zgejsv.f b/SRC/zgejsv.f
index 84c9ea0a..43572816 100644
--- a/SRC/zgejsv.f
+++ b/SRC/zgejsv.f
@@ -39,21 +39,21 @@
 *>
 *> \verbatim
 *>
-* ZGEJSV computes the singular value decomposition (SVD) of a real M-by-N
-* matrix [A], where M >= N. The SVD of [A] is written as
-*
-*              [A] = [U] * [SIGMA] * [V]^*,
-*
-* where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
-* diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
-* [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
-* the singular values of [A]. The columns of [U] and [V] are the left and
-* the right singular vectors of [A], respectively. The matrices [U] and [V]
-* are computed and stored in the arrays U and V, respectively. The diagonal
-* of [SIGMA] is computed and stored in the array SVA.
-*
-*  Arguments:
-*  ==========
+*> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N
+*> matrix [A], where M >= N. The SVD of [A] is written as
+*>
+*>              [A] = [U] * [SIGMA] * [V]^*,
+*>
+*> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
+*> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and
+*> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are
+*> the singular values of [A]. The columns of [U] and [V] are the left and
+*> the right singular vectors of [A], respectively. The matrices [U] and [V]
+*> are computed and stored in the arrays U and V, respectively. The diagonal
+*> of [SIGMA] is computed and stored in the array SVA.
+*>
+*>  Arguments:
+*>  ==========
 *>
 *> \param[in] JOBA
 *> \verbatim
@@ -268,7 +268,6 @@
 *>
 *> \param[out] CWORK
 *> \verbatim
-*> CWORK (workspace)
 *>          CWORK is DOUBLE COMPLEX array, dimension at least LWORK.     
 *> \endverbatim
 *>
@@ -279,7 +278,7 @@
 *>          LWORK depends on the job:
 *>
 *>          1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*>            1.1 .. no scaled condition estimate required (JOBE.EQ.'N'):
+*>            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
 *>               LWORK >= 2*N+1. This is the minimal requirement.
 *>               ->> For optimal performance (blocked code) the optimal value
 *>               is LWORK >= N + (N+1)*NB. Here NB is the optimal
@@ -299,7 +298,7 @@
 *>             (JOBU.EQ.'N')
 *>            -> the minimal requirement is LWORK >= 3*N.
 *>            -> For optimal performance, LWORK >= max(N+(N+1)*NB, 3*N,2*N+N*NB),
-*>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
+*>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF,
 *>               ZUNMLQ. In general, the optimal length LWORK is computed as
 *>               LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZPOCON), N+LWORK(ZGESVJ),
 *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
@@ -491,19 +490,19 @@
 *>
 *> \verbatim
 *>
-* [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
-*     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
-*     LAPACK Working note 169.
-* [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
-*     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
-*     LAPACK Working note 170.
-* [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
-*     factorization software - a case study.
-*     ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
-*     LAPACK Working note 176.
-* [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
-*     QSVD, (H,K)-SVD computations.
-*     Department of Mathematics, University of Zagreb, 2008.
+*> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
+*>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
+*>     LAPACK Working note 169.
+*> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
+*>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
+*>     LAPACK Working note 170.
+*> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
+*>     factorization software - a case study.
+*>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
+*>     LAPACK Working note 176.
+*> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
+*>     QSVD, (H,K)-SVD computations.
+*>     Department of Mathematics, University of Zagreb, 2008.
 *> \endverbatim
 *
 *>  \par Bugs, examples and comments:
@@ -640,7 +639,11 @@
 *
 *     Quick return for void matrix (Y3K safe)
 * #:)
-      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+      IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+         IWORK(1:3) = 0
+         RWORK(1:7) = 0
+         RETURN
+      ENDIF
 *
 *     Determine whether the matrix U should be M x N or M x M
 *