diff --git a/lapack-netlib/SRC/dbdsqr.f b/lapack-netlib/SRC/dbdsqr.f
index 93db95e7a..7d47fa282 100644
--- a/lapack-netlib/SRC/dbdsqr.f
+++ b/lapack-netlib/SRC/dbdsqr.f
@@ -166,7 +166,7 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is DOUBLE PRECISION array, dimension (4*N)
+*> WORK is DOUBLE PRECISION array, dimension (4*(N-1))
*> \endverbatim
*>
*> \param[out] INFO
diff --git a/lapack-netlib/SRC/dbdsvdx.f b/lapack-netlib/SRC/dbdsvdx.f
index 96fdb3d61..10d97a71f 100644
--- a/lapack-netlib/SRC/dbdsvdx.f
+++ b/lapack-netlib/SRC/dbdsvdx.f
@@ -165,7 +165,7 @@
*>
*> \param[out] Z
*> \verbatim
-*> Z is DOUBLE PRECISION array, dimension (2*N,K) )
+*> Z is DOUBLE PRECISION array, dimension (2*N,K)
*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
*> contain the singular vectors of the matrix B corresponding to
*> the selected singular values, with U in rows 1 to N and V
diff --git a/lapack-netlib/SRC/dcombssq.f b/lapack-netlib/SRC/dcombssq.f
new file mode 100644
index 000000000..79f6d95c9
--- /dev/null
+++ b/lapack-netlib/SRC/dcombssq.f
@@ -0,0 +1,92 @@
+*> \brief \b DCOMBSSQ adds two scaled sum of squares quantities.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DCOMBSSQ( V1, V2 )
+*
+* .. Array Arguments ..
+* DOUBLE PRECISION V1( 2 ), V2( 2 )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2.
+*> That is,
+*>
+*> V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq
+*> + V2_scale**2 * V2_sumsq
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in,out] V1
+*> \verbatim
+*> V1 is DOUBLE PRECISION array, dimension (2).
+*> The first scaled sum.
+*> V1(1) = V1_scale, V1(2) = V1_sumsq.
+*> \endverbatim
+*>
+*> \param[in] V2
+*> \verbatim
+*> V2 is DOUBLE PRECISION array, dimension (2).
+*> The second scaled sum.
+*> V2(1) = V2_scale, V2(2) = V2_sumsq.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2018
+*
+*> \ingroup OTHERauxiliary
+*
+* =====================================================================
+ SUBROUTINE DCOMBSSQ( V1, V2 )
+*
+* -- LAPACK auxiliary routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2018
+*
+* .. Array Arguments ..
+ DOUBLE PRECISION V1( 2 ), V2( 2 )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER ( ZERO = 0.0D+0 )
+* ..
+* .. Executable Statements ..
+*
+ IF( V1( 1 ).GE.V2( 1 ) ) THEN
+ IF( V1( 1 ).NE.ZERO ) THEN
+ V1( 2 ) = V1( 2 ) + ( V2( 1 ) / V1( 1 ) )**2 * V2( 2 )
+ END IF
+ ELSE
+ V1( 2 ) = V2( 2 ) + ( V1( 1 ) / V2( 1 ) )**2 * V1( 2 )
+ V1( 1 ) = V2( 1 )
+ END IF
+ RETURN
+*
+* End of DCOMBSSQ
+*
+ END
diff --git a/lapack-netlib/SRC/dgbrfsx.f b/lapack-netlib/SRC/dgbrfsx.f
index fb52d643f..76afb2d6a 100644
--- a/lapack-netlib/SRC/dgbrfsx.f
+++ b/lapack-netlib/SRC/dgbrfsx.f
@@ -308,7 +308,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -344,14 +344,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -359,9 +359,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.0D+0
-*> = 0.0 : No refinement is performed, and no error bounds are
+*> = 0.0: No refinement is performed, and no error bounds are
*> computed.
-*> = 1.0 : Use the double-precision refinement algorithm,
+*> = 1.0: Use the double-precision refinement algorithm,
*> possibly with doubled-single computations if the
*> compilation environment does not support DOUBLE
*> PRECISION.
diff --git a/lapack-netlib/SRC/dgbsvxx.f b/lapack-netlib/SRC/dgbsvxx.f
index 819d20c6d..058b20686 100644
--- a/lapack-netlib/SRC/dgbsvxx.f
+++ b/lapack-netlib/SRC/dgbsvxx.f
@@ -431,7 +431,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -467,14 +467,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -482,9 +482,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.0D+0
-*> = 0.0 : No refinement is performed, and no error bounds are
+*> = 0.0: No refinement is performed, and no error bounds are
*> computed.
-*> = 1.0 : Use the extra-precise refinement algorithm.
+*> = 1.0: Use the extra-precise refinement algorithm.
*> (other values are reserved for future use)
*>
*> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
diff --git a/lapack-netlib/SRC/dgebak.f b/lapack-netlib/SRC/dgebak.f
index 45a86ee57..10a78aa1a 100644
--- a/lapack-netlib/SRC/dgebak.f
+++ b/lapack-netlib/SRC/dgebak.f
@@ -47,10 +47,10 @@
*> \verbatim
*> JOB is CHARACTER*1
*> Specifies the type of backward transformation required:
-*> = 'N', do nothing, return immediately;
-*> = 'P', do backward transformation for permutation only;
-*> = 'S', do backward transformation for scaling only;
-*> = 'B', do backward transformations for both permutation and
+*> = 'N': do nothing, return immediately;
+*> = 'P': do backward transformation for permutation only;
+*> = 'S': do backward transformation for scaling only;
+*> = 'B': do backward transformations for both permutation and
*> scaling.
*> JOB must be the same as the argument JOB supplied to DGEBAL.
*> \endverbatim
diff --git a/lapack-netlib/SRC/dgeesx.f b/lapack-netlib/SRC/dgeesx.f
index 26042a5f9..a08104d3d 100644
--- a/lapack-netlib/SRC/dgeesx.f
+++ b/lapack-netlib/SRC/dgeesx.f
@@ -583,7 +583,9 @@
IF( N.GT.I+1 )
$ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
$ A( I+1, I+2 ), LDA )
- CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
+ IF( WANTVS ) THEN
+ CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
+ END IF
A( I, I+1 ) = A( I+1, I )
A( I+1, I ) = ZERO
END IF
diff --git a/lapack-netlib/SRC/dgejsv.f b/lapack-netlib/SRC/dgejsv.f
index 25ed248d0..a30cfab87 100644
--- a/lapack-netlib/SRC/dgejsv.f
+++ b/lapack-netlib/SRC/dgejsv.f
@@ -82,7 +82,7 @@
*> desirable, then this option is advisable. The input matrix A
*> is preprocessed with QR factorization with FULL (row and
*> column) pivoting.
-*> = 'G' Computation as with 'F' with an additional estimate of the
+*> = 'G': Computation as with 'F' with an additional estimate of the
*> condition number of B, where A=D*B. If A has heavily weighted
*> rows, then using this condition number gives too pessimistic
*> error bound.
@@ -133,7 +133,7 @@
*> specified range. If A .NE. 0 is scaled so that the largest singular
*> value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
*> the licence to kill columns of A whose norm in c*A is less than
-*> DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
+*> DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
*> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
*> = 'N': Do not kill small columns of c*A. This option assumes that
*> BLAS and QR factorizations and triangular solvers are
@@ -230,7 +230,7 @@
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of
*> the left singular vectors, including an ONB
*> of the orthogonal complement of the Range(A).
-*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
+*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
*> then U is used as workspace if the procedure
*> replaces A with A^t. In that case, [V] is computed
*> in U as left singular vectors of A^t and then
@@ -252,7 +252,7 @@
*> V is DOUBLE PRECISION array, dimension ( LDV, N )
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
*> the right singular vectors;
-*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
+*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
*> then V is used as workspace if the pprocedure
*> replaces A with A^t. In that case, [U] is computed
*> in V as right singular vectors of A^t and then
@@ -272,13 +272,13 @@
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (LWORK)
-*> On exit, if N.GT.0 .AND. M.GT.0 (else not referenced),
+*> On exit, if N > 0 .AND. M > 0 (else not referenced),
*> WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
*> that SCALE*SVA(1:N) are the computed singular values
*> of A. (See the description of SVA().)
*> WORK(2) = See the description of WORK(1).
*> WORK(3) = SCONDA is an estimate for the condition number of
-*> column equilibrated A. (If JOBA .EQ. 'E' or 'G')
+*> column equilibrated A. (If JOBA = 'E' or 'G')
*> SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
*> It is computed using DPOCON. It holds
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
@@ -297,7 +297,7 @@
*> triangular factor in the first QR factorization.
*> WORK(5) = an estimate of the scaled condition number of the
*> triangular factor in the second QR factorization.
-*> The following two parameters are computed if JOBT .EQ. 'T'.
+*> The following two parameters are computed if JOBT = 'T'.
*> They are provided for a developer/implementer who is familiar
*> with the details of the method.
*>
@@ -313,8 +313,8 @@
*> Length of WORK to confirm proper allocation of work space.
*> LWORK depends on the job:
*>
-*> If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*> -> .. no scaled condition estimate required (JOBE.EQ.'N'):
+*> If only SIGMA is needed (JOBU = 'N', JOBV = 'N') and
+*> -> .. no scaled condition estimate required (JOBE = 'N'):
*> LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
*> ->> For optimal performance (blocked code) the optimal value
*> is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
@@ -330,7 +330,7 @@
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
*> N+N*N+LWORK(DPOCON),7).
*>
-*> If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
+*> If SIGMA and the right singular vectors are needed (JOBV = '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, DGELQF,
@@ -341,19 +341,19 @@
*> If SIGMA and the left singular vectors are needed
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*> -> For optimal performance:
-*> if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
-*> if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
+*> if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
+*> if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
*> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
*> In general, the optimal length LWORK is computed as
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
-*> Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or
-*> M*NB (for JOBU.EQ.'F').
+*> Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
+*> M*NB (for JOBU = 'F').
*>
-*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
-*> -> if JOBV.EQ.'V'
+*> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
+*> -> if JOBV = 'V'
*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
-*> -> if JOBV.EQ.'J' the minimal requirement is
+*> -> if JOBV = 'J' the minimal requirement is
*> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
*> -> For optimal performance, LWORK should be additionally
*> larger than N+M*NB, where NB is the optimal block size
@@ -369,7 +369,7 @@
*> of JOBA and JOBR.
*> IWORK(2) = the number of the computed nonzero singular values
*> IWORK(3) = if nonzero, a warning message:
-*> If IWORK(3).EQ.1 then some of the column norms of A
+*> If IWORK(3) = 1 then some of the column norms of A
*> were denormalized floats. The requested high accuracy
*> is not warranted by the data.
*> \endverbatim
@@ -377,10 +377,10 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value.
-*> = 0 : successful exit;
-*> > 0 : DGEJSV did not converge in the maximal allowed number
-*> of sweeps. The computed values may be inaccurate.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value.
+*> = 0: successful exit;
+*> > 0: DGEJSV did not converge in the maximal allowed number
+*> of sweeps. The computed values may be inaccurate.
*> \endverbatim
*
* Authors:
@@ -953,7 +953,7 @@
IF ( L2ABER ) THEN
* Standard absolute error bound suffices. All sigma_i with
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
-* agressive enforcement of lower numerical rank by introducing a
+* aggressive enforcement of lower numerical rank by introducing a
* backward error of the order of N*EPSLN*||A||.
TEMP1 = DSQRT(DBLE(N))*EPSLN
DO 3001 p = 2, N
@@ -965,7 +965,7 @@
3001 CONTINUE
3002 CONTINUE
ELSE IF ( L2RANK ) THEN
-* .. similarly as above, only slightly more gentle (less agressive).
+* .. similarly as above, only slightly more gentle (less aggressive).
* Sudden drop on the diagonal of R1 is used as the criterion for
* close-to-rank-deficient.
TEMP1 = DSQRT(SFMIN)
@@ -1294,7 +1294,7 @@
CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
$ WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
CONDR1 = ONE / DSQRT(TEMP1)
-* .. here need a second oppinion on the condition number
+* .. here need a second opinion on the condition number
* .. then assume worst case scenario
* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
* more conservative <=> CONDR1 .LT. DSQRT(DBLE(N))
@@ -1335,7 +1335,7 @@
ELSE
*
* .. ill-conditioned case: second QRF with pivoting
-* Note that windowed pivoting would be equaly good
+* Note that windowed pivoting would be equally good
* numerically, and more run-time efficient. So, in
* an optimal implementation, the next call to DGEQP3
* should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
@@ -1388,7 +1388,7 @@
*
IF ( CONDR2 .GE. COND_OK ) THEN
* .. save the Householder vectors used for Q3
-* (this overwrittes the copy of R2, as it will not be
+* (this overwrites the copy of R2, as it will not be
* needed in this branch, but it does not overwritte the
* Huseholder vectors of Q2.).
CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
@@ -1638,7 +1638,7 @@
*
* This branch deploys a preconditioned Jacobi SVD with explicitly
* accumulated rotations. It is included as optional, mainly for
-* experimental purposes. It does perfom well, and can also be used.
+* experimental purposes. It does perform well, and can also be used.
* In this implementation, this branch will be automatically activated
* if the condition number sigma_max(A) / sigma_min(A) is predicted
* to be greater than the overflow threshold. This is because the
diff --git a/lapack-netlib/SRC/dgelq.f b/lapack-netlib/SRC/dgelq.f
index ece645079..fc14d892f 100644
--- a/lapack-netlib/SRC/dgelq.f
+++ b/lapack-netlib/SRC/dgelq.f
@@ -1,3 +1,4 @@
+*> \brief \b DGELQ
*
* Definition:
* ===========
@@ -17,7 +18,17 @@
* =============
*>
*> \verbatim
-*> DGELQ computes a LQ factorization of an M-by-N matrix A.
+*>
+*> DGELQ computes an LQ factorization of a real M-by-N matrix A:
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a N-by-N orthogonal matrix;
+*> L is an lower-triangular M-by-M matrix;
+*> 0 is a M-by-(N-M) zero matrix, if M < N.
+*>
*> \endverbatim
*
* Arguments:
@@ -138,7 +149,7 @@
*> \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
+*> take them for granted. These details may change in the future, and are not likely
*> 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.
*>
@@ -159,10 +170,10 @@
SUBROUTINE DGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, TSIZE, LWORK
diff --git a/lapack-netlib/SRC/dgelq2.f b/lapack-netlib/SRC/dgelq2.f
index 04aa57fc1..a6c835de4 100644
--- a/lapack-netlib/SRC/dgelq2.f
+++ b/lapack-netlib/SRC/dgelq2.f
@@ -33,8 +33,16 @@
*>
*> \verbatim
*>
-*> DGELQ2 computes an LQ factorization of a real m by n matrix A:
-*> A = L * Q.
+*> DGELQ2 computes an LQ factorization of a real m-by-n matrix A:
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a n-by-n orthogonal matrix;
+*> L is an lower-triangular m-by-m matrix;
+*> 0 is a m-by-(n-m) zero matrix, if m < n.
+*>
*> \endverbatim
*
* Arguments:
@@ -96,7 +104,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -121,10 +129,10 @@
* =====================================================================
SUBROUTINE DGELQ2( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/dgelqf.f b/lapack-netlib/SRC/dgelqf.f
index 834c47168..4b11761f6 100644
--- a/lapack-netlib/SRC/dgelqf.f
+++ b/lapack-netlib/SRC/dgelqf.f
@@ -34,7 +34,15 @@
*> \verbatim
*>
*> DGELQF computes an LQ factorization of a real M-by-N matrix A:
-*> A = L * Q.
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a N-by-N orthogonal matrix;
+*> L is an lower-triangular M-by-M matrix;
+*> 0 is a M-by-(N-M) zero matrix, if M < N.
+*>
*> \endverbatim
*
* Arguments:
@@ -110,7 +118,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -135,10 +143,10 @@
* =====================================================================
SUBROUTINE DGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/dgemlq.f b/lapack-netlib/SRC/dgemlq.f
index bb6b2868f..dea693c24 100644
--- a/lapack-netlib/SRC/dgemlq.f
+++ b/lapack-netlib/SRC/dgemlq.f
@@ -1,3 +1,4 @@
+*> \brief \b DGEMLQ
*
* Definition:
* ===========
@@ -144,7 +145,7 @@
*> \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
+*> take them for granted. These details may change in the future, and are not likely
*> 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.
*>
diff --git a/lapack-netlib/SRC/dgemqr.f b/lapack-netlib/SRC/dgemqr.f
index 8509b13d9..0f7a42233 100644
--- a/lapack-netlib/SRC/dgemqr.f
+++ b/lapack-netlib/SRC/dgemqr.f
@@ -1,3 +1,4 @@
+*> \brief \b DGEMQR
*
* Definition:
* ===========
@@ -144,7 +145,7 @@
*> \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
+*> take them for granted. These details may change in the future, and are not likely
*> 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.
*>
diff --git a/lapack-netlib/SRC/dgeqr.f b/lapack-netlib/SRC/dgeqr.f
index d0a1a18f9..0bff5d1f9 100644
--- a/lapack-netlib/SRC/dgeqr.f
+++ b/lapack-netlib/SRC/dgeqr.f
@@ -1,3 +1,4 @@
+*> \brief \b DGEQR
*
* Definition:
* ===========
@@ -17,7 +18,18 @@
* =============
*>
*> \verbatim
-*> DGEQR computes a QR factorization of an M-by-N matrix A.
+*>
+*> DGEQR computes a QR factorization of a real M-by-N matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -138,7 +150,7 @@
*> \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
+*> take them for granted. These details may change in the future, and are not likely
*> 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.
*>
@@ -160,10 +172,10 @@
SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, TSIZE, LWORK
diff --git a/lapack-netlib/SRC/dgeqr2.f b/lapack-netlib/SRC/dgeqr2.f
index c1e91e9bd..9ce1feb25 100644
--- a/lapack-netlib/SRC/dgeqr2.f
+++ b/lapack-netlib/SRC/dgeqr2.f
@@ -33,8 +33,17 @@
*>
*> \verbatim
*>
-*> DGEQR2 computes a QR factorization of a real m by n matrix A:
-*> A = Q * R.
+*> DGEQR2 computes a QR factorization of a real m-by-n matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a m-by-m orthogonal matrix;
+*> R is an upper-triangular n-by-n matrix;
+*> 0 is a (m-n)-by-n zero matrix, if m > n.
+*>
*> \endverbatim
*
* Arguments:
@@ -96,7 +105,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -121,10 +130,10 @@
* =====================================================================
SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/dgeqr2p.f b/lapack-netlib/SRC/dgeqr2p.f
index 921f79921..9b81ccb33 100644
--- a/lapack-netlib/SRC/dgeqr2p.f
+++ b/lapack-netlib/SRC/dgeqr2p.f
@@ -33,8 +33,18 @@
*>
*> \verbatim
*>
-*> DGEQR2P computes a QR factorization of a real m by n matrix A:
-*> A = Q * R. The diagonal entries of R are nonnegative.
+*> DGEQR2P computes a QR factorization of a real m-by-n matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a m-by-m orthogonal matrix;
+*> R is an upper-triangular n-by-n matrix with nonnegative diagonal
+*> entries;
+*> 0 is a (m-n)-by-n zero matrix, if m > n.
+*>
*> \endverbatim
*
* Arguments:
@@ -97,7 +107,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -124,10 +134,10 @@
* =====================================================================
SUBROUTINE DGEQR2P( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/dgeqrf.f b/lapack-netlib/SRC/dgeqrf.f
index 83d7d8dd7..98666221f 100644
--- a/lapack-netlib/SRC/dgeqrf.f
+++ b/lapack-netlib/SRC/dgeqrf.f
@@ -34,7 +34,16 @@
*> \verbatim
*>
*> DGEQRF computes a QR factorization of a real M-by-N matrix A:
-*> A = Q * R.
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -111,7 +120,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -136,10 +145,10 @@
* =====================================================================
SUBROUTINE DGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/dgeqrfp.f b/lapack-netlib/SRC/dgeqrfp.f
index d182f98c9..5cf4069ed 100644
--- a/lapack-netlib/SRC/dgeqrfp.f
+++ b/lapack-netlib/SRC/dgeqrfp.f
@@ -33,8 +33,18 @@
*>
*> \verbatim
*>
-*> DGEQRFP computes a QR factorization of a real M-by-N matrix A:
-*> A = Q * R. The diagonal entries of R are nonnegative.
+*> DGEQR2P computes a QR factorization of a real M-by-N matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix with nonnegative diagonal
+*> entries;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -112,7 +122,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup doubleGEcomputational
*
@@ -139,10 +149,10 @@
* =====================================================================
SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/dgerfsx.f b/lapack-netlib/SRC/dgerfsx.f
index aafca8d10..495ea1726 100644
--- a/lapack-netlib/SRC/dgerfsx.f
+++ b/lapack-netlib/SRC/dgerfsx.f
@@ -283,7 +283,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -319,14 +319,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -334,9 +334,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.0D+0
-*> = 0.0 : No refinement is performed, and no error bounds are
+*> = 0.0: No refinement is performed, and no error bounds are
*> computed.
-*> = 1.0 : Use the double-precision refinement algorithm,
+*> = 1.0: Use the double-precision refinement algorithm,
*> possibly with doubled-single computations if the
*> compilation environment does not support DOUBLE
*> PRECISION.
diff --git a/lapack-netlib/SRC/dgesc2.f b/lapack-netlib/SRC/dgesc2.f
index 2f01a762f..72d8a38f0 100644
--- a/lapack-netlib/SRC/dgesc2.f
+++ b/lapack-netlib/SRC/dgesc2.f
@@ -90,7 +90,7 @@
*> \verbatim
*> SCALE is DOUBLE PRECISION
*> On exit, SCALE contains the scale factor. SCALE is chosen
-*> 0 <= SCALE <= 1 to prevent owerflow in the solution.
+*> 0 <= SCALE <= 1 to prevent overflow in the solution.
*> \endverbatim
*
* Authors:
@@ -151,7 +151,7 @@
* ..
* .. Executable Statements ..
*
-* Set constant to control owerflow
+* Set constant to control overflow
*
EPS = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' ) / EPS
diff --git a/lapack-netlib/SRC/dgesdd.f b/lapack-netlib/SRC/dgesdd.f
index 926607f98..0218900d2 100644
--- a/lapack-netlib/SRC/dgesdd.f
+++ b/lapack-netlib/SRC/dgesdd.f
@@ -322,7 +322,7 @@
*
IF( WNTQN ) THEN
* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
-* keep 7*N for backwards compatability.
+* keep 7*N for backwards compatibility.
BDSPAC = 7*N
ELSE
BDSPAC = 3*N*N + 4*N
@@ -448,7 +448,7 @@
*
IF( WNTQN ) THEN
* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
-* keep 7*N for backwards compatability.
+* keep 7*N for backwards compatibility.
BDSPAC = 7*M
ELSE
BDSPAC = 3*M*M + 4*M
diff --git a/lapack-netlib/SRC/dgesvdq.f b/lapack-netlib/SRC/dgesvdq.f
new file mode 100644
index 000000000..e495d2bf9
--- /dev/null
+++ b/lapack-netlib/SRC/dgesvdq.f
@@ -0,0 +1,1385 @@
+*> \brief DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGESVDQ + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
+* S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
+* WORK, LWORK, RWORK, LRWORK, INFO )
+*
+* .. Scalar Arguments ..
+* IMPLICIT NONE
+* CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
+* INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
+* INFO
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
+* DOUBLE PRECISION S( * ), RWORK( * )
+* INTEGER IWORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DGESVDQ computes the singular value decomposition (SVD) of a real
+*> M-by-N matrix A, where M >= N. The SVD of A is written as
+*> [++] [xx] [x0] [xx]
+*> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx]
+*> [++] [xx]
+*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOBA
+*> \verbatim
+*> JOBA is CHARACTER*1
+*> Specifies the level of accuracy in the computed SVD
+*> = 'A' The requested accuracy corresponds to having the backward
+*> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
+*> where EPS = DLAMCH('Epsilon'). This authorises DGESVDQ to
+*> truncate the computed triangular factor in a rank revealing
+*> QR factorization whenever the truncated part is below the
+*> threshold of the order of EPS * ||A||_F. This is aggressive
+*> truncation level.
+*> = 'M' Similarly as with 'A', but the truncation is more gentle: it
+*> is allowed only when there is a drop on the diagonal of the
+*> triangular factor in the QR factorization. This is medium
+*> truncation level.
+*> = 'H' High accuracy requested. No numerical rank determination based
+*> on the rank revealing QR factorization is attempted.
+*> = 'E' Same as 'H', and in addition the condition number of column
+*> scaled A is estimated and returned in RWORK(1).
+*> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
+*> \endverbatim
+*>
+*> \param[in] JOBP
+*> \verbatim
+*> JOBP is CHARACTER*1
+*> = 'P' The rows of A are ordered in decreasing order with respect to
+*> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
+*> of extra data movement. Recommended for numerical robustness.
+*> = 'N' No row pivoting.
+*> \endverbatim
+*>
+*> \param[in] JOBR
+*> \verbatim
+*> JOBR is CHARACTER*1
+*> = 'T' After the initial pivoted QR factorization, DGESVD is applied to
+*> the transposed R**T of the computed triangular factor R. This involves
+*> some extra data movement (matrix transpositions). Useful for
+*> experiments, research and development.
+*> = 'N' The triangular factor R is given as input to DGESVD. This may be
+*> preferred as it involves less data movement.
+*> \endverbatim
+*>
+*> \param[in] JOBU
+*> \verbatim
+*> JOBU is CHARACTER*1
+*> = 'A' All M left singular vectors are computed and returned in the
+*> matrix U. See the description of U.
+*> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
+*> in the matrix U. See the description of U.
+*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
+*> vectors are computed and returned in the matrix U.
+*> = 'F' The N left singular vectors are returned in factored form as the
+*> product of the Q factor from the initial QR factorization and the
+*> N left singular vectors of (R**T , 0)**T. If row pivoting is used,
+*> then the necessary information on the row pivoting is stored in
+*> IWORK(N+1:N+M-1).
+*> = 'N' The left singular vectors are not computed.
+*> \endverbatim
+*>
+*> \param[in] JOBV
+*> \verbatim
+*> JOBV is CHARACTER*1
+*> = 'A', 'V' All N right singular vectors are computed and returned in
+*> the matrix V.
+*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
+*> vectors are computed and returned in the matrix V. This option is
+*> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
+*> = 'N' The right singular vectors are not computed.
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of rows of the input matrix A. M >= 0.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of columns of the input matrix A. M >= N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array of dimensions LDA x N
+*> On entry, the input matrix A.
+*> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
+*> the Householder vectors as stored by DGEQP3. If JOBU = 'F', these Householder
+*> vectors together with WORK(1:N) can be used to restore the Q factors from
+*> the initial pivoted QR factorization of A. See the description of U.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER.
+*> The leading dimension of the array A. LDA >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] S
+*> \verbatim
+*> S is DOUBLE PRECISION array of dimension N.
+*> The singular values of A, ordered so that S(i) >= S(i+1).
+*> \endverbatim
+*>
+*> \param[out] U
+*> \verbatim
+*> U is DOUBLE PRECISION array, dimension
+*> LDU x M if JOBU = 'A'; see the description of LDU. In this case,
+*> on exit, U contains the M left singular vectors.
+*> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
+*> case, U contains the leading N or the leading NUMRANK left singular vectors.
+*> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
+*> contains N x N orthogonal matrix that can be used to form the left
+*> singular vectors.
+*> If JOBU = 'N', U is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDU
+*> \verbatim
+*> LDU is INTEGER.
+*> The leading dimension of the array U.
+*> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M).
+*> If JOBU = 'F', LDU >= max(1,N).
+*> Otherwise, LDU >= 1.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*> V is DOUBLE PRECISION array, dimension
+*> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
+*> If JOBV = 'A', or 'V', V contains the N-by-N orthogonal matrix V**T;
+*> If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right
+*> singular vectors, stored rowwise, of the NUMRANK largest singular values).
+*> If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
+*> If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*> LDV is INTEGER
+*> The leading dimension of the array V.
+*> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
+*> Otherwise, LDV >= 1.
+*> \endverbatim
+*>
+*> \param[out] NUMRANK
+*> \verbatim
+*> NUMRANK is INTEGER
+*> NUMRANK is the numerical rank first determined after the rank
+*> revealing QR factorization, following the strategy specified by the
+*> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
+*> leading singular values and vectors are then requested in the call
+*> of DGESVD. The final value of NUMRANK might be further reduced if
+*> some singular values are computed as zeros.
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (max(1, LIWORK)).
+*> On exit, IWORK(1:N) contains column pivoting permutation of the
+*> rank revealing QR factorization.
+*> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
+*> of row swaps used in row pivoting. These can be used to restore the
+*> left singular vectors in the case JOBU = 'F'.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> LIWORK(1) returns the minimal LIWORK.
+*> \endverbatim
+*>
+*> \param[in] LIWORK
+*> \verbatim
+*> LIWORK is INTEGER
+*> The dimension of the array IWORK.
+*> LIWORK >= N + M - 1, if JOBP = 'P' and JOBA .NE. 'E';
+*> LIWORK >= N if JOBP = 'N' and JOBA .NE. 'E';
+*> LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E';
+*> LIWORK >= N + N if JOBP = 'N' and JOBA = 'E'.
+*
+*> If LIWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (max(2, LWORK)), used as a workspace.
+*> On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters
+*> needed to recover the Q factor from the QR factorization computed by
+*> DGEQP3.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> WORK(1) returns the optimal LWORK, and
+*> WORK(2) returns the minimal LWORK.
+*> \endverbatim
+*>
+*> \param[in,out] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK. It is determined as follows:
+*> Let LWQP3 = 3*N+1, LWCON = 3*N, and let
+*> LWORQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
+*> { MAX( M, 1 ), if JOBU = 'A'
+*> LWSVD = MAX( 5*N, 1 )
+*> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
+*> LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
+*> Then the minimal value of LWORK is:
+*> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
+*> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
+*> and a scaled condition estimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left
+*> singular vectors are requested;
+*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left
+*> singular vectors are requested, and also
+*> a scaled condition estimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
+*> singular vectors are requested;
+*> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
+*> singular vectors are requested, and also
+*> a scaled condition etimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R';
+*> independent of JOBR;
+*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
+*> JOBV = 'R' and, also a scaled condition
+*> estimate requested; independent of JOBR;
+*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
+*> full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='N'
+*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
+*> if the full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='N', and also a scaled condition number estimate
+*> requested.
+*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
+*> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
+*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
+*> if the full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='T', and also a scaled condition number estimate
+*> requested.
+*> Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).
+*>
+*> If LWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension (max(1, LRWORK)).
+*> On exit,
+*> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
+*> number of column scaled A. If A = C * D where D is diagonal and C
+*> has unit columns in the Euclidean norm, then, assuming full column rank,
+*> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
+*> Otherwise, RWORK(1) = -1.
+*> 2. RWORK(2) contains the number of singular values computed as
+*> exact zeros in DGESVD applied to the upper triangular or trapeziodal
+*> R (from the initial QR factorization). In case of early exit (no call to
+*> DGESVD, such as in the case of zero matrix) RWORK(2) = -1.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> RWORK(1) returns the minimal LRWORK.
+*> \endverbatim
+*>
+*> \param[in] LRWORK
+*> \verbatim
+*> LRWORK is INTEGER.
+*> The dimension of the array RWORK.
+*> If JOBP ='P', then LRWORK >= MAX(2, M).
+*> Otherwise, LRWORK >= 2
+*
+*> If LRWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit.
+*> < 0: if INFO = -i, the i-th argument had an illegal value.
+*> > 0: if DBDSQR did not converge, INFO specifies how many superdiagonals
+*> of an intermediate bidiagonal form B (computed in DGESVD) did not
+*> converge to zero.
+*> \endverbatim
+*
+*> \par Further Details:
+* ========================
+*>
+*> \verbatim
+*>
+*> 1. The data movement (matrix transpose) is coded using simple nested
+*> DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
+*> Those DO-loops are easily identified in this source code - by the CONTINUE
+*> statements labeled with 11**. In an optimized version of this code, the
+*> nested DO loops should be replaced with calls to an optimized subroutine.
+*> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
+*> column norm overflow. This is the minial precaution and it is left to the
+*> SVD routine (CGESVD) to do its own preemptive scaling if potential over-
+*> or underflows are detected. To avoid repeated scanning of the array A,
+*> an optimal implementation would do all necessary scaling before calling
+*> CGESVD and the scaling in CGESVD can be switched off.
+*> 3. Other comments related to code optimization are given in comments in the
+*> code, enlosed in [[double brackets]].
+*> \endverbatim
+*
+*> \par Bugs, examples and comments
+* ===========================
+*
+*> \verbatim
+*> Please report all bugs and send interesting examples and/or comments to
+*> drmac@math.hr. Thank you.
+*> \endverbatim
+*
+*> \par References
+* ===============
+*
+*> \verbatim
+*> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
+*> Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
+*> 44(1): 11:1-11:30 (2017)
+*>
+*> SIGMA library, xGESVDQ section updated February 2016.
+*> Developed and coded by Zlatko Drmac, Department of Mathematics
+*> University of Zagreb, Croatia, drmac@math.hr
+*> \endverbatim
+*
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*> Developed and coded by Zlatko Drmac, Department of Mathematics
+*> University of Zagreb, Croatia, drmac@math.hr
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2018
+*
+*> \ingroup doubleGEsing
+*
+* =====================================================================
+ SUBROUTINE DGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
+ $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
+ $ WORK, LWORK, RWORK, LRWORK, INFO )
+* .. Scalar Arguments ..
+ IMPLICIT NONE
+ CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
+ INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
+ $ INFO
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
+ DOUBLE PRECISION S( * ), RWORK( * )
+ INTEGER IWORK( * )
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+* .. Local Scalars ..
+ INTEGER IERR, IWOFF, NR, N1, OPTRATIO, p, q
+ INTEGER LWCON, LWQP3, LWRK_DGELQF, LWRK_DGESVD, LWRK_DGESVD2,
+ $ LWRK_DGEQP3, LWRK_DGEQRF, LWRK_DORMLQ, LWRK_DORMQR,
+ $ LWRK_DORMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWORQ,
+ $ LWORQ2, LWORLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
+ $ IMINWRK, RMINWRK
+ LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
+ $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
+ $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
+ DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
+* .. Local Arrays
+ DOUBLE PRECISION RDUMMY(1)
+* ..
+* .. External Subroutines (BLAS, LAPACK)
+ EXTERNAL DGELQF, DGEQP3, DGEQRF, DGESVD, DLACPY, DLAPMT,
+ $ DLASCL, DLASET, DLASWP, DSCAL, DPOCON, DORMLQ,
+ $ DORMQR, XERBLA
+* ..
+* .. External Functions (BLAS, LAPACK)
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ DOUBLE PRECISION DLANGE, DNRM2, DLAMCH
+ EXTERNAL DLANGE, LSAME, IDAMAX, DNRM2, DLAMCH
+* ..
+* .. Intrinsic Functions ..
+*
+ INTRINSIC ABS, MAX, MIN, DBLE, SQRT
+*
+* Test the input arguments
+*
+ WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' )
+ WNTUR = LSAME( JOBU, 'R' )
+ WNTUA = LSAME( JOBU, 'A' )
+ WNTUF = LSAME( JOBU, 'F' )
+ LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA
+ LSVEC = LSVC0 .OR. WNTUF
+ DNTWU = LSAME( JOBU, 'N' )
+*
+ WNTVR = LSAME( JOBV, 'R' )
+ WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' )
+ RSVEC = WNTVR .OR. WNTVA
+ DNTWV = LSAME( JOBV, 'N' )
+*
+ ACCLA = LSAME( JOBA, 'A' )
+ ACCLM = LSAME( JOBA, 'M' )
+ CONDA = LSAME( JOBA, 'E' )
+ ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA
+*
+ ROWPRM = LSAME( JOBP, 'P' )
+ RTRANS = LSAME( JOBR, 'T' )
+*
+ IF ( ROWPRM ) THEN
+ IF ( CONDA ) THEN
+ IMINWRK = MAX( 1, N + M - 1 + N )
+ ELSE
+ IMINWRK = MAX( 1, N + M - 1 )
+ END IF
+ RMINWRK = MAX( 2, M )
+ ELSE
+ IF ( CONDA ) THEN
+ IMINWRK = MAX( 1, N + N )
+ ELSE
+ IMINWRK = MAX( 1, N )
+ END IF
+ RMINWRK = 2
+ END IF
+ LQUERY = (LIWORK .EQ. -1 .OR. LWORK .EQ. -1 .OR. LRWORK .EQ. -1)
+ INFO = 0
+ IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN
+ INFO = -1
+ ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN
+ INFO = -2
+ ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN
+ INFO = -3
+ ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN
+ INFO = -4
+ ELSE IF ( WNTUR .AND. WNTVA ) THEN
+ INFO = -5
+ ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN
+ INFO = -5
+ ELSE IF ( M.LT.0 ) THEN
+ INFO = -6
+ ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
+ INFO = -7
+ ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN
+ INFO = -9
+ ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR.
+ $ ( WNTUF .AND. LDU.LT.N ) ) THEN
+ INFO = -12
+ ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR.
+ $ ( CONDA .AND. LDV.LT.N ) ) THEN
+ INFO = -14
+ ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN
+ INFO = -17
+ END IF
+*
+*
+ IF ( INFO .EQ. 0 ) THEN
+* .. compute the minimal and the optimal workspace lengths
+* [[The expressions for computing the minimal and the optimal
+* values of LWORK are written with a lot of redundancy and
+* can be simplified. However, this detailed form is easier for
+* maintenance and modifications of the code.]]
+*
+* .. minimal workspace length for DGEQP3 of an M x N matrix
+ LWQP3 = 3 * N + 1
+* .. minimal workspace length for DORMQR to build left singular vectors
+ IF ( WNTUS .OR. WNTUR ) THEN
+ LWORQ = MAX( N , 1 )
+ ELSE IF ( WNTUA ) THEN
+ LWORQ = MAX( M , 1 )
+ END IF
+* .. minimal workspace length for DPOCON of an N x N matrix
+ LWCON = 3 * N
+* .. DGESVD of an N x N matrix
+ LWSVD = MAX( 5 * N, 1 )
+ IF ( LQUERY ) THEN
+ CALL DGEQP3( M, N, A, LDA, IWORK, RDUMMY, RDUMMY, -1,
+ $ IERR )
+ LWRK_DGEQP3 = INT( RDUMMY(1) )
+ IF ( WNTUS .OR. WNTUR ) THEN
+ CALL DORMQR( 'L', 'N', M, N, N, A, LDA, RDUMMY, U,
+ $ LDU, RDUMMY, -1, IERR )
+ LWRK_DORMQR = INT( RDUMMY(1) )
+ ELSE IF ( WNTUA ) THEN
+ CALL DORMQR( 'L', 'N', M, M, N, A, LDA, RDUMMY, U,
+ $ LDU, RDUMMY, -1, IERR )
+ LWRK_DORMQR = INT( RDUMMY(1) )
+ ELSE
+ LWRK_DORMQR = 0
+ END IF
+ END IF
+ MINWRK = 2
+ OPTWRK = 2
+ IF ( .NOT. (LSVEC .OR. RSVEC )) THEN
+* .. minimal and optimal sizes of the workspace if
+* only the singular values are requested
+ IF ( CONDA ) THEN
+ MINWRK = MAX( N+LWQP3, LWCON, LWSVD )
+ ELSE
+ MINWRK = MAX( N+LWQP3, LWSVD )
+ END IF
+ IF ( LQUERY ) THEN
+ CALL DGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = MAX( N+LWRK_DGEQP3, N+LWCON, LWRK_DGESVD )
+ ELSE
+ OPTWRK = MAX( N+LWRK_DGEQP3, LWRK_DGESVD )
+ END IF
+ END IF
+ ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
+* .. minimal and optimal sizes of the workspace if the
+* singular values and the left singular vectors are requested
+ IF ( CONDA ) THEN
+ MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWORQ )
+ ELSE
+ MINWRK = N + MAX( LWQP3, LWSVD, LWORQ )
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL DGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ ELSE
+ CALL DGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ END IF
+ LWRK_DGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = N + MAX( LWRK_DGEQP3, LWCON, LWRK_DGESVD,
+ $ LWRK_DORMQR )
+ ELSE
+ OPTWRK = N + MAX( LWRK_DGEQP3, LWRK_DGESVD,
+ $ LWRK_DORMQR )
+ END IF
+ END IF
+ ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
+* .. minimal and optimal sizes of the workspace if the
+* singular values and the right singular vectors are requested
+ IF ( CONDA ) THEN
+ MINWRK = N + MAX( LWQP3, LWCON, LWSVD )
+ ELSE
+ MINWRK = N + MAX( LWQP3, LWSVD )
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL DGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ ELSE
+ CALL DGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ END IF
+ LWRK_DGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = N + MAX( LWRK_DGEQP3, LWCON, LWRK_DGESVD )
+ ELSE
+ OPTWRK = N + MAX( LWRK_DGEQP3, LWRK_DGESVD )
+ END IF
+ END IF
+ ELSE
+* .. minimal and optimal sizes of the workspace if the
+* full SVD is requested
+ IF ( RTRANS ) THEN
+ MINWRK = MAX( LWQP3, LWSVD, LWORQ )
+ IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
+ MINWRK = MINWRK + N
+ IF ( WNTVA ) THEN
+* .. minimal workspace length for N x N/2 DGEQRF
+ LWQRF = MAX( N/2, 1 )
+* .. minimal workspace lengt for N/2 x N/2 DGESVD
+ LWSVD2 = MAX( 5 * (N/2), 1 )
+ LWORQ2 = MAX( N, 1 )
+ MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2,
+ $ N/2+LWORQ2, LWORQ )
+ IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
+ MINWRK2 = N + MINWRK2
+ MINWRK = MAX( MINWRK, MINWRK2 )
+ END IF
+ ELSE
+ MINWRK = MAX( LWQP3, LWSVD, LWORQ )
+ IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
+ MINWRK = MINWRK + N
+ IF ( WNTVA ) THEN
+* .. minimal workspace length for N/2 x N DGELQF
+ LWLQF = MAX( N/2, 1 )
+ LWSVD2 = MAX( 5 * (N/2), 1 )
+ LWORLQ = MAX( N , 1 )
+ MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2,
+ $ N/2+LWORLQ, LWORQ )
+ IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
+ MINWRK2 = N + MINWRK2
+ MINWRK = MAX( MINWRK, MINWRK2 )
+ END IF
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL DGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DGESVD = INT( RDUMMY(1) )
+ OPTWRK = MAX(LWRK_DGEQP3,LWRK_DGESVD,LWRK_DORMQR)
+ IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
+ OPTWRK = N + OPTWRK
+ IF ( WNTVA ) THEN
+ CALL DGEQRF(N,N/2,U,LDU,RDUMMY,RDUMMY,-1,IERR)
+ LWRK_DGEQRF = INT( RDUMMY(1) )
+ CALL DGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DGESVD2 = INT( RDUMMY(1) )
+ CALL DORMQR( 'R', 'C', N, N, N/2, U, LDU, RDUMMY,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DORMQR2 = INT( RDUMMY(1) )
+ OPTWRK2 = MAX( LWRK_DGEQP3, N/2+LWRK_DGEQRF,
+ $ N/2+LWRK_DGESVD2, N/2+LWRK_DORMQR2 )
+ IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
+ OPTWRK2 = N + OPTWRK2
+ OPTWRK = MAX( OPTWRK, OPTWRK2 )
+ END IF
+ ELSE
+ CALL DGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DGESVD = INT( RDUMMY(1) )
+ OPTWRK = MAX(LWRK_DGEQP3,LWRK_DGESVD,LWRK_DORMQR)
+ IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
+ OPTWRK = N + OPTWRK
+ IF ( WNTVA ) THEN
+ CALL DGELQF(N/2,N,U,LDU,RDUMMY,RDUMMY,-1,IERR)
+ LWRK_DGELQF = INT( RDUMMY(1) )
+ CALL DGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_DGESVD2 = INT( RDUMMY(1) )
+ CALL DORMLQ( 'R', 'N', N, N, N/2, U, LDU, RDUMMY,
+ $ V, LDV, RDUMMY,-1,IERR )
+ LWRK_DORMLQ = INT( RDUMMY(1) )
+ OPTWRK2 = MAX( LWRK_DGEQP3, N/2+LWRK_DGELQF,
+ $ N/2+LWRK_DGESVD2, N/2+LWRK_DORMLQ )
+ IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
+ OPTWRK2 = N + OPTWRK2
+ OPTWRK = MAX( OPTWRK, OPTWRK2 )
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ MINWRK = MAX( 2, MINWRK )
+ OPTWRK = MAX( 2, OPTWRK )
+ IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19
+*
+ END IF
+*
+ IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN
+ INFO = -21
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DGESVDQ', -INFO )
+ RETURN
+ ELSE IF ( LQUERY ) THEN
+*
+* Return optimal workspace
+*
+ IWORK(1) = IMINWRK
+ WORK(1) = OPTWRK
+ WORK(2) = MINWRK
+ RWORK(1) = RMINWRK
+ RETURN
+ END IF
+*
+* Quick return if the matrix is void.
+*
+ IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN
+* .. all output is void.
+ RETURN
+ END IF
+*
+ BIG = DLAMCH('O')
+ ASCALED = .FALSE.
+ IWOFF = 1
+ IF ( ROWPRM ) THEN
+ IWOFF = M
+* .. reordering the rows in decreasing sequence in the
+* ell-infinity norm - this enhances numerical robustness in
+* the case of differently scaled rows.
+ DO 1904 p = 1, M
+* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
+* [[DLANGE will return NaN if an entry of the p-th row is Nan]]
+ RWORK(p) = DLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY )
+* .. check for NaN's and Inf's
+ IF ( ( RWORK(p) .NE. RWORK(p) ) .OR.
+ $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN
+ INFO = -8
+ CALL XERBLA( 'DGESVDQ', -INFO )
+ RETURN
+ END IF
+ 1904 CONTINUE
+ DO 1952 p = 1, M - 1
+ q = IDAMAX( M-p+1, RWORK(p), 1 ) + p - 1
+ IWORK(N+p) = q
+ IF ( p .NE. q ) THEN
+ RTMP = RWORK(p)
+ RWORK(p) = RWORK(q)
+ RWORK(q) = RTMP
+ END IF
+ 1952 CONTINUE
+*
+ IF ( RWORK(1) .EQ. ZERO ) THEN
+* Quick return: A is the M x N zero matrix.
+ NUMRANK = 0
+ CALL DLASET( 'G', N, 1, ZERO, ZERO, S, N )
+ IF ( WNTUS ) CALL DLASET('G', M, N, ZERO, ONE, U, LDU)
+ IF ( WNTUA ) CALL DLASET('G', M, M, ZERO, ONE, U, LDU)
+ IF ( WNTVA ) CALL DLASET('G', N, N, ZERO, ONE, V, LDV)
+ IF ( WNTUF ) THEN
+ CALL DLASET( 'G', N, 1, ZERO, ZERO, WORK, N )
+ CALL DLASET( 'G', M, N, ZERO, ONE, U, LDU )
+ END IF
+ DO 5001 p = 1, N
+ IWORK(p) = p
+ 5001 CONTINUE
+ IF ( ROWPRM ) THEN
+ DO 5002 p = N + 1, N + M - 1
+ IWORK(p) = p - N
+ 5002 CONTINUE
+ END IF
+ IF ( CONDA ) RWORK(1) = -1
+ RWORK(2) = -1
+ RETURN
+ END IF
+*
+ IF ( RWORK(1) .GT. BIG / SQRT(DBLE(M)) ) THEN
+* .. to prevent overflow in the QR factorization, scale the
+* matrix by 1/sqrt(M) if too large entry detected
+ CALL DLASCL('G',0,0,SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
+ ASCALED = .TRUE.
+ END IF
+ CALL DLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 )
+ END IF
+*
+* .. At this stage, preemptive scaling is done only to avoid column
+* norms overflows during the QR factorization. The SVD procedure should
+* have its own scaling to save the singular values from overflows and
+* underflows. That depends on the SVD procedure.
+*
+ IF ( .NOT.ROWPRM ) THEN
+ RTMP = DLANGE( 'M', M, N, A, LDA, RDUMMY )
+ IF ( ( RTMP .NE. RTMP ) .OR.
+ $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN
+ INFO = -8
+ CALL XERBLA( 'DGESVDQ', -INFO )
+ RETURN
+ END IF
+ IF ( RTMP .GT. BIG / SQRT(DBLE(M)) ) THEN
+* .. to prevent overflow in the QR factorization, scale the
+* matrix by 1/sqrt(M) if too large entry detected
+ CALL DLASCL('G',0,0, SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
+ ASCALED = .TRUE.
+ END IF
+ END IF
+*
+* .. QR factorization with column pivoting
+*
+* A * P = Q * [ R ]
+* [ 0 ]
+*
+ DO 1963 p = 1, N
+* .. all columns are free columns
+ IWORK(p) = 0
+ 1963 CONTINUE
+ CALL DGEQP3( M, N, A, LDA, IWORK, WORK, WORK(N+1), LWORK-N,
+ $ IERR )
+*
+* If the user requested accuracy level allows truncation in the
+* computed upper triangular factor, the matrix R is examined and,
+* if possible, replaced with its leading upper trapezoidal part.
+*
+ EPSLN = DLAMCH('E')
+ SFMIN = DLAMCH('S')
+* SMALL = SFMIN / EPSLN
+ NR = N
+*
+ IF ( ACCLA ) THEN
+*
+* Standard absolute error bound suffices. All sigma_i with
+* sigma_i < N*EPS*||A||_F are flushed to zero. This is an
+* aggressive enforcement of lower numerical rank by introducing a
+* backward error of the order of N*EPS*||A||_F.
+ NR = 1
+ RTMP = SQRT(DBLE(N))*EPSLN
+ DO 3001 p = 2, N
+ IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002
+ NR = NR + 1
+ 3001 CONTINUE
+ 3002 CONTINUE
+*
+ ELSEIF ( ACCLM ) THEN
+* .. similarly as above, only slightly more gentle (less aggressive).
+* Sudden drop on the diagonal of R is used as the criterion for being
+* close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
+* [[This can be made more flexible by replacing this hard-coded value
+* with a user specified threshold.]] Also, the values that underflow
+* will be truncated.
+ NR = 1
+ DO 3401 p = 2, N
+ IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
+ $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402
+ NR = NR + 1
+ 3401 CONTINUE
+ 3402 CONTINUE
+*
+ ELSE
+* .. RRQR not authorized to determine numerical rank except in the
+* obvious case of zero pivots.
+* .. inspect R for exact zeros on the diagonal;
+* R(i,i)=0 => R(i:N,i:N)=0.
+ NR = 1
+ DO 3501 p = 2, N
+ IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502
+ NR = NR + 1
+ 3501 CONTINUE
+ 3502 CONTINUE
+*
+ IF ( CONDA ) THEN
+* Estimate the scaled condition number of A. Use the fact that it is
+* the same as the scaled condition number of R.
+* .. V is used as workspace
+ CALL DLACPY( 'U', N, N, A, LDA, V, LDV )
+* Only the leading NR x NR submatrix of the triangular factor
+* is considered. Only if NR=N will this give a reliable error
+* bound. However, even for NR < N, this can be used on an
+* expert level and obtain useful information in the sense of
+* perturbation theory.
+ DO 3053 p = 1, NR
+ RTMP = DNRM2( p, V(1,p), 1 )
+ CALL DSCAL( p, ONE/RTMP, V(1,p), 1 )
+ 3053 CONTINUE
+ IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN
+ CALL DPOCON( 'U', NR, V, LDV, ONE, RTMP,
+ $ WORK, IWORK(N+IWOFF), IERR )
+ ELSE
+ CALL DPOCON( 'U', NR, V, LDV, ONE, RTMP,
+ $ WORK(N+1), IWORK(N+IWOFF), IERR )
+ END IF
+ SCONDA = ONE / SQRT(RTMP)
+* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
+* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
+* See the reference [1] for more details.
+ END IF
+*
+ ENDIF
+*
+ IF ( WNTUR ) THEN
+ N1 = NR
+ ELSE IF ( WNTUS .OR. WNTUF) THEN
+ N1 = N
+ ELSE IF ( WNTUA ) THEN
+ N1 = M
+ END IF
+*
+ IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
+*.......................................................................
+* .. only the singular values are requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+*
+* .. compute the singular values of R**T = [A](1:NR,1:N)**T
+* .. set the lower triangle of [A] to [A](1:NR,1:N)**T and
+* the upper triangle of [A] to zero.
+ DO 1146 p = 1, MIN( N, NR )
+ DO 1147 q = p + 1, N
+ A(q,p) = A(p,q)
+ IF ( q .LE. NR ) A(p,q) = ZERO
+ 1147 CONTINUE
+ 1146 CONTINUE
+*
+ CALL DGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU,
+ $ V, LDV, WORK, LWORK, INFO )
+*
+ ELSE
+*
+* .. compute the singular values of R = [A](1:NR,1:N)
+*
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'L', NR-1,NR-1, ZERO,ZERO, A(2,1), LDA )
+ CALL DGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU,
+ $ V, LDV, WORK, LWORK, INFO )
+*
+ END IF
+*
+ ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN
+*.......................................................................
+* .. the singular values and the left singular vectors requested
+*.......................................................................""""""""
+ IF ( RTRANS ) THEN
+* .. apply DGESVD to R**T
+* .. copy R**T into [U] and overwrite [U] with the right singular
+* vectors of R
+ DO 1192 p = 1, NR
+ DO 1193 q = p, N
+ U(q,p) = A(p,q)
+ 1193 CONTINUE
+ 1192 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, U(1,2), LDU )
+* .. the left singular vectors not computed, the NR right singular
+* vectors overwrite [U](1:NR,1:NR) as transposed. These
+* will be pre-multiplied by Q to build the left singular vectors of A.
+ CALL DGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1119 p = 1, NR
+ DO 1120 q = p + 1, NR
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1120 CONTINUE
+ 1119 CONTINUE
+*
+ ELSE
+* .. apply DGESVD to R
+* .. copy R into [U] and overwrite [U] with the left singular vectors
+ CALL DLACPY( 'U', NR, N, A, LDA, U, LDU )
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'L', NR-1, NR-1, ZERO, ZERO, U(2,1), LDU )
+* .. the right singular vectors not computed, the NR left singular
+* vectors overwrite [U](1:NR,1:NR)
+ CALL DGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+* .. now [U](1:NR,1:NR) contains the NR left singular vectors of
+* R. These will be pre-multiplied by Q to build the left singular
+* vectors of A.
+ END IF
+*
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN
+ CALL DLASET('A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL DLASET( 'A',NR,N1-NR,ZERO,ZERO,U(1,NR+1), LDU )
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+* The Q matrix from the first QRF is built into the left singular
+* vectors matrix U.
+*
+ IF ( .NOT.WNTUF )
+ $ CALL DORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U,
+ $ LDU, WORK(N+1), LWORK-N, IERR )
+ IF ( ROWPRM .AND. .NOT.WNTUF )
+ $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
+*
+ ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
+*.......................................................................
+* .. the singular values and the right singular vectors requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+* .. apply DGESVD to R**T
+* .. copy R**T into V and overwrite V with the left singular vectors
+ DO 1165 p = 1, NR
+ DO 1166 q = p, N
+ V(q,p) = (A(p,q))
+ 1166 CONTINUE
+ 1165 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
+* .. the left singular vectors of R**T overwrite V, the right singular
+* vectors not computed
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+ CALL DGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1121 p = 1, NR
+ DO 1122 q = p + 1, NR
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1122 CONTINUE
+ 1121 CONTINUE
+*
+ IF ( NR .LT. N ) THEN
+ DO 1103 p = 1, NR
+ DO 1104 q = NR + 1, N
+ V(p,q) = V(q,p)
+ 1104 CONTINUE
+ 1103 CONTINUE
+ END IF
+ CALL DLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+ ELSE
+* .. need all N right singular vectors and NR < N
+* [!] This is simple implementation that augments [V](1:N,1:NR)
+* by padding a zero block. In the case NR << N, a more efficient
+* way is to first use the QR factorization. For more details
+* how to implement this, see the " FULL SVD " branch.
+ CALL DLASET('G', N, N-NR, ZERO, ZERO, V(1,NR+1), LDV)
+ CALL DGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1123 p = 1, N
+ DO 1124 q = p + 1, N
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1124 CONTINUE
+ 1123 CONTINUE
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+ END IF
+*
+ ELSE
+* .. aply DGESVD to R
+* .. copy R into V and overwrite V with the right singular vectors
+ CALL DLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'L', NR-1, NR-1, ZERO, ZERO, V(2,1), LDV )
+* .. the right singular vectors overwrite V, the NR left singular
+* vectors stored in U(1:NR,1:NR)
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+ CALL DGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL DLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
+ ELSE
+* .. need all N right singular vectors and NR < N
+* [!] This is simple implementation that augments [V](1:NR,1:N)
+* by padding a zero block. In the case NR << N, a more efficient
+* way is to first use the LQ factorization. For more details
+* how to implement this, see the " FULL SVD " branch.
+ CALL DLASET('G', N-NR, N, ZERO,ZERO, V(NR+1,1), LDV)
+ CALL DGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+ END IF
+* .. now [V] contains the transposed matrix of the right singular
+* vectors of A.
+ END IF
+*
+ ELSE
+*.......................................................................
+* .. FULL SVD requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+*
+* .. apply DGESVD to R**T [[this option is left for R&D&T]]
+*
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+* .. copy R**T into [V] and overwrite [V] with the left singular
+* vectors of R**T
+ DO 1168 p = 1, NR
+ DO 1169 q = p, N
+ V(q,p) = A(p,q)
+ 1169 CONTINUE
+ 1168 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
+*
+* .. the left singular vectors of R**T overwrite [V], the NR right
+* singular vectors of R**T stored in [U](1:NR,1:NR) as transposed
+ CALL DGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+* .. assemble V
+ DO 1115 p = 1, NR
+ DO 1116 q = p + 1, NR
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1116 CONTINUE
+ 1115 CONTINUE
+ IF ( NR .LT. N ) THEN
+ DO 1101 p = 1, NR
+ DO 1102 q = NR+1, N
+ V(p,q) = V(q,p)
+ 1102 CONTINUE
+ 1101 CONTINUE
+ END IF
+ CALL DLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+*
+ DO 1117 p = 1, NR
+ DO 1118 q = p + 1, NR
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1118 CONTINUE
+ 1117 CONTINUE
+*
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+ ELSE
+* .. need all N right singular vectors and NR < N
+* .. copy R**T into [V] and overwrite [V] with the left singular
+* vectors of R**T
+* [[The optimal ratio N/NR for using QRF instead of padding
+* with zeros. Here hard coded to 2; it must be at least
+* two due to work space constraints.]]
+* OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
+* OPTRATIO = MAX( OPTRATIO, 2 )
+ OPTRATIO = 2
+ IF ( OPTRATIO*NR .GT. N ) THEN
+ DO 1198 p = 1, NR
+ DO 1199 q = p, N
+ V(q,p) = A(p,q)
+ 1199 CONTINUE
+ 1198 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL DLASET('U',NR-1,NR-1, ZERO,ZERO, V(1,2),LDV)
+*
+ CALL DLASET('A',N,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL DGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1113 p = 1, N
+ DO 1114 q = p + 1, N
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1114 CONTINUE
+ 1113 CONTINUE
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x N1), i.e. (M x N) or (M x M).
+*
+ DO 1111 p = 1, N
+ DO 1112 q = p + 1, N
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1112 CONTINUE
+ 1111 CONTINUE
+*
+ IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU)
+ IF ( N .LT. N1 ) THEN
+ CALL DLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU)
+ CALL DLASET('A',M-N,N1-N,ZERO,ONE,
+ $ U(N+1,N+1), LDU )
+ END IF
+ END IF
+ ELSE
+* .. copy R**T into [U] and overwrite [U] with the right
+* singular vectors of R
+ DO 1196 p = 1, NR
+ DO 1197 q = p, N
+ U(q,NR+p) = A(p,q)
+ 1197 CONTINUE
+ 1196 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL DLASET('U',NR-1,NR-1,ZERO,ZERO,U(1,NR+2),LDU)
+ CALL DGEQRF( N, NR, U(1,NR+1), LDU, WORK(N+1),
+ $ WORK(N+NR+1), LWORK-N-NR, IERR )
+ DO 1143 p = 1, NR
+ DO 1144 q = 1, N
+ V(q,p) = U(p,NR+q)
+ 1144 CONTINUE
+ 1143 CONTINUE
+ CALL DLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV)
+ CALL DGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
+ $ V,LDV, WORK(N+NR+1),LWORK-N-NR, INFO )
+ CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
+ CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
+ CALL DORMQR('R','C', N, N, NR, U(1,NR+1), LDU,
+ $ WORK(N+1),V,LDV,WORK(N+NR+1),LWORK-N-NR,IERR)
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1),LDU)
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ ELSE
+*
+* .. apply DGESVD to R [[this is the recommended option]]
+*
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+* .. copy R into [V] and overwrite V with the right singular vectors
+ CALL DLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL DLASET( 'L', NR-1,NR-1, ZERO,ZERO, V(2,1), LDV )
+* .. the right singular vectors of R overwrite [V], the NR left
+* singular vectors of R stored in [U](1:NR,1:NR)
+ CALL DGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL DLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+ ELSE
+* .. need all N right singular vectors and NR < N
+* .. the requested number of the left singular vectors
+* is then N1 (N or M)
+* [[The optimal ratio N/NR for using LQ instead of padding
+* with zeros. Here hard coded to 2; it must be at least
+* two due to work space constraints.]]
+* OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
+* OPTRATIO = MAX( OPTRATIO, 2 )
+ OPTRATIO = 2
+ IF ( OPTRATIO * NR .GT. N ) THEN
+ CALL DLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL DLASET('L', NR-1,NR-1, ZERO,ZERO, V(2,1),LDV)
+* .. the right singular vectors of R overwrite [V], the NR left
+* singular vectors of R stored in [U](1:NR,1:NR)
+ CALL DLASET('A', N-NR,N, ZERO,ZERO, V(NR+1,1),LDV)
+ CALL DGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. now [V] contains the transposed matrix of the right
+* singular vectors of A. The leading N left singular vectors
+* are in [U](1:N,1:N)
+* .. assemble the left singular vector matrix U of dimensions
+* (M x N1), i.e. (M x N) or (M x M).
+ IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU)
+ IF ( N .LT. N1 ) THEN
+ CALL DLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU)
+ CALL DLASET( 'A',M-N,N1-N,ZERO,ONE,
+ $ U(N+1,N+1), LDU )
+ END IF
+ END IF
+ ELSE
+ CALL DLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU )
+ IF ( NR .GT. 1 )
+ $ CALL DLASET('L',NR-1,NR-1,ZERO,ZERO,U(NR+2,1),LDU)
+ CALL DGELQF( NR, N, U(NR+1,1), LDU, WORK(N+1),
+ $ WORK(N+NR+1), LWORK-N-NR, IERR )
+ CALL DLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV)
+ IF ( NR .GT. 1 )
+ $ CALL DLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV)
+ CALL DGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+NR+1), LWORK-N-NR, INFO )
+ CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
+ CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
+ CALL DORMLQ('R','N',N,N,NR,U(NR+1,1),LDU,WORK(N+1),
+ $ V, LDV, WORK(N+NR+1),LWORK-N-NR,IERR)
+ CALL DLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL DLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+ END IF
+ END IF
+* .. end of the "R**T or R" branch
+ END IF
+*
+* The Q matrix from the first QRF is built into the left singular
+* vectors matrix U.
+*
+ IF ( .NOT. WNTUF )
+ $ CALL DORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U,
+ $ LDU, WORK(N+1), LWORK-N, IERR )
+ IF ( ROWPRM .AND. .NOT.WNTUF )
+ $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
+*
+* ... end of the "full SVD" branch
+ END IF
+*
+* Check whether some singular values are returned as zeros, e.g.
+* due to underflow, and update the numerical rank.
+ p = NR
+ DO 4001 q = p, 1, -1
+ IF ( S(q) .GT. ZERO ) GO TO 4002
+ NR = NR - 1
+ 4001 CONTINUE
+ 4002 CONTINUE
+*
+* .. if numerical rank deficiency is detected, the truncated
+* singular values are set to zero.
+ IF ( NR .LT. N ) CALL DLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N )
+* .. undo scaling; this may cause overflow in the largest singular
+* values.
+ IF ( ASCALED )
+ $ CALL DLASCL( 'G',0,0, ONE,SQRT(DBLE(M)), NR,1, S, N, IERR )
+ IF ( CONDA ) RWORK(1) = SCONDA
+ RWORK(2) = p - NR
+* .. p-NR is the number of singular values that are computed as
+* exact zeros in DGESVD() applied to the (possibly truncated)
+* full row rank triangular (trapezoidal) factor of A.
+ NUMRANK = NR
+*
+ RETURN
+*
+* End of DGESVDQ
+*
+ END
diff --git a/lapack-netlib/SRC/dgesvj.f b/lapack-netlib/SRC/dgesvj.f
index 2cbc5ce0e..cf7aac982 100644
--- a/lapack-netlib/SRC/dgesvj.f
+++ b/lapack-netlib/SRC/dgesvj.f
@@ -90,13 +90,13 @@
*> JOBV is CHARACTER*1
*> Specifies whether to compute the right singular vectors, that
*> is, the matrix V:
-*> = 'V' : the matrix V is computed and returned in the array V
-*> = 'A' : the Jacobi rotations are applied to the MV-by-N
+*> = 'V': the matrix V is computed and returned in the array V
+*> = 'A': the Jacobi rotations are applied to the MV-by-N
*> array V. In other words, the right singular vector
*> matrix V is not computed explicitly, instead it is
*> applied to an MV-by-N matrix initially stored in the
*> first MV rows of V.
-*> = 'N' : the matrix V is not computed and the array V is not
+*> = 'N': the matrix V is not computed and the array V is not
*> referenced
*> \endverbatim
*>
@@ -118,8 +118,8 @@
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> On entry, the M-by-N matrix A.
*> On exit :
-*> If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
-*> If INFO .EQ. 0 :
+*> If JOBU = 'U' .OR. JOBU = 'C' :
+*> If INFO = 0 :
*> RANKA orthonormal columns of U are returned in the
*> leading RANKA columns of the array A. Here RANKA <= N
*> is the number of computed singular values of A that are
@@ -129,9 +129,9 @@
*> in the array WORK as RANKA=NINT(WORK(2)). Also see the
*> descriptions of SVA and WORK. The computed columns of U
*> are mutually numerically orthogonal up to approximately
-*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
+*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
*> see the description of JOBU.
-*> If INFO .GT. 0 :
+*> If INFO > 0 :
*> the procedure DGESVJ did not converge in the given number
*> of iterations (sweeps). In that case, the computed
*> columns of U may not be orthogonal up to TOL. The output
@@ -140,8 +140,8 @@
*> input matrix A in the sense that the residual
*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
*>
-*> If JOBU .EQ. 'N' :
-*> If INFO .EQ. 0 :
+*> If JOBU = 'N' :
+*> If INFO = 0 :
*> Note that the left singular vectors are 'for free' in the
*> one-sided Jacobi SVD algorithm. However, if only the
*> singular values are needed, the level of numerical
@@ -150,7 +150,7 @@
*> numerically orthogonal up to approximately M*EPS. Thus,
*> on exit, A contains the columns of U scaled with the
*> corresponding singular values.
-*> If INFO .GT. 0 :
+*> If INFO > 0 :
*> the procedure DGESVJ did not converge in the given number
*> of iterations (sweeps).
*> \endverbatim
@@ -165,9 +165,9 @@
*> \verbatim
*> SVA is DOUBLE PRECISION array, dimension (N)
*> On exit :
-*> If INFO .EQ. 0 :
+*> If INFO = 0 :
*> depending on the value SCALE = WORK(1), we have:
-*> If SCALE .EQ. ONE :
+*> If SCALE = ONE :
*> SVA(1:N) contains the computed singular values of A.
*> During the computation SVA contains the Euclidean column
*> norms of the iterated matrices in the array A.
@@ -175,7 +175,7 @@
*> The singular values of A are SCALE*SVA(1:N), and this
*> factored representation is due to the fact that some of the
*> singular values of A might underflow or overflow.
-*> If INFO .GT. 0 :
+*> If INFO > 0 :
*> the procedure DGESVJ did not converge in the given number of
*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
*> \endverbatim
@@ -183,7 +183,7 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ
+*> If JOBV = 'A', then the product of Jacobi rotations in DGESVJ
*> is applied to the first MV rows of V. See the description of JOBV.
*> \endverbatim
*>
@@ -201,16 +201,16 @@
*> \param[in] LDV
*> \verbatim
*> LDV is INTEGER
-*> The leading dimension of the array V, LDV .GE. 1.
-*> If JOBV .EQ. 'V', then LDV .GE. max(1,N).
-*> If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
+*> The leading dimension of the array V, LDV >= 1.
+*> If JOBV = 'V', then LDV >= max(1,N).
+*> If JOBV = 'A', then LDV >= max(1,MV) .
*> \endverbatim
*>
*> \param[in,out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (LWORK)
*> On entry :
-*> If JOBU .EQ. 'C' :
+*> If JOBU = 'C' :
*> WORK(1) = CTOL, where CTOL defines the threshold for convergence.
*> The process stops if all columns of A are mutually
*> orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
@@ -230,7 +230,7 @@
*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
*> This is useful information in cases when DGESVJ did
*> not converge, as it can be used to estimate whether
-*> the output is stil useful and for post festum analysis.
+*> the output is still useful and for post festum analysis.
*> WORK(6) = the largest absolute value over all sines of the
*> Jacobi rotation angles in the last sweep. It can be
*> useful for a post festum analysis.
@@ -245,9 +245,9 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
-*> > 0 : DGESVJ did not converge in the maximal allowed number (30)
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
+*> > 0: DGESVJ did not converge in the maximal allowed number (30)
*> of sweeps. The output may still be useful. See the
*> description of WORK.
*> \endverbatim
diff --git a/lapack-netlib/SRC/dgesvxx.f b/lapack-netlib/SRC/dgesvxx.f
index afcd05d8e..21b56f61c 100644
--- a/lapack-netlib/SRC/dgesvxx.f
+++ b/lapack-netlib/SRC/dgesvxx.f
@@ -411,7 +411,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -447,14 +447,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -462,9 +462,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.0D+0
-*> = 0.0 : No refinement is performed, and no error bounds are
+*> = 0.0: No refinement is performed, and no error bounds are
*> computed.
-*> = 1.0 : Use the extra-precise refinement algorithm.
+*> = 1.0: Use the extra-precise refinement algorithm.
*> (other values are reserved for future use)
*>
*> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
diff --git a/lapack-netlib/SRC/dgetc2.f b/lapack-netlib/SRC/dgetc2.f
index 0896a7013..5bf5b890f 100644
--- a/lapack-netlib/SRC/dgetc2.f
+++ b/lapack-netlib/SRC/dgetc2.f
@@ -85,7 +85,7 @@
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
-*> > 0: if INFO = k, U(k, k) is likely to produce owerflow if
+*> > 0: if INFO = k, U(k, k) is likely to produce overflow if
*> we try to solve for x in Ax = b. So U is perturbed to
*> avoid the overflow.
*> \endverbatim
diff --git a/lapack-netlib/SRC/dgetsls.f b/lapack-netlib/SRC/dgetsls.f
index 3b44a40ab..dfc72c8b2 100644
--- a/lapack-netlib/SRC/dgetsls.f
+++ b/lapack-netlib/SRC/dgetsls.f
@@ -1,3 +1,5 @@
+*> \brief \b DGETSLS
+*
* Definition:
* ===========
*
diff --git a/lapack-netlib/SRC/dggesx.f b/lapack-netlib/SRC/dggesx.f
index 47022fbdf..0e57d636e 100644
--- a/lapack-netlib/SRC/dggesx.f
+++ b/lapack-netlib/SRC/dggesx.f
@@ -131,10 +131,10 @@
*> \verbatim
*> SENSE is CHARACTER*1
*> Determines which reciprocal condition numbers are computed.
-*> = 'N' : None are computed;
-*> = 'E' : Computed for average of selected eigenvalues only;
-*> = 'V' : Computed for selected deflating subspaces only;
-*> = 'B' : Computed for both.
+*> = 'N': None are computed;
+*> = 'E': Computed for average of selected eigenvalues only;
+*> = 'V': Computed for selected deflating subspaces only;
+*> = 'B': Computed for both.
*> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
*> \endverbatim
*>
diff --git a/lapack-netlib/SRC/dgsvj0.f b/lapack-netlib/SRC/dgsvj0.f
index 4fd38d37e..318e369fd 100644
--- a/lapack-netlib/SRC/dgsvj0.f
+++ b/lapack-netlib/SRC/dgsvj0.f
@@ -117,7 +117,7 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
+*> If JOBV = 'A', then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then MV is not referenced.
*> \endverbatim
@@ -125,9 +125,9 @@
*> \param[in,out] V
*> \verbatim
*> V is DOUBLE PRECISION array, dimension (LDV,N)
-*> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
+*> If JOBV = 'V' then N rows of V are post-multipled by a
*> sequence of Jacobi rotations.
-*> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
+*> If JOBV = 'A' then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then V is not referenced.
*> \endverbatim
@@ -136,8 +136,8 @@
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of the array V, LDV >= 1.
-*> If JOBV = 'V', LDV .GE. N.
-*> If JOBV = 'A', LDV .GE. MV.
+*> If JOBV = 'V', LDV >= N.
+*> If JOBV = 'A', LDV >= MV.
*> \endverbatim
*>
*> \param[in] EPS
@@ -157,7 +157,7 @@
*> TOL is DOUBLE PRECISION
*> TOL is the threshold for Jacobi rotations. For a pair
*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
-*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
+*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) > TOL.
*> \endverbatim
*>
*> \param[in] NSWEEP
@@ -175,14 +175,14 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> LWORK is the dimension of WORK. LWORK .GE. M.
+*> LWORK is the dimension of WORK. LWORK >= M.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
@@ -1045,7 +1045,7 @@
1993 CONTINUE
* end i=1:NSWEEP loop
-* #:) Reaching this point means that the procedure has comleted the given
+* #:) Reaching this point means that the procedure has completed the given
* number of iterations.
INFO = NSWEEP - 1
GO TO 1995
diff --git a/lapack-netlib/SRC/dgsvj1.f b/lapack-netlib/SRC/dgsvj1.f
index 376682c7f..35a93619f 100644
--- a/lapack-netlib/SRC/dgsvj1.f
+++ b/lapack-netlib/SRC/dgsvj1.f
@@ -61,7 +61,7 @@
*> In terms of the columns of A, the first N1 columns are rotated 'against'
*> the remaining N-N1 columns, trying to increase the angle between the
*> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
-*> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
+*> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parameter.
*> The number of sweeps is given in NSWEEP and the orthogonality threshold
*> is given in TOL.
*> \endverbatim
@@ -147,27 +147,27 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
-*> sequence of Jacobi rotations.
-*> If JOBV = 'N', then MV is not referenced.
+*> If JOBV = 'A', then MV rows of V are post-multipled by a
+*> sequence of Jacobi rotations.
+*> If JOBV = 'N', then MV is not referenced.
*> \endverbatim
*>
*> \param[in,out] V
*> \verbatim
*> V is DOUBLE PRECISION array, dimension (LDV,N)
-*> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
-*> sequence of Jacobi rotations.
-*> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
-*> sequence of Jacobi rotations.
-*> If JOBV = 'N', then V is not referenced.
+*> If JOBV = 'V', then N rows of V are post-multipled by a
+*> sequence of Jacobi rotations.
+*> If JOBV = 'A', then MV rows of V are post-multipled by a
+*> sequence of Jacobi rotations.
+*> If JOBV = 'N', then V is not referenced.
*> \endverbatim
*>
*> \param[in] LDV
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of the array V, LDV >= 1.
-*> If JOBV = 'V', LDV .GE. N.
-*> If JOBV = 'A', LDV .GE. MV.
+*> If JOBV = 'V', LDV >= N.
+*> If JOBV = 'A', LDV >= MV.
*> \endverbatim
*>
*> \param[in] EPS
@@ -187,7 +187,7 @@
*> TOL is DOUBLE PRECISION
*> TOL is the threshold for Jacobi rotations. For a pair
*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
-*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
+*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) > TOL.
*> \endverbatim
*>
*> \param[in] NSWEEP
@@ -205,14 +205,14 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> LWORK is the dimension of WORK. LWORK .GE. M.
+*> LWORK is the dimension of WORK. LWORK >= M.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
diff --git a/lapack-netlib/SRC/dhseqr.f b/lapack-netlib/SRC/dhseqr.f
index 4444b955f..b4fc3af90 100644
--- a/lapack-netlib/SRC/dhseqr.f
+++ b/lapack-netlib/SRC/dhseqr.f
@@ -70,7 +70,7 @@
*> \param[in] N
*> \verbatim
*> N is INTEGER
-*> The order of the matrix H. N .GE. 0.
+*> The order of the matrix H. N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
@@ -87,7 +87,7 @@
*> set by a previous call to DGEBAL, and then passed to ZGEHRD
*> when the matrix output by DGEBAL is reduced to Hessenberg
*> form. Otherwise ILO and IHI should be set to 1 and N
-*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
+*> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim
*>
@@ -100,20 +100,20 @@
*> (the Schur form); 2-by-2 diagonal blocks (corresponding to
*> complex conjugate pairs of eigenvalues) are returned in
*> standard form, with H(i,i) = H(i+1,i+1) and
-*> H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
+*> H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and JOB = 'E', the
*> contents of H are unspecified on exit. (The output value of
-*> H when INFO.GT.0 is given under the description of INFO
+*> H when INFO > 0 is given under the description of INFO
*> below.)
*>
*> Unlike earlier versions of DHSEQR, this subroutine may
-*> explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
+*> explicitly H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1
*> or j = IHI+1, IHI+2, ... N.
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*> LDH is INTEGER
-*> The leading dimension of the array H. LDH .GE. max(1,N).
+*> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim
*>
*> \param[out] WR
@@ -128,8 +128,8 @@
*> The real and imaginary parts, respectively, of the computed
*> eigenvalues. If two eigenvalues are computed as a complex
*> conjugate pair, they are stored in consecutive elements of
-*> WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
-*> WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
+*> WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and
+*> WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in
*> the same order as on the diagonal of the Schur form returned
*> in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
*> diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
@@ -148,7 +148,7 @@
*> if INFO = 0, Z contains Q*Z.
*> Normally Q is the orthogonal matrix generated by DORGHR
*> after the call to DGEHRD which formed the Hessenberg matrix
-*> H. (The output value of Z when INFO.GT.0 is given under
+*> H. (The output value of Z when INFO > 0 is given under
*> the description of INFO below.)
*> \endverbatim
*>
@@ -156,7 +156,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of the array Z. if COMPZ = 'I' or
-*> COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
+*> COMPZ = 'V', then LDZ >= MAX(1,N). Otherwise, LDZ >= 1.
*> \endverbatim
*>
*> \param[out] WORK
@@ -169,7 +169,7 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> The dimension of the array WORK. LWORK .GE. max(1,N)
+*> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient and delivers very good and sometimes
*> optimal performance. However, LWORK as large as 11*N
*> may be required for optimal performance. A workspace
@@ -187,21 +187,21 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0: successful exit
-*> .LT. 0: if INFO = -i, the i-th argument had an illegal
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal
*> value
-*> .GT. 0: if INFO = i, DHSEQR failed to compute all of
+*> > 0: if INFO = i, DHSEQR failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.)
*>
-*> If INFO .GT. 0 and JOB = 'E', then on exit, the
+*> If INFO > 0 and JOB = 'E', then on exit, the
*> remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output
*> value of H.
*>
-*> If INFO .GT. 0 and JOB = 'S', then on exit
+*> If INFO > 0 and JOB = 'S', then on exit
*>
*> (*) (initial value of H)*U = U*(final value of H)
*>
@@ -209,19 +209,19 @@
*> value of H is upper Hessenberg and quasi-triangular
*> in rows and columns INFO+1 through IHI.
*>
-*> If INFO .GT. 0 and COMPZ = 'V', then on exit
+*> If INFO > 0 and COMPZ = 'V', then on exit
*>
*> (final value of Z) = (initial value of Z)*U
*>
*> where U is the orthogonal matrix in (*) (regard-
*> less of the value of JOB.)
*>
-*> If INFO .GT. 0 and COMPZ = 'I', then on exit
+*> If INFO > 0 and COMPZ = 'I', then on exit
*> (final value of Z) = U
*> where U is the orthogonal matrix in (*) (regard-
*> less of the value of JOB.)
*>
-*> If INFO .GT. 0 and COMPZ = 'N', then Z is not
+*> If INFO > 0 and COMPZ = 'N', then Z is not
*> accessed.
*> \endverbatim
*
@@ -261,8 +261,8 @@
*> This depends on ILO, IHI and NS. NS is the
*> number of simultaneous shifts returned
*> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
-*> The default for (IHI-ILO+1).LE.500 is NS.
-*> The default for (IHI-ILO+1).GT.500 is 3*NS/2.
+*> The default for (IHI-ILO+1) <= 500 is NS.
+*> The default for (IHI-ILO+1) > 500 is 3*NS/2.
*>
*> ISPEC=14: Nibble crossover point. (See IPARMQ for
*> details.) Default: 14% of deflation window
@@ -341,8 +341,8 @@
PARAMETER ( NTINY = 11 )
*
* ==== NL allocates some local workspace to help small matrices
-* . through a rare DLAHQR failure. NL .GT. NTINY = 11 is
-* . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
+* . through a rare DLAHQR failure. NL > NTINY = 11 is
+* . required and NL <= NMIN = ILAENV(ISPEC=12,...) is recom-
* . mended. (The default value of NMIN is 75.) Using NL = 49
* . allows up to six simultaneous shifts and a 16-by-16
* . deflation window. ====
diff --git a/lapack-netlib/SRC/dla_gbrcond.f b/lapack-netlib/SRC/dla_gbrcond.f
index e9713c9ca..c9eebcbea 100644
--- a/lapack-netlib/SRC/dla_gbrcond.f
+++ b/lapack-netlib/SRC/dla_gbrcond.f
@@ -141,13 +141,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (5*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/dla_gbrfsx_extended.f b/lapack-netlib/SRC/dla_gbrfsx_extended.f
index 12b2a32a4..282a63f1c 100644
--- a/lapack-netlib/SRC/dla_gbrfsx_extended.f
+++ b/lapack-netlib/SRC/dla_gbrfsx_extended.f
@@ -66,19 +66,19 @@
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
-*> The value is defined by ILAPREC(P) where P is a CHARACTER and
-*> P = 'S': Single
+*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
+*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
-*> = 'X', 'E': Extra
+*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] TRANS_TYPE
*> \verbatim
*> TRANS_TYPE is INTEGER
*> Specifies the transposition operation on A.
-*> The value is defined by ILATRANS(T) where T is a CHARACTER and
-*> T = 'N': No transpose
+*> The value is defined by ILATRANS(T) where T is a CHARACTER and T
+*> = 'N': No transpose
*> = 'T': Transpose
*> = 'C': Conjugate transpose
*> \endverbatim
@@ -270,7 +270,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
diff --git a/lapack-netlib/SRC/dla_gercond.f b/lapack-netlib/SRC/dla_gercond.f
index aa93ca5a4..6f7d70a6a 100644
--- a/lapack-netlib/SRC/dla_gercond.f
+++ b/lapack-netlib/SRC/dla_gercond.f
@@ -123,13 +123,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/dla_gerfsx_extended.f b/lapack-netlib/SRC/dla_gerfsx_extended.f
index 082f810f0..4cb9ef4f9 100644
--- a/lapack-netlib/SRC/dla_gerfsx_extended.f
+++ b/lapack-netlib/SRC/dla_gerfsx_extended.f
@@ -64,19 +64,19 @@
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
-*> The value is defined by ILAPREC(P) where P is a CHARACTER and
-*> P = 'S': Single
+*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
+*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
-*> = 'X', 'E': Extra
+*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] TRANS_TYPE
*> \verbatim
*> TRANS_TYPE is INTEGER
*> Specifies the transposition operation on A.
-*> The value is defined by ILATRANS(T) where T is a CHARACTER and
-*> T = 'N': No transpose
+*> The value is defined by ILATRANS(T) where T is a CHARACTER and T
+*> = 'N': No transpose
*> = 'T': Transpose
*> = 'C': Conjugate transpose
*> \endverbatim
@@ -256,7 +256,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERRS_C is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERRS_C is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERRS_C(i,:) corresponds to the ith
diff --git a/lapack-netlib/SRC/dla_porcond.f b/lapack-netlib/SRC/dla_porcond.f
index 498e707e3..b2f9c4b1e 100644
--- a/lapack-netlib/SRC/dla_porcond.f
+++ b/lapack-netlib/SRC/dla_porcond.f
@@ -113,13 +113,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/dla_porfsx_extended.f b/lapack-netlib/SRC/dla_porfsx_extended.f
index 8c0d6bebd..ece9b00ed 100644
--- a/lapack-netlib/SRC/dla_porfsx_extended.f
+++ b/lapack-netlib/SRC/dla_porfsx_extended.f
@@ -65,11 +65,11 @@
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
-*> The value is defined by ILAPREC(P) where P is a CHARACTER and
-*> P = 'S': Single
+*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
+*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
-*> = 'X', 'E': Extra
+*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] UPLO
@@ -246,7 +246,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
diff --git a/lapack-netlib/SRC/dla_porpvgrw.f b/lapack-netlib/SRC/dla_porpvgrw.f
index 4fe1a1922..8a6f9e1a7 100644
--- a/lapack-netlib/SRC/dla_porpvgrw.f
+++ b/lapack-netlib/SRC/dla_porpvgrw.f
@@ -85,7 +85,7 @@
*> The leading dimension of the array AF. LDAF >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim
diff --git a/lapack-netlib/SRC/dla_syrcond.f b/lapack-netlib/SRC/dla_syrcond.f
index 91d557145..23ed82588 100644
--- a/lapack-netlib/SRC/dla_syrcond.f
+++ b/lapack-netlib/SRC/dla_syrcond.f
@@ -119,13 +119,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/dla_syrfsx_extended.f b/lapack-netlib/SRC/dla_syrfsx_extended.f
index f54d15194..b390600c0 100644
--- a/lapack-netlib/SRC/dla_syrfsx_extended.f
+++ b/lapack-netlib/SRC/dla_syrfsx_extended.f
@@ -67,11 +67,11 @@
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
-*> The value is defined by ILAPREC(P) where P is a CHARACTER and
-*> P = 'S': Single
+*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
+*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
-*> = 'X', 'E': Extra
+*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] UPLO
@@ -255,7 +255,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
diff --git a/lapack-netlib/SRC/dla_syrpvgrw.f b/lapack-netlib/SRC/dla_syrpvgrw.f
index c2e5cb018..5ba03093b 100644
--- a/lapack-netlib/SRC/dla_syrpvgrw.f
+++ b/lapack-netlib/SRC/dla_syrpvgrw.f
@@ -101,7 +101,7 @@
*> as determined by DSYTRF.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim
diff --git a/lapack-netlib/SRC/dla_wwaddw.f b/lapack-netlib/SRC/dla_wwaddw.f
index 99a86c553..4f50540d6 100644
--- a/lapack-netlib/SRC/dla_wwaddw.f
+++ b/lapack-netlib/SRC/dla_wwaddw.f
@@ -36,7 +36,7 @@
*> DLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
*>
*> This works for all extant IBM's hex and binary floating point
-*> arithmetics, but not for decimal.
+*> arithmetic, but not for decimal.
*> \endverbatim
*
* Arguments: